Modified temperature index model for estimating the melt water discharge from debris-covered Lirung Glacier , Nepal

In the Nepalese Himalayas, the complex topography, occurrence of debris covered glaciers, and limited data availability creates substantial difficulties for modelling glacier melt. The proper recognition of melt processes governs the accurate estimation of melt water from glacier dominated systems, even in the presence of debris-covered glaciers. This paper presents a glacier melt model developed for the Lirung subbasin of Langtang valley, which has both a clean glacier area, 5.86 km2, and a debris-covered glacier area, 1.13 km2. We use a temperature index approach to estimate sub-daily melt water discharge for a two week period at the end of monsoon, and the melt factor is varied according to the aspect and distributed to each grid processed from the digital elevation model. The model uses easily available data and simple extrapolation techniques capable of generating melt with limited data. The result obtained from this method provides accurate estimate with an R2 value of 0.89, bias of 0.9% and Nash-Sutcliffe efficiency of 0.86, and suitable in Himalaya where data availability is major issue.


INTRODUCTION
Debris-covered glaciers are common features of the Himalaya and mountain regions throughout the world.The melt rates of debris-covered glaciers vary significantly compared to other types of glacier (e.g.Scherler et al., 2011;Bocchiola et al., 2011).Debris thickness is a controlling factor that influences melt rates of glacier ice; melt is enhanced with debris thicknesses of a few centimetres, but reduced at greater debris thicknesses (Nakawo and Rana, 1998;Scherler et al., 2011).Global temperature rise has resulted in accelerated glacier retreat in most regions (Mingjie et al., 2013).In light of the existing uncertainty in future climatic scenarios, data gaps and glacier mass loss, except in the Karakarom region (Hewitt, 2005), the Himalayan region presents a huge challenge in terms of future water availability and management (Scherler et al., 2011).Modelling is an appropriate tool for places having scarce data (Sivapalan et al., 2003), provided that melt processes can be approximated (Hock, 2003) even in the debris-covered glacier (Bocchiola et al., 2011).
In general, there are two approaches for snow and ice melt modelling.The energy balance approach explicitly models all components of the surface energy balance, and temperatureindexed-modelling considers temperature as the main variable controlling melt (Hock, 2003;Tobin et al., 2013).Numerous studies have attempted to model the melt water discharge using a positive degree day approach (Kayastha et al., 2000a(Kayastha et al., , 2005)).Several studies have incorporated shortwave radiation to improve sub-daily melt totals (Hock, 1999;Pellicciotti et al., 2005).Although the energy balance approach best describes the effects of debris cover on melt totals (Hock, 1999(Hock, , 2003)), input data availability is a significant constraint.Uncertainty in the hydrological models might arise due to errors in input data, inappropriate parameter selection and modelling approach used (Hughes et al., 2010).
In this study, we estimate melt from the debris-covered Lirung Glacier, and compare modelled and observed discharge for a small Himalayan catchment.The model uses a grid-based approach derived from a digital elevation model (DEM) and estimates melt precipitation and evapotranspiration for each grid cell.This paper also tests the effect on modelled discharge by changing the meteorological inputs and degree day assumption.

Study area
Lirung Glacier is a debris-covered glacier in monsoon dominated Langtang Valley, Nepal, located about 60 km north of Kathmandu Valley.Geographically, Lirung Glacier extends from 28.21 o N, 85.52 o E to 28.27 o N, 85.56 o E. The total catchment area is 13.8 km 2 with 1.13 km 2 of debris glacier and 5.86 km 2 of clean glacier area.Elevation within the catchment ranges from 3958 to 7232 m a.s.l.The average slope is 38 o with most of the area lying with an east and south facing aspect.Glacier boundaries in the catchment data are extracted from the International Centre for Integrated Mountain Development (ICIMOD) glacier inventory (Bajracharya et al., 2014).Further processing of glacier extent was carried out using the ASTER DEM 2011 of 30 m resolution.Slope within the catchment is derived from the available DEM and manual delineation is done for debris glacier area based on slope.For the debris-covered terminus, glacier boundaries are further adjusted following a terminus survey conducted with a Garmin 60CSx GPS on 3 December 2013 and considered that the glacier area had not changed from September to December.

