Quantification of different flow components in a high-altitude glacierized catchment (Dudh Koshi, Himalaya): some cryospheric-related issues

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 resource availability and planning the future uses of water in downstream regions. Two of the main issues in the hydrology of high-altitude glacierized catchments are (i) the limited representation of cryospheric processes controlling the evolution of ice and snow in distributed hydrological models and (ii) the difficulty in defining and quantifying the hydrological contributions to the river outflow. This study estimates the relative contribution of rainfall, glaciers, and snowmelt to the Khumbu River streamflow (Upper Dudh Koshi, Nepal, 146 km², 43 % glacierized, elevation range from 4260 to 8848 m a.s.l.) as well as the seasonal, daily, and sub-daily variability during the period 2012–2015 by using the DHSVM-GDM (Distributed Hydrological Soil Vegetation Model – Glaciers Dynamics Model) physically based glacio-hydrological model. The impact of different snow and glacier parameterizations was tested by modifying the snow albedo parameterization, adding an avalanche module, adding a reduction factor for the melt of debris-covered glaciers, and adding a conceptual englacial storage. The representation of snow, glacier, and hydrological processes was evaluated using three types of data (MODIS satellite images, glacier mass balances, and in situ discharge measurements). The relative flow components were estimated using two different definitions based on the water inputs and contributing areas. The simulated hydrological contributions differ not only depending on the used models and implemented processes, but also on different definitions of the estimated flow components. In the presented case study, ice melt and snowmelt contribute each more than 40 % to the annual water inputs and 69 % of the annual stream flow originates from glacierized areas. The analysis of the seasonal contributions highlights that ice melt and snowmelt as well as rain contribute to monsoon flows in similar proportions and that winter outflow is mainly controlled by the release from the englacial water storage. The choice of a given parametrization for snow and glacier processes, as well as their relative parameter values, has a significant impact on the simulated water balance: for instance, the different tested parameterizations led to ice melt contributions ranging from 42 % to 54 %. The sensitivity of the model to the glacier inventory was also tested, demonstrating that the uncertainty related to the glacierized surface leads to an uncertainty of 20 % for the simulated ice melt component.

1 Introduction

The Himalayan mountain range is known for being the water tower of Central and South Asia (Immerzeel et al., 2010). Its high-elevation 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 significant impact on the hydrological regime (Akhtar et al., 2008; Immerzeel et al., 2012; Lutz et al., 2014; Nepal, 2016). Development of tourism is also affecting the accessibility to water during the peak of the tourist season. In the Everest region in Nepal water needs have increased within the past decades due to higher demand in water supply for tourists and hydroelectricity production, leading to water shortages during months with low flows (winter and spring) (McDowell et al., 2013). Understanding the past and present hydrological regimes, and more particularly estimating the seasonal contributions of ice melt, snowmelt, and rainfall to outflows, is thus a key issue for managing water resources within the next decades. Indeed, the quantification of the ice melt contribution enables us to assess the proportion of water currently available which is coming from a long-term accumulation in the glaciers and thus to assess the annual decrease in the basin water storage due to glacier melt. Moreover, knowing the fraction of snowmelt, ice melt, and rainfall to the river outflow, and understanding their hydrological pathways, can give insights into how much water is currently seasonally delayed and how the seasonal outflow and the overall water balance might be impacted in the future when this delay changes or if the ratio of snowfall to rainfall changes (Berghuijs et al., 2014).

Recent studies have estimated present glacier and snowmelt contributions to the outflow in Nepalese Himalayan catchments (e.g., Andermann et al., 2012; Savéan et al., 2015; Ragettli et al., 2015) and simulated future hydrological regimes using glacio-hydrological models (Rees and Collins, 2006; Nepal, 2016; Soncini et al., 2016). Results have demonstrated large differences in the estimates of the contribution of glaciers to the annual outflows of the Dudh Koshi catchment in Nepal, which range from 4 % to 60 % (Andermann et al., 2012; Racoviteanu et al., 2013; Nepal et al., 2014; Savéan et al., 2015).

One of the main sources of uncertainty in modeling the outflow of Himalayan catchments is the representation of cryospheric processes, which control the evolution of ice and snow-covered surfaces in hydrological models. For instance, the representation of the debris-covered glaciers in glacio-hydrological models is a challenge. Debris-covered glaciers represent about 23 % of all glaciers in the Himalaya–Karakoram region (Scherler et al., 2011). The debris layers have been expanding during the last decades due to the glacier recession (Shukla et al., 2009; Bhambri et al., 2011; Benn et al., 2012) and are expected to keep expanding in the near future (Rowan et al., 2015). Since the study of Østrem (1959) it has been known that the debris thickness has a strong impact on the meltwater generation, which means that a good representation of the debris-covered glaciers in glacio-hydrological models is essential for estimating the amount of meltwater generated in glacierized catchments in the Himalayas. Many other cryospheric processes, such as the liquid water storage and transfer through glaciers, snow transport by avalanches or wind, glacial lake dynamics, and snow albedo evolution are either very simplified or not at all represented by the models (Chen et al., 2017). It is therefore important to estimate the impact of such simplified representations of cryospheric processes on modeling results.

Delineation of the glacierized areas is another key entry element in the glacio-hydrological model. Glacier inventories are commonly used as forcing data to delineate glacierized areas in glacio-hydrological modeling studies. There are three global major glacier inventories, i.e., the World Glacier Inventory (Cogley, 2009), GlobGlacier (Paul et al., 2009), and the Randolph Glacier Inventory (Pfeffer et al., 2014), and several regional glacier inventories in the HKH region (ICIMOD, Bajracharya et al., 2010; Racoviteanu et al., 2013), showing substantial differences. These can be due to the definition of the glacierized area itself (Paul et al., 2013; Brun et al., 2017) as well as to the characteristics of the satellite image (date, resolution, spectral properties) used for the delineation (Kääb et al., 2015), and to difficulties related to the interpretation of satellite images for outlying the glaciers, especially when they are debris-covered (Bhambri et al., 2011; Racoviteanu et al., 2013; Robson et al., 2015). Thus, the question of whether the glacier delineation has a significant impact on the model results needs to be addressed.

