Technical note: Uncertainty in multi-source partitioning using large tracer data sets

The availability of large tracer data sets opened up the opportunity to investigate multiple source contributions to a mixture. However, the source contributions may be uncertain and, apart from Bayesian approaches, to date there are only solid methods to estimate such uncertainties for two and three sources. We introduce an alternative uncertainty estimation method for four sources based on multiple tracers as input data. Taylor series approximation is used to solve the set of linear mass balance equations. We illustrate the method to compute individual uncertainties in the calculation of source contributions to a mixture, with an example from hydrology, using a 14-tracer set from water sources and streamflow from a tropical, high-elevation catchment. Moreover, this method has the potential to be generalized to any number of tracers across a range of disciplines.

1 Introduction

Tracer applications have dramatically increased over recent years across a wide range of disciplines (West et al., 2010). Applications in hydrology (Hooper, 2003; James and Roulet, 2006; Kirchner and Neal, 2013), ecology (Phillips and Gregg, 2003; Semmens et al., 2009b), anthropology (Ehleringer et al., 2008), conservation biology (Bicknell et al., 2014), nutrition (Magaña-Gallegos et al., 2018), environmental and ecosystem science (Bartov et al., 2013; Granek et al., 2009), and erosion and sediment transportation (Davies et al., 2018) have been the most prominent. Such widespread use of tracers was mainly facilitated by the availability of analytical techniques that provide highly sensitive, rapid multi-element analysis at a lower cost (Falkner et al., 1995). For example, the use of inductively coupled plasma mass spectrometry (ICP-MS) as one of the most relevant analytical techniques for elemental analysis (Helaluddin et al., 2016) led to the availability and use of large tracers sets (elements) in hydrological studies (Barthold et al., 2017; Belli et al., 2017; Correa et al., 2017; Kirchner and Neal, 2013; Mimba et al., 2017). Trace elements together with water stable isotopes (cavity ringdown laser absorption spectroscopy paved the way: Berman et al., 2009; Lis et al., 2008) as well as physical–chemical water parameters (e.g. electrical conductivity and pH) are now often used to improve understanding of hydro-geochemical cycles, flow pathways and runoff generation in hydrology. Furthermore, mixing models based on tracer mass balance equations are widely applied to identify the dominant sources of a mixture and their contribution dynamics.

In hydrological mixing models, the composition of the stream is assumed to be an integrated mixture of signatures of different sources (Christophersen et al., 1990). The proportional contributions of n+1 sources to the stream can be uniquely determined using n different tracers (Christophersen and Hooper, 1992). Bayesian methods have been developed to identify multiple (>3) sources and compute their contributions to a mixture in a two-dimensional mixing space (Parnell et al., 2010; Stock et al., 2018). In this case, a unique solution is not feasible and a higher uncertainty is attributed to the model (Phillips and Gregg, 2001, 2003). On the other hand, end-member mixing analysis (EMMA) (Hooper, 2003) was developed to use multiple tracers as input and, therefore, allows for a multi-dimensional space that potentially increases the number of identifiable sources (Barthold et al., 2011; Inamdar et al., 2013; Liu et al., 2004). In addition, the use of multiple tracers can avoid bias and subjectivity in the input information. Therefore, EMMA provides a robust and complete conceptualization of catchment functioning and source interactions during runoff generation (Iwasaki et al., 2015). However, despite its benefits, the EMMA approach lacks a formal methodology to assess the uncertainty of multiple end-members (Delsman et al., 2013) and their individual uncertainties in the calculation of source contributions to a stream.

To our knowledge, the uncertainty estimation of source contributions to streams is based on Gaussian error propagation (Genereux, 1998) and was so far only calculated using one or two tracers simultaneously (MixSIAR: Parnell et al., 2010; Phillips and Gregg, 2001; Semmens et al., 2009a). Alternatively, when the number of sources is higher, the uncertainty is usually based on the sum of analytical errors, elevation effects and the spatial variability of end-member concentrations (Uhlenbrook and Hoeg, 2003). Hence, we propose an alternative methodology based on the first-order Taylor series approximation to estimate the uncertainty of individual end-members or sources (e.g. precipitation, soil water, groundwater) to a mixture (e.g. streamflow).

