A Method for Calibrating the Transient Storage Model from the Early and Late-Time Behavior of Breakthrough Curves

Abstract

Solute transport in rivers is controlled by mixing processes that occur over a wide spectrum of spatial and temporal scales. Deviations from the classic advection–dispersion model observed in tracer test studies are known to be generated by the temporary trapping of solutes in storage zones where velocities and mixing rates are relatively small. In this work, the relation between the early and late-time behavior of solute breakthrough curves (BTCs) and the key parameters of the Transient Storage Model (TSM) is analyzed using non-asymptotic approximations of the model equations. Two main slopes are identified corresponding to the rising and decreasing limbs of the BTCs which are linked by specific relationships to transport and storage parameters. The validity of the proposed approximations is demonstrated with both synthetic and experimental data. Consistent with the TSM assumptions, the range of validity of the proposed approximations represents the limit of separability between surface dispersion and transient storage and can be expressed as a function of a nondimensional parameter. The results of this work can help environmental scientists identify solute transport and transient storage parameters and support the design of enhanced field tracer experiments.

1. Introduction

Tracer tests are a widespread technique for characterizing the fate and transport of solutes in streams and rivers. Experiments are typically conducted by injecting a tracer upstream and observing the evolution of the tracer concentration at one or more sections downstream. Reach-averaged transport parameters are then estimated by fitting a one-dimensional model to the observed breakthrough curves (BTCs).

Initial studies on solute transport in rivers were based on Taylor’s theory of longitudinal dispersion, which was found to be incapable of fully representing the skewed and long-tailed BTCs observed in tracer tests. Chatwin proposed that the longitudinal dispersion coefficient could be estimated from the rising limb of the BTCs and argued that the late-time behavior of the BTCs resulted from solute trapping in the viscous sublayer [1]. In rivers, the deviations from the advection–dispersion equation are exacerbated by the temporary storage of solutes in sediment [2,3,4], vegetated zones [5] and lateral pockets [6,7,8], characterized by small velocities and slow mixing rates compared to the transport in the main channel.

Several transport models have been proposed in the literature to represent deviations from Taylor’s dispersion theory, particularly to describe long-term retention phenomena due to exchanges between the surface water and the hyporheic zone [8]. These include the Aggregated Dead Zone model (ADZ) [9,10], general Multi-Rate Mass-Transfer (MRMT) formulations [11], Fractional advection–dispersion models [12], the STIR model [13,14], the Continuous Time Random Walk (CTRW) [15], data-driven models [16], and models representing hyporheic exchange as an exponentially attenuated diffusion process [3,17]. Further advancement in the field of tracer test and solute-transport modelling includes the use of reactive tracers, such as the resazurin-resorufin tracers, for quantifying the exchange with transient storage zones that are metabolically active [18,19,20]. These experimental techniques require the use of more sophisticated solute transport models [21,22,23] and quantification of deficiencies in the mass balance [24,25,26,27].

The Transient Storage Model (TSM) [28] has been widely used to analyze tracer test data in field and laboratory settings. The TSM combines a classical advection–dispersion equation with a first-order mass transfer into a storage zone of finite volume and has been shown to provide reasonably good approximations of tracer test breakthrough curves [28,29,30], although better fits can be obtained by assuming two storage components with different detention time [31,32,33]. A number of works have analyzed the capability of the TSM to consistently represent the dynamics of solute retention in streams [34,35,36].

Whilst these models are often capable of providing a good representation of the observed BTCs, uncertainties can arise in the interpretation of the model parameters. Nordin et al. [37], and more recently Gonzalez-Pinzon et al. [38], showed that the BTCs of tracer tests in rivers exhibit a persistent skewness that is inconsistent with the predictions of traditional 1D solute transport models. Recent studies [32,39] also observed that the transient storage parameters obtained from calibration of a 1D solute transport model with tracer test data can depend on the length of the study reach. Haggerty et al., found that distinct storage mechanisms acting at different timescales can combine to generate a power-law behavior of the tail of the concentration curves [40].

Even though more sophisticated models are available, the TSM is still widely used to investigate solute transport and exchange processes in rivers, thanks to its simplicity, accessibility, and ability to provide useful hints for understanding solute transport [41]. In the present work, a graphical approach is presented for estimating the parameters of the TSM. This is done by (i) deriving specific approximations for the initial rising limb and the final decreasing limb of the concentration curves, and (ii) identifying the conditions which allow separate determination of transport and storage parameters from the shape of the BTCs. The applicability of the newly derived approximations as a method to derive the parameters of the TSM from the independent fit of the rising and decreasing limb of a BTC is demonstrated with a numerical test and with experimental data. The method also provides insights into the limit of separability between surface dispersion and transient storage according to the Transient Storage Model.

