A One‐Dimensional Model for Turbulent Mixing in the Benthic Biolayer of Stream and Coastal Sediments

In this paper, we develop and validate a rigorous modeling framework, based on Duhamel's Theorem, for the unsteady one‐dimensional vertical transport of a solute across a flat sediment‐water interface (SWI) and through the benthic biolayer of a turbulent stream. The modeling framework is novel in capturing the two‐way coupling between evolving solute concentrations above and below the SWI and in allowing for a depth‐varying diffusivity. Three diffusivity profiles within the sediment (constant, exponentially decaying, and a hybrid model) are evaluated against an extensive set of previously published laboratory measurements of turbulent mass transfer across the SWI. The exponential diffusivity profile best represents experimental observations and its reference diffusivity scales with the permeability Reynolds number, a dimensionless measure of turbulence at the SWI. The depth over which turbulence‐enhanced diffusivity decays is of the order of centimeters and comparable to the thickness of the benthic biolayer. Thus, turbulent mixing across the SWI may serve as a universal transport mechanism, supplying the nutrient and energy fluxes needed to sustain microbial growth, and nutrient processing, in the benthic biolayer of stream and coastal sediments.


Introduction
Many physical and biological processes in aquatic ecosystems depend on, or are strongly affected by, turbulent fluid motions at the sediment-water interface (SWI) (Franca & Brocchini, 2015;Grant, Azizian, et al., 2018). In streams, turbulence drives the vertical transport of dissolved constituents through the water column (Hondzo, 1998;O'Connor & Hondzo, 2008;Tomasek et al., 2018) governing the rate at which reactive constituents (e.g., nitrate) can be assimilated and removed by the streambed (Grant, Gomez-Velez, et al., 2018). Stream turbulence also enhances the transport of dissolved and fine particulate material through the benthic biolayer, the upper 5 cm of the streambed where much of the microbial biomass, as well as nutrient and pollutant processing, is concentrated (Battin et al., 2008;Caruso et al., 2017;Dahm et al., 2002;Harvey et al., 2013;Knapp et al., 2017;Tomasek et al., 2018;Trauth et al., 2014;Zarnetske et al., 2011).
Only under highly idealized conditions is it possible to resolve the spatially and temporally complex advection pathways generated by dispersive mixing and turbulent diffusion. Due to its simplicity, tractability, and consistency with scaling approaches, a common alternative is to describe mass transport across the SWI and through the benthic biolayer as a horizontally and temporally (over turbulence time scales) averaged flux gradient process (Voermans et al., 2018): The variables J(y, t) (M L −2 T −1 ) and C s (y, t) (M L −3 ) are the vertical solute flux and interstitial concentration at depth y (which is 0 at the SWI and increases with depth into the streambed, Figure 1a) and time Figure 1. (a) An illustration of how water column turbulence can influence mass transport in the benthic biolayer. In this diagram, the benthic biolayer consists of a flat coarse-grained streambed subject to dispersive mixing and turbulent diffusion by a traveling pressure wave (dashed blue line), a mean velocity boundary layer that crosses the sediment-water interface (envelope of black arrows), and turbulence penetration (red eddies). The vertical mass flux J (y) arising from these phenomena is assumed to follow the flux gradient model (Equation 1a). (b) Turbulent mass transport across the SWI can be measured in the laboratory using closed systems, such as a stirred tank. Two-way coupling across the SWI is indicated by the circular arrows.
The use of Equations 1a and 1b to describe mixing in the benthic biolayer raises three questions. First, given that the mean and turbulent flow fields responsible for mixing across the benthic biolayer decrease with depth into the sediment (Breugem et al., 2006;Pokrajac & Manes, 2009;Roche et al., 2018;Voermans et al., 2017), what is the vertical structure of the effective diffusivity? Second, for a given vertical structure, how well does the flux gradient model (Equation 1a) describe solute transport through the streambed? Third, how do we extrapolate effective diffusivities measured in the laboratory to streams and coastal sediments? Here, we address these three questions in the context of dispersive mixing and turbulent diffusion across a flat sediment bed. Complementary efforts are underway to address mixing in the benthic biolayer by bedform pumping and bedform turnover (cf. Grant et al., 2020).
The paper is organized as follows. In section 2 we demonstrate, through the application of Duhamel's Theorem, that solute concentration in the interstitial fluids of the sediment bed can be represented as a convolution of solute concentration in the water column with the Green's function for mass transport in the streambed (Leij et al., 2000). This leads to a set of explicit solutions for the spatiotemporal evolution of concentrations in the water and sediment columns of a closed system. Notably, these solutions capture the two-way coupling of concentration evolution above and below the SWI, whereby mass transfer out of the streambed alters concentration in the overlying water column which, in turn, alters mass transfer into the streambed, and so on ( Figure 1b). We then derive four Green's functions for two choices of the lower-boundary condition (finite or semi-infinite sediment domain) and three functional forms of the diffusivity profile. In section 3 we demonstrate how this theory can be used to simulate unsteady mass transfer in a closed system, and in section 4 apply it to previously published measurements of turbulent mass transfer across a flat unconsolidated sediment bed in a well-stirred tank. We address the three questions raised above in section 5 and present our conclusions in section 6.

Duhamel's Solution for Turbulent Mixing in the Benthic Biolayer
Averaging over the time scale of turbulence and assuming sediment porosity does not change appreciably through the benthic biolayer (Knapp et al., 2017), mass conservation for a conservative solute in a horizontally uniform system takes the form In this study, we investigate how mass transfer across the SWI is influenced by the variation in effective diffusivity with depth: D eff ( y) = D eff,0 f (y), where D eff,0 is the effective diffusivity at the SWI (henceforth referred to as the surficial effective diffusivity) and f ( y) is a piecewise continuous function for which f(0) = 1. Adopting this functional form for the effective diffusivity, Equation 2a can be rewritten in dimensionless form where the new dependent variable, c s (-), incorporates the solute's initial concentration in the sediment and water columns (C s0 and C w0 , respectively): The constant a (L −1 ) is an inverse depth scale (whose definition will be shown to depend on the choice of diffusivity profile) and t T is a time constant for solute mixing in the benthic biolayer. Given the definition of dimensionless concentration (Equation 2c), the initial condition for Equation 2b becomes At the upper boundary (the SWI, y = 0) we require that the interstitial solute concentration equals the solute concentration in the overlying well-mixed water column, C s (y = 0, t) = C w (t). Expressed in dimensionless form the upper boundary condition becomes The Heaviside step function H t ð Þ (-) in Equation 3b ensures the upper boundary condition is 0 for t ≤ 0 (this detail becomes important for the application of Duhamel's Theorem below). Expression of the upper boundary condition in this way implies that mass transfer across the SWI is limited by the transport of solute within the streambed and not by mixing across the overlying concentration boundary layer (Grant, Azizian, et al., 2018;Grant, Gomez-Velez, et al., 2018). That is, the Biot number (the ratio of time scales for diffusive mixing in the streambed and mass transfer across the concentration boundary layer) is much greater than unity (Incropera et al., 2007).
One of two lower-boundary conditions can be selected, depending on whether the sediment bed is finite (Equation 3d) or semi-infinite (Equation 3e) in extent.
enforces a no-flux boundary condition at the normalized depth d b ¼ ad b (-) where d b is the depth of the sediment bed ( Figure 1). Equation 3e prescribes that, deep within the bed (y→∞), the interstitial concentration is maintained at its initial value.
As documented in the supporting information (Texts S1 and S2), by invoking Duhamel's Theorem (Perez Guerrero et al., 2013) the above system of equations can be solved for any time-varying solute concentration in the overlying water column, any piecewise continuous diffusivity profile, and either a finite or semi-infinite streambed. The solution is a convolution of the dimensionless water column concentration, c w t ð Þ, with a so-called Green's function (Leij et al., 2000;Myers, 1971), G y; t ð Þ (T −1 ), scaled here by the mixing time scale introduced earlier (Equation 2d), G y; t ð Þ ¼ t T G y; t ð Þ: According to Equation 4, dimensionless solute concentration in the interstitial fluid at any depth and time, c s y; t ð Þ, depends on the entire prior history of dimensionless solute concentration in the water column, c w t ð Þ, filtered through the Green's function, G y; t ð Þ. Green's function, in turn, is a fundamental solution to the diffusion equation (Equation 2b) that characterizes the response of solute concentration in the interstitial fluid of the streambed to an impulsive injection of mass at the SWI at is the Dirac delta function. The mathematical form of the Green's function depends on the vertical structure of the diffusivity profile, f y ð Þ, and the lower-boundary condition (either Equation 3d or 3e). Four Green's functions, corresponding to different combinations of the diffusivity profile and lower-boundary condition, are derived in section 2.3.

