Modelling ice-cliff backwasting on a debris-covered glacier in the Nepalese Himalaya

Abstract

Ice cliffs have been identified as a reason for higher ablation rates on debris-covered glaciers than are implied by the insulation effects of the debris. This study aims to improve our understanding of cliff backwasting, and the role of radiative fluxes in particular. An energy-balance model is forced with new data gathered in May and October 2013 on Lirung Glacier, Nepalese Himalaya. Observations show substantial variability in melt between cliffs, between locations on any cliff and between seasons. Using a high-resolution digital elevation model we calculate longwave fluxes incident to the cliff from surrounding terrain and include the effect of local shading on shortwave radiation. This is an advance over previous studies, that made simplified assumptions on cliff geometry and radiative fluxes. Measured melt rates varied between 3.25 and 8.6 cm d⁻¹ in May and 0.18 and 1.34 cm d⁻¹ in October. Model results reproduce the strong variability in space and time, suggesting considerable differences in radiative fluxes over one cliff. In October the model fails to reproduce stake readings, probably due to the lack of a refreezing component. Disregarding local topography can lead to overestimation of melt at the point scale by up to ∼9%.

1. Introduction

Debris-covered glaciers are a common feature in high mountain ranges and constitute a significant share of the total glacierized area in the Hindu Kush–Karakoram–Himalaya (Benn and others, 2012; Bolch and others, 2012; Nuimura and others, 2012). Debris cover over a few centimetres thick is generally understood to act as insulation to the ice mass of the glacier, inhibiting ablation (Oestrem, 1959). Recent remote-sensing studies have, however, failed to identify diminished rates of recession compared to debris-free glaciers (Kääb and others, 2012; Basnett and others, 2013; Gardelle and others, 2013). Despite the fact that these studies refer only to a very recent period of data, they have triggered investigations of the behaviour of debris-covered glaciers at scales larger than that of the numerous point-scale numerical studies carried out until recently. An explanation for high ablation rates at the glacier scale has been attributed to the occurrence of ice cliffs and supraglacial lakes on debris-covered glaciers (Gardelle and others, 2012; Pellicciotti and others, 2015). Although cliffs cover only small portions of the total debris-cover area (<2% for cliffs on Lirung Glacier; Sakai and others, 2000), melt rates >5 cm d⁻¹ during the melt season have been reported and could hence contribute a significant amount of mass loss when extended to the scale of the entire glacier (Sakai and others, 1998; Han and others, 2010; Reid and Brock, 2014). The increased melt on sloping cliffs could be caused by a combination of factors: (1) the low albedo of the cliff ice due to very fine debris originating from the surrounding debris on top of and around the cliffs; (2) the exposure to a high amount of solar radiation when cliffs face south (north) in the Northern (Southern) Hemisphere; and (3) high longwave radiation receipts emanating from the surrounding debris, which, unlike debris-free ice, reaches very high temperatures during the day. While each of these assumptions has been variously suggested in previous studies, none has been tested systematically or in combination with the others.

Debris-covered glaciers are characterized by rough topography (Benn and others, 2012; Pellicciotti and others, 2015). Cliffs also exhibit complex geometry, with strong variability in extension, size, aspect, slope and shape of the vertical profiles (Immerzeel and others, 2014). Measurements of backwasting at single locations on cliffs (inferred from measurements at stakes or at the cliff top) have been useful to provide quantitative estimates of backwasting rates and to qualitatively suggest possible dominant controls (Sakai and others, 1998; Han and others, 2010; Reid and Brock, 2014). However, only a model able to calculate the energy fluxes resulting from the complex cliff topography can properly quantify backwasting rates and their controls. Sakai and others (1998) proposed a basic model on Lirung Glacier, which calculates shortwave and longwave radiation from the sky, as well as turbulent fluxes. In a later publication the shading from the surrounding topography was introduced, together with longwave radiation from the debris opposite the cliff (Sakai and others, 2002). So-called view factors were introduced to define the extent to which the cliff is exposed to the sky and opposite debris. The view factors were derived from a topographical map with contours at 10 m intervals. Han and others (2010) applied this model to a glacier in the Chinese Tien Shan. A low-resolution digital elevation model (DEM) was used to account for the role of the surrounding topography. They also provided a refined calculation of the turbulent fluxes, using the aerodynamic bulk method, while the earlier method was based on empirical constants. Building on Han and others’ (2010) model, Reid and Brock (2014) implemented further improvements to study cliffs melting on Miage glacier, in the European Alps, and used their model to quantify the total cliff contribution at the glacier scale. They improved the calculation of longwave radiation by using the debris surface temperature instead of the air temperature employed by Han and others (2010), and used a higher-resolution DEM to account for local topography. They investigated energy fluxes and melt at different positions along the cliff height and accounted for debris mounds opposite the cliffs that were higher than the stake position itself.

In this paper, we use new data measured on Lirung Glacier in the Nepalese Himalaya, both on the glacier and directly at the cliff’s surface (including radiative fluxes parallel to the cliff surface, a novel measurement). A high-resolution DEM of the glacier was used to provide a much more detailed representation of the complex surface topography than in previous studies. In so doing, we address the following three aims: (1) to analyse the new datasets to gain insights about the variability of observed melt rates and radiative fluxes at the cliff scale, during the day and over the strong seasonality typical of the monsoon climate of the Langtang valley; (2) to improve the model based on these new insights, in particular improve the calculations of radiative fluxes; and (3) to use the improved model to understand backwasting processes, quantify the importance of single fluxes and analyse their spatial variability.

In contrast to earlier studies we have multiple measurements available during the day and can hence check the model’s performance at the sub-daily scale. Our data also include measurements from two distinct seasons, allowing us to judge the model performance with variable climatic input.

2. Study Site and Field Measurements

2.1. Study site

This paper is based on measurements on Lirung Glacier in the Nepalese Himalaya (28°13′57″ N, 85°33′43″ E) over one main melt season, May–October 2013. The glacier lies in the Upper Langtang catchment, which extends over an area of 350 km², ∼30% of which is glacierized (Shiraiwa and Yamada, 1991; Ragettli and others, 2015). The largest glaciers are debris-covered on their tongue, representing ∼25% of the total glacierized area (Pellicciotti and others, 2015; Ragettli and others, 2015). A comprehensive glacier inventory of the catchment lists 14 of the 72 glaciers as debris-covered (Shiraiwa and Yamada, 1991). Lirung Glacier has the greatest elevation range (∼4000–7132 m a.s.l.), with the debris cover extending to 4400 m a.s.l. Its total area is >6 km², with >20% covered in debris (Shiraiwa and Yamada, 1991). The debris-covered tongue has a length of 3.5 km and is on average 500 m wide (Immerzeel and others, 2014). The debris thickness is very heterogeneous, ranging from boulders several metres in height to fine silt and sand. Over most of the tongue the cover is >50 cm deep (Ragettli and others, 2015). In 2013 eight ice cliffs were identified, on two of which measurements were undertaken (Fig. 1). During the 2013 melt season, extensive fieldwork was conducted on the glacier and in the valley to understand the response of debris-covered glaciers to climate.