These issues of the representation of cryospheric processes and of glacier delineation in the hydrological modeling are addressed in the present study by (i) adapting the parameterization of the snow albedo evolution of DHSVM-GDM in order to improve the simulation of the snow cover dynamics; (ii) implementing an avalanche module; (iii) introducing a melting factor for debris-covered glaciers; and (iv) testing the sensitivity of simulated outflows and flow components with respect to these modifications as well as to glacier delineation for three different outlines coming from different glacier inventories. Both in situ measurements and satellite data were used for evaluating the outflow simulations as well as snow cover and glacier evolutions focusing on a small headwater catchment.

Figure 1 Study area: (a) location map of the 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 Daily minimal and maximal air temperatures and daily precipitation measured at the Pyramid station.

There are indeed several ways to define the glacier contribution to runoff (Radić and Hock, 2014): it can be considered either the total outflow coming from glacierized areas, the outflow produced by the glacier itself (snow, firn, and ice melt), or the outflow produced only by the ice melt. How the contributions to the outflow are defined adds further uncertainty to the estimation of the glacier contribution. The definition of the glacial contribution is dependent on the hydrological model (distributed or lumped, representation of glaciers and snow in the model) and cannot always be chosen. In the Dudh Koshi basin, Andermann et al. (2012); Racoviteanu et al. (2013); Savéan et al. (2015) estimated the fraction of the outflow produced by ice melt, whereas Nepal et al. (2014) defined the glacier contribution as the fraction of the outflow coming from glacierized areas. Here, flow components were estimated using two different definitions of the hydrological contributions for assessing their relative contributions to the total water balance. Finally, model results are analyzed at the annual, monthly, daily, and sub-daily scales in order to explain the origin of the water flows and their seasonal and daily variations.

2 Study area

This study focuses on the Pheriche sub-catchment of the Dudh Koshi basin (outlet at coordinates 27.89^∘ N, 86.82^∘ E) located in Nepal on the southern slopes of Mt. Everest in the Sagarmatha National Park (SNP) (Fig. 1). The catchment area is 146 km² and its elevation extends from 4260 to 8848 m a.s.l.

Local climate is mainly controlled by the Indian summer monsoon (Bookhagen and Burbank, 2006) and is characterized by four different seasons: a cold dry winter from December to March with limited precipitation, a warm and moist summer with most of the annual precipitation occurring during the monsoon from June until September, and two transition seasons: the pre-monsoon season in April and May and the post-monsoon season in October and November (Shrestha et al., 2000). At 5000 m, the annual precipitation is around 600 mm and the mean monthly temperature ranges from −8.4 ^∘C in January to 3.5 ^∘C in July, according to temperature and precipitation data from the Pyramid EvK2 station (Fig. 2 and Table 1). The hydrological regime follows the precipitation cycle with high flows during the monsoon season, when most of the annual precipitation occurs, complemented by the melting of snow and ice, and low flows during winter.

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, and (b) GAMDAM (red) and ICIMOD (blue) glacier outlines.

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.

Due to high elevation, vegetation in the catchment is scarce. The basin area is mainly covered by rocks and moraines (43 %) (Bajracharya et al., 2010) and glaciers (43 %) (Racoviteanu et al., 2013). Only 14 % of the basin area is covered by grasslands and shrublands. Glaciers belong to the summer-accumulation type (Wagnon et al., 2013) and are partially fed by avalanches (King et al., 2017; Sherpa et al., 2017); 60 % of the glaciers are located between 5000 and 6000 m a.s.l. Debris-covered glaciers are found at low elevations mainly on the ablation tongues of the glaciers (Fig. 3). According to the Racoviteanu et al. (2013) glacier inventory, debris-covered glaciers represent 30 % of the glacierized area, with smaller melting rates at similar elevations to debris-free glaciers due to the insulating effect of the debris layer (Vincent et al., 2016).

3 Data and model setup

3.1 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 land-cover classification from ICIMOD (Bajracharya, 2014) were used for the soil and land-cover 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, precipitation has been recorded at the Pheriche and Pyramid AWS by two Geonor T-200 sensors designed to measure both liquid and solid precipitation. 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, and shortwave and longwave radiation at Pheriche and Pyramid were provided by the EvK2-CNR stations.

Discharge measurements of the Khumbu River at Pheriche have been obtained using a pressure water level sensor at a 30 min time step since October 2010.

Racoviteanu et al. (2013)(Nuimura et al., 2015)(Bajracharya et al., 2010)

Table 2 Glacier outline characteristics.

The MODImLab algorithm developed by Sirguey et al. (2009) was applied to MODIS reflectance 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 to atmospheric and topographic effects, which makes the snow cover maps more realistic in mountainous areas. Twenty-seven 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).

To describe 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 to different results.

Mass balances estimated by Sherpa et al. (2017) for the clean-ice West Changri Nup and Pokalde glaciers located in the Pheriche basin (Fig. 3) were used as a reference, as well as mean annual glacier mass balances calculated over the Pheriche basin area for the period 2000–2016 by Brun et al. (2017).

3.2 Glacio-hydrological modeling

3.2.1 General description of the model

The DHSVM-GDM (Distributed Hydrological Soil Vegetation Model – Glaciers Dynamics Model) glacio-hydrological model 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; Beckers and Alila, 2004). A glacier dynamics module was recently implemented in DHSVM by Naz et al. (2014) to simulate glacier mass 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 timescales. It uses a two-layer energy and mass balance module for simulating snow cover evolution and a single-layer energy 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 region 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, land cover, soil depth, and ice thickness).

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

3.2.2 Snow albedo parameterization

In the original DHSVM-GDM version, the snow albedo α_s (–) is set to its maximum value ${{\mathit{\alpha }}_{\text{s}}}_{\text{max}}$ (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):

$$\begin{array}{}\text{(1)}& \begin{array}{r}{\mathit{\alpha }}_{\text{s}}={{\mathit{\alpha }}_{\text{s}}}_{\text{max}}({\mathit{\lambda }}_{\text{a}}{)}^{N^{\mathrm{0.58}}}\text{if}{T}_{\text{s}}<\mathrm{0},\\ {\mathit{\alpha }}_{\text{s}}={{\mathit{\alpha }}_{\text{s}}}_{\text{max}}({\mathit{\lambda }}_{\text{m}}{)}^{N^{\mathrm{0.46}}}\text{if}{T}_{\text{s}}\mathit{⩾}\mathrm{0},\end{array}\end{array}$$

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 (^∘C).

