Atmosphere Driven Mass-Balance Sensitivity of Halji Glacier, Himalayas

The COupled Snowpack and Ice surface energy and mass balance model in PYthon (COSIPY) was employed to investigate the relationship between the variability and sensitivity of the mass balance record of the Halji glacier, in the Himalayas, north-western Nepal, over a 40 year period since October 1981 to atmospheric drivers. COSIPY was forced with the atmospheric reanalysis dataset ERA5-Land that has been statistically downscaled to the location of an automatic weather station at the Halji glacier. Glacier mass balance simulations with air temperature and precipitation perturbations were executed and teleconnections investigated. For the mass-balance years 1982 to 2019, a mean annual glacier-wide climatic mass balance of −0.48 meters water equivalent per year (m w.e. a−1) with large interannual variability (standard deviation 0.71 m w.e. a−1) was simulated. This variability is dominated by temperature and precipitation patterns. The Halji glacier is mostly sensitive to summer temperature and monsoon-related precipitation perturbations, which is reflected in a strong correlation with albedo. According to the simulations, the climate sensitivity with respect to either positive or negative air temperature and precipitation changes is nonlinear: A mean temperature increase (decrease) of 1 K would result in a change of the glacier-wide climatic mass balance of −1.43 m w.e. a−1 (0.99m w.e. a−1) while a precipitation increase (decrease) of 10% would cause a change of 0.45m w.e. a−1 (−0.59m w.e. a−1). Out of 22 circulation and monsoon indexes, only the Webster and Yang Monsoon index and Polar/Eurasia index provide significant correlations with the glacier-wide climatic mass balance. Based on the strong dependency of the climatic mass balance from summer season conditions, we conclude that the snow–albedo feedback in summer is crucial for the Halji glacier. This finding is also reflected in the correlation of albedo with the Webster and Yang Monsoon index.

Kropáček et al. [38] determined an annual geodetic MB of (−0.40 ± 0.30) m w.e. a −1 between 2000 and 2013. Shean et al. ([7], dataset: [49]) calculated an annual geodetic MB of (−0.70 ± 0.09) m w.e. a −1 between 2000 and 2018. Both studies use digital elevation model (DEM) differencing but for different periods. As reported by Ye et al. [50], the glaciers in the Naimona'nyi region (∼25 km northwest of Halji glacier) have retreated at least since 1976. The only information about the Halji glacier development pre-millennial are glacier outlines derived from satellite data and the resulting glacier areas as displayed in Figure 1. The glacier area decreased from 3.1 km 2 in 1974 [51] to 2.3 km 2 in 2001 [45]. Between 2001 and 2010, the area was relatively stable with still 2.2 km 2 in 2010 [52], but decreased to 1.9 km 2 by 2018. The 2018 outline was derived from a Sentinel 2 [53] scene found with the Google Earth Engine Digitisation Tool developed by Lea [54]. According to the outline of the Randolph Glacier Inventory 6.0 (RGI6) [45], the glacier has two ice divides (see green lines Figure 1). All other outlines have been adjusted to these divides. 30

Data and Methods
In the following section, the applied physically-based surface energy and mass balance model COSIPY is described. Afterwards, the AWS and the meteorological forcing of COSIPY with the applied downscaling approaches are presented in Section 3.2. The setup for the simulations in this study is summarized in Section 3.3. The applied statistics, the seasonal sensitivity characteristic (SSC) and circulating and monsoon indexes are presented in Section 3.4. An overview of all introduced symbols and constants with their units can be found in Appendix A.