Fig. 1. Map of Lirung Glacier in the Langtang catchment (indicated as the shaded area in the inset map in the top right corner). The black curve indicates the tongue glacier border, while the accumulation area in the melt period is shaded. The white box indicates the area covered by unmanned aerial vehicle (UAV) flights. The slope map extracted from the high-resolution UAV digital elevation model (DEM) is shown enlarged, with the two cliffs investigated marked as C1 and C2.

2.2. Meteorological data

Meteorological data were collected on Lirung Glacier in two field campaigns in May and October 2013. An automatic weather station (AWS Lirung; 4076 m a.s.l.) was installed on the glacier in May (Fig. 1; Table 1) and recorded conditions over the monsoon season until October 2013. Additional meteorological data were available from one off-glacier station to the south of the tongue (AWS Kyanjing; 3857 m a.s.l.; Fig. 1), which has been operating since May 2012. AWS Lirung was installed between the two cliffs (1 and 2; Fig. 1) investigated in this work.

Table 1. Locations and dates of the measurements used in this study. AWS Lirung and AWS Kyanjing indicate the two automatic weather stations on- and off-glacier, respectively. CNR1 is the net radiometer installed on cliff 1 parallel to the cliff surface. Cliffs 1 and 2 are the two cliffs that were monitored for this study, and on which stake readings were taken in May and October 2013

‘May’ and ‘October’ refer to the periods when stake readings were taken in the field, given in Table 1. Premonsoon is defined from the beginning of measurements at AWS Lirung (9 May) until the beginning of the monsoon (14 June). Post-monsoon is defined from the end of the monsoon (20 September) until the end of the measurements (30 November, including the data extension discussed in Section 2.2.2). When referring to the monsoon period we discuss data from 15 June to 19 September. Our division of seasons follows that of Immerzeel and others (2014), except for the end of the monsoon season, which they fixed as 30 September, but we set as 19 September. A clear separation of seasons is not unambiguous, but our definition of the monsoon seems to better define the period with reduced diurnal temperature variability and much higher precipitation rates.

At AWS Lirung the variables shown in Table 2 were measured at 5 min intervals. Incoming longwave radiation was not measured at AWS Lirung and was therefore modelled (Section 2.2.1). A network of temperature sensors with an incorporated logger (referred to as T-Loggers) were installed on the glacier tongue (Petersen and Pellicciotti, 2011). For this study, we use the measurements from a T-Logger located next to AWS Lirung.

Table 2. Characteristics of the sensors installed on AWS Lirung and on the cliff. I↓↑ and L↓↑ are shortwave and longwave radiation, u wind speed, w _d wind direction, T _a and T _s air and surface temperature and RH_a relative humidity. C are sensors at the cliff. K&Z refers to Kipp & Zonen

Due to a technical failure, AWS Lirung provided no measurements after 23 October. We therefore extended the record until 21 November to match the measurements at the cliffs (Table 1), as explained in Section 2.2.2.

A Kipp & Zonen CNR1 net-radiometer (Fig. 2a; Table 2) was placed on cliff 1 where stake 1.1 was located in October (Fig. 3). Because of the rapid backwasting and continuous rockfall, the sensor cannot be maintained for a long period on the cliff, and is only used to validate the modelled fluxes of shortwave radiation (calculated as described in Section 3.1.1).

Fig. 2. (a) The net-radiometer (Kipp & Zonen CNR1) deployed at the cliff surface at a height of 2 m, measuring incoming and reflected shortwave radiation and incoming and outgoing longwave radiation. (b) The AWS deployed parallel to the debris surface.

Fig. 3. Photographs of cliffs 1 and 2 in May (top) and October (bottom) with the respective stake locations, each time taken from a similar position. Notice that, on cliff 1, stake 1.3 in May is at the same location as stake 1.1 in October. On cliff 2, stakes 2.1–2.3 are at approximately the same locations in both seasons.

Because of the assumed importance of low albedo to the cliff ablation, we carried out one measurement of albedo at each stake location on the cliff in October using a handheld luxmeter (Labfacility – LX101). Observed albedo was quite low, with values between 0.04 and 0.29 (Table 3), in accordance with observations of the cliff surface in the field. The ice was covered by a fine dust cover, which becomes part of the surface, as the slurry mix of fine debris and water remains on the cliff surface even at slopes >30°. During both seasons it was observed that meltwater would refreeze on the cliff. The albedo measurements with the CNR1 show the same range (0.1–0.3). These values are lower than the constant value of 0.37 used by Han and others (2010) and correspond more with the constant value of 0.09 used by Reid and Brock (2014). We used the observed albedo as initial estimates for the model optimization described in Section 3.2.

Table 3. Characteristics of the ablation stakes at both cliffs 1 and 2 in May and October and corresponding mean daily melt rates. Albedo was measured in October with a handheld luxmeter. ‘Distance’ refers to the distance of the stake location from top of the cliff

2.2.1. Modelling incoming longwave radiation at AWS Lirung

Incoming longwave radiation for the required period of time was modelled using the Stefan–Boltzmann relationship as

(1)

where ∊ _eff is the effective emissivity of the sky, σ the Stefan–Boltzmann constant and T _a the air temperature (K), measured at the AWS. The sky emissivity was determined with distinct parameterizations for clear sky (Dilley and O’Brien, 1998) and cloudy sky (Unsworth and Monteith, 1975), a combination which has been found to provide the best results in a study comparing different models of longwave radiation (Juszak and Pellicciotti, 2013). The model was calibrated (2013) and validated (2012) with measurements of longwave radiation at AWS Kyanjing for each season separately. Model performance was very good during the dry season, but decreased in the monsoon, when cloud cover was present nearly every day, and the daily cycle could not be reproduced well. The root-mean-square error (RMSE) between modelled and observed incoming long-wave radiation was ∼27 W m⁻², which represents an error of ∼10% for this flux. None of the available models for incoming longwave radiation has been developed for the monsoon climate of South Asia or the Himalaya, hence the use of such a model may not be entirely suitable for the region and future work could be devoted to the development of a parameterization specific for this climate.

2.2.2. Temporal extension of datasets of AWS Lirung in October

To extend the series of observations at AWS Lirung after the period when the AWS failed on 23 October, we used distinct techniques for different meteorological input variables. Clear-sky incoming shortwave radiation was modelled using a non-parametric model (Bird and Hulstrom, 1981; Iqbal, 1983) implemented in the glacio-hydrological model TOPKAPI (Ragettli and others, 2015). Terrain parameters and solar position and the interaction between the solar radiation and the topography were derived with the vector algebra approach proposed by Corripio (2003). The cloud transmittance factors needed to derive all-sky solar radiation were calculated using the approach of Pellicciotti and others (2011), with data from AWS Kyanjing. The entire model is described by Pellicciotti and others (2011) and was applied in simulations of solar radiation in the entire upper Langtang catchment by Ragettli and others (2015).

Air temperature and relative humidity were extrapolated from AWS Kyanjing using lapse rates calculated from the available post-monsoon time series, when the AWS still functioned (20 September–23 October). Since the lapse rate changes during the day when the debris heats up, a lapse rate was determined for each hour of the day. We tested our approach for the post-monsoon period where data were still available. Agreement between lapsed and measured values is very good (with a correlation coefficient, R ², of 0.97 for air temperature and 0.95 for relative humidity). A linear relationship was obtained for wind speed between the two AWSs (R ² = 0.75 for the post-monsoon season) and was used to reconstruct wind speed at AWS Lirung. Surface temperature was derived from air temperature using a linear regression (Foster and others, 2012) for each hour of the day. Comparison of measured with modelled surface temperatures before station failure resulted in R ² = 0.93.

