Surface energy and mass balance of Mera Glacier (Nepal, Central Himalaya) and their sensitivity to temperature and precipitation

Abstract

The sensitivity of glacier mass balance to temperature and precipitation variations is crucial for informing models that simulate glaciers’ response to climate change. In this study, we simulate the glacier-wide mass balance of Mera Glacier with a surface energy-balance model, driven by in situ meteorological data, from 2016 to 2020. The analysis of the share of the energy fluxes of the glacier shows the radiative fluxes account for almost all the energy available during the melt season (May–October). However, turbulent fluxes are significant outside the monsoon (June–September). On an annual scale, melt is the dominant mass flux at all elevations, but 44% of the melt refreezes across the glacier. By reshuffling the available observations, we create 180 synthetic series of hourly meteorological forcings to force the model over a wide range of plausible climate conditions. A +1 (−1)°C change in temperature results in a −0.75 ± 0.17 (+0.93 ± 0.18) m w.e. change in glacier-wide mass balance and a +20 (−20)% change in precipitation results in a +0.52 ± 0.10 (−0.60 ± 0.11) m w.e. change. Our study highlights the need for physical-based approaches to produce consistent forcing datasets, and calls for more meteorological and glaciological measurements in High Mountain Asia.

1. Introduction

The pace of climate warming in High Mountain Asia (HMA) is accelerating (Pepin and others, 2015) and precipitation in these regions exhibits significant heterogeneity and remains insufficiently comprehended (Lutz and others, 2014). Furthermore, the intricate relationship between precipitation and temperature variations proves to be a formidable puzzle. This connection, intricately intertwined with glacier mass balance, poses a challenge for understanding the recent evolution of glaciers in HMA (Oerlemans and Reichert, 2000). While multi-year satellite-based estimates of mass changes allow to map the heterogeneity of glacier mass balance across large scales (e.g. Hugonnet and others, 2021), they need to be complemented by other approaches to further elucidate these patterns of contrasted mass losses. One possible approach is to consider that glacier mass changes are the combination of a change in climate conditions modulated by a glacier sensitivity to these changes (e.g. Oerlemans and Reichert, 2000; Marzeion and others, 2012). Following this approach, Sakai and Fujita (2017) demonstrated that regionally different sensitivity to temperature changes could be the main driver of observed mass losses across Asia. They found that the glacier mass-balance sensitivity to temperature was determined by the general climatology, and in particular, by the summer temperature, the annual range of temperature and the ratio between summer and annual precipitation. However, their approach relies on a number of simplifying assumptions that consider only the climate at the equilibrium line altitude (Ohmura and others, 1992) and the climate data they used have a coarse spatial resolution. There is thus room to improve the methodology they applied, in particular through a better representation of processes responsible for glacier mass losses and gains.

These processes controlling the glacier mass are determined by the surface energy balance (SEB), which is commonly modelled to investigate how glacial mass balance is governed and how sensitive it is to climatic variables (Fujita, 2008; Azam and others, 2014; Fugger and others, 2022). There is a long history of studies that investigated the glacier SEB at various locations and at various temporal and spatial scales to relate atmospheric variables to glacier mass changes (e.g. Oerlemans and Knap, 1998; Favier and others, 2004). Specifically, in Hindu Kush Himalaya (HKH), a number of studies investigated the SEB of glaciers in different climate contexts (Mölg and others, 2012; Huintjes and others, 2015a, 2015b; Zhu and others, 2015, 2018, 2021; Fugger and others, 2022; Arndt and Schneider, 2023). They highlight the different sensitivities to temperature and precipitation in different climate conditions, with the dry and cold (continental) climate that prevails in the northwest margin of HKH being associated with low sensitivities of glacier mass balance to temperature, and the warmer and wetter (oceanic) climate of southeast HKH corresponding to larger sensitivities (e.g. Arndt and Schneider, 2023). The larger sensitivities are associated with the prevalence of surface melt in the surface mass balance. Surface energy-based studies find highly non-linear sensitivity of the mass balance to precipitation, unlike studies based on empirical approaches, such as degree-day modelling (e.g. Wang and others, 2019). This is due to the highly non-linear response of glacier surface mass balance to the albedo effect (e.g. Arndt and Schneider, 2023).

However, in HKH most of the SEB studies have two main limitations: either they were conducted at point scale (Kayastha and others, 1999; Azam and others, 2014; Acharya and Kayastha, 2019; Litt and others, 2019; Mandal and others, 2022), or they used meteorological data from reanalysis products (Arndt and Schneider, 2023). The point-scale modelling of the SEB is limited because the SEB is very sensitive to the surface state of the glacier (ice, snow or debris), and to the distribution of meteorological variables (precipitation, temperature, radiative fluxes, etc.) that vary across the glacier area (Oerlemans and others, 1999). Modelling the SEB of a glacier across its entire area requires distributed measurements of meteorological variables, and measurements of the glacier surface mass balance at multiple locations, including the accumulation area. Unfortunately, such data are seldom available in HKH (Huintjes and others, 2015a, 2015b; Arndt and others, 2021; Oulkar and others, 2022; Srivastava and Azam, 2022). Meteorological variables obtained from reanalysis can be heavily biased, especially if they are not downscaled with local measurements (e.g. Hamm and others, 2020; Khadka and others, 2022). When meteorological variables from reanalysis are used to force a glacier mass-balance model, they first need to be debiased, which is often done by tuning a precipitation correction factor until the glacier mass balance matches observations. While there is usually no alternative, this method is known to be subject of equifinality (e.g. Rounce and others, 2020).

This article presents a glacier-wide SEB analysis of Mera Glacier in the eastern part of Central Himalaya. We applied the ‘COupled Snowpack and Ice surface energy and mass balance model in PYthon’ (COSIPY: Sauter and others, 2020) which has been optimised and evaluated using site-specific measurements (Fig. S1 in the Supplementary material). Among Nepal's monitored glaciers, Mera Glacier stands out for its extensive and continuous meteorological and mass-balance data, making it one of the most comprehensively observed glaciers in the region (Wagnon and others, 2021; Khadka and others, 2022). By integrating field measurements, in situ meteorological data and the SEB model, we aim to enhance our understanding of (1) the physical processes governing the seasonal and spatial variability of the glacier mass balance and (2) the sensitivity of the mass balance to meteorological variables. The findings from this comprehensive study will give us a better understanding of the impact of the on-going climate change on Himalayan glaciers.

2. Study area and climate

2.1 Mera Glacier

Mera Glacier, situated in the eastern part of the Central Himalaya within the Upper Dudh Koshi basin, is a plateau-type debris-free glacier. Encompassing an area of 4.84 km² in 2018, the glacier stretches from an elevation of 6390 m a.s.l. to a minimum of 4910 m a.s.l. (Fig. 1). This north-facing glacier features a gentle slope with a mean inclination of ~16°. At an elevation of ~5900 m a.s.l., the glacier separates into two distinct branches, the Mera branch and the Naulek branch. The Mera branch initially heads north and then curves westwards, while the Naulek branch extends ~2 km towards the northeast. The Mera branch is the largest of the two branches and accounts for ~80% of the glacier's total area.

Figure 1. Map of Mera Glacier showing the network of ablation stakes (blue dots) and accumulation pits (cyan diamonds). The stake location and number are taken from November 2020. The number of stakes vary from year to year, due to total excavation, reinstallation at the original location, snow burial or destruction. The pink stars represent the locations of different AWSs with their respective photos and dates ((a) Khare Geonor, (b) Mera La AWS, (c) AWS-H and (d) AWS-L). The outline of Mera Glacier is from 2018 with a total area of 4.84 km², and the background image was acquired by Sentinel-2 on 24 November 2018. Elevation lines are extracted from the 2012 Pléiades DEM (Wagnon and others, 2021). The inset map gives the location of Mera Glacier in Nepal (black square) and the glacierised areas from RGI6 (shaded blue areas).

2.2 Climate

Like other glaciers in Nepal, Mera Glacier is a summer accumulation type glacier, gaining mass mainly from the summer monsoon (June–September) snowfalls brought by the South Asian monsoon system (Wagnon and others, 2013; Thakuri and others, 2014; Shea and others, 2015a). The glacier experiences most of its accumulation and ablation during the monsoon, which makes it a key season to understand the climatic regime of the glacier (Ageta and Higuchi, 1984). From June to September, the average air temperature measured between 2012 and 2020 at 5360 m a.s.l. on the Naulek branch is 0.3°C, and the average precipitation recorded at 4888 m a.s.l. is equal to 570 mm, with an annual precipitation of 818 mm (Khadka and others, 2022). During this season, warm air masses flow from the Bay of Bengal and bring moisture and precipitation in the Himalaya (Perry and others, 2020). In just a few days, marking the start of the post-monsoon (October–November), generally at the beginning of October, meteorological conditions change abruptly to become dry, sunny and increasingly cold and windy. Very occasionally, this season is marked by the intrusion of typhoons in the Himalaya, which bring large amounts of snowfalls above ~4000 m a.s.l. in just a few days, like in October 2013 and 2014 (Shea and others, 2015a). The winter (December–February) is similar but harsher than the post-monsoon with constantly cold, dry and very windy conditions. At Naulek (5360 m a.s.l.), the average air temperature during this season is −10.4°C. The pre-monsoon starts in March and is characterised by progressively warmer, wetter and less windy conditions until the monsoon is totally installed at the beginning of June. The pre-monsoon is then the second wettest season after the monsoon with approximately one-quarter of the annual precipitation on the glacier (Khadka and others, 2022).

3. Data

3.1 Meteorological data