2. Methods

In the TSM and ADZ models, the river channel is conceptually divided into two mutually interacting domains: a main flow channel with cross-sectional area $A$, in which the bulk flow occurs, and a storage domain of finite cross-sectional area $A_S$ adjacent to the main channel. The cross-sectional average concentration of solute is assumed to be homogeneous in the main channel and in the storage area and is denoted by $C_W$ and $C_S$, respectively. The coupled differential equations governing the evolution of the concentration along the river are:

$$\frac{\partial C_W}{\partial t}+U\frac{\partial C_W}{\partial x}-D_W\frac{{\partial }^2{C}_W}{\partial x^2}=-\frac{A_S}{A}\frac{\partial C_S}{\partial t}$$

(1)

$$\frac{\partial C_S}{\partial t}=\alpha \frac{A}{A_S}\left(C_W-C_S\right)$$

(2)

where $U$ is the advective velocity, $D_W$ is the longitudinal dispersion coefficient in the bulk flow, and $\alpha$ is the exchange rate of solute at the storage zone-bulk flow interface. The model above implies an exponential residence time distribution (RTD) in the dead zones, $\phi \left(t\right)=(1/T_D)e^{-\frac{t}{T_D}}$, with $T_D=\frac{1}{\alpha }\frac{A_S}{A}$ as the mean residence time in the storage zone and can be shown to be equivalent to other, more general residence time formulations under the assumption of a single storage domain with exponential RTD [13,14,40]. Following Davis et al. [42], the solution to Equations (1) and (2) for an instantaneous injection of mass $M$ of conservative tracer at $x=0$ and at $t=0$ can be written as:

$$\begin{gathered} C_W\left(x,t\right)=C_T\left(x,t\right){\mathrm{e}}^{-\alpha t} \\ +e^{-\frac{t}{T_D}}{{\displaystyle \int }}_0^t{C}_T\left(x,\tau \right)e^{-\frac{\tau }{T_D}}\sqrt{\frac{\alpha }{T_D}}{e}^{-\alpha \tau }\sqrt{\frac{\tau }{t-\tau }}\hspace{0.33em}{I}_1\left(2\sqrt{\frac{\alpha }{T_D}}\sqrt{\tau \left(t-\tau \right)}\right)d\tau \end{gathered}$$

(3)

where $I_1$ is the first-order modified Bessel function of the first kind, and

$$C_T\left(x,t\right)=\frac{M}{2A\sqrt{\pi D_{W}t}}{e}^{-\frac{{\left(x-Ut\right)}^2}{4D_{W}t}}$$

(4)

is the solution to Equation (1) when transient storage is disregarded. In the alternative formulation proposed by Davis et al. [10], referred to as the Aggregated Dead Zone (ADZ) model, the longitudinal dispersion term in Equation (1) is dropped, $D_W=0$, and the corresponding solution for an instantaneous mass injection at time $t=0$ is:

$$C_W\left(x,t\right)=\frac{M}{AU}\delta \left(t-\frac{x}{U}\right)\hspace{0.33em}{e}^{-\alpha t}+$$

(5)

$$+\frac{M}{AU}H\left(t-\frac{x}{U}\right)e^{-\frac{t-\frac{x}{U}}{T_D}}{\mathrm{e}}^{-\frac{\alpha x}{U}}\sqrt{\frac{\alpha }{T_D}}\left(\frac{\frac{x}{U}}{t-\frac{x}{U}}\right)\hspace{0.33em}{I}_1\left(2\sqrt{\frac{\alpha }{T_D}}\sqrt{\frac{x}{U}\left(t-\frac{x}{U}\right)}\right)$$

where $\delta \left(t\right)$ is the delta Dirac function and $H\left(t\right)$ is the Heaviside function.

2.1. Approximation for the Rising Limb of the BTC

Equation (3) implies that $C_T{e}^{-\alpha t}$ is the portion of mass $M$ that never entered the storage domain. As a first approximation, this term describes the rising limb of a BTC, observed at a section $x$ following an instantaneous injection of mass $M$ at $x=0$:

$$C_R^{\delta }\left(x,t\right)=\frac{M}{2A\sqrt{\pi D_{W}t}}{e}^{-\frac{{\left(x-Ut\right)}^2}{4D_{W}t}} e^{-\alpha t}$$

(6)