2.3. Topographic data

The two cliffs investigated in this work are shown in Figure 3 and their cross sections in Figure 4. The two cliffs were located close to the AWS (cliff 1 ∼100 m south and cliff 2 ∼200 m north). Slope and aspect measured manually with an inclinometer at the location of each stake are listed in Table 3. Cliff 1 had a maximum height of 13 m and extended over a width of ∼50 m. The undercut and the opening of an englacial conduit suggest the former presence of an ephemeral lake below the cliff. As the cliff backwasted, the western part of cliff 2 (where stake 2.4 was located in May) was buried by debris from above, as were previously exposed areas on the eastern side (not visible). The cliff’s width therefore decreased from 180 m in May to <100 m in October. Cliff 2 had a maximum height of 30 m and was associated with a lake that drains and fills periodically during the year, at least partially from the far side of the cliff. Its lower section was partially overhanging, where most probably the lake has caused accelerated retreat.

Fig. 4. Cliffs 1 (C1) and 2 (C2) in May and October. The stake numbering corresponds to that of Figure 3.

During both field campaigns high-resolution DEMs (20 cm, hereafter referred to as UAV DEM) were obtained from multiple flights over the glacier lower tongue with an unmanned aerial vehicle (UAV) (Immerzeel and others, 2014). To smooth out small errors we resampled these DEMs to 1 m to estimate sky- and debris-view factors. A second DEM was available for the entire catchment, the ASTER GDEM of 30 m resolution (NASA, 2001). It was used to calculate the shading of the incoming solar radiation.

2.4. Melt observations

Bamboo stakes were drilled into the ice at different locations over both cliffs (Fig. 3) and used to determine backwasting rates at the cliff surface. Melt rates were measured perpendicular to the ice surface. Earlier studies (Sakai and others, 1998; Han and others, 2010) measured horizontal backwasting at the crest of the cliff, while Reid and Brock (2014) also measured perpendicular melt.

Measurements were taken daily and, in some cases, twice a day. Slope and aspect were recorded at each cliff site where stakes were located (Fig. 4). In May, the stakes on both cliffs had similar slope values (40–43°; Table 3) with the exception of stake 2.3 on cliff 2, which was set up at a steeper location (51°; Table 3; Fig. 4). In October, slopes were more variable, with a very steep slope at stake 1.5 on cliff 1 and a less steep slope at stake 1.4, on the same cliff. Where possible, we tried to install the stakes at locations with different aspects. Slope and aspect values obtained in the field were compared with those obtained from the UAV DEM, resulting in very similar values.

Observations are shown in Figure 5. Melt rates are considerably higher in May than in October (Fig. 5; Table 3). This can be explained, to a large degree, by differences in shortwave radiation and temperature between the two seasons, with air temperatures higher in May (Fig. 5). Differences between the two cliffs and locations on the same cliff are considerable in the same period. In May, there is large variability especially on cliff 2, with stake 2.4 facing east showing much higher melt than the other stakes. Here melt is highest at stake 2.2 in the centre of the cliff, with the highest and lowest having similar melt rates (Fig. 5a; Table 3). At cliff 1, stake 1.1 has the highest melt (Fig. 5a; Table 3). In October, melt rates on cliff 2 are all similar (probably due to the proximity of the stakes), but variability is important on cliff 1. The two uppermost stakes of the eastern transect (stakes 1.3 and 1.4) have much higher melt rates than those at the other locations, probably due to higher incoming shortwave radiation. Energy fluxes in the two seasons at each stake are analysed in detail in Section 4.2, to explain the high variability evident in the observations.

Fig. 5. Melt readings in (a) May and (b) October at cliffs 1 and 2. Hourly air temperature measurements at the TLogger close to cliff 1 are also plotted for reference. Notice that the length of the horizontal axis is the same in the two panels, to allow comparison of melt rates.

3. Methods

3.1. Cliff energy-balance model

The basis of the model used in this work was presented by Han and others (2010) and further developed by Reid and Brock (2014). Here we briefly summarize the existing model components and provide a more detailed description for parts where improvements to the earlier versions were implemented.

The basic energy-balance model equation for a unit horizontal area is

(2)

where Q _m is heat for vertical ice melt, I _n and L _n are net shortwave and longwave radiative fluxes and H and LE are turbulent sensible and latent heat fluxes (W m⁻²). The conductive heat in the ice and the heat from precipitation are neglected (Reid and Brock, 2014).

The available energy for vertical melt, Q _m, is converted into melt perpendicular to the surface (cm) as

(3)

where Δt is the time step (in this case 1 hour = 3600 s), ρ _i the density of ice (kg m⁻³) and L _f the latent heat of fusion of water (334 kJ kg⁻¹). All fluxes are calculated perpendicular to the surface. The model runs at an hourly resolution and fluxes are calculated as explained below.

3.1.1. Shortwave radiation

The net shortwave radiation incident to a unit of sloped area is calculated as

(4)

where I _s is direct incoming, D _s diffuse and D _t shortwave radiation reflected from the terrain surrounding the cliff (W m⁻²) and α _i is ice albedo. As not all the terrain surrounding the cliff may receive sunlight and not all of its area may point towards the cliff, the value for D _t has to be considered a maximum estimate.

Incoming shortwave radiation is thus calculated as the sum of three terms, all based on the measurements of shortwave radiation at AWS Lirung. I _s is calculated as by Han and others (2010), based on an approach by Garnier and Ohmura (1968), by which the estimated direct normal irradiance measured at AWS Lirung is transferred to the cliff surface using the zenith angle, horizon angle in the direction of the sun, solar declination, latitude and azimuth.

Diffuse shortwave radiation is calculated as

(5)

where V _s is the sky-view factor (described in Section 3.1.2), k _d the diffuse fraction (calculated according to Reindl and others, 1990) and I ₀ is shortwave radiation measured at the AWS (W m⁻²).

Shortwave radiation received from reflection by the terrain surrounding the cliff is calculated as

(6)

where α _t is the albedo from the debris surface, found to be, on average, 0.15 from measurements at the AWS, which is similar to that found by Reid and Brock (2014) but much lower than that of Han and others (2010) (0.12 and 0.24, respectively). It is assumed here that all debris slopes face the ice cliff in a parallel way, while from some of these slopes radiation may be reflected elsewhere. The modelled value of D _t is therefore a maximum estimation.

In this paper we model incoming shortwave radiation as reaching the cliff face perpendicularly. Modelled values do correspond well with measurements at the cliff face during the short period of measurements with the CNR1 (Fig. 6). A slight shift of the modelled flux in comparison with the measured flux is visible in the midday values (with the model simulating the peak radiation ∼1 hour earlier than observations). Beside this, however, the agreement is remarkable and shows the drastic reduction in the flux incident to the sloped north-facing cliff surface, in comparison with the values measured horizontally at the glacier (AWS Lirung). The horizon angle, h, equally necessary for the calculation, is calculated with an enhanced approach described together with the calculations of the longwave flux in Section 3.1.2.

