Providing a non-deterministic representation of spatial variability of precipitation in the Everest region

This paper provides a new representation of the effect of altitude on precipitation that represents spatial and temporal variability in precipitation in the Everest region. Exclusive observation data are used to infer a piecewise linear function for the relation between altitude and precipitation and significant seasonal variations are highlighted. An original ensemble approach is applied to provide non-deterministic water budgets for middle and high-mountain catchments. Physical processes at the soil– atmosphere interface are represented through the Interactions Soil–Biosphere–Atmosphere (ISBA) surface scheme. Uncertainties associated with the model parametrization are limited by the integration of in situ measurements of soils and vegetation properties. Uncertainties associated with the representation of the orographic effect are shown to account for up to 16 % of annual total precipitation. Annual evapotranspiration is shown to represent 26 %± 1 % of annual total precipitation for the mid-altitude catchment and 34%± 3 % for the high-altitude catchment. Snowfall contribution is shown to be neglectable for the mid-altitude catchment, and it represents up to 44 %± 8 % of total precipitation for the highaltitude catchment. These simulations on the local scale enhance current knowledge of the spatial variability in hydroclimatic processes in highand mid-altitude mountain environments.


Introduction
The central part of the Hindu Kush Himalaya region presents tremendous heterogeneity, in particular in terms of topography and climatology. The terrain ranges from the agricultural plain of Terai to the highest peaks of the world, including Mount Everest, over a south-north transect about 150 km long ( Fig. 1).
Two main climatic processes on the synoptic scale are distinguished in the central Himalayas (Barros et al., 2000;Kansakar et al., 2004). First, the Indian monsoon is formed when moist air arriving from the Bay of Bengal is forced to rise and condense on the Himalayan barrier. Dhar and Rakhecha (1981) and Bookhagen and Burbank (2010) assessed that about 80 % of annual precipitation over the central Himalayas occurs between June and September. However, the timing and intensity of this summer monsoon is being reconsidered in the context of climate change (Bharati et al., 2016). The second main climatic process is a west flux that gets stuck in appropriately oriented valleys and occurs between January and March. Regarding high altitudes (> 3000 m), this winter precipitation can occur exclusively in solid form and can account for up to 40 % of annual precipitation (Lang and Barros, 2004) with considerable spatial and temporal variation.
Published by Copernicus Publications on behalf of the European Geosciences Union.
On a large spatio-temporal scale, precipitation patterns over the Himalayan Range are recognized to be strongly dependent on topography (Anders et al., 2006;Bookhagen and Burbank, 2006;Shrestha et al., 2012). The main thermodynamic process is an adiabatic expansion when air masses rise, but, at very high altitudes (> 4000 m), the reduction of available moisture is a concurrent process. Altitude thresholds of precipitation can then be discerned (Alpert, 1986;Roe, 2005). However, this representation of orographic precipitation has to be modulated considering the influence of such a protruding relief .
Products for precipitation estimation currently available in this area, e.g. the APHRODITE interpolation product (Yatagai et al., 2012) and the Tropical Rainfall Measuring Mission (TRMM) remote product (Bookhagen and Burbank, 2006), do not represent spatial and temporal variability in orographic effects at a resolution smaller than 10 km (Gonga-Saholiariliva et al., 2016). Consequently, substantial uncertainty remains in water budgets simulated for this region, as highlighted by Savéan et al. (2015). In this context, groundbased measurements condensed in small areas have been shown to enhance the characterization of local variability in orographic processes (Andermann et al., 2011;Pellicciotti et al., 2012;Immerzeel et al., 2014). However, although the Everest region is one of the most closely monitored areas of the Himalayan Range, valuable observations remain scarce. In particular, the relation between altitude and precipitation is still poorly documented.
The objective of this paper is to provide a representation of the effect of altitude on precipitation that represents spatial and temporal variability in precipitation in the Everest region. The parameters controlling the shape of the altitudinal factor are constrained through an original sensitivity analysis step. Uncertainties associated with variables simulated through the Interactions Soil-Biosphere-Atmosphere (ISBA) surface scheme (Noilhan and Planton, 1989) are quantified.
The first section of the paper presents the observation network and recorded data. The second section describes the model chosen to represent orographic precipitation, including computed altitude lapse rates for air temperature and precipitation. The method for statistical analysis through hydrological modelling is also described. The third section presents and discusses the results of sensitivity analysis and uncertainty analysis.
2 Data and associated uncertainties