A network of automatic weather stations (AWSs) has been installed and gradually expanded since 2012 in the Mera Glacier catchment at different elevations and on various surfaces (Fig. 1). In this study, we mainly use data from two on-glacier AWSs, one located in the ablation area on the Naulek branch at 5360 m a.s.l., AWS-Low (hereafter referred to as AWS-L; 4 years of data between November 2016 and October 2020), and one located in the accumulation area at 5770 m a.s.l., AWS-High (hereafter referred to as AWS-H; 3 years of data between November 2017 and November 2020; Fig. S2). Both AWSs record air temperature (T), relative humidity (RH), wind speed (u), incoming and outgoing longwave (LW_in and LW_out, respectively) and shortwave (SW_in and SW_out, respectively) radiation (Table 1).

Table 1. List of the different AWS operating on Mera Glacier, or in its vicinity, with their elevations, operating periods, list of sensors and associated meteorological variables used as forcing, optimisation or validation data of the SEB model

T, air temperature; RH, relative humidity; u, wind speed; SW_in, incoming shortwave radiation; SW_out, outgoing shortwave radiation; LW_in, incoming longwave radiation; LW_out, outgoing longwave radiation; P _a, atmospheric pressure and P, precipitation. The numbers in brackets indicate the data gap of the variables (second column) or the uncertainty of each sensor provided by the manufacturer (third column).

*Artificially aspired during daytime.