An approximation for this limb of the curve in case of a step injection of tracer of duration $T_S$ at $x=0$ can be derived from the convolution of $C_W$, expressed by Equation (3), with the concentration at the boundary. Neglecting the contribution of the second term of Equation (3), the approximating expression for the rising limb of the BTC is found as

$$\begin{gathered} \begin{array}{r}\hfill C_R^H\left(x,t\right)=\left\{\frac{Q{C}_T\left(x,t\right)}{M}\right\}\times \left\{C_0\hspace{0.33em}\left[H\left(t\right)-H\left(t-T_S\right)\right]\right\}={{\displaystyle \int }}_0^{T_S}\frac{C_{0}Q}{M}{C}_T\left(x,t-\tau \right)d\tau \end{array} \\ =\frac{C_0}{ 2} e^{-\frac{\alpha x}{U^3}\left(U^2-\alpha D_W\right)}\left[\mathrm{erf}\left(\frac{U^{3}t-U^{2}x+2\alpha D_{W}x}{2\sqrt{U^3{D}_{W}x}}\right)-\mathrm{erf}\left(\frac{2\alpha Dx-U^2\left(U{T}_S+x-Ut\right)}{2\sqrt{U^3{D}_{W}x}}\right)\right] \end{gathered}$$

(7)

where $C_R^H$ is the approximation of the rising limb of the BTC, and H(t) is the Heaviside function. The error function, $\mathrm{erf}\left(y\right)$, can be approximated near $y=0$ by its McLaurin series:

$$\mathrm{erf}\left(y\right)=\frac{2}{\sqrt{\pi }}{\displaystyle \sum }_{n=0}^{\infty }\frac{{\left(-1\right)}^n{y}^{2n+1}}{\left(2n+1\right)n!}$$

(8)

By limiting the series expansion to the first two terms and substituting into Equation (7) under the frozen cloud approximation used by Chatwin [1], a linear approximation is found for the rising limb at a distance x from the injection point in the case of a finite-length step injection at $x=0$:

$$C_R^H\left(x,t\right)=m\left(x\right)t+q\left(x\right)$$

(9)

where

$$m\left(x\right)=\frac{C_{0}U{e}^{-\alpha \frac{x}{U}}}{\sqrt{4\pi D_W\frac{x}{U}}}$$

(10)

and

$$q\left(x\right)=\frac{C_0{e}^{-\alpha \frac{x}{U}}}{2}\left[1-\frac{2}{\sqrt{\pi }}\left(\frac{x}{\sqrt{4D_W\frac{x}{U}}}\right)\right]$$

(11)

Equation (9) represents the rising limb of a BTC at a section $x$ resulting from a step concentration $C_0$ at $x=0$ from $t=0$ to the time $t=T_S$. This approximation applies around the point $t=x/U$, where $C_W\left(x,t\right)\simeq C_0/2$.

2.2. Approximation for the Decreasing Limb of the BTC

Equation (5) can be used to obtain an approximation for the decreasing limb $C_D^{\delta }\left(X,t\right)$ of a BTC following an instantaneous injection of mass $M$ at $x=0$. In the early stages of development, the decreasing limb is represented by Equation (5), provided that the term describing the concentration distribution in Equation (3), given by $C_T\left(x,t\right)e^{-\alpha t}$, is negligible for $t\gg 0$. The modified Bessel function of first order and first kind can be approximated by a linear equation $I_1\left(c\right)\simeq c/2$ in the range $0<c<1$ [43], or, more accurately, by the series [44]:

$$\begin{array}{r}\hfill I_{\nu }\left(c\right)\simeq \frac{c^{\nu }}{\nu \Gamma \left(\nu \right)}\left[1+\frac{c^2}{4\left(\nu +1\right)}+\frac{c^4}{32\left(\nu +1\right)\left(\nu +2\right)}\right.\left.+\frac{c^8}{6144{\prod }_{n=1}^4\left(\nu +n\right)}\right]\end{array}$$

(12)

which is valid for $c\le 12$ or $c\le \nu$. Using a second order approximation, thus limiting the series expansion of $I_1$ to the second term of Equation (12), the expression of $I_1$ in Equation (5) can be approximated as follows:

$$I_1\left(2\sqrt{\frac{\alpha }{T_D}}\sqrt{\frac{x}{U}\left(t-\frac{x}{U}\right)}\right)\simeq \sqrt{\frac{\alpha }{T_D}}\sqrt{\frac{x}{u}\left(t-\frac{x}{U}\right)}+\frac{1}{2}{\left(\frac{\alpha }{T_D}\right)}^{3/2}{\left(\frac{x}{U}\right)}^{3/2}{\left(t-\frac{x}{U}\right)}^{3/2}$$