Meteorological station transect
An observation network of 10 stations (Table 1 and Fig. 1) records hourly precipitation (P ) and air temperature (T ) since 2010 and 2014. The stations are equipped with classical rain gauges and HOBO ® sensors for temperature. The stations are located to depict the altitudinal profile of P and T over (1) the main river valley (Dudh Kosi Valley), oriented south-north, and (2) the Kharikhola tributary river, oriented east-west.
To reduce undercatching of solid precipitation, two Geonors ® were installed at 4218 and 5035 m in 2013. Measurements at Geonor ® instrumentation allow us to correct the effect of wind and the loss of snowflakes. Records from four other stations administrated by the Ev-K2-CNR association (www.evk2cnr.org) are also available. Total precipitation, air temperature, atmospheric pressure (AP), relative humidity (RH), wind speed (WS), short-wave radiation (direct and diffuse) (SW) and long-wave radiation (LW) have been recorded with an hourly time step since 2000 at Pyramid station (5035 m a.s.l.). Overall, these 10 stations cover an altitude range from 2078 m to 5035 m a.s.l., comprising a highly dense observation network, compared to the scarcity of ground-based data in this type of environment. The characteristics of the 10 stations are summarized Table 1. The meteorological data provided by the EVK2-CNR association are available by contacting the respective authors. All other data used in this study are freely available through the platform www.papredata.org.
Annual means of temperature and precipitation measured at these stations are presented are Table 2 for the two hydrological years 2014-2015 and 2015-2016. These time series contain missing data periods, which can represent up to 61 % of the recorded period. For stations LUK, NAM, PHA, PAN, PHE and PYR, where relatively long time series are available, gaps were filled with the interannual hourly mean for each variable. For the other stations, gaps were filled with values at the closest station, weighted by the ratio of mean values over the common periods. Time series from 1 January 2013 to 30 April 2016 were then reconstructed from these observations. Two seasons are defined based on these observations and knowledge of the climatology of the central Himalayas: (1) the monsoon season, from April to September, including the early monsoon, whose influence seems to be increasing with the current climate change (Bharati et al., 2014); (2) the winter season, dominated by westerly entrances with a substantial spatio-temporal variability.
Local measurements cannot be an exact quantification of any climatic variables, and they are necessarily associated with errors that follow an random distribution law. In particular, snowfall is usually undercaught by instrumentation (Sevruk et al., 2009). However, since this study focuses most particularly on uncertainty associated with spatialization of local measurements, aleatory errors in measurements will not be considered here.  (Table 3). Uncertainty in discharge is usually considered to account for less than 15 % of discharge (Lang et al., 2006).
Recession times are computed for available recession periods using the lfstat R library (Koffler and Laaha, 2013) with both the recession curves method (World Meteorological Organization, 2008) and the base flow index method (Chapman, 1999). We found recession times for Kharikhola and Tauche catchments of, respectively, around 70 days and around 67 days. Consequently, we consider that there is no interannual storage in either of the two catchments. This hypothesis can be modulated if a contribution of deep groundwater is considered (Andermann et al., 2011). Since these two catchments have no (Kharikhola) or neglectable (Tauche) glacier contribution, we hypothesized that the only entrance for water budgets in these catchments is total precipitation. In this study we used these two catchments as samples to assess generated precipitation fields against observed discharge on the local scale. The hydrological year is considered to start on 1 April, as decided by the Department of Hydrology and Meteorology of the Nepalese Government and generally accepted Savéan et al., 2015).
3 Spatialization methods for temperature and precipitation