There are numerous data gaps in both records, due to AWS failure, power shortage during occasional abundant snowfalls covering the solar panels for instance or sensor breakdowns (Table 1; Fig. 2 and see https://glacioclim.osug.fr/). The largest data gap occurred at AWS-L, when the station fell down from 12 December 2017 to 24 November 2018. To fill these gaps, data from the off-glacier Mera La AWS were used. Linear correlation relationships were established each month between the same variables from the two stations from November 2016 to October 2020, at an hourly time step (Fig. S3 and Table S1). The atmospheric pressure (P _a) measured at Mera La AWS is used at AWS-L without any interpolation, as both AWSs are located <2 km apart at almost the same elevation (Table 1; Fig. 1).

Figure 2. Hourly data from 1 November 2016 to 31 October 2020 of (a) air temperature (T), (b) relative humidity (RH), (c) wind speed (u), (d) incoming shortwave radiation (SW_in), (e) incoming longwave radiation (LW_in) at AWS-L, (f) atmospheric pressure (P _a) at Mera La AWS and (g) precipitation (P) at Khare Geonor. Orange shaded areas indicate data gaps at AWS-L, which have been filled by Mera La AWS data using linear interpolation, and the light blue shaded areas in panel (g) visualise the monsoons.

The Khare Geonor station (4888 m a.s.l.) has been installed on 24–25 November 2016, 472 m lower in elevation and ~3 km northwest of AWS-L. In the Mera catchment, this is the only station of the network recording all-weather precipitation, thanks to a weighing device. The precipitation data have been corrected for undercatch following the method by Førland and others (1996) and Lejeune and others (2007) as a function of wind speed and precipitation phase (liquid or solid) depending on air temperature (see the details in the supplement of Khadka and others, 2022).

3.2 Spatial distribution of meteorological forcings

Given the complexity and heterogeneity of the local climate and terrain, it is a challenging task to distribute point data spatially. For air temperature, relative humidity and incoming longwave radiation, we derived empirical linear relationships from the respective measurements of the two on-glacier AWSs installed with a 410 m difference in elevation. To take advantage of the longest data series without gaps at AWS-H, we calculate the gradients between 1 March 2018 and 28 February 2019, even though a major part of the data from AWS-L is reconstructed over this period. Since AWS-L has a relatively longer and more consistent dataset than AWS-H, we distribute meteorological data from this lower station across the glacier using observed altitudinal gradients of air temperature, relative humidity and incoming longwave radiation (Table S2; Fig. S4). In meteorology, the dew-point temperature gradient is more commonly encountered than the relative humidity gradient. However, in this particular study, the relative humidity gradient is utilised because the COSIPY model is developed based on relative humidity input data. For completeness, we also compute the dew-point temperature gradient from data collected at both stations, which we subsequently convert into a relative humidity gradient. This converted gradient is comparable to the one directly obtained from relative humidity measurements at AWS-L and AWS-H, that we use in our study (Fig. S5 and the corresponding Supplementary text).

The incoming solar radiation has been distributed for each grid following the methods of Sauter and others (2020), already tested and applied on Himalayan glaciers (e.g. Arndt and Schneider, 2023). First the fraction of diffuse radiation (F _diff) is calculated based on Wohlfahrt and others (2016):

(1)$$F_{{\rm diff}} = {\rm e}^{-{\rm e}^{\,p1\;-\;( {\,p2\;-\;p3C_{\rm I}} ) \;}}\;( {1-p4} ) \;\;p4$$

where p1 = 0.1001, p2 = 4.7930, p3 = 9.4758, p4 = 0.2465 are parameters from Wohlfahrt and others (2016) and C _I, for clearness index, is the ratio of incoming solar radiation to maximum incoming solar radiation. F _diff may vary from 0 to 1. Second, the measured incoming shortwave radiation is split into beam (R _b = SW_in(1 − F _diff)) and diffuse (R _d = SW_in × F _diff) radiation. Then, the corrected solar radiation (R _c) is calculated on each grid based on Ham (2005):

(2)$$R_{\rm c} = R_{\rm b}\;\;c_{\rm f} + R_{\rm d}$$

where c _f is the correction factor calculated based on the azimuth and the slope of each grid following Ham (2005).

As our network does not allow to assess precipitation variations with elevation over the Mera Glacier catchment, precipitation amounts are assumed constant all over the catchment and equal to Khare Geonor records. Similarly, wind speed is likely spatially variable due to terrain aspect, roughness and heterogeneity but in first approximation, we had no choice but to consider the wind as constant over Mera Glacier and equal to that at AWS-L. These first-order approximations are discussed in Section 6.3.

3.3 Mass-balance data

Mera Glacier has been monitored since 2007 at least once a year in November and its mass-balance series is one of the longest continuous field-based series of the Himalaya. Its glacier-wide mass balance is obtained annually using the glaciological method based on a network of 16 ablation stakes and five accumulation sites on average (Fig. 1). This glacier-wide mass-balance series has been calibrated with the 2012–18 geodetic mass balance (Wagnon and others, 2021). Over the period 2007–23, the mean corrected glacier-wide mass balance is equal to −0.42 ± 0.23 m w.e. a⁻¹, with only four positive mass-balance years out of 16. Our study period 2016–20 was characterised by constantly negative mass-balance years with a mean glacier-wide value of −0.74 ± 0.18 m w.e. a⁻¹, 2017–18 being the most negative year (−0.92 ± 0.16 m w.e. a⁻¹) and 2019–20 being the least negative (−0.49 ± 0.22 m w.e. a⁻¹) (Table 2). Between 2007 and 2023, the glacier has lost ~10% of its surface area.

Table 2. Glacier-wide mass balance for Mera Glacier, point mass balance at AWS-L (obtained by averaging all stake measurements on the Naulek branch between 5300 and 5380 m a.s.l.) and at AWS-H (obtained by averaging stake measurements close to AWS-H, from 5750 to 5790 m a.s.l.), as well as snow depths in the ablation area (annually measured during field campaigns in November) and in the accumulation area (assumed for the model)

An ice density of 900 kg m⁻³ and measured snow densities were used to compute point mass balances (370 kg m⁻³ for Naulek, and 380 to 430 kg m⁻³ for the accumulation area). The error range for point mass balances is the std dev. of all measurements

*Updated from Wagnon and others (2021).

Point mass balances measured at each stake or at each accumulation site can exhibit significant spatial variability depending on factors such as elevation, slope, aspect and wind redistribution. Table 2 provides the annual and mean values of point mass balances over the study period 2016–20 at the two on-glacier AWSs. For AWS-L, it is computed by using all stake measurements available on the Naulek branch between 5300 and 5380 m a.s.l., using a mean measured snow density of 370 kg m⁻³ and an ice density of 900 kg m⁻³. For AWS-H, point mass balances measured at sites located between 5750 and 5790 m a.s.l. in the vicinity of the station are averaged, using depth-averaged snow densities measured during each field campaign (from 380 to 430 kg m⁻³).

4. Methods

4.1 Model description (COSIPY)

In this study, we use the COSIPY model, which is a python-based coupled snowpack and ice SEB model (Sauter and others, 2020). The model is a 1-D multi-layer discretisation of the snowpack/ice column that resolves the energy and mass conservation, and calculates the surface energy fluxes using input meteorological variables as forcings. For spatially distributed simulations, the point model is run independently at each point of the glacier domain, neglecting the lateral mass and energy fluxes. The model's reliability has been validated across distinct contexts and geographical regions (Sauter and others, 2020; Arndt and others, 2021; Blau and others, 2021). The COSIPY model calculates the energy available for melt (Q _M) for each time step and is expressed as

(3)$$Q_{\rm M} = SW_{{\rm net}} + LW_{{\rm net}} + Q_{\rm S} + Q_{\rm L} + Q_{\rm C} + Q_{\rm R}$$

where SW _net is the remaining net shortwave radiation at the surface after penetration inside the snow/ice surface, LW _net is net longwave radiation and Q denotes the other heat fluxes of different subsequent scripts S: sensible, L: latent, C: sub-surface (called ground-heat flux in Sauter and others, 2020) and R: rain in W m⁻². All the fluxes are positive when directed towards the surface and negative away from the surface. When the surface temperature is at the melting point and Q _M > 0, the excess energy is used to melt. Additionally, COSIPY calculates a subsurface melt, that is calculated from the penetration of the incoming shortwave radiation (Sauter and others, 2020). For the rest of the analysis, we refer to total melt as the sum of surface and subsurface melt.

4.1.1 Model settings

The COSIPY model is used in its default configuration. The turbulent fluxes, Q _S and Q _L, are calculated as

(4)$$Q_{\rm S} = \rho _{{\rm air}}C_{\rm P}C_{\rm S}u( {T-\;T_{\rm S}} ) $$

(5)$$Q_{\rm L} = \rho _{{\rm air}}L_{\rm V}C_{\rm l}u( {q-\;q_{\rm S}} ) $$

where ρ _air is air density (in kg m⁻³); C _P is the specific heat of air at constant pressure (in J kg⁻¹ K⁻¹) and L _V is the latent heat of sublimation/vapourisation (in J kg⁻¹). C _S and C _l are the dimensionless transport coefficients calculated using the bulk method with initial roughness lengths taken from Mölg and others (2012) and further calibrated (see Section 4.2), q is the specific humidity of air (in g kg⁻¹), T _S and q _S are the temperature (in °C) and specific humidity at the surface, respectively. The bulk Richardson number has been used to assess the stability correction.

The snowfall is distinguished from liquid precipitation using a logistic transfer function based on Hantel and others (2000) (Fig. S6):

(6)$$n = \displaystyle{1 \over 2}\;\{ {{\rm tanh}\;[ {( {( {T\;-\;T_0} ) -\;T_{00}} ) \;s_0} ] + 1} \} $$

where n is the fraction of snowfall (1 when it is only snow, 0 if only rain and in between 0 and 1 if this is mixed rain and snow), T ₀ is the melting point (0°C), T ₀₀ is the centre for snow transfer function (in °C) and s ₀ is the spread snow transfer function. The ageing/decay of the snowpack's albedo is then based on Oerlemans and Knap (1998), where it depends on the number of days after the last snowfall. Snow density is another important property of snow, particularly for its liquid water content or thermal conductivity. It is obtained by following Essery and others (2013). The default albedo and densification parameterisations used in this study are described in the Supplementary material (see additional text related to the Method section, pp. 9–10).

4.1.2 Model description and initialisation

We run COSIPY at point scale at AWS-L and in a distributed way over 51 glacierised grid points of 0.003° × 0.003° (0.0984 km²) resolution (corresponding to a total glacier-covered area of 5.01 km²), resampled from Copernicus GLO-30 DEM (https://doi.org/10.5270/ESA-c5d3d65) using the glacier outlines in 2012 (5.10 km²) from Wagnon and others (2021). Additionally, the grid elevations were adjusted downwards by 19 m to align with the elevation of AWS-L. The selection of this particular resolution is a compromise between the computational time and a reasonable representation of the topography. The model is run at an hourly timescale, independently for the 4 years of data, from 1 November of one year until 31 October of the following year over 2016–20, without spin-up time. We impose a temperature of 265.16 K at the glacier sole because the glacier is cold-based (Wagnon and others, 2013). The model is initialised with 600 layers of glacier ice topped by a snowpack whose profile is specified by the user. In our case, this profile is determined by the snow depth (Table 2), with each layer having an optimal snow layer height of 0.1 m. The initial snow depth in the ablation zone is kept closest to the observations made every November, corresponding to the snow depth measurement at AWS-L. Above 5750 m a.s.l., the model is initialised with a snow depth that increases by 20 cm per 100 m with altitude, starting at 50 cm on 1 November at 5750 m a.s.l. (Table 2). This assumption is essential because it ensures that there is always snow in the accumulation zone for the simulations. Even though annual observations in November show that the snow line is always lower than 5750 m a.s.l., which is in line with this assumption, the altitudinal gradient is not verifiable in the field. The snow depth in the accumulation zone is probably greater than that used for the initialisation, as the firn–ice interface is several metres below the glacier surface. However, for modelling purposes, the initial snow depth within this zone is of minimal concern as accumulation consistently outweighs ablation, and the snow depth gradually synchronises with the on-going snowfall.

4.2 Optimisation

SEB models are sensitive to the choice of parameter sets (e.g. Zolles and others, 2019). Consequently, model parameters need to be calibrated, and the model's ability to reproduce the glacier surface mass balance needs to be evaluated. We optimise the model following a multi-objective optimisation procedure using forcing data measured at AWS-L between 1 November 2018 and 31 October 2019, the mass-balance year with the least gaps and the most reliable dataset. The optimisation is performed based on the maximum amount of data available, that is, observed albedo, surface temperature calculated from LW_out using the Stefan–Boltzmann equation with a surface emissivity of 0.99 (Blau and others, 2021), and point mass balance at AWS-L (Table 2).

4.2.1 Parameters

There are many parameters in COSIPY and the eight important ones are listed in Table 3. To identify the most sensitive ones among this set, we conducted 108 manual model runs at a point scale, specifically at AWS-L. In these runs, we alternatively and randomly explored various values for selected sensitive parameters, focusing on five parameters related to albedo and three associated with roughness lengths. This rigorous testing encompassed a plausible range of values, allowing us to qualitatively assess their impact on the model's outcomes. These ranges of albedo parameters and roughness lengths are taken from Mölg and others (2012). As we only have one level of wind speed measurement at AWS-L, roughness lengths cannot be directly calculated at this site. Notably, the roughness lengths exhibited lower sensitivity compared to the albedo parameters, which were identified as the most sensitive (in bold in Table 3). These albedo parameters are known sensitive parameters for energy-balance studies (e.g. Zolles and others, 2019). We then optimise these five parameters following the procedure described in Section 4.2.2 starting from a plausible range of values taken from the literature (Mölg and others, 2012; Zolles and others, 2019). All other parameters are taken from the default settings, except for those listed in bold in Table 3, that are optimised.

Table 3. List of selected parameters used in COSIPY, and manually tested before running our optimisation procedure

In bold are the five most sensitive parameters that are optimised from a plausible range of values (min–max range, taken from the Mölg and others (2012) and used in this study. The investigated range of max–min values for the roughness lengths is also shown.

4.2.2 Multi-objective optimisation

Multi-objective optimisation is a calibration method that enables the possibility of more than one optimal solution and provides a way to evaluate a variety of parameter sets (Yapo and others, 1998; Rye and others, 2012; Zolles and others, 2019). The multi-objective approach can be expressed as

(7)$${\rm minimise\;}\{ {\,f_1( \theta ) , \;f_2( \theta ) , \;\;\ldots , \;\;f_n( \theta ) } \} \;$$

where f ₁(θ), f ₂(θ), …, f _n(θ) are n objective functions of model realisations of parameter sets θ. The optimisation process combines multiple objectives into a single ideal through scalar aggregation. For this, a weighted sum (f _agg(θ)) is applied to find the minimum aggregate of different single objectives (Rye and others, 2012):

(8)$${\rm minimise}\;f_{\rm agg}( \theta ) = \{ {w_1f_1( \theta ) , \;w_2f_2( \theta ) , \;\;\ldots , \;w_n\;f_n( \theta ) } \} \;$$

where w is the weight applied to all single objectives based on their performance and the arguments of the aggregating functions to obtain the Pareto solution (Pareto, 1971). With the multiple single objective, the selection of Pareto solution and multi-objective is more precise. Here, the multi-optimisation is done based on three objective metrics that compare the observation at AWS-L and model outputs:

(9)$$f_1( \theta ) = \left[{\displaystyle{{\mathop \sum \nolimits_{i = 1}^n ( {x_i-\bar{x}} ) ( {y_i-\bar{y}} ) } \over {\sqrt {\mathop \sum \nolimits_{i = 1}^n {( {x_i-\bar{x}} ) }^2} \sqrt {\mathop \sum \nolimits_{i = 1}^n {( {y_i-\bar{y}} ) }^2} }}} \right]^2$$

(10)$$f_2( \theta ) = \displaystyle{1 \over N}\left\vert {\mathop \sum \limits_{i = 1}^N ( {y_i-x_i} ) } \right\vert $$

(11)$$f_3( \theta ) = \vert {( {MB_{{\rm obs}}-MB_{{\rm mod}}} ) } \vert $$

where f ₁ is the coefficient of determination (r ²) between the observed (x _i) and modelled (y _i) values of albedo and surface temperature (T _s), f ₂ is the mean absolute error (MAE) between the observed (x _i) and modelled (y _i) values of both variables and f ₃ is the absolute error (AE) between observed (MB _obs) and modelled (MB _mod) point mass balances.

The simplest way to calibrate the model is to optimise the different objective functions for four study periods (4 years × 3 objective functions = 12 performances). This approach results in a large range of uncertainty with many sensitive parameters for different mass-balance years. Similarly, the set of parameters optimised for one period may not perform better in another period, resulting in higher uncertainty (Soon and Madsen, 2005). However, by simulating Pareto solutions for individual mass-balance years and evaluating the objective functions over the other years, it is possible to select the best set of parameters.

4.3 Evaluation of the model at point and distributed scale

First, we optimise the five sensitive model parameters highlighted in Table 3 at AWS-L for the 2018/19 period. Second, we evaluate the model performance at point scale at AWS-L site for the other three mass-balance years, systematically comparing measured and simulated albedo, surface temperature and point mass balance. Third, we run the model in a spatially distributed way and similarly evaluate it at point scale at AWS-H site, over the 3 years of available data (2017–20). Finally, we also compare simulated mass balance to measured surface mass balance at glacier-wide scale for each year between 2016 and 2020. This allows us to test both the spatial and the temporal transferability of the model optimised parameters.

4.4 Sensitivity analysis

Our study aims to evaluate the sensitivity of the surface mass balance to changes in meteorological forcings. In previous studies, the sensitivity to temperature or precipitation is usually assessed by a constant change in temperature (e.g. ±1–2°C) or a relative change in precipitation (e.g. ±10–30%), keeping all other meteorological variables unchanged at the same time (Kayastha and others, 1999; Mölg and others, 2012; Sunako and others, 2019; Arndt and others, 2021; Gurung and others, 2022; Srivastava and Azam, 2022; Arndt and Schneider, 2023). The main disadvantage of this method, hereafter referred to as the classical method, is that perturbing a single meteorological variable breaks the physical links between the meteorological variables, which is detrimental to the simulation of surface mass balance (e.g. Prinz and others, 2016; Clauzel and others, 2023). To overcome this issue, a number of methods were developed to perturb meteorological data while preserving the link between variables (Sicart and others, 2010; Prinz and others, 2016; Autin and others, 2022). Here we both perturb the meteorological forcings in a classical way and we produce synthetic scenarios, described below (see also Fig. S7 for a flow chart explaining how scenarios are obtained).

4.4.1 Classical scenarios

For the classical sensitivity analysis method, the scenarios are developed by varying temperature by ±1°C and precipitation amount by ±20% for four mass-balance years without changing the other meteorological variables. In total, we produce 16 runs (+1°C, −1°C, +20% and −20% for each of the 4 years) between 2016 and 2020. In order to determine the sensitivity of the mass balance to temperature or precipitation, it is then simply necessary to calculate the average anomaly in the mass balance over the 4 year period for each corresponding perturbation.

4.4.2 Synthetic scenarios

Utilising the initial 4 year dataset, we performed a cyclic selection process, systematically considering each month of the year commencing with November and concluding in October of the subsequent year. During each cycle, we alternately chose data from the unaltered original dataset and data from the warmest, coldest, wettest or driest month within the 4 year period (2016–20). This process allowed us to generate synthetic annual datasets at an hourly resolution. By this way, we create 180 one year-long synthetic meteorological series, exploring a wide range of climatic variability from very warm to very cold conditions, and from very dry to very wet conditions (see details below). For each scenario, for each month, we keep original hourly data and we only shuffle the months from different years with each other. In this way, both the physical integrity between meteorological variables and the weather conditions with respect to the time of the year are preserved.

First, the four most extreme synthetic scenarios are developed by making 1 year-long series of hourly forcing data that contain the most extreme months among the 4 years of data. For instance, the wet scenario (hereafter referred to as We_We) is obtained by combining the wettest (We) months: if November 2017 is the month with the maximum monthly amount of precipitation among the 4 months of November of the 2016–20 series, all hourly forcing data from November 2017 are selected in scenario We_We; then if December 2019 is the wettest of the four December months, then hourly data from December 2019 will be selected as the second month of scenario We_We; etc. until October. In this way the wet scenario We_We combining the hourly data of all meteorological variables from the 12 months with the maximum monthly amount of precipitation is created. Similarly, we combine the hourly data of all meteorological variables from the 12 driest months to create the dry scenario (D_D), or from the 12 warmest months to create the warm scenario (Wa_Wa), or from the 12 coldest months to create the cold scenario (C_C).

Second, starting from these four extreme scenarios referred to hereafter as baseline conditions, we also create 48 additional scenarios by modifying some specific seasons. For instance, we keep the conditions of the wet scenario We, except during one season, let's say the monsoon, where we decide to take the conditions of the warm scenario Wa. The four considered seasons are winter (win), pre-monsoon (pre), monsoon (mon) and melting season (melt). We decide not to consider the post-monsoon which is not critical for the glacier mass balance because there is usually neither any large precipitation nor any large melt but we prefer to introduce a 6 month-long melting season which we suspect to be more critical to control the glacier mass balance. This melting season covers half of the year from May to October, the only months where significant melt is observed in the field. The synthetic scenarios follow the naming convention AB_X, where A is the climate baseline coming from the four extreme scenarios, so A is either We, D, C or Wa, B characterises the conditions of the modified season, that is, B is also either We, D, C or Wa and X refers to the season that has been modified to create the additional scenario, so X is either win, pre, mon or melt. In case of the specific example described above, the scenario is named We_Wa_mon which means that the dataset is 1 year-long, using hourly forcing data from the wet scenario except during the monsoon where the warmest months are selected. A can take four values, B can take three values and X can take four values (the number of values for B cannot be similar to that of A, otherwise we define one of the extreme scenario), resulting in a total of 4 × 3 × 4 = 48 scenarios.

Similarly, we also mix the original unshuffled data (U) from each year of the study period with extreme months. U can be alternatively 2016/17, 2017/18, 2018/19 or 2019/20 respectively referred to as 2017, 2018, 2019 or 2020. These scenarios follow the same naming convention, with U being an additional value for A and B. For example, the scenario 2018_Wa_mon corresponds to the unchanged original data of the year 2017/18, except during the monsoon where data from the warmest months are selected. We obtain here 4 × 4 × 4 = 64 scenarios. Similarly, D_2020_melt corresponds to the driest months, where months of the melting season have been replaced by the original 2019/20 data. Again, we obtain here 4 × 4 × 4 = 64 scenarios. In total, we have 4 + 48 + 64 + 64 = 180 scenarios.

For each of these 180 synthetic annual datasets, we calculate the glacier-wide annual mass balance of Mera Glacier using COSIPY and calculate the difference from the mean annual glacier-wide mass balance simulated by COSIPY using the original 2016–20 dataset. We call this difference the mass-balance anomaly. Similarly, we have an anomaly for each forcing meteorological variable. We can now derive the mass-balance sensitivity to each meteorological variable by fitting a linear regression between the mass-balance anomaly and the anomaly of the variable under consideration. The slope of this linear relationship gives the mass-balance sensitivity to the given variable. The error bars on these sensitivities for ±1°C temperature and ±20% precipitation variation have been calculated using the 99% confidence interval, given by three std dev.s, of this linear regression.

5. Results

5.1 Optimisation and evaluation

5.1.1 Optimisation

The point-scale optimisation strategy at AWS-L results in a well identified group of parameter sets (Pareto solutions) that minimise multiple objective functions simultaneously (Fig. 3). The objective functions of each set of original parameters are largely scattered with the range of r ² and MAE of albedo being 0.06–0.54 and 0.11–0.25, respectively (Fig. 3a). Similarly for surface temperature, r ² and MAE are in-between 0.81–0.92, and 1.37–3.24°C, respectively (Fig. 3d). The range of AE is 0–3.92 m w.e. for the point mass balance at AWS-L (Figs 3b, c, e, f). The 200 solutions closest to the utopia point in terms of how well they perform across all objective functions represent the Pareto solutions (Fig. 3: bold black dots). All these multiple Pareto solutions are almost equally good and plausible. The metrics of the Pareto solutions are less scattered than the original ensemble, with the range of r ² and MAE of albedo being 0.31–0.54 and 0.11–0.15 (Fig. 3a), and the range of r ² and MAE of surface temperature being 0.87–0.92 and 1.58–2.33°C (Table 4). The range AE is 0–1.28 m w.e. for the point mass balance at AWS-L (Table 4).

Figure 3. Solution space for the multi-objective optimisation for the period 1 November 2018 – 31 October 2019. One dot represents results obtained with one set of parameters, and bold black and red dots define the Pareto solution space and optimised solution, respectively. Plots show the scatter plot between (a) 1 − r ² and MAE from albedo comparison, (b) MAE from albedo comparison and AE from mass-balance comparison, (c) AE of mass-balance comparison vs 1 − r ² from albedo, (d) 1 − r ² and MAE from surface temperature comparison, (e) MAE from surface temperature comparison and AE from mass balance and (f) AE of mass balance vs 1 − r ² from surface temperature comparison.

Table 4. Range of different objective function values in the first 200 Pareto solution space for 2018/19 period at AWS-L

Since Pareto solutions perform well over all time and space ranges, the best parameter set among these Pareto solutions is then chosen as the optimised set (Table 3). The r ² and MAE between the observed albedo and that resulting from the selected optimised parameter set are 0.48 and 0.12, respectively; similarly, for surface temperature r ² and MAE are 0.91 and 1.92°C, respectively (Fig. 4), and AE for the point mass balance is 0.17 m w.e. This final optimised set of parameters corresponds to a time scaling factor of 3 days and a depth scaling factor of four centimetres for the albedo model (Oerlemans and Knap, 1998) as well as values of albedo of snow (0.85), firn (0.55) and ice (0.30) similar to the default values used in COSIPY (Table 3).

Figure 4. Mean daily snow albedo (top) and surface temperature (bottom) from observation (Obs., black line) and simulated with COSIPY between 1 November 2018 and 31 October 2019 at AWS-L. r ² and MAE represent the correlation coefficient and MAE between the observed and simulated variables, respectively. The red thick line and the brown thin lines represent the simulated variables using the final optimised parameter set and using all other solution spaces, respectively.

5.1.2 Glacier-wide simulation of COSIPY

The good performance of COSIPY simulations at point scale with the optimal set of parameters does not guarantee the good performance of the distributed simulations that rely on additional hypotheses, such as the meteorological forcing distribution that changes the meteorological forcings even at AWS-L location, because for instance SW_in is re-computed at AWS-L based on the slope and the aspect of the considered gridcell. We thus evaluate the distributed COSIPY simulations over the period 2016–20 with the albedo and surface temperature at AWS-L and AWS-H, and with the glacier-wide mass balance. r ² (MAE) for albedo is 0.32 (0.15) and 0.16 (0.14) for 2016/17 and 2019/20 at AWS-L, respectively (Fig. S8). At AWS-H, over the 3 year period 2017–20, r ² for albedo and surface temperature are 0.50 and 0.92 respectively, and the mean AE for 2016–20 at AWS-H is only 0.15 ± 0.14 m w.e. The surface temperature is always highly correlated with a low bias in both sites (Fig. S8). These metrics are close to the ones from point-scale simulations.

Additionally, we compare the observed and simulated surface point and glacier-wide mass balances. The simulated point surface mass balances match well with the in situ measurements obtained at stakes for all the 4 years (Figs 5 and 6). The location of the equilibrium line altitude is well represented in COSIPY simulations and the general shape of the dependency of the surface point mass balance on elevation is satisfyingly reproduced (Fig. 5). The glacier-wide mass balance from the model is the mean value from all 51 individual gridcells and it is compared to the in situ glacier-wide mass balance taken from Wagnon and others (2021). Over the 4 year period, the mean observed glacier-wide in situ mass balance is −0.74 ± 0.18 m w.e. a⁻¹ (Table 2) and the simulated mass balance is −0.66 m w.e. a⁻¹ with a std dev. of ±0.26 m w.e. a⁻¹. The largest difference between the observed and modelled glacier-wide mass balance occurs in 2019/20, with the simulated mass balance being 0.22 m w.e. less negative than the observed one (Fig. 5).

Figure 5. In situ (blue dots) and simulated at each gridcell (red dots) point mass balances as a function of elevation on Mera Glacier for each year of the 2016–20 period. MB is the glacier-wide mass balance obtained from field measurements (blue text) (Wagnon and others, 2021) and simulated with COSIPY (red text). Also shown are the hypsometries of Mera Glacier used for in situ glacier-wide mass-balance calculations (light blue histograms) and for COSIPY (light brown histograms).

Figure 6. Distributed simulated annual mass balance (MB, in m w.e.) for each year of the study period. Also shown as white circles are the point mass-balance observations (ablation stakes and accumulation pits) with the inside colour corresponding to the respective annual measurement. The glacier outlines (black) are from Wagnon and others (2021).

5.2 SEB and mass-balance components

5.2.1 Seasonal and annual energy-balance components

Figure 7 shows the monthly glacier-wide surface energy and mass-balance components at Mera Glacier for the period 2016–20, and Figs S9–S16 are maps of the glacier, showing the distributed annual energy and mass fluxes for each year of the study period. SW_net is the primary energy source available at the surface throughout the year, the second energy source being the Q _S, that is significant only between November and March. Along the year, SW_in is controlled by the position of the sun responsible for the potential SW_in and also by cloudiness, explaining why it decreases from 274 W m⁻² in the pre-monsoon to 195 W m⁻² during the monsoon (Table 5). Similarly, SW_net decreases from 83 W m⁻² during the pre-monsoon to 61 W m⁻² during the monsoon because of SW_in reduction rather than change in albedo (glacier-wide values of 0.71 during the pre-monsoon and 0.70 during the monsoon). The change in air temperature and water vapour (moisture) is responsible for a strong increase of LW_in from the pre-monsoon (217 W m⁻²) to the monsoon (296 W m⁻²), the only season when LW_in nearly counterbalances the LW_out.

Figure 7. Glacier-wide monthly (a) energy fluxes, (b) mass fluxes and (c) mass balance from November 2016 to October 2020 on Mera Glacier (left panels) and mean monthly annual cycle (right panels). SW_net, net shortwave radiation; LW_net, net longwave radiation; Q _L, latent heat flux; Q _S, sensible heat flux; Q _C, subsurface heat flux; Q _R, rain heat flux; Q _M, available melt energy at the surface; SnowF., solid precipitation; Subl., sublimation; Surf. M., melt at surface; Sub S. M., subsurface melt and Refr., refreezing. Blue shaded areas visualise the monsoons.

Table 5. Annual and seasonal surface energy fluxes (W m⁻²) and their contribution to the total energy intake (Q _in) and outtake (Q _out) over the whole Mera Glacier area, at AWS-L and at AWS-H

The annual values are calculated between 1 November and 31 October of the following year, using all data over the study period. The negative (−) sign indicates the energy loss from the surface.

Winter, Dec–Feb; pre-monsoon, Mar–May; monsoon, Jun–Sep; post-monsoon, Oct–Nov.

The total energy intake at the surface is highest and almost similar during the pre-monsoon (496 W m⁻²) and the monsoon (491 W m⁻²) (Table 5). However, the net all-wave radiation, calculated as the sum of SW_net and LW_net, is 32 W m⁻² during the pre-monsoon and 20 W m⁻² higher during the monsoon. This indicates that the change in cloud cover and atmospheric condition has a relatively minor effect on the total energy absorbed at the glacier surface, but does impact the net all-wave radiation. In the pre-monsoon, this net all-wave radiation is equally compensated by Q _L and Q _M (~−16 W m⁻² each). During the monsoon, LW_net and Q _L are both reduced or close to zero leaving all the energy available for melt with an average energy value of −43 W m⁻². The contributions of Q _S and Q _C are always low during the pre-monsoon and the monsoon, while the Q _R is negligible all the time.

The energy-balance components vary across different glacier areas; they are analysed in the ablation area at AWS-L and in the accumulation area at AWS-H. When considering the annual means, the magnitudes of SW_net and Q _M are higher at AWS-L (93 and −34 W m⁻², respectively) than at AWS-H (61 and −15 W m⁻², respectively). LW_net, Q _L and Q _S exhibit similar annual means throughout the year at both sites. At AWS-L, SW_net remains similar during the pre-monsoon (87 W m⁻²) and the monsoon (88 W m⁻²) because the decrease in SW_in (287 W m⁻² in the pre-monsoon and 201 W m⁻² in the monsoon) is compensated by a decrease in albedo (0.70 in the pre-monsoon and 0.56 in the monsoon). The albedo remains high at AWS-H during the whole year, and there is thus a decrease of SW_net in the monsoon (43 W m⁻²) compared to the pre-monsoon (76 W m⁻²). The variation of LW_net at AWS-L and AWS-H is rather small as the difference between LW_in and LW_out remains similar. Comparing both sites, Q _L remains rather similar during all seasons. Q _M dominates Q _L all year round except during the winter at AWS-L, but at AWS-H, Q _L always dominates Q _M, except during the monsoon. Q _S is significant only during the cold months of the winter and the post-monsoon, with a similar magnitude whichever the location. Q _C is positive and rather small (<5 W m⁻², slightly higher during the post-monsoon) all year round except during the monsoon at AWS-L where it is slightly negative (Table 5).

5.2.2 Seasonal and annual mass-balance components

Table 6 lists the annual and seasonal mass-balance components of Mera Glacier. After the direct accumulation (through snowfalls) and surface melt on Mera Glacier, annual refreezing and sublimation are two major mass-balance components, refreezing being even higher than snowfalls. Indeed, at glacier scale, 44% of the total (surface + sub-surface) melt refreezes annually. The glacier-wide sublimation is −0.15 m w.e. and therefore contributes 23% of the total mass balance or 6% of the ablation terms (total melt + sublimation).

Table 6. Glacier-wide annual and seasonal mass-balance components (mm w.e.) over the total glacier area and at point scale at AWS-L and AWS-H using all data over the study period 2016–20

The negative (−) sign indicates a mass loss from the surface

Looking at seasonal scale, pre-monsoon and monsoon are important seasons in terms of mass-balance processes, as more than 86% of solid precipitation falls and 84% of annual melt happens from March to September. However, post-monsoon is not completely negligible in terms of melt (−0.30 m w.e. or 14% of the annual melt). This melt is almost equal to the surface mass balance (−0.17 m w.e.) due to the limited magnitude of the other processes, and in particular the limited refreezing in the snow-free areas of the glacier. The winter is characterised by limited mass-balance processes, with ~11% of the annual solid precipitation and 2% of the annual total melt (Table 6).

At glacier scale, the total melt (−0.42 m w.e.) and sublimation (−0.05 m w.e.) during the pre-monsoon are balanced by the refreezing (0.30 m w.e.) and snowfall (0.14 m w.e.), leading to a near zero surface mass balance (Table 6). The glacier-wide total melt during the monsoon is −1.42 m w.e., which is higher at AWS-L (−2.15 m w.e.) and lower at AWS-H (−0.99 m w.e.). Depending on the presence and the state of a snowpack (snow depth, density and temperature), meltwater refreezes below the surface. Glacier-wide refreezing is 0.49 m w.e. during the monsoon, it is lower at AWS-L (0.19 m w.e.) where ice is often exposed at the surface, and higher at AWS-H (0.70 m w.e.) where there is always a snowpack with negative temperature. The refreezing preserves 35% (0.49 m w.e) of the total glacier-wide melt during the monsoon; its relative contribution is higher at AWS-H (71%) than at AWS-L (9%).

The annual glacier-wide sublimation is −0.15 m w.e. and is nearly identical at both AWSs (−0.16 and −0.15 mm w.e. at AWS-L and AWS-H, respectively; Fig. S6). Most of the sublimation (93% glacier-wide) happens outside the monsoon, when cold, dry and windy conditions prevail. Wind is not spatially distributed in our simulations, leading to rather homogeneous sublimation across the glacier. There are few exceptions, like a slightly higher sublimation in winter at AWS-L (−0.07 m w.e.) than at AWS-H (−0.05 m w.e.) due to higher roughness length linked to the surface state (exposed ice at AWS-L vs snow at AWS-H) and to the mixing ratio that depends on the air temperature. Due to the lower wind speed, sublimation is insignificant at both AWS sites during the monsoon (Table 6).

5.3 Mass-balance sensitivity to meteorological forcings

5.3.1 Link between meteorological forcing anomalies and mass-balance anomalies

In order to analyse the link between the different input meteorological variables and the outputs from the simulations, we calculate anomalies of each variable from the 180 scenarios by subtracting the mean of the original unshuffled 2016–20 simulations (Fig. 8). From all synthetic scenarios, the magnitude of variation of air temperature is nearly ±1°C, whereas for precipitation, it varies from −35 to +55% annually. The anomalies of relative humidity (±5%) and incoming radiations (±10 W m⁻²) are rather narrow. There are many significant correlations (p < 0.001) between the anomalies of the different variables, suggesting that they are likely to be physically related to each other. On an annual scale, the anomaly of air temperature correlates significantly and positively with the wind speed and the air pressure anomalies, and negatively with the precipitation and relative humidity anomalies. Regarding radiations that are expected to have an impact on the mass balance, the SW_in anomaly correlates significantly and negatively with the relative humidity, the precipitation and the LW_in anomalies. The LW_in anomaly correlates significantly and positively with the relative humidity and the precipitation anomalies, and negatively with the SW_in and wind speed anomalies.

Figure 8. Scatter plot between anomalies of different input variables and glacier-wide mass-balance anomalies for the 180 synthetic runs. Also shown are the Pearson correlation coefficients between the series of annual anomalies. The black lines represent the linear regressions and the grey shaded areas indicate the standard error. The anomalies of each variable are calculated by subtracting the mean of the original unshuffled 2016-20 simulation. (Asterisks represent significance levels, accordingly: *0.05, **0.01, ***0.001.)

The mass-balance anomaly correlates significantly with the anomalies of every meteorological variable except for SW_in (Fig. 8). The correlations between mass-balance anomalies and those of LW_in, atmospheric pressure or wind speed are moderate but significant, positive in case of LW_in, and negative for the other variables. The correlation between mass-balance anomalies and air temperature is significant and highly negative (r = −0.79). Mass-balance anomalies are highly and positively correlated with those of precipitation (r = 0.87) and relative humidity (r = 0.84; Fig. 8).

5.3.2 Mass-balance sensitivity to air temperature and precipitation

From the classical method, we find that perturbing the temperature by +1 (−1)°C leads to a change in glacier-wide mass balance of −0.61 (+0.41) m w.e. (Table 7). A −20 (+20)% change in precipitation leads to a −0.79 (+0.48) m w.e change in glacier-wide mass balance (Table 7). With the synthetic scenario method, we find that a temperature change of +1 (−1)°C leads to a glacier-wide mass-balance change of −0.75 ± 0.17 (+0.93 ± 0.18) m w.e., and a −20 (+20)% change in precipitation results in a mass-balance change of −0.60 ± 0.11 (+0.52 ± 0.10) m w.e. Due to the physical link between variables, and in particular the negative correlation between temperature and precipitation, we find that the sensitivity of mass balance to temperature is significantly higher when calculated from synthetic scenarios than from the classical method, especially in case of cooling. For precipitation, it is significantly reduced in case of a precipitation deficit but almost unchanged in case of an increase.

Table 7. Mass-balance anomalies as compared to the mean of the four 2016–20 years from the classical and synthetic scenarios’ methods

This synthetic scenario approach allows to derive mass-balance sensitivities to any meteorological variable, as long as a significant correlation exists between the anomalies of mass balance and the variable under consideration. In particular, as we find a high correlation between mass balance and relative humidity anomalies, we can also assess the mass-balance sensitivity to this variable: a −4 (+4)% change in RH corresponds to a −1.02 (+1.38) m w.e. change in mass balance. However, we caution on the interpretation of these correlations, as most of the input meteorological variables covary, the correlations may be significant, but they do not show a causal relationship, that needs to be discussed in the light of the knowledge of the processes (see Section 6).

5.3.3 Specific meteorological conditions leading to the most positive/negative mass balances

From the synthetic scenarios, the annual glacier-wide mass balances range from −1.76 to 0.54 m w.e. (Figs 9 and S17), which is a wider range than the historically measured glaciological mass balance since 2007 (min = −0.92 m w.e. in 2017/18 and max = 0.26 m w.e. in 2010/11; Wagnon and others, 2021). The simulated annual mass balances are compared with the original mean annual glacier-wide mass balance of −0.66 m w.e. (Table 6) over the 2016–20 study period, referred to hereafter as the reference year. Most of the scenarios that have warm or dry conditions as a baseline correspond to the first category of scenarios characterised with negative mass balances ranging from −1.76 to −0.81 m w.e. They have a positive SW_net anomaly compared to the reference year (+2 to +17 W m⁻²), associated either with a change in air temperature towards a warming (−0.71 to +1.13°C) or to a decrease in snowfall (0 to −0.29 m w.e.) or to both, resulting in a low glacier-wide albedo (0.53–0.64). With more energy intake, melting is enhanced, and due to the reduced snowfalls, ice is more exposed at the glacier surface favouring runoff, and in turn <46% of this meltwater refreezes, ultimately leading to the most negative glacier-wide mass balances (Fig. S17).

Figure 9. Glacier-wide (a) energy fluxes, (b) mass flux components and (c) mass balance (MB) from (left panels) 12 selected synthetic scenarios (in red, grey, and blue, on the x-axis), as well as (right panels) from the four mean classical scenarios (in black) and the reference year (RY, in green, on the x-axis). The results from the classical scenarios or the reference year have been averaged over the 4 years 2016–20. Based on the MB results, 12 synthetic scenarios are selected (four corresponding to the most negative MBs, four from the middle of the MB set with moderately negative MBs, and four most positive MBs) out of the 180 scenarios (all shown in Fig. S17). The colour code of synthetic scenarios visualises the MB range, from the most negative (red), to the most positive (blue), grey being intermediate and moderately negative. 2019 and 2020 represent the 2018/19 and 2019/20 mass-balance years, respectively. SW_net, net shortwave radiation; LW_net, net longwave radiation; Q _L, latent heat flux; Q _S, sensible heat flux; Q _C, subsurface heat flux; Q _R, rain heat flux; Q _M, available melt energy at the surface; SnowF., solid precipitation; Subl., sublimation; Surf. M., melt at surface Sub S. M., subsurface melt and Refr., refreezing.

The second category of scenarios corresponds to glacier-wide mass balances from −0.80 to −0.25 m w.e. close to that of the reference year. Here, we find scenarios combining a baseline and a seasonal component that would normally lead to opposite mass-balance responses such as dry with wet conditions or warm with cold conditions (e.g. Wa_We_mon, D_We_mon, Wa_We_win, D_C_win). We also find the majority of scenarios that have the unperturbed data as baseline (Figs 9 and S17). The mass balance is affected equally but in an opposite direction by the temperature anomalies (−0.97 and +1.02°C) and the precipitation anomalies (−0.16 to +0.23 m w.e.). In this category, the refreezing ranges from 37 to 56% of total melt, which is close to that of the reference year (44%), and the ranges of SW_net (−8 to +6 W m⁻²) and LW_net (−4 to +3 W m⁻²) anomalies are small (Fig. 8).

The third category corresponds to positive or near-balanced glacier-wide mass balances (>−0.25 m w.e.) mostly produced by scenarios with wet or cold baselines. They have temperature anomalies between −1.00 and +0.52°C and precipitation anomalies between −0.04 and +0.45 m w.e. (Figs 9 and S17). The higher amount of snowfall increases the accumulation, increases the albedo, and in turn decreases the SW_net (−19 to −4 W m⁻²). In addition, the refreezing is high (46–71% of total melt). Overall, the scenarios with a wet year baseline always create a mass balance that is close to balance, and specifically, the highest positive mass balance is produced by the wettest conditions all year round (scenario We_We) (Figs 9 and S17).

For all scenarios, we find that the mass balance is primarily influenced by the baseline conditions, and not by the seasonal variation. We do not find any season that has an influence larger than the other ones (Fig. S17). In particular, when we look at the scenarios of unperturbed data with seasonal variations, we find that winter seems to have as much influence, if not even more influence, than the other seasons on the mass balance (Fig. S17). This result is rather counterintuitive, as most of the mass-balance processes happen in monsoon and pre-monsoon (e.g. Fig. 8).

Regarding the classical scenarios, as expected, both −20% and +1°C scenarios produce negative mass balances, but less extreme than the dry and warm synthetic scenarios. In contrast, both +20% and −1°C classical scenarios are characterised by near-balanced mass balances, far from the positive glacier-wide mass balances obtained with the wet and cold scenarios. The energy and mass fluxes in the classical scenarios are also comparable to those of the synthetic scenarios. The LW_net is similar in all classical scenarios. However, SW_net is 13 W m⁻² higher and refreezing is 25% lower in the −20% precipitation and +1°C scenarios than those in the +20% and −1°C scenarios.

6. Discussion

6.1 Surface energy and mass-balance components of Mera Glacier, and comparison with other similar studies in HKH

It is difficult to compare different glacier surface energy and mass-balance analyses rigorously across the same region because study periods are never similar. Moreover, temporal (multi-annual, annual or seasonal) and spatial (point scale or glacier-wide) resolutions are often different and not comparable (Table 8). In the HKH region, the seasonality of precipitation has a strong impact on the energy and mass-balance components. In the western Himalaya, the winter precipitation dominates the annual accumulation, and sublimation strongly contributes to the ablation processes (Mandal and others, 2022; Oulkar and others, 2022; Srivastava and Azam, 2022). In the Central Himalaya, the glaciers are summer accumulation type with significant longitudinal variability in mean summer temperature, which has the strongest impact on mass-balance sensitivity (Sakai and Fujita, 2017). Still, precipitation, which depends on the monsoon intensity and duration, is a key variable governing the energy and mass balance of glaciers through the albedo effect and its control on the refreezing (Shaw and others, 2022).

Table 8. Comparison of SEB components on different glaciers in HMA, whose location is visible in Figure S18

The various studies already included in Azam and others (2014) are not reported in this table. All the fluxes are in W m⁻².

*Corresponding to Q _S + Q _L here.

The pattern of surface energy and mass balance over the whole Mera Glacier confirms what has already been observed on other glaciers of HKH (Table 8, which is an update of Table 4 of Azam and others, 2014, with the location of glaciers given in Fig. S18). Overall, the radiative fluxes strongly dominate the SEB, and control the amount of energy available for melt (Fig. 7). Between the pre-monsoon and the monsoon, with the gradual establishment of the dense cloud cover typical of this latter wet and warm seasons, incoming shortwave radiation gradually loses intensity, replaced at the same time by increasing incoming longwave radiation, resulting in a constantly high amount of radiative energy available for the glacier. This amount of energy decreases with elevation mainly because albedo is higher in the accumulation zone, where ice is never exposed at the surface. Turbulent fluxes are only significant out of the monsoon, contributing to bring additional energy towards the surface by sensible heat flux (Fig. 7). Contrary to glaciers in the western Himalaya where resublimation occasionally occurs during the monsoon (Azam and others, 2014), on Mera Glacier, the latent heat flux is always negative which means that sublimation is a non-negligible process of mass loss, especially during seasons with strong wind (winter and post-monsoon), like on Zhadang and Parlung No. 4 glaciers on the southeast Tibetan Plateau (Mölg and others, 2012; Sun and others, 2014; Zhu and others, 2018). On Mera Glacier, glacier-wide sublimation accounts for 23% of the glacier-wide mass balance (Table 6), in agreement with previously published values in the Himalaya (e.g. Gurung and others, 2022; Srivastava and Azam, 2022) or lower than that of Sutri Dhaka Glacier in the western Himalaya (55%; Oulkar and others, 2022). However, this value on Mera Glacier is under-estimated because wind speed strongly increases with elevation (Khadka and others, 2022), an effect that has not been taken into account in our study. Refreezing is also an important process, because at annual and glacier-wide scales, as much as 44% of meltwater refreezes on Mera Glacier, with a clear increase of this percentage with elevation (Table 6). This finding is again in agreement with other studies in the region (Stigter and others, 2018; Bonekamp and others, 2019; Kirkham and others, 2019; Veldhuijsen and others, 2022; Arndt and Schneider, 2023). It is noteworthy mentioning that sublimation and refreezing are two important processes for glaciers in this region but they are not included in empirical degree-day models (Litt and others, 2019). Moreover, refreezing is always parameterised in a more or less sophisticated way in physical snowpack models, but field experiments are crucially missing to evaluate the accuracy of such parameterisations.

6.2 Mass-balance sensitivity to different meteorological variables and comparison with other studies in HKH

6.2.1 Mass-balance sensitivity of Mera Glacier to meteorological variables

Estimating the sensitivity of glacier mass balance to a change in temperature and/or precipitation is a classical problem in glaciology (e.g. Ohmura and others, 1992). Different approaches have been implemented in HKH to assess the glacier mass-balance sensitivity at different scales (e.g. Sakai and Fujita, 2017; Wang and others, 2019; Arndt and others, 2021; Gurung and others, 2022). However, the classical approach (perturbing temperature and precipitation by certain values or percentages), which involves perturbing individual meteorological variables while keeping others unchanged, has been criticised due to the interconnectedness of these variables (Nicholson and others, 2013; Prinz and others, 2016). For example, changing air temperature directly impacts the water vapour ratio of the atmosphere, which can ultimately affect the turbulent heat fluxes.

For Mera Glacier, we find strong correlations among the meteorological variables, which makes difficult to decipher their individual effects on mass balance. Notably, the negative correlation between temperature and precipitation anomalies (r = −0.56; slope of the linear relationship: −0.16 m w.e. °C⁻¹) leads to larger mass-balance sensitivities to temperature when this relationship is preserved (Table 7). The mass-balance sensitivity to SW_in is unexpectedly limited (−0.06 m w.e. a⁻¹ (W m⁻²)⁻¹; r = −0.22) due to the role of albedo. Indeed, the correlation between albedo anomalies and mass-balance anomalies is very high (r = 0.97) showing that the surface state, that depends primarily on the amount of snowfalls, is more important than the actual incoming shortwave radiative flux. Despite its direct contribution to the SEB, LW_in correlates positively with mass balance due to its strong correlation with precipitation. One limitation of our approach is that we can virtually find sensitivities of the mass balance to any input variable, as long as a correlation exists. For instance, the strong sensitivity to the relative humidity shown in Figure 8 should not be interpreted as an expected change in mass balance due to a change in relative humidity. Instead it shows that the relative humidity is closely tied to the other meteorological variables, and that meteorological conditions that favour high relative humidity also favour positive mass balances.

Another interesting feature is the asymmetry towards negative values in the sensitivities to temperature with the classical method (Table 7). While continental and sub-continental glaciers typically exhibit less sensitivity to negative temperature changes (e.g. Wang and others, 2019; Arndt and Schneider, 2023), maritime glaciers tend to have symmetrical sensitivity, especially when estimated with degree-day models (Wang and others, 2019). However, with the synthetic scenario approach, we find the opposite asymmetry for temperature sensitivity (Table 7), likely due to the correlations between precipitation and temperature. This leads to higher negative mass balances than the simple +1°C perturbation (Fig. 8), indicative of the maritime climate setting for Mera Glacier. Regarding precipitation, its sensitivity is asymmetric towards negative values with both methods (Table 7). It is clearly due to albedo feedback effects in the model, which are poorly represented in degree-day models that predict a symmetric sensitivity to precipitation (e.g. Wang and others, 2019).

6.2.2 Mass-balance sensitivity comparison with other studies in HKH

The impact of temperature and precipitation changes on mass balance depends on the climate of the region, but it can be attenuated or exacerbated depending on the glacier's morphology and topography (Brun and others, 2023). Contrary to glaciers located in arid cold climates less sensitive to temperature changes (Ohmura and others, 1992), those affected by the Indian summer monsoon, such as Mera Glacier, are sensitive to both temperature and precipitation (Fujita, 2008; Johnson and Rupper, 2020; Arndt and others, 2021). With higher temperature, first less precipitation falls as snow and in turn accumulation is reduced, and second and more important, more shortwave radiation is absorbed through lower albedo leading to enhanced melt (Fujita, 2008). Still, the mass-balance sensitivity to temperature and precipitation varies among different glaciers. The glacier-wide sensitivity of Mera Glacier to changes in temperature and precipitation is of the same order of magnitude as other glaciers in HKH, even though it is noteworthy to mention that these glacier sensitivities have been assessed through the classical method and in turn are not directly comparable (Fig. 10, Table S3).

Figure 10. Location of glaciers where studies of mass-balance sensitivities have been conducted in the Himalaya and Tibetan Plateau regions. Each panel shows the mass-balance sensitivity to temperature and precipitation of each glacier, with the associated reference. Table S3 lists all these glaciers, and provides additional information.

For instance, on Mera Glacier, with the synthetic scenarios approach, a ±1°C temperature perturbation has a greater impact than that of a ±20% precipitation change, which is mostly the case for glaciers in Figure 10, especially those located in Nepal. However, this pattern differs for glaciers in the Indian western Himalaya, which mostly exhibit lower sensitivities, sometimes higher for precipitation than for temperature like in Sutri Dhaka Glacier (Oulkar and others, 2022). Surprisingly, Arndt and Schneider (2023) find extreme sensitivities of glacier-wide mass balances to warming or to increase in precipitation for glaciers in the Central Himalaya (Yala and Halji glaciers) or Nyainqentanglha Range (Zhadang Glacier) compared to what we observe on Mera Glacier, although all these glaciers are subjected to rather humid monsoon dominated conditions.

Our approach with the synthetic scenarios does not allow to investigate whether the sensitivity changes linearly. It is well established that the sensitivity to temperature is non-linear and much higher for larger temperature changes (e.g. Arndt and Schneider, 2023). As we rely only on existing observations, we cannot assess what would be the other meteorological variables in a +2 or +3°C climate setting. Directions to overcome this issue could be to investigate the links between glaciers’ response to synoptic variables (e.g. Mölg and others, 2012; Zhu and others, 2022), to investigate specific monsoon characteristics and their impacts on the mass balance (e.g. Shaw and others, 2022), or to force glacier mass-balance models with downscaled global circulation model outputs that preserve the physical relationships between variables (e.g. Bonekamp and others, 2019; Clauzel and others, 2023). Furthermore, the size of the glacier plays a role on its mass-balance sensitivity. Glaciers with a higher accumulation area, such as Sutri Dhaka, Trambau and Mera glaciers, exhibit a lower sensitivity than glaciers whose accumulation zone is strongly reduced, such as Zhadang and Halji glaciers (Zhu and others, 2018; Sunako and others, 2019; Arndt and others, 2021; Srivastava and Azam, 2022).

6.3 Limitations of our approach

Simulating the distributed surface energy and mass balance of a glacier presents numerous challenges and limitations. One striking example is the relatively large difference between the simulated and observed glacier-wide mass balances for 2019/20, where COSIPY simulated mass balance is 0.22 m w.e. larger than the observed one. While we do not have a definitive explanation for such a discrepancy, we can list a number of sources of errors in our approach. Uncertainties arise primarily from (1) the model's process representation, (2) in situ data, (3) their spatial distribution over the glacier area and (4) the model's initial conditions.

One crucial process in snowpack modelling is the decay of snow albedo over time. COSIPY implements the Oerlemans parameterisation (Oerlemans and Knap, 1998), which has known limitations in certain climate contexts (e.g. Wang and others, 2022; Voordendag and others, 2021). Typically, albedo parameters are fixed or adopted from previous studies using COSIPY or any other energy-balance model. However, here, we optimised these parameters at AWS-L before distributing the meteorological forcings. The snow ageing and depth scaling factors used in this study fall within the range of commonly used factors (Wang and others, 2019; Sauter and others, 2020; Arndt and others, 2021; Sherpa and others, 2023), but differ slightly from those used in other studies (Sauter and others, 2020; Arndt and others, 2021; Potocki and others, 2022). Although the model can accurately predict albedo at a point scale, it is sometimes misrepresented (Fig. S8). The lack of robust parameterisation and uncertainties surrounding the amount of snow deposited on the glacier surface and its redistribution by the wind contribute to this issue.

Additionally, refreezing is another poorly constrained process. In the HKH region, refreezing has primarily been assessed using models (Steiner and others, 2018; Kirkham and others, 2019; Saloranta and others, 2019; Veldhuijsen and others, 2022), rather than snowpack temperature and density measurements, as is done in the seasonal snowpack (e.g. Pfeffer and Humphrey, 1996). Specific experiments should be conducted to evaluate the effects of refreezing on glaciers, particularly to determine whether meltwater percolates below the previous year's horizon and contributes to internal accumulation. Additionally, COSIPY lacks certain processes, such as wind erosion and wind-driven snow densification. These processes can be very important, particularly during the post-monsoon and winter seasons due to the strong winds at high elevations (Brun and others 2023; Litt and others, 2019; Sherpa and others, 2023). During November field campaigns, wind erosion features such as sastrugi are frequently observed in the accumulation area of Mera Glacier.

Distributing meteorological data over a rough terrain is one of the most challenging task. The spatial distribution of the meteorological data based on a single vertical gradient throughout the year is somehow questionable (Section 3.2). Additionally, applying a vertical gradient to distribute meteorological variables weakens the physical links between them, as already discussed in the sensitivity analysis. For the study period and the range of elevation of Mera Glacier, the vertical gradients for temperature, relative humidity and LW_in exhibit high variability (Fig. S3). To derive the gradients, we selected the year with the minimum gaps at AWS-H despite a large portion of the data at AWS-L being reconstructed at that time. Additionally, the T/RH sensor was not artificially ventilated at AWS-H, which introduces some measurement uncertainty. Therefore, we tested a range of gradients to distribute T [−6.5 to −4.2°C km⁻¹] and RH [−25 to 0% km⁻¹]. In alpine environments, LW_in provides large amounts of melt energy, especially in the ablation area (through valley-side walls) and can dominate the energy balance of snow or glacier surfaces. LW_in is highly sensitive to surface melt when the atmosphere is saturated, particularly during the monsoon (Sicart and others, 2010). On Mera Glacier, the daily LW_in gradient as a function of elevation varies greatly from day to day and season to season, with a minimum of −41 W m⁻² km⁻¹ observed during the monsoon. Various LW_in gradients have been tested within the [−40 to 0 W m⁻² km⁻¹] range and optimised to −25 W m⁻² km⁻¹ (Table S2). In conclusion, there is no ideal method for spatially distributing meteorological variables on a complex glacier surface that extends over a large altitudinal range. On Mera Glacier, two on-glacier AWSs enable us to provide reasonable vertical gradients of temperature, relative humidity and longwave incoming radiation. These gradients are highly variable in time and likely in space as well. To maintain simplicity in our modelling approach, we prefer to use a single gradient for these variables instead of using temporally or spatially variable gradients. The use of variable gradients would require a denser observation network than the two AWSs we currently have. Additionally, using different gradients would alter the set of optimised parameters without necessarily affecting the final results.

The depletion of precipitation as a function of elevation in the Upper Khumbu region is still a matter of debate due to the lack of reliable data at high elevations, difficulties in correcting the undercatch of snow, and comparing precipitation records obtained with different devices (Salerno and others, 1994; Perry and others, 2020). On the Mera catchment, there is only one all-weather rain gauge located at 4888 m a.s.l., just below the glacier. The precipitation recorded at this station is considered constant across the glacier surface. The distribution of precipitation is clearly more complex due to the rough topography (Immerzeel and others, 2014) and snow redistribution by wind. We prefer not to apply any vertical precipitation gradient, as it would complicate the modelling and increase equifinality issues, as with other meteorological variables. Similarly, applying a vertical gradient of wind speed based on records at AWS-L and AWS-H is not recommended due to the site-specific nature of the records (Shea and others, 2015a). Therefore, like precipitation, wind speed is assumed to be constant and equal to the wind velocity at AWS-L across the glacier. Given that wind speed mainly affects turbulent fluxes, which are of secondary importance compared to radiative fluxes, any assumptions made only impact the total sublimation and its spatial distribution.

To initialise the model, the snow depth on each gridcell at the start of the simulation, specifically on 1 November, must be known. As field trips typically occur around mid-November each year, and November is a relatively dry month with minimal melting, snow depth initialisation is based on direct observations performed ~15 d later than the initialisation date. However, this method is not entirely error-free as not all gridcells are surveyed and precipitation or wind drift may have occurred between 1 November and the survey date. Such error may have large impacts over the entire simulation period when a surface is initially recognised snow free although it is not, or vice versa.

In COSIPY, between 10 and 20% of SW_net penetrates below the surface and is partly reflected at different depths within the snowpack or the ice (van den Broeke and Bintanja, 1995). This amount of shortwave radiation is not accounted for in the outtake term Q _out of Table 5, which explains why Q _in and Q _out do not exactly compensate each other. Therefore, the energy balance is not perfectly closed in the COSIPY model (e.g. Arndt and others, 2021).

7. Conclusion

The COSIPY energy and mass-balance model was applied to Mera Glacier in the Central Himalaya, Nepal. In situ meteorological datasets were used, recorded both on and off the glacier at different elevations ranging from 4888 m a.s.l. (for precipitation) to 5770 m a.s.l. The data were collected from 1 November 2016 to 31 October 2020 at an hourly time step. The model parameters were optimised at the point scale using data from AWS-L at 5360 m a.s.l. over 2018/19. A multi-objective optimisation was employed, and the albedo ageing and snow depth factors were selected as the most sensitive parameters and then calibrated. The model was validated both at the point scale with multiple AWS measurements (albedo and surface temperature recorded at 5360 and 5770 m a.s.l.) and at the glacier scale with annual point mass-balance measurements obtained at various elevations. The validation indicates that the model effectively simulates the annual glacier-wide mass balance of Mera Glacier, with a simulated mean value of −0.66 m w.e. a⁻¹ for 2016–20 compared to the observed value of −0.74 m w.e. a⁻¹.

The SEB over Mera Glacier is dominated by radiative fluxes, which are responsible for almost all the energy available during the monsoon, the main melting season. From the pre-monsoon to the monsoon, with the increasing cloudiness, incoming shortwave radiation gradually decreases in intensity in favour of incoming longwave radiation thus maintaining a large amount of energy available for melt. Turbulent fluxes are only significant outside the monsoon. The sensible heat is an energy source at the surface whereas the latent heat flux is always negative. Sublimation is therefore an important ablation process, especially during the windy months of the post-monsoon and winter. Annually and at glacier scale, refreezing is a crucial process, because on average 44% of meltwater refreezes, with a large positive gradient with altitude.

To investigate the sensitivity of glacier mass balance to changes in temperature and precipitation, we generated 180 different scenarios by shuffling our 4 year dataset. We aggregated warm, cold, dry or wet months alternatively, depending on the seasons. These scenarios allowed us to explore a wide range of conditions, from very dry and warm to very cold and wet. As a result, the glacier-wide mass balances of Mera Glacier ranged from −1.76 to +0.54 m w.e. a⁻¹. The mass-balance sensitivity to meteorological variables can be quantified from these synthetic scenarios. A temperature change of +1°C (−1°C) results in a change of −0.75 ± 0.17 (+0.93 ± 0.18) m w.e. in glacier-wide mass balance, while a precipitation change of +20% (−20%) results in a change of +0.52 ± 0.10 (−0.60 ± 0.11) m w.e. in mass balance. Compared to the classical approach, the sensitivity of the mass balance is more pronounced with temperature, but not significantly different with precipitation. Similar to other glaciers with summer accumulation, Mera Glacier is highly sensitive to both temperature and precipitation.

To evaluate the mass-balance sensitivity to any meteorological variables, it is advantageous to generate scenarios based on real in situ data. This not only helps to quantify these sensitivities more accurately but also to explore the inter-relationships between variables. Our study demonstrates, for instance, that temperature has a negative correlation with precipitation. Therefore, classical sensitivity approaches that alter temperature and precipitation independently are likely to be biased. It is worth noting that long-term high-quality datasets are necessary to apply such a synthetic method approach. We are lucky enough to have a long-term dataset on Mera Glacier, but we encourage to maintain and develop similar observational networks on other glaciers in HKH, in order to compare glaciers and to assess whether sensitivities obtain locally on a glacier can be extrapolated regionally. Currently, we cannot be certain that this sensitivity is not specific to Mera Glacier.

Like any modelling, our approach has limitations inherent to the model used, the quality of input data, their spatial distribution, and the choice of initial conditions. These limitations are difficult to quantify, but our method allows us to provide an accuracy range for the results based on a 99% confidence interval. A potential next step in this study would be to conduct an uncertainty analysis to assess the weight of all potential errors related to the model and data, as well as to evaluate the equifinality of the results. However, this is beyond the scope of the present study. Nonetheless, such an analysis would be valuable, as many modelling approaches encounter similar issues.