Modelling runoff from a Himalayan debris-covered glacier

Although the processes by which glacial debris mantles alter the melting of glacier ice have been well studied, the mass balance and runoff patterns of Himalayan debris-covered glaciers and the response of these factors to climate change are not well understood. Many previous studies have addressed mechanisms of ice melt under debris mantles by applying multiplicative parameters derived from field experiments, and other studies have calculated the details of heat conduction through the debris layer. However, those approaches cannot be applied at catchment scale because distributions of thickness and thermal property of debris are heterogeneous and difficult to measure. Here, we established a runoff model for a Himalayan debris-covered glacier in which the spatial distribution of the thermal properties of the debris mantle is estimated from remotely sensed multi-temporal data. We applied the model to the Tsho Rolpa Glacial Lake–Trambau Glacier basin in the Nepal Himalaya, using hydro-meteorological observations obtained for a 3.5year period (1993–1996). We calculated long-term averages of runoff components for the period 1980–2007 using gridded reanalysis datasets. Our calculations suggest that excess meltwater, which implies the additional water runoff compared with the ice-free terrain, from the debris-covered area contributes significantly to the total runoff, mainly because of its location at a lower elevation. Uncertainties in runoff simulation due to estimations of the thermal properties and albedo of the debris-covered surface were assessed to be approximately 8 % of the runoff from the debris-covered area. We evaluated the sensitivities of runoff components to changes in air temperature and precipitation. As expected, warmer air temperatures increase the total runoff by increasing the melting rate; however, increased precipitation slightly reduces the total runoff, as ice melting is suppressed by the increased snow cover and associated high albedo. The response of total runoff to changing precipitation is complex because of the different responses of individual components (glacier, debris, and ice-free terrain) to precipitation.