Temperature
In mountainous areas, temperature and altitude generally correlate well linearly, considering a large timescale (Valéry et al., 2010;Gottardi et al., 2012). In the majority of studies based on field observations, air temperature values are extrapolated using the inverse distance weighting method (IDW) (Andermann et al., 2012;Immerzeel et al., 2012;Duethmann et al., 2013;Nepal et al., 2014). An altitude lapse rate θ (in • C km −1 ) is also used to take altitude into account for hourly temperature computation at any point M of the mesh extrapolated by IDW: where T is the hourly temperature, S i the ith station of the observation network, z i the altitude of station S i , z M the altitude of grid point M, and d −1 is the inverse of distance in latitude and longitude. In the Himalayas, seasonal Ragettli et al., 2015) or constant (Pokhrel et al., 2014) altitudinal lapse rates (LRs) are used for temperature. Figure 2 presents seasonal LRs computed from temperature time series at the 10 stations described in Sect. 2.1. The linearity is particularly satisfying for both seasons, even if stations follow differently oriented transects (W-E or N-S orientation). Computed LRs for both seasons are very close to values proposed by Immerzeel et al. (2014) and Heynen Consequently, these values for seasonal LRs will be used in this study. Uncertainties associated with temperature interpolation will therefore be neglected because they have minor impact on modelling compared to uncertainties in precipitation.

Model of orographic precipitation
The complexity of precipitation spatialization methods has been commented on by Barros and Lettenmaier (1993). When orographic effects are not well understood, complex approaches do not necessarily reproduce local measurements efficiently (Bénichou and Le Breton, 1987;Frei and Schär, 1998;Daly et al., 2002). In the central Himalayas, various hydrologic and glaciological studies are based on observation networks to produce a precipitation grid. However, few studies provide precipitation fields on an hourly timescale (Ragettli et al., 2015;Heynen et al., 2016), and precipitation fields on spatial scales lower than 1 km are always obtained using altitude linear lapse rates Nepal et al., 2014;Pokhrel et al., 2014). However, the considered lapse rates are constant in time and/or uniform in space. The spatial and temporal variability in the precipitation is then not represented in these studies. Moreover, the geostatistical co-kriging method has been applied by Gonga-Saholiariliva et al. (2016) for monsoon precipitation interpolation over the Kosi catchment. However, the provided precipitation fields overall underestimate the observations, and this method is shown not to be adequate for the interpolation of solid precipitation. The IDW method is a simple, widely Hydrol. Earth Syst. Sci., 21, 4879-4893, 2017 www.hydrol-earth-syst-sci.net/21/4879/2017/  (Valéry et al., 2010;Gottardi et al., 2012;Duethmann et al., 2013;Nepal et al., 2014). In the French Alps, Valéry et al. (2010) combine the IDW method with a multiplicative altitudinal factor. Precipitation at any point M of the mesh extrapolated by the IDW is given by In Eq.
(2), the altitude effect is represented through the introduction of the altitudinal factor β, defined by Valéry et al. (2010) as the slope of the linear regression between the altitude of stations (in m a.s.l.) and the logarithm of the seasonal volume of total precipitation expressed in millimetres. This method presents the advantage of using an altitudinal factor which can vary in time and space. The spatial and temporal variability in the precipitation is therefore represented in this method. Moreover, the effect of altitude is independently studied and the controlling parameters have physical meaning.