(13)

Equation (13) is valid for $x/U<t<x/U+T_D/\left(2\alpha x/U\right)$. By substituting Equation (13) into Equation (5), the following approximating expression for the decreasing limb of the BTC is obtained:

$$\begin{array}{r}\hfill C_D^{\delta }\left(x,t\right)=\frac{M\alpha \frac{x}{U}}{T_{D}Q}{e}^{-\alpha \frac{x}{U}}{e}^{-\frac{t-\frac{ x}{U}}{T_D}}\left[1+\frac{\alpha }{2T_D}\frac{x}{U}\left(t-\frac{x}{U}\right)\right]\\ \hfill \end{array}$$

(14)

For $x/U<t<x/U+T_D/\left(2\alpha x/U\right)$, another approximation can be derived,

$$\left[1+\frac{\alpha }{2T_D}\frac{x}{U}\left(t-\frac{x}{U}\right)\right]\simeq e^{\frac{\alpha }{2T_D}\frac{x}{U}\left(t-\frac{x}{U}\right)}$$

(15)

which leads to the final approximated expression for Equation (14):

$$C_D^{\delta }\left(x,t\right)=a\left(x\right)e^{n\left(x\right)t}$$

(16)

where:

$$a\left(x\right)=\frac{M\alpha \frac{x}{U}}{T_{D}Q}{e}^{-\frac{\alpha }{2T_D}{\left(\frac{x}{U}\right)}^2-\frac{\alpha T_D-1}{T_D} \frac{x}{U}}$$

(17)

and:

$$n\left(x\right)=\frac{\alpha }{2T_D}\frac{x}{U}-\frac{1}{T_D}$$

(18)

Equation (16) is convenient because it makes it possible to fit the observed BTC at a distance x from the injection point by linear regression in semi-log scale.

Equation (5) describes the behavior of the decreasing limb of a BTC generated by an instantaneous injection of tracer when the effect of the mass transfer into the dead zones is much greater than the effect of longitudinal dispersion at late times. For a step injection of tracer at $x=0$, the expression for the decreasing limb must be suitably modified. In this case, the concentration at the injection section can be expressed as $C_0\left[H\left(t\right)-H\left(t-T_S\right)\right]$. If a solute mass $M$ is injected at $x=0$ with a constant rate over the time $T_S$ starting at $t=0$, the approximated expression of the decreasing limb can be found from the convolution integral,

$$\begin{array}{r}\hfill C_D^H\left(x,t\right)={{\displaystyle \int }}_0^{T_S}\frac{Q{C}_0}{M}a\left(x\right)e^{n\left(x\right)\left(t-\tau \right)}d\tau =\left\{\frac{Q{C}_D^{\delta }\left(x,t\right)}{M}\right\}\times \left\{C_0\left[H\left(t\right)-H\left(t-T_S\right)\right]\right\}\end{array}$$

(19)

The solution to Equation (19) is again an exponential function:

$$C_D^H\left(x,t\right)=b\left(x\right)e^{n\left(x\right)t}$$

(20)

where:

$$b\left(x\right)=\left\{-e^{-\left(\frac{\alpha }{2T_D} \frac{x}{U}-\frac{1}{T_D}\right)\left(T_S+\frac{x}{U}\right)}+e^{-\left(\frac{\alpha }{2T_D} \frac{x}{U}-\frac{1}{T_D}\right)x/U}\right\}\frac{2C_0\alpha }{\alpha \frac{x}{U}-2}\frac{x}{U}{e}^{-\alpha \frac{x}{U}}$$

(21)

Given an observed BTC at a distance x, the coefficients $m\left(x\right)$, $q\left(x\right)$, $b\left(x\right)$ and $n\left(x\right)$ can be determined by performing a linear regression of the rising limb on a linear scale and a linear regression of the decreasing limb on a semi-log scale. Equations (10), (11), (18), and (21) then represent a system of nonlinear equations that can be solved numerically to obtain the unknown transport and storage parameters, $U$, $D_W$, $\alpha$, and $T_D$.

To determine the limit of applicability of the proposed approximation, it is convenient to introduce the non-dimensional time, defined as $t^{*}=tU/x$. Using these non-dimensional quantities, Equation (13) can be written as:

$$\begin{array}{r}\hfill I_1\left(2\sqrt{\frac{\alpha }{T_D}}\frac{x}{U}\sqrt{\left(t^{*}-1\right)}\right)\simeq \sqrt{\frac{\alpha }{T_D}}\frac{x}{U}\sqrt{\left(t^{*}-1\right)}+\frac{1}{2}{\left(\frac{\alpha }{T_D}\right)}^{3/2}{\left(\frac{x}{U}\right)}^3{\left(t^{*}-1\right)}^{3/2}\\ \hfill \end{array}$$

(22)

valid for $1<t^{*}<t_{LIM}^{*}$, where $t_{LIM}^{*}$ is a non-dimensional parameter defined as:

$$t_{LIM}^{*}=1+\frac{T_D{U}^2}{2\alpha x^2}$$

(23)

The parameter $t_{LIM}^{*}$ represents the upper time limit for which Equation (22) approximates the modified Bessel function. The comparison between the approximated and the original Bessel function is reported in Figure 1 for four different Bessel functions, characterized by different values of the parameter $t_{LIM}^{*}$. The lower limit, $t^{*}=1$, corresponds to the non-dimensional advective time $t_{AD}^{*}=t_{AD}/\left(x/U\right)$, where $t_{AD}=x/U$. For the complete list of the symbols used in the equations, refer to

Table A1. Notation list of all the symbols used in the manuscript.

Symbol	Unit	Description
$a$	(kg m−3)	$\mathrm{coefficient} \mathrm{of} \mathrm{the} \mathrm{trendline} \mathrm{for} C_D^{\delta }$
$A$	(m2)	mean flow area
$A_S$	(m2)	transient storage area
$b, b^{*}$	(kg m−3), (-)	$\mathrm{coefficient} \mathrm{and} \mathrm{normalized} \mathrm{coefficient} \mathrm{of} \mathrm{the} \mathrm{trendline} \mathrm{for} C_D^H$
$C_0$	(kg m−3)	concentration at the injection section
$C_R^H$	(kg m−3)	approximation for the rising limb of the BTC
$C_R^{*}$	(-)	dimensionless linear approximation of the rising limb
$C_D^H$$, C_D^{\delta }$	(kg m−3)	approximation for the decreasing limb of the BTC
$C_D^{*}$	(-)	dimensionless exponential approximation of the decreasing limb
$C_S$	(kg m−3)	concentration in the storage area
$C_T$	(kg m−3)	elementary solution of the advection–dispersion equation
$C_W$	(kg m−3)	concentration in the main flow channel
$D_W$	(m2 s−1)	longitudinal dispersion coefficient
$H$	(-)	Heaviside function
$I_1$	(-)	modified Bessel function of the first order and first kind
$m, m^{*}$	(kg m−3 s−1), (s−1)	$\mathrm{slope} \mathrm{and} \mathrm{normalized} \mathrm{slope} \mathrm{of} \mathrm{the} \mathrm{trendline} \mathrm{for} C_R^H$$\mathrm{and} C_R^{*}$
$n$	(s−1)	$\mathrm{exponent} \mathrm{of} \mathrm{the} \mathrm{trendline} \mathrm{for} C_D^H$$\mathrm{and} C_D^{*}$
$q, q^{*}$	(kg m−3), (-)	$\mathrm{intercept} \mathrm{and} \mathrm{normalized} \mathrm{intercept} \mathrm{of} \mathrm{the} \mathrm{trendline} \mathrm{for} C_R^H$$\mathrm{and} \mathrm{for} C_R^{*}$
$Q$	(m3 s−1)	discharge
$t$	(s)	time
$t_{AD}$	(s)	advective time
$T_D$	(s)	mean residence time in the storage area
$T_S$	(s)	time duration of a plateau injection
$t_{LIM}^{*}$	(-)	dimensionless time limit
$U$	(m s−1)	flow velocity
$x$	(m)	longitudinal distance from injection point
$\alpha$	(s−1)	exchange rate
$\delta$	(s−1)	Dirac delta function
$\tau$	(s−1)	dummy variable for time
$\phi$	(s−1)	exponential RTD in the dead zones

in Appendix A.

3. Application

A numerical test was performed to validate the approximations derived in the previous Section 2.1 and Section 2.2. Experimental data were used for calibrating the unknown TSM parameters using the proposed graphical method.

3.1. Numerical Experiment