MODIS albedo images and the albedo measurements from Pyramid and Changri Nup were used to analyze the decrease in 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 3 consecutive days without clouds after the snowfall with the albedo parameterization in DHSVM-GDM. Since the observed values are not well represented by the standard albedo decrease, the parameterization was replaced by Eq. (2), with a decay of the albedo when there is no new snowfall inspired by the ISBA model albedo parameterization (Douville et al., 1995) and with the fresh snow albedo modified as a function of the amount of snowfall:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{\alpha }}_{\text{s}}& =({{\mathit{\alpha }}_{\text{s}}}_{t-\mathrm{1}}-{{\mathit{\alpha }}_{\text{s}}}_{\min })\exp (-cN)+{{\mathit{\alpha }}_{\text{s}}}_{\min }\\ & \mathrm{if}{i}_{\mathrm{snowfall}}=\mathrm{0}\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}},\\ {\mathit{\alpha }}_{\text{s}}& =\max (\mathrm{0.6},{{\mathit{\alpha }}_{\text{s}}}_{t-\mathrm{1}})\\ & \mathrm{if}\mathrm{0}\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}<i_{\mathrm{snowfall}}\mathit{⩽}\mathrm{1}\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}},\\ {\mathit{\alpha }}_{\text{s}}& =\max (\mathrm{0.6},{{\mathit{\alpha }}_{\text{s}}}_{t-\mathrm{1}})+({{\mathit{\alpha }}_{\text{s}}}_{\max }-max(\mathrm{0.6},{{\mathit{\alpha }}_{\text{s}}}_{t-\mathrm{1}}\left)\right){\displaystyle \frac{i_{\mathrm{snowfall}}-\mathrm{1}}{\mathrm{3}-\mathrm{1}}}\\ & \mathrm{if}\mathrm{1}\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}<i_{\text{snowfall}}\mathit{⩽}\mathrm{3}\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}},\\ {\mathit{\alpha }}_{\text{s}}& ={{\mathit{\alpha }}_{\text{s}}}_{\max }\\ & \mathrm{if}{i}_{\mathrm{snowfall}}>\mathrm{3}\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}},\end{array}\end{array}$$

where ${{\mathit{\alpha }}_{\text{s}}}_{t-\mathrm{1}}$ is the albedo from the previous time step, ${{\mathit{\alpha }}_{\text{s}}}_{\min }$ 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 (d⁻¹), and i_snowfall 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):

$$\begin{array}{}\text{(3)}& c=\mathrm{20}\exp (-\mathrm{0.001}Z),\end{array}$$

where Z is the elevation of the cell in m a.s.l.

The new function for the decrease in the snow albedo is also shown in Fig. 4.

3.2.3 Avalanche parameterization

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-elevation cells, where the air temperature remains below 0 ^∘C, and to a deficit of snow in the lower areas, where snowmelt 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 was implemented in DHSVM-GDM. The avalanche model transfers snow to downslope neighbor 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 neighbor cells is larger than 50 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 neighbor cells (between one and four cells). In case of several possible directions downstream, avalanche snow is distributed according to a ratio based on the slope of each of the directions. Within the same time step, the amount of snow in the receiving cells is actualized and the avalanches propagate downslope until the conditions cited above are no longer respected.

3.2.4 Glacier parameterization

Distributed ice thickness is derived from the terrain slope following the method described in Haeberli and Hölzle (1995).

In the original DHSVM-GDM version, glacier melt is instantaneously transferred to the soil surface, which is parameterized as bedrock under glaciers (Naz et al., 2014). This significantly underestimates the transfer time through glaciers. In this study, storage of liquid water inside glaciers was implemented by adding an englacial porous layer between the glacier and the bedrock allowing the liquid water storage within the glacier. Previous studies have shown the good performance of adding a conceptual representation of the storage and drainage in the glaciers within glacio-hydrological models (e.g., Jansson et al., 2003; Hock and Jansson, 2006). The most widely adopted approach is based on a reservoir or a cascade of reservoirs with time-invariant parameters (e.g., Farinotti et al., 2012; Zhang et al., 2015; Hanzer et al., 2016; Gao et al., 2017). Here, storage of liquid water inside glaciers was implemented by adding an englacial porous layer between the glacier and the bedrock allowing the liquid water storage within the glacier. This englacial porous layer has a depth of 2 m and is characterized by a porosity of 0.8 and a hydraulic conductivity (vertical and lateral) of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹ (see Table A2). As in the previously cited studies, the parameters are kept constant through the simulations. They were optimized here according to the constraint of minimizing the differences (in terms of least squares) between the recession shape of the simulated hydrographs and the observed one.

Moreover, 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).

3.2.5 Quantification of the flow components

Quantifying the relative contributions of ice melt, snowmelt, and rainfall in the river discharge at different timescales is a difficult task because hydrological models usually do not track the origin of water during transfer within the catchment (Weiler et al., 2018). There are also different ways of defining the origins of the streamflow. Weiler et al. (2018) lists three types of contributions: (1) contributions from the source areas, i.e. from each class of land cover, (2) contributions from the runoff generation (overland flow, subsurface flow, and groundwater flow), and (3) input contributions (ice melt, snowmelt, and rain).

In this study, two different definitions were used to estimate the hydrological contributions. First, we estimate the contributions of ice melt, snowmelt, and net rainfall to the total water input (definition 1) according to the following water balance equations (all the terms are fluxes expressed in L T⁻¹):

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{Input}=\mathrm{Icemelt}+\mathrm{Snowmelt}+\mathrm{RainNet},\text{(5)}& {\displaystyle }& {\displaystyle \frac{\text{d}{I}_{\mathrm{wq}}}{\text{d}t}}=\mathrm{GlAcc}-\mathrm{Icemelt}-\mathrm{SublIce},\text{(6)}& {\displaystyle }& {\displaystyle \frac{\text{d}{S}_{\mathrm{wq}}}{\text{d}t}}=P_{\text{solid}}-\mathrm{Snowmelt}-\mathrm{SublIce}-\mathrm{GlAcc},\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{RainNet}=P_{\text{liquid}}-E_{\text{int}},\text{(8)}& {\displaystyle }& {\displaystyle \frac{\text{d}\mathrm{Storage}}{\text{d}t}}=\mathrm{Input}-Q-E_{\text{T}},\end{array}$$