Observed relation between altitude and seasonal precipitation
Several studies based on observations (Dhar and Rakhecha, 1981;Barros et al., 2000;Bookhagen and Burbank, 2006;Immerzeel et al., 2014;Salerno et al., 2015) or theoretical approaches (Burns, 1953;Alpert, 1986) have observed that precipitation in the Himalayan Range generally presents a multimodal distribution along elevation. Precipitation is considered to increase with altitude until a first altitudinal threshold located between 1800 and 2500 m, depending on the study, and to decrease above 2500 m. Moreover, the linear correlation of precipitation with altitude is reported to be weak for measurements above 4000 m (Salerno et al., 2015). The decreasing of precipitation with altitude is characterized through various functions (Dhar and Rakhecha, 1981;Bookhagen and Burbank, 2006;Salerno et al., 2015). Nevertheless, the hypothesis of linearity of precipitation (P ) with altitude (z) is often made with a constant  or time-dependent lapse rate (Immerzeel et al., 2014). Gottardi et al. (2012) noted that, in mountainous areas, the hypothesis of a linear relation between P and z is only acceptable over a small spatial extension and for homogeneous weather types. Consequently, we considered altitude lapse rates for precipitation on the seasonal timescale, and we analysed the spatial variability in the relation between P and z.
For this purpose, we chose to regroup the stations into three groups (see Fig. 1): (1) stations with elevation ranging from 2078 to 3022 m, following a west-east transect (Group 1); (2) stations with elevation ranging from 2619 to 3570 m following a south-westerly transect (Group 2); and (3) stations with elevation above 3970 m (Group 3). Figure 3 shows that for Group 1, observed seasonal volumes of precipitation increase globally with altitude at a rate lower than 0.1 km −1 ; for Group 2, seasonal volumes decrease at a rate around −0.3 km −1 ; for Group 3, seasonal volumes decrease at a rate lower than 0.2 km −1 , with a poor linear trend.
The overlapping of altitude ranges between Group 1 and Group 2 highlights that the relation between precipitation and altitude strongly depends on terrain orientation. The difference in seasonal volumes at the BAL (2575 m a.s.l., 3471 mm year −1 ) and MER stations (2561 m a.s.l., 2245 mm year −1 ) (GROUP 1) also results from site effects on precipitation. In summary, β values inferred from local observations mainly express local variability and are not sufficient to establish any explicit relation between precipitation and altitude on the catchment scale. However, for operational purposes, the β factor can be simplified as a multimodal function of altitude within the Dudh Kosi catchment. The β factor is represented as a piecewise linear function of altitude using two altitude thresholds (z 1 and z 2 ) and three altitude lapse rates (β 1 , β 2 and β 3 ): As no deterministic value can be ensured for the five parameters controlling the shape of Eq. (3), an ensemble approach was applied (see Sect. 4) to estimate parameter sets on the scale of the entire Dudh Kosi River basin that are optimally suitable for both the Tauche and the Kharikhola catchments.
4 Sensitivity and uncertainties analysis method 4.1 Overall strategy Saltelli et al. (2006) distinguishes between sensitivity analysis (SA), which does not provide a measurement of error, and uncertainties analysis (UA), which computes a likelihood function according to reference data. SA is run before UA as a diagnostic tool, in particular to reduce variation intervals for parameters and therefore save computation time.
The algorithm chosen for SA was the regional sensitivity analysis (RSA) (Spear and Hornberger, 1980) method. The RSA method is based on the separation of the parameter space into (at least) two groups: behavioural or nonbehavioural parameter sets. A behavioural parameter set is a set that respects conditions (maximum or minimum thresholds) on the output of the orographic precipitation model. Thresholds will be defined for solid and total precipitation in the "Results and discussion" section. The analysis is performed using the R version of the SAFE(R) toolbox, developed by Pianosi et al. (2015). SA and UA are set up as follows (Beven, 2010): 1. First, the parameter space is sampled, according to a given sampling distribution. For each parameter set, hourly precipitation fields are computed at the 1 km resolution using Eq. (2) for both the Tauche and the Kharikhola catchments. Since physical processes strongly differ between the winter and monsoon seasons, we chose to differentiate the altitude correction for the two seasons. Behavioural parameter sets were then selected for each of the two seasons.
2. Then, for each behavioural precipitation field, the ISBA surface scheme, described in the next section, was run separately on Kharikhola and Tauche catchments. The objective function was computed as the difference between simulated and observed annual discharge at the outlet of each catchment. Parameter sets that lead to acceptable discharge regarding observed discharge for the two catchments are finally selected.