The validity of the approximations derived above is tested on an ideal case, by (i) simulating a BTC from a TSM model using previously fixed $A$, $A_S$, $D_W$, and $\alpha$, and (ii) plotting the approximation lines obtained from the Equations (10), (11), (18) and (21). The flow discharge is assumed to be $Q=0.4$ m³ s⁻¹ with stream width $b=5.0$ and depth $d=0.4$ m, respectively. The dispersion coefficient is estimated using Fischer’s formula, $D_W=0.011U^2{b}^2/\left(U^{*}d\right)=0.735$ m²s⁻¹ [45], where the shear velocity is $U^{*}=\sqrt{gdj}=0.037$ m s⁻¹ and $j$ is the energy slope under the assumption of uniform flow and given a Manning coefficient $n=0.05$ m^−1/3s⁻¹. The model output is generated using a single compartment TSM with storage area $A_S=0.1$ m², and exchange rate $\alpha ={10}^{-4}$ s⁻¹ at three cross-sections located at $X_1=500$ m, $X_2=1000$ m and $X_3=1500$ m from the injection point. The boundary condition at $x=0$ is a constant concentration $C_0$ injection of a solute mass $M=192$ g, applied for a period $T_S=480$ s starting at $t=0$ s.

In Figure 2, the dimensionless output concentration $C^{*}\left(x,t\right)=C_W\left(x,t\right)/C_p\left(x,t_p\right)$, where $C_p\left(x,t_p\right)$ is the peak concentration at time $t_p$, is plotted at the three sections mentioned above in both semi-log and linear scales. The approximations for the rising and the decreasing parts of the BTC are plotted according to Equations (9) and (20) in linear and semi-log scale, respectively, where parameters $m^{*}\left(x\right)=m\left(x\right)/C_p\left(x,t_p\right)$, $q^{*}\left(x\right)=q\left(x\right)/C_p\left(x,t_p\right)$, $b^{*}\left(x\right)=b\left(x\right)/C_p\left(x,t_p\right)$ and $n\left(x\right)$ are evaluated from Equations (10), (11), (18) and (21) using the same parameters used for the TSM. The values of $A$, $A_S$, $D_W$, $\alpha$, $m^{*}$, $q^{*}$, $b^{*}$, $n$ and $t_{LIM}^{*}$ are reported in

Table 1. TSM parameters $A$, $A_S$, $D_W$, and $\alpha$ are used for the concentration curves in Figure 2, coefficients $m^{*}$, $q^{*}$, $b^{*}$, and $n$ of the approximation trendlines for the decreasing and rising parts at the three output sections $X_i$, and peak concentration $C_p$ at each section used for normalizing the BTC in Figure 2.

Section	Distance	${\mathit{D}}_{\mathit{W}}$	$\mathit{A}$	${\mathit{A}}_{\mathit{s}}$	$\mathit{\alpha }$	${\mathit{t}}_{\mathit{L}\mathit{I}\mathit{M}}^{*}$	$\mathit{n}$	${\mathit{b}}^{*}$	${\mathit{m}}^{*}$	${\mathit{q}}^{*}$	${\mathit{C}}_{\mathit{p}}$
	(m)	(m2 s−1)	(m2)	(m2)	(s−1)	(-)	(s−1)	(-)	(s−1)	(-)	(g m−3)
$X_1$	500	0.735	2.00	1.0 × 10−1	1.0 × 10−4	1.400	−1.75 × 10−3	4.63 × 101	2.04 × 10−3	−4.32	0.502
$X_2$	1000					1.100	−1.50 × 10−3	2.22 × 103	1.63 × 10−3	−7.27	0.347
$X_3$	1500					1.044	−1.25 × 10−3	2.00 × 104	1.31 × 10−3	−8.96	0.274

for each section.

3.2. Field Experiment

The applicability of the model for evaluating unknown TSM parameters from measured BTCs is demonstrated on experimental data obtained from a real tracer test. The data are taken from the experimental work of Bottacin-Busolin et al. [32], conducted on the Desturo canal, an irrigation channel located in Northern Italy. The studied reach is divided into two sections located at $X_1=262$ and $X_2=567$ m, respectively, from the injection point. The plateau injection of a mass of tracer $M=3.24$ g is $T_S=1470$ s long and the flow discharge is $Q=0.045$ m³ s⁻¹.

