Quantification of different flow components in a high-altitude glacierized catchment ( Dudh Koshi , Nepalese Himalaya )

In a context of climate change and water demand growth, understanding the origin of water flows in the Himalayas is a key issue for assessing the current and future water resources availability and planning the future uses of water in downstream regions. This study estimates the relative contributions of rainfall, glacier and snow melt to the Khumbu River streamflow (Upper Dudh Koshi, Nepal, 146 km2, 43% glacierized, elevation range from 4260 to 8848 m a.s.l.), as well as their seasonal variability during the period 2012-2015, by using the physically based glacio-hydrological model DHSVM-GDM (Distributed 5 Hydrological Soil Vegetation Model Glaciers Dynamics Model). One of the main issues in high elevated and glacierized catchments hydrology is the limited representation of cryospheric processes, which control the evolution of ice and snow, in distributed hydrological models. Here, the impact of different snow and glacier parametrizations was tested by modifying the original DHSVM-GDM snow albedo parametrization, by adding an avalanche module, and by adding a reduction factor for the melt of debris covered glaciers. Results show that this new version of DHSVM improves the simulation of the snow covered 10 area and the glacier mass balances, thus improving the reliability of the overall hydrological simulation. In the presented case study, ice and snow melt contribute each more than 40% to the annual outflow. 69% of the outflow originates from glacierized areas. Our simulations also highlight that winter flows are mainly controlled by the release from the englacial water storage. In general, it is shown that the choice of a given parametrization for the snow and glacier processes has a significant impact on the simulated water balance. The sensitivity of the model to the glaciers inventory was tested, demonstrating that the uncertainty 15 related to the glacierized surface leads to an uncertainty of 20% on the simulated ice melt component.


Introduction
The Himalayan mountain range is known for being the water tower of Central and South Asia (Immerzeel et al., 2010).Its high elevated glaciers and snow cover play an important role in the regional hydrological system (Kaser et al., 2010;Racoviteanu et al., 2013) and provide water resources for the population living in the surrounding countries (Viviroli et al., 2007;Singh et al., 2016;Pritchard, 2017).
In the Hindu Kush-Himalaya (HKH) region climate change is expected to cause shrinkage of the snow and ice cover (Bolch et al., 2012;Benn et al., 2012;Kraaijenbrink et al., 2017).Changes in glacier and snow cover runoff are likely to have a  3 Data and model setup

Database
To describe the topography of the study area, an ASTER DEM originally at 30 m resolution was resampled to a 100 m resolution.The SOTER Nepal soil classification (Dijkshoorn and Huting, 2009) and a landcover classification from ICIMOD (Bajracharya, 2014) were used for the soil and landcover description.Meteorological data were available at hourly time steps from three automatic weather stations (AWS) located at Pangboche (3950 m a.s.l.), Pheriche (4260 m a.s.l.) and Pyramid (5035 m a.s.l.) (Table 1).Since December 2012, the precipitation has been recorded at the Pheriche and Pyramid AWS by two Geonor T-200 sensors designed to measure both liquid and solid precipitations.Data were corrected for potential undercatch following the method used by Lejeune et al. (2007) and Sherpa et al. (2017).Precipitation at the Pangboche station was recorded with a tipping bucket.Air temperature, wind speed, relative humidity short-wave radiation and long-wave radiation measurements at Pheriche and Pyramid were provided by the EvK2-CNR stations.
Discharge measurements of the Khumbu River at Pheriche were obtained using a pressure water level sensor at 1 h time step since October 2010.The MODImLab algorithm developed by Sirguey et al. (2009) was applied to MODIS reflectances data to obtain daily albedo and snow fraction satellite images for the period 2010-2015.We used the Sirguey et al. (2009) algorithm rather than the MOD10A1 500 m snow products because it generates daily regional snow cover images at 250 m resolution and applies corrections on atmospheric and topographic effects which makes the snow cover maps more realistic on mountainous areas.27 cloud free Landsat8 images were used to generate snow maps at 30 m resolution between 1 November 2014 and 31 December 2015.A NDSI (Normalized-Difference Snow Index) threshold of 0.15 was taken to separate snow free and snow covered pixels on Landsat8 data as proposed by Zhu and Woodcock (2012).Daily snow cover maps were then retrieved from the MODImLab snow fraction product: areas with a snow fraction above 0.15 were defined as snow covered areas so that the MODImLab Snow cover area (SCA) matches the Landsat8 SCA on the common dates.For the rest of this study we call MODIS data albedo and snow cover data obtained with the MODImLab algorithm.We also used snow albedo data from in-situ measurements at Pyramid and Changri Nup (Table 1).
For describing the glacierized area in the basin we compared three different glacier outlines available as vector layers for the Khumbu region: the glacier delineation proposed by Racoviteanu et al. (2013) specifically set up for the Dudh Koshi basin; the GAMDAM inventory covering the entire Himalayan range (Nuimura et al., 2015); and the ICIMOD inventory (Bajracharya et al., 2010) (Fig. 3).The three outlines have been derived on different grids, from different datasets at different spatial resolutions and covering different temporal periods (see Table 2), thus leading at different results.was used for simulating outflows at Pheriche.DHSVM is a physically based, spatially distributed model which was developed for mountain basins with rain and snow hydrological regimes (Wigmosta et al., 1994;Nijssen et al., 1997;Wigmosta and Burges, 1997).A glacier dynamics module was recently implemented in DHSVM by Naz et al. (2014)  and mass balance module for glaciers (Andreadis et al., 2009;Naz et al., 2014) and has been applied in a number of studies for snow and cold regions hydrology (e.g., Leung et al., 1996;Leung and Wigmosta, 1999;Westrick et al., 2002;Whitaker et al., 2003;Zhao et al., 2009;Bewley et al., 2010;Cristea et al., 2014;Frans et al., 2015).Distributed meteorological data (air temperature, precipitation, relative humidity, wind speed, and shortwave and longwave incoming radiation) are requested as input, as well as distributed geographical information (elevation, soil type, landcover, soil depth, and ice thickness).