Two-Way Coupling Across the SWI in a Closed System
In typical applications of Duhamel's Theorem, the functional form of the non-homogeneous boundary condition (i.e., the water column concentration, c w t ð Þ) is stipulated in advance. In our case, however, the water and sediment concentrations are fully coupled through mass flux across the SWI. For a closed system, like the stirred tank illustrated in Figure 1b, the change of solute mass in the water column is equal to the rate of mass transfer across the SWI by dispersive mixing and turbulent diffusion: New variables appearing here include the interfacial area, A b , of the streambed and the height of the water column, h w . The streambed porosity, 0 < θ < 1, appears on the righthand side of Equation 5a to capture the abrupt change in area over which mass transport occurs above and below the SWI (Grant et al., 2012). Using the dimensionless variables introduced earlier, Equation 5a can be expressed as follows where the new dimensionless variable, h w (Equation 5c), is a scaled form of the water column depth: ). This two-way coupling can be solved exactly by manipulating the water and sediment mass balance equations in the Laplace domain. As demonstrated in the supporting information (Text S3), the result is a set of fully coupled solutions for solute concentration in the water and sediment columns of a closed system: , and e G is the Laplace transform of the dimensionless Green's function which, as noted earlier, depends on the diffusivity depth profile f y ð Þ and bottom boundary condition. A corresponding set of solutions can be derived without two-way coupling (whereby the diffusion equation's upper boundary condition is maintained at C s y ¼ 0; t ð Þ¼C w0 ): The inverse Laplace transforms in Equations 6a-6d were determined analytically or evaluated numerically using Gaussian quadrature (Graf, 2004).

Laplace Domain Solutions for the Green's Function
In Table 1 we present four Laplace domain solutions for the Green's function given three choices of diffusivity depth profile and two choices of bottom boundary condition (finite or semi-infinite sediment bed) (derivations in the supporting information, Text S4). This analysis therefore provides 16 different solution combinations for concentration in the water and sediment columns of a closed system with (Equations 6a and 6b) or without (Equations 6c and 6d) two-way coupling across the SWI.
Most laboratory (Grant et al., 2012;Marion & Zaramella, 2015;O'Connor & Harvey, 2008) and field (Wörman, 2000) studies of mixing across the SWI adopt the C Profile. However, several studies (Chandler et al., 2016;Nagaoka & Ohgaki, 1990;Roche et al., 2019) have shown that turbulent mixing in the sediment bed decays exponentially with depth, and a recent numerical modeling study concluded that the E Profile is consistent with experimental breakthrough curves measured in the laboratory and field (Bottacin-Busolin, 2019). The C2E Profile captures enhanced mixing at the top of the streambed by extending the surficial effective diffusivity D C2E eff; 0 to a depth y = ℓ t (L) below the SWI (c.f., the analytical model in Roche et al., 2019). The diffusivity profile declines exponentially below this depth, y > ℓ t . The dimensionless form of the enhanced mixing depth is defined as follows, ℓ t ¼ a C2E ℓ t .

Example of the Theory's Application to Mixing Across the SWI in a Stirred Tank
Consider an experiment similar to those described in section 4, in which mass M of a conservative solute is added to the interstitial fluid of the sediment bed in an otherwise solute-free stirred tank. Adopting the notation from Figure 1b, the initial interstitial concentration is C s0 = M/(d b A b θ). At time t = 0 the impeller is turned on, causing the concentration in the overlying water column to rise as solute is turbulently mixed out of the bed. Over time, the water and sediment concentrations will approach a final (well-mixed) equilibrium concentration, C eq : Within the context of our modeling framework, the temporal evolution of solute concentration from its initial (all solute mass in the sediment) to final (well-mixed) state depends on the diffusivity's depth profile, whether two-way coupling across the SWI is considered and whether the sediment bed is modeled as semi-infinite or finite.
If we adopt the C Profile, for example, the theory presented in section 2 leads to the following three solutions for solute concentration in the water column (Text S5): (1) a "null model" for an infinitely deep sediment bed without two-way coupling (Equation 8b); (2) an "infinite bed" model for an infinite sediment bed with two-way coupling (Equation 8c); and (3) a "finite bed" model for a finite sediment bed with two-way coupling (Equation 8d).
The concentration in the water column is unbounded for the null and infinite bed models, because their lower-boundary condition implies that an infinite mass of solute is stored in the sediments. The superscript "C" indicates the solutions are specific to the C Profile. The corresponding set of solutions for interstitial concentration in the sediment bed are as follows: The null model predicts that concentration in the water column increases unboundedly with the square root of time ( Similar patterns are evident for model-predicted concentration in the sediment (Figure 3b). Here we focus on the evolution of interstitial solute concentration in the upper portion of the bed, y C ¼ 0:05. The null model Note. The functions K 0 , K 1 , and K 2 are modified Bessel functions of the second kind. predicts a rapid decrease in concentration initially (as solute in the upper portion of the streambed mixes into the water column) followed by a gradual decline over time (as solute from deeper in the bed mixes upward The theoretical framework above was applied to an extensive set of previously published measurements of turbulent mixing of a conservative tracer (Rhodamine WT) across a flat SWI in a stirred tank (Chandler et al., 2016; "C16" henceforth). C16's experiments are notable because they (1) covered a range of bed shear velocities (u * = 0.01-0.04 m s −1 ), mean grain diameters (d g = 0.15-5.00 mm), and sediment permeabilities (K= 0.18-223 × 10 −10 m 2 ); and (2) simultaneously measured water and sediment column concentrations. The second feature permits a direct comparison of mixing parameters estimated from concentration data collected exclusively above or below the SWI.
C16's experimental methods, and the approach we used for parameter estimation and model performance evaluation, are briefly described here (see Text S6 for details). The sediment column consisted of randomly packed single-sized spherical glass spheres, with a depth d b = 0.2 m and a porosity θ= 0.38-0.39. In all experiments, the initial state was a Rhodamine WT saturated sediment bed (C s0 ≈ 100 ppb) and a Rhodamine WT-free water column (C w0 ≈ 0 ppb) (Table S1). Tracer concentrations were monitored fluorometrically in the water column and at five depths in the sediment column (y = 0.015, 0.049, 0.083, 0.117, 0.151 m) at 0.1 Hz over a period of hours to days. Diffusivity profile parameters were inferred by minimizing the root-mean-square error (RMSE) in nonlinear least squares regression between experimental data and model predictions. The corrected Akaike Information Criterion (AICc), which accounts for the trade-off between model fit and model complexity, was used to rank the performance of the three diffusivity profiles; the top-ranked model has the smallest AICc value (Aho et al., 2014). For model fitting, we used the "infinite bed" model specific to each profile and restricted the experimental time window to periods when Rhodamine WT concentration at the deepest sensor changed by <10% (Table S1). Parameter values for all three diffusivity profiles (C, E, C2E) were inferred from the water column data measured in 20 of C16's experiments (Tables S2, S3, and S5); six experiments were excluded due to missing data or other issues. E Profile parameters were separately estimated from sediment column data measured in the same 20 experiments (Table S4).

Experimental Evaluation of the C Profile
In comparison to C16's experimental data, there is significant bias in the concentrations predicted by the C Profile's infinite bed model (Figure 4a). In this figure, normalized Rhodamine WT measurements are plotted against the square root of dimensionless time, eff; 0 , varies by experiment depending on the inferred value of the effective diffusivity (Table S2); the other two parameters, h w and θ, vary negligibly. Thus, C16's data can be compared directly to a single model prediction for the time evolution of concentration in the water column (Equation 8c, upper graph in Figure 4a) and at two depths (15 and 151 mm) in the sediment column (Equation 9b, lower graph in Figure 4a). For clarity, measured Rhodamine WT concentrations at the three intermediate depths (4.9, 8.3, 11.7 cm) are not included in this figure; Figure 3. The influence of two-way coupling and finite bed depth on the evolution of solute concentration in the (a) water column and (b) interstitial fluids of the sediment bed, assuming solute is initially present only in the interstitial fluids of the sediment bed and the diffusivity profile is constant with depth (C Profile). The normalized bed depth was set to unity for these simulations, d b ¼ 1.
In C16's experiments, the concentration of Rhodamine WT in the water column is proportional to the cumulative mass of Rhodamine WT transferred from the sediment to the water column up to that point in time.
Thus, the C Profile model underestimates and overestimates mass transfer out of the sediment bed at short and long times, respectively (Figure 4a, upper graph). The underlying problem can be diagnosed by comparing model-generated and measured Rhodamine WT concentrations in the sediment bed (lower graph, Figure 4a). The C Profile model underestimates mixing in the surficial portion of the bed at short times (predicted concentrations exceed measured concentrations at 15 mm) and overestimates mixing deeper in the bed at later times (predicted concentrations less than measured concentrations at 151 mm).

Experimental Evaluation of the E Profile
Model bias is reduced substantially when the effective diffusivity decays exponentially with depth ( Figure 4b). The E Profile's infinite bed model was constructed by substituting its Green's function (Equation T3 in Table 1) into the water and sediment solutions for a closed system with two-way coupling (Equations 6a and 6b): The superscript or subscript "E" indicates that these solutions are specific to the E Profile. Because the E Profile model has two unknown parameters (D E eff; 0 and a E ), there is no longer a single curve against which all of C16's data can be compared (as was the case for the C Profile in Figure 4a). For the E Profile (and C2E Profile described below) model-data comparisons must be conducted on an experiment-by-experiment basis. In Figure 4b the comparison is performed for a typical C16 experiment (ID #20110613, red points in the figure). The E Profile infinite bed model reproduces water column measurements of Rhodamine WT concentration (top graph, Figure 4b), although some bias is evident for < 0:05. Sediment concentrations predicted by the E Profile's infinite bed model capture the fast and slow mixing out of the shallow and deep portions of the sediment bed, respectively (lower graph, Figure 4b). For consistency, E Profile model predictions are plotted against the same abscissa, , used for the C Profile model in Figure 4a. The experimentspecific transformation from t C to t E was determined by (1) fitting Equation 10a to Rhodamine WT concentration measured in the water column during experiment ID #20110613 (red points in Figure 4b, a E = 50 ± 1.2 m −1 , D E eff; 0 ¼ 5:6 ± 0:5 ð Þ× 10 −6 m 2 s −1 ), and (2) substituting these inferred parameter values, together with experiment-specific values of the dimensionless water depth ( h w ¼ a E h w =θ ¼ 33 ± 14 ) and the C Profile's effective diffusivity (D C eff; 0 ¼ 3:4 ± 0:2 ð Þ× 10 −7 m 2 s −1 ), into the time transformation,

Experimental Evaluation of the C2E Profile
The C2E Profile's infinite bed model is a near-perfect representation of C16's water column measurements (top graph, Figure 4c). However, compared to the E Profile, the C2E infinite bed model systematically underestimates mixing in the streambed, especially at 15 mm below the SWI (compare lower graphs in Figures 4b  and 4c). These model predictions were constructed by substituting the C2E's Green's function (Equations T4a-T4c in Table 1) into Equations 6a and 6b. C2E model predictions were then plotted against ffiffiffiffi ffi t C p following a modification of the two-step procedure outlined in section 4.3: (1) fitting the C2E model to Rhodamine WT concentration measured in the water column during experiment ID #20110613 (D C2E eff; 0 ¼ 1:5 ± 0:07 ð Þ × 10 −6 m 2 s −1 , ℓ t = 0.04 ± 0.002 m, a C2E = 66 ± 3.1 m −1 ); and (2) substituting these inferred values, along with experiment-specific values of the dimensionless water depth ( h w ¼ a C2E h w =θ ¼ 44 ± 19 ) and the C Profile's effective diffusivity (D C eff; 0 ¼ 3:4 ± 0:2 ð Þ× 10 −7 m 2 s −1 ), into the time transformation, The superscript or subscript "C2E" indicates the variables are specific to the C2E Profile.

Discussion
Here we discuss all 20 of C16's experiments, with the goal of answering the three questions raised in section 1.

How Is the Effective Diffusivity Structured Vertically?
C16's data set allows us to quantitatively compare the performance of the three diffusivity profiles (C, E, and C2E) over a permeability Reynolds number range (0.2-4.3) that spans dispersive mixing and turbulent diffusive regimes. Our earlier conclusion that the E Profile represents a substantial improvement over the C Profile (based on a comparison to Rhodamine WT measurements from Experiment ID #20110613, section 4) extends to the rest of C16's experiments as well ( Figure 5). In all cases, the E Profile's infinite bed model captures a larger fraction of the variance in C16's water column measurements (R 2 > 99.5%, Figure 5a) and has substantially smaller RMSE values (Figure 5b). The E Profile's AICc is also >10 units lower than the C Profile's AICc (Figure 5c) implying that the former model is preferred (Aho et al., 2014;Weijs & Ruddell, 2020).
The C2E Profile's infinite bed model also performs well. Compared to the E Profile, the C2E model has consistently lower RMSE and AICc values (compare crosses and red circles in Figures 5b and 5c) and a slightly improved coefficient of determination (R 2 > 99.8%, Figure 5a). However, these improvements come at the cost of a new parameter (the C2E's inverse depth scale, a C2E ) whose inferred values are poorly constrained (coefficient of variation in excess of 40 for some experiments, see Table S7), highly variable (varying over 1,000-fold, Figure 6a) and, in some cases, not physically meaningful (e.g., the depth over which the effective diffusivity decays is 1/a C2E ≈ 20 μm for the largest value of a C2E indicated in Figure 6a). Inferred values of the C2E Profile's other two parameters (effective diffusivity, D C2E eff; 0 , and depth of constant mixing, ℓ t ) are strongly correlated (R 2 = 0.93 and 0.76) with the effective diffusivity, D E eff; 0 , and decay scale, 1/a E , inferred from the E Profile (Figures 6b and 6c). In summary, of the three profiles trialed in this study, the E Profile appears to be the most parsimonious description of the effective diffusivity's vertical structure.

Is the Flux Gradient Diffusive Model an Accurate Representation of Turbulent Solute Transport Through the Streambed?
While the E Profile's infinite bed model captures a large fraction of the variance in C16's water column measurements (R 2 > 99.5%, Figure 5a), this assessment is based on the same data set that was used for model calibration. A more rigorous test can be stated as follows: Are the same E Profile parameter values obtained when the model is optimized with water column measurements versus when the model is optimized with sediment column measurements? Put another way, can the evolution of solute concentration in the interstitial fluid of the sediment bed be inferred from the evolution of solute concentration in the water column, and vice versa?
The answer is a qualified "yes." Effective diffusivities estimated from C16's water and sediment column data are strongly correlated over a 1,000-fold change in magnitude (Figure 7a, Pearson correlation coefficient, R = 0.867). Values of the inverse depth scale inferred from the water and sediment column data are much less variable and not significantly correlated (Figure 7b), but their respective log means (a E = 10 1.61 ± 0.18 and a E = 10 1.70 ± 0.08 m −1 , respectively) are equal within error (and consistent with the inverse depth scale reported by C16, a = 55 = 10 1.74 m −1 , estimated by dividing the sediment column into a series of layers and computing, with the C Profile model, diffusivities for each layer separately). This inverse depth scale corresponds to an "e-folding depth" (i.e., the depth at which the E Profile's effective diffusivity declines to 1/e ≈ 0.37 of its surficial value) of approximately 2 cm. This depth scale comports with previous field and laboratory estimates for the thickness of the benthic biolayer (2-5 cm) (Battin et al., 2008;Caruso et al., 2017;Dahm et al., 2002;Harvey et al., 2013;Kessler et al., 2013;Knapp et al., 2017;Krause et al., 2017;Tomasek et al., 2018;Trauth et al., 2014;Zarnetske et al., 2011). Hence, turbulent mixing may represent a universal mechanism for delivering the nutrients and energy needed for microbial growth in the benthic biolayer. More generally, the similar parameter values obtained from the water and sediment measurements support the claim that the gradient-flux diffusive model (Equation 1a) is a reasonable representation of turbulent solute transport across the SWI, provided that the diffusivity declines exponentially with depth into the streambed.

Water Resources Research
While effective diffusivities inferred from data collected above and below the SWI are strongly correlated, some bias is evident when the permeability Reynolds number exceeds the threshold for a fully turbulent SWI, Re K > 2 (Voermans et al., 2017) (Figure 7a). One possible explanation is C16's use of fluorometric measurements of tracer at a point in the sediment bed, which contrasts with the modeling assumption that solute concentrations are horizontally averaged (section 2). Chandler (2012) describes how measured concentration in the sediment was sensitive to fluorometer location and that these differences were consistent over time; that is, tracer appeared to mix out of the bed faster on one side of the tank than on the other (Chandler, 2012, p. 173). The authors also document distinct and persistent patterns of mean flow velocity within the tank (Chandler, 2012, p. 118), which would lead to heterogeneous turbulence intensities and corresponding heterogeneous efflux across the SWI at high permeability Reynolds numbers. To the extent that C16's point measurements are not equal to horizontally averaged concentrations, the effective diffusivities inferred from these data may be nonrepresentative. Indeed, Chandler (2012) noted an order of magnitude discrepancy in the time scale over which interstitial concentration declined on opposite sides of the stirred  tank-a discrepancy that could induce order-of-magnitude inaccuracy in estimated sediment side diffusivities (the scale of disagreement seen in Figure 7a for Re K > 2).
An alternative explanation is that the flux gradient model (Equation 1a) is an imperfect descriptor of turbulent mass transfer through the sediment at high permeability Reynolds number-a conclusion supported by the slight reduction in E Profile model performance (i.e., higher RMSE and AICc values and lower R 2 values) with increasing permeability Reynolds number above, Re K > 2 ( Figure 5). The flux gradient description of dispersive mixing across other types of "porous" boundaries, such as vegetation canopies, is known to break down; for example, significant dispersive momentum flux can occur in regions of weak velocity gradient (Poggi & Katul, 2008), the analog of a concentration gradient for mass transfer. However, our estimates of diffusivity above and below the SWI in the dispersive regime (0.01 < Re K < 2) are in general concordance ( Figure 7a). It should also be noted, even at the highest permeability Reynolds numbers evaluated here, the flux gradient diffusive model still explains an overwhelming fraction of the variance in water concentration measurements (R 2 > 0.995)-provided that the diffusivity decays exponentially with depth ( Figure 5a). Therefore, the flux gradient diffusive model is a reasonable representation of turbulent mass transfer in the interstitial fluids of the sediment bed, provided that the vertical structure of the effective diffusivity is correctly specified (e.g., with the E Profile).

Can Laboratory Measurements of Turbulent Mixing Across the SWI Be Extrapolated to Stream and Coastal Sediments?
Translation of our results to the field requires scaling relationships from which the E Profile's two parameters -the surficial effective diffusivity and inverse decay scale-can be estimated. Many studies have reported that effective diffusivities (inferred by fitting the C Profile's null model to flume measurements of turbulent mixing across flat sediment beds) exhibit a quadratic dependence on the permeability Reynolds number, D C eff; 0 ∝ Re 2 K (Grant et al., 2012;O'Connor & Harvey, 2008;Richardson & Parr, 1988;Voermans et al., 2018). The permeability Reynolds number is calculated from the shear velocity, u * , sediment bed permeability, K, and the kinematic viscosity of water, υ. Permeability can be estimated from the grain diameter and porosity of unconsolidated sediments (e.g., using the Kozeny-Carmen equation; Kamann et al., 2007) while the kinematic viscosity of water is determined primarily by temperature (Rumble, 2019). Several methods are available for measuring shear velocity (cf., Johnson & Cowen, 2017) including a force balance approach that provides a spatially averaged value from the stream's hydraulic radius (or depth, h w , for wide streams) and slope, S (-): u * ¼ ffiffiffiffiffiffiffiffiffiffi gh w S p where g is the gravitational acceleration. Thus, if the E Profile's two parameters can be expressed in terms of the permeability Reynolds number, such relationships would allow translation of laboratory measurements to field applications.
When our inferred effective diffusivities are plotted against the permeability Reynolds number, a significant change in slope and intercept (as represented by nonoverlapping 95% confidence intervals) is evident around Re K = 1 (Figure 8a): The power law exponent for the surficial effective diffusivity spans the quadratic dependence noted above, declining from 2.53 ± 0.17 in the dispersive mixing regime (Re K < 1) to 0.99 ± 0.15 in the turbulent diffusive regime (Re K > 1). Note that this scaling relationship is normalized by the molecular diffusion coefficient for Rhodamine WT in water at 21°C (D m = 2.9 × 10 −10 m 2 s −1 ) (Chandler, 2012). Therefore, over the range of permeability Reynolds number captured in C16's study, turbulence enhances mixing by 10 2 to 10 5 above that expected for molecular diffusion alone (Figure 8a).
Equation 11b is our scaling relationship for the E Profile's inverse decay scale, here normalized by a rough estimate for the thickness of the benthic biolayer (Knapp et al., 2017): ℓ B = 2 cm (Figure 8b).
log 10 a E ℓ B ¼ 0:04 ± 0:01; Re K < 1 0 ± 0:02 ð Þ− 0:32 ± 0:08 ð Þ× log 10 Re K ; R 2 ¼ 0:18; The implied turbulent mixing depth (1/a E ) transitions from being roughly equal to the benthic biolayer thickness in the dispersive mixing regime (Re K < 1, 1/a E ≈ ℓ B ) to a weak inverse dependence on the permeability Reynolds number in the turbulent diffusive regime (Re K > 1, a E ℓ B ∝ Re −0:32 ± 0:08 K ). The minimum turbulent mixing depth, 1/a E ≥ ℓ B ≈ 2 cm, is between 4 and 130 times the diameter of the glass spheres that make up the sediment bed (0.150 ≤ d g ≤ 5 mm) and about 20 times the estimated thickness of the Brinkman Layer (the region of enhanced mean velocity at the top of the sediment) δ b ≈ 20 ffiffiffi ffi K p ¼ 0:9 mm (Voermans et al., 2017). It is also about tenfold less than the sediment bed depth (d b = 20 cm), implying that 1/a E is not a proxy for bed depth. These comparisons raise the following question: What is the physical interpretation of the inverse decay scale? Based on a model for mass exchange across the SWI by turbulent dispersive mixing, Higashino et al. (2009) reported that, at depths of around 2 cm and for shear velocities on the lower end of the range employed by C16, the root-mean-square vertical velocity of the interstitial pore fluids are >10% of their value at the SWI (Higashino et al., 2009; Figure 3). Thus, when Re K < 1, the inverse decay scale likely represents the surficial depth over which mass is vigorously mixed by turbulence-induced dispersive mixing. This physical interpretation of a −1 E also holds for cases when turbulent diffusion controls solute transport (i.e., Re K > 1). In this regime, regions of the streambed where solute is rapidly mixed correspond to regions of elevated turbulent shear stresses (Roche et al., 2018). Because turbulent shear stresses propagate deeper into the streambed at elevated Re K (Voermans et al., 2017), rapid solute mixing is expected to extend deeper into the streambed as Re K increases. This expectation is in direct agreement with our observations of a E decreasing with Re K when Re K > 1 (Figure 8b).
Because the above scaling relationships (Equations 11a and 11b) are based on C16's stirred tank experiments, they may not apply to all turbulent environmental flows. For example, while C16's permeability Reynolds number range includes dispersive mixing and turbulent diffusive regimes, their shear velocities are on the low side (u * = 0.01 to 0.04 m s −1 ) for headwater streams (Hall et al., 2009). Field validation of Equations 11a and 11b will be an important next step.

Conclusions
In this paper, we developed and tested a rigorous one-dimensional modeling framework, based on Duhamel's Theorem, for predicting mass transfer across the SWI and through the benthic biolayer of a turbulent stream. The framework allows for depth-varying diffusivity profiles and encodes two-way coupling across the SWI. The theory is applied to previously published measurements of turbulent mixing across a flat sediment bed in a stirred tank (Chandler et al., 2016) to evaluate the performance of three diffusivity depth profiles (C, E, and C2E Profiles). Key findings include (1) the flux gradient diffusive model is a reasonable representation of turbulent mass transfer across the SWI and in the sediment bed, provided that the vertical structure of the effective diffusivity is correctly specified; (2) The experimental data are consistent with an exponentially declining diffusivity profile (i.e., the E Profile); (3) values of the E Profile's two parameters (surficial effective diffusivity at the SWI, D E eff; 0 , and decay depth scale, a E ) vary with the permeability Reynolds number, Re K , providing a direct link between lab results and field-scale applications; (4) the E Profile's dependence on the permeability Reynolds number changes abruptly at Re K = 1, reflecting different modes of mixing below (turbulent dispersive mixing) and above (turbulent diffusion) this threshold; and (5) the effective diffusivity's e-folding depth is concordant with field and laboratory measurements of the benthic biolayer thickness. Therefore, turbulent mixing across the SWI may serve as a universal transport mechanism, supplying the nutrient and energy fluxes needed to sustain microbial growth, and nutrient processing, in the benthic biolayer of stream and coastal sediments.