where $\frac{\text{d}{I}_{\mathrm{wq}}}{\text{d}t}$ and $\frac{\text{d}{S}_{\mathrm{wq}}}{\text{d}t}$ are the variations of the ice and snow storages, GlAcc is the amount of snow that is transferred to the ice layer by compaction on glaciers (Naz et al., 2014), SublIce and SublSnow are the amounts of sublimation from the ice and snow layers, 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 the sum of these contributions (Input) is not equal to the outflow at the catchment outlet Q as it represents all liquid water reaching the soil surface (before infiltration and potential storage in the soils and glaciers $\left(\frac{\text{d}\mathrm{Storage}}{\text{d}t}\right)$ and before evapotranspiration (E_T)).

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 contributing areas (definition 2):

• direct glacier contribution: direct runoff from glacierized areas;

• delayed glacier contribution: resurging meltwater stored inside glaciers;

• direct snow contribution: direct outflow from snow-covered non-glacierized areas;

• direct runoff: direct runoff from areas without snow and glaciers; and

• subsurface and groundwater contribution: resurging water from the soil in non-glacierized areas resulting from infiltrated rainfall, snowmelt, as well as upstream lateral subsurface flows.

These contributions are obtained from the amount of water reaching the soil surface simulated by DHSVM-GDM (see Supplement). On each grid cell, this volume is a mixture of ice melt, snowmelt, and rainfall and can either infiltrate into the soil or produce runoff. Definition 2 combines contributions from source areas (glacierized and non-glacierized areas) and contributions from runoff generation (direct runoff, englacial contribution, and soil contribution).

Figure 5 Definitions of the flow components.

Figure 5 illustrates the two definitions of the different contributions to outflows. Definition 1 allows assessment of the annual impact of glacier melt and snowmelt on the water production, while definition 2 describes the intra-annual routing of the water within the catchment. Moreover, using the two definitions allows us to directly compare the results with other hydrological modeling studies in the Dudh Koshi basin, which have estimated glacier 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 glacier 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 estimated hydrological contributions.

3.3 Experimental setup

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 the 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 result, soil depth outside glacierized areas ranges between 0.5 and 1 m (glaciers are considered to lay on bedrock). 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 modeling, we performed simulations with the four following configurations:

• v0: original DHSVM-GDM snow and glacier parameterization;

• v1: modified snow albedo parameterization;

• v2: modified snow albedo parameterization and avalanche module;

• v3: modified snow albedo parameterization, avalanche module, and melt coefficient for debris-covered glaciers.

All four configurations were run with the Racoviteanu et al. (2013) glacier 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.

Using configuration v3, we also tested the impact of using different glacier outlines (the GAMDAM and ICIMOD inventories were also considered for simulations) and analyzed the sensitivity related to different values of parameters related to the cryospheric processes (see Table B1). Indeed, the debris-covered glacier melt reduction factors 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). A sensitivity analysis of the englacial porous layer parameters (depth, porosity, and hydraulic conductivity) and avalanche parameters was also performed (see Table B1 for the tested parameter values) and the relative impact on the simulated hydrological response was discussed (see Sect. 5.3.2).

The sensitivity to the values of the englacial porous layer parameters (depth, porosity, and hydraulic conductivity), as well as to the values of the avalanche parameters and to the soil depth is also analyzed in the discussion section of the paper.

3.3.1 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 a 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² values higher than 0.75 were retained for precipitation (43 % of the dataset). For smaller R², the lapse rate is considered 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⁻¹. 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⁻¹. We found only 10 days (1 %) showing a temperature inversion with a positive daily lapse rate.

Figure 6 Daily temperature and precipitation lapse rates. Discarded precipitation lapse rates (with a R²<0.75) are represented in orange.

3.3.2 Model evaluation

A multi-criteria evaluation was made considering simulated outflows, SCA, and glacier mass balances. Discharge measurements at Pheriche station were used as a 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 in 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 Pokalde glaciers (Sherpa et al., 2017).

4 Results

4.1 Impact of the snow and glacier parameterizations on the simulated results

This section presents the simulation results obtained with the different configurations of the DHSVM-GDM model (configurations v0, v1, v2, and v3; see Sect. 3.2.5) and the analysis of the impact of the snow and glacier parameterizations on the simulated annual outflow, the daily SCA, and annual glacier mass balances.

Figure 7 Simulated annual hydrological contributions (definition 1) to Pheriche outflow for 3 glaciological years from December 2013 to November 2015.

4.1.1 Annual outflow

Figure 7 represents the annual outflow and flow components (definition 1) simulated with the different model configurations, indicating the impact of each modification of the snow and glacier parameterization on the simulated annual outflow and flow components. Configuration v1 leads to a drastically increased outflow due to an enhanced ice melt component. Implementing the avalanche module (v2) reduces the ice melt component and increases the snowmelt component by 21 %. Configuration v3 including debris-covered glaciers further reduces the ice melt, resulting in a simulated annual outflow close to the observations.

Figure 7 shows that 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 817 mm, which is nearly twice the amount of ice melt obtained with v3 that includes debris-covered glacier melt.

The configuration with all three modifications (v3) gives results similar to the original parameterization of DHSVM-GDM (v0) in terms of glacier mass balance, improving slightly the annual outflow. The ice melt factor for debris-covered glaciers and the avalanches compensate for the increase in ice melt caused by the new snow albedo parameterization, but the modifications implemented in v3 impact the results for the flow components: on average, less ice melt and more snowmelt are generated. Moreover, configuration v3 modifies the seasonal variation of the outflow by increasing winter discharges and reducing monsoon discharges (not illustrated here), which improves the daily NSE and KGE (Table 3).

Table 3 NSE and KGE values calculated on the daily discharges on the period 2012–2015 for each model configuration.

4.1.2 Snow cover dynamics