Fig. 6. Incoming shortwave radiation modelled at stake 1.1 on cliff 1, I _s mod, compared with the measured values from the CNR1, I ₀ meas. Also shown are the measurements at AWS Lirung, I ₀ AWS, that are used to estimate radiation on a slope from the AWS measurements.

Ice albedo, α _i, is assumed to be constant in time at each stake location. Like Han and others (2010) and Reid and Brock (2014), we used a constant value for albedo, but we optimized it with a Monte Carlo approach at each stake (Section 3.2). Han and others used a constant ice albedo of 0.37, while Reid and Brock used a constant mean value of 0.094 obtained from point measurements. In contrast to the earlier studies, here albedo measurements were available in October for each stake location and at one location over a whole day, providing a range of possible values for the modelled albedo. These values (Table 3) were therefore used to define the range of plausible values for the Monte Carlo optimization.

3.1.2. Longwave radiation

Net longwave radiation is calculated as

(7)

where V _l is the sky-view factor for longwave radiation, L _in is the incoming longwave radiation from the atmosphere, assumed equal to that measured at the AWS, L _d the longwave radiation from the surrounding debris and L _o the outgoing longwave radiation from the ice-cliff surface (all in Wm⁻²). Longwave radiation is assumed to be the same on the slope and on the horizontal plane, hence slope and aspect are not considered.

Longwave radiation emitted by the terrain is calculated with the Stefan–Boltzmann law as

(8)

with T _d the measured temperature of the debris (K), ∊ _d the debris emissivity equal to 0.95 and V _d the debris-view factor. Outgoing longwave radiation is determined using the Stefan–Boltzmann relation, with the ice surface temperature, T _i, taken equal to 0°C and the ice emissivity, ∊ _i, equal to 0.97, resulting in a constant outgoing flux over time. Both emissivity values were chosen following Reid and Brock (2014). The assumption T _i = 0°C is a simplification, as it is likely that the ice surface reaches higher temperatures because of the debris on top but could also sink below the freezing point at night. More measurements would help evaluate how sound this assumption is.

Calculation of the portion of the sky and of the debris that are seen from each point of the cliff is crucial for an accurate estimate of both shortwave and longwave radiation. The sky-view factor is the percentage of the overhead area from a point on the cliff that is open to the atmosphere and from which the cliff receives downwelling atmospheric radiation (Fig. 7a). The debris-view factor is the percentage of the view that is obstructed by debris, and hence contributes longwave radiation emitted by the debris-covered glacier surface as well as reflected shortwave radiation from debris (Fig. 7b).

Fig. 7. Determination of the sky-view and debris-view factors for different cases. (a) and (b) show the calculation of the individual horizon angle for the sky-view factor and the debris-view angle for the debris-view factor for a single direction; these views are then aggregated over 360° to determine one sky-view, V _s, and one debris-view factor, V _d, for each location on the cliff. (c) If the local topography (i.e. a debris mound) shades the mountains in the back, the sky-view factor, V _s (light shading), is the same for both shortwave and longwave radiation; the debris-view factor, V _d (darker shading) is determined by the debris mound facing the cliff. (d) When the topography in the distance is visible from the location on the cliff, the sky-view factor, V _s (solid curve and shaded), for shortwave radiation is determined by the mountains, while the sky-view factor for calculation of the longwave radiation from the atmosphere, V _l (dashed curve), is determined by the opposite debris mound.

Reid and Brock (2014) took into account the effect of the surrounding topography by assuming the terrain below the location of the stake on the cliff was opposite a debris mound and the area above exposed to the sky. This way the fraction of the hemisphere exposed to debris could be taken into account for longwave calculations. However, this approach did not take the actual height of the opposite debris mound into account. To account for debris mounds that are higher than the actual location on the cliff, Reid and Brock introduced a parameter representing an additional amount to the debris-view factor, V _d, which they determined through optimization of the model against stake readings. The value was found to vary between 0 and 0.1 for different cliffs. The parameter calculated in this way is specific for individual cliffs and determined through calibration, and might incorporate errors in both stake readings and ablation model.

We therefore improve on this assumption by implementing a model that requires information on the visible horizon from the location on the cliff surface. Horizon angles were calculated from each point of the DEM. The model looks at 360 1° steps in full circle, and hence considers all obstacles that obstruct the view. It takes into account the distance and the elevational difference to each pixel and then calculates the angle to all obstacles within the 1° slice. The maximum for each of the 360 steps is selected and then integrated to determine the sky-view and debris-view factors, respectively (the approach is sketched in Fig. 7 and described by Dozier and Frew, 1990).

The sky-view factor for the incoming longwave radiation from the atmosphere, V _l, is determined from the topography close to the cliff (Fig. 7b). This includes the assumption that longwave radiation from the mountain slopes further away can be treated as atmospheric longwave radiation. Plüss and Ohmura (1997) proposed that the distance to the emitting surface has a large relevance for distances <1 km, in which case additional radiation, depending on surface and air temperature, should be added. They estimated additional radiation of up to 60 W m⁻² in the case of unobstructed skies. The slopes at the study site are relatively close, hence this may play a role. With the data available we are, at this point, not able to make a more accurate statement about the contribution of distance to the emitting surface.

To test the impact of the UAV DEM on model calculations, we also performed model runs with a resampled 30 m DEM.

3.1.3. Turbulent fluxes

Turbulent sensible, H, and latent heat, LE, fluxes are determined with the bulk aerodynamic method used by Han and others (2010) and Reid and Brock (2014). Stability corrections based on the Richardson number (Brock and others, 2007, 2010) are applied. As Reid and Brock noted, surface roughness, z ₀, may need more attention. Han and others used an approach based on the work of Lettau (1969), that uses height and slope of the cliff, as well as wind direction relative to the cliff surface, to calculate z ₀. This implicitly assumes the cliff to be a roughness element and disregards the roughness of the cliff surface itself, where the energy exchange actually takes place. This solution seems problematic, since cliffs cannot be assumed to be a uniform obstacle with a constant height or slope, and wind directions towards the cliff cannot be assumed to be the same as at the AWS. Calculations based on the assumption of Han and others with different realistic values of cliff height, slope, aspect and wind direction lead to maximum values of z ₀ > 1 m. In contrast, literature values for surface roughness rarely exceed 1 cm, and for glacier ice are generally in the range 0.1–5 mm (Brock and others, 2006). The logarithmic model (Andreas, 1987) used to derive the roughness length of temperature and vapour pressure may not be applicable at such a scale. Conversely, if we assume that the main flow direction is parallel to the cliff surface, then the exchange of heat through turbulent fluxes will occur perpendicular to the surface and the logarithmic profile will be valid. We therefore decided to use a commonly used literature value for z ₀ as an initial value (1 mm), as proposed by Pellicciotti and others (2005), which is close to the ideal value found by Reid and Brock (2014) and within the range presented by Brock and others (2006).

A possible source of error is the fact that wind and humidity measurements from the horizontal plane at the AWS are used, which might be different from those on a slope. However, without detailed datasets from the cliff surface, we retain this assumption.

3.1.4. Direction of the fluxes

