Journal topic
Hydrol. Earth Syst. Sci., 24, 3289–3309, 2020
https://doi.org/10.5194/hess-24-3289-2020
Hydrol. Earth Syst. Sci., 24, 3289–3309, 2020
https://doi.org/10.5194/hess-24-3289-2020

Research article 29 Jun 2020

Research article | 29 Jun 2020

# Comparing Bayesian and traditional end-member mixing approaches for hydrograph separation in a glacierized basin

Comparing Bayesian and traditional end-member mixing approaches for hydrograph separation in a glacierized basin
Zhihua He1,a, Katy Unger-Shayesteh3, Sergiy Vorogushyn1, Stephan M. Weise4, Doris Duethmann5,7, Olga Kalashnikova6, Abror Gafurov1, and Bruno Merz1,2 Zhihua He et al.
• 1Hydrology Section, GFZ German Research Centre for Geosciences, Telegrafenberg, Potsdam, Germany
• 2Institute for Environmental Sciences and Geography, University of Potsdam, Potsdam, Germany
• 3International Relations, German Aerospace Center (DLR), Linder Höhe, Cologne, Germany
• 4Department of Catchment Hydrology, Helmholtz Centre for Environmental Research (UFZ), Halle, Germany
• 5Institute of Hydraulic Engineering and Water Resources Management, Vienna University of Technology (TU Wien), Vienna, Austria
• 6Climate, Water and Natural Resources Department, Central Asian Institute for Applied Geosciences (CAIAG), Bishkek, Kyrgyzstan
• 7Department of Ecohydrology, Leibniz-Institute of Freshwater Ecology and Inland Fisheries (IGB), Berlin, Germany

Correspondence: Zhihua He (zhwork3533@163.com)

Abstract

Tracer data have been successfully used for hydrograph separation in glacierized basins. However, in these basins uncertainties of the hydrograph separation are large and are caused by the spatiotemporal variability in the tracer signatures of water sources, the uncertainty of water sampling, and the mixing model uncertainty. In this study, we used electrical conductivity (EC) measurements and two isotope signatures (δ18O and δ2H) to label the runoff components, including groundwater, snow and glacier meltwater, and rainfall, in a Central Asian glacierized basin. The contributions of runoff components (CRCs) to the total runoff and the corresponding uncertainty were quantified by two mixing approaches, namely a traditional end-member mixing approach (abbreviated as EMMA) and a Bayesian end-member mixing approach. The performance of the two mixing approaches was compared in three seasons that are distinguished as the cold season, snowmelt season, and glacier melt season. The results show the following points. (1) The Bayesian approach generally estimated smaller uncertainty ranges for the CRC when compared to the EMMA. (2) The Bayesian approach tended to be less sensitive to the sampling uncertainties of meltwater than the EMMA. (3) Ignoring the model uncertainty caused by the isotope fractionation likely led to an overestimated rainfall contribution and an underestimated meltwater share in the melt seasons. Our study provides the first comparison of the two end-member mixing approaches for hydrograph separation in glacierized basins and gives insight into the application of tracer-based mixing approaches in similar basins.

1 Introduction

Glaciers and snowpack store a large amount of fresh water in glacierized basins, thus providing an important water source for downstream human societies and ecosystems (Barnett et al., 2005; Viviroli et al., 2007; He et al., 2014; Penna et al., 2016). Seasonal meltwater and rainfall play significant roles in shaping the magnitude and timing of runoff in these basins (Rahman et al., 2015; Pohl et al., 2017). Quantifying the seasonal contributions of the runoff components (CRCs), including groundwater, snowmelt, glacier melt, and rainfall, to the total runoff is therefore highly necessary for understanding the dynamics of water resources in glacierized basins under the current climate warming (La Frenierre and Mark, 2014; Penna et al., 2014; He et al., 2015).

The traditional end-member mixing approach (abbreviated as EMMA) has been widely used for hydrograph separation in glacierized basins across the world (Dahlke et al., 2014; Sun et al., 2016a; Pu et al., 2017). For instance, studies in the glacierized catchments of the Italian Alps indicate the successful application of the EMMA to estimate the proportions of groundwater, snow, and glacier meltwater based on water stable isotopes and electric conductivity (EC; e.g., Chiogna et al. 2014, Engel et al. 2016, and Penna et al. 2017). Using EMMA, Li et al. (2014) confirmed significant contributions of snow and glacier melt runoff to total runoff in the Qilian Mountains. Maurya et al. (2011) reported the contribution of glacial ice meltwater to the total runoff in a Himalayan basin on δ18O and EC using a three-component EMMA.

However, uncertainties of CRCs quantified by EMMA in glacierized basins are typically high (Klaus and McDonnell, 2013; Rahman et al., 2015) because of the following reasons. (1) The catchment elevation generally extends over a large range, leading to strong spatial variability in climate forcing (precipitation and temperature) and the tracer signatures of water sources. (2) The number of end-member water sources for runoff is typically high, including snow and glacier meltwater. (3) Water sampling in a high-elevation glacierized catchment is difficult due to logistical limitations that result in small sample sizes for the application of EMMA. The uncertainties of CRCs can be categorized into statistical uncertainty and model uncertainty. Statistical uncertainty refers to the spatiotemporal variability of the tracer signatures, sampling uncertainty, and laboratory measurement error (Joerin et al., 2002). Model uncertainty is determined by the assumptions of the EMMA, which might not agree with the reality in the basin (Joerin et al., 2002; Klaus and McDonnell, 2013). For example, the fractionation effect on isotope ratios caused by evaporation during the mixing process can result in significant errors given the constant tracer assumption in the EMMA (Moore and Semmens, 2008).

The Gaussian error propagation technique has been typically applied along with EMMA to estimate the statistical uncertainty for hydrograph separation, assuming that the uncertainty associated with each source is independent of the uncertainty of other sources (Genereux, 1998; Pu et al., 2013). The spatiotemporal variability of the tracer signatures is estimated by multiplying the t values of the Student's t distribution at the selected significance level with the standard deviations (SDs) of the measured tracer signatures (Pu et al., 2013; Penna et al., 2016; Sun et al., 2016b). Although this approach has been used successfully in various glacierized basins, some recurring issues remain. (1) These include inappropriate estimations of the variability of tracer signatures of water sources when only a few water samples are available (Dahlke et al., 2014). The SD values of the measured tracer signatures likely fail to represent the variability of the tracer signatures of individual water sources across the basin due to the small water sample sizes. (2) The correlations of tracer signatures and the dependence of runoff components are inevitably ignored due to the assumption of the independence of the multiple uncertainty sources. The correlation between δ18O and δ2H of each water source, and the interaction between runoff components, could provide additional constraints on the uncertainty of the quantification of runoff components that are, however, typically ignored in the Gaussian error propagation technique. Furthermore, the model uncertainty caused by the fractionation effect on isotope ratios during the mixing process is also often ignored.

The Bayesian end-member mixing approach (shortened to Bayesian approach) shows the potential of estimating the proportions of individual components to the mixing variable in a more rigorous, statistical way (Parnell et al., 2010). For hydrograph separation, the tracer signatures of the water sources are first assumed to obey specific prior distributions. Their posterior distributions are then obtained by updating the prior distributions with the likelihood derived from water samples. In the last step, the CRCs to the total runoff are estimated based on the balance of the posterior tracer signatures. The posterior distributions of the CRCs are typically estimated in a Markov chain Monte Carlo (MCMC) procedure. In the Bayesian approach, both the statistical and model uncertainties are represented by the posterior distributions of parameters. The parameter uncertainty is estimated based on likelihood observations using MCMC.

Although the Bayesian approach can be applied in cases when the sample sizes are small (Ward et al., 2010), it has rarely been used for hydrograph separation in glacierized basins. To the authors' knowledge, there have been only four studies, including Brown et al. (2006), that conducted the hydrograph separation using a three-component Bayesian approach in a glacierized basin in the French Pyrenees. Furthermore, Cable et al. (2011) quantified the CRCs to total runoff in a glacierized basin in the North American Rocky Mountains. They used a hierarchical Bayesian framework to incorporate the temporal and spatial variability of the water isotope data into the mixing model. Rodriguez et al. (2016) investigated the effects of tracer measurements and mixing model parameters on the quantification of CRCs in a Chilean glacierized basin using an informative Bayesian framework. Recently, Beria et al. (2020) used a classic Bayesian approach to estimate the uncertainty of CRCs at a Swiss Alpine catchment. However, the performance of the Bayesian approach has not been evaluated in comparison to the EMMA. Moreover, the sensitivity of the Bayesian approach to the water sampling uncertainty associated with the representativeness of the water samples caused by the limited sample site and sample size is still not clear. Benefiting from the prior assumptions of changes in isotope signatures during the mixing process, the Bayesian approach bears the potential to estimate the fractionation effect on isotopic signatures (Moore and Semmens, 2008), which has, however, not been investigated either.

In this study, we compare the EMMA and the Bayesian approach for hydrograph separation in a Central Asian glacierized basin using water isotope and EC measurements. In Central Asia, glacierized catchments provide an important fresh water supply for downstream cities and irrigated agriculture. Quantifying the contributions of multiple runoff components to total runoff is important for understanding the dynamics of water resource availability at the regional scale. However, the uncertainty of the quantification of runoff components in the glacierized catchments is particularly large, as mentioned before. Our research question is twofold. First, how do the EMMA and Bayesian approaches compare with respect to the quantification of CRCs? Second, what are the influences of the different uncertainty sources (including variability of the tracer signatures, sampling uncertainty, and model uncertainty) on the estimated CRCs in the two mixing approaches?

