The Role of Anomalous Transport in Long‐Term, Stream Water Chemistry Variability

We investigate the occurrence of anomalous (non‐Fickian) transport in an hydrological catchment system at kilometer scales and over a 36‐year period. Using spectral analysis, we examine the fluctuation scaling of long‐term time series measurements of a natural passive tracer (chloride), for rainfall and runoff. The scaling behavior can be described by a continuous time random walk (CTRW) based on a power‐law distribution of transition times, which indicates two distinct power‐law regimes in the distribution of overall travel times in the catchment. The CTRW provides a framework for assessing anomalous transport in catchments and its implications for water quality fluctuations.

DENTZ ET AL. 10.1029/2023GL104207 2 of 8 water and different fractions of pre-event water, which reside for different times within the catchment subsurface (Hrachowitz et al., 2016). Stable isotopes of H and O in water, and passive chemical tracers, are used to infer travel time distributions of streamflow (Harman, 2015;Rodriguez et al., 2021).
Current approaches to interpret this behavior involve inverse estimation of travel time distributions by least-squares deconvolution (Kirchner, 2019), fitting of forward convolution models, or assignment of empirical age ranked storage as a state variable in combination with "StorAgeSelection (SAS)" functions for streamflow and evapotranspiration, to infer their respective travel time distributions (Harman, 2015;Rinaldo et al., 2015)). Related studies interpret catchment travel times empirically in terms of a single or several gamma distributions (Hrachowitz et al., 2010;Kirchner et al., 2000;Rodriguez & Klaus, 2019;Rodriguez et al., 2021), or a beta distribution (van der Velde et al., 2012).
The purpose of this Letter is to study anomalous chemical transport through a catchment system as reflected in streamflow chemistry data, over both short and long times. To this end, we analyze a unique 36-year data set and extract, for the first time, an indication of CTRW-like anomalous behavior over extensive length and time scales. In this context, we provide a physical picture and quantitative conceptualization of chemical transport through catchments.

Background and Mathematical Development
In the following, we briefly recount relevant background on the conceptual model and spectral representation. We then develop the full CTRW model in the context of catchments which is applied to a unique chemical tracer data set from the Plynlimon experimental catchments in Wales in the next section.

