Scholarly article on topic 'Refined energy-balance modelling of a supraglacial pond, Langtang Khola, Nepal'

Refined energy-balance modelling of a supraglacial pond, Langtang Khola, Nepal Academic research paper on "Earth and related environmental sciences"

Share paper
Academic journal
Annals of Glaciology

Academic research paper on topic "Refined energy-balance modelling of a supraglacial pond, Langtang Khola, Nepal"

Annals of Glaciology 57(71) 2016 doi: 10.3189/2016AoG71A421 29

© The Author(s) 2016. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons. org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.

Refined energy-balance modelling of a supraglacial pond,

Langtang Khola, Nepal

Evan S. MILES,1 Francesca PELLICCIOTTI,2'3 Ian C. WILLIS,1 Jakob F. STEINER,2

Pascal BURI,2 Neil S. ARNOLD1

1 Scott Polar Research Institute, University of Cambridge, Cambridge, UK 2Institute for Environmental Engineering, Swiss Federal Institute of Technology (ETH), Zürich, Switzerland 3 Department of Geography, Northumbria University, Newcastle upon Tyne, UK Correspondence: Evan S. Miles <>

ABSTRACT. Supraglacial ponds on debris-covered glaciers present a mechanism of atmosphere/glacier energy transfer that is poorly studied, and only conceptually included in mass-balance studies of debris-covered glaciers. This research advances previous efforts to develop a model of mass and energy balance for supraglacial ponds by applying a free-convection approach to account for energy exchanges at the subaqueous bare-ice surfaces. We develop the model using field data from a pond on Lirung Glacier, Nepal, that was monitored during the 2013 and 2014 monsoon periods. Sensitivity testing is performed for several key parameters, and alternative melt algorithms are compared with the model. The pond acts as a significant recipient of energy for the glacier system, and actively participates in the glacier's hydrologic system during the monsoon. Melt rates are 2-4 cm d-1 (total of 98.5 m3 over the study period) for bare ice in contact with the pond, and <1 mm d-1 (total of 10.6 m ) for the saturated debris zone. The majority of absorbed atmospheric energy leaves the pond system through englacial conduits, delivering sufficient energy to melt 2612 m3 additional ice over the study period (38.4 m3 d-1). Such melting might be expected to lead to subsidence of the glacier surface. Supraglacial ponds efficiently convey atmospheric energy to the glacier's interior and rapidly promote the downwasting process.

KEYWORDS: debris-covered glaciers, glacier ablation phenomena, glacier hydrology, glacier mass balance


Debris-covered glaciers comprise ~10% of the glacierized area in High Mountain Asia, where many glaciers have been rapidly losing mass in recent years (Benn and others, 2012; Bolch and others, 2012; Pellicciotti and others, 2015). Although the insulating effect of thick debris is known to reduce ablation (0strem, 1959; Scherler and others, 2011; Ragettli and others, 2015), the impact of surface ponds and their associated ice cliffs on the ablation process is much less understood (Sakai and others, 2000; Benn and others, 2001, 2012; Röhl, 2006). Several studies have shown that surface ponds and terminal lakes of debris-covered glaciers are associated with rapid thinning (Basnett and others, 2013; Pellicciotti and others, 2015) and retreat (Benn and others, 2000; Gardelle and others, 2011; Sakai, 2012). A conceptual model of debris-covered glacier response to climate change has been advanced, linking supraglacial pond development to subsequent terminal lake formation (Benn and others, 2000; Reynolds, 2000; Quincey and others, 2007; Benn and others, 2012).

Studies have identified glacier slope, velocity and thinning as controls on supraglacial pond formation and persistence (Benn and others, 2001; Gulley and Benn, 2007; Sakai and Fujita, 2010; Salerno and others, 2012) as these influence the surface's intersection with englacial conduits and crevasses. Observations have highlighted key melt mechanisms associated with ponds, including calving and waterline melting (Reynolds, 2000; Benn and others, 2001; Röhl, 2008), and the rapid backwasting of bare-ice cliffs that are often adjacent

to them (Benn and others, 2001; Han and others, 2010; Reid and Brock, 2014; Steiner and others, 2015).

However, few studies have attempted to quantify the energy exchanges associated with supraglacial ponds, or their effects on glacier ablation. Xin and others (2011) identified two kinetic types of melt: (1) winds may force currents to drive thermo-erosion and notch development near the lake surface and (2) free convection, due to the density/temperature relationship of water, may drive pond circulation and therefore promote melt along the entire water/ice interface. Sakai and others (2000) and Röhl (2008) each adapted empirical relationships from iceberg melt observations to examine subaqueous and waterline melting of ice cliffs, while Lüthje and Pedersen (2006) adapted a method based on free convection to study basal melting of ponds on the Greenland ice sheet. To date, no effort has been made to compare these algorithms or to apply a melt model based on physical principles to supraglacial ponds on debris-covered glaciers. This is the aim of the present study.


Here we present the framework and calculations of the model, while the source data and parameters for analysis are described in the next section.

Physical description of the system

Supraglacial lakes on debris-covered glaciers are complex systems with multiple boundary exchanges of energy and

Fig. 1. Conceptual (a) mass and (b) energy exchanges for the pond control volume.

mass. This paper closely follows the analytical framework of Sakai and others (2000), although improvements have been made to nearly all computations. The approach treats the whole lake as a single control volume which can change in mass and volume (Fig. 1a), as reflected, in practice, by a change in water level, due to adjacent ice-cliff melt or catchment runoff, V\, discharge, Vd, rainfall, \/R, or latent fluxes, \/LE, all with units m3s-1. Three distributed fluxes (W) need to be modelled for energy-balance calculations: (1) at the air/atmosphere surface of the pond's water, Qn; (2) at the saturated debris surface at the pond's base, Qd; and (3) at the subaqueous bare-ice surface, Qi, and all need to be scaled by their respective areas (Fig. 1b). Additionally, energy can be advected into (I) or out of (D) the system via mass transfer, leading to variations in stored energy, AS.

Whole-pond control volume energy balance

Considering the control volume as a reservoir of mass and energy, the internal energy of the pond, 5 (J), that is available for melt can be calculated using

5 = cwpwVp(Tp- 273.15 K) (1)

where cw (J kg-1 K-1) is the specific heat of water, pw (kg m_3) is the density of the water, Tp (K) is the mean pond temperature and \ p (m3) is the current volume of the pond. Changes in this stored energy must be compensated by energy transfer at the control volume boundaries, expressed as

AS = Qn4p + QiA + QdAd + I + D (2)

Here Qn is the net surface flux from the atmosphere, Qi is the energy exchange at subaqueous bare-ice surfaces, Qd is

the energy exchange through subaqueous debris, I (W) is energy advected into the pond by runoff inputs and D (W) is the energy removed from the pond by discharge. Ap, Ai and Ad (m2) are the areas of the pond, ice and debris boundaries with the atmosphere.

Pond energy exchanges

At the water surface, energy exchanges between the pond and atmosphere are given by

Qn = In + Ln + H + LE + Qr

with In and Ln net shortwave and longwave radiative fluxes, H and LE sensible and latent turbulent fluxes, and Qr advected energy due to rainfall (all Wm~2). Most terms are calculated as described by Steiner and others (2015), so here we provide only a brief description of the terms, highlighting the differences between our equations and those of the earlier study.

