Understanding monsoon controls on the energy and mass balance of Himalayan glaciers

. The Indian and East Asian Summer Monsoons shape the melt and accumulation patterns of glaciers in High Mountain Asia in complex ways due to the interaction of persistent cloud cover, large temperature amplitudes, high atmospheric water content and high precipitation rates. While the monsoons dominate the climate of the southern and eastern regions, they progressively lose strength westward towards the Karakoram, where the inﬂuence of Westerlies is predominant. Despite the major role of the monsoon in the Himalayas, a holistic understanding of their inﬂuence on the region’s glaciers is lacking 5 because previous applications of energy- and mass-balance models have been limited to single study sites. In this study, we use a full energy- and mass-balance model and seven on-glacier automatic weather station datasets from different parts of the Himalayas to investigate how monsoon conditions inﬂuence the glacier surface energy and mass balance. In particular, we look at how debris-covered and debris-free glaciers respond differently to monsoonal conditions. The radiation budget mostly controls the melt coming energy due to its dark and cold surface, small conduction distance to the ice and sensitivity to turbulence. The turbulent ﬂuxes ’work for’ the debris-covered glaciers by returning absorbed energy to the atmosphere, and tend to ’work against’ clean-ice and dirty-ice glaciers under monsoonal conditions. We thus expect the mass balance of debris-covered glaciers to react less sensitively to projected future monsoon changes than clean-ice and dirty-ice glaciers.