The ISBA surface scheme
We considered that there was no interannual storage in either of the two subcatchments studied; i.e. the variation in the groundwater content was considered zero from one hydrological year to the other. Consequently, annual simulated discharges were computed as the sum over all grid cells and all time steps of simulated surface flow and simulated subsurface flow. The question of the calibration of flow routing in the catchment was thus avoided. The ISBA surface scheme (Noilhan and Planton, 1989;Noilhan and Mahfouf, 1996) simulates interactions between the soil, vegetation and the atmosphere with a sub-hourly time step (SVAT model). The multilayer version of ISBA (ISBA-DIF) uses a diffusive approach (Boone et al., 2000;Decharme et al., 2011): surface and soil water fluxes are propagated from the surface through the soil column. Transport equations for mass and energy are solved using a multilayer vertical discretization of the soil. The explicit snow scheme in ISBA (ISBA-ES) uses a three-layer vertical discretization of snowpack and provides a mass and energy balance for each layer (Boone and Etchevers, 2001). Snowmelt and snow sublimation are taken into account in balance equations. The separation between runoff over saturated areas (Dunne runoff), infiltration excess runoff (Horton runoff) and infiltration is controlled by the variable infiltration capacity scheme (VIC) (Dümenil and Todini, 1992). The ISBA code is freely available from the respective authors.
The precipitation phase was estimated depending on hourly air temperature readings. Mixed phases occurred for temperatures between 0 and 2 • C, following a linear relation. Other input variables required for ISBA (atmospheric pressure, relative humidity, wind speed, short-and long-wave ra-  diations) were interpolated from measurements at Pyramid station as functions of altitude, using the method proposed by Cosgrove et al. (2003). Short-wave radiation and wind speed are not spatially interpolated and are considered to be equals to the measurements at Pyramid station for the two catchments.