Introduction
Glaciers are considered to play an important role as the water resources for densely populated Asian regions (e.g.Cruz et al., 2007).Recent studies have revealed that the response of glaciers to climate variations varies considerably in Asian highland regions (e.g.Fujita and Nuimura, 2011;Bolch et al., 2012;Yao et al., 2012), and that the response depends partly on the characteristics of the debris mantles on Himalayan glaciers (Scherler et al., 2011).Terminus positions of heavily debris-covered glaciers seem to be insensitive to changes in climate (Scherler et al., 2011), while surface lowering over debris-covered areas seems to be comparable to that in debris-free ablation areas (Nuimura et al., 2011(Nuimura et al., , 2012;;Kääb et al., 2012).It is still unclear whether heterogeneity in climatic forcing or debris cover pattern is responsible for observed temporal variations in glacial melt observed in different Himalayan glacier systems.Experimental studies have revealed that thin debris layers accelerate the melting of underlying ice, whereas thick debris layers suppress melting (e.g.Østrem, 1959;Mattson et al., 1993).Some numerical simulations of conductive heat flux through the debris layer have successfully reproduced patterns of ice melting under the debris layer (e.g.Nicholson and Benn, 2006;Reid and Brock, 2010).However, these heat conduction models cannot be applied to a basin-scale mass balance calculation of debris-covered glaciers because the spatial distributions in debris thickness and thermal conductivity are nearly impossible to measure.On the other hand, some hydrological studies in glacierized catchments containing debris-covered glaciers have parameterized ice melting under the debris layer (e.g.Lambrecht et al., 2011;Immerzeel et al., 2012;Juen et al., 2014).Although these studies have been validated by hydrologic and/or other observational data, continuity in surface conditions over time cannot be guaranteed, especially in systems with rapidly changing glaciers.In addition, the debris-covered surfaces of real glaciers exhibit highly heterogeneous and rugged topography, over which no representative thickness is obtainable.Heat absorption in such rugged topography, which includes ice cliffs and supraglacial ponds, is considered to be one of the significant sources of heat for melting in debris-covered areas (Sakai et al., 2000a(Sakai et al., , 2002)).Therefore, prediction of basin-scale patterns of ice melt on debris-covered glaciers from a simple relationship between debris thickness and ice melting is exceedingly difficult.
To overcome the difficulties discussed above, we have adopted the "thermal resistance" parameter proposed by Nakawo and Young (1982).This parameter is defined as the debris thickness divided by the thermal conductivity of the debris layer and its spatial variations may be obtained from remotely sensed data, such as data obtained from Landsat or ASTER imagery.Nakawo and Rana (1999) used this approach to estimate the distribution of thermal resistance on glaciers from Landsat TM data, and successfully reproduced runoff from the debris-covered Lirung Glacier in the Langtang region of Nepal.Subsequently, Suzuki et al. (2007) demonstrated temporally consistent values of thermal resistance on glaciers in the Bhutan Himalaya, as determined from ASTER data taken on different dates, for which surface temperature and albedo were calibrated using field measurement conducted at the same time as ASTER acquisitions.Zhang et al. (2011) obtained the thermal resistance distribution of a debris-covered glacier in southeastern Tibet and validated the calculated thermal resistance, melt, and runoff with in situ measurements.However, these studies did not evaluate uncertainties in thermal resistance values, or how these affect both the calculated ice melt under the debris and the resulting runoff.In this study, therefore, our goal was to obtain thermal resistance values and to evaluate uncertainties in the values based on ASTER data acquired in different seasons and years.In addition, we establish an integrated runoff model that incorporates variations in surface conditions, such as debris-covered and debris-free glacier surfaces as well as ice-free terrain.Model performance was tested for a catchment with a debris-covered glacier in the Nepal Himalaya.We evaluated and discussed the uncertainties associated with thermal resistance and albedo, and the sensitivity of runoff to meteorological variables.

Location, data and models
Abbreviation, unit and value of all parameters used in this study are summarized in Table 1.

Delineation and classification of the catchment
We chose the Tsho Rolpa Glacial Lake-Trambau Glacier basin located at the head of the Rolwaling Valley, in the east Nepal Himalaya (27.9 • N, 86.5 • E; Fig. 1) as our study site.Tsho Rolpa is one of the largest glacial lakes in the Nepal Himalaya.We delineated the basin using a digital elevation model produced from multi-temporal ASTER data (ASTER-GDEM, 2009;Tachikawa et al., 2011).The basin extends from 4500 to 6850 m a.s.l., with a total area of 76.5 km 2 (Fig. 1a and Table 2).
We divided the surface features of the basin into four categories: debris-covered glacier (debris), debris-free glacier (glacier), ice-free terrain (ground), and lake surface (Tsho Rolpa) to perform the following runoff calculations.Using the clearest available ASTER image acquired in February 2006 (Fig. 1a), we calculated the normalized difference water index (N W ) and normalized difference snow/ice index (N S ) from reflectance of the ASTER sensors (r n ) using the following equations: The N W has been successfully used to delineate glacial lake boundaries in the Himalayas (Fujita et al., 2009).The N S has been used to evaluate snow cover extent in North America (Hulka, 2008).Thresholds of N W and N S are assumed to be 0.42 and 0.94, respectively, to best distinguish the surfaces.Debris-covered surface was visually distinguished from icefree terrain using surface morphology such as rugged relief and ice flow features (Nagai et al., 2013).Steep slope terrain without snow or ice (steeper than 30 • ) was also defined as ice-free terrain.The resulting basin surface category map is shown in Fig. 1b, and the hypsometry (area-altitude profile) based on the ASTER-GDEM is shown in Fig. 2.

Meteorological and hydrological data
Air temperature, solar radiation, relative humidity, and wind speed are required as input variables for models in this study.Meteorological, hydrological, and limnological observations were conducted in the 1990s (Fig. S1 in the Supplement; Yamada, 1998;Sakai et al., 2000b).The observations are used to confirm the plausibility of the gridded data and to validate the calculated runoff, for which we use gridded data as model inputs to examine the long-term mean and seasonal cycle of runoff components.Air temperature, solar radiation, relative humidity, and wind speed are taken from the NCEP/NCAR reanalysis gridded data (NCEP-1, Kalnay et al., 1996).Air temperature at the elevation of the observation site (4540 m a.s.l.) is linearly interpolated from air temperatures at geopotential heights of 500 and 600 hPa; the temperature lapse rate is also obtained from these data.Wind speed at a 2 m height from the surface (U ) is estimated from 10 m wind in the reanalysis data (U 10 ), based on the assumption of    a logarithmic dependence of wind speed on height: where surface roughness (z 0 ) is assumed to be 0.1 m.The ground-based Aphrodite daily precipitation data are used, which have a spatial resolution of 0.5 • × 0.5 • (Yatagai et al., 2009).All variables except for wind speed show significant correlations between gridded and observational data (Fig. S2 in the Supplement).Air temperature shows a partic- ularly high linear correlation, with little bias.The statistical parameters strongly support representativeness of the temperature, though temperature generally shows a good consistency among in situ and reanalysis data because of the climatic seasonality.Solar radiation, relative humidity, and wind speed show less significant or no correlations.Fujita and Ageta (2000) have pointed out that uncertainties in these variables are less important for the mass balance of Tibetan glaciers than those of air temperature and precipitation.Although correlation of pentad (5 day) precipitation is weaker than that of air temperature, it is still significant (Fig. S2 in the Supplement).We therefore use the gridded data for all variables except for precipitation, and compare modelled and observed runoffs to find the best set of calibration coefficients using the Aphrodite precipitation data and elevation corrected precipitation (Sect.3.2).

Thermal resistance
Thermal resistance is defined as debris thickness divided by the thermal conductivity of the debris layer (Nakawo and Young, 1982).Suzuki et al. (2007) established a methodology to obtain the thermal resistance distribution from ASTER and reanalysis climate data.Zhang et al. (2011) confirmed that the distribution of thermal resistance was well correlated with that of debris thickness from in situ measurements over a southeastern Tibetan glacier with a rather gentle and homogeneous debris-covered surface.We obtained the thermal resistance of the debris-covered area from multitemporal ASTER data following their methods.The thermal resistance (R T ) is defined by debris thickness (h) and thermal conductivity (λ) as Assuming no heat storage in the debris layer, no heat conduction into temperate glacier ice, and a linear temperature profile within the debris layer, the conductive heat flux through the debris layer (G d ) is described with the surface temperature (T s ) and the temperature at the interface between debris and ice (T i ), which is assumed to be melting point, as The conductive heat flux from the surface toward the debrisice interface is described as a residual term of the heat balance at the debris surface, according to All components are positive when the fluxes are directed towards the ground.Although turbulent heat fluxes have to be taken into account in the exact heat exchange over the debris surface, Suzuki et al. (2007) demonstrated that these fluxes are negligible because there is only limited mass flux in the low density air at Himalayan high elevation.Clear sky conditions, which are required for satellite data utilization, are also associated with a reduced importance of turbulent heat fluxes, especially of latent heat.We therefore assumed that the turbulent heat fluxes were zero (H S = H L = 0).We can then obtain the thermal resistance at a given point without knowing the debris thickness and thermal conductivity if we know the downward short-wave and long-wave radiation fluxes, the albedo, and the surface temperature.
We selected eight cloud-free images of ASTER level 3A1 data, which is a semi-standard ortho-rectified product available from ERSDAC Japan (Table S1 in the Supplement).Surface albedo is calculated using three visible near-infrared sensors (VNIR; bands 1-3) using the equations described in Yüksel et al. (2008).Surface temperature is obtained from an average of five sensors in the thermal infrared (TIR; bands 10-14) using the formula proposed by Alley and Nilsen (2001).The spatial resolution of the thermal resistance is then constrained by the coarsest resolution of the ASTER TIR sensors (90 m).We utilize NCEP/NCAR reanalysis 6hourly data (Kalnay et al., 1996) for both downward radiation fluxes at the timing closest to ASTER acquisition.

Snowmelt and albedo
Heat balance over the snow surface (Q s ) and daily snowmelt (M s ) are estimated as The turbulent fluxes (H S and H L ) are estimated by bulk formulae as A wetness parameter (τ w ) is assumed to be 1 for the snow and ice surfaces, while it varies over debris-cover.Snow surface albedo on a given day (α day ) is calculated with a scheme proposed by Kondo and Xu (1997), in which an exponential reduction of snow albedo with time after a fresh snowfall is assumed as The number of days after the latest fresh snow date is set to zero (day = 0) when snowfall is greater than 5 mm w.e., and the albedo of firn (α f ) is taken as the minimum snow albedo (0.4).The parameter k depends on air temperature, according to The albedo of the initial fresh snow (day = 0) also depends on air temperature: Surface albedo is affected by the glacier ice or debris surface if the snow layer is thin.According to Giddings and LaChapelle (1961), the penetration of solar radiation into snow is assumed to follow Fick's second law of diffusion with a term for simultaneous absorption, and the surface snow albedo (α s ) over the underlying ice or debris surface Extinction coefficient of snow (K) is assumed to be 30 m −1 (Greuell and Konzelmann, 1994), and albedo of the underlying surface (α b ) is taken to be ice (α i ) or debris (α d ) based on the targets.The albedo of glacier ice (α i ) is assumed to be 0.2 based on our field observations on Asian glaciers (Takeuchi and Li, 2008;Fujita et al., 2011).

Probability of snow and rain
Precipitation across the Himalayan regions takes place mainly during the summer monsoon season so that the precipitation phase (snowfall or rainfall) has to be taken into account.Based on observational reports in Tibet (Ueno et al., 1994;Sakai et al., 2006a), we assume the probability of snowfall (P s ) and rainfall (P r ) to depend on air temperature as follows: P r = P p − P s . (16)

Energy and mass balance of the debris-free glacier
Energy and mass balance over the debris-free glacier surface are calculated in 50 m elevation bands from a model established by Fujita and Ageta (2000), which has successfully calculated the glacier mass balance, equilibrium line altitude and runoff of several Asian glaciers (e.g.Fujita et al., 2007Fujita et al., , 2011;;Sakai et al., 2009aSakai et al., , 2010;;Fujita and Nuimura, 2011;Zhang et al., 2013).The basic equation can be written as Downward long-wave radiation (H LR ) is estimated from an empirical equation using air temperature, relative humidity and the ratio of solar radiation to that at the top of atmosphere based on Glover and McCulloch (1958) and Kondo (1994).
We determine the surface temperature by iterative calculations, in which the conductive heat flux into the glacier ice (G g ) is calculated by changing the ice temperature profile.Daily runoff water (D g ) is obtained as Here refrozen ice in the snow layer (R f ) is obtained from the change in the ice temperature profile when surface water is present.All details are described in Fujita and Ageta (2000) and Fujita et al. (2007).

Energy and mass balance of debris-covered surface
We calculate heat balance at the debris-covered surface using Eq. ( 6), but also take into account the turbulent heat fluxes to give results valid for various weather conditions, such as clear, cloudy, rainy, and snowy conditions, though these were neglected when the thermal resistance was obtained under the clear sky assumption (see Sect. 2.3).We use an alternative bulk coefficient (C d ) for the turbulent heat fluxes over the debris surface in Eqs. ( 9) and ( 10).In addition, we assume that the wetness parameter (τ w ) for the latent heat flux in Eq. ( 10) changes with the thermal resistance (R T ) as because Suzuki et al. (2007) revealed that the debris surface was wet (τ w ≈ 1) when its thickness was thin and became exponentially drier (τ w ≈ 0) with increased thermal resistance in the Bhutan Himalaya.We determine the surface temperature that satisfies Eq. ( 6) by iterative calculation.Once the surface temperature is determined and the heat flux toward the ice-debris interface is positive, the daily melt of ice beneath the debris layer (M d ) and then daily runoff (D d ) generated at a given point are obtained as It is assumed that all heat flux into the debris layer is used to melt ice.Condensation of vapour is also taken into account if available (H L / l e ) though it is generally negligible in many cases.If seasonal snow covers the debris surface, no ice melt beneath the debris layer is assumed until the snow cover completely melts away.Daily runoff in the presence of snow is thus obtained from Eq. ( 21), but with the ice melt beneath the debris-layer (M d ) replaced by the snowmelt (M s ) in Eq. ( 8).
The spatial resolution for the debris-covered surface is 90 m, which is constrained by the ASTER TIR data used to obtain the surface temperature in the thermal resistance calculation (Sect.2.3).

Runoff from ice-free terrain and the lake
Runoff from ice-free terrain is calculated for 50 m elevation bands, based on a simple bucket model proposed by Motoya and Kondo (1999).The potential evaporation rate (E p ) is obtained from the energy balance: Here the albedo of ice-free terrain (α w ) is assumed to be 0.1.Evaporation efficiency (β) depends on soil moisture content: which is expressed as a ratio of water contents (W a ) in the surface storage just below the surface to the maximum water content (W amax ) (Fig. 3).The bulk coefficient (C t ) is parameterized with wind speed (U ) as Runoff from the ice-free terrain (D t ) is obtained when the surface storage is full: If there is snow cover, snowmelt (M s ) is calculated using Eqs.( 7) and ( 8), in which direct liquid condensation is taken into account if available.If there is no snow, evaporated water (βE p ) is reduced from the rainwater value.Water is first used to fill the surface storage capacity (W amax -W a ) in all cases.If there is insufficient water to fill the surface storage, no runoff is generated (D t = 0) and the water content in the next time step (W n ) is given by If evaporation is greater than the sum of rain and water content, evaporated water is constrained by the water in the surface storage (βE p = W a ) and no water content is expected in the next time step (W n = 0).We have little information on the water balance of the Tsho Rolpa Glacial Lake, although the water circulation within the  Motoya and Kondo, 1999).Surface storage is used to calculate energy and water balance of the ice-free terrain (see the Sect.2.3.3).Internal and ground storages are used to calculate final daily runoffs for the individual components such as the debriscovered surface (R d ), the debris-free glacier (R g ), the ice-free terrain (R t ), and the lake (R l ).lake has been thoroughly investigated (Sakai et al., 2000b).Therefore we assumed that precipitation would immediately be removed as runoff from the lake (D l ), giving the maximum runoff without evaporative loss because the outlet is located just below the lake, though the water surface should be a significant source of evaporation.