We illustrate this application using a multi-tracer data set from Correa et al. (2019b), in a three-dimensional space defined by a principal component analysis (PCA). In Correa et al. (2019b), the authors computed the uncertainties but without disclosing any details in the calculation and methodology used. The main objective of this technical note is, therefore, to explicitly describe the mathematical development that allows the calculation of partial derivatives, degrees of freedom and confidence interval limits for each source fraction contribution and, moreover, to provide the code and several examples for their calculation and reproducibility.

2 Uncertainty estimation method development

In this section, the uncertainty estimation method presented in Phillips and Gregg (2001) is expanded for four source contributions to the mixture.

Let 𝒞 represent the following set of sources: A, B, C and D, and mixture M, $\mathcal{C}=\left\{A,B,C,D,M\right\}$. In the following equations, x∈𝒞, $y\in \left\{\stackrel{\mathrm{‾}}{\mathit{\delta }},\stackrel{\mathrm{‾}}{\mathit{\lambda }},\stackrel{\mathrm{‾}}{\mathit{\varphi }}\right\}$ and $z\in \left\{A,M,C\right\}$. x, y and z are variables that belong to the following sets: x to the set of A, B, C, D and mixture M; y to the set of medians of every projected source and mixture in each principal component $\stackrel{\mathrm{‾}}{\mathit{\delta }}\stackrel{\mathrm{‾}}{\mathit{\lambda }}\stackrel{\mathrm{‾}}{\mathit{\varphi }}$ respectively of the sub index used; and z to the set of A, M and C. Furthermore, $f_A,f_B,f_C$ and f_D represent the contribution fraction of sources A, B, C and D respectively to the mixture M.

The data required for this analysis are the median and standard deviations (σ) of each of the sources (A, B, C and D) and the mixture M, projected and expressed in the coordinates of the three-dimensional PCA space. As well as this, the sample size (n) of each source is required. Details of the application are presented in Sect. 3.2.