Figures 8 and 9 compare the simulated SCA and duration obtained with configurations v0, v1, v2, and v3 to data derived from MODIS images. The SCA is strongly overestimated using the original parameterization v0: Fig. 8 shows that after full coverage it does not decrease fast enough compared to MODIS data. Figure 9 demonstrates that the snow cover duration is overestimated for the entire catchment area. This indicates that in the simulations snow does not melt fast enough with the original parameterization. Configuration v1 with the modified snow albedo parameterization (Eq. 2) accelerates the snowmelt and improves the SCA simulation (Fig. 8). The RMSE between the simulated and observed SCA decreases from 29 % using v0 to 14 % using v1 and v2. Figure 9 shows that with configuration v1 in some areas located at high elevation the snow cover duration is underestimated. This bias is rectified in configuration v2 since the avalanche module transfers snow from high-elevation and sloping cells downward and corrects the lack of snow observed with configuration v1 at the edges of the permanent snow cover (Fig. 9). The results for the SCA and snow cover duration using configuration v3 show no difference compared to configuration v2 since only the ice melt rate for debris-covered glaciers is modified.

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

Figure 9 Difference between the mean annual snow cover duration simulated with DHSVM-GDM and derived from MODIS images (in days) for the Pheriche catchment (a–c), with a focus on West Changri Nup (d–f) and Pokalde glaciers (g–i).

4.1.3 Glacier mass balances

Figure 10 compares the simulated mean annual glacier mass balances obtained with the different model configurations with mass balances determined with geodetic methods between 1999 and 2015 (Bolch et al., 2012; Gardelle et al., 2013; Nuimura et al., 2015; King et al., 2017; Brun et al., 2017). These geodetic mass balances range from $-\mathrm{0.67}\pm \mathrm{0.45}$ m w.e. yr⁻¹ (2000–2008) (Nuimura et al., 2015) to $-\mathrm{0.32}\pm \mathrm{0.09}$ m w.e. yr⁻¹ (2000–2015) (Brun et al., 2017).

Our results show that the snow parameterization 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⁻¹ and decreases to −2.02 m w.e. yr⁻¹ with the corrected snow albedo (v1) since the modified snow albedo parameterization accelerates the snowmelt, 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⁻¹. 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 $-\mathrm{0.84}\pm \mathrm{0.14}$ m w.e. yr⁻¹, which is within the intervals of uncertainty and thus in good agreement with geodetic methods.

Figure 10 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.

Figure 11 Annual simulated and measured point mass balances on West Changri Nup (a) and Pokalde (b) glaciers; also shown is the 1 : 1 line.

We also evaluated the mass balance at the point scale. Figure 11 shows the simulated mass balances with parameterizations v0, v1, v2, and v3 versus the observed mass balances of the two debris-free glaciers West Changri Nup and Pokalde measured in situ for the 3 glaciological years (2012–2015). Here, configuration v3 gives the same results as configuration v2 because in configuration v3 only the ice melt rate on debris-covered glaciers is modified. The simulated mass balances vary according to the model configuration. With configuration v0, the model overestimates the point mass balances because of small snowmelt rates (see also Sect. 4.1.2). With configuration v1, the model overestimates the 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. 9). Configuration v2 improves the simulated mass balance by transferring snow due to avalanches on the West Changri Nup glacier and by removing excessive snow accumulation on the Pokalde glacier. For the Pokalde glacier, the mass balances simulated with configuration v3 show a larger variability than the mass balances simulated with configuration v0, but the point mass balances are spread around the diagonal axis, which leads to a bias 10 times smaller (mean bias of 1 m with v0 and 0.1 m with v3).

The results at basin scale and point scale show that the snow parameterization has a strong impact on the simulated glacier mass balance and that the new snow albedo parameterization 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 (including errors depending on the interpolated input fields and errors induced by the representation of slopes and expositions by the DEM) or a lack of representativeness of the measurements.

Figure 12 Simulated annual hydrological contributions to Pheriche outflow for the two definitions of the flow components (definition 1 and definition 2) and for 3 glaciological years from December 2013 to November 2015.

Table 4 Annual hydrological balance simulated with configuration v3 for 3 glaciological years from December 2013 to November 2015.

4.2 Simulated outflow and flow components

This section presents the outflows and flow components simulated in the Pheriche basin during the period 2012–2015 with the modified version of DHSVM-GDM (configuration v3). The simulation results are analyzed using two different definitions of the flow components (definitions 1 and 2; see Sect.3.2.5).

4.2.1 Annual simulated outflow and hydrological contributions

The annual outflows simulated with the new parametrization of the model (configuration v3) are in good agreement with the annual observed outflows since they remain within the 15 % interval of estimated error (Fig. 12 and Table 4).

The results show an inter-annual variability of the flow components. During the period 2013–2015, the ice melt component ranged from 41 % to 50 %, the snowmelt 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 decreased from 2013 to 2015 and explains the decrease in the net rainfall components from 155 mm in 2013 to 88 mm in 2015. The snowmelt component is higher in 2013 because of the 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 melted more rapidly and the glaciers started melting earlier. In contrast, 2015 was a year with a lot of winter snowfall, which delayed the beginning of the glacier melt and explains the lower ice melt component. The losses by evaporation and sublimation are rather constant through the simulation period, ranging from 140 to 150 mm yr⁻¹.

The runoff coefficients (ratio between the annual outflow and annual precipitation) were on average equal to 1.4, which means that a considerable amount of water is withdrawn each year from the catchment through ice melt (eventually in the form of a delayed groundwater flow).

On average, we find that the outflow is mainly produced by meltwater as 46 % of the annual water input is due to ice melt and 41 % to snowmelt (definition 1). The contributions estimated according to definition 2 show the importance of infiltration and subsurface flows in the water balance since more than 40 % of the outflow was coming from water infiltrated in glaciers and more than 20 % from subsurface and groundwater flows generated outside the glacier-covered area.

The choice of the definition of the hydrological components leads to different perceptions of the glacier contribution to the outflow. The glacier contribution to the total outflow is 69 % if the contribution from the entire glacierized area (i.e. contributions of ice melt, snowmelt, and net rainfall) is considered like in definition 2. However, the contribution from ice melt alone, included in definition 1, corresponds to only 46 % of the water input.

4.2.2 Seasonal variations of the flow components

Figure 13 presents the daily simulated discharges simulated with configuration v3 and the flow components estimated with the two different definitions. Daily discharges were well simulated for 2012–2013 and 2014–2015 by the model, with NSE equal to 0.91 and KGE equal to 0.88. However, the outflow is underestimated by the model during the monsoon season in 2014.