Hydro-meteorological data
Water levels at the outlet were measured every 15 minutes using an Ott Pressure Level Sensor (PLS; Fig. 1) and averaged hourly from 21 September to 3 October, 2013.In order to calculate discharge, a rating curve was developed from 30 discharge measurement using a Gurley current meter.Stage is recorded at the outlet of a proglacial pond, which may act as a short-term reservoir delaying the transport time for meltwater to the outlet (Sakai, 1997).A DAVIS Vantage Pro Automatic Weather Station (AWS) was set up on Lirung Glacier at an elevation of 4132 m a.s.l.from 21 September to 3 October, 2013 and measured temperature and precipitation every 5 minutes (Fig. 1).The hourly average temperature and hourly sum of precipitation recorded by the AWS were used as model input data.A permanent AWS in Kyangjing (3862 m) and permanent AWS in Yala (5090 m) (Fig. 1) provided additional meteorological inputs for the experiments described below.To extrapolate the hourly meteorological data (Fig. 3) obtained from the AWS at Lirung Glacier, the 30-m resolution ASTER DEM was used.Collected data are distributed on each grid by using temperature, and precipitation gradients.In this study the precipitation gradient is used as described by (Immerzeel et al., 2014) and the temperature gradient is calculated from station data available from Yala and Kyangjing temperature data.

MODELLING APPROACH
To estimate melt-water discharge a simple water balance equation is considered: (1) where, Q is total discharge (m 3 /s), M(x, t) is the snow and ice melt of location x at time t (m 3 /s), P is the precipitation (mm), ET is the evapotranspiration (mm/h) and Qbase is the base flow (m 3 /s).The lowest flow within the stream with negligible contribution of melt water to discharge is considered as initial base flow and the concave method (Subramanya, 1994) is used to differentiate baseflow diurnally.
In our study, a simple temperature index approach is used to estimate snow and ice melt at an hourly time-step.The general melt equation is given by: where, T(x,t) is the air temperature (ᵒC) and DDF (mm d -1 o C -1 ) is the positive degree day factor used in the model.Three methods are used to estimate the melt water discharge for Lirung Glacier.
The first method specifies a constant degree day factor for snow, ice and debris-covered ice.The second method uses a degree day factor adjusted for the elevation effect (Pradhananga et al., 2014) while the third method uses aspect and elevation modified degree day factors.There is variation of degree day factor with respect to the location of glacier cover.In this model, there is a different degree day factor for debris-covered ice DDFd, clean ice DDFi and snow DDFs.Below 5000 m, DDFs is 5 mm d -1 o C -1 , and DDFi is set to 6 mm d -1 o C -1 .Above 5000 m, DDFs is set to 6.5 mm d -1 o C -1 and DDFi 7.5 mm d -1 o C -1 for ice ablation above 5000 m.These values are calibrated according to the range of degree day factor provided (Pradhananga et al., 2014).For debris covered glacier in the basin, a melt reduction factor of 0.5 is applied (Kayastha et al., 2000b).
In order to represent the variations in melt with respect to aspect and net solar radiation, the degree day factor is also modified following the approach of (Immerzeel et al., 2012): (3) where, DDFm is modified degree day factor with respect to aspect (mm d -1 o C -1 ), Rexp = 0.2 is the factor quantifying the aspect dependence.Due to daily afternoon cloud formation in the Himalaya region the east-facing slopes typically receive more solar radiation than west-facing slopes (Miehe, 1989).In this study, however, east and west facing slopes are considered to have similar degree day factors.
A threshold temperature governs whether precipitation occurs as snow or rain, and we use the transition scheme outline by Kayastha et al. (2000a) to separate snow and rain.
The recession coefficient is used in order to modify the recession flow for this watershed (Martinec et al., 1983): where, K n is the recession coefficient for the n-th hour, Qo is the initial discharge (m 3 /s), Qn is discharge after recession (m 3 /s).The value of K is obtained by solving the equation ( 6), where the constant x = 0.98 and y = 0.36 is obtained by a recession flow plot (Martinec and Rango, 1986): (5) Point evapotranspiration at the on-glacier AWS is estimated using the modified Hargreaves method (Sperna Weiland et al., 2012): where, Ra is the extraterrestrial radiation (MJ m -2 hour -1 ), T is the hourly temperature and TR is the diurnal temperature range ( ᵒ C).The extraterrestrial radiation is estimated based on (Allen et al., 1998).
The model is run several times shuffling the degree-day factor and changing input climate data.We also test the impact of the input climate data and extrapolation method on the modelled discharge by using two different assumptions: (1) a single average temperature gradient, and (2) hourly temperature gradients, calculated from available station data.