At both section $X_1$ and $X_2$, a linear and an exponential trendline are applied to the leading and the decreasing limb, respectively, considering the portion of data in the leading and decreasing limb which maximizes the coefficient of determination R². The fitting lines are reported in Figure 3, showing the observed concentration data, normalized by the peak concentration, and the fitting trendlines in linear scale for the rising limb (panels a and b) and semi-log scale for the decreasing limb (panels c and d). The fitting lines provide the values of $m^{*}\left(x\right)$ and $q^{*}\left(x\right)$ of the linear trendline for the rising part, $C_R^{*}\left(x,t\right)=m^{*}\left(x\right)t+q^{*}\left(x\right)$, and the values of $b^{*}\left(x\right)$ and $n$ of the exponential trendline for the decreasing limb of the curve, $C_D^{*}\left(x,t\right)=b^{*}\left(x\right)e^{n\left(x\right)t}$. The unknown transport parameters $A$, $A_S$, $D_W$ and $\alpha$ are obtained by solving Equations (10), (11), (18), and (21) as a system of equations. Note that, here, the flow rate Q is known and therefore the cross-sectional velocity U can be determined from the cross-sectional area A as U = Q/A, whereas the parameter $T_D$ in Equations (18) and (21) depends on α, A, and A_S, according to the relationship $T_D=\frac{1}{\alpha }\frac{A_S}{A} .$ The resulting transport parameters and $t_{LIM}^{*}$ are reported in the first two rows of

Table 2. Coefficients $m^{*}$ , $q^{*}$, $b^{*}$, $n$ of the approximation trendlines in Figure 3; TSM parameters $A$, $A_S$, $D_W$, $\alpha$ obtained from the trendline approximations (“Approx.”) and from a numerical optimized fit (“Optim.”) of the BTCs; estimation of $t_{LIM}^{*}$ for each section and type of fitting, and peak concentration $C_p$ at each section used for normalizing the BTC in Figure 3.

Section	Distance	Fit	$\mathit{n}$	${\mathit{b}}^{*}$	${\mathit{m}}^{*}$	${\mathit{q}}^{*}$	${\mathit{D}}_{\mathit{W}}$	$\mathit{A}$	${\mathit{A}}_{\mathit{s}}$	$\mathit{\alpha }$	${\mathit{t}}_{\mathit{L}\mathit{I}\mathit{M}}^{*}$	${\mathit{C}}_{\mathit{p}}$
	(m)		(s−1)	(-)	(s−1)	(-)	(m2 s−1)	(m2)	(m2)	(s−1)	(-)	(mg m−3)
$X_1$	262	Approx.	−2.11 × 10−3	98.5	2.21 × 10−3	−2.278	0.279	0.202	2.93 × 10−2	4.00 × 10−4	1.33	48.12
		Optim.	-	-	-	-	0.243	0.200	3.72 × 10−2	5.65 × 10−4	1.21
$X_2$	567	Approx.	−1.80 × 10−3	4252.9	8.86 × 10−4	−2.582	0.090	0.244	1.86 × 10−2	4.55 × 10−4	1.02	42.8
		Optim.	-	-	-	-	0.535	0.270	1.86 × 10−2	1.37 × 10−4	1.16

.

To check the quality of the fitting, the obtained transport parameters are then used to generate a BTC using the TSM. A comparison between the simulated BTCs (dashed-dot blue line) and the observed data (circles) is reported in Figure 3 in linear (panel a and b) and semi-log scales (panel c and d).

The observed BTCs are then numerically fitted with the TSM by using differential evolution as an optimization algorithm. The computed parameters $A$, $A_S$, $D_W$, and $\alpha$ for the optimized fitted BTC are reported in Table 2 for comparison with the parameters obtained using the proposed approximations. The best fitting BTCs are reported in Figure 3 (dashed red line).

4. Results and Discussion

The results of the application presented in Figure 2 show that the goodness of the graphical estimation of the transient storage parameters and longitudinal dispersion coefficient varies at the three sections, suggesting that the method can be applied under certain conditions. In particular, it can be noted that the approximations are good at the closest section $X_1$, while at $X_2$, they start to deviate from the rising and decreasing part of the BTC. At the farthest section, $X_3$, the decreasing part of the BTC is not well fitted by the trendline obtained via Equations (19) and (21). This indicates that this curve does not behave as the exponential function at late time of Equation (20). Similarly, the rising part is no longer described by Equation (9). There are two reasons for this behavior: firstly, the term $C_T=C_W{e}^{-\alpha t}$ in Equation (4) becomes negligible at longer times and for higher exchange rates $\alpha$; secondly, at late times, Equation (13) is not a good approximation of the Bessel function. At late times and large distances from the injection point, the contribution of transient storage to the mixing process becomes more important than longitudinal dispersion, and they become two overlapping processes. This means that the rising limb of the curve no longer represents the pure effect of longitudinal dispersion. At the same time, the decreasing limb of the curve does not behave as an exponential function (Equation (13)) anymore. The dimensionless upper time limit of validity for Equation (22), $t_{LIM}^{*}$, can represent a numerical index for identifying the conditions of applicability of the approximation method proposed in this work. For sections $X_1=500$ m, $X_2=1000$ m, and $X_3=1500$ m, the values of $t_{LIM}^{*}$ are 1.40, 1.10, and 1.04, respectively. These values suggest that transient storage does not produce an exponential signature on the decreasing limb when the upper time limit of validity of Equation (22), $t_{LIM}^{*}$, is close to 1. Although there is no definite value for $t_{LIM}^{*}$, the results suggest that, in general, the exponential signature is not easily distinguishable from the longitudinal dispersion effect when $t_{LIM}^{*}<1.1$.