Parametrization of surfaces
Several products provide parameter sets for physical properties of surfaces on the global scale (Hagemann, 2002;Masson et al., 2003;Arino et al., 2012). However, these products are not accurate enough at the resolution required for this study. The most recent analysis (Bharati et al., 2014;Ragettli et al., 2015) exclusively used knowledge garnered from the literature. To detail the approach, in this study the parametrization was based on in situ measurements. A classification into nine classes of soil/vegetation entities was defined based on Sentinel-2 images at a 10 m resolution (Drusch et al., 2012), using a supervised classification tool of the QGIS Semi-Automatic Classification Plugin (Congedo, 2015). In and around the two catchments, 24 reference sites were sampled during field missions. Data collection included soil texture, soil depth and root depth, determined by augering to a maximum depth of 1.2 m. Vegetation height and structure and dominant plant species were also determined. The results were classified into nine surface types. The nine classes and their respective fractions in the Kharikhola and Tauche catchments are presented Table 4.
Analysis of soil samples showed that soils were mostly sandy (∼ 70 %), with a small proportion of clay (∼ 1 %). Soil depths varied from very thin (∼ 30 cm) at high altitudes to 1.2 m for flat cultivated areas. Forest areas were separated into three classes: dry forests were characterized by high slopes and shallow soils; wet forests presented deep silty  (Congedo, 2015), based on Sentinel-2 images at a 10 m resolution (Drusch et al., 2012). In situ sample points were used to describe the soil and vegetation characteristics of each class. soils (1 m), with high trees (7 m); intermediate forests had moderate slopes and relatively deep, sandy soils. Crop areas presented different soil depths depending on their average slope. In addition, values for unmeasured variables (leaf area index (LAI), soil and vegetation albedos, surface emissivity, surface roughness) were taken from the ECO-CLIMAP1 classification (Masson et al., 2003) for ecosystems representative of the study area. ECOCLIMAP1 provides the annual cycle of dynamic vegetation variables, based both on a surface properties classification (Hagemann, 2002) and on a global climate map (Koeppe and De Long, 1958). The ECOCLIMAP2 product (Faroux et al., 2013) is derived from ECOCLIMAP1 and provides enhanced descriptions of surfaces. However, ECOCLIMAP2 is only available for Europe and therefore is not used in this study. Table 4. Soil and vegetation characteristics of the nine classes defined in the Kharikhola and Tauche catchments; "% KK" and "% Tauche" represent the fraction of each class in the Kharikhola and Tauche catchments. Sand and clay fractions ("% Sand" and "% Clay", respectively), soil depth (SD), root depth (RD), and tree height (TH) are defined based on in situ measurements. The dynamic variables (e.g. the fraction of vegetation and leaf area index) are found in the ECOCLIMAP1 classification (Masson et al., 2003)  5 Results and discussion 5.1 Regional sensitivity analysis The parameter space was sampled using the "All at a time" (AAT) sampling algorithm from the SAFE(R) toolbox (Pianosi et al., 2015). Since no particular information was available on prior distribution and interaction for the five parameters, uniform distributions were considered. The size of parameter samples was chosen according to Sarrazin et al. (2016) (Table 5). The optimization method is highly sensitive to the choice of initial values for the β 1 , β 2 , β 3 , z 1 and z 2 parameters. Several attempts have been made, and the choices presented Table 6 are justified by the following arguments: -Minimum and maximum values for the altitude thresholds z 1 and z 2 are chosen according to both literature review (Barros et al., 2000;Anders et al., 2006;Bookhagen and Burbank, 2006;Shrestha et al., 2012;Nepal, 2012;Savéan, 2014) and observations. The first inquired altitudinal threshold is described in the literature as between 2000 and 3000 m, and the second threshold is described as above 4000 m. These intervals have been enlarged to also test related values.
-Maximum (minimum) values for β 1 (β 2 ) are chosen about 10 times larger than the value computed based on observation. Considering the definition of the beta coefficient, a value greater than 2 km −1 (lower than −2 km −1 ) would lead to a multiplication of precipitation by 1.22 (by 0.82) within 100 m. When applied to the precipitation observed at stations, this would lead to inconsistent precipitation when increasing altitude by 100 m.
-The β 3 coefficient has to be negative because a positive value would lead to unrealistic values at high altitudes. Moreover, the minimum value is chosen to be significantly smaller than the value computed for β 3 based on the observations but also to remain higher than the value computed for β 2 based on the observations.
A behavioural parameter set is a set that leads to an annual amount of total precipitation for both catchments comprised between a minimum and a maximum value. A parameter sets that does not meet these conditions is considered as nonbehavioural. Maximum and minimum conditions on annual total precipitation for a set to be behavioural were chosen according to the annual observed discharge for each of the two catchments. The mean observed discharge for the recorded period was 2043 mm year −1 at the Kharikhola station and 457 mm year −1 at the Tauche station. Annual total precipitation was expected to be greater than the measured annual discharge and lower than annual discharge plus 70 %. These thresholds take into account both the uncertainty in measured discharges and actual evapotranspiration. Based on values proposed in the literature, evapotranspiration is assumed to represent less than 50 % of observed discharge, for both catchments. The minimum and maximum thresholds for both catchments are summarized in Table 7.
The method's convergence (i.e. the stability of the result when the sample size grows) was graphically assessed. The results converged for sample sizes from 1000 samples. Figure 5 shows the cumulative density function (CDF) for behavioural and nonbehavioural parameter sets for the monsoon and winter seasons. Of the 2000 parameter sets sampled, 712 sets verified the chosen minimum and maximum conditions for annual total precipitation and snowfall (i.e. they were behavioural). The sensitivity of the output to each parameter was evaluated by the maximum vertical distance (MVD) between CDF for both behavioural and nonbehavioural parameter sets. Annual total precipitation appeared to be less sensitive to parameters controlling winter precipitation than to parameters controlling monsoon precipitation. This result can be explained by the fact that winter precipitation was less than monsoon precipitation. However, since the applied sampling method does not take into account the Hydrol. Earth Syst. Sci., 21, 4879-4893, 2017 www.hydrol-earth-syst-sci.net/21/4879/2017/  Table 6. Initial ranges considered for the five shape parameters of the altitudinal factor: z 1 , z 2 β 1 , β 2 and β 3 . Ranges are defined based on measurements at stations and on values found in the literature. existing interaction between the five parameters, further analysis for parameter ranking was not significant. The method was necessarily sensitive to the prior hypothesis presented Table 5. In particular, the conditions for a set to be behavioural have a significant impact on the distribution of the behavioural sets. In contrast, increasing the sample size does not affect the output distribution, since minimum size for convergence is reached.