Introduction
High Mountain Asia (HMA) holds the largest ice volume outside the polar regions (Farinotti et al., 2019) and due to the large elevation range and vast geographic extent, HMA glaciers are highly diverse in character and hydro-climatic situation (Yao et al., 2012). Several large-scale weather patterns interact with the region's topography (Bookhagen and Burbank, 2010), caus-25 ing glaciers to contrast in terms of hypsometry (Scherler et al., 2011a) and accumulation and ablation seasonality . The Indian Summer Monsoon dominates the Central Himalaya and the Southeastern Tibetan Plateau during summer, and gradually loses strength moving towards the Karakoram, Pamir and Kunlun ranges in the east, where the influence of Westerlies is particularly strong. A more continental, monsoon-westerlies-influenced regime controls the Central Tibetan Plateau (Yao et al., 2012;Mölg et al., 2014), and the East Asia Monsoon influences the eastern slopes of the Tibetan Plateau 30 Yao et al., 2012). These major modes of atmospheric circulation do not only control the surface processes and runoff regimes of glaciers (e.g. Mölg et al., 2014Mölg et al., , 2012Kaser et al., 2010) but also lead to distinct responses of glaciers to climate change (Yao et al., 2012;Sakai and Fujita, 2017;Scherler et al., 2011b;Kraaijenbrink et al., 2017). For example, mass losses are high throughout most of HMA, and are particularly pronounced on the South-Eastern Tibetan Plateau, while glaciers exhibit a near-balance regime in the Karakoram, Pamir and Kun Lun (Farinotti et al., 2020;Gardelle et al., 2012;35 Brun et al., 2017;Shean et al., 2020).
Accurate glacier mass balance modelling is essential to assess glacier meltwater contribution to mountain water resources, and to predict future glacier states and catchment runoff. Physically-based models of glacier energy-and mass balance represent surface and sub-surface energy fluxes using physical equations to calculate the energy residual, i.e. the energy available for melt, and the glacier runoff. They have provided an understanding about the individual processes controlling the glacier mass 40 balance and the climatic sensitivity of glaciers in their specific hydro-climatic environment (e.g. Mölg et al., 2012). Summeraccumulation type glaciers in HMA experience simultaneous accumulation and ablation and their mass balances are known to be highly sensitive to climatic variability during the monsoon season (Fujita and Ageta, 2000), when warm air temperatures and high moisture influx coincide. Using energy balance modelling for an inter-annual study for the Central Tibetan Zhadang glacier, Mölg et al. (2012) demonstrated that the timing of monsoon onset and the associated albedo variability can change 45 melt-rates substantially in subsequent years. At the same time, they observed a de-coupling of the glacier mass balance from the Indian Summer Monsoon's control during the core monsoon season. Mölg et al. (2014) explain the mass balance variability of Zhadang Glacier as being controlled by both the Indian Summer Monsoon onset and remotely controlled by mid-latitude Westerlies. Combining energy balance with weather forecast modelling, Bonekamp et al. (2019) identify the timing and quantity of snowfalls as the main source of differences in mass-balance regimes between the Shimshal Valley in the Karakoram and 50 been the focus of devoted modelling studies (Giese et al., 2020;Collier et al., 2014). Moreover, the representation of latent heat due to evaporation (Giese et al., 2020;Steiner et al., 2018) and atmospheric stability correction for turbulent fluxes were shown to be important to improve the simulation of sub-debris melt (Reid and Brock, 2010;Mölg et al., 2012). Model implementation, however, remains complex and difficult to validate and transfer. Nevertheless, Nicholson and Stiperski (2020) showed that the turbulent conditions over debris-covered and clean-ice sites are similar enough to be numerically treated in similar ways, i.e. 70 the Monin-Obukhov similarity theory can be leveraged in the same way for the calculation of the turbulent fluxes, but with different parameters.
Observations of glacier surface meteorology, a prerequisite for accurate energy balance modelling, exist only for a few glaciers in HMA, and even fewer records exist for debris-covered glaciers. Direct observations of glacier mass balance in HMA also remain sparse, and remote sensing observations are hindered by heavy cloud cover during the monsoon season. 75 Previous studies explicitly dealing with the imprint of the monsoon on the surface thermal properties of glaciers remained limited to individual clean-ice glaciers in the Central Tibetan Plateau (Mölg et al., 2012(Mölg et al., , 2014. Our main goal is to improve understanding of monsoon controls on glaciers in the Himalayan region by leveraging available automatic weather station (AWS) records in the region. Our specific objectives are: 1) to understand, in detail, the glacier energy and mass balance and related controls at seven Himalayan study sites in a robust and systematic manner; 2) to assess the effects of monsoonal con-80 ditions on debris-covered and clean ice glaciers; 3) to investigate whether these effects are generalisable; and 4) to discuss possible implications of our findings under climate change. We address these objectives by applying the glacier energy and mass balance module of a land surface model suited to both debris-covered and clean-ice glaciers. We apply the model at the point scale of individual AWSs, driven by high-quality in situ meteorological observations that guarantee accurate energy balance simulations, not affected by extrapolation of the meteorological forcing. By identifying and discussing the key surface 85 processes of glaciers and their dynamics under monsoonal conditions, this study promotes their appropriate representation in models of glacier mass balance and the hydrology of glacierised catchments.
2 Study sites and data 2.1 Sites and observations 90 In situ observations from seven on-glacier automatic weather stations in different environments along the climatic gradient of the Himalayas were gathered and are used for forcing and evaluation of the model (Table 2). Our seven study sites are located in the Central and Eastern Himalayas and cover a range of glacier types and local climates (Figure 1, 2 and Table 1). The seven sites include both spring-(24K, Parlung No.4) and summer-accumulation glaciers (all others) as indicated by the proportion of monsoon precipitation to the annual precipitation ( Figure 3). Langtang, Lirung and Yala are neighboring glaciers found in the Langtang Valley (Figure 1), which has an extensive history of glaciological and hydrological reasearch. The Langtang Valley is strongly influenced by the Indian Summer Monsoon (∼ June to October), during which more than 70% of the annual precipitation arrives (Figure 3 and Table 1), while the period from November to May is a drier season Immerzeel et al., 2012). The Langtang Valley has been a site of extensive glaciological (e.g. Fujita et al., 1998;Stumm et al., 2020), meteorological (Immerzeel et al., 2014;Heynen et al., 2016;Steiner and Pellicciotti, 2016;Bonekamp et al., 2019;100 Collier and Immerzeel, 2015) and hydrological (e.g. Ragettli et al., 2015) investigations. On-glacier AWSs were installed during the ablation season on Lirung (2012Lirung ( -2015 and Langtang (2019) glaciers, and year-round on Yala (2012-ongoing) ( Table   2). Both Lirung and Langtang are valley glaciers that have heavily debris-covered tongues, but the tongue of Lirung has disconnected from the accumulation zone ( Figure 2). Yala is a considerably smaller clean-ice glacier, with most of its ice mass located at comparably high elevation. It is oriented to the southwest and has a gentle slope (Fujita et al., 1998) (Figure 2 and   105   Table 1).
North Changri Nup Glacier (hereafter Changri Nup Glacier) is a debris-covered valley glacier located in the Everest region in Nepal ( Figure 1). The southeast-oriented, avalanche-fed glacier discharges into the Koshi River system. The local climate is similar to the one of the Langtang Valley, with 70-80% of precipitation falling during monsoon (Vincent et al., 2016) (Figure   2, 3 and Table 1).

110
24K and Parlung No.4 glaciers are located on the southeastern Tibetan Plateau, feeding water into the upper Parlung Tsangpo, a major tributary to the Yarlung Tsangpo -Brahmaputra River. The summer climate is characterized by monsoonal air masses reaching the Gangrigabu mountain range from the south through the Yarlung Tsangpo Grand Canyon. 24K Glacier is an avalanche fed valley glacier with a debris-covered tongue, located 24 km from the town of Bome . It is small in size, oriented to the northwest and surrounded by shrubland (Figure 1, 2  glacier, which is north-east oriented, considerably larger than 24K and located 130km (south-east) from Bome (Yang et al., 2011) (Figure 1 and Table 1). Full automatic weather stations were installed in the ablation zones of both glaciers in 2016 and subsequent years (Table 2). Hailuogou Glacier, the second-largest of our study sites ( Figure 2) is located on the eastern slope of Mt. Gongga in the easternmost portion of the southeastern Tibetan Plateau (Figure 1). It is located at low elevation and covered with thin and patchy 120 debris, leading to high annual ablation rates ( Figure 2 and Table 1). The local climate influenced by the East Asia Monsoon with typically only 50 to 60% of the annual precipitation arriving during the monsoon period (Figure 1 and 3). The debriscovered tongue is connected to a steep and extensive accumulation zone via a large icefall, but avalanching is the primary mass supply mechanism through the icefall to the valley tongue (Liao et al., 2020), and a dynamic disconnect is expected to occur in the near future. Weather stations were installed at three nearby off-glacier locations and one on-glacier site during 2008, while 125 precipitation was measured at Alpine Ecosystem Observation and Experiment Station of Mt. Gongga, within 1.5 km from the glacier terminus (Table 2).

Climatic and meteorological conditions
Here, we use the monthly averaged ERA5-Land reanalysis data (Muñoz Sabater, 2019) to provide an overview of the long 130 term climatic patterns, e.g. the average monsoonal regime from June through September, and evaluate the representativeness Table 1. Characteristics of the study sites. Planimetric glacier and debris surface areas, mean elevation, slope and aspect were calculated using the updated Randolph Glacier Inventory 6.0 by Herreid and Pellicciotti (2020)   Incoming shortwave radiation ( Figure 3b) shows a clear peak before monsoon onset at all sites. A smaller secondary peak is reached just after the monsoon in October at the Central Himalayan sites, but not at the Eastern Himalayan sites. Interruptions 135 in monsoonal overcast conditions (break periods) seem to be more common at the eastern sites, leading to occasional secondary peaks in incoming shortwave radiation during monsoon.
Average mean monthly 2 m air temperatures have a similar pattern at all study sites (Figure 3a), with a slow increase from January to a peak between July and August, just after peak monsoon, and a steeper decline from post-monsoon into winter. A similar regime is followed by LW ↓ , with highest values reached during the core monsoon ( Figure 3c).

Tethys-Chloris energy balance model
We use the hydrological, snow and ice modules of the Tethys-Chloris (T&C) land surface model (Fatichi et al., 2012;Mastrotheodoros et al., 2020;Paschalis et al., 2018;Botter et al., 2020) to simulate the mass and energy balance of the seven study glaciers. T&C simulates, in a fully distributed manner, the energy and mass budgets of a large range of possible land surfaces, including vegetated land, bare ground, water, snow and ice. Here, we apply the model at the point scale of the AWS locations 160 to simulate the energy fluxes of the underlying surface and subsurface, which can comprise snow, ice and supraglacial debris cover layers, according to the local and dynamic conditions. The melt and accumulation of ice and snow, and the ice melt under debris are also explicitly simulated. The surface energy balances for the three different possible surfaces are for snow, for debris cover, and for ice, where R n [W m −2 ] is the net radiation absorbed by the snow/debris/ice surface, Qv [W m −2 ] is the energy advected from precipitation, Q f m [W m −2 ] is the energy gained or released by melting or refreezing the frozen or liquid water that is held where are the incoming atmospheric and outgoing longwave radiation components, respectively. In this study α is given as an input to the model based on the AWS observations. We prescribe α for all surface types as the daily cumulated albedo, which is the 24 hour sum of SW ↑ divided by the sum of SW ↓ centred over the time of observation (van den Broeke et al., 2004).

Incoming energy with precipitation
For calculating the incoming energy with precipitation, rain is assumed to fall at air temperature (T a ) when positive, with a 190 lower boundary of 0 • C. Snow is assumed to fall at negative T a with an upper boundary of 0 • C. Here, Q v is the energy required to equalize the precipitation temperature with the surface temperature T s and is therefore calculated as where is the density of water and P r,liq [mm], P r,sno [mm] are the rain-and snowfall intensities, respectively.

Phase changes in the snow pack
The snow pack has a water holding capacity Sp wc described in section 3.2.2. The energy flux gained/released by melting/refreezing the frozen/liquid water that is held inside the snow pack is calculated as:

Turbulent energy fluxes
Over snow, debris and ice surfaces, the sensible energy flux is calculated as where is the specific heat of air at constant pressure, ρ a [kg m −3 ] is the density of air. The aerodynamic resistance r ah [s m −1 ], which is also a function of wind speed (W s) is calculated using the simplified solution of the Monin-Obokhov similarity theory (Mascart et al., 1995;Noilhan and Mahfouf, 1996). The roughness lengths of heat (z 0h [m]) and water vapour (z 0w [m]) used in 210 the calculation of the aerodynamic resistance are equal in T&C (z 0h = z 0w ), and (z 0h = z 0w = 0.1z 0m ). The roughness length of momentum (z 0m ) is set to 0.001 m for snow and ice surfaces (Brock et al., 2000), while we optimize it against the surface temperature for debris (see section 3.3).
Correct estimates of the latent energy flux due to water phase changes at the surface are important to accurately model glacier 215 melt, especially under moist conditions (Sakai et al., 2004). Phase changes between the water and gas phase and the resulting energy fluxes are considered over all surfaces. The latent energy is limited by the availability of water in the form of ice and snow, or in the case of a debris surface, by the amount of water intercepted (interception storage capacity is set to 2mm). The latent energy flux is estimated from:

220
where λ s is the latent energy of sublimation defined as λ s = λ + λ f , with λ = 1000 (2501.3 − 2.361 T a ) [J kg −1 ] as the latent energy of vaporisation. q sat is the surface specific humidity at saturation at T s , q a is the specific humidity of air at the measurement height and r aw the aerodynamic resistance to the vapour flux, which we assume equals r ah .

Ground energy flux
The definition of the ground energy flux G [W m −2 ] differs based on the surface type. In the case of snow, it is equal to the 225 energy transferred from the snowpack to the underlying ice or debris surface, where in the assumption of a slowly changing process, G can be approximated with the temperature difference of the previous time step (t-1), which allows to solve for G outside the numerical iteration to find the snow surface temperature of the current time step: where k sno [W K −1 m −1 ] is the thermal conductivity of snow and d sno [m] is the snow depth. For ice in the absence of snow 230 and debris, it is the energy flux from the ice pack to the underlying surface or to the ice at a depth of 2 m: is the temperature of the underlying ice, and d ice [m] is the relevant ice thickness. For debris, which was discretised into eight layers at all debris-covered sites, G is the energy flux conducted into the debris layers. Its calculation is for a given time t and depth z 235 where k d [W K −1 m −1 ] is the debris thermal conductivity (see section 3.3) and T deb (z, t) [°C] is the debris temperature at time t and depth z. G(z, t) can be included in the heat diffusion equation as such: where cv d is the debris heat capacity. Under the assumption of homogeneous debris layers, κ [m 2 s −1 ] as the debris heat diffusivity replaces the term k d cvs and equation (12) can be written as: The heat diffusion equation (13) is solved using iterative numerical methods. This way, the debris temperature profile T deb (z, t) is solved together with G(z, t) at any depth and time. The conductive energy flux at the base of the debris is used to heat the ice and to calculate ice melt once above the melting point.

Precipitation partition
Input precipitation is required to be partitioned into solid P r sno and liquid P r liq precipitation, because of the differing impacts of snow and rain on the energy and mass balance. For this study, the precipitation partition method described by Ding et al. (2014) was implemented into T&C. This parameterisation has been developed specifically for High Mountain Asia based on a 250 large dataset of rain, sleet and snow observations, and does not require recalibration. It determines the precipitation partition based on the wet-bulb temperature, station elevation and relative humidity and allows for sleet events, as a mixture between liquid and solid precipitation. Ding et al. (2014) found the wet-bulb (T wb ) to be a better predictor than T a of the precipitation type, that the temperature threshold between snow and rain increases at higher elevations, and that the probability of sleet is reduced in conditions of low relative humidity.

Water content of the snow, ice and debris layers
The water content of ice is approximated with a linear reservoir model. The liquid water outflow is proportional to the ice pack water content Ip wc [mm w.e.], which is initiated when Ip wc exceeds a threshold capacity, prescribed as 1% of the ice water equivalent (IW E [mm w.e.]). The Ip wc is the sum of ice melt and liquid precipitation, minus the water released from the ice pack. The water released is the sum of the ice pack excess water content plus the outflow from the linear reservoir, given as 260 I out = Ip wc /K ice , where K ice is the reservoir constant which is proportional to the ice pack water equivalent. Unlike within snow packs, Q f m is not accounted for within the ice pack, since water is presumed to be evacuated quickly from the ice due to runoff without refreezing.
The water content of the snow pack Sp wc [mm w.e.] is approximated using a bucket model, in which outflow of water from the snow pack occurs when the maximum holding capacity of the snow pack is exceeded. Following the method of Bélair 265 et al. (2003), the maximum holding capacity of the snow pack is based on SW E, a holding capacity coefficient and rho sno .
Snowmelt plus liquid precipitation, minus the water released from the snow pack gives the current Sp wc . If T sno is greater than 0°C then the snow pack water content is assumed to be liquid, whereas otherwise it is assumed frozen.
For supraglacial debris, both observations and methods for modelling its water content are lacking. We thus use a simplified we assume debris to have a dynamic interception storage, which can hold a maximum of 2 mm water at all debris-covered sites and can be refilled by snowmelt or liquid precipitation. The evaporative flux from the debris and the latent energy flux of evaporation is therefore limited by this interception storage.

Glacier mass balance 275
The mass balance calculation of snow and ice is rather similar, so they will be described together here. Calculations are performed for snow if there is snow precipitation during a timestep or the modelled SW E at the surface is greater than zero.
Net input of energy to the snow or ice pack will increase its temperature, and after the temperature has been raised to the melting point, additional energy inputs will result in melt. The change in the average temperature of the ice or snowpack dT is controlled using: Where dt is the time step is IW E or SW E before melting and limited to a maximum of 2000 mm, assumed to be the water equivalent mass exchanging energy with the surface. Energy inputs into an iso-thermal ice/snow pack The water equivalent mass of the snow/ice pack after melting W E(t) [mm w.e.] is updated conserving the mass balance following: Here E = λE/λs [mm] is the sublimation from ice and snow. The snow density is assumed to be constant with depth and calculations are performed assuming one single snow pack layer. The snow density evolves over time using the method proposed 290 by Verseghy (1991) and improved by Bélair et al. (2003). In this parameterisation the snow density increases exponentially over time due to gravitational settling and is updated when fresh snow is added to the snowpack. Two parameters are required in this scheme, ρ M 1 sno and ρ M 2 sno [kg m −3 ], which represent the maximum snow density under melting and freezing conditions, respectively. The depth of the ice pack can be increased through the formation of ice from the snow pack (ice accumulation), which is prescribed to occur if the snow density increases to greater than 500 kg m −3 (a density associated with the firn to ice 295 transition) and at a rate of 0.037 mm h −1 (Cuffey and Paterson, 2010). The density of ice is assumed constant with depth and equal to 916.2 kg m −3 .

Debris parameters
A major challenge in physically based mass-balance modelling of debris-covered glaciers is the assignment of appropriate debris properties. Besides the debris thickness, which was measured at the AWS location, the thermal conductivity k d , the roughness length z 0 of the debris surface, the surface emissivity d , the debris volumetric heat capacity cv d and and the debris density ρ d have to be assigned. While the latter three can be quantified using literature values, there is more uncertainty about k d and z 0 , two highly variable quantities that are difficult to measure in the field, but which EB models are highly sensitive to.
We thus choose to optimize them, since our primary requirement is an accurate representation of the energy and mass balance: (1) in a first step, we optimize k d simulating only the conduction of energy through the debris during snow free conditions, with 305 the LW ↑ -derived surface temperature T s,LW as an input, the ice melt as the target variable and the Nash-Sutcliffe Efficieny N SE [−] as performance metric.
(2) Next, we run the whole energy balance model and optimize z 0 for snow-free conditions, with all required meteorological inputs, and the optimal k d from step (1), while comparing modelled T s against T s,LW , using N SE as performance metric. The resulting parameters are given in Table 3.

310
We calculate the uncertainty associated with all energy and mass balance components by performing 10 3 Monte Carlo simulations for each study site at the AWS location. We vary three debris parameters (k d ,z 0m , d ), debris thickness h d , as well as six measured model input variables: air temperature T a , the vapor pressure at reference height ea[−], SW ↑ , SW ↓ , LW ↓ , the total precipitation before partition P r, and the wind speed W s. Measured outgoing shortwave radiation SW ↑ was included into the Monte Carlo set, as it determines our input α, as discussed in Section 3.1.1. While the parameter uncertainty range was defined 315 based on previously published values for debris (e.g. Yang et al., 2017;Rounce et al., 2015;Evatt et al., 2015;Reid and Brock, 2010;Nicholson and Benn, 2006;Rowan et al., 2020;Lejeune et al., 2013;McCarthy, 2018), the debris thickness measurement uncertainty was given with a range of 1cm and the range for the meteorological inputs was set based on the respective sensor uncertainties (see Table 4). All uncertainties were equally distributed around the standard parameter values and observations. Uncertainties are given as one standard deviation of the error of the Monte Carlo runs against the 320 standard run.

Modelled mass balance
The model accurately reproduces the measured surface height change (ablation and accumulation) at both debris-covered and clean-ice glaciers (Figure 4). The maximum uncertainties associated with each model output ranges from ±4% (Parlung No.4, 325 Figure 4f) to ±15% (Yala, Figure 4c). Where Ultrasonic Depth Gauge (UDG) records are available (Lirung, Langtang, Yala, Changri Nup), the deviations of the simulations from the observations stay within the uncertainty range (Figure 4a-d). We decided to not consider the UDG record from Changri Nup after a large August snowfall, as variables describing the surface state (e.g. α, LW ↑ ) following this event indicate a discontinuous snow cover at the AWS location, while the UDG, which is some meters away from the AWS, shows continuous snow cover with depths at the order of tens of centimeters. This discrep-  (Figure 5f). Similarly, intermittent snow cover protects the ice from melting at the two highest sites (Yala and Changri Nup) during the summer months (Figure 5c and d). Sublimation from ice and snow represents a very small share of the total mass losses, and ranges 345 from 0.01% (Lirung, Figure 5a) to 1.2% (Changri Nup, Figure 5d). It mostly occurs under dry conditions during premonsoon at the highest sites (Changri Nup, Yala). Debris cover is clearly another important control of total mass losses: The presence of debris at least 10 cm thick (Changri Nup, Lirung, Langtang and 24K) causes average melt rates to be comparatively low (8.2, 8.5, 6.1 and 19.4 mm d −1 , Figure 5a,b,d and e), while at the 'dirty-ice' site (Hailuogou, with thin and patchy debris, approximated to be 1 cm), the melt rate is higher (42.7 mm d −1 , Figure 5g) over the simulation period, peaking in mid-July.

350
Relatively high average melt rates are also reached at the clean ice sites Parlung No.4 (34.6 mm d −1 , Figure 5f) and Yala (16.6 mm d −1 , Figure 5c), despite the high altitude of the latter.

Modelled energy balance 355
The main energy source on all glaciers and during all seasons is SW net (Figure 6). At all sites and over the entire modelling period, LW out is the largest energy sink, but it remains comparably stable in magnitude between the seasons ( Figure A8). The turbulent heat fluxes (H and LE), and dQ (G on debris) act as energy sinks on average at all sites. G can however act in both ways, as an energy source when warming the glacier, and as an energy sink when warming the air. In our definition G sums up all types of conductive energy fluxes in the snow-debris-ice column. Part of the energy is used for heating the snow, ice 360 or debris-ice surface layer until melt occurs (dG, Table 5). The relative contributions of the individual sinks strongly depends on the surface type: LW net is an important sink at the clean-ice sites, while the turbulent fluxes remain small there. On the debris-covered glaciers with debris >10 cm thick, where LW net plays a relatively small role, the turbulent fluxes act as major energy sink. Where debris is thin, they instead act as additional energy sources.
In the following, we consider how energy fluxes change when moving from pre-monsoon to monsoon. Note that the sign of 365 this change depends on the sign of the original flux, which is negative when the flux acts as an energy sink (energy away from the surface) and positive when the flux acts as a source (energy towards the surface). As an example, a negative change in LE, when acting as a sink, means an increase in its magnitude.
The seasonal imprint of the monsoon on the energy balance is very clear at all study sites on average and diurnally: There is a strong reduction in the magnitude of SW ↓ during monsoon across sites (with decreases ranging between −41.8 and

370
−135 W m −2 at the seven sites, Figure 7 and  Table 5) as a consequence of ephemeral snow cover ( Figure 4). LW ↓ on the other hand increases at all sites (by +15.7 to +57.0 W m −2 ), and LW ↑ slightly increases in magnitude as well (by −1.0 to −13.3 W m −2 ) (Figure 7 and Table 5).

375
At the debris-covered sites, the reduction in R net from pre-monsoon to monsoon is balanced by a reduction in the magnitude of H (with changes of +15.6 to +67.7 W m −2 between sites), which acts as a major energy sink during both sub-seasons ( Figure   6, 7 and Table 5). LE, which plays a small role as a sink compared to H during pre-monsoon, increases in magnitude at the debris-covered sites (with changes ranging from −2.7 to −24.4 W m −2 , Figure 7 and Table 5), and acts as a considerable heat sink during the monsoon. As a result of these contrasting changes, dQ remains fairly similar between the pre-monsoon and 380 monsoon at the debris covered sites (with changes ranging between +1.3 and −3.3 W m −2 ) (Figure 7a,c,e,g and Table 5). This is partly a consequence of a compensation between increased/reduced melt rates before/after noon ( Figure A10). An exception is the thin-debris site Hailuogou, where dQ increases in magnitude (by −26.2 W m −2 ), mostly as a result of warmer nights ( Figure A10).
At the clean-ice sites, LW ↑ is initially larger than LW ↓ in the pre-monsoon, but the two become almost equal while adjusting to 385 monsoon conditions, causing the LW net to be close to zero (LW ↓ : +40.7 to +57.0 W m −2 , and LW ↑ : −3.1 to −11.7 W m −2 ), 12.7 W m −2 , Table 5). The turbulent fluxes, however, remain comparably small in magnitude at the clean ice sites ( Table 5).

Turbulent fluxes at debris-covered sites and their controls
The turbulent fluxes LE and H are important heat sinks on debris-covered glaciers and heat sources on the dirty-ice glacier, while they remain small on the clean-ice glaciers ( Figure 6). Their magnitude largely determines the energy that is conducted 400 into the debris and used for ice melt. We observe a decrease in the magnitude of H and an increase in the magnitude of LE going from pre-monsoon to monsoon at the debris-covered glaciers, while both increase as sources of energy at the dirty ice glacier (Figure 7).
The predictive power of the temperature gradient between surface and air gT [°C −1 ] and W s as well as their combination for determining H and LE were assessed using a univariate polynomial regression model for the single predictors and a multiple 405 polynomial regression model using both variables, and both models had 2 degrees of freedom. H is largely controlled by the temperature gradient between surface and air (g T ) on debris-covered glaciers, which explains between 75 and 99% of the variability of H ( Figure A9), and g T decreases during monsoonal conditions by −0.05 to −1.44 • C m −1 (Table 6). Here it is clear that a smaller temperature gradient between surface and air during the monsoon reduces the magnitude of H as a sink. Wind emerges as the most important control of H at the dirty-ice glacier, explaining up to 89% of variability ( Figure   410 A9). The mean magnitude of W s increases at that site from 1.23 m s −1 in pre-monsoon to 2.10 m s −1 in monsoon (Table 6).
A cold surface in combination with a wind-enhanced turbulence and fast turnover of warm air mass results in H becoming a potent heat source on the dirty-ice glacier. From a physical perspective, the same holds for clean-ice glaciers experiencing simultaneously high wind speeds and high air temperatures. This is, however, not observed on Yala and Parlung No.4, where air temperatures stay comparably low and average wind speed decreases slightly in monsoon compared to pre-monsoon (Table   415 6).
Neither RH, g T , or W s on their own explain the variability of LE across sites ( Figure A9). LE however increases consistently from pre-monsoon to monsoon together with an increase in the duration of moisture availability at the surface of debris-covered glaciers (by 21 to 63%, Table 6). In fact, evaporation and the related LE release tends to be water-limited during pre-monsoon, but energy-limited during monsoon (Figure 8). This implies that the availability of moisture is driving the increase of LE from 420 pre-monsoon to monsoon. However, on the dirty-ice glacier, LE acts as a heat source across sub-seasons and switches sign during monsoon at the clean-ice sites. This happens when T a is consistently higher than T s , causing condensation when the warm air touches the cold glacier-, or thin debris-surface.    . Budyko diagram with the 5-day mean potential evaporation rate during snow-free conditions (Epot) relative to the mean available intercepted water (In) on the x-axis, and the actual evaporation rate during snow-free conditions (Eact) relative to In on the y-axis. Only debris-covered glaciers, where LE is a sink term on average, are shown. Table 6. Mean cloud-cover fraction (ccf ), relative humidity (RH), temperature gradient between surface and air (gT ), wind speed (W s) and the percentage of time during which the debris is modelled to hold intercepted water (In) for each site and season, also indicating percent changes between the sub-seasons.

All glaciers
Overcast conditions caused by monsoon increase LW ↓ and with that the magnitude of T a at all sites (Figures 6 and 7). SW ↓ on the other hand is reduced at all glacier surfaces due to the reflection and scattering of persistent, heavy clouds. While the changes of these two primary sources of energy partly counteract each other in terms of energy supply to glacier surfaces during monsoon, we find that net incoming radiation decreases in magnitude everywhere. Another major control on the energy and mass balance of all glaciers are the snowcover dynamics, which in turn are driven by the precipitation seasonality and the partition of precipitation into rain and snow. For example, in the case of Parlung No.4, the onset of glacier melt was delayed until well after monsoon onset, until all snow had disappeared ( Figure 5). After snow has melted out, ephemeral snowcover from monsoonal precipitation increases surface albedo and raises SW ↑ , protecting the ice and suppressing melt rates throughout the 440 summer. This is especially relevant at high elevation sites (Yala, Changri Nup). An interruption of the monsoon at 24K occurred in August 2016, possibly associated with an El Niño event (Kumar et al., 2006). During this interruption the energy balance returned to a pre-monsoonal regime of clearer skies, more pronounced diurnal temperature amplitudes, low precipitation rates and lower relative humidity ( Figure A5), resulting in higher melt rates during that period (Figure 6e). Our analysis shows that some effects of monsoon are common for all surface types, while the presence or absence of debris and its thickness control 445 how the incoming energy is absorbed and transmitted to the ice.

Debris-covered glaciers
Monsoon conditions in combination with debris cover affects the energy balance in a way that stabilizes melt rates over the ablation period, despite increasing air temperatures typical of monsoon ( Figure 6). Enhanced melt during the night and morning in monsoon, resulting from higher T a , is partly offset during the cooler afternoon hours ( Figure A10) when the wet air masses 450 usually arrive in this region, bringing intense cloud-cover and precipitation. While H decreases as a consequence of a smaller average temperature gradient between surface and air, more latent energy is released from the wet debris surface, which shifts from a water-limited process during pre-monsoon to an energy-limited process during monsoon (Figure 8). It has been known from studies at individual sites that debris cover protects the ice by returning energy to the atmosphere in the form of turbulent fluxes  and that the turbulent fluxes can be a major component in the energy balance during both dry and 455 wet conditions (Steiner et al., 2018). Accounting for the debris water content through the inclusion of a simple interception storage has allowed us to identify the importance of the latent heat flux in the cooling effect of evaporation from debris during monsoon, a process that has often been neglected in previous modelling studies. We have been able to quantify how cloud overcast and additional moisture modify the energy balance of debris-covered glaciers, and especially the turbulent fluxes, to result in a melt-equalizing effect.

Clean-ice glaciers
In contrast to debris-covered glaciers, when clean-ice glaciers are snow-free and the ice has been heated to the melting point, almost all net radiation leads directly to ice melt, while the turbulent flux contribution remains small. When entering the monsoon period, the latent heat flux switches sign, changing from sublimation to condensation, which adds energy to the 465 surface instead of removing it. This behaviour has previously been observed at sites with a 'southern influence' (Yang et al., https://doi.org/10.5194/tc-2021-97 Preprint. Discussion started: 10 May 2021 c Author(s) 2021. CC BY 4.0 License.

Dirty-ice glacier
At our site with thin debris, the effects that we have observed for debris-covered and clean-ice sites combine to create a 470 melt-enhancing effect that becomes particularly potent during monsoon: the dark debris surface absorbs almost 90% of SW ↓ in the case of Hailuogou. With a short conduction length, the energy influx leads almost directly to melt. The cold surface favours condensation and a strong temperature gradient between the surface and the air (g T pre-monsoon:−2.08, monsoon: −2.61 • C m −1 ). Driven by higher wind speeds ( that thin debris causes higher melt rates than at both clean-ice sites and sites with thicker debris cover (Reznichenko et al., 2010;Östrem, 1959;Reid and Brock, 2010;Fyffe et al., 2020), especially in humid environments (Evatt et al., 2015), e.g. the location of Hailuogou glacier.

Sensitivity of seasonal flux changes to elevation and debris thickness
Our results are derived from simulations at one location (AWS) on each glacier. To understand how representative they are of 480 conditions across the glacier ablation zone at each site, and across the possible range of debris thicknesses in particular, we conduct a sensitivity experiment at each site. We re-run the model synthetically varying the AWS elevation to represent the range of elevation of each glacier ablation zone by applying a T a lapse rate of 0.6°C/100m and, for the debris-covered sites, by varying also the debris thickness in the range 10-80 cm (for ranges and steps see Table 7). Using the station-measured, accumulated albedo is not appropriate during this experiment, due to changing snow conditions with varying elevation. We 485 therefore include the parameterisation introduced by Ding et al. (2017) for modelling α. From the resulting range of EB flux outputs, we calculate the range of expected changes for the entire ablation zone when moving from pre-monsoon to monsoon (∆-range). This allows us to place our results in the context of the changes that can be expected over the entire ablation zone, given its elevation span and debris thickness variability. Figure 9 shows that even accounting for the range of conditions across each glacier ablation area, the pattern of pre-monsoon to monsoon difference in flux components, and importantly dQ, remain 490 similar for debris-covered sites: The ∆-range of dQ stays within the uncertainty range, with the exception of Langtang, where the unrealistic combination of relatively thin debris and low elevation causes high dQ ∆-range. This lends confidence to the results obtained at the individual AWS locations. Although we adjusted forcing data for elevation in this exercise, we could not represent the effects of variable debris thicknesses in modifying 2 m meteorological variables (Steiner and Pellicciotti, 2016;Shaw et al., 2016). This comes with the assumption that surface-atmosphere interactions are negligible compared to 495 the altitudinal patterns and temporal changes. While this might be acceptable at thicker debris sites, it is more questionable at Hailuogou, where the observations were taken above thin and cold debris. However, also at this site, the ∆-range ends up to be small ( 5 W m −2 ) and close to zero when debris between 10 and 80 cm thickness is applied artificially. Table 7. Ranges of elevations and debris thicknesses used for the sensitivity runs, including the glacier terminus elevation (min), the AWS elevation (AWS) and the upper debris limit on debris-covered glaciers or to the approximated ELA elevation on clean-ice glaciers (max). We also show the range of debris thicknesses hm modelled for debris-covered glacier sites. All combinations of elevations and debris thicknesses were used.  to debris thickness varied from 10 and 80 cm combined with the elevation range of debris-covered ablation zone of the individual glaciers;

Limitations
Our modelling approach allowed us to constrain the energy and mass balance at the AWS location in the best possible way 500 for the different glacier surface types. Two debris-specific parameters (k d and z 0 ) were optimised in two separate steps at the debris sites (see section 3.3) and α was derived from measured SW ↓ and SW ↑ and given as an input at all sites. While this is not a feasible approach at the scale of entire glaciers, for application at other sites, or for future simulations, as it requires specific observations of radiation and ablation, both choices were made here in order to constrain the energy balance with observations to the greatest degree possible. Requisite data (e.g. thermistor strings, wind towers, high-resolution topography) 505 to estimate these two parameters with established methods were not uniformly available across our study sites. These observational methods themselves can be sensitive to the measurement period (e.g. Nicholson and Benn, 2006) and are subject of ongoing research (e.g. Rounce et al., 2015;Miles et al., 2017;Rowan et al., 2020;Quincey et al., 2017), so we concluded that the calibration approach here was best suited to the application at all sites. In this respect, all optimized values fall within the expected range based on prior energy-balance studies on debris-covered glaciers Rounce et al., 2015;Evatt 510 et al., 2015;Reid and Brock, 2010;Nicholson and Benn, 2006;Rowan et al., 2020;Lejeune et al., 2013;McCarthy, 2018), but consistent measurements of debris parameters remain a community research priority to ensure robust modelling of sub-debris melt. To represent evaporation and the latent heat flux from the debris, we assigned an interception storage of 2 mm water to the debris surface, as, to the authors' knowledge, very few measurements have been collected of the debris layer's interception capacity and moisture content. This is important because the thermal properties of debris vary with 515 moisture fluctuations inside the debris (Nicholson and Benn, 2006;Collier et al., 2014;Steiner et al., 2018), so the moisture retention properties of the surface debris layer may play an important role in regulating energy transfer. However, these thermal properties have been found to be stable during the core monsoon months (Rowan et al., 2020), so we kept k d constant in time at each site in our model. Additional, targeted investigations should examine the possible role of water storage and mobility within debris layers (Giese et al., 2020;Collier et al., 2014). 520

Implications
Our results show that the surface type plays a large role in modulating a glaciers' response in the seasonal transition from premonsoon to monsoon: we find a melt-enhancing effect of thin debris, driven by enhanced energy uptake attributable to low albedo and turbulent energy fluxes. We also find that, on clean-ice glaciers, LE switches seasonally to become a source instead of a sink of energy to the surface, while net radiation determines melt. Importantly, we additionally identify a melt-equalizing 525 effect of debris cover under monsoon conditions at a number of sites. Monsoon-influenced, summer-accumulation glaciers (such as Langtang, Lirung, Yala, and Changri Nup) have been previously shown to be especially vulnerable to warming due to a decrease in accumulation and an enhancement of ablation with lowering albedo (Fujita and Ageta, 2000), and our results confirm that net shortwave radiation is the key control on monsoon-period melt rates for clean ice glaciers. Our results also emphasize that the longevity of pre-monsoon snowcover into the monsoon period is a key control on melt rates, supporting 530 past findings that the strength and timing of the monsoon onset has a profound impact on small mountain glaciers (Mölg et al., 2014(Mölg et al., , 2012 through the phase change of precipitation in the transition to monsoon conditions (Fujita and Ageta, 2000;Ding et al., 2017;Zhu et al., 2018).
The distinct responses of the surface energy balance of debris-covered and clean-ice glaciers to the pre-monsoon to monsoon transition potentially accumulate to large mass balance biases when not adequately taken into account in models. It is thus im-portant to also consider how this seasonal transition may evolve under changing climate, and what the consequences for glacier mass balance might be. All future model simulations agree on continued warming during the 21st century over High Mountain Asia, together with a strengthening of elevation dependent warming (Palazzi et al., 2017) and increases in moisture availability.
An analysis on the ensemble estimates of regionally downscaled CMIP5 projections (CORDEX) for the Himalayas (Sanjay et al., 2017) shows that total summer precipitation is projected to increase for 2036-2065 (2066-2095)  and spatial distribution of precipitation across High Mountain Asia (Kadel et al., 2018;Sanjay et al., 2017), which is likely a result of complex and poorly understood drivers of past monsoon changes (Saha et al., 2014;Saha and Ghosh, 2019), coarse resolution of the baseline products and topographic variability in the region (Sanjay et al., 2017). A slight shift towards an 545 earlier monsoon onset of <5 days over the coming century together with an increasing shift towards a later retreat by 5 to 10 days (mid-century) and 10 to 15 days (end-century) might increase the length of the monsoon period, with stronger lengthening in the Eastern Himalaya (Moon and Ha, 2020).
The prospect of warmer temperatures together with increased precipitation would (1) cause a shift in the precipitation partition from snow to rain in the monsoon, resulting in snow cover shifting to higher elevations and increasing total melt; (2) potentially 550 lead to an increase in early spring snowfall, which would delay the onset of ice melt; (3) increase the likelihood of ephemeral monsoonal snow cover but move it to higher elevations, thus leaving more of the lower ablation zones exposed; (4) increase the wetbulb temperature together with humidity to result in a reduction of the solid fraction of precipitation during monsoon.
Overall it is likely that glacier ablation zones will be exposed for longer periods under future climate due to a net decrease of the snow covered duration, with a resulting increase in total ablation. A lengthening of the monsoon into autumn, on the other 555 hand, (Moon and Ha, 2020) would somewhat offset warmer air temperatures with regards to the late-season melt for all glacier types.
Warmer and wetter monsoonal conditions, including increased cloudiness, are likely to result in an overall increase of melt rates on clean-ice and dirty-ice glaciers. This is because (1) they are more directly controlled by R net , which is likely to increase in magnitude; (2) the turbulent fluxes towards cold surfaces are also likely to increase in magnitude, and they tend to 'work 560 against' these types of glaciers. In contrast, the turbulent fluxes 'work for' debris-covered glaciers, and the melt-equalizing effect of debris under monsoon would likely remain in place. With these components summing up to have an overall protective effect on debris-covered glaciers, they are likely to resist the projected changes in the monsoonal summer longer into the future.
Previous studies suggested that the mass balance of DCGs might be less sensitive to climate warming than clean-ice glaciers (e.g. Anderson and Mackintosh, 2012;Wijngaard et al., 2019;Mattson, 2000). Here we confirm this hypothesis and suggest 565 that this is even more the case under monsoonal conditions.
Simplified methods of glacier melt calculation (e.g. relying only on temperature or temperature and shortwave radiation) may integrate some of these processes through calibration (e.g. Ragettli et al., 2016). However, in a study on two clean-ice glaciers (Yala and Mera), Litt et al. (2019) were able to transfer Temperature Index and Enhanced Temperature Index models between sites and years only during pre-monsoon. During monsoon, the transfer failed due to site-specific and inter-annually variable cloudiness.

Future work
By modelling seven different study glaciers across the Himalayan arc, we have been able to identify important monsoon effects on the energy and mass balance of debris-covered and clean-ice glaciers. In the context of the large spatial and inter-annual 575 variability in the climate, including the monsoon strength and timing, our sample size remains small due to data availability.
The timing and quantity of spring snow-cover and monsoonal snowfalls have large effects on the energy and mass balance, and these variables can be highly variable between years. A future analysis could leverage a greater number of complete AWS records, and possibly multi-year records, in order to extend some of our findings. Despite the advanced representation in Tethys-Chloris, the surfaces of glaciers, and especially of debris-covered glaciers, are more complex systems than indicated at the AWS 580 location. Future work should be invested into spatially distributing forcing data and parameters necessary to run energy-balance models in a distributed framework at glacier and catchment scales. Studies employing distributed energy balance models could, for example, test the melt-equalizing effect of debris at the glacier scale and its overall effect on the catchment runoff. This could be combined with generating spatially consistent forcing data using high-resolution weather modelling as in Bonekamp et al. (2019) or Potter et al. (2020), and expanded to larger domains. Realistic representations of the glacier surface, including 585 distributed debris thickness, and supraglacial features, such as ice cliffs and surface ponds, should also be established to provide well-constrained water budgets and improved representations of glacierised environments in land-surface models under present and future environmental conditions.

Conclusions
By modelling the energy and mass balance of seven glaciers in the Himalayas at on-glacier weather stations, we identify and 590 explain the main effects on the energy and mass balance caused by the Indian and East Asian summer monsoons. Heavy cloud cover, liquid precipitation, wind speed, the presence and thickness of debris cover, and elevation shape the energy and mass balances of the Himalayan glaciers during the ablation season. The timing of snow melt-out and the presence of ephemeral monsoonal snowcover play a particularly important role, especially at high elevations. We highlight key pre-monsoon to monsoon changes in energy fluxes, distinct for clean-ice and debris-covered glacier sites: the melt of clean-ice glaciers is primarily 595 radiation-driven at any point during the ablation season, and strongly influenced by albedo variations. The latent heat flux can initially be a small sink and turn into a source during the core monsoon. Debris cover can act in two ways: (1) Once a debris layer of a few centimetres is established, it returns most of the heat it absorbs back to the atmosphere via longwave and turbulent heat fluxes. The sensible heat flux reduces during core monsoon, but the latent heat flux removes energy from the surface due to evaporation of liquid water. The turbulent fluxes readjust in this way to monsoonal temperature, wind and moisture regimes, 600 maintaining a nearly constant melt rate over the entire ablation period. Sensitivity analyses of our energy-balance model shows that these findings hold for large portions of the debris-covered glacier surfaces.
(2) When it is thin, debris amplifies the in-