If the system is composed of Eq. (1)

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{ccccccccc}{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_A{f}_A& +& {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_B{f}_B& +& {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_C{f}_C& +& {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_D{f}_D& =& {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_M\\ {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_A{f}_A& +& {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_B{f}_B& +& {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_C{f}_C& +& {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_D{f}_D& =& {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_M\\ {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_A{f}_A& +& {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_B{f}_B& +& {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_C{f}_C& +& {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_D{f}_D& =& {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_M\\ f_A& +& f_B& +& f_C& +& f_D& =& \mathrm{1}\end{array}\right.\end{array}$$

and has solution¹ for f_Af_Bf_Cf_D>0, the contribution fractions take the following form:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& f_A={\displaystyle \frac{\left({\mathrm{\Phi }}_M-{\mathrm{\Delta }}_M\right)\left({\mathrm{\Lambda }}_C-{\mathrm{\Delta }}_C\right)-\left({\mathrm{\Lambda }}_M-{\mathrm{\Delta }}_M\right)\left({\mathrm{\Phi }}_C-{\mathrm{\Delta }}_C\right)}{\left({\mathrm{\Phi }}_A-{\mathrm{\Delta }}_A\right)\left({\mathrm{\Lambda }}_C-{\mathrm{\Delta }}_C\right)-\left({\mathrm{\Lambda }}_A-{\mathrm{\Delta }}_A\right)\left({\mathrm{\Phi }}_C-{\mathrm{\Delta }}_C\right)}}\\ & ={\displaystyle \frac{\mathrm{Num}}{\mathrm{Den}}},\\ & f_C={\displaystyle \frac{\left({\mathrm{\Delta }}_M-{\mathrm{\Lambda }}_M\right)-\left({\mathrm{\Delta }}_A-{\mathrm{\Lambda }}_A\right)f_A}{\left({\mathrm{\Delta }}_C-{\mathrm{\Lambda }}_C\right)}},\\ & f_B={\mathrm{\Delta }}_M-\left({\mathrm{\Delta }}_C{f}_C+{\mathrm{\Delta }}_A{f}_A\right),\\ & f_D=\mathrm{1}-\left(f_C+f_B+f_A\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\mathrm{\Delta }}_x={\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x-{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_D}{{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_B-{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_D}},{\mathrm{\Lambda }}_x={\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x-{\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_D}{{\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_B-{\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_D}},\\ & {\mathrm{\Phi }}_x={\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x-{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_D}{{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_B-{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_D}}.\end{array}\end{array}$$

The partial derivatives of Eq. (2) are given by

$$\begin{array}{rl}& {\displaystyle \frac{\partial f_A}{\partial y_x}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{Den}}^{\mathrm{2}}}}\left[\right[\left({\mathrm{\Lambda }}_C-{\mathrm{\Delta }}_C\right)\left({\displaystyle \frac{\partial {\mathrm{\Phi }}_M}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_M}{\partial y_x}}\right)\\ & +\left({\mathrm{\Phi }}_M-{\mathrm{\Delta }}_M\right)\left({\displaystyle \frac{\partial {\mathrm{\Lambda }}_C}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_C}{\partial y_x}}\right)\\ & -\left({\mathrm{\Phi }}_C-{\mathrm{\Delta }}_C\right)\left({\displaystyle \frac{\partial {\mathrm{\Lambda }}_M}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_M}{\partial y_x}}\right)\\ & -\left({\mathrm{\Lambda }}_M-{\mathrm{\Delta }}_M\right)\left({\displaystyle \frac{\partial {\mathrm{\Phi }}_C}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_C}{\partial y_x}}\right)]\mathrm{Den}\\ & -[\left({\mathrm{\Lambda }}_C-{\mathrm{\Delta }}_C\right)\left({\displaystyle \frac{\partial {\mathrm{\Phi }}_A}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_A}{\partial y_x}}\right)\\ & +\left({\mathrm{\Phi }}_A-{\mathrm{\Delta }}_A\right)\left({\displaystyle \frac{\partial {\mathrm{\Lambda }}_C}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_C}{\partial y_x}}\right)\\ & -\left({\mathrm{\Phi }}_C-{\mathrm{\Delta }}_C\right)\left({\displaystyle \frac{\partial {\mathrm{\Lambda }}_A}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_A}{\partial y_x}}\right)\\ & -\left({\mathrm{\Lambda }}_A-{\mathrm{\Delta }}_A\right)\left({\displaystyle \frac{\partial {\mathrm{\Phi }}_C}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_C}{\partial y_x}}\right)\left]\mathrm{Num}\right],\end{array}$$

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\displaystyle \frac{\partial f_C}{\partial y_x}}={\displaystyle \frac{\mathrm{1}}{{\left({\mathrm{\Delta }}_C-{\mathrm{\Lambda }}_C\right)}^{\mathrm{2}}}}\left[\right[\left({\displaystyle \frac{\partial {\mathrm{\Delta }}_M}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Lambda }}_M}{\partial y_x}}\right)\\ & -\left({\displaystyle \frac{\partial {\mathrm{\Delta }}_A}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Lambda }}_A}{\partial y_x}}\right)f_A-\left({\mathrm{\Delta }}_A-{\mathrm{\Lambda }}_A\right){\displaystyle \frac{\partial f_A}{\partial y_x}}]\left({\mathrm{\Delta }}_C-{\mathrm{\Lambda }}_C\right)\\ & -\left({\displaystyle \frac{\partial {\mathrm{\Delta }}_C}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Lambda }}_C}{\partial y_x}}\right)\left[\left({\mathrm{\Delta }}_M-{\mathrm{\Lambda }}_M\right)-\left({\mathrm{\Delta }}_A-{\mathrm{\Lambda }}_A\right)f_A\right]],\\ & {\displaystyle \frac{\partial f_B}{\partial y_x}}={\displaystyle \frac{\partial {\mathrm{\Delta }}_M}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_C}{\partial y_x}}{f}_C-{\mathrm{\Delta }}_C{\displaystyle \frac{\partial f_C}{\partial y_x}}-{\displaystyle \frac{\partial {\mathrm{\Delta }}_A}{\partial y_x}}{f}_A\\ & -{\mathrm{\Delta }}_A{\displaystyle \frac{\partial f_A}{\partial y_x}},\\ & {\displaystyle \frac{\partial f_D}{\partial y_x}}=-{\displaystyle \frac{\partial f_C}{\partial y_x}}-{\displaystyle \frac{\partial f_B}{\partial y_x}}-{\displaystyle \frac{\partial f_A}{\partial y_x}}.\end{array}\end{array}$$