An important change here to the treatment of Reid and Brock (2014) is the trigonometrical adjustment of the fluxes. Reid and Brock multiplied radiative fluxes by the cosine of β (where β is the ice-cliff slope (°)) to determine the equivalent vertical melt. Turbulent fluxes were calculated on a horizontal plane, hence assumed to cause vertical melt. Subsequently, the resulting energy for melt, Q _m (Eqn (3)), was multiplied by the cosine of β, to calculate the equivalent melt in the perpendicular direction. Since radiative fluxes are modelled as acting perpendicular to the slope and hence result in melt parallel to the slope, the respective adjustment for slope is superfluous and can be discarded. Since the logarithmic profile for wind also holds on a slope, the turbulent fluxes can equally be assumed to act perpendicular to the slope.

3.2. Monte Carlo analysis and model optimization

Model parameters that were not measured but can be constrained within a range of plausible values were optimized using a Monte Carlo approach (as in Ragettli and Pellicciotti, 2012; Heynen and others, 2013). The parameters optimized in this way are the albedo of ice (_i, for net shortwave radiation), albedo of the debris (_t, for shortwave radiation reflected from the surrounding terrain), the emissivities of ice and albedo (∊ _i and ∊ _t, for outgoing and debris longwave radiation, respectively) and the surface roughness of the ice cliff (z ₀, for turbulent fluxes). One million model runs with randomly selected parameter combinations from an initial prescribed range were carried out. The initial parameter ranges are listed in Table 4. The RMSE of the model output against the stake readings was determined for each parameter combination and the combination with the best fit was used to calculate melt and energy fluxes at each stake. Optimization was carried out separately for May and October and separately for each stake. No stake readings were available during the monsoon.

Table 4. Optimal parameters obtained with the Monte Carlo analysis in May (top) and October (bottom). Stake numbers are as in Figure 3. α _i and α _t are albedo and ∊ _i and ∊ _t emissivity values for ice and debris, respectively. z ₀ is the surface roughness. In the bottom rows the RMSE is given for the model performance. In parentheses are initial parameter ranges and the initial value. In October the measured values for α _i were used as the initial value (Table 3)

The initial value for ice albedo was chosen as 0.15 and the possible range limited to 0.01–0.35, corresponding to measured values from the ice cliffs. The initial value for debris albedo was taken to be 0.15, corresponding to a mean value obtained from radiation measurements at the AWS, ranging between 0.1 and 0.2. Here the range was limited to 0.01–0.3. Emissivity values were 0.97 for ice and 0.95 for debris (Reid and Brock, 2014), with ranges 0.95–0.99 and 0.93–0.97, respectively. Surface roughness was initially set to 1 mm and varied between 1 and 5 mm (Brock and others, 2006).

3.3. Sensitivity analysis

Sensitivity analysis was conducted in two steps. First, we investigated the Monte Carlo runs, which reflect a global sensitivity of the model to the parameters used in the optimization, since all parameters are changed at the same time. The 100 best runs from the Monte Carlo analysis (i.e. resulting in the lowest RMSE) were selected and the parameters from these runs plotted for each stake in both seasons, to evaluate their spread within the best realizations (e.g. Finger and others, 2011). Parameters with a large spread are those that the model is not sensitive to, as any of those values lead to a high model performance. On the contrary, parameters with a small spread (indicated by the confidence interval) are those the model is sensitive to.

In a second step, a well-established one-at-a-time (OAT) sensitivity analysis was conducted for a selection of model variables. The method is described by McCuen (1973) and has been used for sensitivity analysis of glacier models (Ragettli and Pellicciotti, 2012; Heynen and others, 2013). We test the model sensitivity to variables that might be affected by errors. Some of these were measured in the field but may vary on relatively small scales (β, aspect, α _i). Others were not measured at the cliff surface but at AWS Lirung over debris, and we assume them to be similar over the cliff (u, T _a). The view factors are dependent on the accuracy of the DEM. Since the use of this high-resolution DEM was an improvement over earlier studies, the sensitivity of the model to the view factors is also investigated. Longwave radiation was not measured on the glacier but modelled at AWS Lirung. This may introduce further uncertainties. Each variable was individually varied in a range of ±20% in 5% increments. This change was compared with the change in mean daily melt. The slope of the fitted curve around the optimal provides the sensitivity (Ragettli and Pellicciotti, 2012; Heynen and others, 2013). A steeper slope corresponds to larger sensitivity of the model to the variable. A negative (positive) slope corresponds to an increase (decrease) of melt in correspondence to a decrease in the variable. This OAT method can only be applied if the parameters are optimal and their interdependence small.

4. Results

4.1. Parameter optimization and optimal model outputs

Figure 8 shows Monte Carlo simulations for May for a stake where a global minimum was found (stake 1.1 in May) and one where it is not reached (stake 2.4 in May). At four of seven stakes in May the optimization was successful (Fig. 8). However it was only at stake 2.4 where the model completely failed to reproduce observed melt; at the two other stakes (1.3 and 2.2) the model is very close. This discrepancy can be due to either a model failure or an error in the field measurement at stake 2.4 (slope, aspect and melt readings might be affected by errors). The reading at stake 2.4 is much higher than all the others, which supports it being a measuring error.

Fig. 8. One million model realizations with different randomized parameter combinations at two stake locations in May. The triangle denotes the result with the initial parameter setting, the diamond the optimized result. In red are the 100 best results according to the RMSE. Model optimization worked where a global minimum was reached (as at stake 1.1) and did not work where it failed to reach a minimum (as at stake 2.4).

Table 4 shows the optimized parameters at all stakes. The results point to considerably higher albedo values at cliff 1 than at cliff 2 in May. This corresponds to visual observations, since cliff 2 had a thicker dust cover and black stripes covering the whole cliff face (Fig. 3). An apparent exception at cliff 1 is stake 1.5 in October, which was located below an overhanging piece of ice and showed different properties to all other stakes.

Cumulative melt calculated with the initial and optimized model parameters is shown in Figure 9 for May, and for two stakes in October in Figure 10. In May, the optimized model can reproduce ablation readings accurately at all stakes, except stake 2.4 on cliff 2. The initial parameter sets overestimate melt at cliff 1 but underestimate it at cliff 2, most likely because of different albedo values.

Fig. 9. Cumulative melt calculated with the initial and optimal model parameters at all stakes in May. Stake positions are shown in Figure 3. Note that stake 2.4 has a different vertical axis scale.

Fig. 10. Cumulative melt calculated with the initial and optimal model parameters at two example stakes in October. Stake positions are shown in Figure 3. Stake 1.2 is an example where the model fails due to the unconsidered refreezing process. Stake 1.3 is located at the top of the cliff where refreezing may play a smaller role.

In October the model overestimates measured melt considerably and at all stakes except 1.3, 1.4 and, to a lesser extent, 2.2 (Fig. 10). In agreement with the lower temperature and solar radiation, the model indeed simulates much lower average melt values than in May. However, stake readings are very low at a majority of stakes (Fig. 10; Table 3). The optimization does not reach a global minimum for any but two stakes and we assume that a process is at play that is currently not reproduced by the model. From visual observations in the field, we suggest that this might be refreezing of meltwater at the cliff surface, and investigate this in Section 5. The model parameter optimization was thus carried out only against stake readings within the first two days, for which sub-daily readings were available and where a cumulative error may be negligible. It is encouraging that the model can reproduce the daily cycle of melt, which is clear from the fit with the initial stake measurements (Fig. 10). One earlier study has collected sub-daily measurements and could show agreement with a simple melt model (Benn and others, 2001). The reproduction of the daily melt cycle was, however, not discussed in the earlier modelling studies (Han and others, 2010; Reid and Brock, 2014) and we believe the model is well able to perform at this temporal scale.