Bucket model calculating river runoff
All water generated over the debris-covered part (D d ), debris-free snow or ice (D g ), ice-free terrain (D t ) and the lake (D l ) is added to the river system through two types of storage, internal and ground storage (Motoya and Kondo, 1999).A schematic diagram that also includes the surface storage for ice-free terrain is shown in Fig. 3.The surface water inflow (F a ), which is made up of the individual surface water inflows (D d , D g , D t , D l ), is added to the internal storage.Outflow from the internal storage (F b ) will occur and be directly added to the final runoff when the volume of water stored (W b ) exceeds the maximum capacity (W bmax , 500 mm) according to Leakage from the internal storage (F c ) is simultaneously calculated as We here assume that 30 % of the internally stored water will be lost in a day (k b = 0.3).Part of the leakage from the internal storage will be directly added to the final runoff and the rest will flow into the ground storage (W c ).There is no limit on the capacity of the ground storage.Leakage from the ground storage (F d ) is given by This flow, which is assumed to be 3 % of the ground storage (k c = 0.03) will form the continuous basal flow of the river system.We obtain the final runoff (F f ) as The fraction (r c ) is assumed to be 0.8.The final runoff can be calculated for individual runoffs from the debris-covered surface (R d ), the debris-free glacier (R g ), the ice-free terrain (R t ), and the lake (R l ).We summarize these runoffs by considering a debris grid with 90 m resolution and the hypsometry of the debris-free glacier surface and ice-free terrain in 50 m elevation bands (Fig. 2).