It is trivial that

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\displaystyle \frac{\partial {\mathrm{\Delta }}_z}{\partial w_x}}=\mathrm{0},w\in \left\{\stackrel{\mathrm{‾}}{\mathit{\lambda }},\stackrel{\mathrm{‾}}{\mathit{\varphi }}\right\};{\displaystyle \frac{\partial {\mathrm{\Lambda }}_z}{\partial w_x}}=\mathrm{0},w\in \left\{\stackrel{\mathrm{‾}}{\mathit{\delta }},\stackrel{\mathrm{‾}}{\mathit{\varphi }}\right\};\\ & {\displaystyle \frac{\partial {\mathrm{\Phi }}_z}{\partial w_x}}=\mathrm{0},w\in \left\{\stackrel{\mathrm{‾}}{\mathit{\delta }},\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right\},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial {\mathrm{\Delta }}_z}{\partial {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}}={\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_B-{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_D\right)}^{-\mathrm{1}}\left\{\begin{array}{ll}\mathrm{1}& z\in \left\{A,C,M\right\}\mathrm{and}x=z\\ -{\mathrm{\Delta }}_z& z\ne B\mathrm{and}x=B\\ {\mathrm{\Delta }}_z-\mathrm{1}& z\ne D\mathrm{and}x=D\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right.,\text{(7)}& {\displaystyle \frac{\partial {\mathrm{\Lambda }}_z}{\partial {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}}={\left({\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_B-{\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_D\right)}^{-\mathrm{1}}\left\{\begin{array}{ll}\mathrm{1}& z\in \left\{A,C,M\right\}\mathrm{and}x=z\\ -{\mathrm{\Lambda }}_z& z\ne B\mathrm{and}x=B\\ {\mathrm{\Lambda }}_z-\mathrm{1}& z\ne D\mathrm{and}x=D\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right.\end{array}$$

and

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial {\mathrm{\Phi }}_z}{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}}={\left({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_B-{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_D\right)}^{-\mathrm{1}}\left\{\begin{array}{ll}\mathrm{1}& z\in \left\{A,C,M\right\}\mathrm{and}x=z\\ -{\mathrm{\Phi }}_z& z\ne B\mathrm{and}x=B\\ {\mathrm{\Phi }}_z-\mathrm{1}& z\ne D\mathrm{and}x=D\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right..\end{array}$$

For example, for f_A we have

$$\begin{array}{rl}& {\displaystyle \frac{\partial f_A}{\partial {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{Den}}^{\mathrm{2}}}}\left[\right[{\displaystyle \frac{\partial {\mathrm{\Delta }}_M}{\partial {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}}\left({\mathrm{\Phi }}_C-{\mathrm{\Lambda }}_C\right)\\ & -{\displaystyle \frac{\partial {\mathrm{\Delta }}_C}{\partial {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}}\left({\mathrm{\Phi }}_M-{\mathrm{\Lambda }}_M\right)]\mathrm{Den}-[{\displaystyle \frac{\partial {\mathrm{\Delta }}_A}{\partial {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}}\left({\mathrm{\Phi }}_C-{\mathrm{\Lambda }}_C\right)\\ & -{\displaystyle \frac{\partial {\mathrm{\Delta }}_C}{\partial {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}}\left({\mathrm{\Phi }}_A-{\mathrm{\Lambda }}_A\right)\left]\mathrm{Num}\right],\\ & {\displaystyle \frac{\partial f_A}{\partial {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{Den}}^{\mathrm{2}}}}\left[\right[{\displaystyle \frac{\partial {\mathrm{\Lambda }}_C}{\partial {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}}\left({\mathrm{\Phi }}_M-{\mathrm{\Delta }}_M\right)\\ & -{\displaystyle \frac{\partial {\mathrm{\Lambda }}_M}{\partial {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}}\left({\mathrm{\Phi }}_C-{\mathrm{\Delta }}_C\right)]\mathrm{Den}-[{\displaystyle \frac{\partial {\mathrm{\Lambda }}_C}{\partial {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}}\left({\mathrm{\Phi }}_A-{\mathrm{\Delta }}_A\right)\\ & -{\displaystyle \frac{\partial {\mathrm{\Lambda }}_A}{\partial {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}}\left({\mathrm{\Phi }}_C-{\mathrm{\Delta }}_C\right)\left]\mathrm{Num}\right],\end{array}$$

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\displaystyle \frac{\partial f_A}{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{Den}}^{\mathrm{2}}}}\left[\right[{\displaystyle \frac{\partial {\mathrm{\Phi }}_M}{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}}\left({\mathrm{\Lambda }}_C-{\mathrm{\Delta }}_C\right)\\ & -{\displaystyle \frac{\partial {\mathrm{\Phi }}_C}{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}}\left({\mathrm{\Lambda }}_M-{\mathrm{\Delta }}_M\right)]\mathrm{Den}-[{\displaystyle \frac{\partial {\mathrm{\Phi }}_A}{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}}\left({\mathrm{\Lambda }}_C-{\mathrm{\Delta }}_C\right)\\ & -{\displaystyle \frac{\partial {\mathrm{\Phi }}_C}{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}}\left({\mathrm{\Lambda }}_A-{\mathrm{\Delta }}_A\right)\left]\mathrm{Num}\right].\end{array}\end{array}$$

Using Eq. (9), the first-order Taylor series approximation (Taylor, 1982) for the variance (σ²) of f_A evaluated at the mean can be calculated by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathit{\sigma }}_{f_A}^{\mathrm{2}}=\sum _x{\left({\displaystyle \frac{\partial f_A}{\partial {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}^{\mathrm{2}}+\sum _x{\left({\displaystyle \frac{\partial f_A}{\partial {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}^{\mathrm{2}}\\ & +\sum _x{\left({\displaystyle \frac{\partial f_A}{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}^{\mathrm{2}}=\sum _y\sum _x{\left({\displaystyle \frac{\partial f_A}{\partial y_x}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{y_x}^{\mathrm{2}}.\end{array}\end{array}$$

To calculate γ_A (the Satterthwaite, 1946, approximation for the degrees of freedom), we define $f_{A{y}_x}=c_A{\left(\frac{\partial f_A}{\partial y_x}\right)}^{\mathrm{2}}$, where c_A is a scale constant that relates $f_{A{y}_x}$ with the respective derivative. It means that f_A with respect to y_x can be a scalar multiple of the derivative value.

In this case, we get

$$\begin{array}{}\text{(11)}& {\mathit{\gamma }}_A={\displaystyle \frac{{\left(\sum _y\sum _x{f}_{A{y}_x}{\mathit{\sigma }}_{y_x}^{\mathrm{2}}\right)}^{\mathrm{2}}}{\sum _y\sum _x\frac{{\left(f_{A{y}_x}{\mathit{\sigma }}_{y_x}^{\mathrm{2}}\right)}^{\mathrm{2}}}{n_{y_x}-\mathrm{1}}}}.\end{array}$$

Note that whatever the value of c_A is, Eq. (11) leads to

$${\mathit{\gamma }}_A={\displaystyle \frac{{\left(\sum _y\sum _x{\left(\frac{\partial f_A}{\partial y_x}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{y_x}^{\mathrm{2}}\right)}^{\mathrm{2}}}{\sum _y\sum _x\frac{{\left({\left(\frac{\partial f_A}{\partial y_x}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{y_x}^{\mathrm{2}}\right)}^{\mathrm{2}}}{n_{y_x}-\mathrm{1}}}},$$

and if we set $f_{A{y}_x}^{\ast }={\left(\frac{\partial f_A}{\partial y_x}\right)}^{\mathrm{2}}$ then the numerator of the last equation can be replaced by ${\left({\mathit{\sigma }}_{f_A}^{\mathrm{2}}\right)}^{\mathrm{2}}$. In other words, we can use Eq. (10) and the derivatives Eq. (9) to estimate the value of γ_A, resulting in $f_{A{y}_x}=c_A{f}_{A{y}_x}^{\ast }$. Of course, it is required that c_A is constant with respect to y_x. Then,

$$\begin{array}{}\text{(12)}& {\mathit{\gamma }}_A={\displaystyle \frac{{\left({\mathit{\sigma }}_{f_A}^{\mathrm{2}}\right)}^{\mathrm{2}}}{\sum _y\sum _x\frac{{\left({\left(\frac{\partial f_A}{\partial y_x}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{y_x}^{\mathrm{2}}\right)}^{\mathrm{2}}}{n_{y_x}-\mathrm{1}}}}.\end{array}$$

Let $w\in \mathcal{C}\setminus \left\{A\right\}$. The first-order Taylor series approximation for the variance of f_w can be calculated by (as above Eq. 10)

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& {\mathit{\sigma }}_{f_w}^{\mathrm{2}}=\sum _x{\left({\displaystyle \frac{\partial f_w}{\partial {\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_x}^{\mathrm{2}}+\sum _x{\left({\displaystyle \frac{\partial f_w}{\partial {\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_x}^{\mathrm{2}}\\ & +\sum _x{\left({\displaystyle \frac{\partial f_w}{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_x}^{\mathrm{2}}=\sum _y\sum _x{\left({\displaystyle \frac{\partial f_w}{\partial y_x}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{y_x}^{\mathrm{2}}.\end{array}\end{array}$$

If we construct γ_w as γ_A, we get

$${\mathit{\gamma }}_w={\displaystyle \frac{{\left(\sum _y\sum _x{f}_{w{y}_x}^{\ast }{\mathit{\sigma }}_{y_x}^{\mathrm{2}}\right)}^{\mathrm{2}}}{\sum _y\sum _x\frac{{\left(f_{w{y}_x}^{\ast }{\mathit{\sigma }}_{y_x}^{\mathrm{2}}\right)}^{\mathrm{2}}}{n_{y_x}-\mathrm{1}}}},$$

where $f_{w_{y_x}}=c_w{f}_{w_{y_x}}^{\ast }$ and $f_{w_{y_x}}^{\ast }={\left(\frac{\partial f_w}{\partial y_x}\right)}^{\mathrm{2}}$ with c_w constant with respect to y_x, then we finally get

$$\begin{array}{}\text{(14)}& {\mathit{\gamma }}_w={\displaystyle \frac{{\left({\mathit{\sigma }}_{f_w}^{\mathrm{2}}\right)}^{\mathrm{2}}}{\sum _y\sum _x\frac{{\left({\left(\frac{\partial f_w}{\partial y_x}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{y_x}^{\mathrm{2}}\right)}^{\mathrm{2}}}{n_{y_x}-\mathrm{1}}}}.\end{array}$$

The upper and lower confidence interval limits for each end-member fraction can be calculated using partial derivatives and the 95 % confidence intervals constructed as follows:

$$\begin{array}{}\text{(15)}& f_w\pm t_{\mathrm{0.05}{\mathit{\gamma }}_w}{\mathit{\sigma }}_{fw},\end{array}$$

where t_0.05,γ is the Student's t for α=0.05 (two-tailed) (Walpole et al., 2017) and γ degrees of freedom related with σ_fw.

Figure 1 Three-dimensional mixing space generated using stream data, where the median of end-members are projected. U1 represents 59.6 % of the variance, U2 19.7 % and U3 7.4 % (From PCA); RF: rainfall; AN: Andosols; HS: Histosols; SW: spring water; M: the median of stream data (mixture).

Table 1 Median and standard deviation (SD) of end-members and stream projected in three-dimensional space for the study period 2013–2014.

^∗ Coordinates of end-members and stream (mixture) medians in three-dimensional
space (U1, U2 and U3). n represents the sample size.

Figure 2 Boxplots of end-members projected in the three-dimensional mixing space for the study period 2013–2014. The y axis represents the coordinates of the mixing space and the x axis the principal components U1, U2 and U3 (the central bar in the box represents the median; notches represent the 95 % confidence intervals; whiskers 1.5 times the interquartile range and circles represent outliers). SW: spring water; HS: Histosol; AN: Andosol; RF: rainfall.

Table 2 Median and standard deviation (SD) of end-members and stream projected in three-dimensional space considering 50 % of the data sets.

The example (1) considers the initial 50 % and (2) the remaining 50 % of the sample sets. ^∗ Coordinates of end-members and stream (mixture) medians in three-dimensional space (U1, U2 and U3). n represents the sample size.

3 Application

3.1 Study site and data

This methodology was tested using data from a high-elevation (3500–3900 m a.s.l.) tropical catchment (7.53 km²) located in southern Ecuador (3^∘4^′38^′′ S, 79^∘15^′30^′′ W). The mean annual precipitation for this study site is 1300 mm (Padrón et al., 2015), the mean annual discharge is 860 mm yr⁻¹. The catchment is of a volcanic origin and dominated by volcanic Histosol (24 %) and Andosol (72 %) soils (Quichimbo et al., 2012), both with a high percentage of organic matter content (450 and 310 g kg, respectively) (Quichimbo et al., 2012) and large water-holding capacities (Buytaert et al., 2006). Histosols are primarily located at the valleys and covered by cushion plants, while Andosol soils are predominant on the hillslopes under a cover of tussock grass. Nearly saturated conditions of the riparian zone are observed year-round, and a spring is located in the north-western part of the catchment. Streamwater samples from five nested streams were collected weekly from March 2013 to April 2014 (n=257) and at a higher frequency during experimental campaigns. Additionally, bi-weekly water samples from four potential end-members, rainfall (RF), soil water from Andosols (AN) and Histosols (HS), and spring water (SW) (n∼30, for each end-member) were collected. A multi-tracer (14 tracers) data set of conservative tracers was obtained from each water sample (Na, Mg, Al, Si, K, Ca, Rb, Sr, Ba, Ce, V, Y, Nd) at the Institute for Landscape Ecology and Resource Management of the Justus Liebig University using an ICP-MS (Agilent 7500ce, Agilent Technologies) and the electrical conductivity (EC) was measured in situ. More detailed information on the study site and data set can be found in Correa et al. (2017, 2019b).

3.2 Uncertainty estimation of water source contributions

Using the classic EMMA approach (Christophersen and Hooper, 1992), data from end-members SW, HS, AN, RF and stream M were projected into a three-dimensional space (Correa et al., 2019b) and presented in Fig. 1. The resulting median and standard deviation of end-members and stream coordinates are shown in Table 1. Furthermore, Fig. 2 shows the distribution of projected samples from individual end-members in the PCA coordinates.

Table 3 Median and standard deviation (SD) of end-members and stream projected in three-dimensional space including artificial outliers.

The example (3) considers outliers included at the positive extreme of the data set of each source and (4) outliers included at the negative extreme. ^∗ Coordinates of end-members and stream (mixture) medians in three-dimensional space (U1, U2 and U3). n represents the sample size.

Table 4 Median and enlarged standard deviation (SD) of end-members and stream projected in three-dimensional space.

The example (5) considers 3 times the standard deviation of the original data set and (6) 5 times the standard deviation of the original data set. ^∗ Coordinates of end-members and stream (mixture) medians in three-dimensional space (U1, U2 and U3). n represents the sample size.

The uncertainty range of each of the four end-member contributions to the stream was determined using the above developed Eq. (15) based on the first-order Taylor series approximation (Eq. 14) (MATLAB code in Correa et al., 2019a). The f_w gives the proportion of w in M and ${\mathit{\sigma }}_{f_w}^{\mathrm{2}}$ gives the variances of w. The upper uncertainty limit was computed as $f_w+t_{\mathrm{0.05},\mathit{\gamma }}{\mathit{\sigma }}_{fw}$ and the lower limit as $f_w-t_{\mathrm{0.05},\mathit{\gamma }}{\mathit{\sigma }}_{fw}$. This procedure was applied to all end-members. The resulting uncertainty estimates for each source end-member are shown in Table 5.

Note that the set of sources A, B, C and D used for the development of the equations are represented here by SW, HS, AN and RF in this specific order. U1, U2 and U3 represent the principal components PC1, PC2 and PC3, respectively.

Table 5 Uncertainty of individual end-member contributions to the stream and Satterthwaite (1946) approximation for the degrees of freedom calculated for the study period 2013–2014.

3.3 Sensitivity of water source uncertainty to input data

From the above-mentioned data set, we have generated six examples to assess the sensitivity of the uncertainty calculation to the source sample size, the artificial inclusion of outliers (upper and lower extremes) and the increased standard deviations of the source data sets.

Table 6 Uncertainty of individual end-member contributions to the stream and Satterthwaite (1946) approximation for the degrees of freedom computed considering 50 % of the data sets.

The example (1) was computed considering the initial 50 % and (2) the remaining 50 % of the sample sets.

Table 7 Uncertainty of individual end-member contributions to the stream and Satterthwaite (1946) approximation for the degrees of freedom computed after including artificial outliers.

The example (3) was computed after including outliers at the positive extreme of the data set and (4) including outliers at the negative extreme.

Table 8 Uncertainty of individual end-member contributions to the stream and Satterthwaite (1946) approximation for the degrees of freedom computed with enlarged standard deviations.

The example (5) was computed considering 3 times the standard deviation of the original data set and (6) 5 times the standard deviation of the original data set.

The first example considers 50 % of the samples (collected in the first half of the monitoring period) from each source. The median, standard deviation and sample size are input data (Table 2) to calculate the uncertainty bands (Table 6).

The second example considers the remaining 50 % of the samples and was similarly executed (Table 2).

In the third example, outliers were artificially included at the upper positive end of data sets for each source at each coordinate, respectively. The outliers consisted of twice the maximum positive value of the observed data (Table 3).

Using the same criteria, the negative extremes were included in the fourth example (Table 3).

Sources affected by dispersed data clouds were taken into account by an increase in the standard deviation. We considered two cases, the first, in example five, increasing 3 times the value of the standard deviation of the initial data set (Table 4) and finally, increasing the standard deviation 5 times for the sixth example (Table 4).

The results of this analysis are presented in Tables 6–8. In examples 1 and 2 the sample size reduction from 24 to 12 and 13 samples respectively (Table 6) had a minimal effect (less than 3 %) on the calculation of the uncertainty ranges compared to the original complete set (Table 1). The fractions of source contributions did not experience changes. The inclusion of outliers affected the values of the medians at levels of the second decimal (Table 3) concerning the median of the initial data (Table 2). However, the standard deviations increased in a range of 1.2 to 2.5 times the original value for AN and HS, and more for RF (2.5 to 10.5) and drastically more for SW (4 to 20 times wider). These variations were reflected in the widening (1 % to 12 %) of uncertainty bands for all existing cases (Table 7) in comparison with those calculated from the original data set (Table 5). Furthermore, the widening of the standard deviations to 3 and 5 times their initial values resulted in an increase in the range of uncertainty between 2 % and 22 % for the first case and between 5 % and 37 % for the second case. For the latter, the minimum limit of the uncertainty range was reached in all the reported cases. The large number of samples used in these exercises was reflected in high degrees of freedom.

4 Summary and remarks

Our methodology was developed to calculate the contribution of sources to the mixture and its associated uncertainty (based on multiple tracer sets) and was shown to be effective in real application cases. The application of the method reflected that the calculations of the uncertainty ranges of multiple source contributions to a mixture do not experience significant changes with sample size reduction or inclusion of outliers. Rather, it shows marginally different results by incorporating standard deviations from widely dispersed data.

The methodology, based on Phillips and Gregg (2001) combined with EMMA applications (Hooper, 2003) presents high potential for use as an alternative method to the simple sum of analytical errors (Uhlenbrook and Hoeg, 2003) or the Bayesian approach (Parnell et al., 2010; Stock et al., 2018). We provide a tool to close the gap in studies of mixing processes when a larger number of source contributions (>3) and related uncertainty estimates are needed for a more complete conceptualization (Iwasaki et al., 2015).

The MATLAB code provided and the illustrative examples facilitate the understanding of the methodology and promote future scientific applications. We are confident that the use of this methodology will help the scientific community that is increasingly using large tracer sets in its research to obtain results supported by a sound uncertainty analysis.