The simulated total water input (i.e. the sum of snowmelt, ice melt, and net rainfall) is always higher than the simulated outflow at the catchment outlet before the monsoon season (from February to June) and lower during the post-monsoon and winter seasons. This is mainly due to glacier meltwater stored inside the glaciers during the pre-monsoon and monsoon seasons and continuing surging during winter, as well as to changes in the soil water storage (Figs. 13b and 14b).

Figure 14 shows the mean monthly flow components averaged over the simulation period. From February to May–June, the water input is entirely controlled by snowmelt and ice melt (snowmelt between 50 % and 60 %, ice melt between 40 % and 48 %; Fig. 14a). The net rainfall, snowmelt, and ice melt absolute contributions are at their maxima in July and August during the monsoon season. During these 2 months, 24 % of the runoff is generated by net rainfall, 37 % by snowmelt, and 38 % by ice melt. From October to January, the runoff is produced by ice melt (up to 80 % in December) and snowmelt (between 20 % and 30 %). Groundwater and englacial water represent a significant fraction of the monthly outflow as they contribute more than 50 % of the outflow during the monsoon season and can contribute up to 90 % during winter (Fig. 14b). Direct contributions from glacierized areas, snow areas, and direct runoff are highest during the monsoon season, when the englacial and soil storage is saturated.

Figure 13 Daily discharges and flow components simulated with configuration v3: (a) production of ice melt, snowmelt, and net rainfall (note that the sum of the flow components represents the total water input and is not equal to the discharge at the catchment outlet; see definition 1, Sect. 3.2.5), (b) relative difference between the simulated glacier liquid water content and the inter-annual mean, and (c) hydrological contributions to the outflow (definition 2, Sect. 3.2.5). Observed discharges are represented by the black line with a 15 % interval of error.

Figure 14 Average monthly contributions to the water input (definition 1, Sect. 3.2.5) (a) and hydrological contributions (definition 2, Sect. 3.2.5) (b) simulated with configuration v3 for the years 2012–2015.

4.2.3 Diurnal cycle

Figure 15 presents the diurnal cycles of precipitation and hydrological components averaged for each considered season (winter, pre-monsoon, monsoon, and post-monsoon) obtained with configuration v3. During winter, pre-monsoon, and post-monsoon, the observed outflow is rather constant during the day, with a weak peak around noon when the temperature is at its maximum. During this period, almost all of the precipitation is in the form of snowfall leading to no direct response for the outflow. The peak around noon can be explained by snowmelt 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 causing a peak in the discharge around midnight.

The model simulates ice and snowmelt during day time with a maximum at noon as expected. Except for the monsoon season, it seems to simulate accurately the baseflow during night without melt production: the discharge is rather controlled by the release of the glacier and soil storage. However, the model simulates a peak of discharge around 14 h originating mostly from glacierized areas, 2 h after the maximum of ice and snowmelt, which does not correspond to observed discharges. At daily and longer timescales the water balance is correctly simulated. However, at a sub-daily scale the model responds too quickly to the snowmelt and ice melt production.

Figure 15 Mean hourly precipitation, discharge, and flow components simulated with configuration v3 and averaged for the winter (DJFM), pre-monsoon (AM), monsoon (JJAS), and post-monsoon (ON) seasons. Note the different y axis scales for each season.

5 Discussion

5.1 Simulation of the discharges and flow components

Soncini et al. (2016) studied flow components in the Pheriche catchment for the period 2013–2014 and estimated an annual ice melt component of 55 % and a snowmelt component of 20 % of the annual outflow. The ice melt components are thus quite similar in terms of relative contributions to outflow, which is not the case for the snowmelt components. We think that the main reason for 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 (Mimeau et al., 2019).

Concerning the seasonal contributions to the outflow, our results are consistent with the results from Soncini et al. (2016), who found a main contribution of snowmelt during the pre-monsoon season, mixed contributions of rainfall, snowmelt, and ice melt during the monsoon season, and mixed contributions of snowmelt and ice melt during the post-monsoon and winter seasons. The studies of Ragettli et al. (2015) and Racoviteanu et al. (2013) concerning the Upper Langtang and the Dudh Koshi basin, respectively, showed that most of the winter outflow surges from soil, channel, surface, and englacial storage changes, which is consistent with our results. However, the estimated flow components presented in this study, particularly the soil and englacial contributions, strongly depend on the model setup. Figure 13 shows that the main part of the soil infiltrated water resurges within a day, whereas liquid water can be stored for several months within the glaciers. This difference between the response of the soil storage and the englacial storage results from the soil and glacier parameterization (see the sensitivity analysis in Sect. 5.3.2).

At hourly scale, the results show that the model cannot represent the diurnal cycle of the outflow correctly as the simulated hydrological response is anticipated. 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 outflows which are not represented by the model. This shows that the current representation of the glacier storage in DHSVM-GDM does not allow us to reproduce accurately the diurnal variations of discharge, and further studies are needed in order to improve the model.

Overall, the comparison between the two definitions of the hydrological contributions shows that contributions must be explicitly specified in order to allow inter-comparison between models, especially for catchments with a large glacierized area. Moreover, the use of two different definitions allows us to get complementary information on the origin of the outflow (processes at the origin of the runoff, types of flow generation, contributive zones). A perspective to improve the quantification of the hydrological contributions to the outflow is to track the ice melt, snowmelt, and rainfall component pathways in the model as suggested in Weiler et al. (2018). This would enable us to quantify the fractions of the three components contributing to subsurface and groundwater flow, which is not possible with the current version of DHSVM-GDM.

5.2 Representation of the cryospheric processes in the model

One of the main difficulties in glacio-hydrological modeling is to correctly simulate both river discharges, snow cover dynamics, and glacier mass balances. The results showed that two different representations of the cryospheric processes in the model (v0 and v3) can lead to similar simulated annual outflows but different estimations of the ice melt and snowmelt contributions. This is particularly true for the glaciological year 2014–2015, when the ice melt contribution decreases from 59 % with v0 to 41 % with v3 and the snowmelt contribution increases from 29 % to 47 % (Fig. 7). This can be explained by the fact that 2014–2015 was a year with a high amount of snowfall (Fujita et al., 2017); therefore, the representation of snow processes in the model has a larger impact on the simulated runoff production than the 2 other years. This 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 in the flow component estimation.