COSIPY
COSIPY is of medium complexity within the range of available surface energy balance (SEB) and MB models [42]. It is physically based and combines SEB processes with an adaptive (non-equidistant) subsurface scheme. The possibility of the identification of important atmospheric MB drivers, the easy implementation on High-Performance Computing Clusters (HPCCs) and the modular structure (maximum traceability) are some of the key points focused on during the development of the model. The open-source model is written in Python 3. The source code is freely available on a git repository (https://github.com/cryotools/cosipy, accessed on 13 March 2021). We use version COSIPY v1.3 (https://doi.org/10.5281/zenodo.3902191, accessed on 13 March 2021) in the distributed setup. The temporal and spatial resolution of the simulations is scalable, depending on the temporal resolution of forcing variables and the spatial resolution of the applied DEM. The SEB is defined as the sum of all energy fluxes at the surface: where Q M is available melt energy, Q SWin is incoming shortwave radiation, α is snow/ice albedo, Q LWin is incoming longwave radiation, Q LWout is outgoing longwave radiation, Q H is sensible heat flux, Q E is latent heat flux, Q G is glacier heat flux and Q R sensible heat flux of rain. Q SWin and Q LWin are input parameters of COSIPY whereby Q LWin can be parametrized (see Equations (14) and (15), [42]). Cloud cover fraction is needed as input in the latter case. The decay of α is calculated after Oerlemans and Knap [57] with the parameters presented and studied by Mölg et al. [58] in HMA. All other terms on the right-hand side of Equation (1) depend on the surface temperature T s . The resulting nonlinear equation is solved iteratively with an optimization algorithm. Due to physical constraints, T s ≤ melting point temperature T m must be fulfilled and therefore energy surplus results in Q M > 0 when T s equals T m . If T s < T m , there is no available energy at the surface, i.e., Q M = 0. All energy fluxes are of positive (negative) algebraic sign towards (directed away from) the surface and presented in W m −2 .
The SEB is linked through T s as an upper Neumann boundary condition [59] to the heat equation of the subsurface module. Furthermore, surface melt and rain percolates through the snowpack and can result in refreezing within the snow layers. A small part of net shortwave radiation Q SWnet = Q SWin (1 − α), penetrates the uppermost layers [60], warms these layers and result in subsurface melt, if the additional energy would otherwise result in inconsistent layer temperature T l >T m . Moreover, parametrizations for the densification of the snow [61,62] and the development of the snow/ice surface roughness [58] were implemented.
COSIPY is a one-dimensional point model that calculates in its distributed setup the MB at each point on the glacier and neglects exchange of lateral mass and energy transport. Processes at the base of the glacier are also neglected, which is why in the following, the simulated MB is referred to as the climatic mass balance b clim in accordance with Cogley et al. [63]. It is the sum of mass changes through accumulation and ablation and therefore links climate variability to glacier-mass gain or loss [64]. It is calculated as follows [63]: where surface accumulation is c s f c , surface ablation is a s f c , internal accumulation is c i , internal ablation is a i , surface mass balance is b s f c and internal mass balance is b i . Accumulated snowfall SF c and deposition (resublimation) of water vapor result in c s f c . COSIPY does not consider direct sublimation during snowfall and processes associated with snowdrift. Therefore, frozen precipitation equals SF c within COSIPY. Surface ablation consists of sublimation and surface melt resulting from Equation (1). Internal accumulation equals refreezing within the snowpack and a i equals subsurface melt through penetrating radiation. Mass balance components are positive (negative) if they are a mass gain (loss) in m w.e. In accordance with Cogley et al. [63], a capital B denotes the glacier-wide MB (sum of all point balances) with the same suffixes than b. The annual climatic mass balance is b clim,a . The cumulative climatic mass balance is b clim,cum . When we refer to the massbalance year (MB-year) of b clim,a and other variables, we speak of the hydrological year, starting after the approximated end of the preceding ablation season on 1 October and ending on 30 September of the following year, whereby the year of the included January is used as name of the MB-year in accordance with Cogley et al. [63]. For further information on applied parametrizations, physical principles and technical infrastructure of COSIPY, please read Sauter et al. [42].

Automatic Weather Station Measurements
An AWS [65] was installed in December 2013 by the Chair of Climatology (Technische Universität Berlin) and was replaced by a new AWS in 2018 (see Figure 1). The AWS is located adjacent to the Halji glacier at an elevation of 5359 m a. s. l. The data are continuously transmitted via Iridium satellite telemetry to a server of the Chair of Climatology at Technische Universität Berlin, Germany. Observations from 20 April 2018 until 2 November 2019 are used to statistically downscale atmospheric reanalysis data. Due to technical problems with total precipitation measurements, the available period starts in that case on 2 August 2018. The measured variables, measuring instruments and nominal accuracies are summarized in Table 1. Data of the first installed AWS from 1 December 2013 until 17 April 2018 were additionally used to evaluate downscaled surface pressure p s f c , air temperature at 2 m T 2 and relative humidity at 2 m RH 2 .

Downscaling ERA5-L
Required COSIPY forcing variables are derived from the atmospheric reanalyses ERA5-L produced by ECMWF within the Copernicus Climate Change Service. ERA5-L-a subversion of ERA5-provides ERA5 surface variables on a ∼9 km grid, instead of the original ∼31 km ERA5 grid. After producing ERA5, the Hydrology revised Tiled ECMWF Scheme for Surface Exchanges over Land (H-TESSEL) of the ERA5 integrated forecast system (IFS CY41R2) was re-executed with spatially interpolated (to a 9 km grid) ERA5 variables to produce ERA5-L [66]. Air temperature at 2 m T 2 and dewpoint temperature at 2 m T d,2 are the only COSIPY forcing variables that contain new information in comparison to ERA5. All other variables are also H-TESSEL forcing variables and therefore contain no new information [41,66,67]. We used hourly data from January 1981 until April 2020 from the grid cell, which represents the area between 30.25°N, 81.45°E and 30.35°N, 81.55°E. Halji glacier lies completely within this grid cell. ERA5-L data can be downloaded free of charge from the climate data store (https://cds.climate.copernicus.eu/cdsapp#!/dataset/ reanalysis-era5-land?tab=form, accessed on 12 December 2020).
The preprocessing steps to start COSIPY simulations can be divided into two steps. An overview of the applied approaches is provided in Table 2. First, we downscale the ERA5-L data to the location of the installed AWS at the glacier. Second, the downscaled forcing variables at the location of the AWS are interpolated with a DEM to the applied distributed fields on the glacier. The used ERA5-L grid cell has a model elevation of 5154 m a. s. l., which lies 205 m below the elevation of the AWS. For surface pressure p s f c , we apply the barometric formula for both the downscaling to the elevation of the AWS and the interpolation to the distributed fields [69]: where p 0 is atmospheric pressure at reference height, a thermal gradient, T 0 air temperature at reference height, h height difference, M average molar mass of air, g gravitational acceleration and R gas constant (see Appendix List of constants). We observed that the seasonal amplitude of ERA5-L's T 2 is greater than the measured T 2 , resulting in a warm bias of temperature in summer and a cold bias in winter temperature. Therefore, we applied a quantile-mapping approach [70] to downscale the data to the AWS. For generation of the distributed T 2 fields, we use a lapse rate of −0.6 K (100 m) −1 . The lapse rate is calculated from the long-term mean (1981 to 2019) of T 2 of 121 ERA5-L grid cells around Halji glacier (see Appendix A Figure A2).
Derived ERA5-L relative humidity at 2 m RH 2 is calculated from saturation water vapor of T 2 and T d,2 (see [71,72]). We use a constant value for all glacier grid points (GGPs) for the distributed RH 2 fields on the glacier. For the two radiation components Q SWin and Q LWin we do not use downscaling from raw ERA5-L to the location of the AWS. However, a radiation model after Wohlfahrt et al. [68] is used to calculate the distributed input fields of Q SWin on the glacier. The model takes slope, aspect, timestamp, latitude, longitude of the GGPs and Q SWin of ERA5-L as input and calculates Q SWin for each GGPs. Q LWin is taken constant in space for all GGPs without interpolation. To derive wind speed at 2 m U 2 , first wind speed at 10 m U 10 has to be calculated from u and v component of 10 m wind, which is described on the ERA5 website (https://confluence.ecmwf.int/pages/viewpage.action? pageId=133262398, last access: 12 December 2020). Afterwards, we use the logarithmic wind profile to calculate U 2 [73]: where z 0 is surface roughness. We use a value of 2.12 mm for z 0 , which is the mean between the upper and lower boundary condition (firn: 4 mm [74] and fresh snow: 0.24 mm [75]) in the surface roughness parametrizations of COSIPY. We use a constant value for all GGPs for distributed U 2 fields on the glacier. Figure 2 presents the comparison between measured and ERA5-L downscaled variables p s f c , T 2 and specific humidity at 2 m SH 2 . For comparison, we use SH 2 instead of RH 2 , the forcing variable required by COSIPY, to exclude differences that are solely based on differences in T 2 and p s f c . The course of the year of all three measured variables can be reproduced by ERA5-L, which is displayed in the left panels. This is supported by the scatterplots in the middle panels and the statistics of the applied linear least-square regression with all coefficients of determination r 2 ≥ 0.9 (p s f c : 0.99, T 2 : 0.90, SH 2 : 0.93), small root mean square errors (RMSE, p s f c : 0.5 hPa, T 2 : 2.3 K, SH 2 : 0.7 g kg −1 ), small mean bias errors (MBE, p s f c : 0.18 hPa, T 2 : <0.01 K, SH 2 : 0.03 g kg −1 ) and regression line slopes close to 1 (p s f c : 1.1, T 2 : 0.94, SH 2 : 1.18). All results are statistically significant with pvalue < 0.01. The cumulative distribution functions in the right panels reveal that the downscaling approach in case of p s f c and T 2 result in a substantially improved agreement with AWS measurements. A decrease in MBE from 14.3 to 0.2 hPa (RMSE: from 14.3 to 0.5 hPa) in case of p s f c and from 0.14 to <0.01 K (RMSE: from 3.4 to 2.3 K) in case of T 2 could be achieved. In case of SH 2 , almost no improvement (MBE: 0.05 to 0.03 g kg −1 ) can be observed. The reason might be that RH 2 mean of the measured values and the corresponding raw ERA5-L values differ only by 0.3 % (Measured: 77.1 %, ERA5-L: 76.8 %). We observed a distinct difference when comparing ERA5-L U 2 and measured U 2 . The hourly mean of measured U 2 is 5 times the ERA5 U 2 . Hourly values of U 2 are compared in Figure 3. Figure 3a shows the difference when comparing the data. It was not possible to apply a complex wind modeling approach within the scope of this study, also because of the lack of a highly resolved DEM. Therefore, we applied a scale factor of 5 to ERA5-L U 2 , which is presented in Figure 3b. The cumulative distribution function in Figure 3c reveals the improvement in the distribution of scaled ERA5-L U 2 data in comparison with measured data. This is supported by the decrease of the MBE from 3.74 to < 0.01 m s −1 .  In the period from August 2018 to November 2019, measured total precipitation TP was twice the amount of ERA5-L TP. We used this as a scale factor ( Figure 4). The scaling results in a decrease of MBE from 0.12 mm to 0.02 mm h −1 .

COSIPY Simulations
As static data, a glacier outline and DEM are required to start COSIPY simulations in the distributed setup. We used a time-invariant reference outline from RGI6 [45] (see Figure 1) for the simulation period October 1981-September 2019. Further, we used the 1 arcsecond global product [44] of the Shuttle Radar Topography Mission (SRTM) as DEM.
With the elevation information of all GGPs, we used the Geospatial Data Abstraction Library (GDAL) [76] to calculate all slopes and aspects of the GGPs needed as input for the radiation module of Wohlfahrt et al. [68]. Furthermore, we aggregated the DEM with GDAL to spatial resolutions between ∼30 and ∼1500 m. The applied spatial resolutions in arcseconds and ∼m and the number of resulting GGP is summarized in Table 3. We created the distributed COSIPY input fields and started simulations for the whole available ERA5-L period from January 1981 to April 2020 at an hourly resolution for all spatial resolutions. We then calculated annual deviations of B clim,a to the 30 m resolution simulation in order to identify the best trade-off between computational power and reasonable results. The deviations are presented in Figure 5.
All deviations until the resolution of ∼100 m are within a tolerable range (±0.02 m w.e. a −1 ). Therefore, for all further simulations in this study, we used a spatial resolution of ∼100 m with the resulting 248 GGPs representing the entire Halji glacier. For all distributed simulations, we used the HPCC of the Climate Geography lab of Humboldt-Universität zu Berlin, Germany. A simulation with ∼100 m resolution and 248 GGPs has a runtime of less than two hours for the whole ERA5-L period from January 1981 to April 2020.

Sensitivity Studies and Large Scale Teleconnections
Statistics: To investigate the influence of the forcing variables on the interannual variability of B clim,a , we calculate non-parametric Spearman's rank correlation coefficients r s [77] between B clim,a and the different forcing variables and the intermediate variables Q SWnet , α and net longwave radiation Q LWnet . We define two significance levels with a 95 % and 99 % confidence interval and access the significance with the two-sided p-value of the Spearman's rank correlation [77]. Furthermore, we start simulations with ±0.5, ±1.0, ±2.0 K T 2 and ±5, ±10, ±20 % TP perturbations and analyzed the resulting annual deviations of B clim,a to the original run.
Seasonal Sensitivity Characteristic (SSC): Oerlemans and Reichert [43] proposed calculating SSC to describe the dependency of B clim of a glacier to the local seasonal climate in a uniform and structured manner. The SSC of a glacier is a 2 × 12 matrix quantifying the B clim sensitivity to temperature and precipitation perturbations each in the 12 months of the year. For their calculation in the first step, the reference COSIPY simulation was adjusted so that the mean B clim,a is zero. Therefore, the mean forcing temperature is adjusted in the range of ±2 K. For the Halji glacier, we applied a T 2 offset of −0.37 K (see Appendix A Figure A3) to arrive at a zero B clim,a for the period October 1981-September 2019. In the second step, 24 COSIPY simulations were forced with monthly T 2 perturbations of ±0.5 K and 24 COSIPY simulations were forced with monthly TP perturbations of ±10 %. The resulting differences of B clim,a between the positive and negative perturbation for each month are the 12 temperature values and the 12 precipitation values of the SSC. These 24 values are displayed as a bar plot. For further information, including the equations of the concept, please see Section 2, especially Equations (3) and (4) in Oerlemans and Reichert [43] and Section 3b in Reichert et al. [78].
Indexes: Finally, we calculate r s between B clim,a and 22 common teleconnection indexes. The temporal resolution and coverage, short and long name, reference and download location of the indexes are summarized in Table A1 in the Appendix A. In doing so, we test the relationship between COSIPY-simulated Halji glacier B clim and atmospheric drivers.

Results
The   On annual basis the forcing variables T 2 , U 2 and TP have the highest correlation with all r s > 0.4. Whereby, T 2 and U 2 by themselves are highly significantly correlated with an annual r s of 0.69 (not displayed in Table 4, p-value < 0.01). No significant annual correlation could be found between T 2 and TP, and U 2 and TP. On a monthly basis, the October temperature at the beginning of the MB-year and the February, June and July temperature mainly determine the interannual B clim,a variability. Interestingly, the monthly sum of TP in September before the MB-year starts explains the interannual variability with r s = 0.43 (p-value < 0.01). The high statistical significance between RH 2 in June and B clim,a results to some degree from the correlation (r s = −0.37, p-value < 0.05, not displayed in Table 4) between T 2 and RH 2 . A statistically significant correlation of B clim and the forcing variables p s f c and Q SWin could be found neither on annual nor monthly time scales. However, a significant correlation (p-value < 0.05) with the monthly mean of Q LWin in February, May and June is present. Whereby, T 2 and Q LWin are by themselves highly significant correlated in February (r s = 0.47, p-value < 0.01), May (r s = 0.74, p-value < 0.01) and June (r s = 0.82, p-value < 0.01). The intermediate variables SF c , Q SWnet and α are superior in explaining interannual MB variability. The annual SF c is mainly a result of the TP and T 2 variability. The interannual variance of α is the dominating driver (r s = 0.8, p-value < 0.01) of the B clim long-term variability. The correlation of Q SWnet is in the same range because Q SWnet and α are the strongest coupled variables with an annual r s of 0.98 (not displayed in Table 4, p-value < 0.01).

Sensitivity to Climate Forcing Input Variables
The monthly deviations in percent (absolute values in Appendix A Figure A5) [43]. The B clim,a sensitivity to T 2 is dominated by the months June to August and partly September. The peak in summer is also visible in B clim,a sensitivity to monthly TP perturbations while it is not that pronounced as the T 2 sensitivity. The B clim,a sensitivity to November TP perturbations is the lowest. B clim,a is for both T 2 and TP most sensitive to July perturbations. The sum of the T 2 response to July and August is significantly greater than the response to all other months combined. In the case of TP, it is the sum of the response to July, August and September, which is greater than the sum of the response to all other months combined.

Index Correlations
As presented in Section 3.4, we analyzed the response of B clim,a to 22 circulation and monsoon indices. Analogous to the forcing variables, we calculated the non-parametric r s between B clim,a and the indices. With exception of the Webster and Yang Monsoon index (WYM) and the Polar/Eurasia index (POL), no significant correlations could be found with any of the indexes. Figure 10 displays B clim,a and the two indexes (annual values) that show significant correlations (p-value < 0.05) with B clim,a . The WYM is calculated from the zonal wind shear between 850 and 200 hPa [79] between 0°and 20°N and 40°and 110°E. It is a broad-scale monsoon index differentiating between strong and weak monsoon seasons [79]. . Seasonal sensitivity characteristic (SSC) after Oerlemans and Reichert [43]. Red bars are the dependence of annual glacier-wide climatic mass balance B clim,a on monthly temperature perturbations of 1 K and blue the dependence of B clim,a on monthly total precipitation perturbation of 10 %.
The POL links the Arctic polar vortex activity and mid-latitude circulation over the Asian continent [80,81] and is one of three modes describing anomalies that are derived from orthogonally rotated component analysis of the geopotential height at 700 hPa north of 15°N in the northern hemisphere [82,83]. The annual r s (see Table 5) is 0.48 (p-value < 0.01) with POL and −0.42 (p-value < 0.05) with WYM. On a monthly bases, the highest values could be found between June values of POL and B clim,a with r s = 0.56 (p-value < 0.01) and February with r s = 0.41 (p-value < 0.05). The highest correlation of the seasonal index WYM and B clim,a is the combined June, July and August value with a r s of −0.44 (p-value < 0.01).

Uncertainties
Model uncertainties can be separated in uncertainties stemming from forcing variables and uncertainties from the model itself.

Forcing Data Uncertainties
The statistical downscaling of the variables p s f c , T 2 and RH 2 to the elevation of the AWS could be realised to a high degree of precision. The AWS is not located on the glacier and can therefore not account for possible T 2 dampening through melting glacier surface (e.g., [84]). With the current setup, we cannot assess this effect. At an elevation of 5350 m a. s. l. the area around the AWS is temporally snow-covered even in summer and therefore T 2 can be damped as well. The lower summer AWS T 2 values compared to ERA5-L T 2 might be related to this effect. In the case of U 2 , the comparison between the ERA5-L and AWS values revealed a large difference. We scaled ERA5 U 2 to the location of the AWS. When it comes to local wind patterns, the location of the AWS is not representative for the whole glacier and affected by small scale features. Without a highly resolved DEM, we cannot achieve a better representation of U 2 . We have planned to implement a snowdrift scheme (e.g., [85]) to account for effects of drifting and blowing snow in future studies. The effects are redistribution of snow, changed albedo patterns, direct sublimation of drifting and blowing snow and modified vertical temperature and moisture profiles [85,86]. This will require highly resolved wind fields that can be used to improve the U 2 representation. With the current setup, it is not possible to access the uncertainties of Q LWin . We used the radiation model after Wohlfahrt et al. [68] for the generation of the distributed Q SWin fields, which corrects the input fields, mainly dependent on slope and aspect of the GGPs. Olson and Rupper [87] and Olson et al. [88] have shown that radiation modeling is most important for valley glaciers with steep surrounding terrain, for steep glaciers with high elevation differences and glaciers at high latitudes. The Halji glacier located at low latitude contradicts all these conditions as a relatively flat glacier with uninterrupted skyline. Therefore, the influence of the radiation modeling is assumed to be of minor importance.
As described in Section 3.2.2 TP from ERA5-L does not contain any new information in comparison with TP from ERA5. Taking TP from an atmospheric dataset with ∼31 km spatial resolution results in the largest uncertainty when it comes to the forcing variables especially in alpine regions with complex terrain and at high evaluations above 5000 m a. s. l. (e.g., [89][90][91][92]). Immerzeel et al. [91] used an inverse approach to find the amount of precipitation, which is needed to sustain the measured glacier MBs. Precipitation in extreme cases in the upper Indus basin is assumed to be up to ten times higher than previously assumed according to their results. We scaled ERA5-L TP to the cumulative sum of measured TP by the AWS. We applied a constant scale factor of two to the ERA5-L data. After applying the scale factor, the timing of ERA5-L TP follows measured TP with some differences (see Figure 4). In follow-up studies, spatially higher-resolved atmospheric datasets (e.g., the High Asia Refined analysis version 2-2 m [92,93]) should be used to reduce uncertainties stemming from TP.
The applied DEM and its resolution have an influence on the model results as well. The first reason is the influence on the radiation module [87,94]. Furthermore, the coarser the resolution, the less the resulting GGPs are representative for the different elevation ranges of the glacier. This is the reason for large and unsystematic deviations of resolutions with less than 30 GGPs (∼300 m) for Halji glacier compared to the resolution of 30 m (see Figure 5). With the presented setup, the deviations of B clim,a to the ∼30 m simulation are up to ∼100 m resolution in a tolerable range (see Section 3.3).

COSIPY Simulation Uncertainties
COSIPY is a point model with no lateral exchange of energy and mass. In its distributed setup, at each GGP all processes are calculated without any exchange with neighboring grid cells. Furthermore, over a period of 40 years, ice-dynamics play an important role. Nevertheless, due to its very low flow velocities (see Appendix A Figure A1), mass balance response of the Halji glacier can be analyzed concerning meteorological forcing only. An approximated mean flow velocity of 2.3 m a −1 ( [47], dataset: [48]) would result in 90 m ice movement over the whole period of 40 years. However, for further long-term integration of glacier change ice-dynamics will have to be integrated into the modeling framework. Snowdrift and direct sublimation of falling snow are further processes that are not resolved by the model. The model has been among others applied to the Zhadang glacier (see Section 5.1 in Sauter et al. [42]) and the Urumqi Glacier No. 1 (see Thiel et al. [95]), both located in HMA. The simulations in both cases revealed reasonable results. To assess the uncertainties stemming from all parameters and constants a full Monte Carlo simulation would have to be executed for the Halji glacier, which is not feasible in this study due to the enormous computational demand. A statistical error would affect single years, but not a long term cumulative value because the error would decrease according to the central limit theorem [96]. A systematic error would indeed affect the cumulative value of B clim . Due to the significant correlation between α and B clim,a and experiences from Mölg et al. [58], Sauter et al. [42] and Thiel et al. [95], the Oerlemans and Knap [57] α parametrization is crucial for model output. We used the parameters of Mölg et al. [58]. For a sensitivity test, we varied the α values for fresh snow, firn and ice within the uncertainties of their study. The fresh snow α variance results in ±0.5 m w.e. a −1 feedback, the firn variance in a ±0.43 m w.e. a −1 feedback and the ice variance in a ±0.17 m w.e. a −1 feedback. The test shows the high sensitivity to the albedo parametrization. Therefore, COSIPY simulated b clim always has to be evaluated against remote-sensing based geodetic mass balance or mass balance derived from the glaciological method through direct observations on the ground. If B clim lies within a reasonable range compared to the independent evaluation data, we understand the main purpose of applying COSIPY in identifying drivers of seasonal and interannual mass balance variability and in advancing process understanding of interactions between atmospheric drivers and B clim . These analyses are feasible and justified when overall calculated B clim falls within a reasonable range. This has also been demonstrated in other studies using medium complexity energy and mass balance models similar to COSIPY (e.g., [97][98][99][100][101][102]).

Glacier-Wide Climatic Mass Balance
Mean simulated B clim,a of the MB-years 1982 to 2019 is −0.48 m w.e. a −1 . Model simulations can be compared to ice volume changes derived by remote sensing studies of Kropáček et al. [38] and Shean et al. [7] on the mass budget of the Halji glacier. Those values are summarized in Table 6. Corresponding simulated B clim,a falls within the range of uncertainty in the case of the study by Kropáček et al. [38]. In the case of Shean et al. [7], the COSIPY simulated result is less negative than their estimate. Overall, the simulated negative B clim falls within a reasonable range and is in line with observed reduction in area, presented by the various outlines in Figure 1, the general glacier retreat in HMA (e.g., [3][4][5][6]) and in particular in western Nepal (e.g., [19,103,104]). A possible explanation for the more negative result of the Shean et al. [7] study are the simulated extreme negative MB-years 2016, 2017 and 2018. which are included in their study period, while they are not included in the Kropáček et al. [38] study. Table 6. Comparison of simulated annual glacier-wide climatic mass balance B clim,a (m w.e. a −1 ) with other studies using digital elevation model differences at the Halji glacier. The latter separates TP in liquid and solid precipitation and determines together with U 2 the density of fresh snow within COSIPY. A high correlation (r s = 0.54, p-value < 0.01) is found between the annual variability of SF c and B clim,a . The combination of SF c and surface melt patterns determines the variability of α, which has the highest correlation with B clim,a with r s = 0.8 (p-value < 0.01) on annual bases. August is the second last month in the MB-year but still B clim,a is extremely sensitive to α variability in August with r s = 0.91 (p-value < 0.01). Monthly deviations to their long-term mean of α in June, July and August reveal their strong coupling to monthly B clim variability in the same months. July, August and September are the months with the highest monthly variability of B clim in comparison with their long-term mean. Although monthly T 2 variability in the same months is smallest compared to the rest of the year this variability is nonetheless of major importance. These results are backed by the SSC, which reveal that peak in T 2 and TP sensitivity occurs during summer. In comparison with the SSC of the six glaciers from different mountain ranges in the world and including only Abramov glacier from HMA presented in Oerlemans and Reichert [43], Halji glacier is the only one with this behavior. Moreover, the difference between the sensitivity to December and November compared to September, August, and July is greater than for any glacier studied by Oerlemans and Reichert [43]. However, the applied mass balance model is not the same and differences resulting solely from the sensitivities of the applied models cannot be evaluated.

Mass-Balance
When looking at monthly mean deviations to their long-term mean in Figure 7, a possible explanation of the equilibrium state between 2008 and 2015 might be an increased TP amount, and therefore increased SF c between June and September, resulting in a positive albedo feedback (see [28]). This feedback can partly be linked to the high correlation of the WYM and B clim,a . The annual value of the WYM is ten years in a row negative between 2005 and 2015. In contrast, there are only three consecutive negative years of the WYM from 1984 to 2004. The WYM is a measure for the large-scale monsoon circulation intensity [79,105]. Besides the high correlations between WYM and B clim,a , a high correlation with r s = −0.45 (p-value < 0.01) between June values of WYM and α highlights again the importance of α.
The correlations between B clim and WYM and B clim and the forcing variables emphasize the importance of timing and strength of the monsoon for the Halji glacier. The correlations among the forcing variables T 2 , RH 2 , U 2 and Q LWin in June are further indications of the importance of the monsoon. A strong monsoon with an unstable atmosphere is accompanied by higher T 2 , higher U 2 , higher Q LWin and decreased RH 2 . Remarkably, other than the WYM, none of the other monsoon (e.g., Indian Summer Monsoon index, Australian Monsoon Index) or monsoon related (e.g., Indian Ocean Dipol index, Southern Oscillation Index) indices revealed a significant correlation within this study. The analysis of possible reasons is beyond the scope of this study and would be speculative at this point.
A couple of studies have shown the connection between POL and precipitation patterns (e.g., [82,106,107]), whereby they focus on the central and whole Tibetan Plateau, North China or winter precipitation patterns in the Karakoram or western Himalaya. In contrast, in case of Halji glacier, the high correlation of POL in June and B clim,a might rather be related to T 2 than to TP with r s = −0.47 (p-value < 0.01) between June values of T 2 and POL.
All results of this study reveal that a high amount of Halji glaciers interannual B clim,a variability can be accounted to the variability of T 2 and TP in summer, which can be at least partly be due to monsoon patterns. This finding is in accordance with the literature (e.g., [19,108]). As a result, Halji glacier is extremely sensitive to temperatures changes. An increase in T 2 of 0.5 K results in a B clim,a of −1.23 m w.e. a −1 (−0.75 m w.e. a −1 difference compared to reference simulation), which mainly affects B clim during summer. An increase in summer temperature has three main negative effects on the glacier mass budget for that kind of glacier [25]: (1) Enhanced melt, (2) lowered ratio between solid and liquid precipitation, which results in less surface accumulation, and (3) lower values of α as a result of a decrease in SF c , which enlarges Q SWnet and results in further enhancement of melt. The latter two effects are denoted as the snow-albedo feedback [26][27][28].

Conclusions
We simulated a B clim,a of −0.48 m w.e. a −1 between October 1981 and September 2019 for the Halji glacier in northwestern Nepal. Given the lack of direct mass balance observations, simulations were compared to geodetic mass balances derived from two remote sensing studies [7,38]. The simulation results are in a reasonable range compared to both of these remote sensing studies. The simulation also reveals high interannual variability of B clim,a . All results revealed the importance of the monsoon and T 2 and TP in summer for the variability of B clim,a . Only the combination of effects and resulting variability of α can explain the pronounced B clim variability in this season. The peak in summer of the calculated seasonal sensitivity characteristic back these findings. The variability of forcing variables in winter has a minor influence on B clim,a .
Comparison of ERA5-L derived COSIPY forcing variables with AWS measured variables revealed that most input data related uncertainties arise from the amount of TP and that downscaling and scaling procedures are necessary in order to obtain reliable model forcing data. A decrease of 20 % in TP results in an annual glacier-wide climatic mass balance of −1.77 m w.e. a −1 . Therefore, there is an urgent need in spatially higher resolved atmospheric datasets, to serve as forcing for long-term runs of mass balance models.
Presented approaches and statistics revealed that downscaling of T 2 , p s f c and RH 2 is straight-forward given that some observations in the field are available. Furthermore, we recommend to investigate the influence of the applied spatial resolution on B clim for each study site individually in order to identify the best trade-off between computational power and reasonable results.
COSIPY proves to be a powerful tool to identify climatic drivers of seasonal and interannual mass balance variability and to improve the process understanding of glacier responses to atmospheric forcing.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following acronyms are used in this manuscript:     Figure A3. Cumulative climatic mass balance of temperature adjusted simulation of Halji glacier from 1981 to 2020. We apply a temperature offset of −0.37 K to reach the long-term zero cumulative climatic mass balance which is needed for the computation of the seasonal sensitivity characteristic according to Oerlemans and Reichert [43]. 30 Figure A4. Differently resolved spatial grids over Halji glacier based on the glacier outline from Randolph Glacier Inventory 6.0 [45] and with a satellite image map as background from [55]. Colors within the inset map represent elevation [56].  , air temperature at 2 m (c), accumulated snowfall (d), incoming shortwave radiation (e) and incoming longwave radiation (f). The color bar is arranged so that deviations that act positively (negatively) on glacier-wide climatic mass balance are shown in blue (red).