2. METHODS
Candidate RTDs were evaluated for different streambed depths in three steps. First, a particle tracking technique was used to generate empirical RTDs associated with bedform pumping through a hyporheic zone of various depths. Second, each of the four candidate analytical distributions (EXP, GAM, LN, and FR) were fit to the empirical RTDs, and parameter sets inferred. Finally, for each condition used to generate the empirical RTDs the four candidate analytical distributions were ranked relative to goodness of fit. Details for these three steps are described next.
2.1 Numerical Generation of the Empirical RTDs
To mimic the advective flow field associated with hyporheic exchange through stationary bedforms we adopted the analytical two-dimensional laminar flow model published by Packman et al. (2000), which is based on earlier analytical solutions of hyporheic exchange through streambeds by Elliott and Brooks (1997) and Vaux (1968). These models posit a sinusoidal pressure variation over the sediment-water interface (mimicking the static and dynamic pressure variations that develop on the surface of streambeds in a turbulent overlying flow (Cardenas et al., 2008), isotropic and homogeneous hydraulic conductivity, constant sediment porosity and fluid density, a so-called “Toth domain” for the upper boundary (Frei et al., 2019; Tóth, 1962), and an impermeable lower boundary at depth \(d_{b}\) below the surface:
\(u^{*}\)= -cos( \(x^{*}\)) [tanh( \(\text{db}^{*}\)) sinh( \(y^{*}\))+cosh( \(y^{*}\))], (1)
\(v^{*}\)= -sin( \(x^{*}\)) [tanh( \(\text{db}^{*}\)) cosh( \(y^{*}\))+sinh( \(y^{*}\))], (2)
\(h_{m}=0.28\ \frac{U^{2}}{2g}\ \left\{\par \begin{matrix}\left(\frac{\frac{H}{d}}{0.34}\right)^{\frac{3}{8}}\text{\ \ }\frac{H}{d}\ \leq 0.34\\ \left(\frac{\frac{H}{d}}{0.34}\right)^{\frac{3}{2}}\text{\ \ \ \ }\frac{H}{d}\ \geq 0.34\\ \end{matrix}\right.\ \) (3)
\(u_{m}=k\ K_{c}h_{m}\ \tan h(\text{db}^{*})\ \) (4)
\(u^{*}=\frac{u}{u_{m}}\), (5)
\(v^{*}=\frac{v}{u_{m}},\) (6)
\(k=\frac{2\pi}{\lambda}\), (7)
\(x^{*}=k\ X\), (8)
\(y^{*}=k\ Y,\) (9)
\(d_{b}^{*}=k\ d_{b},\) (10)
\(t^{*}=\frac{ku_{m}t}{},\) (11)
here, \(u\) and \(v\) are the horizontal and vertical Darcy fluxes, respectively, \(H\) is the bedform height, \(d\) is the stream depth,\(K_{c}\) is the hydraulic conductivity, \(u_{m}\) is the maximum Darcy flux at the bed surface, \(X\) and \(Y\) are the horizontal and vertical coordinates,\(\ d_{b}^{*}\) is the relative sediment depth, \(\lambda\)is the bedform wavelength, and \(k\) is the wavenumber of the bedforms. The vertical coordinate, \(Y\), is centred at the SWI (\(Y=0\)) and oriented upward; i.e., depth into the bed corresponds to negative values of \(Y\). The range of \(d_{b}^{*}\) was chosen between \(0.1\ \)and\(20\) (Supporting Information) consistent with published experimental studies performed with dune-like bedforms (Elliott & Brooks, 1997; Marion et al., 2002; Packman et al., 2000b, 2004; Packman & MacKay, 2003; Rehg et al., 2005; Ren & Packman, 2004).
The volumetric water flux predicted by equation (\(2\)) varies sinusoidally with horizontal distance along the SWI, forming well-defined upwelling and downwelling regions that are fully characterized by a repeating unit cell, one of which occurs over the domain, \(x^{*}\in\ [-\pi,\pi]\) (see figure 1b in Grant et al. (2020)). Thus, the RTD associated with Packman et al.’s hyporheic exchange flow field can be fully characterized by tracking the arrival times of particles released in the downwelling zone of a single unit cell. Accordingly, we released \(10000\) particles in the downwelling zone (\(0\leq\ x^{*}\leq\frac{\pi}{2})\) of the unit cell centered on \(x^{*}=0\), with a flux weighting scheme that added particles in proportion to the local downwelling flux (to assure that roughly the same number of particles entered the hyporheic zone along every streamline). The RTD for different choices of \(d_{b}^{*}\) in Packman et al.’s flow model was obtained by fixing the length of each particle step within the sediment domain (\({s}^{*}=ks)\) to be\({5*10}^{-3}\). Then, the i -th time step (\({t}^{*}\)) was calculated as:
\({t}_{i}^{*}=\frac{{s}^{*}}{\sqrt{{u_{i}^{*}}^{2}+{v_{i}^{*}}^{2}}}\), (12)
where \(u_{i}^{*}\) and \(v_{i}^{*}\) denote the velocity components at the particle location at the end of i -th step. For each particle the time at the end of the i -th step (\(t_{i}^{*}\)) within the sediment bed domain (\(y^{*}<0)\) is:
\(t_{i}^{*}=t_{i-1}^{*}+{t}_{i}^{*},\ \ \ i=1,\ldots,\ N\), (13)
where \(t_{0}^{*}=0\), N is the number of steps undertaken by each particle and the corresponding horizontal and vertical particle displacements (\({x}_{i}^{*}\ \)and \({y}_{i}^{*}\), respectively) of the i -th step are:
\({x}_{i}^{*}=\ u_{i}^{*}\frac{{t}_{i}^{*}}{}\) (14)
\({y}_{i}^{*}=\ v_{i}^{*}\frac{{t}_{i}^{*}}{}\ \)(15)
where is the sediment effective porosity. For each \(d_{b}^{*},\)cumulative distribution function (CDF) and probability density function (PDF) forms of the RTD were calculated from the observed residence times of the 10000 particles.
2.2 Analytical Distribution Parameter Inference.
Separate particle tracking RTDs were generated for \(125\) dimensionless bed depths ranging from \(d_{b}^{*}=0.1\) to \(20\). Each of these RTDs was fit to the four analytical distributions (EXP, GAM , LN, and FR, see Table \(1\)) described earlier using Maximum Likelihood Estimation (Mathematica, Wolfram). EXP is characterized by a decreasing monotonic function and parametrized by a single parameter (\(\rho\)). The parameters for the GAM distribution, \(\alpha\) and \(\beta\), control the shape and scale characteristics of the distribution, respectively. The parameters of LN, \(\mu\) and σ, determine the mean of the log-transformed random variable and its standard deviation, respectively. The FR distribution is a three-parameter distribution (shape parameter \(s\), scale parameter \(q\), and location parameter\(m\)), but for the sake of parsimony and for consistency with the other distributions considered in this study, the FR distribution with two parameters was chosen by fixing \(s=1\), as assumed by Grant et al. (2020). To determine which of these distributions best represents the empirical RTD at each dimensionless depth \(d_{b}^{*}\), the candidate analytical distributions were ranked using the Kolmogorov-Smirnov (KS) test, which measures the difference in the CDF between the particle tracking RTD and the assumed distribution:
\begin{equation} KS=\sup_{t}\left|\hat{F}\left(t^{*}\right)-F\left(t^{*}\right)\right|\text{\ \ }\left(16\right)\nonumber \\ \end{equation}
Where \(\sup_{t}\) is the supremum over \(t\),\(\hat{F}\left(t^{*}\right)\) is the CDF form of the empirical RTD for a certain \(d_{b}^{*}\), and \(F\left(t^{*}\right)\) is the CDF form of the candidate analytical distribution. A set of regression formulae for each distribution was then prepared to correlate the analytical distribution parameters with \(d_{b}^{*}\), using the “Curve Fitting” tool in Matlab.