Snow albedo parametrization
In the original DHSVM-GDM version, the snow albedo α s [-] is set to its maximum value α smax (to be fixed either by calibration or from observed albedo values), after a snowfall event and then decreases with time according to the following equations (Wigmosta et al., 1994): Where N is the number of days since the last snowfall, λ a [-] and λ m [-] correspond to 0.92 and 0.70 for the accumulation season and the melt season, respectively, and T s is the snow surface temperature [ MODIS albedo images and the albedo measurements from Pyramid and Changri Nup were used to analyse the decrease of snow albedo with age in various locations of our study area.Figure 4 compares the observed albedo decay as a function of time for snow events with at least three consecutive days without clouds after the snowfall with the albedo parametrization in DHSVM-GDM.Since the observed values are not well represented by the standard albedo decrease, the parametrization was replaced by Eq. 2, with a decay of the albedo when there is no new snowfall inspired by the ISBA model albedo parametrization (Douville et al., 1995) and with the fresh snow albedo modified as a function of the amount of snowfall: Where α st−1 is the albedo from the previous time step, α smin is the minimal snow albedo of 0.3 (estimated using the mean minimal albedo values observed at the station and on MODIS images), N is the number of days since the last snowfall, c is the coefficient of the exponential decrease [days Where Z is the elevation of the cell in m a.s.l.
The new function for the decrease of the snow albedo is also shown in Fig. 4.

Avalanches parametrization
The transport of snow by avalanches is not represented in the original version of DHSVM-GDM.The absence of avalanches in the model can lead to an unrealistic accumulation of snow in steep high elevated cells, where the air temperature remains below 0 • C, and to a deficit of snow in the lower areas, where snow melt occurs during the melting season.The simulated water balance directly depends on the snow cover, thus not considering avalanches can lead to significant errors.In order to address these discrepancies, an avalanche module inspired by Wortmann et al. (2016) was implemented in DHSVM-GDM.
The avalanche model transfers snow to downslope neighbour cells under the following conditions: if the terrain slope is steeper than 35 • and the amount of dry snow water equivalent (total snow water equivalent minus liquid water content) is higher than 30 cm: 5 cm of snow water equivalent remains in the cell and the rest is removed by avalanches; if the terrain slope is less steep than 35 • but the difference in snow water equivalent compared to the downslope neighbour cells is larger than 5 cm: 95 % of the difference is removed by avalanches.
The transfer of snow by avalanches is based on the surface runoff routing in DHSVM-GDM: at every time step starting from the highest cell of the DEM to the lowest, each cell can transfer snow to its closest downslope neighbour cells (between 1 and 4 cells).Within the same time step, the amount of snow in the receiving cells is actualised and the avalanches propagate downslope until the conditions cited above are no longer respected.

Glacier parametrization
The distributed ice thickness is derived from the terrain slope following the method described in Haeberli and Hölzle (1995).
Since the standard DHSVM-GDM model does not take into account the impact of the debris layer on melting of the glaciers, the insulating effect of the debris layer is not represented.Here, we implemented a reduction factor for ice melt generated in grid cells with debris-covered glaciers (see Sect. 3.3).
Moreover, in the original DHSVM-GDM version glacier, melt water is instantaneously transferred to the surface of the soil layer which is parameterized as bedrock under glaciers (Naz et al., 2014).This significantly underestimates the transfer time through glaciers.In this study we modified the soil parametrization in order to include the liquid water storage by increasing the soil depth to 2 m under all glaciers, in order to include the liquid water storage and the delayed englacial transfer.This also implies to change the values of three soil parameters under glaciers: the vertical and the lateral conductivities, as well as the porosity (Table A2).These were fixed by optimizing the recession shape of the hydrographs.

Quantification of the flow components
Quantifying the relative contributions of ice melt, snow melt and rainfall in the river discharge at different time scales is a difficult task because hydrological models usually do not track the origin of water particles during transfer within the catchment Weiler et al. (2018).
In this study, two different definitions were used to estimate the hydrological contributions.First, we estimate the contributions of ice melt (V icemelt ), snow melt (V snowmelt ), and net rainfall (V rainN et ) to the total production of runoff (V runof f ) (definition 1) according to the following equations: Where V glAcc is the amount of snow that is transferred to the ice layer by compaction on glaciers (Naz et al., 2014), S ice and S snow are the amounts of sublimation from the ice and snow layers, dIwq dt and dSwq dt are the variations of the ice and snow storages, P solid and P liquid are the amounts of solid and liquid precipitation, and E int is the amount of evaporation from intercepted water stored in the canopy.It is worth noting that at daily or monthly time scale, the sum of these contributions may not be equal to the outflow at the catchment outlet as liquid water can be stored by or evacuated from the soil or glaciers.
In order to evaluate the seasonal components of the outflow at the catchment's outlet, we also define the hydrological contributions as fractions of the outflow coming from the different contributive areas (definition 2): direct glacier contribution: direct runoff from glacierized areas, delayed glacier contribution: resurging melt water stored inside glaciers, direct snow contribution: direct outflow from snow covered non-glacierized areas, direct runoff: direct runoff from areas without snow and glaciers, subsurface and groundwater contribution: resurging water from the soil in non-glacierized areas resulting from infiltrated rainfall, snow melt, as well as upstream lateral subsurface flows.
Figure 5 illustrates the two definitions of the different contributions to outflows.Definition 1 allows assessing the annual impact of glaciers and snow cover on the water production, while Definition 2 describes the intra-annual routing of the water within the catchment.Moreover, using the two definitions allows to directly compare our results with other hydrological modelling studies in the Dudh Koshi basin, which have estimated glaciers contributions either from effective ice melt (Savéan et al., 2015;Ragettli et al., 2015;Soncini et al., 2016) or runoff from glacierized areas (Immerzeel et al., 2012;Nepal et al., 2014).Further, we assessed the impact of the definition of hydrological components on the estimated glaciers contribution.
Flow components were estimated for the period 2012-2015 at annual scale, on the basis of the glaciological year (from 1 December to 30 November), as well as monthly, daily, and sub-daily scales, in order to have a better understanding of the seasonal variation of the hydrological contributions.

Experimental set-up
Simulations were run with a 1 h time step and a spatial resolution of 100 m for the period from 1 November 2012 to 27 November 2015 corresponding to the period with most available meteorological and discharge data.A soil depth map was derived from the DEM using the method proposed in the DHSVM-GDM documentation (Wigmosta et al., 1994).As a results, soil depth outside glacierized areas ranges between 0.5 and 1 m.Under the glaciers, the soil depth was set to 2 m.All parameter values retained for the simulations (with no calibration) are summarized in Appendix A.
In order to test the impact of the representation of the cryospheric processes on the hydrological modelling, we performed simulations with the four following configurations: -v0: original DHSVM-GDM snow and glacier parametrization; Using configuration v3, we also tested the impact of using different glaciers outlines (the GAMDAM and ICIMOD inventories were also considered for simulations).The debris-covered glacier melt reduction factor estimated in Konz et al. (2007), Nepal et al. (2014) and Shea et al. (2015) are respectively equal to 0.3, 0.33 and 0.47.Thus, values between 0.3 and 0.5 were also considered (in addition to the reference of 0.4) in order to evaluate the sensitivity of the model to the debris covered glacier reduction factor.

Model forcing
Meteorological data from the Pheriche and Pyramid stations (Table 1) were spatialized over the basin by an inverse distance interpolation method.Altitudinal lapse rates of precipitation and temperature were calculated at 1 h time step from data collected at Pangboche (3950 m a.s.l.), Pheriche (4260 m a.s.l.) and Pyramid (5035 m a.s.l.) (Fig. 6).Only significant lapse rates with R 2 values higher than 0.75 were retained for precipitation (43% of the dataset).For smaller R 2 , the lapse rate is considered as not significant and thus set to 0.
In this study, the precipitation lapse rates show a large seasonal variability with daily lapse rates ranging from -41 to 9 mm km -1 .Precipitation decreases with elevation during the monsoon season and increases with elevation in winter: during the simulation period, we found 450 days (40%) with no precipitation, 83 days (8%) with a strictly negative lapse rate and 165 days (15%) with a strictly positive lapse rate.Concerning temperatures, daily lapse rates range from -0.009 to +0.006 • C m -1 .
We found only 10 days (1%) showing a temperature inversion with a positive daily lapse rate.

Model evaluation
A multi criteria evaluation was made considering simulated outflows, SCA and glacier mass balances.Discharge measurements at Pheriche station were used as reference for the evaluation of simulated outflows.A 15% confidence interval was retrieved as representative of the uncertainty of measured discharge.Nash-Sutcliffe Efficiency (NSE) (Nash and Sutcliffe, 1970) and Kling-Gupta Efficiency (KGE) (Gupta et al., 2009) were chosen as objective functions and applied to daily discharges.The simulated SCA was evaluated in comparison to daily SCA derived from MODIS images.Because a large number of MODIS images suffer from cloud coverage, we only compared the simulated and observed SCA during days with less than 5% of cloud cover on the catchment.The simulated glacier mass balances were evaluated at basin scale by a comparison with published regional geodetic mass balances and at local scale using available stake measurements on the clean ice West Changri Nup and the Pokalde glaciers (Sherpa et al., 2017).shows that after full coverage it does not decrease fast enough compared to MODIS data.Figure 8 demonstrates that the snow cover duration is over-estimated for the entire catchment area.This indicates that in the simulations snow does not melt fast enough with the original parametrization.Configuration v1 with the modified snow albedo parametrization (Eq.2) accelerates the snow melt and improves the SCA simulation (Fig. 7).The RMSE between the simulated and observed SCA decreases from 29% using v0 to 14% using v1 and v2. Figure 8 shows that with configuration v1 in some areas located at high elevation the 10 snow cover duration is underestimated.This biais is rectified in configuration v2 since the avalanche module transfers snow from high elevated and sloping cells downward and corrects the lack of snow observed with configuration v1 at the edges of the permanent snow cover (Fig. 8).The results for the SCA and snow cover duration using the configuration v3 show no difference compared to the configuration v2 since only the ice melt rate for debris covered glaciers is modified.
Our results show that the snow parametrization has a significant impact on the simulated glacier mass balance.The mass balance simulated with v0 is on average -0.82 m w.e.yr -1 and decreases to -2.02 m w.e.yr -1 with the corrected snow albedo (v1) since the modified snow albedo parametrization accelerates the snow melt which leads to more uncovered ice and stronger glacier melt.The avalanche module (v2) adds snow on glaciers and increases the accumulation and, thus, reduces the glacier melt to -1.69 m w.e.yr -1 .Nevertheless the mass balance remains too negative compared to geodetic mass balances which suggests that the model produces too much ice melt.The implementation of debris-covered glaciers (v3) gives a mean annual glacier mass balance of -0.84 ±0.14 m w.e.yr -1 which is within the intervals of uncertainty and thus in good agreement with geodetic methods.
It is worth notify that only the order of magnitude of the simulated and geodetic mass balances can be compared because the considered areas are not exactly the same in the different studies and the considered time periods differ.Indeed, four of the geodetic mass balances are estimated for the Khumbu-Changri glacier, while mass balances from this study and from Brun  We also evaluated the mass balance at the point scale.Figure 10 shows the simulated mass balances with parametrizations v0, v1, v2 and v3 versus the observed mass balances of the two debris-free glaciers West Changri Nup and Pokalde measured insitu.Here, the configuration v3 gives the same results than the configuration v2 because the configuration v3 only modifies the ice melt rates on debris-covered glaciers.As previously, the simulated mass balances vary according to the model configuration and the variability between the different configurations is larger than the spatial and inter-annual variability.
With configuration v0, the model overestimates the point mass balances because of low snow melt rates (see also section 4.1.1).With configuration v1, the model simulates too much ice melt on the West Changri Nup glacier due to a lack of accumulation in the western part of the catchment and a too strong accumulation on the Pokalde glacier (Fig. 8).The configuration v2 improves the simulated mass balance by transfering snow due to avalanches on the West Changri Nup glacier and by removing exceeding snow accumulation on the Pokalde glacier.
These results show that the snow parametrization has a strong impact on the simulated glacier mass balance and that the new snow albedo parametrization and the avalanching module clearly improve the simulated glacier mass balance on debris-free glaciers.Nevertheless, regarding point mass balances, the agreement is far from being perfect, due either to simulation errors

Annual outflow and flow components
Figure 11 represents the annual outflow and flow components simulated with the different model configurations, indicating the impact of each modification of the snow and glacier parametrization on the simulated annual outflow and the glacier contribution to the runoff.Configuration v1 leads to a drastically increased outflow and ice melt component.This increase of ice melt is due to the acceleration of the snow melt compared to v0.Implementing the avalanche module (v2) reduces the ice melt component since more snow is maintained on the glaciers and also increases the snow melt component because snow that was originally accumulated at high elevation is transported downwards by avalanches.Configuration v3 including debris covered glaciers further reduces the ice melt, resulting in a simulated annual outflow close to the observations.This shows that the modification of one specific hydrological process (here, the representation of avalanches) can have a significant impact on the simulated hydrological response of the catchment and requires improving other processes (here, considering specific representation of debris-covered glaciers).
The configuration with all three modifications (v3) gives results similar to the original parametrization of DHSVM-GDM (v0) in terms of glaciers mass balance, improving slightly the annual outflow.The ice melt factor for debris covered glaciers and the avalanches compensate the increase of ice melt caused by the new snow albedo parametrization, but the modifications implemented in v3 impact the results for the flow components: on average, less ice melt and more snow melt are generated.This is particularly true for the glaciological year 2014-2015, when the ice melt contribution decreases from 60% with v0 to 16 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-34Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 9 February 2018 c Author(s) 2018.CC BY 4.0 License.41% with v3 and the snow melt contribution increases from 29% to 47%.This can be explained by the fact that 2015 is a particularly snowy year (Fujita et al., 2017), which highlights the importance of a correct representation of snow processes in the model.This also shows the need to use as much validation data as possible to assess the coherence between the ice, snow and hydrological processes and reduce the uncertainty on the flow components estimation.Moreover, the configuration v3 modifies the seasonal variation of the outflow by increasing winter discharges and reducing monsoon discharges (not illustrated here) which improves the monthly KGE (0.9 with configuration v3 versus 0.83 with configuration v0).
Configuration v2, which does not consider the debris-covered glaciers, overestimates the outflow at Pheriche with a mean bias of +32% compared to the annual observed outflow.Without the debris layer, the ice melt component represents 61% of the annual outflow, which is nearly twice the amount of ice melt obtained with v3 that includes debris-covered glacier melt.When the coefficient for debris covered glaciers melt varies between 0.3 and 0.5 (Table 3), it modifies the simulated annual outflow by ± 7% and the ice melt flow component ranges from 42 to 50%.
Using configuration v3, the simulated annual outflow always stays within the interval of confidence of the observed outflow.The losses by evaporation and sublimation are rather constant trough the simulation period ranging between 140 and 150 mm/yr.Results show an inter-annual variability of the hydrological contributions to the outflow.During the period 2013-2015, the ice melt component ranges from 41 to 50%, the snow melt component from 37 to 47% and the net rainfall component from 12 to 16%.These variations are related to the meteorological annual variability.The amount of rainfall decreases from 2013 to 2015 and explains the decrease of the net rainfall components from 155 mm in 2013 to 88 mm in 2015.The snow melt component is higher in 2013 because of warmer pre-monsoon and monsoon seasons.The ice melt component is mainly controlled by the amount of winter snowfall.In 2014 a low amount of snowfall was observed, so the snowpack melts more rapidly and the glaciers start melting earlier.In contrast, 2015 is a year with a lot of winter snowfall, which delays the beginning of the glaciers melt and explains the lower ice melt component.
On average, we find that 69% of the annual outflow comes from glacierized areas and that 46% of the annual outflow is due to ice melt and 41% to snow melt (Fig. 11).Soncini et al. (2016) studied flow components in the Pheriche catchment on the period 2013-2014 and estimated an annual ice melt component of 55% and a snow melt component for 20%.The ice melt components are thus quite similar in terms of relative contributions to outflow, which is not the case for the snow melt components.We think that the main reason of such a difference is that we use different precipitation input.Indeed, precipitation data are measured here by Geonor sensors, while in Soncini et al. (2016) precipitation data come from tipping buckets.At the Pyramid station, where both sensors are installed, the Geonor sensor measures 60% more precipitation than the tipping bucket over the period 2013-2015 and the main differences are in terms of solid precipitation (309 mm of mean annual snowfall measured by the Geonor sensor versus 83 mm measured by the tipping bucket, which is known to badly perform with solid precipitation).
The choice of the hydrological components definition leads to different perceptions of the glacier contribution to the outflow (Fig. 11a and 11b).When considering configuration v3, the glacier contribution to the outflow is 69% if we consider the contribution from the entire glacierized area (i.e.contributions of ice melt, snow melt and net rainfall) and 46% if we consider only the contribution from ice melt.Likewise, the snow components with definition 1 and 2 are equal to 41% and 6% respectively,  3) which clearly shows the importance of infiltration and subsurface flows in the water balance.The comparison between the two definitions of the hydrological shows that contributions must be explicitly specified in order to allow inter-comparison between models, especially for catchment with a large glacierized area.

Sensitivity to the glacier outline
The three inventories used in this study result in very different estimates of the glacierized area: between 43% and 26% of the Pheriche basin with the inventories proposed by Racoviteanu et al. (2013) and GAMDAM (Fig. 3).Table 4 presents the mean annual results of glacier mass balance, outflow, and flow components for the configuration v3 using the three inventories.
The GAMDAM inventory leads to a more negative glacier mass balance than the two others inventories with -1.17 m w.e.yr -1 against -0.84 and -0.89 m w.e.yr -1 for the Racoviteanu et al. (2013) and ICIMOD inventories.This is due to smaller glacier accumulation areas in the GAMDAM inventory.The amount of snowfalls collected over those areas is lower, leading to more negative mass balances: glaciers receive less snowfall for accumulation, which lowers the mass balance value.Concerning the simulated outflow and flow components, the GAMDAM and ICIMOD inventories lead to fewer ice melt than the Racoviteanu et al. ( 2013) inventory due to their smaller surfaces in ablation areas, which leads to a smaller simulated annual outflow.
From these results we estimate an uncertainty of 20% (407 mm with the Racoviteanu et al. (2013) inventory versus 327 mm with the ICIMOD inventory, cf.Table 4) on the simulated annual ice melt volume related to the glaciers outline.The glacier outline mainly affects the simulated outflow during the monsoon season, when the ice melt contribution to the outflow is more important and leads to an uncertainty of 8% (154 mm with the Racoviteanu et al. (2013) inventory versus 141 mm with the ICIMOD inventory) on the monthly discharges during monsoon season.This result shows that the choice of the glacier inventory as an input data of the glacio-hydrological model affects the simulations and is crucial for the results.Here, the Racoviteanu et al. (2013) inventory gives the best results in terms of glacier mass balance and the smallest bias with respect to the annual outflow.As its area is significantly larger than the other inventories, it gives the largest amount of ice melt.This potentially compensates a lack of precipitation due to a poor knowledge of the precipitation distribution over the catchment, specifically in the areas above 5000 m a.s.l. which constitute more than three quarters of the total area and for which no observations exist.It is worth noting that the glacier mass balances obtained with the Racoviteanu et al. (2013) and ICIMOD inventories are very similar, which shows that glacier mass balances alone are not sufficient for model evaluation: a consistent mass balance can lead to errors on the simulated glacier contributions and total outflow.It also shows that in a region with a large uncertainty in the precipitation forcing data (Savéan et al., 2015;Eeckman et al., 2017) and on the delimination of glacierized areas, calibrating glacier parameters on the mass balance can lead to inaccurate estimations of the total ice melt volume.

Simulated seasonal and diurnal flows and contributions
In this section we analyse the hydrological contributions to the outflow simulated with configuration v3, which gives the best results for snow cover area, glacier mass balance, and discharges (as previously pointed out).Figure 12  winter seasons.This is mainly due to glacier melt water stored inside the glaciers during pre-monsoon and monsoon seasons and continuing surging during winter, as well as to changes in the soil water storage (Fig. 12b and 13b).
Figure 13 show the mean monthly flow components averaged over the simulation period.From February to May/June, the runoff production is mainly controlled by snow melt (between 50 and 60%) and by ice melt (between 40 and 48%) (Fig. 13a).
The rainfall, snow melt and ice melt absolute contributions are all at their maxima in July and August during the monsoon season.During these two months, runoff is generated at 24% by net rainfall, 37% by snow melt, and 38% by ice melt.From October to January, the runoff is mainly produced by ice melt (up to 80% in December) and snow melt.This is consistent with the results from Soncini et al. (2016), who found a main contribution of snow melt during the pre-monsoon season, mixed contributions of rainfall, snow melt and ice melt during the monsoon season and mixed contributions of snow melt and ice melt during post-monsoon and winter season.Figure 13b shows that liquid water stored inside glaciers contributes to more than 30% of the outflow during the monsoon season and can contribute up to 78% during winter.This corresponds well to the studies of Ragettli et al. (2015) and Racoviteanu et al. (2013) concerning the Upper Langtang and the Dudh Koshi basin respectively, who found that most of the winter outflow surges from soil, channel, surface, and englacial storage changes.The soil contribution includes snow melt and rainfall that infiltrate in the upper soil layers and contribute rapidly to the outflow.
It represents a significant share of the monthly outflow with a maximum between March and May (more than 30%) and a minimum in post-monsoon and winter (less than 20%).The response of the soil storage is faster than the englacial storage as the main part of the soil infiltrated water resurges within a day, whereas liquid water can be stored for several months within the glaciers.Direct contributions from glacierized areas, snow areas, and direct runoff are highest during the monsoon season, when the englacial and soil storage is saturated.the observed outflow is rather constant during the day, with a weak peak around noon when the temperature is at its maximum.Almost all of the precipitation is in form of snowfall, leading to no direct response for the outflow.The peak around noon can be explained by snow melt or the melting of small frozen streams.During the monsoon season, there is a strong diurnal cycle of the precipitation with a maximum occurring during late afternoon or at night and a discharge peak occurring around midnight.
The model simulates ice and snow melt during day time with a maximum at noon as it was expected.Except for the monsoon season, it seems to simulate accurately the baseflow during night when there is no melt production; the discharge is then controlled by the release of the glaciers and soil storage.Unfortunately the model simulates a peak of discharge coming mostly from glacierized areas around 14 h, two hours after the peak of ice and snow melt, which does not correspond to the

Conclusions and perspectives
In this study we used a distributed physically-based glacio-hydrological model (DHSVM-GDM) to simulate the outflow of a small catchment in the Everest region and to estimate the different contributions to flows, which can be useful for water resources and water-related risks management.
One of the main difficulties for hydrological modelling of highly glacierized catchments is to correctly simulate at the same time the hydrograph, the dynamics of the snow cover, and the glacier mass balances.In this study, an improvement of the parameterization of cryospheric processes in DHSVM-GDM was proposed in order to better represent ice melt under debris covered glaciers, avalanches, the storage and transport of melt water within glaciers.Snow cover dynamics has been shown to play an important role for the simulated water balance.The albedo parametrization proposed in this study and the implemented avalanche module enabled to simulate the snow cover spatial distribution and the glacier mass balances more accurately than in the original version of DHSVM-GDM by increasing the glacier accumulation and reducing ice melt.However, only after the consideration of the insulating effect of the debris layer on glaciers the overestimation of more than 30% of the annual outflow is rectified.
Among the different snow and glacier parametrizations that were tested, the most satisfactory (configuration v3) shows major contributions from glaciers and snow to the outflow with 46% of the annual outflow produced by ice and 41% by snow melt.
Winter flows are mainly controlled by the release of englacial water storage (up to 78% in December), which corroborates other studies (Racoviteanu et al., 2013;Ragettli et al., 2015).We estimate an uncertainty related to the ice melt reduction factor by debris (ranging from 0.3 to 0.5) of ±0.14 m w.e.yr -1 for the annual glacier mass balance and ±7% for the annual outflow.The glacier inventories used to outline the glacierized areas have also an important impact on the simulation results.
The three inventories tested in this study give estimations of the glacierized area ranging from 26 to 43% of the basin area and a corresponding uncertainty of 20% for the ice melt production.
The runoff coefficient (ratio between the annual outflow and annual precipitation) is on average equal to 1.4, which means that a considerable amount of water is withdrawn every year from the catchment through ice melt (eventually in the form of a delayed groundwater flow).Thus, if the precipitation regime (in terms of both intensity and phase) does not change within the next decades, the access to water resources is likely to be reduced, especially during the fall and the winter seasons, as the glaciers outflow will decrease due to glaciers shrinkage, even without taking into account climate warming.By the way, this study also reminds that glacial and snow contributions to the streamflow must be clearly defined, as considering Previous results are already useful for decision makers and stakeholders, but also indicate possible forthcoming works for increasing the simulations reliability and reducing uncertainties, especially at short time steps.Indeed, at daily and longer scales, the different hydrological components seem to be well reproduced by the model.However, an analysis of the diurnal cycle (Figure 14) showed that DHSVM-GDM responds too rapidly to the ice melt production and that the representation of the water storage within the glaciers needs to be improved.Further improvements should be based on studies that analyze the mechanisms of glaciers drainage systems in the Khumbu region and their influence on glaciers outflow (e.g., Gulley et al., 2009;Benn et al., 2017).These studies show that englacial conduits and supraglacial channels, ponds and lakes play a key role in the response of glaciers: DHSVM-GDM could thus be upgraded by implementing a parameterization of such systems and delay the response of glacierized areas, as successfully proposed, for instance, in the model developed by Flowers and Clarke (2002).Other processes such as supraglacial ponds and ice cliffs melting, transport of snow by wind or variation of temperature in the ice pack are not considered in DHSVM-GDM and their impact on the hydrological modelling should also be studied.
Another limitation of our model lies in the application of a uniform reduction factor for ice melt under debris covered glaciers.In order to have a more realistic representation of the debris, the reduction factor could be spatially distributed, at least following elevation or slope exposition, and eventually taken as time-variant.As an example, Ragettli et al. (2015) considered a distributed debris thickness in their glacio-hydrological model and obtained a mean reduction of ice melt under debris of 84%.
A main issue is related to the availability of data for validating the glacio-hydrological modelling parametrizations and outputs.The lack of in-situ measurements at high elevations and the uncertainty related to the available data prevent from assessing the performance of the model in simulating the different cryospheric processes we consider.They only allow to evaluate integrated variables such as the annual glacier mass balance, the seasonal snow cover area dynamics and the outflow at the catchment outlet with significant uncertainties that impact the estimation of the different flow components.It is worth noting that the snow cover distribution evaluation is particularly challenging on the Pheriche catchment, as clouds cover more than 50% of the catchment during more than 150 days per year on average (and almost all the time during the monsoon season).
Finally, a major source of uncertainty lies in the lack of meteorological data at high elevation (because of the inaccessibility on the terrain) and in the measurement errors when observations are available, due to extreme meteorological conditions.Indeed, precipitation is known to be underestimated due to the difficulty of measuring solid precipitation with rain gauges (Wagnon et al., 2013).Precipitation fields provided by different atmospheric models and satellites show also high discrepancies in this region of the Himalayas (Andermann et al., 2011;Palazzi et al., 2013;Ceglar et al., 2017) : another perspective of this study is to test the sensitivity of the model to different precipitation forcing data sets (in-situ, reanalysis, and satellite) and analyze the impact of different precipitation amounts and spatial distributions on the simulated discharges and flow components.

Figure 1 .
Figure 1.Study area : (a) Location map of Pheriche catchment (black) in the Sagarmatha National Park (green) in Nepal.Characteristics of the meteorological stations are summarized in Table 1.(b) Hypsometric curve of the Pheriche catchment.

Figure 2 .
Figure 2. Daily minimal and maximal air temperature and daily precipitation measured at the Pyramid station.
Glacier outlines characteristicsMass balances estimated bySherpa et al. (2017) for the clean-ice West Changri Nup and Pokalde glaciers located in the Pheriche basin (Fig.3) were used as reference, as well as mean annual glacier mass balances calculated over the Pheriche basin area for the period 2000-2016 byBrun et al. (2017).

Figure 3 .
Figure 3. Glacier outlines in the Pheriche catchment.(a) Clean glaciers and debris-covered glaciers from Racoviteanu et al. (2013) and location of the clean ice West Changri Nup and Pokalde glaciers (b) GAMDAM (red) and ICIMOD (blue) glacier outlines to simulate glacier mass 6 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-34Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 9 February 2018 c Author(s) 2018.CC BY 4.0 License.balance and the runoff production in catchments with glaciers, thus extending the application to ice dominated hydrological regimes.The resulting DHSVM-GDM simulates the spatial distribution and the temporal evolution of the principal water balance terms (soil moisture, evapotranspiration, sublimation, glacier mass balance, snow cover, and runoff) at hourly to daily time scales.It uses a two-layer energy and mass balance module for simulating snow cover evolution and a single layer energy

Figure 4 .
Figure 4. Original and modified parametrization of the snow albedo evolution in DHSVM-GDM and comparison with observed albedo data (2010-2015) in Pheriche, Pyramide and Changri Nup.

Figure 5 .
Figure 5. Definitions of the flow components.

-
v1: modified snow albedo parametrization; -v2: modified snow albedo parametrization and avalanche module; -v3: modified snow albedo parametrization, avalanche module and melt coefficient for debris covered glaciers.All four configurations were run with the Racoviteanu et al. (2013) glaciers outline.Concerning the melt of the debris-covered glaciers, we use a reduction factor of 0.4 as estimated by Vincent et al. (2016) from a study on uncovered and debris covered areas of the Changri Nup glacier.

Figure 6 .
Figure 6.Daily temperature and precipitation lapse rates.Discarded precipitation lapse rates (with a R 2 <0.75) are represented in orange.
Figures 7 and 8 compare the simulated snow cover area (SCA) and duration obtained with the configurations v0, v1, v2, and v3 to data derived from MODIS images.The SCA is clearly overestimated using the original parametrization v0: Figure75

Figure 7 .
Figure 7.Comparison of the MODIS SCA and the simulated daily SCA with the four modelling configurations (v0, v1, v2 and v3) for the Pheriche catchment.

Figure 8 .
Figure 8. Difference between the mean annual snow cover duration simulated with DHSVM-GDM and derived from MODIS images (in days) for the Pheriche catchment (top panels), with a focus on West Changri Nup (medium panels) and Pokalde glaciers (bottom panels).
et al. (2017) represent the mean mass balance for all the glaciers located in the Pheriche basin.Moreover, the mean annual glacier mass balance is estimated here from a three years simulation run (2012-2015), whereas the geodetic mass balances 14 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-34Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 9 February 2018 c Author(s) 2018.CC BY 4.0 License.were estimated for longer intervals (5 to 15 years) covering different periods and the inter-annual variability is thus not taken into account.

Figure 9 .
Figure 9. Mean annual glacier mass balances simulated with configurations v0, v1, v2 and v3.The error bar for configuration v3 represents the uncertainty related to the debris layer coefficient melt varying between 0.3 and 0.5.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-34Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 9 February 2018 c Author(s) 2018.CC BY 4.0 License.(including errors depending on the interpolated input fields and errors induced by the representation of slopes and expositions by the DEM) and/or from a lack of representativeness of the measurements.

Figure 10 .
Figure 10.Annual simulated and measured point mass balances on West Changri Nup (left panel) and Pokalde (right panel) glaciers; also shown is the 1:1 line.

Figure 11 .
Figure 11.Simulated annual hydrological contributions to Pheriche outflow for the two definitions of the flow components: (a) definition 1, (b) definition 2 and for 3 glaciological years from 12/2013 to 11/2015.
presents the daily simulated discharges and the flow components estimated with the two different definitions (runoff production vs. contributive areas).Daily discharges are well simulated in 2013 and 2015 by the model, with NSE equal to 0.9 and KGE equal to 0.86.However, outflow is under-estimated by the model during the monsoon season in 2014.The simulated total runoff is higher than the outflow at the catchment outlet before the monsoon season (from February to June) and lower during post-monsoon and 20 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-34Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 9 February 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 12 .
Figure 12.Daily discharges and flow components simulated with configuration v3: (a) production of ice melt, snow melt and net rainfall (note that the sum of the flow components represent the total runoff and is not equal to the discharge at the catchmant outlet, see definition 1 at section 3.2.5)(b) hydrological contributions to the outflow (definition 2 at Sect. 3.2.5).Observed discharges are represented by the black line with a 15% interval of uncertainty.

Figure 13 .
Figure 13.Monthly share of the runoff production (a) and hydrological contributions (b) simulated with configuration v3.Water balance contributions are introduced at section 3.2.5.

Figure 14 .
Figure 14.Mean hourly precipitation, discharge and flow components by seasons (DJFM: winter, AM: pre-monsoon, JJAS: monsoon, ON: post-monsoon) simulated with configuration v3.Water balance contributions are introduced at Sect. 3.2.5.Note that the y-axis scale is different for each season.

Figure 14
Figure 14 presents the different diurnal cycles of precipitation and hydrological components for each considered season (winter, pre-monsoon, monsoon, post-monsoon) obtained with configuration v3.During winter, pre-monsoon and post-monsoon, 10 observed discharges.At the daily scale (and longer time scales) the water balance is correctly simulated.However, at a sub-24 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-34Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 9 February 2018 c Author(s) 2018.CC BY 4.0 License.daily scale the model responds too quickly to the snow and ice melt production.Irvine-Fynn et al. (2017) found that on the Khumbu glacier the presence of supraglacial ponds buffers the runoff by storing diurnally more than 20% of the discharge.This could explain the longer transfer time observed on the measured outflow which are not represented by the model.This also shows that the current representation of the glacier and soil storage in DHSVM-GDM does not allow to reproduce accurately the diurnal variations of discharge and further studies are needed in order to determine the causes of such a shortcoming.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-34Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 9 February 2018 c Author(s) 2018.CC BY 4.0 License.glacial contribution as the total outflow from the glacierized area or as outflow produced by ice melt can lead to very different estimations.

Table 1 .
Location of measurements.T air temperature, P precipitation, WS wind speed, RH relative humidity, SWin incoming shortwave radiation, SWout outgoing shortwave radiation, LWin incoming longwave radiation.
-1 ], and i snowf all the snowfall intensity [mm/h].Since the observed decrease is dependent on elevation, the coefficient c is calculated as a function of elevation according to Eq. 3:

Table 3 .
Annual hydrological balance for configuration v3 and sensitivity to the debris covered melt reduction factor ranging from 0.3 to 0.5 (from left to right), for 3 glaciological years from December 2013 to November 2015.

Table 4 .
Mean annual glaciers mass balance (MB), outflow and flow components simulated with different glaciers inventories (configuration v3)