Conceptual Model and Spectral Representation
For clarity, we define a hydrological catchment as a three-dimensional control volume comprising a sloping soil layer, with local undulations, underlain by a (usually porous and/or fractured) rock layer. Within this volume, the streamlines generally converge into a river network, hence ideally, the entire set of surface and subsurface runoff components feed an outflowing stream, with the catchment being treated as an input-output system. Rainfall may either run off the surface (Zehe & Blöschl, 2004), or infiltrate into partially saturated, near-surface soils, where it may feed evapotranspiration and/or further percolate into saturated soil-rock formations below (Kirchner & Allen, 2020). Each of these regions is structurally heterogeneous over a range of length scales and subject to anomalous transport (Beven & Germann, 1982;Jackisch et al., 2017;Wienhofer et al., 2009). Anomalous transport can arise from the wide range of residence times that are caused by the time-dependent mixtures of rainfall, surface runoff and pre-event water fractions, and the varying concentrations of chemical species that they contain.
To quantify chemical transport behavior in such a complex system, we consider the effective travel time distribution response to a pulse of rainfall over the entire area of a catchment. Thus chemical and isotopic tracers in precipitation enter the catchment almost simultaneously at each point on its surface. For simplicity, but maintaining a realistic picture, we assume that a stream runs through the catchment. This stream acts as a collector-a line sink that drains the catchment-for the tracer. Typically, continuous tracer concentration measurements (tracer arrivals) are made at a control (gauging) point downstream, which yields a time series of tracer concentrations. As the travel time in the stream is very short compared to the travel time in the catchment, we consider in the following the streambank at x = 0 as the point of reference. The travel time distribution from a (pulse) source at distance x to the streambank is denoted by Φ(x, t). The concentration c(t) of a chemical tracer at the streambank is then with r(x, t) the rainfall input of a chemical tracer at position x in a hillslope of length L, and w(x) is a weighting function proportional to the fraction of the drainage area that lies at distance x from the stream (Kirchner et al., 2001). The drainage area is assumed to be regular with w(x) = 1. We note that explicit quantification of evapotranspiration, and wet and dry periods, which affect water and chemical outflows in both space and time, is far from straightforward. Here, these effects are included in an implicit sense: the 36-year measurement period reflects all inputs and outputs to the system, including those of evapotranspiration.
Chemical arrivals are sampled downstream, and can be considered an "instantaneous" integration of all chemical tracer arrivals, from all catchment pathways, along the entire length of the stream. The travel times within the stream can be neglected, because stream velocities are generally much faster than downslope subsurface transport velocities within the catchment. By integrating over all precipitation tracer inputs that reach the stream, we can determine the total flux into the stream. This flux is the first-passage time distribution (overall travel time distribution) at the downstream point of measurement. Given a uniform rainfall distribution over Ω, it follows that r(x, t) ≈ r(t), and we can hence define for the overall travel time distribution h(t) in the entire catchment: Note that the time integral over the travel time distribution Φ(x, t) integrates to one. Kirchner et al. (2000Kirchner et al. ( , 2001 originally analyzed the then-available ∼10 years of time series measurements from the Plynlimon catchments, comparing, in particular, time series of measurements of chloride concentration in rainfall to the chloride tracer concentrations in the (outlet) collector Hafren stream. In this context, a spectral analysis of the concentration signal in the stream was first performed; the analysis considers the Fourier transform where f is the frequency with units of inverse time and i denotes the imaginary unit. Thus, large frequencies correspond to short times (high frequency fluctuations), while small frequencies correspond to long times (low frequency fluctuations).
The power spectrum S c (f) is defined by where the bar denotes the complex conjugate. The power spectra S h (f) of h(t) and S r (f) of r(t) are defined analogously. S c (f) can be written in terms of S h (f) and S r (f) as Thus, S h (f) is obtained by division of the observed spectra of the concentration and rainfall signals. In the following, we refer to S h (f) as the transfer function.
Spectral analysis showed the rainfall chloride concentrations to exhibit an approximate white noise spectrum, while streamflow concentrations displayed a fractal 1/f 2α scaling spectrum (Kirchner et al., 2000), which is consistent with a Gamma distribution of overall travel times h(t). Kirchner et al. (2001) showed that the exponent α = 1/2 can be obtained from an inverse Gaussian travel time distribution Φ(x, t), which assumes that transport along the catchment can be characterized by a mean flow velocity and a macrodispersion coefficient. Recently, matrix diffusion typical of fractured aquifers was invoked as a possible explanation of the observed power-law scalings (Rajaram, 2021).

Continuous Time Random Walk Model
In this letter, we adopt the CTRW framework to integrate the impact of broad distributions of mass transfer times, typical of transport in heterogeneous media (Berkowitz et al., 2006), on the distribution of overall travel times in catchments. The approach is used to analyze and interpret the extensive 36-year record (1983-2019) of chloride tracer concentration measurements from three catchments in Plynlimon, Wales (Kirchner & Neal, 2013;Kirchner et al., 2000;Neal et al., 2011Neal et al., , 2012. The Upper Hafren, Lower Hafren, and Tanllwyth catchments have areas of 1.22, 3.47, and 0.916 km 2 , respectively, and can be abstracted to be of roughly rectangular shape (longer in the upstream to downstream direction); see Neal et al. (2011) for a detailed description and map of the catchments. The measurements were processed using a spectral analysis of the concentration signal in precipitation and at the stream outlet, using the methods described in detail in Kirchner and Neal (2013); see Supporting Information S1 for additional information.
In the development here, the catchment is assumed to be composed of hydrological "compartments" of characteristic length ℓ. Each compartment j is characterized by a distribution ψ j (t) of transition times (as distinct from "travel times," which are defined differently in Section 2.1), which encodes the different retention and transport processes on the compartment scale. The compartments are assumed to be statistically independent and of uniform length ℓ. In this framework, the transport of a solute parcel through the catchment is given by The transition times, τ j , are distributed according to ψ j (t). The travel times t n from the nth compartment toward the stream at the distance nℓ are then given by Thus, the distribution Φ n (t) of t n is given by the n − fold convolution of the ψ j (t), which reads in Fourier space as Note that this formulation accounts for the variability of travel time distributions along the hillslope. For simplicity, here, we assume that ψ j (t) = ψ(t) is compartment-independent. This is a reasonable assumption when the bulk of the water flow is in the subsurface, and when the spatial subsurface properties can be considered statistically stationary from watershed boundaries to the streambank.
The overall travel time distribution h(t) results from the recharge and flow from N compartments toward the outflow. Using Equation 9 for ψ j (t) = ψ(t) in the Fourier transform of the discrete version of Equation 2, we can write (see Supporting Information S1 for details) This geometric sum can be evaluated explicitly, which gives The transition times τ n are assumed to be broadly distributed, characterized by the power law ψ(t) ∼ t −1−β with 0 < β < 1. We use here the one-sided stable distribution (Nolan, 2018;Uchaikin & Zolotarev, 1999) with 0 < β < 1, τ ℓ = a 0 ℓ 1/β the time scale that indicates the onset of the power law, and a 0 a scale parameter. The inverse Fourier transform of ψ*(f) decays as t −1−β for t ≫ τ ℓ and decreases exponentially fast toward 0 for t ≪ τ ℓ . Furthermore, transition times are shifted toward smaller values with decreasing τ ℓ , that is, transitions are on average faster. These features are illustrated in Figure 1. Note that the stable distribution Equation 12 is an operational assumption. In general, the transition time distribution is characterized by a cut-off that corresponds to the largest mass transfer time over a characteristic heterogeneity scale (Dentz et al., 2004). In the present context, this time scale would correspond, for example, to diffusion times across low permeability regions. Particularly with increasing depth, subsurface preferential pathways for fluid flow and tracer transport in heterogeneous and 10.1029/2023GL104207 5 of 8 fractured rock formations can become increasingly ramified, with small fractures tending to be compressed and less conductive and less weathered host rock having lower permeability. We thus assume here that the cut-off time is of the order of or larger than the largest observation time of about 30 years, which corresponds to characteristic molecular diffusion lengths in rock of ∼0.3-0.7 m.
The typical compartment size ℓ may be estimated from the characteristic length scale of the distribution of physical aquifer heterogeneity. We assume here that this scale is much smaller than the hillslope length ℓ ≪ L and consider in Equation 11 the scale limit N = L/ℓ → ∞ such that Nℓ = L. In this scale limit, we obtain for h*(f) where we defined τ L = τ ℓ (L/ℓ) 1/β . The corresponding h(t) scales as see also Scher et al. (2002) and the Supporting Information S1. The transfer function ℎ = ℎ * ( )ℎ * ( ) corresponding to Equation 13 is * The detailed derivations of Equations 13-15 can be found in the Supporting Information S1.

Results and Discussion
Fits to the data for each of the three catchments, using the explicit Equation 15, are shown in Figures 2-4. For comparison, corresponding fits of the spectra of the advection-dispersion equation (ADE)-based solution, with underlying Fickian transport are also shown. Details on the fitting method and performance indicators are given in the Supporting Information S1.
The data for the power spectra S h (f) of the three catchments show similar behaviors. At intermediate and high frequencies, 1/yr < f < 1000/yr, the data have a power-law tangent. With decreasing frequencies, they asymptote toward 1. The two models show relatively good fits for f > 10/yr. For low frequencies, the ADE-based spectra converge faster toward 1 than the CTRW-based spectra, for which we assume that the cut-off scale is larger than the largest observation time.
The ADE-based spectra are fit by the time scale τ c = D/v 2 between 1 and 7 years, which is required to extend the power-law tangent f −1 to sufficiently low frequencies and fit the data. As outlined by Kirchner et al. (2001), this implies that the dispersion lengths L D = D/v are similar to the lengths of the hillslopes. However, as noted by Gelhar and Axness (1983), "on the order of 10-100 dispersion lengths of downstream displacement will be required before" the macrodispersion concept is strictly applicable.
The relatively low, best-fit values of β in the CTRW model indicate chemical transport in a highly dispersive, heterogeneous medium, as noted from the spectra by Kirchner et al. (2001) and Kirchner and Neal (2013), with a  particularly broad range of tracer transition times through the entire catchment (Berkowitz et al., 2006). This is supported by hydrogeological information on the three catchments, which consist largely of acidic soils that overlie fractured mudstone and shale bedrock (Kirchner & Neal, 2013). However, relatively few data points are available for small frequencies, that is, large observation times. Furthermore, we find that Equations 13 and 11 with Equation 12 are indistinguishable for N ≥ 300. This implies that for a catchment length of ∼1 km, the characteristic transition length would be approximately 3 m. Note that because we focus on the scale limit Equation 13, the characteristic length scale of the compartments is much smaller than the length L of the hillslope; as such, the results are independent of the number of compartments.
To further illustrate the above, Figure 5 shows the catchment residence time distributions corresponding to the fitted spectrum of S h based on the CTRW, and ADE models (see Supporting Information S1), for Upper Hafren as shown in Figure 2. The residence time distribution corresponding to the CTRW model Equation 15 scales according to the behaviors given in Equation 14. Note that the power-law behavior at large times relies on the operational assumption that the cut-off scale of ψ(t) is of the order of or larger than the largest observation time. The ADE-based model exhibits t −1/2 -scaling at short times and exponential decay toward zero at large times.
The small Péclet numbers associated with the ADE formulation imply that a stable macrodispersion coefficient (as discussed by Gelhar and Axness (1983)) has not yet emerged after this transport distance, which implies the existence of preferential flow at all scales. Similarly, by its formulation, the CTRW analysis indicates that preferential flow is occurring at all length and time scales.

Conclusion
We investigated the occurrence of anomalous transport in a hydrological catchment, at large spatial scales and over long times. Using spectral analysis, we show that the behavior of long-term (36 years) time series measurements of a natural passive tracer (chloride), in rainfall and runoff in a hydrological catchment, can be described by a CTRW approach that accounts for anomalous chemical transport behavior. The ADE spectra shown here are likewise indicative of anomalous transport, because they imply that subsurface permeability is heterogeneous at  all scales, up to the flowpath length itself (Kirchner et al., 2001). Natural flow systems are not characterized by heterogeneity only at small scales, and thus do not fit the Fickian paradigm. The CTRW model is founded on the presence of a broad distribution of travel times, which may be characterized by a power-law regime that determines both the short and long time behaviors. CTRW provides a modeling framework to account for the impact of broad distributions of travel times on large scale solute transport. The ability to quantify such behavior opens an important path to analyzing chemical transport in hydrological catchments.

Data Availability Statement
No new datasets were generated in this study. The data on which this article is based are available in Kirchner et al. (2000), Neal et al. (2011), Neal et al. (2012, 2013a, 2013b and Norris et al. (2017a), Norris et al. (2019a).