The model performance is expressed by the RMSE determined between the stake reading and the model output (Table 4). While the RMSE for the stake location where the model failed in May (stake 2.4) is high, it is on average 1.23 cm for the other stakes, and 0.29 cm for the initial stake readings in October. The RMSE amounts to ∼2% of the total melt in May.

4.2. Energy fluxes at the cliff surface

Figure 11 shows the radiative and turbulent fluxes calculated with the optimal parameters at each stake. A strong variability is evident between the two seasons, between the two cliffs in each season and between the locations on the same cliff.

Fig. 11. Mean hourly energy fluxes at all stake locations during the period of stake readings in the pre-monsoon and post-monsoon seasons. L _in − L _o is the difference between incoming atmospheric longwave radiation and outgoing longwave radiation from the ice. L _d is longwave radiation emitted from the surrounding terrain. I _s is direct shortwave radiation from the sky, D _s and D _t diffuse and terrestrial shortwave radiation. H and LE are turbulent fluxes (sensible and latent heat).

Incoming shortwave radiation is a major source of energy in May, with a relatively important contribution from its diffuse component (Table 5; Fig. 11). However, it varies between the stakes. On cliff 1 in May, direct incoming shortwave radiation is much higher at stake 1.3 than at stakes 1.1 and 1.2, as stake 1.3 is more exposed to the sky because of its central position on the cliff. At stake 1.1, however, the shortwave radiation from the terrain is much higher, resulting in higher values of total incoming shortwave radiation. In the same season, values of radiative fluxes on cliff 2 are similar at stakes 2.1 to 2.3, but considerably different at stake 2.4 both in magnitude and distribution. The highest amount of direct incoming shortwave radiation is received by stake 2.1, which is located on top of the cliff and hence exposed longest to the sun. This amount decreases consistently for the lower stakes (stakes 2.2 and 2.3). At stake 2.4 it is conspicuously low, which results in very low melt rates overall. Note that here the model considerably underestimates observed melt. Diffuse shortwave radiation and shortwave radiation from the terrain are of the same order of magnitude as at the other stakes. A large part of the shortwave radiation flux (in the afternoon hours and at nearly all stakes during the whole day in October) is made up of its diffuse component and shortwave radiation reflected from the opposite terrain. The peak of direct shortwave radiation is in May around 11:00, after which it drops relatively quickly to zero just before 16:00. This corresponds to observations from the field, where even though the glacier surface was still illuminated, ice cliffs facing north were already shaded. It is also apparent that shortwave radiation from the terrain increases for stakes that are located further down the cliff face (stakes 1.3 and 2.3) or face more debris because of their aspect (stake 2.4). Direct shortwave radiation at the north-facing stakes (all except stake 2.4) increases sharply once the sun has reached the overhead position over the cliff itself.

Table 5. Mean values of energy fluxes at each stake during the measurement period in May (8–20 May; top) and October (23 October–21 November; bottom). I _s is direct incoming and D _s diffuse shortwave radiation, D _t shortwave radiation reflected from the surrounding terrain and I _ref reflected from the ice. L _in is incoming longwave radiation from the atmosphere, L _d from the surrounding terrain and L _o outgoing longwave radiation. LE and H are latent and sensible heat fluxes

Longwave radiation emitted by the surrounding debris and reaching the cliff is a major flux. It is able to compensate the net longwave flux, incoming from the atmosphere and outgoing from the cliff (L _in − L _o). This is a striking result and the first quantification of the role played by the surrounding debris in the energy balance. The longwave flux from the terrain has a diurnal cycle associated with debris temperature and varies from stake to stake, but remains a major source of energy, of the same order of magnitude as the incoming shortwave radiation. It is evidently a key flux in the energy balance of the cliffs. Longwave radiation from the debris increases from top to bottom of the cliff (compare stake 1.1 at the top with stake 1.2 at the bottom in May; stakes 2.1 and 2.3 in May; stakes 1.3 to 1.5 in October in Table 5), since lower stakes are more exposed to the surrounding debris.

If the slopes of the surrounding mountains were taken into account in addition to the longwave radiation from debris, its contribution would probably rise considerably, as Plüss and Ohmura (1997) provide an estimate of up to 60 W m⁻² for the case of unobstructed skies. At the same time however, the contribution of longwave radiation from the atmosphere would decrease, as the overhead fraction exposed to the sky would decrease.

Turbulent fluxes play a minor role in both dry seasons. Sensible heat is present all day, while latent heat only contributes to the energy balance in the post-monsoon season.

In October, when the sun is much lower, incoming radiative fluxes are considerably smaller. Direct shortwave radiation hardly ever reaches the cliff surface at some stakes, and shading plays a major role. The close surrounding topography or the cliff itself may shade the ice surface from the sun for a long time given the lower sun angles. For example, stake 1.5 is shielded from the sun due to its location at the bottom of the cliff with a slope of >70° and located below an ice roof receives no direct sunlight at all. This corresponds to observations from the field. In contrast, stakes 1.3 and 1.4, which are located on the upper part of a cliff face with a relatively shallow slope and high exposure to the sky, receive high values of radiation for a very short period. At these two stakes we observed the highest melt rates during October, while stake 1.5 situated only a few metres below had a much lower mean daily melt rate (Table 5). The stakes at cliff 2 are equally shaded by the cliff itself, as the cliff is very high and a debris mound rises just to the south of the ice face.

The balance of atmospheric longwave radiation minus outgoing longwave from the cliff (L _in − L _o) is more negative in the post-monsoon season, due to increased outgoing radiation (higher ice emissivity) and, more importantly, due to lower incoming longwave radiation. Longwave radiation from the surrounding debris does not entirely compensate the negative longwave balance, but offsets it strongly, and peaks at midday when the temperature of the debris is highest. In the post-monsoon season the net longwave balance resulting from the three terms is negative. Mean modelled melt rates at the stakes (Table 5) correspond to the observations (Table 3). Sakai and others (1998) reported higher modelled values for the same glacier in the premonsoon and monsoon seasons.

4.3. Sensitivity analysis

Figure 12 shows the confidence interval for the parameters of the 100 best runs. The model is not sensitive to debris emissivity or surface roughness. It is sensitive to ice emissivity only at all stakes of cliff 1 in October. Here the outgoing longwave radiation plays a significant role relative to other fluxes, also since incoming shortwave radiation is limited. Similarly, the model is sensitive to debris albedo in October, since shortwave radiation reflected from the terrain plays a relatively more important role with hardly any direct radiation reaching the ice surface, due to lower solar angles. The model is most sensitive to ice albedo in both seasons. This is expected, since it determines the amount of shortwave radiation absorbed by the cliff, and the shortwave radiation flux is a major source of energy even when direct shortwave radiation does not reach the cliff (as is the case in October).