Shortwave radiation

Incoming shortwave radiation is determined based on individually modelled direct, Is, diffuse, Ds, and debris-reflected, Dt, radiation components (section 3.1.1 of Steiner and others, 2015), incorporating view factors determined using the lake's position and a 0.2 m resolution digital elevation model (Immerzeel and others, 2014). The pond's albedo is estimated using the empirical relation previously applied by Sakai and others (2000),

= 0.78 es

which is determined from Tsho Rolpa Lake, where a is albedo and ds is solar elevation. The net shortwave flux is then calculated as

In = (Is + Ds + Dt)(1 - ,

Longwave radiation

The net longwave radiation is calculated as a combination of atmospheric, Lin, and debris, Ld, sources, each scaled by view factors (section 3.1.2 of Steiner and others, 2015) and radiation emitted by the pond surface, Lo:

Ln = Lin + Ld — Lo Lo is calculated using the Stefan-Boltzmann law

Lo = Tws4 (7)

with the pond's surface temperature (Tws (K)) and emissivity ew = 0.95 (Sakai and others, 2000).

Turbulent fluxes and rain

Turbulent fluxes are determined using the bulk aerodynamic method with atmospheric stability correction, as implemented by Han and others (2010), Reid and Brock (2014) and Steiner and others (2015). The only changes made are associated with the application to a water surface rather than debris, resulting in a saturated film of air at the water temperature and the use of an appropriate surface roughness value for a water surface (z0 = 5 x 10_4 m from literature values ranging between 2.7 x 10_s and 1 x 10_3 m). The specific energy flux associated with rainfall is estimated using the rainfall rate, qr (ms-1), and the air temperature, Ta (°C):

Qr — cwPwTaqr

Subaqueous subdebris melt

Energy transfer through the saturated debris zone at the pond bottom is assumed to occur only via conduction (Sakai and others, 2000; Röhl, 2008). Like Röhl (2008), we calculate the bulk thermal conductivity, kd (J kg-1 K-1), as a combination of the values for rock, kr, and water, kw, scaled by porosity,

kd = kr(1 - kw^ (9)

The energy flux through this layer, Qd (Wm-2), can then be calculated based on the pond bottom water temperature, Twb (°C), the ice temperature, Ti (assumed to be at the freezing point 0°C), and the thickness of the saturated debris zone, dd (m):

kd(Twb - Ti)

The melt rate, v^ (ms-1), is the energy flux divided by the latent heat of fusion for water, Lf (J kg-1):

vd = 17

Subaqueous bare-ice melt

Several different methods have been applied to model subaqueous melt rates, vi (ms-1). Sakai and others (2000) applied an empirical relation determined from iceberg melt observations (Weeks and Campbell, 1973). This method assumed a strong forced-convection scenario, based on average pond temperatures, Tw (°C), and a contact length, xi (m), but the pond's vertical velocity at the bare-ice contact was fixed at uw = 0.02 m s-1:

vi = 7.14 x 10-

(Tw - Ti)

Röhl (2008) implemented a different empirical relation, which depended entirely on water temperature (here converted to ms-1):

1.8 x 10-2(Tw + 1.8)15

vi =----— (13)

i 24 x 60 x 60

However, the equation was initially formulated based on the freezing point of sea water at —1.8°C (Russell-Head, 1980), so a more appropriate equation would implement the 0°C freezing point of fresh water instead. Lüthje and Pedersen (2006) applied free-convection algorithms optimized for horizontal plates (Linden, 2002; Taylor and Feltham, 2004) to estimate subaqueous melt rates on the Greenland ice sheet:

rcw (ßgKl2 Y/3(T _T )4/3

Lf V Ц ) ('w Ti)

where r = 0.1 is a dimensionless constant, g = 9.81 ms-2 is gravitational acceleration and is the temperature-dependent thermal conductivity of water (Wm-1 K-1).

Our approach takes a step back to the driving causes of melt at the ice/water interface. While wind-driven water currents leading to forced convection can develop for ponds with sufficient fetch (Sakai, 2012), many ponds are much smaller and vertical water velocities are negligible (Xin and others, 2011). Instead, natural convection can occur, where temperature-dependent density differences drive a convec-tive current. Freshwater density peaks at 4°C, so for temperatures below this point, colder parcels will tend to rise, driving circulation. For this analysis, suspended sediment

concentrations (SSC) are assumed to be near-constant (i.e. the pond is well mixed), so only thermal differences drive circulation. Observed SSC values for similar ponds are 10400 mg L-1 (Bhatt and others, 2007; Takeuchi and others, 2012), and SSC variations of >100 mg L-1 have a sufficiently strong effect on density to drive stratification for much larger ponds (Chikita and Joshi, 2000). The case of natural convection along a vertical plate was examined and linearized by Churchill and Chu (1975), whose approach is applied here to the case of water near the freezing point.

This first requires evaluation of the Grashof number, Gr, representing the ratio between buoyancy and viscosity:

gDl3(Pws - P0)

where Di (m) is the mean vertical subaqueous ice/water contact length (characteristic length for this geometry), pws and p0 are the water density (kg m-3) at the lake surface and freezing point, respectively, which together drive convection, and pf and u are the density and kinematic viscosity evaluated at Tf, the mean of the surface and freezing-point temperatures (°C). The Prandtl number, Prf (the ratio between kinematic viscosity and thermal diffusivity), is evaluated at Tf to determine the Rayleigh number, Ra:

Ra = GrPrf

which is used to determine whether the flow is turbulent and to calculate the Nusselt number, Nu:

С 0.68 + r °.67 Ra1/4

1 _i_ I 0.492

l + l Prf

) 9/16

0.825 + r °.387 Ra1/6

1 _L Í 0.492

l + \ Prf

105 < Ra < 109

Ra > 109

Finally, the mean coefficient of heat transfer, hi, can be evaluated and applied to determine the heat flux, Qi (Wm-2), and melt rate, vi (ms-1), due to free convection:

Nu k„

Qi = hi (Tp- Ti) Qi

Mass balance of the pond

The mass balance of the pond according to Figure 1a is given by:

AV = V + VLE + VR + MWi + MWd- Vd (21)

where the observed changes in volume, AV, are accounted for by inflows, Vi, generated melt, MWi, MWd, exchanges of vapour, VLE, or rain, \/R, at the pond's surface, or outflows, Vd. The volumetric rate of evaporation or condensation can be calculated directly from the magnitude of the latent energy flux, QLE:

VLE = (22)

cwpws ' ws

The rain input is simply the precipitation rate scaled by lake

Fig. 2. Study site on Lirung Glacier and station locations used in this study, with estimated pond depths indicated. The three bodies of water are hydraulically connected through the thick debris cover, and are modelled jointly.

area, Al,

Vr = qA and the outflows are given by

\d = D

cwpws Tp

Two fluxes remain unspecified at the control volume boundaries (Fig. 1): mass and energy brought into the system by inflows, \Vi and I, and mass and energy removed from the system by discharge, Vd and D. Our field measurements indicate surface water temperatures of <0.1°C for all surveyed flows emerging at the glacier surface or at the terminus, suggesting that a negligible amount of energy is added to the system by inflows (I = 0J). This enables the simultaneous solution of the mass and energy equations, as the outflow discharge energy, D, may be calculated directly from Eqn (2), and assuming that this discharge removes water at the average pond temperature, Tp, the volumetric discharge, Vd, can be estimated from Eqn (24). Finally, the flow of water into the pond, \Vi, can be inferred from Eqn (21).


The above framework has substantial data requirements and has many physically meaningful parameters that are difficult to constrain, and cannot be calibrated. As a result, we adopt a diagnostic modelling approach to identify the importance of different lake processes. This section describes the data collected and physical parameter selections for driving the standard model configuration.

Study site and field instrumentation

To develop and test the model, a supraglacial pond was instrumented on Lirung Glacier in the Langtang Valley of Nepal (85°33'43" N, 28°13'57" E; Fig. 2). Lirung Glacier is situated in the Upper Langtang catchment, which is ~30% glacierized (Ragettli and others, 2015), and where the glaciers exhibit debris cover on tongues totalling 25% of the total glacierized area (Shiraiwa and Yamada, 1991). Lirung Glacier covers a very large elevation range (40407180 m a.s.l.), and has very heterogeneous debris cover over ~20% of its 6 km2 total area (Shiraiwa and Yamada, 1991; Immerzeel and others, 2014).

Climate in the Langtang Valley is characterized by a strong monsoon, with ~70% of the annual precipitation falling between June and September (Immerzeel and others, 2011). This period also contains the highest daily maximum and minimum temperatures experienced at the site, and consequently the highest rates of ablation occur during these months. Very cold and dry conditions are dominant during the winter, while the pre- and post-monsoon transition periods experience intermediate temperatures and occasional atmospheric instabilities (Ragettli and others, 2015).

The study pond at 4070 m a.s.l. was monitored beginning on 8 May 2013, equipped with a HOBO water-level logger to measure water pressure and pond bottom water temperature (just above the saturated-debris zone). The sensor was cast into the pond, approximately equidistant from the shore and ice cliff (Fig. 2). The pond's drainage led to the subaerial exposure of the instrument on 15 July 2013 (Fig. 3), so the study period is day-of-year (DOY) 130-197, encompassing the late pre-monsoon and the early monsoon

Fig. 3. The study pond in (a) May and (b) October 2013, showing ice-cliff retreat and a decreased water level after the monsoon.

of 2013. For May-October 2014 a float was installed with a surface temperature sensor and tethered to the pressure transducer cable to investigate temperature variability within the pond. An automatic weather station (AWS Lirung; 4076ma.s.l.) was installed 80 m to the south of the pond to monitor atmospheric conditions at the site, recording incoming shortwave radiation, Is, 2 m air temperature, Ta (shielded and vented), relative humidity, rH, and wind speed, ua, and direction. All variables were measured at 5 min intervals and aggregated to 1 hour values. The barometric sensor from an off-glacier AWS 2.3 km from the study pond at Kyanjing village (AWS Kyanjing; Fig. 2) was adjusted to the pond altitude using the ideal gas law, then applied to correct the HOBO pressure signal for conversion to at-sensor depth according to the hydrostatic equation. Precipitation data were taken from the tipping bucket gauge at AWS Kyanjing. Lirung Glacier was observed in May and October 2013 using an unmanned aerial vehicle (UAV), from which a detailed digital elevation model (DEM) and orthophoto were developed (Immerzeel and others, 2014) covering the study pond.

Pond geometry

Converting the pressure records to lake volume requires a depth/area/volume relationship. Basic observations with a sonar transducer were limited, due to the hazard of rockfall from the small pond's ice cliff, but measurements of pond depths increased linearly approaching the ice cliffs. This suggests a vertical or undercut subaqueous ice cliff, in agreement with observations of notch formation and overhangs at the base of ice cliffs on other debris-covered glaciers (Benn and others, 2001; Röhl, 2008). Considering that a pond in a supraglacial depression likely once filled the entire depression, the current exposed terrain may be similar to the pond basal topography. The slope of Lirung Glacier's surface depressions was investigated by identifying all closed depressions >50m2 in the UAV DEM and collecting the slope statistics for each depression and the glacier as a whole. The glacier-wide mean slope for such depressions was 53.2% (28°), while the mean value for the study pond's depression was 67.8% (34°).

The depth beneath water surface at the time of DEM acquisition was then extrapolated from the pond shore (excluding the ice-cliff edge) using the 53.2% slope to produce an estimate of the pond depth (Fig. 2), which was added to the pond elevation to produce a grid of subaqueous surface altitudes. Finally, to develop a depth/area/volume

curve, the minimum altitude of the pond bottom was determined and incremented by 0.01 m, determining at each step the submerged debris area, pond volume, ice-cliff contact length and mean ice-cliff depth. A time series of water level altitudes was determined using the hydrostatic equation and the pond's surface elevation at the time of UAV acquisition, and time series of geometric properties were determined from the depth/area/volume curve.

Pond temperature estimation

To apply the model, several measures of water temperature are needed. The pond bottom temperature, Tb, is measured directly and assumed to be uniform, and we prescribe an ice temperature, Ti, at the ice/water interface of 0°C for the subaqueous ice cliff and the saturated debris zone. The mean water temperature for the pond, Tw, the water surface temperature, Tws, and film temperature, Tf, are all necessary for the energy-balance calculations.

Field observations in this study and those of other investigators (Sakai and others, 2000; Röhl, 2008) indicate very small water surface horizontal velocities for small ponds, suggesting that pond overturning is minimal, other than at the ice interface and near inlet and discharge mass exchanges. This limited evidence of widespread convection suggests that the pond's average temperature and surface temperature could be estimated from the temperature gradient within the pond. Our sensors showed strong temperature fluctuations at the pond's surface and the frequent occurrence of an inversion layer (Fig. 4), similar to that observed by Xin and others (2011) in timing and gradient magnitude, although the pond was much shallower during the 2014 observation period than the 2013 model period here. Mean observed hourly gradients between the water surface and pond bottom were applied with the 2013 monsoon observed pond bottom temperatures and pond depths, to generate surface and average temperatures for the study period (Fig. 4). The film temperature, Tf, is derived according to Churchill and Chu (1975) as the mean of the surface and freezing-point temperatures.

Sensitivity analysis

Due to a paucity of calibration and validation data for this study, several sensitivity tests are performed to consider the effects of key assumptions and parameter choices. For each test, the resulting subaqueous melt and estimated discharge were determined. The test-run titles are in italic, and model outcomes are summarized in Table 1.

Fig. 4. Observed hourly water column temperature gradients (positive downward), coloured by depth at the time of measurement for May-October 2014. The mean hourly gradient is shown as the solid black curve, while the ±1 standard deviations are denoted by +. Depth is according to the pressure transducer at the time of the joint surface and pond bottom temperature measurements. Night-time temperatures exhibit a strong inversion, which switches due to diurnal heating at the surface.

Configuration tests

This study attempts to represent pond geometry in a more realistic manner than previous analyses, but is limited by the assumption of a constant subaqueous slope of 53.2%. This is varied between 40% (21.8°) and 80% (38.7°) to consider maximum subaqueous melt quantities based on reasonable subaqueous geometries.

The observed temperature gradients in the pond exhibit hourly clustering over the 2014 monsoon season, but these measurements were made across a narrow range of pond depths. The pond temperature is a critical variable for most

components of the model, so four distinct temperature distributions were used to consider a range of possible outcomes: the mean temperature gradient based on the 2014 observations, the mean gradient plus and minus one standard deviation, and a uniform pond temperature as measured at the pressure transducer (Tw = T|b = Tws).

In addition to the free-convection approach presented in this paper (Eqns (15-20)), three other subaqueous bare-ice melt algorithms were evaluated (Sakai and others, 2000; Lüthje and Pedersen, 2006; Röhl, 2008), as in Eqns (12-14). The Sakai algorithm (Eqn (12)) was tested with three assumed vertical pond velocities (runs sakai, sakai0.01 and saka/0.04), and the Röhl algorithm (Eqn (13)) was tested with freezing points based on sea water (rohl) orfresh water (rohl0C), while the Lüthje algorithm (luthje; Eqn (14)) is independent of parameter assumptions, though all three algorithms would also be strongly affected by assumed temperature gradients.

The discharge is assumed to be at the mean pond temperature, but drainage points could occur at the surface or base of the pond. Using these temperature estimates will bound the estimated discharge flow required to balance the energy equation (Eqn (2)). Inputs to the pond are modelled to be at the freezing point, and therefore do not add melt-available energy to the pond system. Field measurements indicate supraglacial runoff commonly <0.1°, but the potential effect of this assumption is tested by assuming inlet temperatures of 0.1°C and 0.25°C, encompassing the range of values observed by Takeuchi and others (2012).

Parameters varied

Key parameters affecting the energy exchanges were varied within a range of literature values. Table 1 indicates the parameter values used in the Standard model run and the alternative values used for the sensitivity tests. The

Table 1. Principal results of model runs. - indicates a result no different to the Standard run. MWi refers to the cumulative subaqueous bare-ice melt volume, MWj refers to the cumulative subaqueous subdebris melt volume, Vd is the mean pond discharge, vi and vJ are the subaqueous ice and debris mean melt rates and Vdmax is the peak pond discharge

Run title Parameter Run value Standard value MWi vi MWd vJ Vd Vdmax

m3 m d-1 m3 m d-1 m3 s-1 m3 s-1

Standard - - - 98.5 0.029 10.6 3.4 ; <10-4 0.028 0.24

slope040 basal slope 40% 53.2% 78.8 - 9.9 - 0.026 0.24

slope0678 basal slope 67.8% 53.2% 117.1 - 11.2 - 0.030 0.25

slope080 basal slope 80% 53.2% 131.9 - 11.6 - 0.031 0.25

Tminuslsig VT obs ß — G obs ß 97.6 0.029 - - - 0.32

Tpluslsig VT obs ß + G obs ß 107.6 0.032 - - 0.030 0.27

Tlb VT 0; T = Tlb obs ß 94.1 0.028 - - 0.023 0.16

sakai Sakai algorithm - - 79.3 0.024 - - - 0.24

sakai0.01 Sakai, uw 0.01 ms-1 0.02 ms-1 45.5 0.014 - - - 0.24

sakai0.04 Sakai, uw 0.04 ms-1 0.02 ms-1 138.1 0.042 - - 0.029 0.24

rohl Röhl algorithm - - 329.3 0.098 - - 0.023 0.23

rohlOC Röhl, Tmelt 0°C -1,8°C 92.1 0.028 - - 0.026 0.24

luthje Lüthje algorithm - - 97.9 0.030 - - 0.026 -

TdTlb Td Tlb TP - - - - 0.024 0.15

TdTws Td Tws TP - - - - 0.076 2.00

TiO.1 Ti 0.1°C 0°C - - - - 0.032 0.35

Ti0.25 Ti 0.25°C 0°C - - - - 0.039 0.99

dd0.5 dd 0.5 m 1.5 m - - 31.9 1.00 x10-3 0.029 -

dd3 dd 3m 1.5 m - - 5.3 1.67 x 10-4 - -

keffSakai keff 0.4 J K-1 m-1 s-1 1.28 J K-1 m-1 s-1 - - 3.3 1.04 x 10-4 - -

keffO.565 keff 0.565 J K-1 m-1 s-1 1.28 J K-1 m-1 s-1 - - 4.7 1.47 x 10-4 - -

keff2 keff 2 J K-1 m-1 s-1 1.28 J K-1 m-1 s-1 - - 16.6 5.22 x 10-4 - -

Fig. 5. Subsets of modelled time series for the 2013 pre-monsoon (a-e) and monsoon (f-j). The surface energy balance is the key exchange of energy, dominated by net shortwave, In, and latent, E, fluxes (a, f). Stored energy, S, diurnally fluctuates by 300% of its base value (b, g). Subaqueous subdebris melt, v^, is low, depending entirely on parameter selection and pond temperature (c, h). Timing of subaqueous bare-ice melt, v., varies greatly, based on algorithm choice and pond temperature (d, i). Discharged energy, D, balances the energy budget for the pond (e, j). Heavily overcast conditions prevailed during DOY 166-169, resulting in continuously positive net longwave radiation, Ln, but a reduction in net surface energy inputs (f), leading directly to a decrease in modelled discharge, D (j). During this period, pond temperatures declined, resulting in lower modelled v^ and v, (h, i).

conductivity of the saturated debris layer, keff, is varied between the values for water (kw = 0.565 J K-1 m-1 s-1) and rock (kr = 2.0J K-1 m-1 s-1) given by Röhl (2008), while the value for permafrost used by Sakai and others (2000) is also tested (keff = 0.4J K-1 m-1 s-1). The assumed subaqueous debris thickness, dd, is also adjusted to 0.5 and 3.0 m, the limits of debris thickness observed at ice-cliff exposures near the study site. For the surface energy balance, emissivity, ew, and surface roughness, z0 (m), are varied among common literature values for water. The albedo model used by Sakai and others (2000), normally evaluated for the date and hour and based on the pond's geographic location, is instead replaced by fixed values at the extremes of on-site observations (a = 0.08, 0.12).

RESULTS Standard run

Based on the subsurface geometry estimate and the pressure transducer data, the pond contained 1250 m3 of water and had a surface area of 650 m2 at the beginning of the study period. The water level was initially stable then drained slowly, lowering ~2.3 m in 50 days before the sensor was exposed subaerially. At this time, the pond's volume was estimated to be 200 m3, with a surface area of 400 m2. Observed pond bottom temperatures fluctuated 0.8-3.0°C, with peaks at midday, while calculated pond surface temperatures varied between 3.5°C at midday and freezing at night. Consequently, the melt-available energy stored within the pond rises dramatically from baseline values

~4 x109J at night to peak values ~10 x109J at midday (Fig. 5b and g).

The diurnal peaks in stored energy are supplied by the residual of the pond surface energy balance (Eqn (3)), which often peaks above 1000 Wm~2 (Fig. 5a and f). The peak shortwave balance is commonly >800Wm~2, with the longwave balance fluctuating between 50 and 50 W m 2. While sensible and rain energy fluxes have a minimal effect (peaks of 40 and 10 W m 2, respectively), the latent surface flux commonly peaks at 200 Wm~2. The latent flux also switches roles: in the dry pre-monsoon (Fig. 5a), it is an energy sink as the pond surface evaporates, but condensation is prevalent in the wet monsoon, when it is an energy source (Fig. 5f).

Modelled subaqueous subdebris melt rates (Eqns (9) and (10)) were very low, between 2 x 10_4 and 6 x 10~4md-1 (Fig. 5c and h), while melt rates for subaqueous bare ice (Eqns (15-19)) were considerably higher, between 0.01 and 0.06 m d-1 (Fig. 5d and i). Scaled by the areas of these surfaces, the cumulative subaqueous melt over the period of record was 8.3 m3 (0.12 m3d-1) for saturated debris and 98.5 m3 (1.4 m3 d-1) for bare ice (Standard in Table 1). Both quantities peak at midday, when both the pond bottom and water surface temperatures are highest.

The sum of all surface energy fluxes (Eqn (2)) has peak values up to 8 x 105 W. This is dominated by the surface energy balance, which is an order of magnitude higher than the rate of change in stored energy. The excess energy (Fig. 5e and j) is accounted for by the pond's discharge, which has an average value of 0.028 m3s-1. The calculated discharge

peaks at 0.1-0.3 m3s-1 in the late morning, then slowly decreases through the afternoon. The calculated influx nearly matches the discharge, as mass changes from melt, evaporation/condensation and rain are all much smaller.

Sensitivity analyses

A summary of the results of the Standard model run and the sensitivity analyses is presented in Table 1.


The application of different slope values resulted in a linear change in the pond's calculated initial and final volumes, ranging from 950 to 160 m3 (run slope040, 40% slope) and 1800 to 800 m3 (run slope080, 80% slope). While the slope is important for the entire energy budget of the pond, it is particularly important as it determines the pond's subaqueous surface areas, directly controlling the cumulative melt. It also determines the ice-face depth, which is the critical length for free convection, so it affects the melt rate at the ice face, vf. This removes more or less energy from the

pond, which has a small effect on mean discharge, Vd|. As shown in Table 1, of the four configurations, slope080 produced the highest average and peak discharges (0.031 and 0.25 m3s-1, respectively) and run slope040 produced the lowest average and peak discharges (0.026 and 0.24 m3s-1, respectively). Slope also has a small effect on the surface energy balance by altering Tws (Table 1).

Pond temperature

The different pond temperature scenarios had no effect on subaqueous subdebris melt, which is calculated based on the pond bottom temperature. They had a only a small effect on the subaqueous bare-ice melt, as this depends on the pond's maximum temperature. Of the three scenarios, run Tlb produced the lowest estimate of cumulative melt (94.1m3) and melt rate (0.028md-1), and run Tpluslsig produced the highest cumulative melt (107.6 m3) and melt rate (0.032 md-1). The temperature scenarios had a substantial effect on the calculated energy stored within the pond, resulting in high estimated discharges to accommodate the high rates of change within the energy reservoir for the high gradient case, and a very low average discharge for the no-gradient case (Table 1).

Melt models

The Sakai algorithm was evaluated for uw = 0.01, 0.02 and 0.04ms-1, giving 45.5, 79.3 and 138.1m3 of meltwater, respectively (Table 1). The Röhl algorithm estimates 329.3 m3 of meltwater, but adapting the equation to fresh water's 0°C freezing point gives 92.1 m3 of meltwater. The algorithm of Lüthje and Pedersen (2006) generates 97.9 m3 of meltwater. The different algorithms produce different average discharge values, based on the amount of energy removed from the pond by the melt process. Excluding the Röhl algorithm as applied to salt water, the Lüthje algorithm produces the lowest average (0.026 m3) and peak (0.24 m3s-1) discharge estimates. The Sakai algorithm with uw = 0.04 ms-1 generates the highest average (0.029 m3s-1) and peak (0.24 m3s-1) discharge estimates, closely followed by the standard run (Table 1).

Inflow and discharge temperature

Two alternative assumptions were considered for the

temperature of the pond's discharge, compared with

Td = Tw for the Standard run. For Td = Tlb (run TdTlb), the average pond discharge is slightly reduced at 0.024 m3s-1, but the peak discharge estimate is halved to 0.15 m3 s-1. The opposite is true for Td = Tws (run TdTws), which estimates a high average discharge of 0.076 m3s-1 and peak discharge of 2.00 m3s-1. These assumptions do not affect the melt calculations.

Two distinct scenarios were assessed for the temperature of inflows, set to 0°C for the Standard run. For Ti = 0.1°C (run Ti0.1) the slight increase in pond energy must be accommodated by a moderate increase in pond discharge to an average value of 0.032 m3 s-1 and maximum value of 0.35 m3s-1. For Ti = 0.25°C (run Ti0.25), the increase is substantial, to an average value of 0.039 m3 s-1 and maximum value of 0.99 m3s-1. Beyond Ti = 0.25°C, Ti approaches Td and the model cannot shunt stored heat with the discharge, producing unrealistic estimates.


The subaqueous subdebris melt is determined from two unconstrained variables, the debris thickness, dd, and the effective saturated debris thermal conductivity, keff. Changing these parameters resulted in subdebris melt estimates ranging from 44% to 300% of the Standard results (runs dd0.5 to keff2, Table 1), while keff = 0.4J K-1 m-1 s-1 (as in Sakai and others, 2000) gives only 31.2% of the standard result (run keffSakai). However, varying these parameters had almost no effect on discharge calculations, because sub-debris melt accounts for a very small portion of the pond's energy budget. Three parameters controlling the pond surface energy balance were also varied: the water emissiv-ity, ew (affects Lo), roughness height, z0 (affects H and E), and albedo, a (affects In). However, none of these parameters can affect the subaqueous melt rates in the model, and literature values altered discharge results by <5%.

DISCUSSION Geometric considerations

Modelled melt rates are dependent on the geometry assumed (runs slope040, slope0678 and slope080 in Table 1), but pond geometry is difficult to constrain. Few systematic pond-depth readings exist; notable exceptions are Benn and others (2000, 2001), Röhl (2008) and Thompson and others (2012), who report bathymetric surveys for much larger ponds. Basic surveys of several ponds were also reported by Sakai and others (2000), but pond depths were estimated. Sakai (2012) established empirical estimates of volume and maximum depth for terminal lakes in the Himalaya, all of which were much larger than our study pond. Their empirical relations produce smaller volumes and depths than our model run slope040. More work is needed to understand the geometry of supraglacial ponds, which has a strong influence on melt rates.

Our study pond appears in figure 7B of Immerzeel and others (2014), who observed the surface downwasting and velocity for Lirung Glacier for the 2013 monsoon. As our pressure transducer was exposed subaerially before the second UAV flight, the later DEM and lake extent (145.6 m2) cannot be used directly for calibration. Since the lake drained partially between May and October 2013, formerly subaqueous topography was exposed with a mean slope of 41%. The October 2013 DEM was analysed to develop a rating curve, which was directly compared with the 53.2%

slope-derived rating curve, with topography truncated at the October 2013 pond altitude for comparability. For the relevant range of altitudes, the 53.2% slope assumption produces pond volume, surface area and debris area estimates in close correspondence with the revealed terrain in October 2013.

Temperature distribution

Modelled subaqueous bare-ice melt rates are also dependent on the assumed temperature distribution in the pond (model runs Tminuslsig, Tpluslsig and Tlb in Table 1). Note that Tlb is measured, so these model runs calculate identical subdebris melt rates. The temperature profile influences the energy reservoir calculations (Eqn (1)), and consequently has a strong effect on the discharge estimate. Using the observed 2014 monsoon gradients gives night-time pond surface temperatures at the freezing point during the 2013 monsoon for model runs Standard and Tpluslsig, which was not observed in the 2014 monsoon. This is due to the different pond depths during the monsoon periods of 2013 (~2 m) and 2014 (~1 m). A more realistic approach would take the variable air and pond temperature into account, while a computational fluid-dynamics approach would provide the most accurate assessment of temperature distributions and energy exchanges, but at significant computational cost. Our study pond showed similar mean daily surface and pond bottom temperatures to observations by Sakai and others (2000), who ran their model with observed pond bottom temperatures. This approach (model run Tlb) only slightly reduces our calculated melt rates. Few distributed observations of supraglacial pond temperatures have been made (Röhl, 2008), and the melt models are very sensitive to estimated mean temperatures, so additional data are needed to constrain melt estimates and understand pond circulation patterns.

Melt estimates

The subaqueous melt values calculated in this study are lower than those observed on other lakes on debris-covered glaciers, probably due to the small size and particular location of the lake. Benn and others (2001) reported windgenerated currents driving thermo-erosional melt rates between 0.168 and 0.648 md-1, but for a supraglacial lake of 52 500 m2 compared with 600 m2 for our pond. Röhl (2006) identified high rates of thermo-erosion at the water level, although the subaqueous ice cliff was not affected substantially. Sakai and others (2009) identified lake fetch as a key feature in determining wind-driven thermo-erosional melt, which is identified as a criterion for calving in supraglacial lakes. Our pond had average water temperatures of 1-1.5°C, and had a fetch of 20 m in the direction of the dominant up-glacier winds. Following their figures 5-7, Sakai and others (2009) estimate a melt rate of <0.08 md-1 during the monsoon. Our modelled melt rates are even lower, but the pond is hidden from up-glacier winds by an overhanging ice cliff (Fig. 3), making the site suitable for free convection. Most ponds previously studied have had higher temperatures, 2-7°C (Sakai and others, 2009; Xin and others, 2011), which would produce much higher melt rates for all the model runs, and Röhl (2008) suggests a relationship between size and pond temperature. The combined effects of small pond size, shading from the sun and low wind exposure result in lower temperatures and lower thermo-erosional rates for our study pond.

Following the 2013 monsoon period and the associated pond drainage, evidence of waterline melting that had occurred during the pond drawdown was revealed (Fig. 6). These step cuts were several centimetres deep, and occasionally exhibited a lower sill, as observed by Röhl (2006), but were only apparent for 1.85 m above the observed water level in October, which was 3.0 m lower than in May. Lowering during the earlier monsoon period may have cut similar notches, and if so these must have melted subaerially since. It appears that these formed during periods of relative water level stability concurrent with peak melt rates (Fig. 6), melting up to 3 cm in 12 hour periods according to the Standard run, twice the mean melt rate over the study period.

The results summarized in Table 1 show estimated melt in close agreement between melt algorithms, at ~98.5 m3 subaqueous ice melted during the study period. As shown in Figure 5d, the timing of melt differs greatly among the models. Our model (Standard run) predicts a strong diurnal cycle in the melt rate, but a gradual decline as overturning weakens. The other melt models predict much higher peak melt rates during the day but lower melt rates at night. Consequently, the distinct parameterizations for subaqueous ice melt are suited for different scenarios. The model developed in this study is reliable for the free-convection case, where a pond has minimal currents and a monitored temperature gradient, but is likely to underestimate melt if significant currents are present. Performance of the forced-convection algorithm applied by Sakai and others (2000) depends entirely on the selected value of uw. This study shows that uw < 0.02 ms-1 estimates less energy transfer than a simple free-convection model. However, the same algorithm shows great promise for water bodies with longer fetch (Sakai and others, 2009).

The empirical algorithm applied by Röhl (2008) was determined for sea water and probably overestimates melt values significantly. When corrected for fresh water, the algorithm produces melt estimates in close agreement with our model. As it depends solely on an estimated average pond temperature, it is ideal for scaling or for estimating localized melt rates (e.g. for studying thermo-erosional notch development).

Similarly, the Lüthje and Pedersen (2006) algorithm is driven solely by average water temperature, but was developed for natural convection over a horizontal clean-ice surface. Critically, the algorithm assumes a well-mixed water column, which does not hold for our pond or others with thorough temperature measurements (Xin and others, 2011). However, it exhibited good agreement with the free-convection model for our pond's geometry, but is likely to be inaccurate for ponds affected by wind. There is need for a unifying model that performs well under both kinetic settings, estimating wind-induced thermo-erosional melting while also calculating theoretical minimum melt rates based on free convection.

Inferred hydrologic activity

The results for pond discharge and influx indicate that even with a fairly stable water level, the study pond plays an active role in the hydrologic system, in agreement with results of Sakai and others (2000) and Benn and others (2001). The pond is a significant recipient of atmospheric energy, does not seem to increase substantially in temperature, and can only use a small proportion of excess heat to

Fig. 6. (a) Observed notching at the study site during monsoon 2013. (b) Water level decline over the study period, highlighting the water level at the time of notch development. (c) Melt rates corresponding to notch development.

drive local subaqueous bare-ice or subdebris melt. The calculated discharge values may be high based on geometry, temperature and parameter assumptions, but the sensitivity analysis indicates only a small change in mean discharge for most of these variables, so the high discharges are likely to be realistic.

The excess energy is advected from the pond by discharge averaging 0.028 m3s-1 (Standard run), suggesting a low residence time in the pond of 9.2 hours. Field observations in 2014 identified a very small surface outlet meandering to an englacial channel opening that had been exposed due to surface thinning and ice-cliff backwasting, but in 2013 this passage must have been blocked. Assuming discharge occurs at the surface temperature (run TdTws), the low estimated surface temperatures at night produce discharge estimates of up to 2.00 m3s-1. This scenario is unlikely, considering that the daily peak flow at the Lirung Glacier outlet for this and prior periods was 2.5-3 m3s-1 (Bhatt and others, 2007; Ragettli and others, 2015). Instead, the discharge most likely exits in an inefficient manner, accounting for the slow drainage rates, in contrast to the fast lake drainage mechanism envisaged by Gulley and Benn (2007). Inefficient, slow drainage may be associated with flow, through accumulated debris blocking a cut-and-closure conduit, which has backed up water to a temporarily stable higher level (Gulley and Benn, 2009).

Furthermore, the influx water is not sourced locally. The pond's catchment area of 14 800 m2 would have to down-waste 11.2 m on average to supply adequate inflows. The adjacent ice cliff was studied by Buri and others (2015), but modelled melt from this source only accounts for 1-4% of the inflows estimated by this study, while this ice cliff certainly provides the greatest melt signal within the pond's basin (Immerzeel and others, 2014). These clues point to the importance of interactions between englacial conduits and the glacier's surface in determining the development, role and eventual drainage of supraglacial ponds.

Cliff/lake system propagation

For an ice cliff in a combined lake system to stably backwaste, the lake's subaqueous backwasting rate must exceed the subaerial horizontal melt rate, whether via melt or calving (Sakai, 2012). Our study pond demonstrated thermo-erosional notching at the waterline (Fig. 6), in spite of limited pond fetch and low temperature, supporting the assertion by Röhl (2006) that subaqueous cliff melting is controlled additionally by geometry, water fluctuations and debris supply. In our study, subaqueous ice is estimated to have backwasted at an average rate of 2.91 cm d-1, comparable to the 3.25-5.65 cm d-1 (May) and 0.180.23 cm d-1 (October) observed at the adjacent 40-51° cliffs (Steiner and others, 2015).

As suggested by several studies (Sakai and others, 2000; Benn and others, 2001), a significant role of the ponds is to convey atmospheric energy to the glacier's interior. For this study, most energy inputs to the pond are accounted for by the pond's discharge to the glacier's englacial and subglacial conduits, where it is likely to cause rapid melting (Gulley and Benn, 2007; Röhl, 2008). Observations of near-0°C water temperatures at the glacier's terminus outlet and at englacial conduit emergence points indicate that all of the discharge thermal energy is lost to melting. If this is true, our study estimates a total of 2612 m3 of melting in the interior of the glacier solely due to this small pond (over 68 days at a rate of 38.4 m3 d-1), far outweighing the locally caused melt. With an average area of 496 m2 across the study period, this is the equivalent of 5.3 m ablation attributable to the pond's area. As noted by Sakai and others (2000) and Benn and others (2001), this contributes to the formation of new cliff/ lake systems by causing down-glacier conduit collapse and blockage.


Although promising for understanding pond-related energy exchanges, there are limitations to this diagnostic approach.

Temperature gradients have been assigned based on limited observations, and no attempt has been made to model energy dynamics within the pond itself. Pond edge effects are not thoroughly treated in the model. Pond-induced calving is neglected, and although not observed at the field site, this is generally the mechanism of fastest pond expansion (Benn and others, 2001). Rockfall is neglected in our model, but was shown to be a minimal energy input by Sakai and others (2000), and it also displaces volume. Due to a lack of observations, a limited representation of the pond's saturated debris base is used, requiring improved understanding in terms of composition and energy fluxes. In this implementation the model is also dependent on source data and empirical relationships developed outside the study pond.


This study advances energy-balance modelling efforts for supraglacial ponds on debris-covered glaciers, then applies many model configurations for a small pond on Lirung Glacier during May-October 2013 to understand the importance of unconstrained properties and the likely range of melt values. Notably, the pond/atmosphere surface inputs large amounts of energy, with exceptionally high latent fluxes. The net surface energy balance dominates the pond's energy fluxes by an order of magnitude for a variety of parameter choices. The excess energy can only be accounted for by the pond's discharge, and is likely to contribute to substantial amounts of englacial or subglacial melt, the equivalent of 5.3 m local ablation for the ponded area. Therefore, ponds seem to be able to convey a large amount of energy into the glacier interior, demonstrated by this study for a relatively small pond. This distal melt may lead to conduit collapse and the formation of additional cliff/lake systems (Benn and others, 2012).

The study tested several model configurations to calculate subaqueous melt. Subaqueous bare-ice melt was estimated to occur at an average rate of 2.91 cm d-1 for a total volume of 98.5 m3, in the correct range of values to match the adjacent cliff's backwasting. This result is in close agreement with algorithms used by Sakai and others (2000) and Lüthje and Pedersen (2006), and an adaptation of the algorithm used by Röhl (2008). Moreover, the result is also in good agreement with modelled and observed back-wasting of the adjacent ice cliff (Steiner and others, 2015), a precondition for the two systems to occur together (Sakai, 2012). The subaqueous melt algorithms are expected to be suitable for distinct applications based on the dominant kinetic regime driving pond-associated melting (Xin and others, 2011). Subaqueous subdebris melt is unconstrained in the model, but a sensitivity analysis indicates that it plays a minor role at this site in terms of energy loss and melt production (10.6 m3 total) unless the debris is very thin, pond basal temperatures are higher than those observed or convection occurs in the saturated debris layer.

A combination of field measurements and physically based modelling has enabled us to identify several important processes associated with a small supraglacial pond on a debris-covered glacier that are likely to be relevant to other ponds in similar settings. First, water inflow to the pond is not only sourced within the pond's immediate catchment area, but involves a significant supraglacial or englacial input from up-glacier. Second, outflow discharge occurs slowly, suggesting an inefficient outflow channel,

perhaps one blocked by debris. Third, the pond is an active component of the entire glacier's hydrologic system, with up to 10% of the glacier's total discharge passing through it, and resulting in a high overturning rate. Taken together, our results support those of others suggesting that supraglacial ponds are both an important indicator of, and provide a key feedback mechanism for, a debris-covered glacier's response to climate change (Benn and others, 2001, 2012; Sakai and Fujita, 2010). Recent pronounced thinning on Lirung Glacier means that the glacier surface is rapidly approaching base level. This means that the surface regularly intersects former englacial conduits, supplying water to locations where it may be impounded. The surface ponds then absorb atmospheric energy and convey it to the interior of the glacier, leading to englacial conduit enlargement and collapse, further basin formation and the creation of new ponds.


The authors are grateful to Joe Shea and the International Centre for Integrated Mountain Development, Kathmandu, Nepal for providing the Kyanjing AWS data and to Walter Immerzeel for use of the UAV DEM. We thank Markus Holzner for discussions about fluid dynamics and heat transfer. We are indebted to the field assistants and Nepali support staff led by Tek Rai. We thank reviewers Katrin Röhl and Doug Benn, our Scientific Editor, Joe Shea, and the Chief Editor for excellent comments which improved the manuscript. This research was enabled by PhD studentship funding from the Gates Cambridge Trust. Fieldwork was supported by the USAID (United States Agency for International Development) High Mountain Glacier Watershed Programs Climber-Scientist Grant (CCRDCS0010), Swiss National Science Foundation project UNCOMUN (SNF 200021L146761), Trinity College, Cambridge, the B.B. Roberts Fund and the Philip Lake and William Vaughn Lewis Fund.


Basnett S, Kulkarni AV and Bolch T (2013) The influence of debris cover and glacial lakes on the recession of glaciers in Sikkim Himalaya, India. J. Glaciol., 59(218), 1035-1046 (doi: 10.3189/ 2013JoG12J184)

Benn DI, Wiseman S and Warren CR (2000) Rapid growth of a supraglacial lake, Ngozumpa Glacier, Khumbu Himal, Nepal. IAHS Publ. 264 (Workshop at Seattle 2000 - Debris-Covered Glaciers), 177-185 Benn DI, Wiseman S and Hands KA (2001) Growth and drainage of supraglacial lakes on debris-mantled Ngozumpa Glacier, Khumbu Himal, Nepal. J. Glaciol., 47(159), 626-638 (doi: 10.3189/172756501781831729) Benn DI and 9 others (2012) Response of debris-covered glaciers in the Mount Everest region to recent warming, and implications for outburst flood hazards. Earth-Sci. Rev., 114(1-2), 156-174 (doi: 10.1016/j.earscirev.2012.03.008) Bhatt MP, Masuzawa T, Yamamoto M and Takeuchi N (2007) Chemical characteristics of pond waters within the debris area of Lirung Glacier in Nepal Himalaya. Limnology, 66(2), 71-80 (doi: 10.4081/jlimnol.2007.71) Bolch T, Fujita K, Scheel M, Bajracharya S and Stoffel M (2012) The state and fate of Himalayan glaciers. Science, 306(6079), 310-314 (doi: 10.1126/science.1215828) Buri P, Pellicciotti F, Steiner JF, Miles ES and Immerzeel WW (2015) A grid-based model of backwasting of supraglacial ice

cliffs over debris-covered glaciers. Ann. Glaciol., 57(71) (see paper in this issue) (doi: 10.3189/2016AoG71A059) Chikita K and Joshi SP (2000) Hydrological and thermal regimes in a supra-glacial lake: Imja, Khumbu, Nepal Himalaya. Hydrol. Sci. J., 45(4), 507-522 (doi: 10.1080/02626660009492353) Churchill SW and Chu HHS (1975) Correlating equations for laminar and turbulent free convection from a vertical plate. Int. J. Heat Mass Transfer, 18(11), 1323-1329 (doi: 10.1016/0017-9310(75)90243-4) Gardelle J, Arnaud Y and Berthier E (2011) Contrasted evolution of glacial lakes along the Hindu Kush Himalaya mountain range between 1990 and 2009. Global Planet. Change, 75(1-2), 47-55 (doi: 10.1016/j.gloplacha.2010.10.003) Gulley J and Benn DI (2007) Structural control of englacial drainage systems in Himalayan debris-covered glaciers. J. Glaciol., 53(182), 399-412 (doi: 10.3189/002214307783258378) Gulley J and Benn DI (2009) A cut-and-closure origin for englacial conduits in uncrevassed regions of polythermal glaciers. J. Glaciol., 55(189), 66-80 (doi: 10.3189/002214309788608930) Han H, Wang J, Wei J and Liu S (2010) Backwasting rate on debris-covered Koxkar glacier, Tuomuer mountain, China. J. Glaciol., 56(196), 287-296 (doi: 10.3189/002214310791968430) Immerzeel WW, Beek LPH, Konz M, Shrestha AB and Bierkens MFP (2011) Hydrological response to climate change in a glacierized catchment in the Himalayas. Climatic Change, 110(3-4), 721-736 (doi: 10.1007/s10584-011-0143-4) Immerzeel WW and 6 others (2014) High-resolution monitoring of Himalayan glacier dynamics using unmanned aerial vehicles. Remote Sens. Environ., 150, 93-103 (doi: 10.1016/j.rse. 2014.04.025)

Linden PF (2002) Convection in the environment. In Batchelor GK, Moffatt HK and Vorster MG eds. Perspectives in fluid dynamics: a collective introduction to current research. Cambridge University Press, Cambridge, 289-343 Lüthje M and Pedersen LT (2006) Modelling the evolution of supraglacial lakes on the West Greenland ice-sheet margin. J. Glaciol., 52(1 79), 608-61 8 (doi: 10.3189/ 172756506781828386) 0strem G (1959) Ice melting under a thin layer of moraine, and the existence of ice cores in moraine ridges. Geogr. Ann., 41(4), 228-230

Pellicciotti F, Stephan C, Miles E, Immerzeel WW and Bolch T (2015) Mass balance changes of debris-covered glaciers in the Langtang Himal in Nepal between 1974 and 1999. J. Glaciol., 61, 373-386 (doi: 10.3189/2015JoG13J237) Quincey DJ and 6 others (2007) Early recognition of glacial lake hazards in the Himalaya using remote sensing datasets. Global Planet. Change, 56(1-2), 137-152 (doi: 10.1016/j.gloplacha. 2006.07.013)

Ragettli S and 9 others (2015) Unraveling the hydrology of a Himalayan catchment through integration of high resolution in-situ data and remote sensing with an advanced simulation model. Adv. Wat. Resour., 78, 94-111 (doi: 10.1016/j.advwatres.2015. 01.013)

Reid TD and Brock BW (2014) Assessing ice-cliff backwasting and its contribution to total ablation of debris-covered Miage

Glacier, Mont Blanc massif, Italy. J. Glaciol., 60(219), 3-13 (doi: 10.3189/2014JoG13J045) Reynolds JM (2000) On the formation of supraglacial lakes on debris-covered glaciers. IAHS Publ. 264 (Workshop at Seattle 2000 - Debris-Covered Glaciers), 153-161 Röhl K (2006) Thermo-erosional notch development at freshwater-calving Tasman Glacier, New Zealand. J. Glaciol., 52(177), 203-213 (doi: 10.3189/172756506781828773) Röhl K (2008) Characteristics and evolution of supraglacial ponds on debris-covered Tasman Glacier, New Zealand. J. Glaciol., 54(188), 867-880 (doi: 10.3189/002214308787779861) Russell-Head DS (1980) The melting of free-drifting icebergs. Ann.

Glaciol., 1, 119-122 Sakai A (2012) Glacial lakes in the Himalayas: a review on formation

and expansion processes. Global Environ. Res., 16, 23-30 Sakai A and Fujita K (2010) Formation conditions of supraglacial lakes on debris-covered glaciers in the Himalaya. J. Glaciol., 56(195), 177-181 (doi: 10.3189/002214310791190785) Sakai A, Takeuchi N, Fujita K and Nakawo M (2000) Role of supraglacial ponds in the ablation process of a debris-covered glacier in the Nepal Himalayas. IAHS Publ., 264, 119-130 Sakai A, Nishimura K, Kadota T and Takeuchi N (2009) Onset of calving at supraglacial lakes on debris-covered glaciers of the Nepal Himalaya. J. Glaciol., 55(193), 909-917 (doi: 10.3189/ 002214309790152555) Salerno F and 6 others (2012) Glacial lake distribution in the Mount Everest region: uncertainty of measurement and conditions of formation. Global Planet. Change, 92-93, 30-39 (doi: 10.1016/ j.gloplacha.2012.04.001) Scherler D, Bookhagen B and Strecker MR (2011) Spatially variable response of Himalayan glaciers to climate change affected by debris cover. Nature Geosci., 4(3), 156-159 (doi: 10.1038/ ngeo1068)

Shiraiwa T and Yamada T (1991) Glacier inventory in the Langtang

Valley, Nepal Himalayas. Low Temp. Sci., 50, 47-72 Steiner JF, Pellicciotti F, Buri P, Miles ES and Immerzeel WW (2015) Modeling ice-cliff backwasting on a debris-covered glacier in the Nepalese Himalaya. J. Glaciol., 61, (doi: 10.3189/ 2015JoG14J194)

Takeuchi N, Sakai A, Shiro K, Fujita K and Masayoshi N (2012) Variation in suspended sediment concentration of supraglacial lakes on debris-covered area of the Lirung Glacier in the Nepal Himalayas. Global Environ. Res., 16, 95-104 Taylor PD and Feltham DL (2004) A model of melt pond evolution on sea ice. J. Geophys. Res., 109(C12), C12007 (doi: 10.1029/ 2004JC002361)

Thompson SS, Benn DI, Dennis K and Luckman A (2012) A rapidly growing moraine-dammed glacial lake on Ngozumpa Glacier, Nepal. Geomorphology, 145-146, 1-11 (doi: 10.1016/ j.geomorph.2011.08.015) Weeks WF and Campbell WJ (1973) Icebergs as a freshwater

source: an appraisal. J. Glaciol., 12(65), 207-233 Xin W, Shiyin L, Han H, Jian W and Qiao L (2011) Thermal regime of a supraglacial lake on the debris-covered Koxkar Glacier, southwest Tianshan, China. Environ. Earth Sci., 67(1), 175-183 (doi: 10.1007/s12665-011-1490-1)