Age-ranked hydrological budgets and a travel time description of catchment hydrology

Age-ranked hydrological budgets and a travel time description of catchment hydrology Riccardo Rigon1, Marialaura Bancheri1, and Timothy R. Green2 1Dipartimento di Ingegneria Civile Ambientale e Meccanica, Universitá degli Studi di Trento, Italy 2USDA-ARS, Water Management and Systems Research Unit, Fort Collins, Colorado USA Correspondence to: Rigon,R. , Dipartimento di Ingegneria Civile Ambientale e Meccanica, Universitá degli Studi di Trento, Via Mesiano, 77, 38123, ITALY. (riccardo.rigon@unitn.it)

associated compounds. First, van der Velde et al. (2012) made the SSF a function of the residence time PDFs using actual time, rather than using the "injection time". Subsequently, Harman (2015b) reformulated the SSF to be a function of the watershed storage and actual time.
These were valuable advances to the theory, but the literature remains obscured by different termi-25 nologies and notations, as well as model assumptions that are not fully explained. Thus there remains a need for theoretical developments that are clearly explained and developed using a consistent set of notations. Questions also remain about how to apply the theory of age-ranked distributions in terms of the model form and parameter estimation. Harman (2015b) noted the importance of selecting an appropriate SSF, but until very recently ( (Harman, 2015a)) there was no proposed method for 30 selecting the form of an SSF and estimating it from available data. Selection of a SSF for a given watershed remains a topic of importance, because it should not be imposed arbitrarily.
Here, we explore some complex cases of the SSF that are consistent with the theoretical advances and can be estimated from available data in some watersheds. In the following sections, the theory to date is reviewed and synthesized into a framework with consistent notation. This alone will ad-35 vance the reader's understanding of how to solve the travel-time distribution problem. The SSF is also defined within the theoretical framework, and the concepts of forward and backward PDFs are fully explored. These conceptual developments are followed by improved methods for selecting the appropriate form of a SSF and estimating its parameters. Guidance for hierarchical approaches to parameter estimation is given based on available data. Finally, the proposed framework and methods 40 are illustrated using data from experimental watersheds.

Definitions of age-ranked quantities
Residence time, travel time and life expectancy of water and associated constituents flowing through watersheds are three related quantities whose meaning is well defined by the following equation:

where T [T] ([T] means time units) is the travel time, t [T] is the actual time measured by a clock, τ
[T] is the injection time (i.e., the time in which a certain amount of water enters the control volume) and ι [T] is the exit time (i.e., the time in which a certain amount of water exits the control volume).
Based upon these definitions, T r := t − τ [T] is the so called residence time, or the age of water entered at time τ , and L e := ι − t [T] is the life expectancy of the same water molecules which are 50 inside of the control volume.
Consider, for example, a control volume as the one shown in figure 1. Its (bulk) water budget is written as: J(t) is the precipitation, usually a given (measured) quantity, while the discharge and the actual evapotranspiration, Q(t) and AE T (t), are modeled. Common simple estimates for the two latter quantities are: where λ [T ] and b are the parameters of the non-linear reservoir model, S max is the maximum water storage and E(t) is the potential ET, temporal function of the radiation inputs and atmospheric conditions. Assuming that radiation and various parameters used to model Q and AE T are given, 65 eq.(2) can be solved and S(t) obtained. If b = 1 the budget is a linear ordinary differential equation, and its solution is analytical as in Coddington and Levinson (1955); otherwise, the solution can be obtained through an appropriate numerical solver (e.g. Butcher (1987)).
Being interested in knowing the age of water we need to consider a more general set of equations.
Assume that the water storage S(t) can be decomposed in its sub-volumes s(t, τ ) [L 3 T −1 ] which 70 refer to water injected into the system at time τ . Thus: where the initial time t = 0 comes before any input into the control volume. Analogously, Q(t) [L 3 T −1 ] is the discharge out of the control volume, and q(t, τ ) [L 3 T −2 ] is the part of the discharge exiting the control volume at time t composed of water molecules that entered at time τ : Finally, let J(t) [L 3 T −1 ] denote the input to the control volume. This input can have an "age", and therefore, it can be defined All these bivariate functions of t and τ , s(t, τ ), q(t, τ ), and ae t (t, τ ) are null for t < τ and can present a derivative discontinuity at the origin (t = τ ) . Given the above definitions, we can rewrite the water budget as a set of age-ranked budget equations:

85
These equations were introduced first by van der Velde et al. (2012) and named by Harman (2015b).

