1. Introduction
The pace of climate warming in High Mountain Asia (HMA) is accelerating (Pepin and others, Reference Pepin2015) and precipitation in these regions exhibits significant heterogeneity and remains insufficiently comprehended (Lutz and others, Reference Lutz, Immerzeel, Shrestha and Bierkens2014). 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, Reference Oerlemans and Reichert2000). 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, Reference Hugonnet2021), 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, Reference Oerlemans and Reichert2000; Marzeion and others, Reference Marzeion, Jarosch and Hofer2012). Following this approach, Sakai and Fujita (Reference Sakai and Fujita2017) 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, Reference Ohmura, Kasser and Funk1992) 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, Reference Fujita2008; Azam and others, Reference Azam2014; Fugger and others, Reference Fugger2022). 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, Reference Oerlemans and Knap1998; Favier and others, Reference Favier, Wagnon, Chazarin, Maisincho and Coudrain2004). Specifically, in Hindu Kush Himalaya (HKH), a number of studies investigated the SEB of glaciers in different climate contexts (Mölg and others, Reference Mölg, Maussion, Yang and Scherer2012; Huintjes and others, Reference Huintjes, Neckel, Hochschild and Schneider2015a, Reference Huintjes2015b; Zhu and others, Reference Zhu2015, Reference Zhu2018, Reference Zhu2021; Fugger and others, Reference Fugger2022; Arndt and Schneider, Reference Arndt and Schneider2023). 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, Reference Arndt and Schneider2023). 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, Reference Wang, Liu, Shangguan, Radić and Zhang2019). This is due to the highly non-linear response of glacier surface mass balance to the albedo effect (e.g. Arndt and Schneider, Reference Arndt and Schneider2023).
However, in HKH most of the SEB studies have two main limitations: either they were conducted at point scale (Kayastha and others, Reference Kayastha, Ohata and Ageta1999; Azam and others, Reference Azam2014; Acharya and Kayastha, Reference Acharya and Kayastha2019; Litt and others, Reference Litt2019; Mandal and others, Reference Mandal2022), or they used meteorological data from reanalysis products (Arndt and Schneider, Reference Arndt and Schneider2023). 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, Reference Oerlemans1999). 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, Reference Huintjes, Neckel, Hochschild and Schneider2015a, Reference Huintjes2015b; Arndt and others, Reference Arndt, Scherer and Schneider2021; Oulkar and others, Reference Oulkar2022; Srivastava and Azam, Reference Srivastava and Azam2022). Meteorological variables obtained from reanalysis can be heavily biased, especially if they are not downscaled with local measurements (e.g. Hamm and others, Reference Hamm2020; Khadka and others, Reference Khadka2022). 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, Reference Rounce, Hock and Shean2020).
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, Reference Sauter, Arndt and Schneider2020) 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, Reference Wagnon2021; Khadka and others, Reference Khadka2022). 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 km2 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 km2, 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, Reference Wagnon2021). 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, Reference Wagnon2013; Thakuri and others, Reference Thakuri2014; Shea and others, Reference Shea2015a). 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, Reference Ageta and Higuchi1984). 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, Reference Khadka2022). During this season, warm air masses flow from the Bay of Bengal and bring moisture and precipitation in the Himalaya (Perry and others, Reference Perry2020). 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, Reference Shea2015a). 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, Reference Khadka2022).
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 (LWin and LWout, respectively) and shortwave (SWin and SWout, 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; SWin, incoming shortwave radiation; SWout, outgoing shortwave radiation; LWin, incoming longwave radiation; LWout, 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 (SWin), (e) incoming longwave radiation (LWin) 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 (Reference Førland1996) and Lejeune and others (Reference Lejeune2007) 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, Reference Khadka2022).
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 (Reference Sauter, Arndt and Schneider2020), already tested and applied on Himalayan glaciers (e.g. Arndt and Schneider, Reference Arndt and Schneider2023). First the fraction of diffuse radiation (F diff) is calculated based on Wohlfahrt and others (Reference Wohlfahrt2016):
where p1 = 0.1001, p2 = 4.7930, p3 = 9.4758, p4 = 0.2465 are parameters from Wohlfahrt and others (Reference Wohlfahrt2016) 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 = SWin(1 − F diff)) and diffuse (R d = SWin × F diff) radiation. Then, the corrected solar radiation (R c) is calculated on each grid based on Ham (Reference Ham2005):
where c f is the correction factor calculated based on the azimuth and the slope of each grid following Ham (Reference Ham2005).
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, Reference Wagnon2021). Over the period 2007–23, the mean corrected glacier-wide mass balance is equal to −0.42 ± 0.23 m w.e. a−1, 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−1, 2017–18 being the most negative year (−0.92 ± 0.16 m w.e. a−1) and 2019–20 being the least negative (−0.49 ± 0.22 m w.e. a−1) (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−3 and measured snow densities were used to compute point mass balances (370 kg m−3 for Naulek, and 380 to 430 kg m−3 for the accumulation area). The error range for point mass balances is the std dev. of all measurements
*Updated from Wagnon and others (Reference Wagnon2021).
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−3 and an ice density of 900 kg m−3. 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−3).
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, Reference Sauter, Arndt and Schneider2020). 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, Reference Sauter, Arndt and Schneider2020; Arndt and others, Reference Arndt, Scherer and Schneider2021; Blau and others, Reference Blau, Turton, Sauter and Mölg2021). The COSIPY model calculates the energy available for melt (Q M) for each time step and is expressed as
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, Reference Sauter, Arndt and Schneider2020) and R: rain in W m−2. 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, Reference Sauter, Arndt and Schneider2020). 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
where ρ air is air density (in kg m−3); C P is the specific heat of air at constant pressure (in J kg−1 K−1) and L V is the latent heat of sublimation/vapourisation (in J kg−1). 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 (Reference Mölg, Maussion, Yang and Scherer2012) and further calibrated (see Section 4.2), q is the specific humidity of air (in g kg−1), 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 (Reference Hantel, Ehrendorfer and Haslinger2000) (Fig. S6):
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 0 is the melting point (0°C), T 00 is the centre for snow transfer function (in °C) and s 0 is the spread snow transfer function. The ageing/decay of the snowpack's albedo is then based on Oerlemans and Knap (Reference Oerlemans and Knap1998), 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 (Reference Essery, Morin, Lejeune and Ménard2013). 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 km2) resolution (corresponding to a total glacier-covered area of 5.01 km2), resampled from Copernicus GLO-30 DEM (https://doi.org/10.5270/ESA-c5d3d65) using the glacier outlines in 2012 (5.10 km2) from Wagnon and others (Reference Wagnon2021). 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, Reference Wagnon2013). 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, Reference Zolles, Maussion, Peter Galos, Gurgiser and Nicholson2019). 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 LWout using the Stefan–Boltzmann equation with a surface emissivity of 0.99 (Blau and others, Reference Blau, Turton, Sauter and Mölg2021), 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 (Reference Mölg, Maussion, Yang and Scherer2012). 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, Reference Zolles, Maussion, Peter Galos, Gurgiser and Nicholson2019). 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, Reference Mölg, Maussion, Yang and Scherer2012; Zolles and others, Reference Zolles, Maussion, Peter Galos, Gurgiser and Nicholson2019). 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 (Reference Mölg, Maussion, Yang and Scherer2012) 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, Reference Yapo, Gupta and Sorooshian1998; Rye and others, Reference Rye, Willis, Arnold and Kohler2012; Zolles and others, Reference Zolles, Maussion, Peter Galos, Gurgiser and Nicholson2019). The multi-objective approach can be expressed as
where f 1(θ), f 2(θ), …, 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, Reference Rye, Willis, Arnold and Kohler2012):
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, Reference Pareto1971). 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:
where f 1 is the coefficient of determination (r 2) between the observed (x i) and modelled (y i) values of albedo and surface temperature (T s), f 2 is the mean absolute error (MAE) between the observed (x i) and modelled (y i) values of both variables and f 3 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, Reference Soon and Madsen2005). 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, Reference Kayastha, Ohata and Ageta1999; Mölg and others, Reference Mölg, Maussion, Yang and Scherer2012; Sunako and others, Reference Sunako, Fujita, Sakai and Kayastha2019; Arndt and others, Reference Arndt, Scherer and Schneider2021; Gurung and others, Reference Gurung2022; Srivastava and Azam, Reference Srivastava and Azam2022; Arndt and Schneider, Reference Arndt and Schneider2023). 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, Reference Prinz, Nicholson, Mölg, Gurgiser and Kaser2016; Clauzel and others, Reference Clauzel2023). To overcome this issue, a number of methods were developed to perturb meteorological data while preserving the link between variables (Sicart and others, Reference Sicart, Hock, Ribstein and Chazarin2010; Prinz and others, Reference Prinz, Nicholson, Mölg, Gurgiser and Kaser2016; Autin and others, Reference Autin, Sicart, Rabatel, Soruco and Hock2022). 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 ABX, 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_Wamon 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_Wamon 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_2020melt 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 2 and MAE of albedo being 0.06–0.54 and 0.11–0.25, respectively (Fig. 3a). Similarly for surface temperature, r 2 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 2 and MAE of albedo being 0.31–0.54 and 0.11–0.15 (Fig. 3a), and the range of r 2 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 2 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 2 from albedo, (d) 1 − r 2 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 2 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 2 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 2 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, Reference Oerlemans and Knap1998) 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 2 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 SWin 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 2 (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 2 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 (Reference Wagnon2021). Over the 4 year period, the mean observed glacier-wide in situ mass balance is −0.74 ± 0.18 m w.e. a−1 (Table 2) and the simulated mass balance is −0.66 m w.e. a−1 with a std dev. of ±0.26 m w.e. a−1. 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, Reference Wagnon2021) 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 (Reference Wagnon2021).
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. SWnet 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, SWin is controlled by the position of the sun responsible for the potential SWin and also by cloudiness, explaining why it decreases from 274 W m−2 in the pre-monsoon to 195 W m−2 during the monsoon (Table 5). Similarly, SWnet decreases from 83 W m−2 during the pre-monsoon to 61 W m−2 during the monsoon because of SWin 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 LWin from the pre-monsoon (217 W m−2) to the monsoon (296 W m−2), the only season when LWin nearly counterbalances the LWout.
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). SWnet, net shortwave radiation; LWnet, 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−2) 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−2) and the monsoon (491 W m−2) (Table 5). However, the net all-wave radiation, calculated as the sum of SWnet and LWnet, is 32 W m−2 during the pre-monsoon and 20 W m−2 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−2 each). During the monsoon, LWnet 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−2. 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 SWnet and Q M are higher at AWS-L (93 and −34 W m−2, respectively) than at AWS-H (61 and −15 W m−2, respectively). LWnet, Q L and Q S exhibit similar annual means throughout the year at both sites. At AWS-L, SWnet remains similar during the pre-monsoon (87 W m−2) and the monsoon (88 W m−2) because the decrease in SWin (287 W m−2 in the pre-monsoon and 201 W m−2 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 SWnet in the monsoon (43 W m−2) compared to the pre-monsoon (76 W m−2). The variation of LWnet at AWS-L and AWS-H is rather small as the difference between LWin and LWout 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−2, 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−2) 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 SWin anomaly correlates significantly and negatively with the relative humidity, the precipitation and the LWin anomalies. The LWin anomaly correlates significantly and positively with the relative humidity and the precipitation anomalies, and negatively with the SWin 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 SWin (Fig. 8). The correlations between mass-balance anomalies and those of LWin, atmospheric pressure or wind speed are moderate but significant, positive in case of LWin, 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, Reference Wagnon2021). 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 SWnet anomaly compared to the reference year (+2 to +17 W m−2), 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. SWnet, net shortwave radiation; LWnet, 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_Wemon, D_Wemon, Wa_Wewin, D_Cwin). 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 SWnet (−8 to +6 W m−2) and LWnet (−4 to +3 W m−2) 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 SWnet (−19 to −4 W m−2). 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 LWnet is similar in all classical scenarios. However, SWnet is 13 W m−2 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, Reference Mandal2022; Oulkar and others, Reference Oulkar2022; Srivastava and Azam, Reference Srivastava and Azam2022). 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, Reference Sakai and Fujita2017). 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, Reference Shaw2022).
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 (Reference Azam2014) are not reported in this table. All the fluxes are in W m−2.
*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, Reference Azam2014, 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, Reference Azam2014), 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, Reference Mölg, Maussion, Yang and Scherer2012; Sun and others, Reference Sun2014; Zhu and others, Reference Zhu2018). 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, Reference Gurung2022; Srivastava and Azam, Reference Srivastava and Azam2022) or lower than that of Sutri Dhaka Glacier in the western Himalaya (55%; Oulkar and others, Reference Oulkar2022). However, this value on Mera Glacier is under-estimated because wind speed strongly increases with elevation (Khadka and others, Reference Khadka2022), 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, Reference Stigter2018; Bonekamp and others, Reference Bonekamp, de Kok, Collier and Immerzeel2019; Kirkham and others, Reference Kirkham2019; Veldhuijsen and others, Reference Veldhuijsen2022; Arndt and Schneider, Reference Arndt and Schneider2023). 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, Reference Litt2019). 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, Reference Ohmura, Kasser and Funk1992). Different approaches have been implemented in HKH to assess the glacier mass-balance sensitivity at different scales (e.g. Sakai and Fujita, Reference Sakai and Fujita2017; Wang and others, Reference Wang, Liu, Shangguan, Radić and Zhang2019; Arndt and others, Reference Arndt, Scherer and Schneider2021; Gurung and others, Reference Gurung2022). 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, Reference Nicholson, Prinz, Mölg and Kaser2013; Prinz and others, Reference Prinz, Nicholson, Mölg, Gurgiser and Kaser2016). 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−1) leads to larger mass-balance sensitivities to temperature when this relationship is preserved (Table 7). The mass-balance sensitivity to SWin is unexpectedly limited (−0.06 m w.e. a−1 (W m−2)−1; 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, LWin 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, Reference Wang, Liu, Shangguan, Radić and Zhang2019; Arndt and Schneider, Reference Arndt and Schneider2023), maritime glaciers tend to have symmetrical sensitivity, especially when estimated with degree-day models (Wang and others, Reference Wang, Liu, Shangguan, Radić and Zhang2019). 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, Reference Wang, Liu, Shangguan, Radić and Zhang2019).
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, Reference Brun2023). Contrary to glaciers located in arid cold climates less sensitive to temperature changes (Ohmura and others, Reference Ohmura, Kasser and Funk1992), those affected by the Indian summer monsoon, such as Mera Glacier, are sensitive to both temperature and precipitation (Fujita, Reference Fujita2008; Johnson and Rupper, Reference Johnson and Rupper2020; Arndt and others, Reference Arndt, Scherer and Schneider2021). 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, Reference Fujita2008). 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, Reference Oulkar2022). Surprisingly, Arndt and Schneider (Reference Arndt and Schneider2023) 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, Reference Arndt and Schneider2023). 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, Reference Mölg, Maussion, Yang and Scherer2012; Zhu and others, Reference Zhu2022), to investigate specific monsoon characteristics and their impacts on the mass balance (e.g. Shaw and others, Reference Shaw2022), 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, Reference Bonekamp, de Kok, Collier and Immerzeel2019; Clauzel and others, Reference Clauzel2023). 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, Reference Zhu2018; Sunako and others, Reference Sunako, Fujita, Sakai and Kayastha2019; Arndt and others, Reference Arndt, Scherer and Schneider2021; Srivastava and Azam, Reference Srivastava and Azam2022).
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, Reference Oerlemans and Knap1998), which has known limitations in certain climate contexts (e.g. Wang and others, Reference Wang2022; Voordendag and others, Reference Voordendag, Réveillet, Macdonell and Lhermitte2021). 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, Reference Wang, Liu, Shangguan, Radić and Zhang2019; Sauter and others, Reference Sauter, Arndt and Schneider2020; Arndt and others, Reference Arndt, Scherer and Schneider2021; Sherpa and others, Reference Sherpa2023), but differ slightly from those used in other studies (Sauter and others, Reference Sauter, Arndt and Schneider2020; Arndt and others, Reference Arndt, Scherer and Schneider2021; Potocki and others, Reference Potocki2022). 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, Reference Steiner2018; Kirkham and others, Reference Kirkham2019; Saloranta and others, Reference Saloranta2019; Veldhuijsen and others, Reference Veldhuijsen2022), rather than snowpack temperature and density measurements, as is done in the seasonal snowpack (e.g. Pfeffer and Humphrey, Reference Pfeffer and Humphrey1996). 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 Reference Brun2023; Litt and others, Reference Litt2019; Sherpa and others, Reference Sherpa2023). 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 LWin 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−1] and RH [−25 to 0% km−1]. In alpine environments, LWin 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. LWin is highly sensitive to surface melt when the atmosphere is saturated, particularly during the monsoon (Sicart and others, Reference Sicart, Hock, Ribstein and Chazarin2010). On Mera Glacier, the daily LWin gradient as a function of elevation varies greatly from day to day and season to season, with a minimum of −41 W m−2 km−1 observed during the monsoon. Various LWin gradients have been tested within the [−40 to 0 W m−2 km−1] range and optimised to −25 W m−2 km−1 (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, Reference Salerno1994; Perry and others, Reference Perry2020). 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, Reference Immerzeel, Petersen, Ragettli and Pellicciotti2014) 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, Reference Shea2015a). 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 SWnet penetrates below the surface and is partly reflected at different depths within the snowpack or the ice (van den Broeke and Bintanja, Reference van den Broeke and Bintanja1995). 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, Reference Arndt, Scherer and Schneider2021).
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−1 for 2016–20 compared to the observed value of −0.74 m w.e. a−1.
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−1. 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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2024.42.
Data
The data utilised in this study are accessible through the GLACIOCLIM database (https://glacioclim.osug.fr/Donnees-himalaya). The model outputs are available at https://doi.org/10.5281/zenodo.10053093.
Acknowledgements
This work has been supported by the French Service d'Observation GLACIOCLIM (part of IR OZCAR), the French National Research Agency (ANR; iFROG project, grant No. ANR-21-CE01-0012), and a grant from Labex OSUG@2020 (Investissements d'avenir ANR10-LABX56). This work would not have been possible without the International Joint Lab Water-Himal (principal investigators D. Shrestha and P. Wagnon) supported by IRD and all the efforts from people in the field: porters, students, and helpers who are greatly acknowledged here. A. Khadka is grateful to IRD and Perpetual Planet programme supported by the National Geographic Society and Rolex for providing financial support for his PhD. He had benefited from the support of the International Centre for Integrated Mountain Development (ICIMOD), which is funded in part by the governments of Afghanistan, Bangladesh, Bhutan, China, India, Myanmar, Nepal and Pakistan. We appreciate the Nepal Mountaineering Association, which provides free permits to conduct fieldwork on Mera peak. We also thank Anselm Arndt for the support in understanding and modifying the initial condition on the COSIPY. Also, we thank Simon Gascoin, Jakob Steiner, Argha Banejee, Arindan Mandal and one anonymous referee whose comments greatly helped to improve the work.
Author contributions
P.W. initiated and D.S. supported the monitoring programme on Mera Glacier. P.W., F.B. and A.K. jointly designed the study. A.K. performed the analysis, under the supervision of F.B. and P.W., with the support of T.C.S. A.K. prepared all the figures and wrote the first draft of the manuscript. A.K., F.B. and P.W. jointly developed the discussion and interpretation of the results. All authors participated in some field campaigns, and contributed to the analysis and to the paper writing.
Open access
