2.1 Definitions

The general theoretical framework relies on the works by Botter et al. (2011), van der Velde et al. (2012), Harman (2015) and Benettin et al. (2015b). Here, we consider a typical hydrologic system with precipitation J(t) as input and evapotranspiration ET(t) and streamflow Q(t) as outputs. The total system storage is obtained as $S\left(t\right)=S_{\mathrm{0}}+V\left(t\right)$, where S₀ is the initial storage in the system and V(t) is the storage variations obtained from the hydrologic balance equation $\text{d}V/\text{d}t=J-\text{ET}-Q$.

The system state variable is the age distribution of the water storage. Indeed, at any time t, the water storage is comprised of precipitation inputs that occurred in the past and that have not left the system yet. Each of these past inputs can be associated with an age T, representing the time elapsed since its entrance into the watershed. Hence, at any time t the storage is characterized by a distribution of ages p_S(T,t). Similarly, discharge and evapotranspiration fluxes are characterized by age distributions p_Q(T,t) and p_ET(T,t), respectively. Each water parcel in storage can also be characterized by its solute concentration C_S(T,t), which in case of an ideal tracer is equal to the concentration of precipitation upon entering the catchment C_J(t−T). Tracer concentration in streamflow is indicated as C_Q. A useful, transformed version of the storage age distribution is the rank storage S_T, which is defined as $S_T(T,t)=S\left(t\right){\int }_{\mathrm{0}}^T{p}_{\text{S}}(\mathit{\tau },t)d\mathit{\tau }$ and represents the volume in storage younger than T at time t.

The key element of the formulation is the SAS function, which formalizes the functional relationship between the age distribution of the system storage and that of the outflows. Different forms have been proposed to express the SAS function directly as a function of age or as a derived distribution of the storage age distribution, (e.g., absolute, fractional or ranked SAS functions; see Harman, 2015). For numerical convenience, here SAS functions are expressed in terms of cumulative distribution functions (CDFs) of the rank storage, for both discharge (Ω_Q(S_T,t)) and evapotranspiration (Ω_ET(S_T,t)). Namely, Ω_Q(S_T,t) is, at any time t, the fraction of total discharge which is produced by S_T(T,t). Hence, it is equal to the fraction of discharge younger than T. The corresponding probability density functions are indicated as ω_Q(S_T,t) and ω_ET(S_T,t). The main model variables are illustrated in Fig. 1.

Figure 1Conceptual illustration of the main variables of the theoretical formulation. Precipitation volumes are represented through colored circles, with darker colors indicating the older precipitations with respect to current time t. Due to transport and mixing processes, precipitation volumes are retained in the catchment storage and released to streamflow (a). Both the catchment storage and its outfluxes are characterized by a distribution of ages (b, c). For example, the youngest water (age t−t₄, light blue color) accounts for 8∕20 of the storage and 3∕8 of streamflow. By cumulating such distributions one gets the rank storage S_T(T,t) and the cumulative discharge age distribution P_Q(T,t) (d, e, red lines). The relationship between S_T(T,t) and P_Q(T,t) is quantified by the SAS function Ω_Q(S_T,t) (f).

Download

2.2 The age master equation

The age ME (Botter et al., 2011) can be seen as a hydrologic balance applied to every parcel of water stored in the catchment. Two different equations can be formulated that describe the forward-in-time or the backward-in-time process (Benettin et al., 2015b; Calabrese and Porporato, 2015; Rigon et al., 2016). Here, we focus on the backward form, as it is the most convenient to model solute concentration in streamflow. The backward form of the ME can be written in a number of equivalent forms that have been proposed in the literature (e.g., Botter et al., 2011; van der Velde et al., 2012; Harman, 2015). Here, we employ the cumulative version, which has a less intuitive physical interpretation but a better suitability to numerical implementation. The complete set of equations reads as follows.