The results also showed that the modification of one specific hydrological process (here, the representation of the snow albedo evolution) can have a significant impact on the simulated hydrological response of the catchment and requires improvement of other processes (here, considering specific representation of avalanches and debris-covered glaciers).

Further modifications of the model could also lead to different model results, and it is also not excluded that different model errors are compensating for each other. For example, the results showed that the original model version leads to a correct simulation of the river discharges because the non-representation of the insulation effect for debris-covered glaciers on the ice melt was compensated for by the incorrect representation of the snow albedo decrease. Due to the complexity of the model and the represented processes, no guarantee can be given that similar compensating effects still occur in the model. In this study, the validation of the model output was extended beyond the annual river discharge to discharges at different timescales, the snow-covered area, and glacier mass balances in order to validate the simulations of the snow cover, glacier melt, and discharges separately. The results demonstrate that the new version of the model performs well for all three signals. Moreover, the new parameterizations of the snow albedo and ice melt under debris were based on observed data (MODIS and in situ albedo measurements, and coefficient for ice melt under debris from Vincent et al., 2016) and do not result from a calibration in order to avoid compensation effects. Therefore, it is very likely that the new implementation improved the quality of the represented processes.

The results presented in this study also indicate potential future improvements to increase the reliability and to reduce the uncertainty of the simulations at short time steps. While at daily and longer timescales the different hydrological components seem to be well reproduced by the model, the analysis of the diurnal cycle (Fig. 15) shows that DHSVM-GDM responds too rapidly to the ice melt production, with too high diurnal peak discharges. This is probably related to the use of constant parameters in the parameterization of the englacial porous layer for glacier storage. Taking into account the seasonal variation of the efficiency of the englacial drainage system appears necessary to simulate the diurnal flow cycle correctly (Hannah and Gurnell, 2001). Therefore, further improvements should be based on studies analyzing the mechanisms of glacier drainage systems in the Khumbu region and their influence on glacier outflows (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 modeling should also be studied.

The avalanche routine implemented in this study is simplified and only considers four directions for the snow redistribution. A perspective of this study is to improve the representation of the avalanches in DHSVM-GDM by considering eight directions for the snow redistribution and considering other parameters such as the age of the snow cover, the snow density, and the type of land cover, as was proposed in Frey and Holzmann (2015).

5.3 Model sensitivities and uncertainties in forcing data

5.3.1 Sensitivity to the glacier outline

The three inventories used in this study result in very different estimates of the glacierized area: between 43 % and 24 % of the Pheriche basin with the inventories proposed by Racoviteanu et al. (2013) and GAMDAM (Fig. 3). Table 5 presents the average annual glacier mass balances, outflows, and flow components for configuration v3 using the three inventories. The GAMDAM inventory leads to a more negative glacier mass balance than the two other inventories, with −1.17 m w.e. yr⁻¹ compared to −0.84 and −0.89 m w.e. yr⁻¹ 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 melts than the Racoviteanu et al. (2013) inventory due to their smaller areas in ablation zones, which leads to a smaller simulated annual outflow. From these results we estimate an uncertainty of 20 % (430 mm with the Racoviteanu et al. (2013) inventory versus 345 mm with the ICIMOD inventory; cf. Table 5) in the simulated annual ice melt volume related to the glacier's 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) in the monthly discharges during the monsoon season. This result shows that the choice of the glacier inventory as an input datum of the glacio-hydrological model contributes to the uncertainty in the simulation 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 for a lack of precipitation due to 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, but the amounts of simulated ice melt are different, which shows that a consistent mass balance can lead to errors in the simulated glacier contributions and total outflow.

Racoviteanu et al. (2013)

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

5.3.2 Sensitivity to model parameter values

DHSVM-GDM is a physically based model with a large number of parameters which are difficult to define, particularly in a study area with little information on catchment characteristics and limited validation data. In order to have an indication of how sensitive the results are to the chosen parameter values, 27 additional simulations were performed with the new version of DHSVM-GDM (configuration v3). The model sensitivity concerning the melt coefficient on debris-covered glaciers, the soil depth, the parameters of the englacial porous layer, and the avalanche parameters (see Table B1) was tested.

Figure 16 Sensitivity of the simulated daily discharge to different parameter values related to cryospheric processes. Black line represents the reference simulation (configuration v3), grey area shows the range (minimum and maximum values) of the 27 simulations with changing parameters (see Table B1), and observed discharges are represented in red.

Figure 16 shows the impact of the parameter values on the simulated daily discharge. The overall variability ranges from 0.1 up to 7.4 m³ s⁻¹ and is larger during the monsoon season (3.8 m³ s⁻¹ on average in June). Overall, the parameter values have a low impact on the simulated daily discharge but a higher impact on the estimated relative contributions to the outflow (see Table B2). NSE values remain close to the reference simulation (around 0.91), while snowmelt and ice melt contributions (definition 1) vary from 34 % to 44 % and from 42 % to 54 %, respectively. Concerning contributing areas (definition 2), the parameter values mainly impact the distribution between direct and delayed runoff from glacierized areas.

A major limitation to the estimation of the contributions from the runoff generation (overland flow, subsurface and groundwater flow, and englacial flow) is the representation of the englacial flows in the model. The performed simulations show that the parameter values of the porous englacial layer have a strong impact on the results (Table B2, sim. 3 to 11). While it has no impact on the total annual contributions from glacierized and non-glacierized areas, it impacts the partitioning between direct and delayed glacier contributions ranging from 34 % to 20 % (direct) and 35 % to 47 % (delayed), with the tested range of parameter values. The overall performance of the simulation of the daily discharge is affected by the choice of the parameter values of the porous englacial layer (NSE values ranging between 0.84 and 0.9). Indeed, an increase in the glacier storage capacity leads to a rather important delay of the outflow (as there is more infiltration simulated during the pre-monsoon and monsoon seasons). It is worth noting that we kept constant parameters for the englacial porous layer in our simulations, which is clearly a simplification, at least considering seasonally evolving flow networks within the glacier. A perspective would be to consider seasonal parameters values in order to reduce the uncertainty regarding direct and delayed runoff generation, but data to validate englacial storage and glacier outflows are currently not available, making it very difficult to adapt the parameter values to the study area.

Another limitation of our model lies in the application of a uniform reduction factor for ice melt under debris-covered glaciers. Table B2 shows the sensitivity of the model to the ice melt reduction factor on debris-covered glaciers (sim. 1 and 2). When the reduction factor varies between 0.3 and 0.5, the simulated annual outflow is modified by ±7 % and the mean ice melt flow component ranges from 42 % to 50 %. 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 %.

We also tested the model sensitivity with respect to the avalanche parameters by changing the values of the thresholds for the slope (35, 40, and 45^∘) and the snow height (15, 30, and 60 cm) for steep terrains and for the snow water equivalent difference between two neighboring cells (25, 50, and 100 cm) in the case of a flat terrain (parameters used for these simulations are presented in Tables B2 and B1, simulations 12 to 27). The results confirmed that higher thresholds for triggering avalanches and their propagation lead to less accumulation on the downslope parts of the glaciers and, thus, to more overall ice melt. Nevertheless, the results also show that the simulated outflow is not very sensitive to the avalanche parameterization since we derived a variation of 3 % of the mean annual outflow. However, the modification of the avalanche parameters leads to an uncertainty of 12 % in the estimation of the ice melt and snowmelt components: the simulated contribution of ice melt to the water input varied between 46 % and 52 % and the contribution of snowmelt varied from 36 % to 41 %.

5.3.3 Validation and forcing data uncertainties

A main issue is related to the availability of data for validating the glacio-hydrological modeling parameterizations and outputs. The lack of in situ measurements at high elevations and the uncertainty related to the available data prevent us from assessing the performance of the model in simulating the different cryospheric processes we consider. They only allow us 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. For instance, concerning the validation of the simulated glacier mass balances, only the order of magnitude of the simulated and geodetic mass balances used in this study can be compared because the considered areas are not the same in the different studies, and the considered time periods differ as well. Indeed, four of the geodetic mass balances were derived for the Khumbu–Changri glacier, while mass balances from this study and from Brun et al. (2017) represent the mean mass balance for all glaciers located in the Pheriche basin. Moreover, the mean annual glacier mass balance is estimated here for the 3 simulated years (2012–2015), whereas the geodetic mass balances correspond to longer (5–15 years) as well as earlier periods beginning between 1999 and 2002 (see Fig. 10) and do not take into account the inter-annual variability of the glacier mass balances. Nevertheless, in Fig. 10, the variability of the simulated glacier mass balances is much larger than the variability of the geodetic mass balances, showing the significant impact of the snow and glacier parameterization on the simulation results. It is also worth noting that the snow cover distribution evaluation is particularly challenging in the Pheriche catchment, as clouds cover more than 50 % of the catchment during more than 150 d yr⁻¹ 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. In Mimeau et al. (2019), we tested the sensitivity of the model to different precipitation forcing datasets (in situ, reanalysis, and satellite) and analyzed the impact of different precipitation amounts and spatial distributions on the simulated discharges and flow components.

6 Conclusions

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 streamflows, which can be useful for water resources and water-related risk management. Different parameterizations of the snow and glacier processes in the model were compared in order to evaluate the uncertainty related to the representation of cryospheric processes in the simulated water balance. Simulated SCA was compared with MODIS images and calculated glacier mass balances with local in situ and geodetic measurements.

Results showed that the representation of the cryospheric processes in the model has a significant impact on the simulated outflow and flow components and that the improvement of one process' parameterization does not necessarily improve the overall simulation. Major contributions from glaciers and snow to the outflow were found, with 46 % of the annual water input produced by ice melt and 41 % by snowmelt in control simulation v3. Winter flows are mainly controlled by the release of englacial and soil water storage (up to 78 % in December in control simulation v3), which corroborates other studies (Racoviteanu et al., 2013; Ragettli et al., 2015). These results are sensitive to both the chosen cryospheric parameterizations (v0 to v3 setups) and their relative parameter values (27 simulations added to control v3), in a similar order of magnitude. For instance, 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⁻¹ for the annual glacier mass balance and ±7 % for the annual outflow. Given the different tested parameter values, the uncertainty on the annual ice melt and the snowmelt contributions ranges from −8 % to +17 % and from −17 % to +7 % with respect to the control contributions, respectively. When looking closely at the glacier contributions, these range from −23 % to +38 % and from −20 to +14 % for the direct and delayed contributions, respectively. Finally, the glacier inventories used to outline the glacierized areas also have 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.

This study also reminds us that glacial and snow contributions to the streamflow must be clearly defined, as considering glacial contribution to be the total outflow from the glacierized area or outflow produced by ice melt can lead to very different estimations.

Appendix A: Parameter values used in DHSVM-GDM

Brutsaert (2005)Brutsaert (2005)Andreadis et al. (2009)Brock et al. (2006)L'hôte et al. (2005)L'hôte et al. (2005)Singh (2001)Vincent et al. (2016)

Table A1 Global parameter values used for the glacio-hydrological simulation with DHSVM-GDM.

Clapp and Hornberger (1978)Niu et al. (2005)Morel-Seytoux and Nimmo (1999)Rawls et al. (1982)Rawls et al. (1982)Meyer et al. (1997)Meyer et al. (1997)Clapp and Hornberger (1978)Burns (2012)Burns (2012)Burns (2012)

Table A2 Soil and glacier parameter values used for the glacio-hydrological simulation with DHSVM-GDM.

Wigmosta et al. (1994)Wigmosta et al. (1994)Wigmosta et al. (1994)Dickinson et al. (1991)Wigmosta et al. (1994)Wigmosta et al. (1994)

Table A3 Vegetation parameter values used for the glacio-hydrological simulation with DHSVM-GDM.

Appendix B: Sensitivity analysis

Table B1 Parameter values used for the sensitivity analysis. The first line corresponds to the reference simulation (configuration v3) and simulations 1 to 27 refer to the sensitivity analysis simulations (parameters different from the reference simulation are highlighted).

Table B2 Nash–Sutcliffe coefficient (NSE) of simulated daily discharge and relative contributions to annual outflow (definitions 1 and 2) for the 27 sensitivity analysis simulations.