Backward and forward approaches
"Backward" and "forward" are well known concepts in the study of travel time distributions. They were firstly introduced by Niemi (1977), Cornaton and Perrochet (2006), to cite few, and recently refined by Benettin (2015). Benettin (2015), in particular, related the concept of backward to the 90 travel times analysis and forward to the life expectancy analysis. However, according to us, these previous works didn't fully disclose the inner meaning of the two concepts. In fact, in our theory, the probabilities are defined as backward when they "look" in time to the history of water molecules and forward when they "look" in time till their exit from the control volume. According to the previous statements, we can define a backward travel time probability, which is conditioned to t and "looks" 95 backward to τ , and a forward travel time probability, which is conditioned to τ and "looks" forward 4 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-210, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 May 2016 c Author(s) 2016. CC-BY 3.0 License.
to t. In the same way, we can define a backward life expectancy probability, which is conditioned to ι and "looks" backward to t, and a forward life expectancy probability, which is conditioned to t and "looks" forward to ι. All these concepts will be better clarified in the following sections.

100
Based on the previous definitions, it is easy to define the pdfs of the residence time, travel time and evapotranspiration time. In particular, the pdf of residence time conditional on the actual time t, p(T r |t), can be defined as: where "≡" means equivalence, and ":=" a definition. Benettin (2015) denoted p(T r |τ ) as ← − p (T r , t) 105 but since this probability density is conditional to the actual time, standard probability notation is clear and unambiguous.
It is evident that this probability is time variant, since the integral in equation (1) that gives S(t) stops at the actual time t.
The pdf of travel time for water exiting the control volume as discharge, p Q (t−τ |t), can be defined 110 as: Eventually, the pdf of travel time for water exiting the control volume as water vapor, p ET (t − τ |t), can be defined as:

115
It is also possible to define the mean age of water for any of the two outlets, which is given by: for i ∈ {Q, E T }, which is a function of the sampling time.
After the above definitions, the age-ranked equation (9), can be rewritten as: 120 when a single "new water" injection of mass is considered, and the bulk quantities S(t), Q(t), AE T (t) are known as soon as the bulk water budget, equation (2), is solved. The travel time probabilities on the right side of (14) are not known. Consequently Botter et al. (2011) introduced a storage selection function, ω(t, τ ) [-], for each of the outputs, so that: Therefore equation (14), after the proper substitutions, becomes: Once assigned the ω(t, τ ) values on the basis of some heuristic, as in Botter et al. (2011), equation

130
(17) represents a linear ordinary differential equation and can be solved exactly as: where : and p(0|t) is the initial condition. Figures 2 and 3 show, respectively, a representation of P (T r |t), 135 as function of t (so it does not integrate to one), and of p(T r |t), obtained by considering three different injection times, named generically τ 1 , τ 2 and τ 3 , and assuming ω Q (t, τ ) = ω ET (t, τ ) = 1.
In particular, what is evident from figure 2 is that for those time intervals in which J(t) = 0, 140 and, therefore, Figure 4 shows the evolution of p(T r |t) with the actual time t, obtained for the same three injection times. The integral of the area under the curves, in this case, is not equal to 1 and the three functions shown in the figure are not pdfs.

Forward Probabilities
Consider again the budget age-ranked equation (9) in its integral form: It can be rewritten as a probability conditional to τ : having defined: and , as shown in figure 5, varies (with t), as expected, between 0 and 1 and has density: It can be observed instead that: and not assumed apriori, their complete shape is known only at t → ∞. For any finite time the actual knowledge we have, is better represented in Figure 6, which shows that the progress of the three curves P , F and G is unknown for future times. What is necessary to normalize F and G and define therefore probabilities function is the asymptotic value of the partition coefficient among the two fluxes: Then, it is easy to show that: are the forward probabilities density function of discharges and evapotranspiration, which properly normalize to 1 when integrated over t. The missing knowledge of Θ at any finite time, obviously does not prevent to know the actual state of the system, which is obtained by solving the budget 175 equation. More information and details on this partition coefficient are provided in Appendix A.

Niemi's relation
As a results of definitions made in sections (4) and (5) there exist two relations involving q(t, τ ), i.e. equations (11) and (30), and ae T (t, τ ), i.e. equations 12 and 31. Equating the correspondent two expression, it results: and: The above relations are known in literature as Niemi's relations or formulas, after Niemi (1977) and also reworked by Botter et al. (2010).

