Benchmark tests for separating n time components of runoff with one stable isotope tracer

A validation of the recently introduced iterative extension of the standard two-component hydrograph separation method is presented. The data for testing this method are retrieved from a random rainfall generator and a rainfall-runoff model composed of linear reservoirs. The results show that it is possible to reconstruct the simulated event water response of a given random model input by applying the iterative separation model and using a single stable isotope tracer. The benchmark model also covers the partially delayed response of event water so that a situation can be simulated in which pre-event water is rapidly 5 mobilized. It is demonstrated how mathematical constraints, such as an ill-conditioned linear equation system, may influence the separation of the event water response. In addition, it is discussed how the volume weighted separated event water response can serve as an estimator for a time-varying backward travel time distribution.


Introduction
Over the past few decades, the separation of storm hydrographs using stable isotope tracers has become a standard method for 10 investigating runoff generation processes in catchment hydrology. It is an approach that relies on steady-state mass balance equations, whereby the early pioneering work was accomplished in the late 1960s and 1970s (Pinder andJones (1969), Dincer et al. (1970), Martinec et al. (1974), and Fritz et al. (1976), Sklash and Farvolden (1979)), over the years, the methodology has been progressively expanded and adapted to the given challenges and tasks in the field. On the catchment scale, for instance, stable isotope tracers (deuterium, oxygen-18, or tritium) have been combined with geochemical tracers, such as a sense of the reliability of the method.
3. A consideration of how the results of the iterative separation approach can be related to the description of travel times in catchment hydrology, which is expected to be one of the favorite applications, but has not been analyzed in detail yet.