Annual simulated water budgets
The precipitation fields generated using each behavioural parameter set were used as input data within the ISBA surface scheme. The simulations over the Tauche and Kharikhola catchments were run separately over the 1 January 2013 to 31 March 2016 period, on an hourly timescale. The 2013-2014 hydrological year was used as a spin-up period, and the results were observed for the 2014-2015 and 2015-2016 hydrological years. To overcome the issue of calibrating a flowrouting module, the simulated discharge were aggregated on an annual timescale and compared to annual observed discharge at the outlet (Q obs ). Figure 6 presents boxplots obtained for the 712 behavioural parameter sets for the terms of the annual water budget, i.e. liquid and solid precipitation, discharge, and evapotranspiration. The dashed line represents Q obs for each catchment. The mean annual volumes of simulated variables were also computed for each parameter set in 2014-2015 and 2015-2016, and the intervals of uncertainty associated with simulated annual volumes are provided. This method highlights the propagation of uncertainties associated with Table 7. Maximum and minimum condition on total precipitation for a parameter set to be behavioural, for the Kharikhola and Tauche catchments. Annual total precipitation was expected to be greater than the measured annual discharge plus 20 % and lower than annual discharge plus 50 %.

Minimum Maximum
Kharikhola 2043 3473 mm year −1 Tauche 457 777 mm year −1 the representation of orographic effects toward the simulated terms of annual water budgets. Table 8 presents the mean value, standard deviation and relative standard deviation for all of the ISBA-simulated variables for the Kharikhola and Tauche catchments for 2014-2015 and 2015-2016. The annual actual evapotranspiration accounted for 26 % of annual total precipitation for Kharikhola and 34 % for Tauche. In comparison, evapotranspiration was estimated at about 20, 14 and 53 % of total annual precipitation, respectively, by Andermann et al. (2012 and Savéan et al. (2015) over the entire Dudh Kosi Basin, and Ragettli et al. (2015) estimated it at 36.2 % of annual total precipitation for the upper part of the Langtang Basin.
Annual snowfall volume for Kharikhola was a neglectable fraction of annual total precipitation (∼ 1 %), and it was around 44 % for Tauche. Annual snowfall was estimated at, respectively, 15.6 and 51.4 % of annual total precipitation by Savéan et al. (2015) (entire Dudh Kosi River basin) and Ragettli et al. (2015) (upper part of the Langtang Basin).
Moreover, this statistical approach shows that the only uncertainties associated with representation of the orographic effect result in significant uncertainties in simulated variables. These uncertainties account for up to 16 % for annual total precipitation, up to 25 % for annual discharge and up to 8 % for annual actual evapotranspiration. Uncertainty in annual snowfall is quantified at 16 % for a high-mountain catchment and up to 32 % for a middle-mountain catchment. These uncertainty intervals are essentially conditioned by model structure and parametrization, and these results indicate that simulated water budgets provided by modelling studies must necessarily be associated with error intervals.

Toward optimizing parameter sets with bias in annual discharge
Going further into the simulation results, the hydrological cycle was inverted, in order to use observed discharge to optimize the relation between precipitation and altitude, as presented for mountainous areas by Valéry et al. (2009). Precipitation fields were then constrained on the local scale according to simulated discharges. Annual bias in discharge was computed for each catchment as the absolute value of the ratio between the observed and simulated annual dis-   Figure 7 presents the scatter plot of the distributions of bias in annual discharge for the Kharikhola and Tauche catchments. The Pareto optima, minimizing bias in annual discharge for both catchments, were computed using the R rPref package (Roocks and Roocks, 2016). For example, the first 10 Pareto optima were selected among the 712 behavioural parameter sets considered. The values of parameters for the winter and monsoon seasons for the 10 first optimum sets are summarized in Table 9. For the 10 parameter sets selected, the altitudinal threshold z 1 was located between 2010 and 3470 m a.s.l. during the monsoon season and between 2287 m a.s.l. and 3488 m a.s.l. during winter. The second altitudinal threshold z 2 was located between 3709 and 6167 m a.s.l. during monsoon and between 3734 and 6466 m a.s.l. during winter. Altitudes found for z 1 were globally higher than altitudes proposed in the literature for the second mode of precipitation (between 1800 and 2400 m a.s.l., as described in Sect. 3.2.2). Since these values were calibrated on the local scale, according to groundbased measurements, they can be considered to accurately represent the local variability encountered in the Tauche and Kharikhola catchments. Moreover, values for an altitudinal threshold of precipitation located above 4000 m a.s.l. were proposed.

Ensemble of hourly precipitation fields on the Dudh Kosi River basin
Observed precipitation at measuring stations was then interpolated on an hourly timescale over the Dudh Kosi River   basin at the 1 km spatial resolution. The method given by Eq.
(2) is applied, using shape parameters for the altitudinal factor selected Table 9. The average annual volumes of computed total precipitation ranged between 1365 and 1652 mm, and annual snowfall volumes ranged between 89 and 126 mm, on average over the 2014-2015 and 2015-2016 hydrological years. These values are consistent with other products available for the area. In particular, Savéan (2014) showed that the APHRODITE (Yatagai et al., 2012)   period. Different relations between altitude and annual precipitation are then represented. The steeper the slope, the more variation in precipitation. This has to be considered in light of the physical properties of convection at such high altitudes.

Conclusions
The main objective of this paper is to provide a representation of the effect of altitude on precipitation that represents spatial and temporal variability in precipitation in the Everest region. A weighted inverse distance method coupled with a multiplicative altitudinal factor was applied to spatially extrapolate measured precipitation to produce precipitation fields over the Dudh Kosi Basin. The altitudinal factor for the Dudh Kosi Basin is shown to acceptably fit a piecewise linear function of altitude, with significant seasonal variations. A sensitivity analysis was run to reduce the variation interval for parameters controlling the shape of the altitudinal factor. An uncertainty analysis was subsequently run to evaluate an ensemble of simulated variables according to observed discharge for two small subcatchments of the Dudh Kosi Basin located in mid-and high-altitude mountain environments. Non-deterministic annual water budgets are provided for two small gauged subcatchments located in high-and midaltitude mountain environments. This work shows that the only uncertainties associated with representation of the orographic effect account for about 16 % for annual total precipitation and up to 25 % for simulated discharges. Annual evapotranspiration is shown to represent 26 % ± 1 % of annual total precipitation for the mid-altitude catchment and 34% ± 3 % for the high-altitude catchment. Snowfall contribution is shown to be neglectable for the mid-altitude catchment, and it represents up to 44 % ± 8 % of total precipitation for the high-altitude catchment. These simulations on the lo-cal scale enhance current knowledge of the spatial variability in hydroclimatic processes in high-and mid-altitude mountain environments.
This work paves the way to produce hourly precipitation maps extrapolated from ground-based measurements that are reliable on the local scale. However, additional criteria would be needed to provide a single optimum parameter set for the altitudinal factor that would be suitable for the entire Dudh Kosi River basin. For example, snow cover areas simulated on a scale larger than the two catchments could be compared to available remote products (Behrangi et al., 2016). Independent measurements of precipitation could also be used to constrain the ensemble of precipitation fields.
Moreover, since observations are made over a very short duration and contain long periods with missing information, the results are limited to the 2014-2015 and 2015-2016 hydrological years and to the Dudh Kosi River basin. In addition, this study focuses only on one source of uncertainty in the measurement-spatialization-modelling chain, whereas sensitivity analysis should include all types of uncertainty (Beven, 2015;Saltelli et al., 2006). A more complete method would include epistemic uncertainty in model parameters and aleatory uncertainty in input variables in the sensitivity analysis (Fuentes Andino et al., 2016).
Data availability. The ISBA model implemented within the Surfex platform is freely available at www.umr-cnrm.fr/surfex/. The SAFE(R) toolbox is freely available at www.safetoolbox.info/. The data provided by the EV-K2-CNR association are freely available at www.evk2cnr.org. For collaboration reasons, the observation data provided by the PRESHINE project will be freely available in December 2018. These data will be distributed through www.papredata.org. The codes used for this work are implemented in R language and they are freely available at www.papredata.org/ codes.asp.