Fig. 12. Sensitivity of the model to the parameters used in the model optimization, shown by plotting all parameter values for the 100 best model runs. The box shows the 0.25–0.75 confidence interval, the cross shows the mean value. The boxplots are for all stakes in May (1.1–2.4) and October (1.1–2.3) from left to right. In red are stakes where the Monte Carlo analysis did not lead to an optimal solution.

The results of the OAT sensitivity analysis are shown for the three seasons separately (Table 6). The model is highly sensitive to slope and aspect, as well as the debris-view factor, V _d. A change in aspect affects the exposure to the sun greatly and aspect is therefore always a dominant parameter. It is slightly less important during the monsoon, when the solar angle is larger, and hence melt is less affected by a change in the direction of view. The same is true for slope, though with a lower overall sensitivity. The standard deviation, however, shows that the sensitivity to aspect is variable between the stakes, while less variable for slope. As discussed above, longwave radiation plays an important role in the overall energy balance. This explains the high sensitivity values for the debris-view factor, as well as the incoming longwave radiation itself. Melt is relatively insensitive to surface temperature, at least for the range considered of 10%, which corresponds to the observed range of measurements between the different temperature sensors located on the glacier. The changes in air temperature, important for the turbulent fluxes, have little impact on the overall melt.

Table 6. Sensitivity of the model to slope, aspect, albedo of ice, α _i, the debris- and sky-view factors, V _d and V _s, the air temperature, T _a, incoming longwave radiation, L _in, and surface temperature of the debris, T _s. The values are the respective mean slopes determined at all stakes (cm d⁻¹ 10%⁻¹). In parentheses is the standard deviation between these stakes. A high standard deviation points to a variable that is heterogeneous in space. In bold are the four variables to which the model is most sensitive in the indicated season

5. Discussion

5.1. Importance of accurate, high-resolution topography

Melt values vary greatly between the two seasons and this variability is driven by radiative fluxes, which also explain the strong variability between stakes. Since longwave and shortwave radiation fluxes are strongly affected by the topography of the cliff and of the surrounding glacier surface, we tested the effect of a coarser DEM on modelled radiative fluxes. For this, we replaced the accurate, high-resolution UAV DEM with a coarser resolution (30 m) DEM obtained from resampling the original UAV DEM and recalculated the energy balance at each stake.

The incoming shortwave radiation, incoming longwave radiation from the atmosphere and the longwave radiation emitted by the surrounding debris and reaching the cliff are shown in Figure 13 for three stake locations in May. The coarse DEM is not able to reproduce the shape of the cliff accurately. Locations that are on the cliff slope (hence shaded from the sun by the cliff) may be exposed earlier in a coarse representation of the topography. Equally, it cannot account for the local topography around the cliff, which defines the amount of longwave radiation emitted by the debris.

Fig. 13. Incoming radiative fluxes at ablation stakes (top row: stake 2.1; middle row: stake 2.3; bottom row: stake 1.2) in the pre-monsoon (grey) and post-monsoon (black) seasons, modelled with a high-resolution DEM from the UAV (solid curves) and with a low-resolution resampled DEM (dashed curves). The locations correspond to stake numbering in May. I _s (left column) is incoming shortwave radiation, L _in (middle column) is direct incoming longwave radiation from the atmosphere and L _d (right column) is longwave radiation emitted from the surrounding debris.

DEM resolution does not have a strong effect on incoming shortwave radiation, except partly at stake 1.2 on cliff 1 and stake 2.3 on cliff 2, where use of the low-resolution DEM results in an overestimation of incoming shortwave radiation in the early morning hours, when shading is not correctly reproduced. Here the cliff actually still shades the stake location, which only the high-resolution DEM accounts for. However, DEM selection has an important impact on the longwave flux from the atmosphere and the longwave flux from the surrounding debris, evident at all stakes in Figure 13. The longwave flux from the atmosphere increases when the coarser-resolution DEM is used, while the flux from the debris decreases significantly, as the debris complex topography surrounding the cliffs is not reproduced by the low-resolution DEM. For both DEMs all three fluxes are lower in October, because of the lower angle of the sun and lower air and surface temperatures.

During the pre-monsoon season, using the low-resolution DEM results in a mean overestimation of daily melt of 5.6%, with a standard deviation between stakes of σ = 4.3, which increases to 6.2% (σ = 3.1) in the monsoon and 8.8% (σ = 10.6) post-monsoon, when the sun is at its lowest angle. The relatively high standard deviation between stakes before and after the monsoon indicates that this effect is dependent on the location, namely cliff aspect, height, slope and the surrounding topography.

The self-shading of the cliff itself and an accurate representation of the surrounding topography has been neglected in earlier studies, due to the lack of a high-resolution DEM. A comparison of results with the two DEMs suggests that self-shading of the cliff is important, especially for north-facing cliffs and for locations lower down on the cliff (Fig. 13). This supports the hypothesis proposed by Sakai and others (2002) that north-facing cliffs persist while south-facing cliffs disappear quickly.

5.2. Overestimation of melt in October and refreezing

Although the model is able to reproduce the lower melt values post-monsoon, it still overestimates them to a considerable degree. An exception are two stakes located at the top of cliff 1, which are subject to considerably more shortwave radiation than the others (stakes 1.3 and 1.4). In the field, refreezing was observed on the ice surface in mid- afternoon, when the sun no longer reached the cliff surface. This may be surprising at first. However, because of the considerable dust cover at the ice surface the runoff is not clear meltwater, but rather a slurry mix of debris and water (Fig. 14b). Water that was previously observed in its liquid state slowed down considerably and became solid after a while. The two prerequisites for a refreezing process are the availability of meltwater and a negative energy balance, resulting in negative temperatures necessary to refreeze. To estimate the former is difficult both from field observations and from modelling, since measuring runoff at the cliff slope accurately is virtually impossible. The available energy for refreezing can be estimated simply using Eqn (3) with the assumption that the energy necessary for refreezing an amount of water is the opposite to that required to melt it. A look at the total energy available for melt (Q _m, Fig. 14c) shows that the energy balance becomes negative between 15:00 and 16:00 in the post-monsoon season. We therefore calculated refreezing from the negative energy between 15:00 and 19:00 by solving Eqn (3). Accounting for refreezing in this manner led to a good agreement of the model with the stake readings at all stakes. The resulting refreezing rates at all stakes is 0.05–0.15 cm h⁻¹, with a mean of 0.09 cm h⁻¹ for the considered hours.

Fig. 14. (a) Cumulative melt at stake 1.1 obtained with the model with optimized parameters and with inclusion of a simple calculation of refreezing as explained in Section 5.2. Including refreezing during the post-monsoon campaign corrected the earlier model offset at all stakes (the example shown here is stake 1.1). (b) Refreezing at the cliff was visible in the field when the water debris mix stopped flowing in the late afternoon. (c) The total energy for melt becomes negative in October as early as 16:00–17:00 and remains so until 09:00 the next day.

It is noteworthy that according to the model results inclusion of refreezing was not necessary at stakes 1.3 and 1.4, where the model worked with the initial set-up. Both are at the top of the cliff and are exposed to longer sunshine hours. Receiving higher radiative fluxes, melt from these locations is likely to run off before it refreezes and no additional meltwater is supplied from above. In this case, the first prerequisite for the refreezing process, availability of meltwater, is not likely.

