Unifying Advective and Diffusive Descriptions of Bedform Pumping in the Benthic Biolayer of Streams

Many water quality and ecosystem functions performed by streams occur in the benthic biolayer, the biologically active upper (~5 cm) layer of the streambed. Solute transport through the benthic biolayer is facilitated by bedform pumping, a physical process in which dynamic and static pressure variations over the surface of stationary bedforms (e.g., ripples and dunes) drive flow across the sediment‐water interface. In this paper we derive two predictive modeling frameworks, one advective and the other diffusive, for solute transport through the benthic biolayer by bedform pumping. Both frameworks closely reproduce patterns and rates of bedform pumping previously measured in the laboratory, provided that the diffusion model's dispersion coefficient declines exponentially with depth. They are also functionally equivalent, such that parameter sets inferred from the 2D advective model can be applied to the 1D diffusive model, and vice versa. The functional equivalence and complementary strengths of these two models expand the range of questions that can be answered, for example, by adopting the 2D advective model to study the effects of geomorphic processes (such as bedform adjustments to land use change) on flow‐dependent processes and the 1D diffusive model to study problems where multiple transport mechanisms combine (such as bedform pumping and turbulent diffusion). By unifying 2D advective and 1D diffusive descriptions of bedform pumping, our analytical results provide a straightforward and computationally efficient approach for predicting, and better understanding, solute transport in the benthic biolayer of streams and coastal sediments.