The initial condition $S_{T_{\mathrm{0}}}$ indicates some initial distribution of the rank storage at time 0. Note that to ensure that p_S, p_Q and p_ET are distributions over the age domain $(\mathrm{0},+\mathrm{\infty })$, the SAS functions must verify the condition ${\mathrm{\Omega }}_Q(S_T\to S(t),t)={\mathrm{\Omega }}_{\text{ET}}(S_T\to S(t),t)=\mathrm{1}$. This condition, however, is automatically verified as the SAS functions were defined as CDFs.

The solution of Eq. (1) gives the rank storage S_T(T,t), from which the discharge age distributions p_Q(T,t) can be obtained as

where P_Q(T,t) is the cumulative distribution of p_Q(T,t) and $P_Q(T,t)={\mathrm{\Omega }}_Q(S_T,t)$ by definition. Stream solute concentration C_Q(t) follows from

The same reasoning applies to the age distributions and concentration of the evapotranspiration flux.

2.3 The SAS functions

As explained in Sect. 2.1, SAS functions are CDFs over the finite interval [0,S(t)]. A simple class of probability distributions that is suitable to serve as an SAS function is the power-law distribution (Queloz et al., 2015b; Benettin et al., 2017b), which takes the form

The parameter $k\in (\mathrm{0},+\mathrm{\infty })$ controls the affinity of the outflow for relatively younger or older water in storage. Specifically, k<1 [k>1] implies affinity for young (old) water, whereas the case k=1 represents random sampling, i.e., outfluxes select water irrespective of its age. k can be conveniently made time variant (e.g., dependent on the system wetness) to account for possible changes in the properties of the system (see van der Velde et al., 2015; Harman, 2015). Equation (6) also requires knowledge of the initial storage in the system S₀, which can be difficult to estimate experimentally and it is often treated as a calibration parameter. When using power-law SAS functions for both Q and ET, the system only requires three calibration parameters: k_Q, k_ET and S₀. Different classes of probability distributions can be used to have more flexibility in the SAS function shape, e.g., the beta (van der Velde et al., 2012; Drever and Hrachowitz, 2017) or the gamma (Harman, 2015; Wilusz et al., 2017) distributions. Such functions can be more difficult to implement numerically, but they are usually available in software libraries.

2.4 The special case of well-mixed or random sampling

In case all the outflows remove the stored ages proportionally to their abundance, the outflow age distributions become a perfect sample (or random sample, RS) of the storage age distribution. The SAS functions in this case assume the linear form ${\mathrm{\Omega }}_Q(S_T,t)={\mathrm{\Omega }}_{\text{ET}}(S_T,t)=S_T(T,t)/S\left(t\right)$ and Eq. (1) has the analytical solution (Botter, 2012)

Equation (7) can be seen as a generalization of the linear reservoir equation to fluctuating storage. Indeed, in the special case of a stationary system, where $J=Q+\text{ET}$ and the ratio J∕S is a constant c, Eq. (7) takes the simple form $p_{\text{S}}\left(T\right)=c\exp (-cT)$.

3.1 Problem discretization

Equation (1) does not have an exact solution, except for the particular case of randomly sampled storage (Sect. 2.4). Thus in general a numerical implementation is required. Following the approach by Queloz et al. (2015b) and Harman (2015), the partial differential equation (Eq. 1) is first converted into a set of ordinary differential equations using the method of characteristics. Indeed, along a characteristic line of the type $t=T+t_{\mathrm{0}}$, Eq. (1) simplifies into an ordinary differential equation in the single variable T:

with initial conditions $S_T(\mathrm{0},t_{\mathrm{0}})=\mathrm{0}$. In this context, reformulating the problem along characteristic lines means following the variable $S_T(T,T+t_{\mathrm{0}})$, i.e., the fraction of storage younger than the water input entered in t₀. This can be equally interpreted as the amount of water storage entered after time t₀. The solution $S_T(T,T+t_{\mathrm{0}})$ starts from the value 0, corresponding to the initial time t₀. Then, as time (and age) grows, $S_T(T,T+t_{\mathrm{0}})$ increases when precipitation J introduces younger water into the system and decreases when outfluxes Q and ET withdraw water younger than T. Water entered after t₀ gradually replaces the water entered before t₀ and for very large T the solution reaches (asymptotically) the total storage in the system, as no water that had entered before t₀ is still present in the system.