The paper is organized as follows: details on the study basin and water sampling are introduced in Sect. 2; assumptions of the two mixing approaches are described in Sect. 3; Sect. 4 estimates the CRCs and the corresponding uncertainties; and the discussion and conclusion finalize the paper in Sects. 5 and 6, respectively.

2 Study area and data

## 2.1 Study area

Located in Kyrgyzstan, Central Asia, the Ala-Archa basin drains an area of 233 km2 (Fig. 1), and glaciers cover around 17 % of the basin area. The elevation of the study basin extends from 1560 to 4864 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, and the elevation range of the glacierized area extends from 3218 to 4857 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, with about 76 % located between 3700 and 4100 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ The Golubin glacier has an area of  5.7 km2 and extends over an elevation range from 3232 to 4458 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 1). Both the elevation range and the mean elevation (3869 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) of the Golubin glacier are close to those of the entire glacierized area (mean elevation is 3945 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). The Golubin glacier represents about 14.4 % of the entire glacierized area, while its elevation range covers around 95.6 % of the entire glacier range. The annual mean precipitation and air temperature measurements taken at the Baitik meteorological station in Kyrgyzstan during 2012–2017 were 538 mm yr−1 and 7.2 C, respectively. The mean daily streamflow during 2012–2017 was about 6.3 m3 s−1 (Fig. S1 in the Supplement). The seasonal dynamics of runoff in the river play an important role in the water availability for downstream agricultural irrigation. The generation of snow and glacier melt runoff generally shows the largest effect on the runoff seasonality (Aizen et al., 2000, 2007). In particular, the snowmelt runoff mainly occurs in the warm period from early March to mid-September, and the glacier melt (referring to ice melt in our study basin) typically generates runoff from the high-elevation areas from July to September (Aizen et al., 1996; He et al., 2018, 2019). We subsequently defined three runoff-generation seasons as follows. The cold season from October to February, in which the streamflow is fed mainly by groundwater and, to a smaller extent, by snowmelt and rainfall; snowmelt season from March to June, in which the streamflow is chiefly fed by snowmelt and groundwater and, additionally, by rainfall; and glacier melt season from July to September, in which the streamflow is fed by significant glacier melt and groundwater together with rainfall and snowmelt.

Figure 1Study area of the Ala-Archa basin (derived from the World Topographic Map by © ESRI) and the Golubin glacier, including the locations of the water sampling points.

Two meteorological stations (Fig. 1), i.e., Alplager (at an elevation of 2100 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) and Baitik (at an elevation of 1580 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), were set up in the basin in the 1960s to collect precipitation and temperature data. The Ala-Archa hydrological station was set up at the same site as the Baitik meteorological station to collect the daily average streamflow data from the 1960s onwards. The dynamics of glacier mass balance and snowpack in the accumulation zone were surveyed in summer field campaigns throughout 2012–2017. The daily precipitation, temperature, and streamflow measured at the basin outlet during 2012–2017 are presented in Fig. S1.

## 2.2 Tracer data

Since July 2013, local station operators have collected weekly stream water samples from the river channel close to the Alplager and Baitik meteorological sites using pure 50 mL high-density polyethylene (HDPE) bottles (He et al., 2019). The slightly varied sampling time was at around noon every Wednesday. Precipitation samples were collected during 2012–2017 at four sites across the basin (Fig. 1). At the Alplager and Baitik meteorological sites, the precipitation samples were first collected from fixed rain collectors (immediately after the rainfall/snowfall events) and were then accumulated in two indoor rain containers over one month. The mixed water in the containers were then sampled for isotopic analysis every month. The indoor rain containers were filled with thin mineral oil layers for monthly precipitation accumulation and stored in cold places. Additionally, two plastic rain collectors (PALMEX, as in Gröning et al., 2012), specifically designed for isotopic sampling and to prevent evaporation, were set up at elevations of 2580 and 3300 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ to collect precipitation in high-elevation areas (Fig. 1). Precipitation samples were collected monthly from these two rain collectors during the period from May to October when the high-elevation areas were accessible.

Glacier meltwater was sampled during the summer field campaigns in each year from 2012 to 2017. Samples of meltwater flowing on the Golubin glacier in the ablation zone and at the glacier tongue were collected by pure 50 mL HDPE bottles and then stored in a cooling box (Fig. 1; the elevation of the sampling sites ranges from 3280 to 3805 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). We only collected glacier meltwater samples from the Golubin glacier due to the logistic limitations of the remaining glacierized area. Snow samples were collected from early March to early October during 2012–2017, as the sampling sites are generally inaccessible due to the heavy snow accumulation in the remaining months. The elevation of the multiple snow sampling sites ranges from 1580 to 4050 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 1). The whole snow profile at each sampling site was collected by drilling a 1.2 m pure plastic tube into the snowpack. The snow in the whole tube were then collected in plastic bags and stored in a cooling box. After all the snow in the plastic bags melted out, the mixed snow meltwater samples were then collected with pure HDPE bottles. Groundwater samples were collected from a spring draining to the river using pure HDPE bottles from March to October during 2012–2017 (Fig. 1; 2400 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). The spring is located at the foot of a rocky hill, around 60 m away from the river channel.

All samples were stored at 4 C and then delivered to the laboratory at the Helmholtz Centre for Environmental Research (UFZ) in Halle, Germany, by air. Isotopic compositions of water samples were measured using laser-based infrared spectrometry (Triple Isotope Water Analyzer (LGR TIWA 45), Los Gatos Research, Inc.; Picarro L1102-i, Picarro, Inc.). A correction procedure has been carried out to minimize the effects of drifts and sample-to-sample memory following the Laboratory Information Management System (LIMS) for lasers 2015 developed by Coplen and Wassenaar (2015). The measurement precisions of both the LGR TIWA 45 and Picarro L1102-i for δ18O and δ2H are ± 0.25 and ± 0.4 ‰, respectively, after the calibration against the common standard of the Vienna Standard Mean Ocean Water (VSMOW). We used the HI9813 portable meter (with pH, electrical conductivity (EC), and total dissolved solids (TDS) measurement functions; Hanna Instruments) to measure the EC values of water samples with a measurement precision of 0.1 µs cm−1. The EC data have been widely used for hydrograph separation due to their easy use and quick measurement. While the EC of water sources is not a conservative tracer when transporting along the subsurface path, in our case this may only have a small effect on the application of hydrograph separation. The measured EC values of water sources (rainfall, snow, and glacier melt) primarily label the surface direct runoff that has weak interaction with mineral soil. The EC indicator measured from the spring water is assumed to be the mean EC value of the groundwater contributing to the streamflow because the elevation of the sampled spring is close to the mean elevation of the basin, and the areas of the regions above and below the spring are very close. Abnormal isotopic compositions caused by evaporation and abnormal EC values caused by impurities were discarded. We used threshold values to identify abnormal values of δ18O and EC, which are defined as values located more than 5 % away from the sample clusters. For δ18O, sample values higher than 5 ‰ were excluded. For EC, sample values higher than 210 µs cm−1 were excluded. Tracer data of individual water sources at the sampled date are presented in Fig. S1.

3 Methodology

The hydrograph separation is carried out in each of the three seasons (i.e., cold season, snowmelt season, and glacier melt season). Water samples collected in the period from 2012 to 2017 are split into each of the three seasons for the hydrograph separation. The CRCs estimated by the mixing approaches refers to the mean contributions in each of the three seasons during the period from 2012 to 2017. The mixing approaches applied for the hydrograph separation in each season are summarized in Table 2. Considering that the groundwater and snowmelt samples were rarely collected in the cold season, we used all available groundwater and snowmelt samples from the three seasons for hydrograph separation in the cold season. Tracer signatures of rainfall are assumed to be the same as the measured tracer signatures of precipitation samples in all three seasons.

## 3.1 Traditional end-member mixing approach (EMMA)

The main assumptions of EMMA include the following points (Kong and Pang, 2012). (1) The tracer signature of each runoff component is constant during the analyzed period. (2) The tracer signatures of the runoff components are significantly different to each other. (3) Tracer signatures are conservative in the mixing process. In the cold and snowmelt seasons, a three-component EMMA method (EMMA_3; Table 2) is used. Because the precision of δ18O (± 0.25 ‰) measured in the lab is higher than that of δ2H (± 0.4 ‰), and both are strongly correlated, the EMMA_3 is based on δ18O and EC. In the glacier melt season, both the EMMA_3 and the four-component EMMA (EMMA_4; Table 2) are used. In the EMMA_3, glacier melt and snowmelt are assumed to be one end member because of their similar tracer signatures. In the EMMA_4, glacier melt and snowmelt are treated separately as two end members, and δ18O and δ2H are used as two separate tracers. The following equations (Eqs. 1–5) are used to estimate CRCs (f1−3) and the corresponding uncertainty in the EMMA_3 (Genereux, 1998).