Distributions of thermal resistance and albedo
We calculated the distribution of thermal resistance from eight ASTER images (Figs. 4 and 5).Some images showed a plausible distribution of thermal resistance (Fig. 4) but a fragmented distribution was obtained in winter images (Fig. 5).
Because the ice-debris interface is assumed to be at the melting point temperature in the calculation of thermal resistance (T i in Eq. 5), it may not be possible to calculate the thermal resistance under cold winter conditions.We therefore obtain an average distribution of the thermal resistance from the four plausible distributions as shown in Fig. 1b.Where calculations were not possible for the debris-covered part, as shown by grey shading in Fig. 1b and at higher elevations, zero thermal resistance is assumed, implying a debris-free glacier.We assumed that the topographical features remained unchanged through the simulation period 1979-2007.
Comparisons of individual thermal resistances against the average show some degree of variability (Fig. 6a).A linear regression of standard deviation against the average suggests that the thermal resistance has an uncertainty of 30 % (Fig. 6c).We simultaneously obtain a distribution of surface albedo, which is required to calculate the thermal resistance and complete the energy mass balance model of the debriscovered surface.Although one image taken in October 2004 shows rather large scatter (Fig. 6b), the uncertainty in albedo expressed as a standard deviation is of a similar level to that of thermal resistance (Fig. 6d).We evaluate the influences of these uncertainties on runoff from the debris-covered surface later (Sect.4.1).

Validation
A 1-year cycle of the calculation runs from 1 October to 30 September of the next year.We first conducted a 4-year calculation from 1 October 1992 to 30 September 1996, and compared the results with the observed runoff at the outlet of the Tsho Rolpa Glacial Lake (shown as a yellow cross in Fig. 1a).Because the reanalysis air temperature represents the observations well (Fig. S2 in the Supplement), we seek the best set of precipitation ratios relative to the Aphrodite precipitation and elevation gradient of precipitation to produce the best estimate of total runoff.We calculated both the root mean square error (D RMS ) and the Nash-Sutcliffe model efficiency (E N ) of the simulation against the observed runoff (Nash and Sutcliffe, 1970).We found that the best estimation was obtained along an isoline of the precipitation ratio of 74 % against the original Aphrodite precipitation averaged over the whole basin (Fig. 7).We adopt 55 % as the precipitation ratio and 35 % km −1 as the elevation gradient of precipitation for the subsequent analysis (thin dashed lines in Fig. 7) based on the comparison of precipitation (Fig. S2 in the Supplement), and elevation gradients of precipitation observed in another Himalayan catchment (Seko, 1987;Fujita et al., 1997).Daily runoff is well reproduced for the 3 hydrological years (Fig. 8).We also performed the calculation using gap-filled meteorological variables without assuming an elevation gradient of precipitation, for which the original observed data were used where available.We obtained  1a), as a function of precipitation ratio (horizontal axis) against the original Aphrodite precipitation and elevation gradient of precipitation (vertical axis).We adopt 55 % as the precipitation ratio and 35 % km −1 as the elevation gradient of precipitation for subsequent analysis (thin dashed lines).The thick dashed line denotes the 74 % precipitation ratio isoline for the whole basin.
similar values of D RMS and E N at the precipitation ratio of 75 %.This implies that reanalysis gridded data are useful to drive the models if the temperature representativeness is sufficiently good and precipitation data are calibrated accordingly.
In Fig. 8, the model overestimates the runoff at the beginning of the melting season.This discrepancy could be caused by model settings in which the generated meltwater was immediately put into the internal storage (Fig. 3).At the beginning of melting season, meltwater could be retained within the snowpack (Gao et al., 2012) or internal channels of glacier.In addition, the lake could be a strong buffer to cause runoff delay when lake level rose.