We discretize time and age using the same time steps $\mathrm{\Delta }T=\mathrm{\Delta }t=h$, resulting in $T_i=i\cdot h$ and $t_j=j\cdot h$, with $i,j\in \mathbb{N}$ and we use the convention that the discrete variables T_i and t_j refer to the beginning of the time step. To simplify the notation, square brackets are used to indicate the numerical evaluation of a function and the indexes i and j are used for T_i and t_j, respectively. For example, f[i,j] indicates the numerical evaluation of function f(T_i,t_j). The conventions used for the discretization are illustrated in Fig. 2. For numerical convenience and because real-world data often represent an average over a certain time interval, all fluxes (J, Q, ET) are considered as averages over the time step h (e.g., $J\left[j\right]=\mathrm{1}/h{\int }_{jh}^{(j+\mathrm{1})h}J\left(\mathit{\tau }\right)\text{d}\mathit{\tau }$). As a consequence, storage variations obtained from a hydrologic balance are linear during a time step and each value refers to the beginning of the time step.

Figure 2Illustration of the conventions used to discretize the time domain. Time steps have a fixed length h (e.g., 12 h) and each time step j starts in $t_j=j\cdot h$. The numerical evaluation of a function f at time t_j is indicated as f[j].

Download

To solve Eq. (8), we implement a forward Euler scheme. This explicit numerical scheme is intuitive and fast to solve, and its numerical accuracy is shown to be satisfactory for many hydrologic applications (see model verification, Sect. 5.1). By terming $\mathrm{\Omega }[i,j]=\mathrm{\Omega }\left(S_T\right[i,j],t_j)$, the discretized problem becomes

for $i,j\in [\mathrm{0},N]$, with N indicating the number of time steps in the simulation, and boundary condition $S_T[\mathrm{0},j]=\mathrm{0}$. In a pure forward Euler scheme, this boundary condition implies that $\mathrm{\Omega }[\mathrm{0},j]=\mathrm{\Omega }(\mathrm{0},t_j)=\mathrm{0}$, meaning that no input can be part of an output during the same time step. This can be a limitation for catchment applications, where “event” water is often not negligible and it can bear important information on catchment form and function. For this reason, in Eq. (9) we use a modified Ω^* defined as