Introduction
The movement of water into and out of the hyporheic zone, or "hyporheic exchange," occurs over a wide range of spatial (and temporal) scales, from >10 km (>1 year) to <1 m (<1 hr) Gomez-Velez & Harvey, 2014;Wörman et al., 2007). This >10 3 range of temporal and spatial scales raises trade-offs-relative to residence times, reaction times, and exchange rates-that can influence the hyporheic zone's ability to process nutrients and other pollutants (Harvey et al., 2013). For example, Gomez-Velez et al. (2015) evaluated the residence time/exchange rate trade-off for aerobic respiration and denitrification in the Mississippi River Network, calculating for each reach a so-called reaction significance factor (RSF) (Harvey et al., 2013). In the RSF framework, more nutrients are removed when hyporheic zone residence times are comparable to reaction times and the uptake length is short compared to the reach length (i.e., the RSF is large). These authors found that the smallest scales of hyporheic exchange are the most important for nutrient processing in streams, with RSFs consistently larger for vertical exchange over submerged ripples and dunes (length scales of the order of 10 0 m) compared to lateral exchange over larger geomorphic features such as river bars and meandering banks (length scales of the order of 10 2 to 10 3 m). This conclusion, which is based on physical arguments, is reinforced by findings that microbial biomass and nitrification and denitrification potential tend to be concentrated in the upper 5 cm of the streambed, a region of the hyporheic zone known as the "benthic biolayer" (Caruso et al., 2017;Knapp et al., 2017;Tomasek et al., 2018). Collectively, these results underscore the importance of elucidating physical mechanisms responsible for hyporheic exchange at the scale where nutrient transformations primarily occur, that is, in the benthic biolayer.
At the scale of the benthic biolayer, one important driver of hyporheic exchange is bedform pumping, which occurs when dynamic and static pressure variations over the surface of bedforms (e.g., ripples and dunes) drive flow across the sediment-water interface (SWI) in spatially isolated upwelling and downwelling zones (Azizian et al., 2015(Azizian et al., , 2017Cardenas et al., 2008;Elliott & Brooks, 1997a, 1997bFleckenstein et al., 2010;Grant et al., 2012Grant et al., , 2014Thibodeaux & Boyle, 1987) (Figure 1a). Since its discovery in 1987 (Thibodeaux & Boyle, 1987), analytical models have been proposed to describe bedform pumping and its influence on stream water quality (reviewed in Boano et al., 2014). Generally, these models can be grouped depending on whether they conceptualize bedform pumping as an advective or diffusive process. Advective models are notable for their relatively faithful representation of the laminar flow fields generated by bedform pumping (Elliott & Brooks, 1997a, 1997b). An advantage of diffusive models is their ability to incorporate multiple Shown are bedform morphology, streamlines through the sediment, pressure variation over the surface of the SWI, and upwelling and downwelling zones (upward and downward facing arrows). (b) A simplified analytical model of this process, the bedform pumping model (BPM), assumes a sinusoidal pressure head variation over a Toth domain representation of a rippled bed (see main text). Streamlines through the sediment consist of repeating identical (or mirror image) unit cells (a single unit cell is indicated by the two red vertical dotted lines). (c) Blow-up of the unit cell extending from x ¼ −π=2 to x ¼ π=2. Each streamline in the unit cell is uniquely identified by where it enters the sediment at x ¼ x 0 , 0 < x 0 < π=2. Shown are the solute concentrations at the entrance, C w (t), and exit, C w t − τ x 0 ð Þ ð Þ , locations of a single streamline, where τ x 0 ð Þ is the residence time of the streamline that enters the sediment at x ¼ x 0 . mechanisms for mass transport across the SWI (i.e., not just bedform pumping) including molecular diffusion, turbulent diffusion, and turbulent dispersion (Chandler et al., 2016;Bottacin-Busolin, 2017;Grant, Azizian, et al., 2018;Grant, Gomez-Velez, & Ghisalberti, 2018;Grant et al., 2020;Roche et al., 2019;Voermans et al., 2017Voermans et al., , 2018. As commonly implemented, both types of analytical models rely on multiple assumptions that limit their practical utility: (1) Solute concentration in the overlying water column is assumed constant in time; (2) two-way coupling across the SWI-whereby mass transfer out of the streambed alters mass concentration in the overlying water column which, in turn, alters mass transfer into the streambed, and so on-is not accounted for; (3) diffusive mixing in the streambed is constant in depth, while the interstitial flow field generated by bedform pumping decays exponentially (Elliott & Brooks, 1997a); and (4) some published diffusive models fail to account for the finite porosity of the streambed, a violation of mass balance that can bias estimates of the diffusivity downwards by a factor of 10 (Grant et al., 2012).
In this paper we derive two parallel analytical frameworks, one advective and the other diffusive, that collectively address the model limitations noted above. The paper is organized as follows. In section 2 we review a canonical 2D analytical model for advective bedform pumping originally developed by Brooks (1997a, 1997b) (section 2.1), show that its residence time distribution (RTD) closely follows the extreme value Fréchet distribution (section 2.2), and derive from this result a set of fully coupled solutions for the evolution of solute concentration in the water and sediment columns of a closed system (section 2.3).
In section 3 we derive a parallel 1D diffusive analytical framework for bedform pumping (section 3.1), show how the choice of a diffusivity profile (constant or exponentially declining) leads to different Green's function (Leij et al., 2000) representations of mass transport in the streambed (section 3.2), and then derive from these Green's functions a set of fully coupled solutions for the evolution of solute concentration in the water and sediment columns of a closed system (section 3.3). We test these models against previously published measurements of unsteady solute transport across artificial and natural bedforms in a recirculating flume (section 4). Discussion of these results is presented in section 5 and conclusions in section 6.

Canonical Solution by Elliott and Brooks
A canonical 2D advective model of bedform pumping (originally solved by Vaux, 1968, and expanded on by Elliott & Brooks, 1997a, 1997b, hereafter referred to as the bedform pumping model (BPM), adopts the Toth approach (Toth, 1962) of approximating hyporheic exchange across a rippled streambed with a sinusoidal pressure head variation over a rectangular domain (cf. Frie et al., 2019) (Figure 1b). The wavelength λ (L) of the pressure wave corresponds to the wavelength of the bedform, and the trough and peak of the pressure wave correspond to where the velocity boundary layer detaches (at the bedform crest) and reattaches (on the lee side of the bedform), respectively (Cardenas & Wilson, 2007a, 2007bSawyer & Cardenas, 2009) (all variables in this paper are defined in Table 1). If the hydraulic conductivity K h (L T −1 ) and porosity θ (-) of the streambed are constant, Darcy's law and the continuity equation can be jointly solved to yield the BPM's well-known formulae for the two-dimensional pressure head distribution and velocity field in the interstitial pores of the hyporheic zone (Equations R1-R5, Figure 1b) (Elliott & Brooks, 1997a, 1997b. As documented in Supporting Information (Text S1), if the sediment bed is initially solute free (at t adv ¼ 0) and the solute in question is conservative (i.e., inert and does not adsorb to sediments), the average advective flux, J adv t adv ð Þ(M L −2 T −1 ), of mass into the streambed can be represented as a convolution over all past water column concentrations, C w t adv ð Þ (M L −3 ): Advective model's normalized horizontal and vertical coordinates (-) X K h , θ Sediment hydraulic conductivity (L T −1 ) and porosity (-) X X u m , h m Maximum Darcy flux (L T −1 ) and pressure head amplitude (L) X h SWI x ð Þ Pressure head variation at the sediment-water interface (L) Normalized pressure head amplitude and normalized Darcy fluxes along horizontal and vertical coordinates (-) X u k k Modulus of the normalized Darcy flux (-) X τ Residence time of a water parcel transiting through the hyporheic zone (T) X C w , C s , t Solute concentration (M L −3 ) in the water column and interstitial fluids of the sediment bed, and calendar and mass (M) in the water column, volume of water (L 3 ) above the sediment bed and in recirculating flume pipes, bed area (L 2 ), bed depth (L), water depth (L), and stream velocity (L T −1 ) Advective model's normalized calendar time and residence time (-), and normalized starting position (on the sediment-water interface) of a water parcel transiting through the hyporheic zone (-) X F RTD τ ð Þ Fraction of water circulating through the hyporheic zone with normalized residence time, τ, or younger (-) Fraction of water circulating through the hyporheic zone with residence within dτ of τ (-) X Fraction of water circulating through the hyporheic zone with residence times within dlog 10 τ of log 10 τ (-) X Normalized solute concentrations in the water column, in the sediment column, and in both the water and sediment columns under equilibrium conditions (-) X X T = h w π/u m , T ¼ T=t T Dimensional (T) and normalized (-) timescale for water in the overlying water column to undergo hyporheic exchange Laplace transform and inverse Laplace transform X X s adv ¼ s adv t T Normalized Laplace transform variable (advective model) (-) X β , μ Fréchet distribution parameters (-) X γ , g Exponent (-) in Fehlman's formula for the pressure head amplitude, and gravitational acceleration (L T −2 ) X y diff ¼ ay, a Normalized depth into the streambed (-) and the inverse depth scale over which diffusivity decays (L −1 ) X D eff (y) Depth dependence of the effective diffusivity ( Normalized critical time at which mass transfer across the SWI deviates from a square root dependence on time (-) Normalized form of the effective water depth (-) X Green's function normalized by the timescale for diffusive mixing in the streambed (-) Green's functions for the exponentially declining and constant diffusivity models (-) X m E t diff ð Þ, m C t diff ð Þ Normalized mass transferred out of the sediment bed after an impulsive input, assuming an exponentially or constant diffusivity profile (-) Laplace domain solution of the Green's function for diffusive mixing in the streambed (-) X Permeability Reynolds number (-) depends on the shear velocity u * (L T −1 ), sediment permeability (L ), and kinematic viscosity (L 2 T −1 ) X X f D ¼ 8u 2 * =V 2 Darcy-Weisbach friction factor (-) X X

10.1029/2020WR027967
Water Resources Research The function f RTD τ ð Þ (-) is the probability density function (PDF) form of the BPM's RTD, defined such that the quantity f RTD τ ð Þdτ is the fraction of water circulating through the hyporheic zone with dimensionless residence times in the range τ to τ þ dτ . The variable u m appearing on the right-hand side of Equation 1a is the maximum Darcy flux of water across the SWI, and calendar time and residence time (t adv ¼ t=t T and τ ¼ τ=t T , respectively) have been scaled by an advective timescale for the transport of solute through a bedform: t T = λθ/πu m (BPM variables defined in Figure 1).

The Fréchet Distribution and the BPM's RTD
For any choice of the dimensionless residence time τ, numerical evaluation of the BPM's RTD requires two steps. First, the dimensionless starting position, x 0 τ ð Þ(-), of the streamline in the unit cell with dimensionless residence time τ (see Figure 1c) is obtained by numerically solving the implicit expression for x 0 τ ð Þ (Equation 1c). This estimate of x 0 τ ð Þ is then substituted into the RTD formula (Equation 1b) to obtain the fraction of flow leaving the hyporheic zone with that dimensionless residence time. Because hyporheic zone residence times vary over multiple orders of magnitude, it is convenient to divide the unit area under the RTD into evenly spaced logarithmic increments of dimensionless residence time (Azizian et al., 2017): The cumulative distribution function (CDF) form of the RTD appearing in Equation 2b, F RTD τ ð Þ (-), is defined as the fraction of water circulating through the hyporheic zone with dimensionless residence time of τ or younger; the PDF and CDF forms of the RTD are related in the usual way: f RTD τ ð Þ ¼ dF RTD τ ð Þ=dτ. As demonstrated in the Supporting Information (Text S2), our definition of F RTD τ ð Þ is mathematically consistent with the one derived by Elliott and Brooks (hereafter, EB) in their original publication of the BPM (Elliott & Brooks, 1997a).
The BPM's RTD spans a thousand-fold change in dimensionless residence times, from τ < 0:1 to τ > 100 (black curves in Figures 2a and 2b). The BPM's RTD was compared to five common analytical distributions (Fréchet, Pareto, Log Normal, Gamma, and Exponential). These analytical probability distributions were optimized by randomly sampling the BPM's RTD 10,000 times and employing maximum likelihood estimation to infer distribution parameter values from these realizations (Ang & Tang, 2007) (see Table 2 for mathematical definition of each distribution and optimized parameter values). The BPM's RTD is well described by both the Fréchet and Pareto distributions, reasonably well described by the Log Normal distribution, and poorly described by the Gamma and Exponential distributions (colored curves in Figure 2). The remarkable similarity between the BPM's RTD and the Fréchet distribution-a heavy-tailed extreme value distribution (Kotz & Nadaraj, 2000)-has not, to our knowledge, been noted in the literature. More commonly, the power law or Pareto distribution is adopted to represent hyporheic exchange (Bottacin-Busolin & Marion, 2010). However, the three-parameter version of the Pareto distribution was required to obtain a reasonable match to the BPM's RTD, and, even then, the Pareto distribution ranked second behind the (two-parameter) Fréchet distribution (see Kolmogorov-Smirnov ranking in Table 2). The Log Normal distribution, which is sometimes used to model residence times in the hyporheic zone (e.g., Azizian et al., 2017;Wörman et al., 2002), underpredicts the RTD's heavy tail but is otherwise comparable (green and black curves, Figure 2b). The Gamma distribution has been used to represent the RTD of water parcels transiting through hillslopes (Kirchner et al., 2000;Leray et al., 2016) while the Exponential RTD underpins the transient storage model, a popular hyporheic exchange modeling framework (Knapp & Kelleher, 2020). Based on the results in Figure 2, these last two distributions should not be used to represent the BPM's RTD. For the analysis that follows, we adopted the optimized Fréchet distribution in place of the BPM's RTD for three reasons: (1) The Fréchet distribution is parsimonious and closely matches the BPM's actual RTD (Table 2 and Figure 2); (2) this approach side steps the numerical challenges associated with the two-step process required to solve the BPM's RTD (see Equation 1b and discussion thereof); and (3) the Laplace transform of the Fréchet distribution can be computed analytically, which simplifies the mass balance analysis described next.

Unsteady Solute Concentration in the Water Column of a Closed System
An example of bedform pumping in a closed system is the recirculating flume setup illustrated in Figure 3a.
A mass M of a conservative solute is added to the water column of a solute-free recirculating flume at time t = 0. After a short mixing period, the concentration in the water column is approximately C 0 = M/V w where V w is the volume of water above the sediment bed and in the recirculating pipes. At this point in time, the second term on the right-hand side of Equation 1a is negligible (because no solute has yet passed through the hyporheic zone and returned to the stream), and therefore, the BPM predicts that the initial flux of solute into the bed should be J 0 = C 0 u m /π. With increasing elapsed time (t > 0), the solute concentration in the overlying water column declines (Figure 3b), the integral term in Equation 1a becomes progressively larger in magnitude (as solute in the streambed returns to the stream), and the net flux across the SWI asymptotically approaches zero ( Figure 3c). In practice, if the experiment runs long enough, the water column solute concentration will approach an equilibrium value, C eq (M L −3 ), reflecting a well-mixed final state in which the solute concentration is the same everywhere in the overlying water column and the interstitial fluids of the sediment bed: New variables in Equation 3 include the effective water depth h w (L) (equal to the volume of water V w (M L −3 ) in the overlying water column and recirculating pipes divided by the surface area A b (L 2 ) of the bed, h w = V w /A b ), sediment bed depth d b (L), and porosity θ (-). As the BPM assumes that the streambed is infinitely deep, this analytical model never achieves equilibrium.
One challenge associated with deriving a solution for the recirculating flume problem illustrated in Figure 3a is the two-way coupling of solute concentrations above and below the SWI. This two-way coupling is evident when mass balance is performed over the recirculating flume's water column: Here, the water column concentration has been scaled by its initial concentration at t = 0 (C w ¼ C w =C 0 ), and the new variable T = h w π/u m represents a characteristic timescale for all water in the overlying water column and recirculating pipes to undergo hyporheic exchange. Two-way coupling manifests mathematically as a dependence of the time rate of change of the water column solute concentration (left-hand side of Equation 4) on the entire past history of water column solute concentration filtered through the hyporheic zone's RTD (convolution integral on the right-hand side of Equation 4). In addition to providing an elegant interpretation of two-way coupling, the convolution representation of hyporheic exchange flux permits an analytical solution to the overall mass balance problem. This is because the Laplace transform of a convolution of two variables is equal to the product of their respective Laplace transforms (Graff, 2004). Thus, after applying the Laplace transform to Equation 4, solving for the solute concentration in the water column becomes a simple algebraic exercise: Here, the variable T ¼ T=t T is a dimensionless timescale for hyporheic zone processing of water above the streambed, s adv ¼ st T is a dimensionless form of the Laplace transform variable, s (T −1 ), and the symbol Figure 3. (a) A conceptual diagram of a recirculating flume experiment, in which streamflow over the top of stationary bedforms induces bedform pumping and hyporheic exchange. The water column has an "effective depth" h w equal to the total volume V w of water above the SWI and in the recirculation pipes divided by the area of the SWI A b . (b) In a typical step-change experiment, the concentration of a conservative solute in the overlying water column C w (t) is increased suddenly to C 0 at time t = 0 and is everywhere equal to zero for t < 0. Mixing across the SWI causes C w (t) to decline toward an equilibrium value.
(c) The advective mass flux across the SWI J adv (t) jumps to J 0 at time t = 0 and then declines toward zero over time.
L −1 [·] denotes the inverse Laplace transform which, in practice, is solved numerically (see section 4.1 and Mathematica codes in Supporting Information). The Laplace transform of the Fréchet distribution can be computed analytically by applying the Right Shift rule (Graff, 2004) where K 1 (·) is the modified Bessel function of the second kind and u is a dummy integration variable: Because the Fréchet distribution parameters (β, μ) were previously optimized (see Table 2), we can infer from Equation 5a that C w t ð Þ depends on a single dimensionless parameter, T ¼ T=t T . The two timescales appearing in this dimensionless parameter depend on physical characteristics of the recirculating flume as follows: Therefore, implementation of this analytical solution requires knowledge of the bedform wavelength λ, streambed porosity θ, streambed hydraulic conductivity K h , the effective depth of the water column h w , and the half amplitude of the pressure head variation h m . With the exception of h m , these parameters can be measured or estimated; for example, bedform wavelength and height can be estimated from field measurements of flow depth, flow velocity, and median grain size using bedform stability diagrams and hydraulic criteria (cf. Zheng et al., 2018) while hydraulic conductivity can be estimated using standard relationships (such as the Kozeny-Carman equation, McCabe et al., 2010) or machine learning algorithms (Stewardson et al., 2016). To estimate h m , the widely cited empirical formula proposed by EB (based on pressure measurements over a triangular bedform reported by Fehlman, 1985) can be employed: H=d 0:34 γ The variables V (L T −1 ) and d (L) represent, respectively, the average velocity and depth of the overlying stream, H (L) is bedform height, and the empirical exponent is taken as either The value of the multiplicative constant (0.28) on the right-hand side of Equation 7 can be adjusted depending on the height-to-wavelength ratio of the bedform (Fox et al., 2014;Shen et al., 1990).

Unsteady Interstitial Solute Concentration in the Sediment Column of a Closed System
A corresponding analytical solution can be derived for the spatiotemporal (2D) evolution of solute concentration in the interstitial fluids of the sediment bed. The solution is premised on the idea that the interstitial concentration of a conservative solute at any dimensionless time t adv is equal to the concentration that was present in the water column some location-dependent dimensionless residence time ago: C s x; y adv ; t adv ð Þ¼ Þ . New variables appearing here include the dimensionless interstitial solute concentra- is the solute mass per unit volume of interstitial fluid (i.e., as opposed to per bulk sediment volume) and τ x; y adv ð Þ¼τ x; y adv ð Þ=t T is a dimensionless form of the location-dependent residence time function, τ x; y adv ð Þ(T), defined as the time it takes interstitial water parcels to travel from the SWI to any x; y adv ð Þposition in the sediment bed (Azizian et al., 2015) (derivation in Supporting Information, Text S3): The unsteady solution for the interstitial concentration of a conservative solute in the streambed directly follows from this last result, where the time-dependent solute concentration in the water column, C w t adv ð Þ, is given by Equation 5a: It should be noted that Equation 8a is valid only within the bounds of the unit cell illustrated in Figure 1c (i.e., −π=2 < x < π=2). Outside of the unit cell, the equation must be translated, with or without reflection, using the following substitution rule for the dimensionless horizontal coordinate: x→ −1 ð Þ n x − nπ ð Þ , where the integer n is given by n ¼ R x=π ½ and the function R[·] rounds to the nearest positive or negative integer value. Finally, a solution for the location of the concentration front in the sediment bed at any dimensionless time t adv can be obtained by substituting t adv for τ on the left-hand side of Equation 8a, and numerically solving the resulting implicit expression for y adv given x, or vice versa. The implicit solution for the concentration front also applies to locations outside of the unit cell (i.e., for x < −π=2 or x > π=2) after translation with or without reflection, using the substitution rule presented above for x.

Diffusive Model of Hyporheic Exchange by Bedform Pumping
Advective models, like the BPM, are premised on the idea that pore-scale advection dominates the transport of solutes in the hyporheic zone. Over the years, researchers have also explored the possibility of employing 1D diffusive models to describe hyporheic exchange, generally, and bedform pumping, in particular (O'Connor & Harvey, 2008). EB, for example, argued that a diffusivity for bedform pumping should take the form of a dispersion coefficient, E ≈ 0.04λu m /θ (L 2 T −1 ), where u m /θ is a characteristic pore-scale velocity associated with the BPM (see Equation R5 in Figure 1) and the mixing length scale is the bedform wavelength, λ (from Equation 11 in Elliott & Brooks, 1997b). Applied to mass transfer in recirculating flumes, the constant diffusivity model predicts that the water column solute concentration, and the penetration depth of solute into the streambed, should both scale with the square root of time (Elliott & Brooks, 1997a, 1997b. In their recirculating flume experiments, EB found that mass transfer across the SWI followed the predicted square root dependence until a transition time of around, t c ≈ 8λθ/u m . Afterward, measured mass transfer rates fell below those predicted by the constant diffusivity model. Similarly, Marion and Zaramella (2005) reported that constant diffusivities inferred from recirculating flume studies decline as the timescale over which hyporheic exchange is measured increases.
From a mechanistic perspective, all these problems with the 1D constant diffusivity model can be rationalized by noting that, as time increases, mass transfer across the SWI slows dramatically as the relative contribution of deeper streamlines to bedform pumping increases. We hypothesize that this effect can be represented by requiring the dispersion coefficient to decline exponentially with depth, consistent with the exponentially declining velocity field that underpins hyporheic exchange by bedform pumping (see Equations R1-R4 in Figure 1) (cf. Bottacin-Busolin, 2019).

Governing Equations for Diffusive Bedform Pumping in a Closed System
Bedform pumping generates concentration fields in the interstitial fluids of the sediment bed that are at least two-dimensional (e.g., for artificially shaped triangular bedforms in the laboratory) and more often three-dimensional (e.g., in natural streams). However, if the goal is to predict average rates of mass transfer over, for example, a stream reach or a recirculating flume, knowledge of the two-and three-dimensional flow and subsurface concentration fields are not required. Thus, for many applications, mass transport and mixing by bedform pumping in the benthic biolayer can be conceptualized as an unsteady one-dimensional diffusion problem, for which the horizontally averaged vertical flux, J diff (y,t) (M L −2 T −1 ), of solute through the sediment is described by a flux-gradient process where the mixing coefficient, or effective diffusivity D eff (y) (L 2 T −1 ), varies with depth y in the benthic biolayer: Grant et al. (2020) demonstrated that Equation 9a is a reasonable descriptor of vertical solute transport by turbulent dispersion and turbulent diffusion through the benthic biolayer of a flat streambed, provided that the diffusion coefficient declines exponentially with depth through the sediment bed. In this paper we hypothesize that a similar result applies to bedform pumping, but with the effective diffusivity replaced by an exponentially declining dispersion coefficient, D eff (y) = E(y) = E 0 e −ay . New variables include a surficial dispersion coefficient, E 0 (L 2 T −1 ), and an inverse length scale, a (L −1 ); the former characterizes the dispersion coefficient at the SWI, while the latter characterizes the depth over which the dispersion coefficient decays. These two parameters are emergent properties of the two-and three-dimensional flow and concentration fields that drive bedform pumping; that is, they are determined by spatial correlations between the time-averaged vertical component of the velocity field and the local mean solute concentration (Voermans et al., 2017). The corresponding one-dimensional mass balance equation can be written as follows: Equation 9b equates the accumulation of mass in a differential horizontal slice of the sediment bed (left-hand side) to the vertical diffusive transport (right-hand side) of a conservative (non-reactive and non-adsorbing) solute (Incropera et al., 2007). The coordinate y increases with depth into the streambed, and its origin (at y = 0) is positioned at the horizontal plane of the SWI.
Substituting the proposed exponential form for the dispersion coefficient into Equation 9b and assuming streambed porosity θ does not vary substantially over the vertical dimension of the benthic biolayer (approximately 5 cm) (Knapp et al., 2017), we arrive at the following mass balance equation for interstitial solute transport in the sediment bed: In Equation 9c, calendar time, t diff ¼ t=t E , has been normalized by a characteristic timescale for dispersive mass transport through the benthic biolayer, t E = 1/a 2 E 0 , depth into the streambed has been normalized by the diffusivity's inverse length scale, y diff ¼ ay, and the interstitial solute concentration has been normalized by the initial concentration in the overlying water column, C s ¼ C s =C 0 (same as for the BPM, see section 2.3). By analogy to the BPM, we also assume that the streambed is initially solute free (Equation 9d), solute concentration drops off to zero deep in the streambed (Equation 9e), and the interstitial solute concentration at the top of the streambed equals the solute concentration in the overlying water column (Equation 9f) where H t ð Þ (-) is the Heaviside function (included here to satisfy the requirements of Duhamel's theorem described later): In writing Equation 9f we have assumed that the interstitial concentration at the SWI is equal to the solute concentration in the overlying water column, which implies that mass transfer into the streambed is not rate limited by convective mass transfer across the concentration boundary layer above the streambed; that is, the dimensionless Biot number, the ratio of timescales for diffusive mixing in the streambed and convective mass transfer across the turbulent boundary layer above the streambed, is much greater than unity (Incropera et al., 2007). For a closed system with a well-mixed water column, like the recirculating flume illustrated in Figure 3a, mass balance over the water column takes the following form: In this equation, the change of solute mass in the overlying water column and recirculation system of the flume (left-hand side) equals the mass transfer rate across the SWI by bedform pumping (represented 10.1029/2020WR027967

Water Resources Research
here as a dispersive process, right-hand side). Streambed porosity θ appears on the right-hand side of the equation to account for the abrupt change in area over which solute mass transport occurs above and below the SWI (Grant et al., 2012). Expressing Equation 10a using the dimensionless variables introduced above for the diffusion equation, we obtain Equation 10b where h w is a normalized form of the effective water depth.

Duhamel's Theorem and Green's Functions
As detailed in Grant et al. (2020), we can link the mass balance equations above (Equation 10b) and below (Equation 9c) the SWI and thereby account for two-way coupling across the SWI, using Duhamel's theorem, an analytical approach for solving the diffusion equation in cases where the forcing function at one boundary is a piece-wise continuous function of time (Guerrero et al., 2013). Duhamel's theorem allows us to express the evolution of interstitial solute concentration in the sediment bed as a convolution of the water column concentration C w t diff ð Þand a so-called auxiliary function C A s y diff ; t diff ð Þwhere v is a dummy integration variable (Guerrero et al., 2013): The auxiliary function is a solution to the same system of equations (Equations 9c-9f), but with the inhomogeneous term replaced by the Heaviside function (cf. Equations 9f and 11e): Substituting the coordinate transformation, ξ ¼ e y diff , into Equation 11b (Yates, 1992) and solving the system of equations in the Laplace domain, we arrive at the following analytical solution for the auxiliary function where s diff ¼ st E is a dimensionless form of the Laplace transform variable and K 1 [·] is the modified Bessel function of the second kind: Duhamel's theorem (Equation 11a) can also be expressed as a convolution of the dimensionless water column concentration, C w t diff ð Þ, and a so-called Green's function, G y diff ; t diff ð Þ(T −1 ), which is normalized here by the dispersive mixing timescale t E , G y diff ; t diff ð Þ¼t E G y diff ; t diff ð Þ : Substituting Equation 12 into Equation 13b yields a Green's function specific to the exponentially declining diffusivity profile (denoted here by the subscript "E"): The similarity between Equation 13a and the convolution integral derived earlier for the BPM (Equation 1a) is striking. In both cases, the water column concentration is convolved with a function (Green's function in the case of the 1D diffusive model and an RTD function in the case of the 2D advective model) that captures the response of the streambed to an impulsive injection of mass into the SWI at dimen- The Green's function above (Equation 13c) is specific for the choice of an exponentially decaying diffusivity profile. For the same set of initial and boundary conditions, a second Green's function can be derived for a constant diffusivity profile (the so-called signaling problem, see Gorenflo & Mainardi, 1998) these two Green's functions, we can interrogate how the choice of a diffusivity profile (i.e., exponentially declining or constant) influences the temporal scaling of mass remaining in the sediment bed following an impulsive input of mass at the SWI at time t diff ¼ 0 . This is because, for t diff > 0 , the upper boundary condition for these two Green's functions is zero, C w t diff > 0 ð Þ¼δ t diff > 0 ð Þ¼0, and therefore, the normalized solute mass in the sediment bed, m t diff ð Þ ¼ ∫ ∞ 0 G y diff ; t diff ð Þ dy diff , will diffuse back into the water column after its injection at time t diff ¼ 0 (the subscripts "E" and "C" denote exponentially declining and constant diffusivity profiles, respectively): As might be expected based on the discussion at the beginning of section 3, the constant diffusivity model predicts that solute mass remaining in the sediment bed declines inversely with the square root of time (Equation 14b). The exponentially declining diffusivity model (Equation 14a) exhibits this same 1= ffiffiffiffiffiffiffi t diff p scaling initially but falls off more rapidly after t ¼ t=t E ≈ 0:3 (Figure 4). This result can be rationalized by noting that, for an exponentially declining diffusivity profile, deeper portions of the streambed are relatively inaccessible to solute injected at the SWI at t diff ¼ 0 and consequently contribute little to the release of stored mass at long times. The similarity between the scaling behavior illustrated in Figure 4 and the scaling behavior described earlier for mass transfer across the SWI in recirculating flumes (see preamble to section 3) is the first indication that our overarching hypothesis -that bedform pumping can be represented by an exponentially decaying diffusivity -may be valid.

Solute Concentration in the Water and Sediment Columns of Closed System
From the results presented above, a set of explicit solutions can be derived for solute concentration in the water column and interstitial fluids of a closed system (Grant et al., 2020):  Table 1 for variable definitions).
These analytical solutions are written in terms of the Laplace transform of the Green's function and its derivative, which, in this context, are tailored to the choice of diffusivity profile. For an exponentially declining diffusivity profile, they are as follows: For a constant diffusivity profile, these two functions can be computed directly from the solution to the signaling problem introduced earlier:

EB's Bedform Pumping Data Set and Model Optimization
To test the parallel 2D advective and 1D diffusive analytical frameworks derived above, we turned to one of the first published recirculating flume experiments specifically designed, along the lines of Figure 3a, to investigate the unsteady transfer of a conservative solute across the SWI by bedform pumping (Elliott & Brooks, 1997b). EB's experiments were conducted with stationary bedforms (either artificial triangular bedforms or natural ripples), a non-adsorbing and stable fluorescent dye (Lissamine), and under various flow velocities (8.6 to 13.2 cm s −1 ), water depths (1.14 to 2.54 cm), and shear velocities (1.3 to 3 cm s −1 ) (Experiment IDs 8,9,12,(14)(15)(16)(17). The sediment bed, which ranged in depth from 12.5 to 22.0 cm depending on the experiment, consisted of medium or fine-grained unconsolidated sand with measured hydraulic conductivity K h = 0.11 and 0.0079 cm s −1 and porosity θ = 0.325 and 0.295, respectively (a summary of experimental conditions is included in the Supporting Information, Table S1). Published over 20 years ago, EB's study remains one of the few where the evolution of dye concentration is followed both above and below the SWI-a feature we will take advantage of here.
Experimental evaluation of our analytical models was carried out in two steps. First, we fit the 2D advective and 1D diffusive models for C w (t) (Equations 5a and 15a, respectively) to EB's measurements of dye concentration in the water column over time.  Table S1). In the second step, parameter values inferred from the water column studies were used to predict the movement of dye through the interstitial fluids of the streambed over time (Equations 8b and 15b). These model predictions were compared to observations of the dye front in the sediment bed, which EB recorded by periodically marking the location of the leading edge of the dye plume on the side of their flume (the wall of the flume was transparent, and dye was visualized with a hand-held UV light).

Evaluation of Model-Predicted Water Column Solute Concentrations
Across all seven experiments, the 2D advective model (Equation 5a) and 1D diffusive model (with an exponentially declining diffusivity profile, Equation 15a) closely conform to EB's time series measurements of dye concentration in the water column (R 2 > 0.9998 for both models, Figures 5a and 5b, also see Tables S2 and S3 in the Supporting Information). For comparison, water column concentrations for Experiment #17 were also simulated with the constant diffusivity model; this involved substituting Equation 15f into Equation 15a and adopting the superficial diffusivity, E 0 , inferred from fitting the exponentially declining diffusivity model to the same data set (Table S3). The constant diffusivity model also conforms to measurements of the water column concentration until around ffiffi t p ≈ ffiffiffiffiffi 30 p min 1/2 ; thereafter, the constant diffusivity model seriously underpredicts observed concentration measurements (blue curve, Figure 5b). Water column concentrations predicted by the constant diffusivity model decline approximately linearly when plotted against ffiffi t p , consistent with EB's observations (see preamble to section 3) and the ffiffi t p scaling of the constant diffusivity's Green's function (see Equation 14b and Figure 4).
We can also evaluate the advective and diffusive models based on how well their inferred parameter values reproduce values expected based on theory or measurements. For example, values of the half-amplitude pressure head inferred from the advective model (ranging from h m = 0.04 to 0.57 mm) are similar (roughly factor of 2 or better) to values estimated from EB's empirical formula (Equation 7) (ranging from h m = 0.09 to 0.31 mm, Table S2). Likewise, values of the effective water depth inferred from the advective model (ranging from h w = 8.8 to 16.7 cm) are similar (roughly factor of 2 or better) to values estimated from reported flume water volume (excluding interstitial fluid) and streambed area (h w = V w /A b = 11.3 to 12.5 cm) (Table S2). Deviations between inferred and predicted (or measured) values of h m and h w do not necessarily imply that the model-inferred values are incorrect. For example, the half-amplitude pressure head values predicted by Equation 7 are only approximately correct (Fox et al., 2014;Shen et al., 1990). Measurement errors associated with flume water volume and bed surface area (which may be difficult to define, given the undulatory nature of the SWI with bedforms) also contribute uncertainty and bias to experimental estimates of h w .
A more rigorous assessment of the inferred parameter values can be framed as follows: Can parameter values inferred from the 2D advective model be translated directly into parameter values for the 1D diffusion model and vice versa? To answer this question, we equated EB's proposed formula for a bedform pumping dispersion coefficient to the diffusivity model's surficial dispersion coefficient: E 0 ≈ 0.04λu m /θ. Substituting the BPM's solution for the maximum Darcy flux (u m , Equation R5 in Figure 1), this formula predicts that the diffusive model's surficial dispersion coefficient is directly proportional to the advective model's half-amplitude pressure head, E 0 ≈ 0.08πK h h m /θ. When values of h m inferred from the 2D advective model are substituted into this formula, the predicted values of E 0 are highly correlated with values of E 0 inferred from the 1D diffusion model (Pearson's correlation coefficient, R = 0.99) (open black circles, Figure 5c). Adjusting the equation's pre-factor to correct the bias evident in the figure, we arrive at the following regression relationship between the 2D advective and 1D dispersive descriptions of bedform pumping (blue filled triangles, Figure 5c): E 0 ≈ 0:133πK h h m =θ; 0:08 ≤ K h mm s −1 À Á ≤ 1:1; 0:042 ≤ h m mm ð Þ≤ 0:11; 0:295 ≤ θ ≤ 0:325 (16a) Likewise, the inverse decay length scale, a, inferred from the 1D diffusion model is highly correlated (R = 0.95) with the inverse of the average bedform wavelength ( Figure 5d): The inequalities indicate the range of EB's parameter values over which each regression was trained. Equations 16a and 16b provide a direct link between our 2D advective and 1D diffusive descriptions of bedform pumping, such that a parameter set for one can be directly translated into a parameter set for the other. The implication is that these two descriptions of bedform pumping are, in fact, functionally equivalent, provided that the limitations with existing analytical models outlined earlier (water column concentration constant in time, diffusivity is constant in depth, two-way coupling across the SWI neglected, and the sediment bed's finite porosity neglected) are properly addressed, as they have been in this study. Because Equations 16a and 16b are calibrated with data from EB's study alone, they are per force restricted to a limited range of flow and streambed conditions (indicated by the inequalities above). A meta-analysis is underway to evaluate the predictive power of these equations beyond the set of experiments analyzed here.

Evaluation of Model-Predicted Interstitial Solute Concentrations
We can also evaluate the advective and diffusive models relative to their ability to predict the unsteady transport of dye plumes through the interstitial fluids of the sediment bed. The progression of one such plume beneath an artificial triangular ripple (EB's Experiment #9) is reproduced in Figures 6a-6d (thick dashed black curves). The dye plume penetrated to a depth of about 8 cm in the first 75 min but required an additional 575 min to progress downward another 4 cm. To compare these observations with the advective model solution, the BPM's coordinate system must first be aligned with EB's triangular ripple.
To this end we used the parameter values estimated from the water column optimization study of Experiment #9 (see Table S2) to predict (with Equation 8b) the interstitial dye concentration in the sediment bed at t = 75 min, coinciding with EB's first dye front measurement. The model's horizontal coordinate was then adjusted to align the left and right edges of the observed and predicted dye fronts. Finally, the model's vertical coordinate was adjusted so that the top of the (flat) model domain is equidistant between the crest and trough of the triangular bedform (final alignment is shown in Figure 6a). After making these adjustments, the advective model's predictions for the downward migration of the dye plume over time closely agree with EB's observations of the dye front at t = 150, 320, and 650 min (Figures 6b-6d).
Two-way coupling is evident in the model-predicted interstitial concentration field. Predicted dye concentrations are elevated along the front of the plume because water parcels at the front moved into the sediment bed near time t = 0 when dye concentration in the water column was near its initial (maximum) value, C w = C 0 . As time progresses, dye concentration in the water column declines and the trailing edge of the dye plume, which consists of younger water parcels, becomes less concentrated. This pattern-high dye concentration along the plume's front and low concentration along the plume's trailing edge-is particularly striking for the simulation at t = 150 min ( Figure 6b). Eventually, the plume's concentration field takes on a more uniform appearance as older water parcels (with higher dye concentrations) return to the stream along slow moving streamlines (Figure 6d).
Thus far, we have found little difference between our 2D advective and 1D diffusive models of bedform pumping. One aspect where these two models differ substantially is their respective concentration depth profiles (Figure 6e). The diffusivity model's depth profiles are convex in shape and characterized by a diffuse concentration front that becomes increasingly smeared out over the vertical extent of the streambed with increasing time. By contrast, the advective concentration depth profiles (generated by horizontally averaging the two-dimensional concentration fields appearing in Figures 6a-6d) are convex and characterized by persistent and very sharp concentration fronts. Integration of these advective and diffusive depth profiles over depth (at any fixed time) indicates that both models transfer roughly the same solute mass into the streambed (within 5%, see Figure S1, Supporting Information). Thus, while the two models appear to be functionally equivalent with respect to the transport of a conservative solute (as in EB's experiments), their very different concentration depth profiles could lead to divergent predictions for the transport of reactive solutes through the benthic biolayer. This begs the question: Which of these two profiles is more representative of natural systems?
To answer this question, we turned to recirculating flume experiments EB conducted with stationary natural ripples. These experiments entailed operating the flume under high flow conditions (to induce sediment transport and ripple formation) and then lowering the flow velocity (to immobilize the bedforms and conduct he dye exchange experiments). Not surprisingly, dye plumes generated by natural ripples are variable with respect to their horizontal extent and the depth to which dye penetrates the streambed (Figure 7a). This variability, which arises from variations in bedform geometry (i.e., height H and wavelength λ) and the three-dimensional nature of natural ripples, can be formally analyzed using spectral methods (e.g., Stonedahl et al., 2010). However, if the goal is to obtain bulk estimates for the downward progression of solute through the benthic biolayer over time, both the 2D advective and 1D diffusive analytical solutions derived in this study perform remarkably well (compare Figures 7a-7c). The sharp dye fronts predicted by the 2D advective model are comparable to patterns of dye penetration beneath "average" bedforms; for example, the two dye plumes located 18 to 30 cm along the horizontal axis (Figure 7a). The smeared-out dye fronts predicted by the 1D diffusive model, on the other hand, may be more representative of the concentration profile one would obtain by horizontally averaging the interstitial concentration field across all bedforms (this hypothesis could not be tested with EB's data set because these authors recorded the time evolution of concentration fronts, not concentration fields). What these results imply for reactive solute transport through the benthic biolayer is an interesting topic for future study.

Discussion
The functional equivalence of the analytical 2D advective and 1D diffusive frameworks derived here implies that their application can be tailored to the problem at hand. The advective model is a relatively faithful  Table S1).
representation of the two-dimensional interstitial flow fields associated with bedform pumping. Consequently, this framework will be useful in cases where knowledge of flow paths through the benthic biolayer, and their associated Darcy fluxes and residence times, is required. The removal of stream-borne particles in the benthic biolayer by deep bed filtration, for example, requires detailed information about the interstitial flow field. This is because, as particles move through the streambed, their filtration rate depends on the local flow velocity (through the contact efficiency η 0 [-], see Tufenkji & Elimelech, 2004) which varies continuously along a streamline (see Equations R1-R5 in Figure 1). Another example is the spatial zonation of interstitial oxygen concentration beneath bedforms, including the formation of so-called anoxic chimneys in upwelling zones (Azizian et al., 2015;Kessler et al., 2012Kessler et al., , 2013. This biogeochemical zonation, which arises from the coupling between in-bed redox reactions and bedform pumping of electron donors and acceptors, can impose significant constraints on important streambed functions, such as coupled nitrification-denitrification (Kessler et al., 2013). Because the advective model's flow field is quantitatively linked to bedform geometry and stream flow (i.e., bedform height and wavelength, as well as stream depth and velocity, see Equation 7), physicochemical (e.g., particle filtration) and biogeochemical (e.g., nutrient transformation) functions of the benthic biolayer can be tied directly to geomorphic processes, such as the adjustment of bedform morphology to changes in land use and flow regime, for example, as a result of urbanization (Harvey et al., 2012).
On the other hand, a strength of the 1D diffusive model is its ability to combine multiple mechanisms for mass transport across the SWI. As noted earlier, mixing across a flat SWI can be characterized by an effective diffusivity, D eff , that incorporates three transport mechanisms (Grant et al., 2012;O'Connor & Harvey, 2008;Richardson & Parr, 1988;Voermans et al., 2018): (1) (Voermans et al., 2017(Voermans et al., , 2018. For turbulent mass transfer across a flat SWI, and accounting for the exponential decay of diffusivity with depth, the surficial effective diffusivity exhibits different permeability Reynolds number scaling behavior in the dispersive (D eff; 0 ∝ Re 2:5 K ,0.01 < Re K < 1) and turbulent diffusive (D eff; 0 ∝ Re 1 K , Re K > 1) regimes (Grant et al., 2020). Our formula linking 2D advective and 1D dispersive descriptions of bedform pumping (Equation 16a) implies that dispersive mixing by bedform pumping also increases with the permeability Reynolds number, E 0 ∝ Re 2 K . This last result can be demonstrated by substituting into Equation 16a definitions for the Darcy-Weisbach friction factor, f D ¼ 8u 2 * =V 2 (-) (Sabersky & Acosta, 1989), and the streambed permeability, K = υK h /g (McCabe et al., 2010), and noting that the friction factor will likely scale with the fraction of the water depth taken up by a bedform, f D ∝ (H/d) γ . Thus, the one-dimensional diffusive framework captures solute mixing through the benthic biolayer by molecular diffusion, turbulent dispersion, turbulent diffusion, and bedform pumping.
Another advantage of the 1D diffusive model is that it can be readily modified to account for groundwater recharge or discharge (through the addition of an advective term to Equation 9b) and the inclusion of biogeochemical reaction networks, such as the Monod kinetic expressions associated with nitrogen cycling in streambeds (respiration, ammonification, nitrification, and denitrification, Azizian et al., 2017). While surface-groundwater exchange can be factored into the 2D advective model as well (cf. Boano et al., 2008), doing so invalidates a key requirement of the BPM analytical solution namely, that the x component of the velocity is everywhere constant along a streamline (see Text S1 in Supporting Information). While outside of the scope of this paper, it is also interesting to note that, at sufficiently high celerity, bedform migration tends to reduce the complexity of the interstitial concentration fields, in effect transforming the two-and three-dimensional concentration fields associated with bedform pumping across stationary bedforms into simple one-dimensional vertical concentration gradients (e.g., of interstitial oxygen concentrations, Wolke et al., 2020) that may also be amenable to analysis with a 1D-diffusive modeling framework.
A benefit of analytical models (compared to numerical simulations) is the relative ease with which they can be implemented, and the physical insights afforded by expressing the quantity of interest (e.g., interstitial solute concentration or mass flux across the SWI) as an explicit function of key system variables. An obvious limitation is that their derivation entails simplifying assumptions that may not be valid in practice. Key assumptions associated with the 2D advective and 1D diffusive modeling frameworks derived here include the following: (1) The interstitial flow field underlying bedform pumping is assumed to be steady state, although solute concentrations in the water column and interstitial fluids of the streambed can vary with time while fully accounting for two-way coupling across the SWI; (2) for uniform flow over a homogeneous sediment bed, the stream gradient drives a stream parallel flow in the lower part of the bed, that is, an underflow. Underflow, which is not accounted for in the 2D advective and 1D diffusive modeling frameworks developed here, can reduce the depth to which advective flow paths extend into the bed (thereby reducing mass transfer at long times, Bottacin-Busolin & Marion, 2010) and possibly alter the effective diffusivity's vertical structure. On the other hand, the general concordance of our two models with Elliott and Brooks' data set ( Figure 5) suggests that this process may exert a second-order effect on measured exchange rates; (3) the advective model is premised on a simplified two-dimensional representation of bedform pumping, while in reality the interstitial flow fields responsible for hyporheic exchange across the benthic biolayer are three-dimensional and temporally varying (e.g., due to bedform migration, Zheng et al., 2018); (4) the solutions are specific to closed systems (such as recirculating flumes) while most hyporheic exchange problems of practical interest are open systems (such as streams and coastal sediments); (5) the transport properties of the benthic biolayer (such as porosity and hydraulic conductivity) are assumed constant while in practice hydraulic conductivity varies spatially at the scale of both single and multiple geologic facies (Tonina et al., 2016) and with time by colmation (Stewardson et al., 2016) and bioclogging (Caruso et al., 2017); and (6) as already noted, our results assume that the solute in question is conservative (i.e., non-reactive and does not absorb to the porous matrix) while most hyporheic exchange problems of practical interest involve reactive solutes or particles. These limitations can be addressed, to varying degrees, within the context of our analytical frameworks, and efforts to do so are currently underway.

Conclusions
In this paper we derived two parallel analytical frameworks, one advective and the other diffusive, that together relax many of the assumptions that limit the practical utility of presently available analytical models for bedform pumping. Both frameworks allow the water column concentration to vary with time while accounting for the two-way coupling of solute concentrations above and below the SWI; the 1D diffusive framework additionally allows the mixing rate, or diffusivity, to vary with depth through the sediment bed. When applied to previously published measurements of bedform pumping in a recirculating flume (Elliott & Brooks, 1997a), we find that both analytical frameworks closely reproduce average patterns and rates of hyporheic exchange, provided that the 1D diffusion model's diffusivity declines exponentially with depth. Practical application of these two frameworks can be tailored to the problem at hand, depending on whether detailed knowledge of the interstitial flow fields and associated Darcy fluxes and residence times is required (2D advective model) or solute transport across the SWI is subject to multiple transport mechanisms, not just bedform pumping (1D diffusive model). Because the advective framework is grounded in a physical description of bedform pumping, it explicitly accounts for how changes in stream flow and sediment transport (e.g., associated with urbanization) influence bedform geometry (wavelength and height), the half amplitude of the pressure head variation, and the hydraulic conductivity of the sediment bed. The 1D diffusivity framework, on the other hand, lumps these geomorphic processes into a surficial dispersion coefficient and an inverse decay length scale that can be directly calculated from the aforementioned advective model parameters (see Equations 16a and 16b). These formulae also predict that the surficial dispersion coefficient for bedform pumping increases with the dimensionless permeability Reynolds number, consistent with diffusivities measured for turbulent exchange across flat streambeds (Grant et al., 2020;Voermans et al., 2018) and streambeds with bedforms (Grant, Gomez-Velez, & Ghisalberti, 2018;Grant et al., 2012;O'Connor & Harvey, 2008). Efforts are currently underway to extend these analytical solutions to open systems (e.g., stream networks), bedform turnover, unsteady flows, and the nonlinear reactions that drive nutrient cycling in the benthic biolayer of streams.

Data Availability Statement
Data are available in Elliott and Brooks (1997b).