Long-term averages
We further calculated the average value of each component in the long term to understand the present condition of the basin.We calculated daily runoff in the period 1979-2007 (28 hydrological years) and then obtained the seasonal cycle (Fig. 9) and annual average (Table 2).We assumed that the geometry and surface condition of the basin were unchanged in the calculation, though expansion of the glacial lake should have supplied excess meltwater in the runoff by calving of the glacier front.Runoff contribution and seasonal cycle show that runoff from the debris-covered surface accounts for more than half of the total runoff (55 %).Comparing area ratio and runoff contribution, the ice melt beneath debris cover supplies significant excess water to the total runoff, which implies the additional water runoff compared with the ice-free terrain (Table 2).This is clearly shown in terms of runoff depth.Both annual averages (Table 2) and seasonal cycle (Fig. 9b) suggest that the debris-covered area yields runoff depth approximately seven times greater than from the debris-free glacier surface or ice-free terrain.Both runoff depths from the debris-free glacier surface and icefree terrain are slightly less than that due to precipitation because of evaporative loss.The similar runoff depths of debrisfree glacier surface and ice-free terrain suggest that the entire debris-free glacier is in a state of balanced budget (Table 2).

Uncertainty caused by thermal resistance and albedo of debris cover
In earlier work on thermal resistance, Suzuki et al. (2007) calibrated surface temperature and albedo with field measurements performed at the same time as ASTER acquisitions in the Bhutan Himalaya.In this study, however, we have no in situ data for calibration so that the reanalysis data are used without adjustment.Our thermal resistances have greater scatter than those of Suzuki et al. (2007;Fig. 5) probably because uncalibrated data were used.On the other hand, Zhang et al. (2011) obtained the thermal resistance distribution of a debris-covered glacier in southeastern Tibet from a single ASTER image.They validated the thermal resistance and melt calculations with their in situ measurements of debris thickness and melt rate.However, these studies did not evaluate how the uncertainty of thermal resistance affects the calculated melt under the debris or runoff.We therefore calculated the influence of uncertainties in thermal resistance and albedo (Fig. 6c and d).At some points there are insufficient data to calculate the standard deviation of thermal resistance, so for such points we use an estimate from the linear regression (Fig. 6c).The standard deviation of the albedo was obtained for all points so that the regression curve was not used (Fig. 6d).We obtain the runoff anomaly associated with changed parameter (dv) from the control calculation in Sect.3.3 by averaging positive and negative cases according to where R, v and δv denote calculated runoff, the parameter used for the control calculation (thermal resistance or albedo) and its standard deviation, respectively.Changes in albedo (R T ave , dα) or thermal resistance (dR T , α ave ) reduce the debris runoff (Table 3).Uncertainty due to albedo (−8 %) has a slightly larger effect than that due to thermal resistance (−5 %).The simultaneous change of both parameters (dR T , dα) results in an additive impact on the debris runoff (−13 %).+δα) suggest that the uncertainty due to albedo has more effect on the debris runoff than that due to thermal resistance.
In the absence of calibration data, the use of multiple ASTER images to derive thermal resistance and albedo gives a runoff uncertainty within 8 %.

Uncertainty caused by other parameters
Because the value of parameters was used, as the original studies proposed in this study (Table 1), relevant uncertainties were masked by the corrected precipitation.Here we performed a sensitivity analysis by changing the parameter by ±30 %, which was equivalent to the degree of uncertainty in the thermal resistance and albedo of debris cover.We averaged anomalies using Eq. ( 33) and then expressed by percentage to the control values (Table 4).The Nash-Sutcliffe model efficiency coefficient (E N ) was obtained against the control calculation.Albedo of firn (α f ) seems to alter the calculated runoff significantly, while those of glacier ice (α i ) or ice-free terrain (α w ) are less influential.In addition to a wider range of change in the firn albedo (±0.12), rather than those in glacier ice (±0.06) and ground (±0.03) by the 30 % perturbation, timing of disappearance of snow cover could be significantly altered by the setting of firn albedo.Bulk coefficient for the snow-ice surface (C s ) is slightly more influential than those for the debris surface (C d ) and for the ice-free terrain (C t ), probably because of the same reason mentioned above.Assumption of wetness for the debris surface (τ w ) was tested not by changing 30 %, but by two extremes; no water (τ w = 0) and water surface (τ w = 1).Among these two extreme settings, uncertainty due to the wetness parameter on debris runoff is less than half of those due to thermal resistance and albedo (Table 3).Although parameter settings in the bucket model (W bmax , k b , k c , r c ) seem not to affect the runoff amount, a low value of the Nash-Sutcliffe model efficiency coefficient suggests that change in the fraction for leakage from the internal storage (r c ) could alter the seasonal cycle of runoff.As a whole, boundary conditions such as thermal resistance and albedo are key parameters in the runoff modelling for the Himalayan debris-covered glacier.