The application of this approach to the experimental data confirmed its validity for $t_{LIM}^{*}>1.1$. The results show that when the approximations for the leading and decreasing limbs of the observed BTCs are applied, the BTC generated by the TSM with the graphically estimated parameters for section $X_1=262$ m (dashed-dotted lines in Figure 3a,c) optimally fits the observed data. At this section, the upper nondimensional time limit of validity for Equation (22) is $t_{LIM}^{*}=1.33>1.10$. The set of parameters obtained with the proposed graphical method can provide a good first approximation of the parameters obtained by the numerically calibrated parameters. Their relative percentage errors are: +15% for $D_W$, +1% for $A$, −21% for $A_S$, and −29% for $\alpha$. For the mean residence time in the storage zone, $T_D$, evaluated according to Equation (3), the relative error between the graphical parameter estimation and the numerically calibrated model is +10%. In contrast, at section $X_2=567$ m, where $t_{LIM}^{*}=1.02$, the graphical approximation does not lead to a well-fitted BTC (dashed-dotted lines in Figure 3a,b): the time of the peak concentration is poorly caught and both the leading and decreasing limbs are not well approximated. At this section, the model calibration based on a global optimization algorithm yields a very different set of parameters, as shown in Table 2. The relative errors between the graphically obtained and the numerically optimized parameters are −89% for $D_W$, −9.5% for $A$, −0% for $A_S$, −230% for $\alpha$ and −67% for $T_D$.

The graphical method presented in this work thus provides a simple preliminary way to evaluate transient storage parameters in field experiments. The empirical rule $t_{LIM}^{*}>1.1$ can be also expressed in terms of the advective time $t_{AD}$ according to the following relations:

$$t_{LIM}^{*}=1+\frac{T_D{U}^2}{2\alpha X^2}=1+\frac{A_s}{2{\alpha }^{2}A}\frac{1}{t_{AD}^2}>{1.1}^{ }$$

(24)

$$t_{AD}<\sqrt{\frac{1}{0.1}\frac{A_s}{2A{\alpha }^2}}$$

(25)

which are more likely satisfied when the BTC is observed at relatively short distances from the injection.

5. Conclusions

Approximations were derived for the early and late time behavior of tracer test breakthrough curves for both slug and plateau injections. The approximations assumed that the transport of a solute in a river can be described by an advection-dispersion first-order mass transfer equation. Consistent with the model, the tracer BTCs measured at a section downstream of the injection point exhibited two main slopes, associated with the initial rising part of the concentration curve and the final decreasing part after the concentration peak. The parameters of the model could be inferred directly from these two slopes using the analytical relations presented in this work. The results showed that the graphical estimation of transient storage parameters and longitudinal dispersion coefficient could be made if two conditions were satisfied: first, when the two mechanisms did not produce an overlapping effect on the rising limb of the BTC and, secondly, when the condition of approximation of the decreasing limb as an exponential function was satisfied. However, a specific condition applies to the distance from the injection point for the validity of these approximations, which represents a limit of separability between longitudinal dispersion and first-order mass transfer in immobile domains. As a general empirical rule, this limit can be expressed by the condition $t_{LIM}^{*}<1.1$, where $t_{LIM}^{*}$ is a nondimensional time given by Equation (23). When this condition is not satisfied, longitudinal dispersion and transient storage become two undistinguishable mixing processes that cannot be uniquely separated by fitting the model to a tracer BTC. The condition above also provides a lower limit for the timescales of solute retention that can be inversely modelled from a tracer experiment. This limit should therefore be considered in the design of tracer tests and subsequent analyses to provide sensitive parameter estimates and identify sources of uncertainty. The results of this work provide a twofold contribution to current modelling and experimental approaches for field tracer studies: (1) by providing a simple method for estimating the surface and subsurface storage parameters; and (2) by identifying a temporal limit for the separation of surface and transient storage processes which can be used for the design of improved field tracer studies.