where e[j] is an estimate of the youngest water stored in the system at the end of time step j. Such an estimate is obtained here as $e\left[j\right]=max(\mathrm{0},J[j]-Q[j\left]{\mathrm{\Omega }}_Q\right[\mathrm{0},j-\mathrm{1}]-\text{ET}[j\left]{\mathrm{\Omega }}_{\text{ET}}\right[\mathrm{0},j-\mathrm{1}]$, i.e., it is a water balance for current precipitation input using the SAS functions evaluated at a previous time step. The classic Euler scheme is returned if e[j]=0. This modification of the classic numerical scheme only affects the behavior of the youngest age in the system and it is a simple and efficient way to account for transport of event water. The accuracy of this numerical scheme is evaluated in Sect. 5.1.

3.2 Numerical routine

The model solves Eq. (9) by implementing an external for-loop on j (i.e., on the chronological time) and an internal for-loop on i (i.e., on the ages). This means that during one time step j all the characteristic curves (Eq. 9) are updated by one time step. The internal loop is implemented using vector operations. The vector length is indicated as n_j and it depends on the number of age classes (which is also the number of characteristic curves) that are included in the computations at time j (see Sect. 3.3). At any time step, the two fundamental operations to solve the discretized ME are

• compute ${\mathrm{\Omega }}_Q^{*}[i,j]$ and ${\mathrm{\Omega }}_{\text{ET}}^{*}[i,j]$ using Eq. (10)

• compute S_T[i,j] using Eq. (9) for $i\in [\mathrm{1},n_j]$.

To compute the model output, further operations are required. In particular,

• update $C_S[i,j]=C_J[i-j]$, which is valid for conservative solutes entering through precipitation;

• compute $p_Q[i,j]\cdot h={\mathrm{\Omega }}_Q[i,j]-{\mathrm{\Omega }}_Q[i-\mathrm{1},j]$;

• compute $C_Q\left[j\right]={\sum }_{i=\mathrm{1}}^{n_j}{C}_S[i,j]\cdot p_Q[i,j]\cdot h$.

Starting from these basic routines, many additional operations can be implemented to, for example, characterize the nonconservative behavior of solutes or to compute some age distribution statistics.

3.3 Additional numerical details

A first issue that the model needs to take into account is that age distributions are defined over an age domain $[\mathrm{0},+\mathrm{\infty })$, meaning that the rank storage is made of an infinite number of elements where the oldest elements typically represent infinitesimal stored volumes. To have a finite number of elements in the computations, an arbitrary old fraction of rank storage can be considered as a single undifferentiated volume of “older” water. This allows us to merge a high number of very little residual volumes into a single old pool. Note that the term old should be used carefully as its definition depends on the particular system under consideration and it may differ depending on the characteristic timescales of the solute used to infer water age (Benettin et al., 2017a). The old pool is defined here as the volume $S_T(T,t)>S_{\text{th}}$, where S_th is a numerical parameter that can be fixed for each different application. S_th also defines the age T_th, corresponding to $S(T=T_{\text{th}},t)=S_{\text{th}}$, which indicates the oldest age that is computed individually. Numerically, the parameter S_th controls the number n_j of age classes (or equivalently rank storage volumes) that are taken into account in the computations. S_th should be chosen so that the number of elements used in the computations remains small but the numerical accuracy is not compromised. It can be convenient to define a nondimensional threshold $f_{\text{th}}\in [\mathrm{0},\mathrm{1}]$ such that S_th=f_th S(t). In this case, a value f_th=0.9 means that the old pool comprises the oldest 10 % of the water storage. When f_th=1 no old pool is taken into account. Once a storage element is merged to the old pool, its individual age and concentration properties cannot be retrieved, but the mean properties of the old pool like the mean solute concentration are preserved.

A second, connected problem regards the initial conditions of the system, i.e., the unknown storage age distribution and solute concentration to be used at the beginning of the calculations. In the absence of information, the initial storage can be considered as one single old pool; hence the initial number of age classes n₀ is equal to 1. Once computations start, new elements are introduced and accounted for in the balance, reducing the impact and the influence of the initial conditions. The old pool gets progressively smaller (and vector length n_j larger) until it reaches the stationary value defined by S_th. An initial spinup period can be used to initialize the ME balance and reduce the size of the initial old water pool. This is particularly indicated when modeling solutes with long turnover times like tritium. The influence of the initial conditions decreases with time, but given the long timescales that may characterize transport processes, it is likely never completely exhausted. This has little impact on the output concentration but it limits the maximum computable age to the time elapsed since the start of the simulation.

The computational time of a simulation can be reduced by not accounting for zero-precipitation inputs as they have no influence in the balance but increase the number of operations required at each time step. In such a case, however, the position of an element in the vector does not correspond with its age anymore and age has to be counted separately. To keep the model intuitive, we decided to not remove zero-precipitation inputs.

5.1 Model verification

Here we evaluate the numerical accuracy of the model in computing the solution of the age ME (i.e., the rank storage S_T) and streamflow concentration C_Q. The numerical model is first evaluated by comparing our modified Euler solution (Eq. 9) to a numerical implementation of the analytic solution (Eq. 1). This comparison is only possible for the case of random sampling (see Sect. 2.4), as no analytic solution is usually available for other transport schemes. Then, the comparison is made for other shapes of the SAS function, approximating the true solution with a higher-order implementation of Eq. (8). As in Sect. 4, comparisons are made on the example dataset, using daily average fluxes and the sinusoidal tracer input concentration.

Figure 5Solute concentration (C_Q) time series obtained from power-law SAS functions with parameter S₀=1000 and parameter $k\in [\mathrm{0.2},\mathrm{3.0}]$, using a 24 h time step (a). The time series are rather different, as they are progressively more lagged and damped for increasing values of k. The difference with the higher-order solution forms the residual time series (b, same scale as a). Residuals are overall limited and they do not cumulate during the 8-year simulation.

Download

For the RS comparison, the analytic solution was obtained by implementing Eq. (7) at a daily scale, considering that fluxes are piecewise constant while the storage is piecewise linear during the time step. The numerical solution for the RS was obtained by setting both Ω_Q and Ω_ET as power laws with parameters $k_Q=k_{\text{ET}}=\mathrm{1}$. The numerical model was run for eight different aggregation time steps h: 1, 2, 3, 4, 6, 8, 12 and 24 h. For each run, the resulting streamflow concentration and one rank storage (corresponding to the end of day 2745) were used for comparison with the analytic solution. Models were run for 8 years using 4 years of spinup. To allow direct comparisons across different aggregation time steps, streamflow concentrations were extracted at the end of each day, resulting in eight different time series (one per h) of 2920 elements. The time series were then normalized by the mean and standard deviation of the analytic solution. A time series of model errors on streamflow concentration (err${}_{C_Q}$) was finally obtained from the difference between the analytic and the numerical (normalized) solutions. The rank storage was evaluated on the entire age domain every 24 h (again, to allow comparisons across different time steps). To avoid comparisons between cumulative functions, the rank storage was used to compute the storage age probability density function p_S (see Sect. 2.1). The errors on p_S were obtained from the difference between the analytic and the numerical solutions. In this case, the error time series (err${}_{p_{\text{S}}}$) consists, for each of the eight aggregation time steps, of 2745 elements. For additional comparisons, the performance of our numerical implementation (EF^*) was compared to the classic implementation of the forward Euler scheme (EF, i.e., Eq. (10) with e[j]=0). Results are obtained for four different values of the initial storage S₀: 300, 500, 1000 and 2000 mm. The standard deviations of err${}_{p_{\text{S}}}$ and err${}_{C_Q}$ are shown in Fig. 4 as a function of the aggregation time step. The EF and EF^* implementations almost have the same error on p_S, indicating that accounting for the event water does not have a major impact on the overall solution of the age ME. However, as different ages do not contribute equally to streamflow, the event water can have a larger impact on streamflow concentration. This is evident in the performance on err${}_{C_Q}$, where the modified EF^* implementation is about 1 order of magnitude more accurate than the classic Euler scheme. The error is on average smaller than 10⁻² the variance of the C_Q signal, which is lower than most measurement errors. The performance on err${}_{C_Q}$ also shows that the errors tend to grow with decreasing values of the mean storage, i.e., when the storage gets depleted (or filled) faster. The error of the EF^* scheme shows a good stability. This is not surprising as the RS case resembles a linear reservoir (see Sect. 2.4) with a coefficient c approximately equal to the mean ratio between the fluxes and the storage 〈J∕S〉 during a time step. The stability condition for the Euler forward scheme in the case of a linear reservoir requires that $c<\mathrm{2}/h$ (no fast decay). In typical hydrologic applications, fluxes are usually much smaller than the storage; hence $\langle J/S\rangle \ll \mathrm{1}/h$ and the EF solution is stable.

Results show that the numerical implementation of the ME is satisfactory for the RS solution in terms of both accuracy and stability. However, solutions other than the RS case may be more challenging owing to the nonuniform age selection played by the outflows. For this reason, we tested power-law SAS functions (Eq. 6) with different values of the exponent k: 0.2, 0.3, 0.5, 0.7, 1, 1.2, 1.5, 2 and 3. The same exponent was used each time for both Ω_Q and Ω_ET. The model was run with a fixed initial storage S₀=1000, for the same timespan and aggregation time steps as in the RS case, and the performance was again evaluated in terms of err${}_{p_{\text{S}}}$ and err${}_{C_Q}$. Given the lack of analytical solutions, we approximated the true solution by using a higher-order implementation (built-in MATLAB solver “ode113”; Shampine and Reichelt, 1997) for Eq. (8). An example of the C_Q time series obtained from the different values of k for h=24 h is reported in Fig. 5. The C_Q time series are rather different, being progressively more lagged and damped for increasing values of k. Although the residual with respect to the higher-order solution can occasionally be up to 1.3 mg L⁻¹, it is on average very low compared to the signal. Thus in this case the accuracy of the model is satisfactory even for h=24 h. Note that for this dataset, the parameters of the SAS function (k=0.2 and S₀=1000) imply that 30 % of the input, on average, becomes output during the same day. The residuals are overall low and do not accumulate during the 8-year simulation, suggesting that even the 24 h simulation is stable. The performance on C_Q was further evaluated in the same way as for the RS case: we normalized the concentration signals and obtained the error time series err${}_{C_Q}$ from the difference with the higher-order solution. Similarly, we computed the errors err${}_{p_{\text{S}}}$ with respect to the higher-order solution for simulation day 2745.

The standard deviations of the errors are shown in Fig. 6 for different values of k and aggregation time steps. The errors on p_S grow for increasing preference of the SAS functions for the younger stored volumes (lower values of k). This indicates that the young water preference is a more challenging numerical condition for the solution of the age ME. This behavior is to be mostly attributed to the errors on the youngest waters in storage. Although we use a modified version of the EF scheme to account for the presence of event water in the outflows (Eq. 10), this approximation has some limitations. In particular, the youngest age in storage (e[j]) is quantified through the SAS function from previous time step; thus it may give rise to errors at the onset of intense storm events. The interpretation of the behavior of the error on C_Q (Fig. 6b) is less straightforward as the errors on the solution p_S can be amplified in various ways by the different SAS functions. Errors appear not too dissimilar for k in the range 0.5–1.2 and they are all reduced by 1 order of magnitude moving from daily to hourly time steps. The more extreme age selections (i.e., k≤0.3 and k≥2) tend to result in higher errors, although the error magnitude remains low (less than 10⁻² the signal variance) and the solution is stable.

Figure 6Numerical errors on the storage age distribution (a) and on streamflow concentration (b) as a function of the aggregation time step. The error time series are summarized through their standard deviation. Each plot shows the model performance for several shapes of the SAS function, parameterized as a power-law distribution with parameter k (Eq. 6). The color code is the same as in Fig. 5. The random-sampling case (i.e., k=1) is indicated in black and it is equivalent to the curves featuring S₀=1000 in Fig. 4.

Download

These examples suggest that the behavior of the system can be interpreted using a (nonlinear) reservoir analogy. Each individual water parcel can be seen as a depleting reservoir that decreases in time owing to the particular outflow removal (Eq. 8). This removal is mediated by the SAS functions; thus it can become large corresponding to high values of ω(S_T,t), potentially leading to an unstable fast decay. The depletion pattern of the reservoir is rather complex as it is nonlinear and it changes at every time step, but it suggests that very pronounced age selections should be considered carefully and checked for potential numerical instabilities. Note that for illustration purposes the effects of the two power-law SAS function parameters k and S₀ were presented separately (Figs. 4 and 6), but they should be considered together as lower storage values may enhance the selection of younger or older waters and increase the numerical errors. The model was tested here for several shapes of the SAS functions on a realistic hydrochemical dataset. Although every dataset is different and it would be impossible to perform a model verification valid for all applications, these results provide some first guidelines as to where the explicit numerical implementation may become critical.

5.2 Model applicability, limitations and perspectives

The model is based on a catchment-scale approach, so it only requires catchment-scale fluxes like precipitation, discharge and evapotranspiration. These fluxes can often be measured (or modeled in the case of ET) without the need for a full hydrologic model. Moreover, the pure SAS function approach implies that, differently from previous approaches (e.g., Bertuzzo et al., 2013; Benettin et al., 2015a), the transport equations which are solved in the model are completely decoupled from the way fluxes were obtained. This notably reduces the number of involved parameters and it simplifies the applicability of the model to different datasets and contexts. Although more research is needed to classify the expected shapes of the SAS functions based on measurable catchment properties, one can quickly obtain first-order evaluations of solute transport by using SAS functions already tested in the literature (e.g., van der Velde et al., 2015; Harman, 2015; Queloz et al., 2015b; Benettin et al., 2017b; Wilusz et al., 2017) and a reasonable choice of the initial storage S₀.

The use of an explicit numerical scheme has the potential of greatly reducing the computational times. Short aggregation time steps are generally recommended, especially when testing the affinity for younger storage volumes (e.g., Eq. (6) with parameter k<0.3), but in case larger time steps (e.g., h=24 h) prove satisfactory, the model can typically run in less than a second on a normal computer. The short computational times make the use of calibration techniques easier and the model structure is directly compatible with the DREAM (Vrugt et al., 2009; ter Braak and Vrugt, 2008) calibration packages. The model can be made faster by not considering the zero-precipitation times but, as explained in Sect. 3.3, this improvement is currently not implemented to keep the model more intuitive.

The model is based on a catchment-scale formulation of transport processes; thus it cannot provide spatial information unless the system is partitioned into a series of spatial compartments (e.g., Soulsby et al., 2015). Even in this case, one would need to know the fluxes to and from each compartment, hence losing one of the main advantages of the general SAS approach. The catchment-scale nature of the formulation also implies that SAS functions have a conceptual character and they cannot be determined directly from physical properties of the system. Their general shape, however, can be traced back to elementary advection–dispersion processes (Benettin et al., 2013) and the mechanistic basis for time-variable SAS functions has recently been highlighted (Pangle et al., 2017).

Although the numerical accuracy of the computations has to be evaluated for each different application, Sect. 5.1 provides some first guidelines to cases in which the numerical accuracy may not be satisfactory. Systems whose storage is quickly depleted by the fluxes are prone to inaccuracies and instabilities. This can happen, for instance, if the system storage is small compared to the fluxes and the SAS functions have a very strong preference for some storage portions. In such cases, higher-order schemes may become desirable. The model package already provides a higher-order solution to Eq. (8) (obtained through the MATLAB built-in function “ode113”) that can help evaluate the numerical accuracy of the results.

The codes implemented in the tran-SAS package can be used to simulate the transport of conservative solutes through a catchment. This represents a first step towards the modeling of large-scale solute transport. Simple reactive transport equations can be easily implemented in the main model routine (Sect. 3.2) using effective formulations that integrate biogeochemical processes across the catchment heterogeneity (Rinaldo and Marani, 1987). Being based on a travel time formulation of transport, the model is obviously not suited to simulating the circulation of solutes for which the chronology of the inputs and the age of water are irrelevant. For a number of cases of interest, however, both the time of entry into the catchment and the residence time of water within the catchment storage may play an important role in the transport process. Many such examples have been addressed in the literature using a catchment-scale approach, including the case of nitrate export from agricultural catchments (Botter et al., 2006; van der Velde et al., 2012), solutes influenced by evapoconcentration effects (Queloz et al., 2015b), pesticide transport (Bertuzzo et al., 2013; Lutz et al., 2017) and solutes produced by mineral weathering (Benettin et al., 2015a). The provided codes are designed to be easy to understand, so that they can be easily customized by the user and adapted to different contexts and applications. The next step is then to adapt the model to real-world problems, for which solutes' nonconservative behavior has to be taken into account.