Effects of debris cover, lake and glacier
It is clear that the debris-covered area supplies significant excess meltwater to the total runoff (Table 2 and Fig. 9).Elevation profiles of debris-covered and debris-free surfaces suggest comparable runoffs from both surfaces (Fig. 10a), so that the significant excess meltwater may be attributed to the lower elevation of the debris-covered area (Fig. 2).To evaluate whether the excess meltwater is generated by accelerated melt due to thinner and darker debris cover or by the lower elevation of the debris-covered area, we performed a sensitivity calculation by assuming no debris cover (the no debris assumption in Table 3).Compared with the control calculation, a debris-free surface would yield slightly less water (−3 %), implying that the significant excess water is generated mainly by the lower elevation of the debris-covered area, and is slightly increased by the acceleration effect of thin and dark debris.This is the opposite of that estimated for the Lirung Glacier in the Langtang region, Nepal (Nakawo and Rana, 1999), where the debris cover significantly suppressed the melting of ice underneath.In fact, the regional distribution of thermal resistance suggested that the debris cover over the Trambau Glacier was thinner than that over the other Note: "N/A" denotes that no specific value is available because of given boundary conditions or parameterization."-" denotes that parameter setting does not affect result.Influence of wetness parameter was obtained by two extreme cases (0 and 1).Upper bound of fraction for the leakage from the internal storage was set at 1.0.
glaciers in the Khumbu and neighbouring regions (Suzuki et al., 2007).Suzuki et al. (2007) pointed out that not only the Trambau Glacier but also other glaciers having glacial lakes tended to have thinner debris-cover in terms of thermal resistance.Sakai and Fujita (2010) also demonstrated that glacial lakes in the Nepal and Bhutan Himalayas were found at the termini of debris-covered glaciers, at which thinning since the Little Ice Age was greater than 50 m and for which the slope was gentler than 2 • .Although it is still unclear whether the thinner debris triggered the glacial lake formation or thick debris portion has turned into the glacial lake, the topographical setting should affect the debris-covered area through debris supply (Nagai et al., 2013).Recent studies have revealed that the thinning rates of debris-covered surfaces were comparable to those of debris-free surfaces in the Himalayas (Nuimura et al., 2011(Nuimura et al., , 2012;;Kääb et al., 2012).Because the surface thinning of a glacier is affected not only by surface melting but also by dynamics (the degree of compressive flow in the ablation area), we cannot naively attribute the significant thinning of the debris-covered surface to the comparable melt of debris-covered ice (Sakai et al., 2006b;Berthier and Vincent, 2012).Nevertheless, the significant melting of ice under the debris layer would be one of the reasons for the significant thinning of debris-covered areas in the Himalayas.Another surface feature of the basin is the Tsho Rolpa Lake, one of the largest glacial lakes in the Nepal Himalayas.In our calculation, the size of the lake is assumed to be constant, but the lake has expanded since the 1950s at a rate of 0.03 km 2 year −1 (Yamada, 1998;Komori et al., 2004;Sakai et al., 2009b).We therefore performed another sensitivity calculation without a lake.The thermal resistance (0.0151 m −1 K −1 W), albedo (0.230) and elevation (4560 m a.s.l.) over the lake are taken from values at the lowermost part of the debris-covered area (approximately 500 m from the glacier terminus).Runoff from the debris-covered surface and the lake significantly increases by 22 % from the control (Table 3).This suggests that ice located at a lower elevation is the main source of excess meltwater in the basin, and the meltwater might have decreased with expansion of the lake.
Glaciers are recognized as water resources in the Asian highlands.The disappearance of glaciers is projected to result in severe depletion of river water that will threaten human life (e.g.Cruz et al., 2007).In this regard, the contribution of glacier meltwater to river runoff has been evaluated in a number of studies (e.g.Immerzeel et al., 2010;Kaser et al., 2010).However, considering that precipitation could still be expected over the terrain after the glaciers disappear, the reduction in glaciers would not directly result in such a severe depletion of river runoff.A future runoff projection for a Himalayan catchment demonstrated that increased precipitation and seasonal snowmelt would compensate for the decrease in glacier meltwater (Immerzeel et al., 2012).We therefore simply assumed that runoff depth over the ice-free terrain (941 mm) was applicable to the whole basin and then evaluated the runoff under the no ice environment (Table 5).Because the excess meltwater is added to the control total runoff, the no ice assumption results in a significant runoff reduction of 43 %.Although an increase in evaporation, which is expected under the climatic conditions, resulting in the disappearance of a glacier, is not taken into account, river water will still be available from the basin.In terms of future water availability, more uncertainty will be caused by projected changes in precipitation.

Sensitivities
To understand how the basin consisting of debris-covered glaciers responds to changes in climatic variables such as air temperature and precipitation, we calculated the runoff  sensitivities by altering the annual air temperature and precipitation, by ±0.1 • C for air temperature or ±10 % for precipitation from the control conditions.We assumed no topographical change here, though glacier extent and surface features would change as responses to long-term climate change.The runoff anomaly is obtained in the same way as uncertainty (Sect.4.1): by averaging positive and negative cases for which signs are taken into account (Eq.33).
Warmer air temperature significantly increases the melting of ice over both debris-covered and debris-free surfaces and thus total runoff, while increased evaporative water loss over the ice-free terrain is negligible (Table 6).Elevation profiles of the response to the warming show no significant difference between debris-covered and debris-free surface at a given elevation (Fig. 10b).The doubled sensitivity of runoff depth over the debris-covered area should be attributed to its lower elevation as discussed in Sect.4.2 (Table 6).Because the debris-free glacier surface mainly consists of the high accumulation zone reaching to above 6000 m a.s.l., warming will have a limited impact overall.On the other hand, an increase in precipitation will potentially prolong the duration of high albedo snow cover, which suppresses the absorption of solar radiation, and will result in a runoff reduction from ice, while runoff from the ice-free terrain and the lake will increase with precipitation (Table 6).These opposing responses compensate for each other and thus result in a smaller influence on total runoff than that caused by warming.The elevation profiles of response show greater sensitivity over the debrisfree glacier than over the debris-covered ice, which may be caused by different albedo settings for glacier ice and debris (Fig. 10b).These sensitivities to changes in air temperature and precipitation cannot simply be compared directly.Considering the standard deviations of air temperature (0.47 • C) and precipitation (93 mm, 9.4 %) for the period 1979-2007 (28 years), we obtain the variability in runoff associated with the variability in climatic variables (Table 6).Variability of the total runoff caused by air temperature variability is 23 times greater than that caused by precipitation variability, though this result could be variable if we used a different time span.The small changes applied above simply result in a linear response, but we further tested runoff sensitivity to precipitation by changing the precipitation over a wider range, from 40 to 200 % of that used in the control calculation.Runoffs from the ice-free terrain and the lake respond linearly to changing precipitation in proportion to their areas, while those from debris-covered and debris-free surfaces respond non-linearly (Fig. 11a).In particular, a deficit of precipitation will yield extreme ice melt because it gives a dark surface without snow cover (blue line in Fig. 11a).Glacier runoff will become stable under conditions of extreme humidity because of compensation between suppressed ice melting and increased rain water.Summing components with different sensitivities results in a complicated total runoff response (black line in Fig. 11a).This suggests that present climatic and topographic conditions of the target basin have the smallest sensitivity to changing precipitation.If the precipitation regime changes significantly for a long period of time, runoff would respond more significantly than under the present regime, though the glacier extent would also change with time.Seasonal cycles of runoff under the two extreme conditions (50 and 200 % of the control) show impressive responses (Fig. 11b).A wetter climate simply increases the runoff during the humid monsoon season, which is affected by precipitation seasonality (blue line in Fig. 11b) while a drier climate significantly alters the seasonal cycle of runoff (orange line in Fig. 11b).Reduced precipitation will accelerate ice melting in the spring as the dark ice surface uncovered by high albedo snow.Although such an effect may not be obvious in the sensitivity obtained by changing precipitation over a small range (±10 % for instance), a change in precipitation could potentially alter the seasonality of runoff, which is important for regional water availability (e.g.Kaser et al., 2010).

Conclusions
We have developed an integrated runoff model, in which the energy-water-mass balance is calculated over different surfaces such as debris-covered surface, debris-free glacier and ice-free terrain.To take into account the effect of debris on ice melt, we adopted an index of thermal resistance, defined as debris thickness divided by debris conductivity, which could be obtained from thermal remote sensing data and reanalysis climate data.Using multiple ASTER data taken on different dates, we obtained distributions of thermal resistance and albedo for the Trambau Glacier in the Nepal Himalaya that were not calibrated by in situ observational data.Both thermal resistance and albedo had uncertainties of approximately 30 % (Fig. 6), and we calculated that these uncertainties could translate into a runoff uncertainty of approximately 8 % (Table 3).Our calculation, which was validated with observational runoff data in the 1990s (Fig. 8), showed that meltwater from both debris-covered and debris-free ice bodies contributed to more than half of the total present runoff (Table 2).In particular, the debris-covered ice supplied significant excess meltwater to the total runoff (Fig. 9).A sensitivity analysis, in which no debris cover was assumed, suggested that the excess meltwater was attributable mainly to the lower elevation location and less importantly to the acceleration effect of thinner debris cover (Table 3) because the elevation profiles of runoff from debris-covered and debris-free surfaces were comparable (Fig. 10a).
Sensitivity analysis showed that change in precipitation affects runoffs from the ice (both debris-covered and debrisfree) and the ice-free terrain in opposing directions.Increased precipitation suppresses ice melting through the high albedo of snow cover, whereas runoff from ice-free terrain increases with precipitation.The two effects compensate each other, so that the response of the total runoff is smaller than that for changes in air temperature (Table 6).However, the potential response to change in precipitation could be complicated for a large perturbation (Fig. 11a).In particular, a deficit of precipitation could alter the seasonal cycle of runoff (Fig. 11b).It is also noted that responses of glacier extent and/or debris distribution have to be taken into account for a longer timescale, though the static condition was assumed in this study.
Other studies calculating heat conduction through a debris layer have accurately reproduced the melt rate of ice beneath the debris mantle if its thickness and conductivity are known (Nicholson and Benn, 2006;Reid and Brock, 2010).Even for a single glacier, however, the distributions of debris thickness and thermal conductivity are unobtainable because of the heterogeneously rugged surface of debris-covered areas.In this regard, our approach using thermal resistance is a practical solution to calculate the ice melting under the debris cover on a large scale, such as a basin or region because the concept of thermal resistance involves coexistences of debris with various thickness, ice cliffs and supra-glacial ponds (Nakawo et al., 1993).The assumption of a linear temperature profile within the debris layer may cause large uncertainty in both deriving thermal resistance and calculating ice melt.In particular, this linear approximation is unrealistic when the debris layer is too thick.Further research is required to understand whether this method can be applied to thicker debris layers or whether any modifications are required.Apart from debris processes, settings for precipitation (ratio to reanalysis data and elevation gradient) will be the main source of uncertainty.In particular, precipitation would decrease with elevation in extremely higher and thus colder environments.Mass balance data from such high elevations enable us to gain more insight on hydrology in the Himalayan catchment.
The Supplement related to this article is available online at doi:10.5194/hess-18-2679-2014-supplement.

Figure 1 .
Figure 1.Tsho Rolpa Glacial Lake-Trambau Glacier basin.(a) Catchment (green line) and (b) categorized surface features with thermal resistance of debris cover.Inset box shows locations of Kathmandu (KTM), Mt.Everest (EV), and the study site.Average thermal resistance is superimposed on the debris-covered area where available.Yellow cross in (a) denotes the location at which meteorological and hydrological observations were conducted in the 1990s.The background image is ASTER data taken in February 2006.

Figure 3 .
Figure 3. Schematic diagram of the bucket model used in this study (modified afterMotoya and Kondo, 1999).Surface storage is used to calculate energy and water balance of the ice-free terrain (see the Sect.2.3.3).Internal and ground storages are used to calculate final daily runoffs for the individual components such as the debriscovered surface (R d ), the debris-free glacier (R g ), the ice-free terrain (R t ), and the lake (R l ).

Figure 4 .
Figure 4. Distributions of thermal resistance for individual ASTER scenes on the Trambau Glacier, which we used to generate the averaged thermal resistance for the calculations (Fig. 1b).

Figure 5 .Figure 6 .
Figure 5. Distributions of thermal resistance in individual ASTER scenes on the Trambau Glacier, which were not used in the runoff calculations.

Figure 7 .
Figure 7. Root mean square difference (D RMS , colour shading) and Nash-Sutcliffe model efficiency coefficient (E N , contour lines) of the model performance for the Tsho Rolpa Glacial Lake-Trambau Glacier basin calculated for the period 1993-1996 (location shown as yellow cross in Fig.1a), as a function of precipitation ratio (horizontal axis) against the original Aphrodite precipitation and elevation gradient of precipitation (vertical axis).We adopt 55 % as the precipitation ratio and 35 % km −1 as the elevation gradient of precipitation for subsequent analysis (thin dashed lines).The thick dashed line denotes the 74 % precipitation ratio isoline for the whole basin.

Figure 8 .
Figure8.Observed and calculated runoff at the outlet of the Tsho Rolpa Glacial Lake-Trambau Glacier basin for the period 1993-1996 (location shown as yellow cross in Fig.1a).

Figure 9 .
Figure 9. Seasonal cycles of (a) daily runoff and (b) daily runoff depth of the Tsho Rolpa Glacial Lake-Trambau Glacier basin calculated for the period 1979-2007 (28 years).Shading denotes interannual variability obtained for the same period.The interannual variability of runoff depth from the lake is not shown for better visibility in (b).Also shown is the seasonal cycle of precipitation averaged for the whole basin (thin red lines).

Figure 10 .
Figure 10.Elevation profiles of (a) annual runoff depth over debrisfree glacier (Glacier) and debris-covered (Debris) surfaces, and of (b) responses under conditions of warming (air temperature increase 0.1 • C) and wetting (precipitation increase 10 %) of the Trambau Glacier calculated for the period 1979-2007 (28 years).Shading and error bars in (a) denote interannual variability for the same period.

Figure 11 .
Figure 11.(a) Response of annual runoffs to changing precipitation ration against the control condition and (b) seasonal cycles of total runoff (thick lines) and precipitation (thin dotted lines) in the two extreme cases of the Tsho Rolpa Glacial Lake-Trambau Glacier basin calculated for the period 1979-2007 (28 years).Response in (a) is described by the anomaly with respect to the control calculation.

Figure S1 .
Figure S1.Meteorological variables observed at the outlet (the cross in Fig. 1a) of the Tsho Rolpa Glacial Lake-Trambau Glacier basin for the period 1993-1996.T a , P p , R S , h r and U denote air temperature, precipitation, solar radiation, relative humidity and wind speed, respectively.The smoothed line in the solar radiation plot denotes radiation calculated at the top of the atmosphere.

Table 1 .
Abbreviation, unit and value of parameters used in this study.

Table 2 .
Area, ratio of area, annual runoff, runoff contribution, and runoff depth of the Tsho Rolpa Glacial Lake-Trambau Glacier basin for different surface types.Errors represent interannual variability calculated for the period 1979-2007 (28 years).

Table 3 .
Uncertainty in runoff due to thermal resistance (R T ) and albedo (α), and topographical assumptions of the debris-covered area of the Trambau Glacier.Influences are obtained by averaging anomalies from positive and negative cases (Eq.33 in Sect.4.1).Combinations of changes in both parameters in different directions were also tested (+δR T and −δα, −δR T and +δα).Errors represent interannual variability calculated for the period 1979-2007 (28 years).
* Control variables for the no lake assumption are the summations of debris and lake.

Table 4 .
Uncertainty in runoff due to parameters used in the model of the Tsho Rolpa Glacial Lake-Trambau Glacier basin.Influences are obtained by changing each parameter by ±30 %, averaged anomalies from positive and negative cases (Eq.33), and then expressed as a percentage of the control values (Table1).Nash-Sutcliffe model efficiency coefficient (E N ) is obtained against the control calculation.

Table 5 .
Annual runoff and runoff depth associated with the presence of ice in the Tsho Rolpa Glacial Lake-Trambau Glacier basin.Errors represent interannual variability calculated for the period 1979-2007 (28 years).

Table 6 .
Sensitivities of annual runoff (million m 3 ) and runoff depth (mm) (parentheses) associated with changes in air temperature (dT a , 0.1 • C) and precipitation (dP p , 10 % or 103 mm) for the Tsho Rolpa Glacial Lake-Trambau Glacier basin.Also shown are sensitivities associated with interannual variability (δ, standard deviation) of air temperature and precipitation calculated for the period 1979-2007 (28 years).