While lack of a refreezing component in the model is likely to be responsible for the overestimation of melt by the model in October, radiative cooling at night might also partly account for the overestimation obtained, as part of the energy currently used for melt is needed to raise the surface temperature to the melting point before melting can take place. Finally, it should be noted that the October simulations from 23 October onwards are obtained with meteorological input extrapolated from the Kyanjing AWS, a fact that might affect our simulations. Despite this, we believe that refreezing is an important process that should be taken into account in future energy-balance model development.

When refreezing is neglected during the cold period, it would lead to considerable overestimation of total cliff backwasting, and efforts should be invested to take it into account in models. The importance of refreezing for debris-covered glaciers was also noted by Lejeune and others (2013), who suggested that discrepancies between observed and modelled values in their energy-balance model over a debris-covered glacier may be due to a neglected refreezing process, although they did not quantify this amount. Literature dealing with refreezing exists, but is largely concerned with refreezing over clean ice in the High Arctic and refreezing in the firn layer (Oerlemans, 1991; Huybrechts and De Wolde, 1999; Wright and others, 2007; Reijmer and Hock, 2008; Pelt and others, 2012; Reijmer and others, 2012). All publications lack data to test these models, a deficiency that may be challenged in the future with more measurements for ice cliffs.

5.3. Modelling backwasting during the monsoon

During the monsoon season no field measurements are available. However, since previous work has suggested that a large portion of total melt occurs during this season, we ran the model during the monsoon to provide a first estimate of energy fluxes and melt rates. Stakes installed in May had melted out at the end of the pre-monsoon season and the new stakes were installed only in October. We thus assumed the same stake locations and topographic characteristics as in May for these simulations. This is a simplifying assumption and therefore results should be regarded with some caution, as the cliffs’ topography might have changed due to backwasting.

Figure 15 shows the energy fluxes for the same stake positions as in May. While radiative fluxes from direct shortwave radiation are smaller than in May (Fig. 11), owing largely to increased cloud cover, shortwave radiation reflected from the terrain, D _t, and longwave radiation from the surrounding debris, L _d, are relatively more important than during the dry seasons. The net total longwave flux is less negative, since incoming longwave radiation is considerably higher due to clouds. In contrast to the dry seasons, the latent heat flux is positive and more important in the monsoon. As a result of lower incoming shortwave radiation combined with higher longwave radiation and turbulent fluxes, melt rates are similar to those pre-monsoon (Table 7). The mean values at the stakes, 3.5–5.2 cm d⁻¹ (Table 7), are lower than reported by Sakai and others (1998) for the same glacier in this period. Maximum values are, however, >9 cm d⁻¹ on cloud-free days.

Table 7. Mean values of energy fluxes at each stake during the monsoon season. I _s is direct incoming and D _s diffuse shortwave radiation, D _t shortwave radiation reflected from the surrounding terrain and I _ref reflected from the ice. L _in is incoming longwave radiation from the atmosphere, L _d from the surrounding terrain and L _o outgoing longwave radiation. LE and H are latent and sensible heat fluxes

Fig. 15. Mean hourly energy fluxes during the monsoon season, modelled at the stake locations from May. L _in − L _o is the difference between incoming atmospheric longwave radiation and outgoing longwave radiation from the ice. L _d is longwave radiation emitted from the surrounding terrain. I _s is direct shortwave radiation from the sky, D _s and D _t diffuse and terrestrial shortwave radiation. H and LE are turbulent fluxes (sensible and latent heat).

6. Conclusions

In this paper, we have improved an existing model to assess ice-cliff backwasting on a debris-covered glacier using a dataset of field measurements collected during two seasons with different climatic properties. For both seasons a high-resolution DEM was available for the debris-covered surface. Use of the DEM in combination with improved algorithms for calculation of shortwave and longwave radiative fluxes allowed a very accurate representation of the complex topographic effects typical of debris-covered glaciers. The model was optimized with a Monte Carlo-type approach. We calculated rates of backwasting for the two seasons and calibrated the model by comparison with the observations at ablation stakes drilled in the cliffs. We used calculations of the energy fluxes to explain the observed variability in melt rates. We also looked at the model sensitivity to both model parameters and input variables. Our main conclusions are summarized below.

Observed melt rates at numerous stakes over two cliffs exhibited large variability, which is difficult to explain using single variables (e.g. slope or aspect). It is likely to be a combination of several factors and their interaction. Mean melt rates were in the range 3.25–8.6 cm d⁻¹ in May and 0.18–1.34 cm d⁻¹ in October. With lower absolute radiative fluxes and temperature, melt rates are considerably lower in the post-monsoon season than pre-monsoon.

The model includes a new, accurate representation of the sky- and debris-view factors and surrounding topography. In this form, it was able to reproduce stake readings well in May and October (for the initial period of multiple readings). In October, when multiple stake readings a day were available the model can reproduce the daily cycle. However, it fails to reproduce the very low melt rates at lower stakes in October. This is probably because the model does not include refreezing at the ice surface. A simple model adaptation that includes refreezing can reproduce the stake readings and results in refreezing rates of on average 0.09 cm h⁻¹ during the afternoon hours when refreezing was observed.

Modelled melt rates during the monsoon are slightly lower than pre-monsoon. While solar radiation is markedly lower due to the increased cloud cover, longwave fluxes from the surrounding terrain as well as from the atmosphere increase.

The influence of the local topography (shading the ice surface from incoming solar radiation and exposing it to more longwave radiation and reflected shortwave radiation from the terrain) at a small scale was shown to be important. It becomes more important during times of the year when the solar angle is low, when, while the glacier surface may be reached by the sun’s rays, north-facing steep cliff faces are in the shade longer in the morning and earlier in the afternoon. This is expressed by a high sensitivity of the model to the topographic variables of slope and especially aspect, as well as the sky- and debris-view factors.

A high-resolution DEM of the glacier surface enabled us to differentiate between the radiative fluxes from the sky and the surrounding terrain in greater detail, allowing us to explain the highly variable melt rates between different cliff locations. While shortwave radiation is the dominant forcing for melt, especially in the pre-monsoon season, longwave radiation from the surrounding terrain plays an important role, contributing between 43 W m⁻² (post-monsoon) and 87 W m⁻² (during the monsoon) to the energy balance, and causes the net longwave radiation balance in the pre-monsoon and monsoon periods to become slightly positive during the day. Areas on top of a cliff are more exposed to direct radiation from the sky, while locations further down, or those covered by surrounding terrain, may not receive any direct influx at all. This deficit, however, is to some degree offset by radiation from the terrain. We showed that disregarding this local topography and using a low-resolution DEM resulted in an overestimate of mean daily melt at the point scale of 5.6–8.8%, depending on the season.

The sensitivity of ice-cliff backwasting to the heterogenic topographic properties of the ice cliff itself and the surrounding glacier surface suggests that only an energy-balance model that can take these factors into account can accurately reproduce melt on such surface forms. Neglecting the local topography could lead to an overestimate of melt from a debris-covered glacier on a distributed scale. The model applied in this paper is still a point-scale model, in the sense that it is applied at each stake separately, despite the fact that all calculations of radiative fluxes are grid-based. A next step would be to implement a grid-based model that can calculate ablation from the entire cliff. A major challenge still, however, is the extrapolation of the meteorological variables other than shortwave radiation at the grid scale.