$\begin{array}{}\text{(1)}& & \left\{\begin{array}{ll}\mathrm{1}={f}_{\mathrm{1}}+{f}_{\mathrm{2}}+{f}_{\mathrm{3}},& \text{for water balance}\\ A={A}_{\mathrm{1}}\cdot {f}_{\mathrm{1}}+{A}_{\mathrm{2}}\cdot {f}_{\mathrm{2}}+{A}_{\mathrm{3}}\cdot {f}_{\mathrm{3}},& \text{for water tracer A}\\ B={B}_{\mathrm{1}}\cdot {f}_{\mathrm{1}}+{B}_{\mathrm{2}}\cdot {f}_{\mathrm{2}}+{B}_{\mathrm{3}}\cdot {f}_{\mathrm{3}},& \text{for water tracer B}\end{array}\right\\text{(2)}& & {f}_{\mathrm{1}}=\frac{A{B}_{\mathrm{2}}-A{B}_{\mathrm{3}}+{A}_{\mathrm{2}}{B}_{\mathrm{3}}-{A}_{\mathrm{2}}B+A{}_{\mathrm{3}}B-{A}_{\mathrm{3}}{B}_{\mathrm{2}}}{{A}_{\mathrm{1}}{B}_{\mathrm{2}}-{A}_{\mathrm{1}}{B}_{\mathrm{3}}+{A}_{\mathrm{2}}{B}_{\mathrm{3}}-{A}_{\mathrm{2}}{B}_{\mathrm{1}}+{A}_{\mathrm{3}}{B}_{\mathrm{1}}-{A}_{\mathrm{3}}{B}_{\mathrm{2}}}\text{(3)}& & {f}_{\mathrm{2}}=\frac{A{B}_{\mathrm{3}}-A{B}_{\mathrm{1}}+{A}_{\mathrm{1}}B-{A}_{\mathrm{1}}{B}_{\mathrm{3}}+A{}_{\mathrm{3}}{B}_{\mathrm{1}}-{A}_{\mathrm{3}}B}{{A}_{\mathrm{1}}{B}_{\mathrm{2}}-{A}_{\mathrm{1}}{B}_{\mathrm{3}}+{A}_{\mathrm{2}}{B}_{\mathrm{3}}-{A}_{\mathrm{2}}{B}_{\mathrm{1}}+{A}_{\mathrm{3}}{B}_{\mathrm{1}}-{A}_{\mathrm{3}}{B}_{\mathrm{2}}}\text{(4)}& & {f}_{\mathrm{3}}=\frac{A{B}_{\mathrm{1}}-A{B}_{\mathrm{2}}+{A}_{\mathrm{1}}{B}_{\mathrm{2}}-{A}_{\mathrm{1}}B+A{}_{\mathrm{2}}B-{A}_{\mathrm{2}}{B}_{\mathrm{1}}}{{A}_{\mathrm{1}}{B}_{\mathrm{2}}-{A}_{\mathrm{1}}{B}_{\mathrm{3}}+{A}_{\mathrm{2}}{B}_{\mathrm{3}}-{A}_{\mathrm{2}}{B}_{\mathrm{1}}+{A}_{\mathrm{3}}{B}_{\mathrm{1}}-{A}_{\mathrm{3}}{B}_{\mathrm{2}}},\end{array}$

where the subscripts 1–3 refer to the three runoff components (i.e., groundwater, snowmelt/meltwater, and rainfall), and A1A3 (B1B3) refer to the mean δ18O (EC) values of runoff components. A and B stand for the mean δ18O and EC values of the stream water. The mean isotope and EC values of precipitation are calculated as the monthly precipitation weighted average values. Similarly, the mean isotope and EC values of stream water are calculated as the weekly streamflow weighted average values.

Assuming the uncertainty of each variable is independent of the uncertainty in others, the Gaussian error propagation technique is applied to estimate the uncertainty of the CRCs (f1−3) using the following equation (Genereux, 1998):

$\begin{array}{}\text{(5)}& {W}_{{f}_{i}}=\sqrt{\begin{array}{l}{\left(\frac{\partial {f}_{i}}{\partial {A}_{\mathrm{1}}}{W}_{{A}_{\mathrm{1}}}\right)}^{\mathrm{2}}+{\left(\frac{\partial {f}_{i}}{\partial {A}_{\mathrm{2}}}{W}_{{A}_{\mathrm{2}}}\right)}^{\mathrm{2}}+{\left(\frac{\partial {f}_{i}}{\partial {A}_{\mathrm{3}}}{W}_{{A}_{\mathrm{3}}}\right)}^{\mathrm{2}}+\\ {\left(\frac{\partial {f}_{i}}{\partial A}{W}_{A}\right)}^{\mathrm{2}}+{\left(\frac{\partial {f}_{i}}{\partial {B}_{\mathrm{1}}}{W}_{{B}_{\mathrm{1}}}\right)}^{\mathrm{2}}+{\left(\frac{\partial {f}_{i}}{\partial {B}_{\mathrm{2}}}{W}_{{B}_{\mathrm{2}}}\right)}^{\mathrm{2}}+\\ {\left(\frac{\partial {f}_{i}}{\partial {B}_{\mathrm{3}}}{W}_{{B}_{\mathrm{3}}}\right)}^{\mathrm{2}}+{\left(\frac{\partial {f}_{i}}{\partial B}{W}_{B}\right)}^{\mathrm{2}},\end{array}}\end{array}$

where fi stands for the contribution of a specific runoff component, and W is the uncertainty of the variable specified by the subscript. For the uncertainty of tracer signatures (WAi and WBi), we multiply the SD values of the measured tracer signatures with t values from the Student's t value table at the confidence level of 95 %. The degree of freedom for the Student's t distribution is estimated as the number of water samples for each water source minus one. Analytical measurement errors are not considered in this approach, which are, however, minor compared to the uncertainty generated from tracer variations (Penna et al., 2017; Pu et al., 2017). The lsqnonneg function in MATLAB is used to solve Eqs. (1)–(4), which solves the equations in a least squares sense, given the constraint that the solution vector f has nonnegative elements. The EMMA_4 uses the equations similar to Eqs. (1)–(5). The values of δ18O and δ2H are typically correlated for each water source. However, the coefficients representing the correlation between δ18O and δ2H (typically calculated as the deuterium excess values) vary among the water sources in a glacierized catchment, thus providing a basis for the EMMA_4 to quantify four runoff components. When quantifying four runoff components using three tracers, four conservative equations of water volume, EC, δ18O, and δ2H are used (similar to Eq. 1). The contributions of runoff components (f), and the partial derivatives used to calculate the uncertainty, are solved from the four conservative equations using MATLAB. However, the solutions are too lengthy to show in the text.

## 3.2 Bayesian mixing approach

The Bayesian approaches applied for each season are summarized in Table 2. Similar to the EMMA, we apply a three-component Bayesian approach to all seasons and, additionally, a four-component Bayesian approach in the glacier melt season. The three-component Bayesian approach has two types. The Bayesian_3_OHcor approach considers the correlation between δ18O and δ2H, whereas the Bayesian_3_OHind approach assumes independence. The four-component Bayesian approach also has two types, namely the Bayesian_4_OHcor, which considers the correlation, and Bayesian_4_OHind, which assumes independence between δ18O and δ2H. A Kolmogorov–Smirnov test has been carried out for both isotope and EC tracers of all water sources before the application of Bayesian approaches. The tracer data of runoff components (i.e., rainfall, snowmelt, groundwater, and glacier melt) pass the normal distribution test at significance levels of p values > 0.3, apart from the EC data of glacier melt. The small glacier melt sample size for the EC measurement probably provides insufficient data for the distribution test. The tracer data of stream water also fail to pass the normal distributions test, which is partly caused by the extreme isotope and EC values (see Figs. S1a and S1b). Thus, the prior assumptions for the Bayesian approaches are listed (similar to Cable et al. 2011) as follows. In approaches considering the correlation between δ18O and δ2H, the prior distributions of δ18O and δ2H of runoff components are assumed to be bivariate normal distributions with means and precision matrix as μ18O, μ2H, and Ω, respectively (Eq. 6a). The precision matrix (Ω, i.e., the inverse of the covariance matrix) for the two isotopes is assumed to be Wishart prior (Eq. 6b). When assuming independence between δ18O and δ2H, the prior distributions of δ18O (δ2H) from runoff components are assumed to be normal distributions with means and variance of μ18O and λ18O (μ2H and λ2H; Eqs. 6c and d). The mean values of the isotopes of runoff components (i.e., μ18O and μ2H) are further estimated by independent normal priors (Eq. 7; Cable et al. 2011), which is assumed to consider the spatial variability of μ18O and μ2H.

$\begin{array}{ll}& \left\{\begin{array}{lr}\left[\begin{array}{l}\mathit{\delta }{}^{\mathrm{18}}\mathrm{O}\\ \mathit{\delta }{}^{\mathrm{2}}\mathrm{H}\end{array}\right]\sim \text{Multi}\mathit{_}\text{normal}\left(\left[\begin{array}{l}\mathit{\mu }{}^{\mathrm{18}}\mathrm{O}\\ \mathit{\mu }{}^{\mathrm{2}}\mathrm{H}\end{array}\right],\mathbf{\Omega }\phantom{\rule{0.33em}{0ex}}\right)& \text{(6a)}\\ \mathbf{\Omega }\sim \text{Wishart}\phantom{\rule{0.25em}{0ex}}\left(\mathrm{2},\mathbf{V}\right)& \text{(6b)}\\ \mathit{\delta }{}^{\mathrm{18}}\mathrm{O}\sim \text{Normal}\phantom{\rule{0.25em}{0ex}}\left(\mathit{\mu }{}^{\mathrm{18}}\mathrm{O},\mathit{\lambda }{}^{\mathrm{18}}\mathrm{O}\right)& \text{(6c)}\\ \mathit{\delta }{}^{\mathrm{2}}\mathrm{H}\sim \text{Normal}\phantom{\rule{0.25em}{0ex}}\left(\mathit{\mu }{}^{\mathrm{2}}\mathrm{H},\mathit{\lambda }{}^{\mathrm{2}}\mathrm{H}\right)& \text{(6d)}\end{array}\right\\\ & \left\{\begin{array}{lr}\mathit{\mu }{}^{\mathrm{18}}\mathrm{O}\sim \text{Normal}\phantom{\rule{0.25em}{0ex}}\left(\mathit{\gamma }{}^{\mathrm{18}}\mathrm{O},\mathit{\sigma }{}^{\mathrm{18}}\mathrm{O}\right)& \text{(7a)}\\ \mathit{\mu }{}^{\mathrm{2}}\mathrm{H}\sim \text{Normal}\phantom{\rule{0.25em}{0ex}}\left(\mathit{\gamma }{}^{\mathrm{2}}\mathrm{H},\mathit{\sigma }{}^{\mathrm{2}}\mathrm{H}\right),& \text{(7b)}\end{array}\right\\end{array}$

where λ18O, γ18O, and σ18O (λ2H, γ2H and σ2H) are parameters used to describe the normal priors of δ18O and μ18O (δ2H and μ2H; see Table 3), which are estimated by likelihood observations. V is a 2 × 2 unit positive definite matrix, and “2” stands for the degree of freedom in the Wishart prior distribution.

The priors of EC values of runoff components are assumed to be normal distributions (Eq. 8a), with mean ε and variance τ. Similarly, the spatial variability of the mean EC values of runoff components (ε) is assumed to follow a normal distribution with mean θ and variance ω (Eq. 8b). Moreover, τ, θ and ω are parameters estimated by likelihood observations (Table 3).

$\begin{array}{ll}\text{(8a)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{EC}\sim \text{Normal}\phantom{\rule{0.25em}{0ex}}\left(\mathit{\epsilon },\mathit{\tau }\right)\text{(8b)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathit{\epsilon }\sim \text{Normal}\phantom{\rule{0.25em}{0ex}}\left(\mathit{\theta },\mathit{\omega }\right).& \left\{\begin{array}{}\\ \end{array}\right\\end{array}$

The prior distributions of stream water are calculated in two steps. First, the prior distributions of δ18O, δ2H, and EC of stream water are assumed to be the same as those of runoff components in Eqs. (6) and (8a). Second, the mean isotopes (μ18O and μ2H) and EC (ε) of stream water are constrained by a mixing model (Eq. 9a and b), which estimates the isotope and EC mean values of stream water by multiplying the contribution of each runoff component (fi) with the corresponding mean isotope and EC values of each runoff component (Eq. 9a).

$\left\{\begin{array}{ll}{\left[\begin{array}{l}\mathit{\mu }{}^{\mathrm{18}}\mathrm{O}\\ \mathit{\mu }{}^{\mathrm{2}}\mathrm{H}\\ \mathit{\epsilon }\end{array}\right]}_{\text{stream water}}& \\ ={\sum }_{i=\mathrm{1}}^{N}{f}_{i}\cdot {\left[\begin{array}{l}\mathit{\mu }{}^{\mathrm{18}}\mathrm{O}\\ \mathit{\mu }{}^{\mathrm{2}}\mathrm{H}\\ \mathit{\epsilon }\end{array}\right]}_{\text{runoff component}\phantom{\rule{0.25em}{0ex}}i}& \text{(9a)}\\ \mathbit{f}\sim \text{Dirichlet}\phantom{\rule{0.25em}{0ex}}\left(\mathbit{\alpha }\right)& \text{(9b)}\\ \mathbit{\alpha }=\mathbit{\rho }+\mathbit{\psi }& \text{(9c)}\\ \left[\mathbit{\rho },\mathbit{\psi }\right]\sim \text{Multi}\mathit{_}\text{normal}\phantom{\rule{0.25em}{0ex}}\left(\mathbit{\beta },\mathbf{\Omega }\right),& \text{(9d)}\end{array}\right\$

where N is the number of runoff components. The contribution vector (f) is represented by a Dirichlet distribution with an index vector α (Eq. 9b), in which the sum of the contributions of all runoff components (fi) equals one. The index vector α is estimated by two variable vectors, namely ρ and ψ (Eq. 9c), considering the temporal and spatial variability in the CRCs (Cable et al. 2011). ρ and ψ are assumed to be bivariate normal distribution with means and precision matrix β and Ω (Eq. 9d). Moreover, β is a parameter vector estimated by likelihood observations (Table 3).

The initial value ranges of parameters that need to be estimated in Eqs. (6)–(9) are summarized in Table 3. The posteriors of parameters describing the spatial variability of tracer signatures in Eqs. (7) and (8b) are first estimated by the mean tracer signatures of runoff components measured at different spatial locations. Parameters describing the overall variability of tracer signatures in Eqs. (6) and (8a) are then constrained by likelihood observations of tracer signatures from all water samples at different times and locations. The posterior distributions of CRCs (f) are estimated by Eq. (9), based on the posterior tracer signatures of runoff components and measured tracer signatures from stream water samples. The posteriors of parameters and contributions are estimated by the R software package RStan. We ran four parallel Markov chain Monte Carlo (MCMC) chains with 2000 iterations for each chain. The first 1000 iterations were discarded for warm-up, generating a total of 4 × 1000 samples for the calculation of the posterior distributions. Uncertainties are presented as the 5–95 percentile ranges from the iterative runs. The parameter values are assumed to follow uniform prior distributions within their value ranges to initialize the MCMC procedure.

Note that the four-component approaches (EMMA_4, Bayesian_4_OHcor, and Bayesian_4_OHind) are developed in our study to investigate the following two questions. (1) Is the EMMA able to quantify four runoff components just by using δ18O, δ2H, and EC? (2) Does the correlation between δ18O and δ2H help to reduce the uncertainty of the quantification of runoff components? The correlation between δ18O and δ2H is ignored in Bayesian_4_OHind. Thus, we used independent prior distributions for δ18O and δ2H of each water source. In Bayesian_4_OHcor, the posterior parameters describing the correlation between δ18O and δ2H vary among the water sources, thus providing a basis for the quantification of four runoff components using four mixing equations of tracer signatures, which is similar to Eq. (9).

## 3.3 Effects of the uncertainty of the meltwater sampling

Due to limited accessibility, meltwater samples are typically difficult to collect in high-elevation glacierized areas. Often only a few water samples are available to represent the tracer signatures of meltwater generated from the entire glacierized area. Hence, the representativeness of collected meltwater samples implies an additional uncertainty source in the hydrograph separation.

We thus define three virtual sampling scenarios to evaluate the effects of meltwater sampling on the EMMA and Bayesian mixing approaches. Scenario I is used to evaluate the effects of the sample size of meltwater in which four groups of meltwater sample are tested. The four sample groups have the same mean value and SD of δ18O or EC but different sample sizes. Mean and SD values of δ18O or EC are calculated for all available meltwater used samples in each group. Scenario II is used to evaluate the effects of the sampled mean value of δ18O (or EC) of meltwater. The four sample groups have the same sample size and SD but different mean values of δ18O (or EC). Scenario III is used to investigate the effects of SD values of sampled δ18O (or EC). The four sample groups have the same sample size and mean tracer signature but different SD values. Because meltwater is particularly difficult to collect, and is the dominant runoff component in the glacier melt season, we investigated the effects of the meltwater sampling uncertainty on the mixing approaches in this season. For the water samples of other runoff components and stream water, we used all the available measurements in the glacier melt season for the three virtual scenarios and kept the same sample characteristics. We only investigated the effects of sampling uncertainty in the glacier melt season because of the following reasons. (1) Runoff in the glacier melt season contributes the largest part to the annual runoff in our study basin. Accurate quantification of each runoff component in this season is extremely important for understanding the dynamics of water availability in the study area. (2) There are more meltwater samples available in this season (15 snowmelt samples and 23 glacier melt samples) than in the snowmelt season (only 15 snowmelt samples; Table 1), thus providing a good observation database for the investigation.

Table 1Tracer signatures measured from water samples in three seasons.

CV stands for coefficient of variation.

Table 2Mixing approaches applied for hydrograph separation in different seasons.

Table 3Parameters used for prior distributions in the Bayesian approaches.

## 3.4 Effects of water isotope fractionation on hydrograph separation

The water sources for runoff, such as rainfall and meltwater, are subject to evaporation before reaching the basin outlet – especially in summer. However, the isotopic composition of stream water was measured at the basin outlet, and the contributions of runoff components are also quantified for the total runoff at the basin outlet. After the long routing path from the sampled sites to the basin outlet, the isotopic compositions of rainfall and meltwater mixing at the basin outlet could be different from those measured at the sampled sites, which is caused by the evaporation fractionation effect. To consider the changes in the isotope signatures of water sources caused by the fractionation effect during the mixing process, we set up two modified Bayesian approaches, i.e., Bayesian_3_OHcor_Frac and Bayesian_4_OHcor_Frac (Table 2). The fractionation effect on the estimated CRCs is quantified by comparing two Bayesian scenarios. In the first scenario (using Bayesian_3_OHcor and Bayesian_4_OHcor), the isotopic compositions of water sources at the basin outlet are assumed to be the same as those measured from the sample sites even though the water sources have suffered evaporation before reaching the basin outlet (using Eqs. 6–9). In the second scenario (using Bayesian_3_OHcor_Frac and Bayesian_4_OHcor_Frac), the evaporation fractionation effect on the isotopic compositions of water sources is considered, and the mixing of water tracers for stream water is represented by Eq. (10). We modify the mean values in Eq. (9a) using the fractionation factors of ξ18O and ξ2H. The priors for ξ18O and ξ2H are assumed to be bivariate normal distributions in Eq. (11).

$\begin{array}{ll}& {\left[\begin{array}{l}\mathit{\mu }{}^{\mathrm{18}}\mathrm{O}\\ \mathit{\mu }{}^{\mathrm{2}}\mathrm{H}\end{array}\right]}_{\text{stream water}}\\ \text{(10)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}=\sum _{i=\mathrm{1}}^{N}{f}_{i}\cdot {\left[\begin{array}{l}\mathit{\mu }{}^{\mathrm{18}}\mathrm{O}+\mathit{\xi }{}^{\mathrm{18}}\mathrm{O}\\ \mathit{\mu }{}^{\mathrm{2}}\mathrm{H}+\mathit{\xi }{}^{\mathrm{2}}\mathrm{H}\end{array}\right]}_{\text{runoff component}\phantom{\rule{0.25em}{0ex}}i}\text{(11)}& & \left[\begin{array}{l}\mathit{\xi }{}^{\mathrm{18}}\mathrm{O}\\ \mathit{\xi }{}^{\mathrm{2}}\mathrm{H}\end{array}\right]\sim \text{Multi}\mathit{_}\text{normal}\phantom{\rule{0.25em}{0ex}}\left(\left[\begin{array}{l}\mathit{\eta }{}^{\mathrm{18}}\mathrm{O}\\ \mathit{\eta }{}^{\mathrm{2}}\mathrm{H}\end{array}\right],\mathbf{\Omega }\right),\end{array}$

where η18O and η2H are parameters describing the mean values of the changes in isotopes caused by the fractionation effect. Ω is the inverse of the covariance matrix defined in Eq. (6b). The parameters in Eqs. (6)–(11) are then reestimated, using the MCMC procedure, by the measurements of tracer signatures. In particular, parameters describing the prior distributions of isotopic compositions at the sample sites in Eqs. (6) and (7) are estimated by the likelihood observations of isotope signatures of runoff components. The fractionation factors ξ18O and ξ2H are estimated by the likelihood observations of isotope signatures of stream water.

4 Results

## 4.1 Seasonality of tracer signatures

Tracer measurements from all the water samples are summarized in Table 1 and Fig. 2 (see also Fig. S1). The mean values in Table 1 indicate that precipitation is most depleted in heavy water isotopes (18O and 2H) among the water sources in the cold season. In the melt seasons, snow and glacier meltwater show the most depleted heavy isotopes. The EC values are highest in groundwater in all seasons and are followed by stream water and precipitation. Among the water sources, snowmelt and glacier melt tend to have the lowest EC values. Figure 2 shows that the slope of the local meteoric water line (LMWL) is lower than that of the global meteoric water line (GMWL). The δ18O of precipitation and snowmelt range from 22.82 ‰ to 1.51 ‰ and from 17.31 ‰ to 6.95 ‰, respectively. The isotopic composition of glacier meltwater is more depleted than that of groundwater and stream water. Stream water shows a similar isotopic composition to groundwater. Three samples from the stream water are far below the LMWL, which is likely caused by the evaporation effect.

Figure 2Isotope signatures of water samples from the three seasons in the Ala-Archa basin.

CV values in Table 1 and boxplots in Fig. 3a–f show that the δ18O and δ2H of precipitation generally show the largest variability in all seasons, followed by the isotopes of snowmelt. Groundwater and stream water show the smallest CV values for δ18O in all three seasons. The stream water presents the lowest CV value for EC in all seasons, followed by the groundwater. The snowmelt EC shows high CV values in the snowmelt and glacier melt seasons, which may be attributed to variable dust conditions at the sampling locations (from the downstream gauge station to the upper glacier accumulation zone). The highest CV value of EC for glacier melt indicates large variability in the glacier melt samples (also see Fig. 3g–i). This is because the glacier melt water samples were collected from a rather clean location (EC value is only 1.5 µs cm−1) and also a relatively dusty location (EC value is 33.4 µs cm−1).

Figure 3(a–i) Boxplots of tracer signatures in three seasons. (j–l) δ18O–EC mixing space of the various water sources in the three seasons. The solid lines indicate the ranges of tracer signatures measured from water samples.

For each water source, except groundwater, the tracer signatures show a significant seasonality (Table 1 and Fig. 3). In particular, the δ18O and δ2H of precipitation are most depleted in the cold season and reach the highest values in the glacier melt season, which is partly caused by the seasonality of the temperature. Stream water shows higher values of δ18O and EC in the cold season when groundwater dominates the streamflow and has lower values in the melt seasons when meltwater has a dominant contribution. Snowmelt has a lower EC value in the glacier melt season than in the cold and snowmelt seasons. In the cold and snowmelt seasons, some snowmelt samples also have EC values as low as those in the glacier melt season. The snow samples in the glacier melt season were only collected from the accumulation zone of the glacier, thus resulting in a small variability in the EC values. The snowpack in the accumulation zone is accumulated by fresh snow in the snowy period (summer-type accumulation glacier).This leads to low EC values in the snowmelt samples. The tracer signature of groundwater is relatively stable across the seasons.

Figure 3j–l show the δ18O–EC mixing space of runoff components in the three seasons. The ranges of solid lines indicate the minimum and maximum tracer values of individual water samples. In the cold season, the δ18O and EC values of stream water are very close to those of groundwater (Fig. 3j), whereas the snowmelt and precipitation tracer signatures show much more of a difference. These results indicate the dominance of groundwater on streamflow during the cold season. In the snowmelt and glacier melt seasons (Fig. 3k–l), the stream water samples are clearly located within the triangle formed by the samples of runoff components. The tracer signatures of glacier meltwater and snowmelt water are similar. The precipitation samples are further away from the stream water samples when compared to the meltwater and groundwater samples. The stream water samples are located nearly in the middle between the meltwater and groundwater samples. This indicates that the contribution of rainfall to total runoff is the smallest, and the contributions of meltwater and groundwater are similar in the melt seasons.

## 4.2 Contributions of runoff components estimated by the mixing approaches

Table 4 and Fig. 4 compare the CRCs estimated by the mixing approaches. In the cold season (Fig. 4a), the EMMA_3 estimated the mean contributions of groundwater and snowmelt as 83 % and 17 %, respectively. The mean contribution of rainfall is zero. The mean contributions of groundwater, snowmelt, and rainfall were estimated as 86 % (87 %), 13 % (12 %), and 1 % (1 %) by the Bayesian_3_OHind (Bayesian_3_OHcor) approach. As shown in Fig. 3j, the tracer signature of stream water in this season is close to that of groundwater, while obviously different from that of rainfall. Meanwhile, the stream water samples are outside of the triangle formed by the runoff components, leading to the zero contribution of the rainfall estimated by the EMMA_3.

Table 4Contributions of runoff components (CRCs) estimated by the different mixing approaches (percentage, %). The ranges (%) show the difference between the 95 % and 5 %  percentiles. SD values refer to the standard deviations.

Figure 4Contributions of runoff components (CRCs) to total runoff estimated by different mixing approaches in three seasons. The Bayesian_3_OHind and Bayesian_3_OHcor were applied in the cold and melt seasons (a–c), and the Bayesian_4_OHind and Bayesian_4_OHcor were applied in the glacier melt season (d). The horizontal lines in the boxes refer to the median contributions, and the whiskers refer to the 95 % and 5 % percentiles.

In the snowmelt season (Fig. 4b and Table 4), the EMMA_3 estimated the mean contributions of groundwater, rainfall, and snowmelt as 44 %, 36 %, and 20 %, respectively. The Bayesian_3_OHind estimated similar mean CRCs to EMMA_3, whereas the Bayesian_3_OHcor delivered a lower contribution of snowmelt (32 %). When treating the glacier melt and snowmelt as one end member (i.e., meltwater) in the glacier melt season (Fig. 4c), the EMMA_3 estimated the mean contributions of groundwater, meltwater, and rainfall as 45 %, 46 %, and 9 %, respectively. The Bayesian_3_OHind and Bayesian_3_OHcor estimated a lower contribution of groundwater (43–44 %) and a higher contribution of rainfall (11 %) compared to EMMA_3. The ranges and SD values of CRCs in Table 4 indicate the uncertainty of the estimates associated with the corresponding mixing approaches, showing that the EMMA_3, followed by Bayesian_3_OHind, produced the highest uncertainty of CRCs in all the three seasons. The Bayesian_3_OHcor reduced the uncertainty slightly when compared to Bayesian_3_OHind, as it benefits from the consideration of the correlation between δ18O and δ2H.

When treating glacier melt and snowmelt as two separate end members in the glacier melt seasons (Fig. 4d), the EMMA_4 failed to separate the hydrograph in the glacier melt season given the large uncertainty range in the contributions of snowmelt and rainfall (0–100 %). The tracer signatures of snow and glacier meltwater are rather close to each other, which violates the second assumption of the EMMA (see Sect. 3.1). In contrast, the Bayesian_4_OHcor and Bayesian_4_OHind estimated the shares of glacier melt and snowmelt as 25–24 % and 21–25 %, respectively. Considering the significant snow cover area in September in the study basin (He et al. 2018, 2019), the contribution of snowmelt in the glacier melt season should be higher than zero. Again, the Bayesian_4_OHcor produced smaller uncertainty ranges and SD values for the contributions of groundwater and meltwater compared to Bayesian_4_OHind and EMMA_4 (Table 4).

The posterior distributions of tracer signatures estimated by the Bayesian_4_OHcor in the glacier melt season are compared with the measured histograms of tracer signatures in Fig. 5. The Bayesian_4_OHcor generally produced similar distributions of water isotopes to the measured distributions in terms of the similar mean values. The estimated posterior SD values of the water isotopes are smaller than SD values of the measurements. This can be explained by the incorporation of prior distributions by the Bayesian_4_OHcor, which reduces the variability of water isotopes. The posterior SD values for EC of water sources are also smaller than the measured SD values. However, the posterior distributions of EC show some deviations from the distributions of measured EC (Fig. 5k–o), which is partly due to the very small sample sizes (see Table 1). The comparison between the posterior distributions of tracer signatures estimated by the Bayesian_3_OHcor and the measured distributions in the other seasons generally shows a similar behavior (not shown for brevity).

Figure 5Posterior distributions of tracer signatures estimated by the Bayesian_4_OHcor in the glacier melt season. Measurement refers to the distributions of tracer signatures from the water samples. (a–e) Distributions of δ18O; (f–j) distributions of δ2H; and (k–o) distributions of EC.

The Bayesian_4_OHind estimated similar posterior distributions of tracer signatures to the Bayesian_4_OHcor (except the glacier melt isotopes; Fig. 6), with similar mean tracer values. It is noted that the Bayesian_4_OHcor estimated smaller SD values for most water sources compared to the Bayesian_4_OHind (e.g., Fig. 6f, g, i, and j). Benefiting from the prior information and the consideration of the correlation between δ18O and δ2H, the Bayesian_4_OHcor tended to produce the smallest variability in the posterior tracer signatures among all the mixing approaches (Figs. 5 and 6), thus resulting in the smallest uncertainty for CRCs (Fig. 4d). Figure 7 compares the correlation between δ18O and δ2H of the measured tracers and the posterior estimates by Bayesian approaches. The Bayesian_4_OHcor reproduced the correlation between δ18O and δ2H well in comparison to the measured data, whereas the Bayesian_4_OHind failed to capture the correlation.

Figure 6Comparison of the posterior distributions of tracer signatures estimated by the Bayesian approaches with Bayesian_4_OHcor and without Bayesian_4_OHind considering the correlation between δ18O and δ2H in the glacier melt season.

Figure 7Correlation between posterior δ18O and δ2H estimated by the Bayesian_4_OHcor and the Bayesian_4_OHind approaches in the glacier melt season.

## 4.3 Uncertainty of hydrograph separation caused by sampling uncertainty of meltwater

Figure 8 shows the sensitivity of the Bayesian_3_OHcor and EMMA_3 approaches to the sampled δ18O of meltwater in the glacier melt season. The mean CRCs quantified by the two mixing approaches shows minor sensitivity to the sample size (Scenario I). However, the uncertainty ranges of contributions tend to decrease with increasing sample sizes, especially for EMMA_3. When assuming only two meltwater samples, the EMMA_3 resulted in very large uncertainty ranges (0–100 %; Fig. 8d), due to the very wide confidence interval for the SD at a sample size of two. The mean contributions of groundwater and meltwater estimated by the two mixing approaches decrease with the increasing mean δ18O of the adopted meltwater sample (Scenario II), while the estimated contribution of rainfall increases with the increasing mean δ18O (Fig. 8k). Variations in the mean CRCs quantified by EMMA_3 are larger than those estimated by the Bayesian_3_OHcor. Using EMMA_3, both the mean contributions of groundwater and meltwater declined by 9 % with the assumed increase of the mean δ18O (Fig. 8e and h), and the contribution of rainfall increased by 17 %. Using Bayesian_3_OHcor, the reduction of the contributions of groundwater and snowmelt are 4 % and 7 %, respectively, and the increase of the contribution of rainfall is only 11 % (Fig. 8k). In Scenario III, the uncertainty ranges of CRCs (especially for rainfall; Fig. 8l) increase with the increasing SD of the sampled δ18O. Again, the increases in the uncertainty ranges estimated by EMMA_3 tend to be larger than those estimated by the Bayesian_3_OHcor. The sensitivity of the mixing approaches to the sampled EC values of the meltwater are similar to the sensitivity of the sampled δ18O (not shown).

Figure 8Sensitivity of the CRCs estimates to the sample size (Scenario I), the mean (Scenario II), and standard deviation (Scenario III) of δ18O of meltwater samples in the glacier melt season. Red boxes show the contributions estimated by the Bayesian_3_OHcor, and the blue boxes refer to the contributions estimated by the EMMA_3.

## 4.4 Effect of isotope fractionation on the hydrograph separation

The changes of δ18O caused by the fractionation effect (referring to ξ18O in Eq. 10) during the mixing process are estimated in Fig. 9a–c. The fractionation has the smallest effect on the δ18O of groundwater, while having the largest effect on the δ18O of rainfall. On average, the δ18O of rainfall increased by around 2.8 ‰ through fractionation in all the three seasons. The CRCs estimated by the Bayesian_3_OHcor_Frac and Bayesian_4_OHcor_Frac are compared with those estimated by the Bayesian_3_OHcor and Bayesian_4_OHcor in Fig. 9d–f, respectively. The mean contribution of groundwater estimated by the Bayesian_3_OHcor_Frac in the cold season is around 9 % lower than that estimated by the Bayesian_3_OHcor (Fig. 9d), while the mean contributions of snowmelt and rainfall are 3 % and 5 % higher, respectively. The reduction of the groundwater contribution should be attributed to the increased contributions of snowmelt and rainfall caused by the fractionation effect. In the snowmelt season, the mean contributions of groundwater and rainfall are 1 % and 7 % lower, respectively (Fig. 9e), while the mean contribution of snowmelt estimated by the Bayesian_3_OHcor_Frac is 8 % higher. In the glacier melt season, the mean contributions of groundwater and meltwater estimated by the Bayesian_4_OHcor_Frac are higher than those estimated by the Bayesian_4_OHcor (Fig. 9f) and are compensated by the 6 % lower contribution of rainfall.

Figure 9Effects of isotope fractionation on the estimates of CRCs in the Bayesian approach for the three seasons. (a–c) Estimated changes in δ18O of runoff components caused by the fractionation effect; (d–e) CRCs estimated by the Bayesian_3_OHcor and the Bayesian_3_OHcor_Frac; and (f) CRCs estimated by the Bayesian_4_OHcor and the Bayesian_4_OHcor_Frac.

The fractionation effect also produced visible changes in the posterior distributions of δ18O and δ2H of runoff components (Fig. 10 shows the example in the glacier melt season). The mean isotopic compositions of runoff components are increased by the fractionation effect. The SD values of the posterior isotopes estimated by the Bayesian_4_OHcor_Frac tend to be higher than those estimated by the Bayesian_4_OHcor due to the increased parameter space in the prior assumptions (Eq. 11), thus leading to the larger uncertainty ranges in the contributions of glacier melt and snowmelt (Fig. 9f). As expected, the estimates of posterior distributions of isotopic compositions of stream water are less sensitive to the fractionation effect of runoff components (Fig. 10e and j). The fractionation also has minor effects on the estimates of posterior distributions of EC values (Fig. 10k–o).

Figure 10Effects of isotope fractionation on the estimated posterior distributions of tracer signatures of water sources in the glacier melt season.

5 Discussion

## 5.1 Uncertainty of the contributions of runoff components

The EMMA estimated similar CRCs but with a larger uncertainty than the Bayesian approaches. The reasons for this are twofold. First, the EMMA estimated the uncertainty ranges of CRCs using the standard deviations (SD) of the measured tracer signatures. SD values are likely overestimated in this study due to the small sample sizes (i.e., low number of water samples) and thus represent the variability of the tracer signatures of the corresponding water sources across the basin insufficiently. Due to the limited accessibility of the sample sites caused by snow cover, the water samples of meltwater and groundwater are often collected sporadically. The small sample size and strong variability in sampled tracer signatures likely led to a large SD value in the measurement. Second, the EMMA assumes that the uncertainty associated with each water source is independent of the uncertainty of other water sources (Eq. 5), which increases the uncertainty ranges for CRCs.

In contrast, by updating the prior probability distributions, the Bayesian approaches estimated a smaller variability of tracer signatures in the posterior distributions when compared to the measured tracer signatures. The posterior distributions were sampled continuously from the assumed initial value ranges by the MCMC runs, thus reducing the sharp changes and yielding lower variability for the tracer signatures. Moreover, the uncertainty ranges for CRCs were quantified using Eqs. (6)–(10) instead of calculating independently as in the EMMA. Additionally, the assumed prior distributions of tracer signatures and the CRCs take the correlation between the tracer signatures and the dependence between the runoff components into account, thus resulting in smaller uncertainty ranges (Soulsby et al., 2003). For example, the Bayesian approaches that considered the correlation between δ18O and δ2H generally estimated smaller uncertainty ranges for CRCs compared to those that did not consider this correlation.

The Gaussian error propagation technique is only capable of considering the uncertainty of CRCs resulting from the variation in the tracer signatures (Uhlenbrook and Hoeg, 2003). The uncertainty of CRCs that originated from the sampling uncertainty of meltwater was then investigated in separate virtual sampling experiments. The EMMA produces large uncertainty ranges and SD values for CRCs in the glacier melt season when the meltwater sample size is rather small. The mean CRC quantified by the EMMA relies more heavily on the mean tracer values of the sampled meltwater, since the mean tracer values have been used directly in Eqs. (1)–(4), compared to the mean CRC estimated by the Bayesian approach.

The EMMA assumes that the tracer signature of each runoff component is constant during the mixing process; thus, it is unable to estimate the uncertainty of CRCs caused by the isotope fractionation effect. The virtual fractionation experiments, using the modified Bayesian approaches, show that the isotope fractionation could increase the contribution of snowmelt by 8 % and reduce the contribution of rainfall by 7 % in the snowmelt season. We assume that the mean CRCs estimated by the Bayesian approaches that consider the isotope fractionation are more plausible – despite the larger uncertainty ranges. Along the flow path from the source areas to the river channel, the isotopic compositions of meltwater and rainfall are likely increased by the evaporation fractionation effect – especially in the warm seasons. The increased isotopic compositions of meltwater and rainfall during the routing process need to be considered in the mixing approaches for hydrograph separation.

In general, the uncertainty of CRCs is visibly caused by the spatiotemporal variability in the tracer signatures, the water sampling uncertainty, and the isotope fractionation during the mixing process. The uncertainty caused by the water sampling of meltwater tends to be smaller than the uncertainty caused by the variations of the tracer signatures in both the EMMA and Bayesian mixing approaches. This is consistent with the findings that the SD values of the tracer measurements of water samples are the main uncertainty sources for the quantification of CRCs (Schmieder et al., 2016, 2018). The Bayesian approach tends to be superior for narrowing the variability of posterior tracer signatures benefitting from the prior assumptions and the consideration of the dependence between tracer signatures and runoff components when compared to EMMA.

## 5.2 Limitations

The representativeness of the water samples is one of the limitations of this study. The groundwater was only sampled from a single spring located at an elevation of 2400 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, which is rather close to the average altitude of the entire river network in the study basin (2530 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). We thus assume that the measured isotopic composition of the spring water represents the mean isotopic composition of the groundwater feeding the river in the basin (see also He et al., 2019). Collecting samples from a few springs to represent the groundwater end member has been proposed before (e.g., Ohlanders et al., 2013 and Mark and McKenzie, 2007), as the accessibility and availability of more potential springs are hampered. Again, for the snow and glacier meltwater samples, we assume that meltwater occurring at similar elevations has similar tracer signatures (He et al., 2019). The sampled elevation ranges from 1580 to 4050 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ and matches the elevation range where meltwater mainly occurs in the basin (from 1580 to 3950 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). Considering that the isotopic compositions of meltwater are particularly dependent on the elevation, the sampled meltwater could represent meltwater that originated from the primary melting locations in the entire basin. The sampled sites thus bear the potential to provide tracer signatures of the major meltwater generated in the basin.

We split the entire sampling period (from 2012 to 2017) into three seasons, i.e., cold season, snowmelt season, and glacier melt season, due to the low availability of the water samples in each year. By concentrating water samples in the three seasons, we increased the sample sizes of each runoff component for each season, thus increasing the ability of water samples to represent the spatiotemporal variability of seasonal tracer signatures. We used all available groundwater and snowmelt samples from the three seasons for hydrograph separation in the cold season due to the rather low number of samples collected in the cold season. This likely leads to overestimated contributions from groundwater and snowmelt in the cold season. However, the overestimation of the groundwater contribution is probably small because the tracer signatures of groundwater generally show small seasonal variability. The estimated contributions of snowmelt in the cold season are a bit higher than the contribution modeled by He et al (2018) during winter months of December, January, and February; these are still reasonable when considering that the cold season includes October and November when the snow is more prone to melting.

The assumptions of the mixing approaches lead to another limitation of this study. The EMMA assumes the tracer signatures of water sources are constant during the mixing process, which is a common assumption for the practical application of EMMA. It thus fails to consider the uncertainty originating from the changes of tracer signatures. In the Bayesian approach, we assumed normal prior distributions for the tracer signatures of water sources and Dirichlet prior distribution for the CRCs based on the literature (Cable et al., 2011). To refine the description of the temporal and spatial variability of the CRCs in the Dirichlet distribution, more hydrological data relating to the runoff processes in the basin are required. We acknowledge that the estimated CRCs could be strongly affected by the assumptions of prior distributions. However, testing the effects of the prior assumptions goes beyond the scope of this study. We assume that collecting more water samples from various locations and at different times for each water source could improve the estimation of tracer signature distributions.

6 Conclusions

This study compared the Bayesian end-member mixing approach with a traditional end-member mixing approach (EMMA) for hydrograph separation in a glacierized basin. The contributions of runoff components (CRCs) to the total runoff were estimated for three seasons, i.e., cold season, snowmelt, and glacier melt seasons. The mean CRCs estimated by the two mixing approaches are similar in all the three seasons. The uncertainties of these contributions, caused by the variability of tracer signatures, water sampling uncertainty, and isotope fractionation, were evaluated as follows:

1. The Bayesian approach generally estimates smaller uncertainty ranges of CRCs in comparison to the EMMA. Benefiting from the prior assumptions of tracer signatures and CRCs, and from the incorporation of the correlation between tracer signatures in the prior distributions, the Bayesian approach reduced the uncertainty. The Bayesian approach jointly quantified the uncertainty ranges of CRCs. In contrast, the EMMA estimated the uncertainty of the contribution of each runoff component independently, thus leading to higher uncertainty ranges.

2. The estimates of CRCs in EMMA tend to be more sensitive to the sampling uncertainty of meltwater when compared to those in the Bayesian approach. For small sample sizes (e.g., two), EMMA estimated very large uncertainty ranges. The mean of the CRCs quantified by EMMA is also more sensitive to the mean value of the tracer signature of the sampled meltwater than those values estimated by the Bayesian approach.

3. Ignoring the isotope fractionation during the mixing process likely overestimates the contribution of rainfall and underestimates the contribution of meltwater in the melt seasons. The EMMA currently used is unable to quantify the uncertainty of CRCs caused by the isotope fractionation during the mixing process due to the underlying assumptions.

Code availability
Code availability.

RStan code for the Bayesian end-member mixing approach is available at https://doi.org/10.5281/zenodo.3897266 (He, 2020).

Supplement
Supplement.

Author contributions
Author contributions.

ZH, KUS, and SV conceptualized this research. ZH, KUS, SMW, OK, and AG collected the data. ZH, KUS, and SV developed the methodology that was used. The original draft was compiled by ZH, SV, and DD. All authors contributed to writing the review and the editing of the paper.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Financial support
Financial support.

This paper has been funded by the German Federal Ministry for Science and Education (GlaSCA-V; grant no. 88 501) and the Volkswagen Foundation (GlaSCA; grant nos. 01DK15002A and B).

The article-processing charges for this open-access publication were provided by a Helmholtz Association German Research Centre.

Review statement
Review statement.

This paper was edited by Markus Weiler and reviewed by three anonymous referees.

References

Aizen, V. B., Aizen, E. M., and Melack, J. M.: Precipitation, melt and runoff in the northern Tien Shan, J. Hydrol., 186, 229–251, https://doi.org/10.1016/S0022-1694(96)03022-3, 1996.

Aizen, V., Aizen, E., Glazirin, G., and Loaiciga, H. A.: Simulation of daily runoff in Central Asian alpine watersheds, J. Hydrol., 238, 15–34, https://doi.org/10.1016/S0022-1694(00)00319-X, 2000.

Aizen, V. B., Kuzmichenok, V. A., Surazakov, A. B., and Aizen, E. M.: Glacier changes in the Tien Shan as determined from topographic and remotely sensed data, Global Planet. Change, 56, 328–340, https://doi.org/10.1016/j.gloplacha.2006.07.016, 2007.

Barnett, T. P., Adam, J. C., and Lettenmaier, D. P.: Potential impacts of a warming climate on water availability in snow-dominated regions, Nature, 438, 303–309, https://doi.org/10.1038/nature04141, 2005.

Beria, H., Larsen, J. R., Michelon, A., Ceperley, N. C., and Schaefli, B.: HydroMix v1.0: a new Bayesian mixing framework for attributing uncertain hydrological sources, Geosci. Model Dev., 13, 2433–2450, https://doi.org/10.5194/gmd-13-2433-2020, 2020.

Brown, L. E., Hannah, D. M., Milner, A. M., Soulsby, C., Hodson, A. J., and Brewer, M. J.: Water source dynamics in a glacierized alpine river basin (Taillon-Gabietous, French Pyrenees), Water Resour. Res., 42, W08404, https://doi.org/10.1029/2005WR004268, 2006.

Cable, J., Ogle, K., and Williams, D.: Contribution of glacier meltwater to streamflow in the Wind River Range, Wyoming, inferred via a Bayesian mixing model applied to isotopic measurements, Hydrol. Process., 25, 2228–2236, https://doi.org/10.1002/hyp.7982, 2011.

Chiogna, G., Santoni, E., Camin, F., Tonon, A., Majone, B., Trenti, A., and Bellin, A.: Stable isotope characterization of the Vermigliana catchment, J. Hydrol., 509, 295–305, https://doi.org/10.1016/j.jhydrol.2013.11.052, 2014.

Coplen, T. B. and Wassenaar, L. I.: LIMS for Lasers 2015 for achieving longterm accuracy and precision of δ2H, δ17O, and δ18O of waters using laser absorption spectrometry, Rapid Commun. Mass Sp., 29, 2122–2130, https://doi-org.cyber.usask.ca/10.1002/rcm.7372, 2015.

Dahlke, H. E., Lyon, S. W., Jansson, P., Karlin, T., and Rosqvist, G.: Isotopic investigation of runoff generation in a glacierized catchment in northern Sweden, Hydrol. Process., 28, 1383–1398, https://doi.org/10.1002/hyp.9668, 2014.

Engel, M., Penna, D., Bertoldi, G., Dell'Agnese, A., Soulsby, C., and Comiti, F.: Identifying run-off contributions during melt-induced run-off events in a glacierized alpine catchment, Hydrol. Process., 30, 343–364, https://doi.org/10.1002/hyp.10577, 2016.

Genereux, D.: Quantifying uncertainty in tracer-based hydrograph separations, Water Resour. Res., 34, 915–919, https://doi.org/10.1029/98WR00010, 1998.

Gröning, M., Lutz, H. O., Roller-Lutz, Z., Kralik, M., Gourcy, L., and Pöltenstein, L.: A simple rain collector preventing water reevaporation dedicated for δ18O and δ2H analysis of cumulative precipitation samples, J. Hydrol., 448–449, 195–200, https://doi.org/10.1016/j.jhydrol.2012.04.041, 2012.

He, Z. H., Parajka, J., Tian, F. Q., and Blöschl, G.: Estimating degree-day factors from MODIS for snowmelt runoff modeling, Hydrol. Earth Syst. Sci., 18, 4773-4789, https://doi.org/10.5194/hess-18-4773-2014, 2014.

He, Z. H., Tian, F. Q., Gupta, H. V., Hu, H. C., and Hu, H. P.: Diagnostic calibration of a hydrological model in a mountain area by hydrograph partitioning, Hydrol. Earth Syst. Sci., 19, 1807–1826, https://doi.org/10.5194/hess-19-1807-2015, 2015.

He, Z., Vorogushyn, S., Unger-Shayesteh, K., Gafurov, A., Kalashnikova, O., Omorova, E., and Merz, B.: The Value of Hydrograph Partitioning Curves for Calibrating Hydrological Models in Glacierized Basins, Water Resour. Res., 54, 2336–2361, https://doi.org/10.1002/2017WR021966, 2018.

He, Z., Unger-Shayesteh, K., Vorogushyn, S., Weise, S. M., Kalashnikova, O., Gafurov, A., Duethmann, D., Barandun, M., and Merz, B.: Constraining hydrological model parameters using water isotopic compositions in a glacierized basin, Central Asia, J. Hydrol., 571, 332–348, https://doi.org/10.1016/j.jhydrol.2019.01.048, 2019.

He, Z. H.: Bayesian end-member mixing approach, Zenodo, https://doi.org/10.5281/zenodo.3897266, 2020.

Joerin, C., Beven, K. J., Iorgulescu, I., and Musy, A.: Uncertainty in hydrograph separations based on geochemical mixing models, J. Hydrol., 255, 90–106, https://doi.org/10.1016/S0022-1694(01)00509-1, 2002.

Klaus, J. and McDonnell, J. J.: Hydrograph separation using stable isotopes: Review and evaluation, J. Hydrol., 505, 47–64, https://doi.org/10.1016/j.jhydrol.2013.09.006, 2013.

Kong, Y. L. and Pang, Z. H.: Evaluating the sensitivity of glacier rivers to climate change based on hydrograph separation of discharge, J. Hydrol., 434, 121–129, https://doi.org/10.1016/j.jhydrol.2012.02.029, 2012.

La Frenierre, J. and Mark, B. G.: A review of methods for estimating the contribution of glacial meltwater to total watershed discharge, Prog. Phys. Geog., 38, 173–200, https://doi.org/10.1177/0309133313516161, 2014.

Li, Z. X., Feng, Q., Liu, W., Wang, T. T., Cheng, A. F., Gao, Y., Guo, X. Y., Pan, Y. H., Li, J. G., Guo, R., and Jia, B.: Study on the contribution of cryosphere to runoff in the cold alpine basin: A case study of Hulugou River Basin in the Qilian Mountains, Global Planet. Change, 122, 345–361, https://doi.org/10.1016/j.gloplacha.2014.10.001, 2014.

Mark, B. G. and Mckenzie, J. M.: Tracing increasing tropical Andean glacier melt with stable isotopes in water, Environ. Sci. Technol., 41, 6955–6960, https://doi.org/10.1021/es071099d, 2007.

Maurya, A. S., Shah, M., Deshpande, R. D., Bhardwaj, R. M., Prasad, A., and Gupta, S. K.: Hydrograph separation and precipitation source identification using stable water isotopes and conductivity: River Ganga at Himalayan foothills, Hydrol. Process., 25, 1521–1530, https://doi.org/10.1002/hyp.7912, 2011.

Moore, J. W. and Semmens, B. X.: Incorporating uncertainty and prior information into stable isotope mixing models, Ecol. Lett., 11, 470–480, https://doi.org/10.1111/j.1461-0248.2008.01163.x, 2008.

Ohlanders, N., Rodriguez, M., and McPhee, J.: Stable water isotope variation in a Central Andean watershed dominated by glacier and snowmelt, Hydrol. Earth Syst. Sci., 17, 1035–1050, https://doi.org/10.5194/hess-17-1035-2013, 2013.

Parnell, A. C., Inger, R., Bearhop, S., and Jackson, A. L.: Source Partitioning Using Stable Isotopes: Coping with Too Much Variation, PLoS ONE, 5, e9672, https://doi.org/10.1371/journal.pone.0009672, 2010.

Penna, D., Engel, M., Mao, L., Dell'Agnese, A., Bertoldi, G., and Comiti, F.: Tracer-based analysis of spatial and temporal variations of water sources in a glacierized catchment, Hydrol. Earth Syst. Sci., 18, 5271–5288, https://doi.org/10.5194/hess-18-5271-2014, 2014.

Penna, D., van Meerveld, H. J., Zuecco, G., Fontana, G. D., and Borga, M.: Hydrological response of an Alpine catchment to rainfall and snowmelt events, J. Hydrol., 537, 382–397, https://doi.org/10.1016/j.jhydrol.2016.03.040, 2016.

Penna, D., Engel, M., Bertoldi, G., and Comiti, F.: Towards a tracer-based conceptualization of meltwater dynamics and streamflow response in a glacierized catchment, Hydrol. Earth Syst. Sci., 21, 23–41, https://doi.org/10.5194/hess-21-23-2017, 2017.

Pohl, E., Gloaguen, R., Andermann, C., and Knoche, M.: Glacier melt buffers river runoff in the Pamir Mountains, Water Resour. Res., 53, 2467–2489, https://doi.org/10.1002/2016WR019431, 2017.

Pu, T., He, Y. Q., Zhu, G. F., Zhang, N. N., Du, J. K., and Wang, C. F.: Characteristics of water stable isotopes and hydrograph separation in Baishui catchment during the wet season in Mt. Yulong region, south western China, Hydrol. Process., 27, 3641–3648, https://doi.org/10.1002/hyp.9479, 2013.

Pu, T., Qin, D. H., Kang, S. C., Niu, H. W., He, Y. Q., and Wang, S. J.: Water isotopes and hydrograph separation in different glacial catchments in the southeast margin of the Tibetan Plateau, Hydrol. Process., 31, 3810–3826, https://doi.org/10.1002/hyp.11293, 2017.

Rahman, K., Besacier-Monbertrand, A. L., Castella, E., Lods-Crozet, B., Ilg, C., and Beguin, O.: Quantification of the daily dynamics of streamflow components in a small alpine watershed in Switzerland using end member mixing analysis, Environ. Earth Sci., 74, 4927–4937, https://doi.org/10.1007/s12665-015-4505-5, 2015.

Rodriguez, M., Ohlanders, N., Pellicciotti, F., Williams, M. W., and McPhee, J.: Estimating runoff from a glacierized catchment using natural tracers in the semiarid Andes cordillera, Hydrol. Process., 30, 3609–3626, https://doi.org/10.1002/hyp.10973, 2016.

Schmieder, J., Hanzer, F., Marke, T., Garvelmann, J., Warscher, M., Kunstmann, H., and Strasser, U.: The importance of snowmelt spatiotemporal variability for isotope-based hydrograph separation in a high-elevation catchment, Hydrol. Earth Syst. Sci., 20, 5015–5033, https://doi.org/10.5194/hess-20-5015-2016, 2016.

Schmieder, J., Garvelmann, J., Marke, T., and Strasser, U.: Spatio-temporal tracer variability in the glacier melt end-member – How does it affect hydrograph separation results, Hydrol. Process., 32, 1828–1843, https://doi.org/10.1002/hyp.11628, 2018.

Soulsby, C., Petry, J., Brewer, M. J., Dunn, S. M., Ott, B., and Malcolm, I. A.: Identifying and assessing uncertainty in hydrological pathways: a novel approach to end member mixing in a Scottish agricultural catchment, J. Hydrol., 274, 109–128, https://doi.org/10.1016/S0022-1694(02)00398-0, 2003.

Sun, C. J., Chen, Y. N., Li, W. H., Li, X. G., and Yang, Y. H.: Isotopic time series partitioning of streamflow components under regional climate change in the Urumqi River, northwest China, Hydrolog. Sci. J., 61, 1443–1459, https://doi.org/10.1080/02626667.2015.1031757, 2016a.

Sun, C. J., Yang, J., Chen, Y. N., Li, X. G., Yang, Y. H., and Zhang, Y. Q.: Comparative study of streamflow components in two inland rivers in the Tianshan Mountains, Northwest China, Environ. Earth Sci., 75, 727, https://doi.org/10.1007/s12665-016-5314-1, 2016b.

Uhlenbrook, S., and Hoeg, S.: Quantifying uncertainties in tracer-based hydrograph separations: a case study for two-,three- and five-component hydrograph separations in a mountainous catchment, Hydrol. Process., 17, 431–453, https://doi.org/10.1002/hyp.1134, 2003.

Viviroli, D., Durr, H. H., Messerli, B., Meybeck, M., and Weingartner, R.: Mountains of the world, water towers for humanity: Typology, mapping, and global significance, Water Resour. Res., 43, W07447, https://doi.org/10.1029/2006wr005653, 2007.

Ward, E. J., Semmens, B. X., and Schindler, D. E.: Including Source Uncertainty and Prior Information in the Analysis of Stable Isotope Mixing Models, Environ. Sci. Technol., 44, 4645–4650, https://doi.org/10.1021/es100053v, 2010.