RESULTS AND DISCUSSION
The observed and modelled discharges simulated are shown in (Fig. 2).All three methods capture the diurnal variability in streamflow at the Lirung Glacier outlet.However, the aspect-elevation modified degree day factors provide better results compared to classical methods and the elevation adjusted method.For all three methods, the R 2 value is >0.80 but there is a clear difference in the Nash-Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) as well as mean bias.The result obtained using the modified temperature index method is efficient and could be further explored and modified compared to the classical and elevation threshold methods.The aspect and elevation modified temperature index method has the greatest accuracy of the three melt modelling methods, with R 2 = 0.89, a small mean bias (2.4%) and the highest NSE (0.86).For the elevation-modified degree day factors, R 2 is 0.84, mean bias is 3.7 % and NSE is 0.71, while the classical approach has R 2 = 0.85, mean bias = 5.1 % and NSE = 0.66.
Since a large part of Lirung Glacier has a south and east facing aspect, the aspect-elevation modified method provides the most accurate result.For other methods, such consideration is neglected which caused a greater mean bias from the observed discharge.Discharge estimated using the elevation modified degree day factor provides more efficient results than the classical method but still has a higher mean bias and lower NSE.

Uncertainty analysis
To examine uncertainty in the modelled discharge that arises from meteorological forcing, the model is run for Lirung Glacier using the meteorological dataset available from the Kyangjing station for the same period.All model parameters are based on the aspect variable degree day approach for the first model run.For the second model run, hourly lapse rate obtained from the Lirung-Kyangjing station and Yala-Kyangjing is used.Discharge estimated from the Kyangjing and Yala datasets demonstrates the importance of the input data assumptions.There is either underestimation or over-estimation of discharge when single or hourly lapse rate is used (Fig. 3(a)).There is significant underestimation of peak discharge using the Kyangjing data and Yala data with a single mean temperature gradient and discharge estimated using hourly lapse rate from Yala Kyangjing data.However, there is significant overestimation of discharge when temperatures are extrapolated using the hourly lapse rate computed from Lirung-Kyangjing dataset.The difference in melt water discharge obtained from different extrapolation methods is due to the location of the AWS and the difference in type of land-use.The Kyangjing station is located in a more exposed location in Langtang valley and has a different climate to the Lirung station which is installed on the glacier, while Yala station is located near the Yala Glacier which is a clean glacier.Temperature typically decreases with the elevation but there is a diurnal fluctuation in temperature gradient for the highly elevated regions (Fig. 3(b)).The hourly temperature gradient for Lirung-Kyangjing station is completely different to Yala-Kyangjing station (Fig. 3(c)) which causes the difference in estimated discharge (Fig. 3(a)).In addition to temperature extrapolations, there are other parameters resulting in uncertainty in model simulations.The debris thickness, rating curve, glacier mapping, pond near the outlet, and evapotranspiration estimate can cause additional uncertainty.However, our results demonstrate that in situ meteorological observations are invaluable for providing accurate estimates of glacier melt.

CONCLUSION
This study presents a melt modelling approach that can be used in glacierized river basins with limited data and is useful for the Himalayan region.A simple temperature index model modified based on elevation and aspect to account for differences in solar energy available for glacier melt, provides the most accurate estimate of discharge.We further demonstrate that in situ data are invaluable for estimating glacier melt, as extrapolation from nearby stations using constant temperature gradients results in significant errors.
Debris-covered glaciers present a unique challenge to glacier melt modelling.However, our study suggests that cumulative melt in catchments with both clean and debris-covered ice can be successfully modelled using relatively simple approaches.Future studies on the energy balance of debris-covered glaciers should help improve parameterizations examined in this study.

Fig. 1
Fig. 1 Location map of the study area showing location of AWS and PLS.

Fig. 2
Fig. 2 Time series of simulated discharge vs observed discharge: (a) discharge obtained by varying degree day with aspect and elevation threshold, (b) discharge obtained using elevation threshold for changing degree day factor, and (c) classic degree day method.

Fig. 3
Fig. 3 Uncertainty analysis: (a) estimated discharge using single and hourly temperature gradient from different stations; (b) air temperature at Lirung, Kyangjing and Yala stations; and (c) hourly temperature gradient.