185
Dividing, for instance (32), for the total volume of water: and observing that: can be considered the marginal pdf of the injection times, and: which has the form of the well known Bayes theorem. More than being of some practical utility, this shows that the interpretation of the backward and forward probabilities as conditional ones is fully 195 consistent. On the other hands, this reveals that the joint probability of T r and t is: Following what written in section 5 there should be a working Niemi's relation for any finite time t, which does not require the knowledge of the asymptotic value Θ(τ ). This can be easily derived after having defined: and From these definitions, it is trivially: and, therefore, for discharges and for evapotranspiration.
These relations become useful when the backward probabilities are the known quantities (up to time t) and can be used to obtain the forward functions, f and g. As a by product, the SSF and the forward functions are also shown to be related, because, for instance, for the discharges is, for any 215 time t:  In the control volume, we can conceptually denote the subsets of the storage which contains the water molecules expected to exit at time ι as: Analogous to what was done before, we can observe that the quantity 225 has the structure of a probability density function once integrated over all ι-s, and it is reasonable to call it the probability density of storage-life expectancy for particles in the control volume at time t.
Based on equation (1) for any t: where indicates convolution among the probability density functions.

230
However, p ι (ι − t|t) can also be related to the forward probabilities discussed in the previous section. In fact, it can be observed that the probability of storage-life expectancy satisfies the following relation with the age-ranked forward quantities: (49) The integral spans the time interval up to t because we are considering the storage at this time. In (49) 235 the first equality says that the life-storage at time time is equal to the water injected at time τ which is expected to exit as discharge or evapotranspiration at time ι, integrated over all τ s. This integral is not effectively known, at time t, because, what is happening between time t and ι is unknown, and so the pdfs (as in Figure 6), unless they are specified from some educated guess, as made in the last section of this paper. Then, it follows: Thus, either as a convolution (i.e. as in equation (48)) or as related to forward probabilities, (i.e. as in equation (50)), the relation between the storage-life expectancy and the previously introduced backward and forward probabilities, is mediated by an integral.

Passive and reactive solutes 245
The formalism developed in section 2 to 6 is applicable, in principle to any substance, say indicated by a superscript i. Therefore we have a bulk budget equation for substance i, and age-ranked budget and which represent trivial extensions of equations (2) and (9), where for making the illustration simpler, we have neglected evapotranspiration, which will be re-introduced eventually. However, if the substance is diluted in water, it is usually treated as concentration in water (either in term of mass, moles or volume). Because we have various terms in the equations, concentrations are possibly as 255 many as the terms that appear. In this case, three: for the concentration in storage; for concentration in input; for discharges. The latter is actually the one which is usually covered by literature, since it is the one measured at the outlet of a control volume/catchment. For the solute discharge, in fact, it is usually assumed the validity of an integral expression like: with the usual interpretation of the symbols, and where the i has been dropped from the probability distribution function, assuming that a passive solute moves like water does. Dividing (56) by the water discharge, it is obtained: and, finally, applying the Niemi's formula: water and their load of substance. The bulk substance budget can instead be written as: and the missing concentration C i S (t) can be easily estimated with the help of (53) since S(t) is also known.
However, the age-ranked formalism can be used to understand a little more about the processes 280 in action. Starting from the quantities that appear in equation (52), the backward probability can be defined as: and analogous definitions (e.g. equation 11) can be given for the discharge and the inputs, such as to obtain, after the appropriate substitutions: (and usually guessed) is the SSF ω Q (t, τ ). However, (61) and (17) can be seen as two coupled equations, in p(t−τ |t) and ω Q (t, τ ) and we can conclude that the SSF cannot be arbitrarily imposed, but viceversa, derived.
From a practical point of view there could be some obstacles in the correct determination of the SSF, because, the distribution of the input of the substance can be unknown. In this case (61)  Finally, reactions that transform the substance i into another or fix it to the matrix, can be seen as a further output. If, for instance, we assume that the rate of decreasing of a substance due to chemical 305 reactions follows a first order kinetic, independent from τ , then this output (to be subtracted from, for instance, equation 61), could be: where k 1 and k 2 are suitable reaction's constants and s i eq represents an equilibrium storage. Whilst more complex type of reactions can be envisioned, this type of reaction (or sink term), being linear, 310 do no alter the essential traits of the theory described above.

An example of the other way around
With the scope to further clarify the formalism, we assume in this section that the forward pdfs introduced in the previous sections are assigned. To this scope we use the concept of linear reservoir, which has a long history in surface hydrology, Dooge (2003).

315
Let's consider first just only one outflow.
The bulk equation for the water budget of a single linear reservoir is: where it has been assumed, for simplicity, that J(t) = n τ =1 R τ , i.e. that the precipitation is accounted as a sequence of instantaneous impulses at different times τ s. It is also, by definition 320 of the linear reservoir: where λ [T] is the mean travel time in the reservoir. If this is the case, the age-ranked water budgets can be written as: where it is Equation (65), after integration over τ reduces to equation (63). By definition, it is s(t, τ ) = 0 for t < τ and the solution, for t > τ is well known as:

330
The equivalent solution, for S(t) gives: and the backward probability can be written, then as: If, and only if, R τ = const the probability simplifies, and it is time invariant, i.e. dependent only 335 on the residence time T r = t − τ . Please, notice that, in this case we did not appeal to equation (17) to estimate the backward probability but we could use directly definitions in equation (69).
Because discharge is just linearly proportional to the storage, it is easy to show that p q (t − τ |t) = p(t−τ |t) and, therefore, in this case, ω(t, τ ) = 1. This shows that the linear reservoir case, where for all injection times the mean residence time is equal (to λ), the SSF function is necessarily unitary.

340
However, a more general case, can be set if the mean residence time is a function of τ , meaning that equation (65) can be modified into: and its solution for t > τ is the same as (67), but with λ muted into λ τ . However, due to the dependence of λ τ on the injection time, the SSF is not anymore a constant, being equal to: This seems to suggest that imposing the characteristics of the pdf could completely determine the ω Q (t, τ ). Viceversa as already known, assigning ω Q (t, τ ) from some heuristic, obviously, would determine a mean residence time dependence on the injection time.
Non trivial ω(t, τ ) can also derive from assuming as a model for discharge a sequence of linear 350 reservoir, as in the so called Nash model, Dooge (2003). Without entering in details, a sequence of linear reservoirs implies that just le last reservoir maintain a linear relation between storage and outflow. Instead a nonlinear relationship exists between the whole storage and the same outflow, implying also a nonlinear SSF.
Using non-linear reservoirs does not allow to obtain semi-analytical results, but the fact suitably 355 tuning the parameters of each age-ranked equation that control the mean residence time affects the form of the SSF cannot change, as is also suggested by arguments below.
Other aspects come into play when the outputs are multiple. Expanding the previous linear case to include evapotranspiration, the bulk equation, under liner hypothesis becomes: where, the further assumption made is that the actual evapotranspiration is equal to: with a linear dependence on the soil water content, as for instance in Rodriguez-Iturbe et al. (1999).
The equations of water budget for the generations becomes then: where the bivariate dependence of ae(t, τ ) on the actual time and the injection time can be justified by arguing that, being the water of different ages not perfectly mixed in the control volume, plants roots sample water of different ages in different modes, according to spatial arrangements. Since the above equation (74) remains a linear ordinary differential equation, it is exactly solvable, and: where: and: Notably, soon as the outflows terms are expressible as a function of the storage multiplying the 375 age-ranked storage: the problem remains linear and analytically solvable. The quantity µ(t, τ ) is usually called age and mass-specific output rate, Calabrese and Porporato (2015). Solving equation (74) it is not even necessary to show that: The latter condition is regained if and only if aet(t, τ ) = aet(t), i.e. it depends only on the current time (which is a condition which requires the perfect mixing of aged waters). In fact, in case a dependence on τ remains, then, trivial algebra says that: first of all, to clarify the notation and unify concepts between previous related works. This was necessary to the obtain understanding of the theoretical framework, which was in some aspects still unclear. The theory in terms of age-ranked storages and fluxes was reworked to obtain a form of the 395 master equation, which allows an easy computation of the backward pdfs.
The relationship between the backward and forward formulation was clarified better defining and discussing the role of the partition coefficient between the two outputs, discharge and evapotranspiration. The importance of a correct estimate of the partitioning coefficient is a key point in the description of the watershed processes, as explained in the appendix. 400 Niemi's relationship was rederived usign our new definitions, obtaining the Bayes theorem. The consistency of the interpretation of the backward and forward pdfs as conditional ones was demonstrated.
To complete the theory, the life expectancy pdf was also defined in two different ways, through the convolution among the residence time and travel time pdfs and in relation with the forward pdfs.

405
Some aspects connected with the predictability of life expectancy were singled out and discussed.
The extension of the theory to any passive substance diluted in water, showed how the SSF func- Finally the abstract theory of age-ranked reservoirs was analyzed through the use of linear reservoirs, which hopefully clarifies the meaning and utility of SSFs for travel time analysis.
Appendix A: Symbols, Acronyms, and Notation A1 The partition coefficient Θ 415 Θ(τ ) has been introduced to complete the algebra of probabilities, in presence of more that one outflow. However studying it is important by itself, because it summarizes the relevant element of hydrologic fluxes partition.
The first plot in figure 7 shows a time-series of Θ(τ ) values obtained from a single injection time, using data from the Posina River generated from the simulation of the hydrological budget

A2 Reproducible research
In order interested researchers can replicate or extend our results, our codes are made available at https://github.com/geoframecomponents. Instructions for using the code can be found at: http:// geoframe.blogspot.com. All the material, with further information, is also linked at http://abouthydrology. blogspot.com/search/label/Residence%20time.