Separation of n time components 70
Consider a control volume, for instance, a catchment in a river basin, with the following bulk water balance: where S(t) is the time evolution of the water storage, J(t) is the precipitation, ET (t) is the evapotranspiration, and Q t (t) is the total stream discharge. Let c be a conservative isotope tracer with the bulk mass balance: 75 where Q t , c t , and c J are the only measured physical quantities. In addition, there are the time points t 0 , t 1 , . . . , t n and time intervals [t 0 , t 1 [, [t 1 , t 2 [, . . . , [t n−1 , t n [ that describe the start and end of n rainfall-runoff events along the time axis. Furthermore, there is a semantic time measure with the intervals e (event) and p (pre-event) that can be moved across the rainfall-runoff events, whereas the interval p is the range of all intervals just before interval e. The stream discharge Q t during event e is composed of water of the current rainfall-runoff event and of the prior rainfall-runoff events, such that where c e and c p are bulk tracer concentrations in the event and pre-event components, respectively.
According to Hoeg (2019), the separation of n time components in stream discharge Q t is based on an iterative balance of catchment input and output tracer mass flows along the time axis. Hereby, the volume flows of precipitation J and evapotranspiration ET and change in water storage S are not explicitly involved. Instead, it is postulated that the pre-event water fraction 85 of each rainfall-runoff event is entirely composed of the event water and pre-event water of the last event. Therefore, for a first backward iteration, we have: Q p = Q e-1 + Q p-1 c p Q p = c e-1 Q e-1 + c p-1 Q p-1 (4) 3 https://doi.org/10.5194/hess-2021-213 Preprint. Discussion started: 10 May 2021 c Author(s) 2021. CC BY 4.0 License.
In the literature, the following criteria are mentioned in relation to the standard two component separation (Sklash and Farvolden (1979), Buttle (1994)), which also apply to the iterative separation model (5): 1. The isotopic content of the event component is significantly different from that of the pre-event component. 3. The pre-event component maintains a constant isotopic signature in space and time, or any variations can be accounted for.
The criteria mentioned above represent a subset of the criteria usually belonging to investigations that consider also the flow paths of the separated runoff components (Klaus and McDonnell (2013)).

100
I would like to add another criterion (Criterion 4) that is usually implicitly considered and simply demands that both the event water Q e and pre-event water Q p cannot be less than zero or larger than the total runoff Q t . Given the equations above, this is the case if the tracer concentration in total runoff c t is always between that of the event water c e and pre-event water c p .
In the context of separation model (5), it is required for all the backward iterations τ that c e < c t < c p ∨ c p < c t < c e (6) significant impact on the results of the hydrograph separation, as I will discuss later.
Based on the knowns Q t , c t and c e , c e-1 . . . c e-τ and c p , c p-1 . . . c p-τ , we can iteratively derive for τ ∈ N backward iterations the following solutions: In addition, from the linear equation system (5), we get the following composition of the stream discharge Q t : 2.2 Separated event water response as an estimator for the time varying backward travel time distribution 120 The iterative separation model (5) can be used to trace event water over a much longer period after the initial event that is being. For the obtained result, I use the term separated event water response to emphasize the fact that the traced response relates to the water of exactly one rainfall event. Conceptually, it can be related to the more commonly used term travel time distribution. Basically, it can be a time-varying approximation for it.
The travel time distribution is the response or breakthrough of an instantaneous, conservative tracer addition over the entire 125 catchment area. It is a probability distribution that can be derived analytically based on the physical assumptions of the system under investigation (Maloszewski and Zuber (1982)). By applying a convolution integral, it can balance the tracer inputs and outputs of equation (2), as follows: and sources of water in a single characteristic (McGuire and McDonnell (2006)). Assuming steady-state conditions, travel time distributions are often interpreted and applied as time invariant, for instance, as a mean over the period under investigation.
In addition, travel time distributions are usually inferred using lumped parameter models that simplify the description of a spatially distributed catchment behavior. Of course, there is evidence and knowledge to the contrary (Hrachowitz et al. (2010), 135 McDonnell et al. (2010), Botter (2012), Heidbüchel et al. (2013)).
For non-stationary systems, it makes sense to distinguish between two types of probability density functions: the forward travel time distribution and the backward travel time distribution (Niemi (1977)). The forward travel time distribution, − → h (ϕ, t in ), is the probability distribution of the travel times ϕ conditional on the injection time t in of a volume flow (e.g., precipitation). The backward travel time distribution, ← − h (ϕ, t), is the probability distribution of the travel times ϕ conditional on 140 the exit time t of a volume flow (e.g., discharge). When the system is in a steady state (constant input/output fluxes), then the forward and backward travel time distributions collapse into a single probability density function (Niemi (1977), Botter et al. (2011); otherwise, the following relation (Niemi's theorem) applies: where θ(t in ) is a partition function that describes the fraction of rainfall J(t in ) that ends up as runoff Q t . In addition, an age 145 function ω Q can be defined that describes the ratio between the number of water particles with an age in the interval (ϕ, ϕ+dϕ) sampled by Q t at time t and the amount of particles with the same age stored in the control volume at that time: where h( , t) is the probability distribution of the residence times of the water particles stored within the control volume at time t.

150
The age function ω Q (ϕ, t) is an interesting quantity in the sense that the tracer concentrations c e , c e-1 , . . . , c e-τ and c p , c p-1 , . . . , c p-τ from equation system (5) can be represented as a function of the same. For instance, the pre-event concentration c p for the event e can be calculated as follows: In addition, the tracer concentrations from the iterative separation model (5) are initially subject to the following assumption: Therefore, it implicitly follows that the age function ω Q (ϕ, t) should remain stable during the subsequent rainfall-runoff events e, e + 1, . . . , e + τ , which requires piece wise steady-state conditions, and that is indeed the core message of Criterion 2 and Criterion 3.
The basic procedure to reconstruct the event water response has already been demonstrated in Hoeg (2019), where the single 160 components Q e , Q e-1 . . . Q e-τ were arranged in the following way: Assume we are interested in the event water contribution of event 1 during events 2, 3, and 4. For this case, I can arrange one after the other: the event water Q e of event 1, the last event water Q e-1 of event 2, the next-to-last event water Q e-2 of event 3, and the next-to-next-to-last event water Q e-3 of event 4, as illustrated in Figure A1.
When referring to the rainfall events J 1 , J 2 , . . . , J τ +1 and the related catchment responses Q 1 , Q 2 , . . . , Q τ +1 , I can define 165 the separated event water response as follows: This is the volume flow of water that entered the catchment at interval [t 1 , t 2 [ with rainfall J 1 and that appears in the stream discharge Q during interval [t 1 , t 1+τ +1 [, as a result of the rainfall events J 1 , J 2 , . . . , J τ +1 and the related catchment responses Q 1 , Q 2 , . . . , Q τ +1 .

170
Furthermore, I postulate that the volume-weighted function of the time-varying separated event water response H [ti,ti+1[ can be considered as an approximation of the time-varying backward travel time distribution, [ regarding all water molecules that entered the catchment (the control volume of system (1) regarding the semantic time intervals e, e + 1, . . . , e + τ on τ backward iterations.
Later in the results, which are given in section 3, I systematically compare the separated event water response with the simulated event water response of the rainfall-runoff model introduced in section 2.4 for different test scenarios. Two of these (section 3.3) are based on time-varying travel time distributions. The separated event water response is related to the left hand side 180 of Niemi's theorem, that is, equation (13), whereas the simulated event water response is related to the right hand side of it.

Condition number and error estimation
For the analysis of field data, the sensitivity of the model plus the input errors of the known variables are usually included in one measure. For instance, Genereux (1998) and Uhlenbrook and Hoeg (2003) demonstrated this based on analytic expressions for the case of uncorrelated known variables and assumed uncertainties, that is, a classical Gaussian error propagation. Others, 185 for example, Kuczera and Parent (1998), Joerin et al. (2002), Weiler et al. (2003), Iorgulescu et al. (2007), and Segura et al. (2012), approximated the expected values based on designed field scenarios and the law of large numbers, which is better known as the Monte Carlo method. By contrast, the condition number allows for the focus to be solely placed on the sensitivity of the model. The iterative separation model (5) can also be represented in the form of a linear equation system: where A is a n × n coefficient matrix. It is a sparse matrix (band matrix) with entries along and around the main diagonal; x is a column vector with n unknown runoff components; b is a column vector with n entries and that is equal to zero except for the first two entries. See appendix B for the case of τ = 3.
For the benchmark tests in section 3, I am interested in a measure that reflects the sensitivity of the relative error of the 195 solution x to the changes on the right side b. In addition, the sensitivity of the solution x to changes in matrix A is required.
This measure is known as the condition of matrix A (Turing (1948), Stoer (1994), Higham (1995)). Applying the matrix norm , this is defined as follows: For the linear equation system (20), the following estimation for the relative error in x can be performed (Stoer (1994)): which means that in the case of a large condition number, small disturbances in b can lead to larger errors in x. A similar estimation can be derived for small disturbances in the coefficients of matrix A in relation to the relative error in x: ∆x For the interpretation of the results in section 3, the following properties of the condition number are useful (Deuflhard and 205 Hohmann (  In general, an equation system with a low condition number is said to be well-conditioned, while an equation system with a high condition number (much larger than one) is said to be ill-conditioned. Note that the relative error of the complete solution x, which is subject to the vector norm , is estimated, not the relative error of each unknown, which is a single entry of vector x. In addition, the condition number does not depend on the size of the input error, as is the case with the Gaussian error propagation, and for the simulations presented in section 3, no input errors of the known variables are available.

215
In section 3, I show that a well-conditioned equation system that is an analog of (5) leads to more reliable separation results in many cases. Vice versa, an ill-conditioned system at time step t might question the applicability of the iterative separation method, regardless of other factors.
Because the relative error of a single component in solution vector x is not quantified by the condition number, it should be not used for the analysis of field data. For this purpose, the Gaussian standard error is the more suitable measure. Nevertheless, 220 the condition number, which is analogous to the Gaussian standard error, can be related to the first order partial derivatives of the linear equation system (20); for instance, it can be shown that 8 https://doi.org/10.5194/hess-2021-213 Preprint. Discussion started: 10 May 2021 c Author(s) 2021. CC BY 4.0 License.

Rainfall generator and rainfall-runoff model
For the benchmark tests, I use synthetic data generated by a simple rainfall-runoff model that consists of two linear reservoirs, 225 as shown in Figure A2. A linear reservoir is one in which the storage is directly proportional to the outflow (Dooge (1959)) where the constant k, which has the dimension of time, is the average delay time (mean residence time) resulting from an inflow into the reservoir. It has the draining process (Zoch (1934)) 230 and impulse response function from which the outflow rate is obtained by using the convolution integral. I have designed the model such that a fraction η of the outflow Q U from the upper reservoir contributes directly to the streamflow Q, and the complementary fraction 1 − η recharges the lower reservoir, which contributes to the streamflow at rate Q L : A larger η value would result in an immediate model response to rainfall (flashy scenario), whereas a lower η value results in a relaxed model response to rainfall (damped scenario). The model operates at hourly time steps over a period of four years.
The upper reservoir is initialized with a storage water level S U of 0.1 mm at time t 0 . It has an initial tracer concentration δ 18 O of -12.5‰. The mean residence time (k U ) of the upper reservoir is set to five days, whereas the lower reservoir has a mean 240 residence time (k L ) of three months. The tracer mass flows work under the assumption that both reservoirs are well-mixed.
The above model set-up is actually similar to the benchmark model presented in Kirchner (2016) and Kirchner (2019), but it is of a purely stateless and linear nature. In section 3, I make use of these properties by applying the superposition principle to exactly calculate the runoff response for each single rainfall event. In the context of the current investigation, I call it the simulated event water response.

245
In addition, the model linearity allows for applying an instantaneous delayed impulse of parts of the precipitation to simulate a rapid mobilization of pre-event water, that is, for instance, Q e-1 or Q e-2 , during subsequent events. To achieve this, a fraction α of the precipitation series J is applied to the rainfall-runoff module at delayed time t − t 0 ( Figure A2). In this scenario, the stream discharge is as follows: The delayed fraction α is not a constant value; rather, it is related to the volume of the subsequent rainfall event, such that the delayed rain impulse cannot exceed a specified part (α max ) of the subsequent event: It has been shown in several studies (e.g., Sklash and Farvolden (1979), Buttle (1994), Hinton et al. (1994), Hoeg et al. (2000) Iorgulescu et al. (2007)) that catchments can store water for days, weeks, or months, but then, the catchments will release the 255 water in minutes or hours in response to rainfall inputs. The benchmark model shown in Figure A2 does not describe the complex hydrodynamic processes that lead to this catchment behavior. Nevertheless, in section 3, I show that the iterative separation model (5) is principally able to detect a delayed response of event water and correctly assign it to the corresponding prior rainfall-runoff event.

260
In the following, I begin with a controlled sequence of precipitation and tracer inputs to demonstrate the basic behavior of the rainfall-runoff model shown in Figure A2 and the ability of the iterative separation model (5) to correctly assign the model response to the corresponding model inputs, which are expressed by a sufficient match of the simulated and separated event water response. In all test scenarios, three backward iterations (τ = 3) are calculated. In this way, an almost complete reconstruction of the event water response is possible.

265
In the second section, a random sequence of precipitation and tracer inputs is considered. Hereby, a flashy and a damped scenario is examined. Finally, a delayed response from parts of the event water can be validated.
For the tracer concentrations in precipitation c J , I choose an δ 18 O interval between -20.0‰ and -5‰. Throughout each rainfall event, the values for c J are kept constant and there are no evapotranspiration processes, interception effects, altitude effects, or similar effects. In addition, the linear storages shown in Figure A2 are assumed to be well-mixed, so that applies to all simulations. The simulated δ 18 O values for the total runoff are taken as pre-event water compositions c p right before the hydrograph rises. All results are accompanied by an error estimator. Here, the condition number of the linear equation system (section 2.3) that belongs to the iterative separation model (5)  For all the test scenarios, the deviation between the total mass input and simulated output, as well as for the tracer mass input and simulated output, is below 0.5%. In the first simulation, I start with a constant isotope signature in precipitation of -10‰. Rainfall is assumed at a rate of 4 mm/h for exactly 24 hours every 1.5 months. This interval corresponds to half the mean residence time of the lower reservoir (k L ).

Controlled rainfall and 18 O input
The rainfall-runoff model operates in the damped scenario.
the precipitation (see Figure A3). Here, Criterion 1 of the iterative separation model 5 is increasingly violated. As shown in Figure A4, the condition number (cond(A)) and deviation between the simulated and separated event water response (∆) grow continuously up to high values: more than 2x10 6 for the condition number and 4x10 −3 mm/h for the deviation. The separated event water responses have values below 0.08 mm/h and appear to be smooth during the first events, but there are increasing (numerical) instabilities in the last events. This is visible in Figure A5, in which the mean deviation between the simulated event 290 water response and separated event water response is shown in relation to the components Q e , Q e-1 , Q e-2 , and Q e-3 . There is an initial mean deviation of above 10% related to the initial water storage level S U at time t 0 . Obviously, it is not possible to make an exact separation at the beginning of the simulation. Then, until event no. 10, there is a small mean deviation of about 0.01% and a final mean increase up to 10% for the last events. After event 16, the separation increasingly falls apart. Until event 32, there are absolute deviations from the simulation on a scale of 1x10 8 mm/h.

Alternating 18 O signature in precipitation
In the second simulation, I continue with an alternating isotope signal in precipitation between -20‰ and -5‰. Again, rainfall occurs at a rate of 4 mm/h for exactly 24 hours every 1.5 months. The rainfall-runoff model operates in the damped scenario.
Because of the alternating 18 O signature in precipitation, the concentration in the total runoff cannot approach that of the precipitation (see Figure A6). As mentioned, the simulated δ 18 O values for the total runoff are taken as pre-event water 300 compositions c p right before the hydrograph rises. As a result, Criterion 1 is always met quite well. As shown in Figure A7, the condition number (cond(A)) and deviation between the simulated and separated event water response (∆) immediately drop down to small values: less than 100 for the condition number and 3x10 −6 mm/h for the deviation. The separated event water responses have values below 0.08 mm/h and look smooth for all the events. This is visible also in Figure A8, in which the mean deviation between the simulated event water response and separated event water response is shown in relation to the separated 305 components Q e , Q e-1 , Q e-2 , and Q e-3 . Again, there is an initial mean deviation of above 10% related to the initial water storage level S U at time t 0 . Then, from event 2, there is a small mean deviation of less than 0.05%, and a final mean decrease down to 0.01% for the last events. The separated event water responses approach the simulated ones very well, with absolute deviations from the simulation on a scale below 3x10 −6 mm/h.

Random rainfall and 18 O input 310
In the third and fourth simulations, I continue with a random isotope signal in precipitation between -20‰ and -5‰. Rain occurs at a maximum rate of 8 mm/h and maximum of three days for every one to two months. The simulation covers a period of four years, the mean annual precipitation is 1254 mm, and the mean annual runoff is 1252 mm.

Flashy scenario
In the flashy scenario, the tracer concentrations for the discharge instantly follow the tracer concentrations in the precipitation 315 (see Figure A9). Criterion 1 is not fulfilled if the tracer concentration in precipitation is similar to the current tracer concentration in the total runoff, and this happens at exactly two time steps: t = 8493 h (event no. 9) and t = 32340 h (event no. Figure A10, the condition number (cond(A)) and deviation between the simulated and separated event water response (∆) shows increased values right after t = 8493 h (up to 1583 for the condition number, 0.056 mm/h for the deviation) and t = 32340 h (up to 3441 for the condition number, 0.053 mm/h for the deviation). The separated event water responses 320 have values below 0.08 mm/h and look smooth for all the events. Except for the events 6-9 (∆ = 1121%), and events 28-31 (∆ = 214%), the separated event water responses approach the simulated ones quite well, as shown in Figure A11, with a mean deviation between 0.03% and 16%. Very often, there is an increasing deviation for each further backward iteration, that is, ∆ e <= ∆ e-1 <= ∆ e-2 <= ∆ e-3 , with values between 1x10 −5 % and 16%.

325
In the damped scenario, the tracer concentrations in discharge less quickly follow the tracer concentrations in precipitation as in the flashy scenario, as shown in Figure A12. Again, Criterion 1 is not fulfilled if the tracer concentration in precipitation is similar to the current tracer concentration in total runoff, and this happens at exactly two time steps: t = 8493 h (event 9) and t = 32340 h (event 31). As shown in Figure A13, the condition number (cond(A)) and the deviation between the simulated and separated event water response (∆) show increased values right after t = 8493 h (up to 10920 for the condition number, (∆ = 115%), and events 28-31 (∆ = 25%), the separated event water responses approach the simulated ones better than in the flashy scenario, as shown in Figure A14, with a mean deviation between 0.003% and 1.7%. Very often, there is an increasing deviation for each further backward iteration, that is, ∆ e <= ∆ e-1 <= ∆ e-2 <= ∆ e-3 , with values between 1x10 −6 % and 1.7%.

335
By comparing 100 flashy scenarios and 100 damped scenarios (Table A1), the above observation can be confirmed. In the damped scenarios, the separated event water responses approach the simulated ones better than in the flashy scenarios, showing a mean of the median values in the deviation of 0.15 % for the damped scenario and 1.56 % for the flashy scenario.

Random rainfall and 18 O input and delayed response of event water
In the fifth and sixth simulations, I continue with a random isotope signal in precipitation between -20‰ and -5‰. Here, 340 rainfall occurs at a maximum rate of 8 mm/h and a maximum of three days for exactly every three months. The simulation covers a period of four years, the mean annual precipitation is 615 mm, and the mean annual runoff is 615 mm. The model operates in the damped scenario.
In addition, as introduced in section 2.4 and illustrated in Figure A2, a delayed response of the event water is considered, whereas the delay time corresponds exactly to the time interval between the rain events (three months), meaning that the 345 impulse of each rain event and delayed impulse of the previous event coincide exactly. With this approach, I would like to simulate a situation in which part of the pre-event water is mobilized by the current rain event.
3.3.1 Large delayed fraction: α max = 0.3 I start with a large delayed fraction of α max = 0.3, which means that a maximum of 30% of the precipitation impulse is shifted by exactly one event interval. Like in the damped scenario in section 3.2.2, the tracer concentrations in discharge less quickly 350 follow the tracer concentrations in precipitation as in the flashy scenario, as shown in Figure A15. The strongly delayed mobilization of the event water can bring the rainfall-runoff model closer to a violation of Criterion 1 and Criterion 4; this occurs, for example, during event 9, that is, between steps t = 17280 h and t = 19440 h. Here, the tracer concentrations in total runoff c t are just between that of the event water c e and the pre-event water c p . As shown in Figure A16, the condition number (cond(A)) and the deviation between the simulated event water response and separated event water response (∆) show  Figure A20 and Figure A21, with a mean deviation between 2.5% and 7.1%. The deviations of the single components ∆ e , ∆ e-1 , . . . , ∆ e-3 , 370 vary between between 0.2% and 14.7%.
By comparing 200 damped scenarios (3200 events) for a range of delayed fractions α max between 1% and 30% ( Figure   A22), the above observations are put into a more positive perspective. For at least half of the events (the median), the separated event water responses approach the simulated ones quite well for the complete range of α max values, with a maximum deviation below 14.0%, and a mean deviations below 1.1%. For all calculated deviations, the mean of medians is below 0.3% and the 375 mean of minimum values is below 1.1x10 −6 %.

Design and applicability of the rainfall-runoff model
One objective of the current study was to validate the iterative approach of the system (5) based on a simple rainfall-runoff model and synthetic-generated precipitation data. For this purpose, I have chosen a linear approach that enables a fairly accu-380 rate, reliable, and fast calculation of the simulated event water response based on the convolution integral.
Furthermore, I have designed the model, such that it can simulate critical situations that, on the one hand, affect the mathematical prerequisites when applying the iterative separation approach and, on the other hand, may play a significant role in field investigations, for instance, a violation of Criterion 1 or Criterion 4. Both situations critically affect the applicability of the iterative hydrograph separation and are often mentioned in tracer-based field studies (Klaus and McDonnell (2013)). Nev-385 ertheless, they can be detected or monitored rather easily: the former based on of the condition number and the latter based on equation (6).
A significant difference in the isotopic content of the event and pre-event water component (Criterion 1) is not given, for instance, if rain falling on the catchment with an isotope signature similar to that in the channel network, as it is simulated in section 3.2 for the flashy scenario ( Figure A9) and the damped ( Figure A12) scenario, with immediately increased values of 390 the condition number and deviation from the simulated event water response for exactly τ + 1 events ( Figure A10 and Figure   A13). In the damped scenario, the condition numbers and percentage deviations are lower than in the flashy scenario (Figure (A11) and Figure A14). The reason for this is obvious: in the damped scenario, differences such as |c e − c p | are larger than in the flashy scenario.
Because the pre-event water concentration is determined before an event, the isotope composition in the total runoff may not 395 represent a valid mixture of the event water and pre-event water component (Criterion 4) if the proportion of rapidly mobilized pre-event water is relatively large. This is simulated in section 3.3.1 for the damped ( Figure A15) scenario during event 9. Here, the large fraction (30%) of rapidly mobilized pre-event water influences the isotope signature in the total runoff so strongly that the condition of equation (6) is just met to some extent. A further increase of the delayed fraction α max would let the separated event water fraction exceed 100%, whereas the pre-event water fraction would fall below 0%, that is, show negative values.

400
If the isotope composition of the event component does not sufficiently differ from the pre-event component, a violation of Criterion 4 is accompanied by immediately increased values of the condition number and deviation from the simulated event water response for exactly τ + 1 events ( Figure A16 and Figure A17).
Conversely, the iterative separation method performs better for medium and smaller fractions of rapidly mobilized pre-event water (α max ), as shown in section 3.3.2. Here, Criterion 4 is consistently fulfilled, and the iterative separation method can 405 reconstruct the complete event water response, as shown in Figure A20 and Figure A21. Nevertheless, for at least half of the events (the median), a good match between the simulated and separated event water responses can be expected for all considered values of α max ( Figure A22). This is an important capability of the iterative equation system (5) (2018)), which would affect the assumption made in Criterion (3). Furthermore, a constant isotopic signature for each simulated event is given as requested in Criterion (2), but under field conditions, intra-storm variabilities of the isotopic composition in precipitation (Kendall and McDonnell (1998)), complex patterns of throughfall depth and water isotopic composition because of interception processes (Saxena (1986), Herbstritt et al. (2019)), along with spatial or temporal variabilities because of altitude effects (Hoeg et al. 420 (2000), Schmieder et al. (2018)) may violate the same.

Applicability of the condition number as error estimator
Initially introduced for the analysis of rounding-off errors (Turing (1948)) and the numerical solvability of linear equation systems (Mallock and Darwin (1933)), the condition of a linear equation system can be regarded as a measure of the sensitivity of the solution of a linear system to the disturbances in the data and as a measure of the sensitivity of the matrix inverse to the 425 disturbances in the matrix. The condition number bounds the level of uncertainty inherent in the solution before a numerical algorithm is applied (Stoer (1994), Higham (1995)). Hence, it helps to evaluate the mathematical constraints that arise directly from the addressed physical model, which is a hydrograph separation model.
In section 3, I have demonstrated that a well-conditioned equation system analogous to (5) leads to more reliable hydrograph separation results in general, whereas an ill-conditioned system at a certain time step might impede the applicability of 430 the iterative hydrograph separation method for the current and subsequent events, here independent of other factors such as uncertainties in the field data. As shown in equations (22) and (23), the condition number stretches an upper bound for the relative error in the solution vector, which might limit the applicability for analyzing the field data. However, because an ill-conditioned system will lead to larger gradients in a Jacobian matrix and to potentially higher errors in a Gaussian error propagation, every tracer-based hydrograph separation should be accompanied by this kind of error analysis. For the standard two component 435 hydrograph separation, this is a common, though sometimes neglected, practice (Bansah and Ali (2017), Kirchner (2019)).
Regarding the iterative hydrograph separation method, Jacobian entries for two backward iterations (τ = 2) are published in Hoeg (2019). For arbitrary n-component systems, the exact analytic expressions are sometimes hard to retrieve or compute; in these cases, a Monte Carlo analysis should be a valid fallback. case. There is currently no hydrological concept or method that can directly solve the problem shown in section 3.3.1. In fact, it is Niemi's theorem (equation (13)) that needs to be balanced here.

Event water response as an estimator for time-varying backward travel time distributions
In section 3, I have demonstrated that the event water response can be interpreted as the breakthrough of an isotope impulse because of a rainfall-runoff event over τ subsequent rainfall-runoff events. Based on equation (18), I have postulated that the volume weighted separated event water response can serve as an estimator for the backward travel time distribution over the 455 current and τ subsequent rainfall-runoff events, implying that recent water fractions, such as Q e-1 , Q e-2 . . . Q e-τ , can be used to estimate the time-varying structure of the catchment storage.
Regarding the performed flashy and damped benchmark tests (section 3.2), the above postulate could be very well and easily confirmed, provided that Criterion (1) is sufficiently met. By contrast, regarding the performed delayed fraction benchmark tests (section 3.3), the above postulate could be confirmed only for the small and medium delayed fraction scenario, provided that 460 Criterion (1) and Criterion (4) are sufficiently met. One reason is that with an increasing proportion of rapidly mobilized preevent water, such as Q e-1 , the pre-event water concentration c p may decreasingly represent the bulk pre-event concentration of the considered control volume. Please note that in all scenarios shown in section 3, the pre-event concentrations are taken from the stream discharge concentration c t right before the hydrograph rises and are assumed to be constant over the entire event.
Here, intra-storm time-varying values instead of constant concentrations would probably lead to a better match between the 465 simulation and separation. This is demonstrated in Figure A23, where the values c e (e) and c e-1 (e + 1) are inversely calculated based on the iterative equations system (5) and the simulated event water response, that is, the response of the rainfall-runoff model to a specific rainfall event, as follows: It can be seen that for strong delayed fractions of the precipitation impulse (with a maximum of 30%), the event water con-470 centrations of the first subsequent event c e-1 (e + 1) are not constant at this point, as it was assumed in the benchmark tests in section 3.3, which is in contrast to the event water concentration of the event related to the precipitation impulse c e (e). Thus, in relation to the bulk mass balance equations (1) and (2), for scenarios with a large fraction of rapidly mobilized pre-event water, it does make sense to express the event water concentration of the prior event c e-1 (e + 1) as a function of the tracer concentration in the water storage S(t) to better represent the bulk pre-event concentration of the considered control volume. For 475 the strong delayed fractions scenario based on the rainfall-runoff model in section 2.4, the simulated event water response can be separated by applying the event concentrations c e-1 (e + 1) that arise from equation (32). However, to satisfy the complete scenario, an optimization regarding all simulated event water responses would be necessary to find the intra-storm varying tracer concentrations that provide the best results.
What about field scenarios in which time dependent variables occur at all time scales and a simulated event water response Q e , Q e-1 . . . Q e-τ , and Q p-τ , to connect the separation results closer to intra-storm related catchment travel times.
In addition, the rate of change of event water Q e and recent water Q e-1 , Q e-2 . . . Q e-τ and their contributions to discharge with storage can also be characterized by the age function ω Q , as shown in equation (14). The age function describes the ratio between the number of water particles with an age in a specific time interval sampled in the stream discharge at time t and the amount of particles with the same age stored in the control volume at that time. If we allow event water to be represented by 490 water with age close to the semantic time interval e, that is, t ∈ e, then value ω Q (e, t) scales with the proportion of event water.
It is the water delivered to the catchment outlet instantaneously, without spending time in storage. On the contrary, ω Q (e−1, t) refers to the water delivered to the catchment outlet delayingly, by staying in the catchment storage over the entire event e − 1.
Then, ω Q (e − 2, t) refers to the water delivered to the catchment outlet delayingly by staying in the catchment storage over the entire events e − 2 and e − 1, and so forth. Regarding a non-stationary problem like the strong delayed fraction scenario  (9), this could be achieved by finding a correlation between (c p-τ +1 − c p-τ ) and (c e-τ − c p-τ ) for a longer observation period. Subgroups of correlations could be retrieved from this, for example, with a seasonal dependency or a relation to antecedent moisture. For the correlation function: we finally would get the approximation for Q e-τ depending on the observations in the stream only, as follows: Figure A24 shows the corresponding correlation functions for the scenarios presented in section 3.

520
A validation of the recently introduced iterative extension of the standard two-component hydrograph separation method was presented. Data were retrieved from a random rainfall generator and a rainfall-runoff model composed of linear reservoirs. The results showed that it is possible to reconstruct the simulated event water response for a given model input by using a single stable isotope tracer. It was demonstrated, how mathematical constraints, such as an ill-conditioned linear equation system, may influence a separation of the event water response. It was considered and discussed how the volume weighted separated 525 event water response can serve as an estimator for the time-varying backward travel time distribution. The benchmark model also covered the partially delayed response of event water so that a situation was simulated in which pre-event water is rapidly mobilized. Related to this, it was shown that in the case of large fractions of rapidly mobilized pre-event water, a time-varying determination of the pre-event water tracer concentrations to better represent the entire control volume becomes a key aspect.
The results obtained above should provide certainty about the capabilities of the iterative hydrograph separation method.

530
Of course, further validation scenarios with different parameter sets are conceivable. Non-linear rainfall-runoff models can be involved. In addition, further tests with experimental field data would be desirable. The effort will be worth it because the iterative hydrograph separation method is easy to use and requires no calibration effort. It provides immediate results directly from observational data with low coding and low computational effort, and it can be connected with further methods for more detailed analyses and special use cases.  The iterative separation model (5) can be represented in the form of a linear equation system for τ = 3, as follows:        Figure A17. Simulation with the damped scenario and random 18 O strongly delayed input -Mean deviation between the simulated event water response and separated event water response for each event in % -∆ e : Deviation related to Q e , ∆ e-1 : Deviation related to Q e-1 ,∆ e-2 : Deviation related to Q e-2 ,∆ e-3 : Deviation related to Q e-3 , ∆: The mean of all deviations.