A Unifying Model for Turbulent Hyporheic Mass Flux Under a Wide Range of Near‐Bed Hydrodynamic Conditions

Existing models for estimating hyporheic solute mass flux often require numerous parameters related to flow, bed, and channel characteristics, which are frequently unavailable. We performed a meta‐analysis on existing data set, enhanced with high Reynolds number cases from a validated Computational Fluid Dynamics model, to identify key parameters influencing effective diffusivity at the sediment water interface. We applied multiple linear regression to generate empirical models for predicting eddy diffusivity. To simplify this, we developed two single‐parameter models using either a roughness or permeability‐based Reynolds number. These models were validated against existing models and literature data. The model using roughness Reynolds number is easy to use and can provide an estimate of the mass transfer coefficient for solutes like dissolved oxygen, particularly in scenarios where detailed bed characteristics such as permeability might not be readily available.


Introduction
At the sediment-water interface (SWI), sediment oxygen demand (SOD) functions as an oxygen removal flux, transporting dissolved oxygen (DO) from the water to the sediment.This balances the penetration of DO caused by near-bed turbulence with the DO consumed by sediment and benthic chemical processes within the bed (Boudreau & Jorgensen, 2001;Gundersen & Jorgensen, 1990;Jørgensen & Revsbech, 1985;Mackenthun & Stefan, 1998;Nezu & Nakagawa, 1993).Accurate modeling of DO dynamics and oxygen mass transfer at the SWI is crucial for understanding nutrient cycling in riverine systems (Boano et al., 2014;Motta et al., 2010;Waterman et al., 2009Waterman et al., , 2011Waterman et al., , 2016) ) This motivated our present work on unified models to predict the mass flux of solutes like oxygen under different flow and bed conditions.
In smooth-wall hyporheic flows, the mass transfer coefficient (K L ) is a function of diffusivity (D eff ) and diffusive layer thickness (δ DL ): K L = D eff /δ DL at the SWI.K L is affected by flow shear velocity (u * ), Reynolds number, bed roughness, and the momentum exchange due to hyporheic flow.O'Connor et al. (2009) summarized data for K L from the literature collected over hydrodynamically smooth beds (Arega & Lee, 2005;Dade, 1993;Hondzo et al., 2005;O'Connor & Hondzo, 2008;Shaw & Hanratty, 1977;Steinberger & Hondzo, 1999) and expressed the dimensionless mass transfer coefficient (K L+ ) as a function of a temperature-dependent Schmidt number (S c ): , with ranges for α and β being 0.052 to 0.164 and 0.704 to 0.67, respectively.K L+ has Reynolds number dependence at low and moderate Reynolds number flows (Grant et al., 2018;Grant & Marusic, 2011) and reaches a self-similar plateau value for large enough Reynolds numbers (Shaw & Hanratty, 1977;Steinberger & Hondzo, 1999).
High roughness height (Lai et al., 1994;Nikora, Goring, & Biggs, 2002, Nikora, Koll, et al., 2002), bed permeability (Jiménez, 2004;Perry et al., 1969;Raupach et al., 1991;Richardson & Parr, 1988;Wu et al., 2019), and bed forms (Elliott & Brooks, 1997;Marion et al., 2002;Packman et al., 2000Packman et al., , 2004;;Packman & MacKay, 2003;Rehg et al., 2005;Ren & Packman, 2004;Tonina & Buffington, 2007) can also enhance hyporheic momentum exchange and thus the mass transfer by increasing the shear stress at the SWI.Han et al. (2018) summarized K L+ values from experimental studies over rough beds and found that rough bed K L+ values can be up to two orders of magnitude higher than those over smooth beds (Han et al., 2018;Inoue & Nakamura, 2011;Nagaoka & Ohgaki, 1990).D. J. O'Connor (1984) studied the transfer coefficient for openchannel flows with smooth and rough beds and proposed an analytical equation for the transitional regime.δ DL values are difficult to be estimated for the case of rough wall.Thus, different length scales may be used as flow characteristic length scales at the SWI of rough walls.Nagaoka and Ohgaki (1990) adopted a pore scale restricted mixing length B = 2ϕ 2 3(1 ϕ) D, D is particle diameter or D 50 if particle size is not uniform in the bed, for the estimation of D eff at the SWI.Manes et al. (2012) and Voermans et al. (2017) proposed flow length scales associated with characteristic turbulent eddy size across the SWI δ * p ) and the depth of turbulent shear penetration in the bed (δ p ) respectively.Another relevant length scale is the vertical location of the inflection point in the mean velocity profile δ (where dU 2 dz 2 = 0 and d < U> dz is maximum) which also corresponds to the position of the SWI defined with respect to the top of the sediments (Voermans et al., 2017).
In our work, we focused on extending the analysis previously done for the scaling and modeling of effective diffusivity, estimating mass transfer coefficient, and developing new simple relationships for the estimation of these parameters in practice.We used a computationally efficient method, Improved Delayed Detached Eddy Simulations (IDDES), to extend the range of previously reported Re k and Re * in ranges that measurements are not available to date.The numerical results were compared and validated against data from the literature.We also performed a reanalysis and MLR using data from the literature combined with our numerical results to reexamine the relationships between D eff /D m and inner and outer-/bulk parameters of the bed and the flow.Finally, novel unifying single-parameter models for the prediction of mass transfer coefficient using roughness or permeability scales were proposed, based on Re * and Re k .The proposed models can accurately predict data from field and laboratory conditions from the literature.The developed model can be applied as a boundary condition for more detailed biochemical models that could potentially compute the transport of solutes such as oxygen within the sediment, particularly under mass-transfer limited condition.

Definitions
The physical parameters influencing the mass transfer flux of a solute at the SWI (J S ) is shown in Figure 1a.These parameters include the bulk and near-bed hydrodynamics (U b , u * , k s , H w ) as well as the SWI and porous sediment bed characteristics (K, ϕ, H b ).The total/effective mass flux comprises molecular J MD S ) , dispersive J DIS S ) , and turbulent J T S ) fluxes, which can be expressed as (Voermans et al., 2018b): The quantities mentioned above utilize Reynolds and spatial decompositions: ψ = ψ + ψ′ and ψ = < ψ > + ψ, in which the variable ψ is represented by overbarred, primed, bracketed, and tilded quantities.These indicate the time-averaged, time-fluctuating, spatially averaged, and spatial-fluctuating quantities, respectively (Lopez & Garcia, 1997;López & García, 1998, 2001;Nikora, Koll, et al., 2002, 2007;Voermans et al., 2018b).In this study, we performed IDDES simulations on theoretical cases, employing surrogate beds composed of closely packed, monodisperse spheres, as in previous research examining hyporheic mass exchange (Manes et al., 2009;Stoesser et al., 2007;Wu et al., 2019) (Figures 1a and 1b). Figure 1b illustrates typical normalized instantaneous streamwise velocity (u x ) normalized by u * , while mean streamwise velocity profile is depicted in Figure 1c.The inflection point (δ) (Figure 1d) in the mean streamwise velocity profile signifies the extrusion of bed-penetrating eddies and the virtual "interface" of the bed (Manes et al., 2012;Voermans et al., 2018b).The length scale (δ) is preferred over penetration length scale (δ p ) and total height (δ p* ) (all depicted in Figures 1c-1e) due to its independence from the need for endoscopic measurement techniques (Blois et al., 2014) or refractive index matching (Voermans et al., 2017), even though porous bed permeability still influences δ.Later in the paper, we capitalize on the benefits of δ to create an easy-to-use predictor for hyporheic mass exchange under a variety of flow and bed conditions.The values and comparison of different length scales with Voermans et al. (2018b) are summarized in SI6 in Supporting Information S1.
Figure 1e shows typical shear stresses as they are computed by the deployed numerical method: where τ ν = μd < u >/dz is the viscous stresses, τ RS is the sum of resolved and modeled Reynolds stresses, and τ form is the form-induced stresses (Nikora et al., 2007).The τ total profile can be approximated as (Garcia, 2008): with τ total = 0 at z = H w .τ total (z bed ) is estimated as the integral of viscous and drag forces over the top hemispheres (A is the area of the bed above z = 0).

Computational Fluid Dynamics Solver
Using the Detached Eddy Simulation (DES) approach, we expanded our data set via 3D hydrodynamic simulations of hyporheic boundary layer flows with OpenFOAM.We employed the PIMPLE algorithm for incompressible 3D Navier-Stokes equations and applied the Spalart-Allmaras Improved Delayed Detached Eddy Simulation (SA-IDDES) for modeling ν eff .The equations of motion were solved assuming that the flows were both incompressible and isothermal.This model was validated using benchmark hyporheic flow cases, compared against Smagorinsky LES model simulations and literature data (Lian et al., 2021;Manes et al., 2009).We performed detailed analysis to establish result independence considering domain size, boundary conditions, mesh resolutions, and integration time (see SI1 in Supporting Information S1 for details).(Voermans et al., 2017).Thus, the similarity relations in Ghisalberti (2009)

Hyporheic Mass Flux and Effective Diffusivity
The effective diffusivity (D eff ) can be used to parameterize hyporheic exchange under the assumption that the mass transfer to the sediment bed can be modeled by Fick's First and Second Law, which can be expressed for homogeneous porous mediums as (O'Connor et al., 2009;O'Connor & Harvey, 2008): The mass transfer flux at the SWI which is a boundary conditions for the above equation is defined as: In Figure 1a, we defined the mass transfer flux (J S ) which under equilibrium should balance the solute consumed within the sediment bed by chemical processes (C w , C s are solute concentration in water and sediment).This mass transfer flux includes the effect of molecular (D m ), dispersive (D dis ) and turbulent (D t ) diffusivities (Voermans et al., 2018b).For the case of smooth beds, D eff /δ DL are typically used for the computation of the flux.In our

Parameterization for D eff
We followed a similar approach as O'Connor and Harvey ( 2008) and Grant et al. (2012) for parameterizing the effective diffusivity (D eff ).Specifically, we used Buckingham's Pi theorem to create dimensionless groupings of the controlling independent variables, as demonstrated in O' Connor and Harvey (2008) and Grant et al. (2012): In the equation above, there are nine variables (n = 9) and two primary dimensions (L, T ), m = 2. Consequently, a predictive equation for effective diffusivity should have a maximum of seven (n-m = 7) non-dimensional groups.
For our analysis, we propose that a normalized D eff /D m can be predicted as a function of the following dimensionless parameters: Equation 7 has the potential to account for bulk and near sediment bed hydrodynamic effects (Re bulk , Re * ), SWI exchange (Re k , Pe k , ϕ), and flume facility/computational-domain size dependencies (Re H w ,Re Hb ) .
Assuming a power law model (Grant et al., 2012) for modeling the dependence of the equation on the above parameters, we can write: We used MLR to develop a model based on the available data from the literature.This data set integrates field and flume data from earlier studies (Grant et al., 2012;O'Connor & Harvey, 2008;Voermans et al., 2018b) along with our numerical results.The data set details are supplied in SI3 in Supporting Information S1 and this is the foundation for model development.In testing the possible models, we experimented with 255 different combinations of seven parameters, along with a constant α.The procedure of MLR necessitates that the dependent variable be defined as functions of a group of independent variables, which must not be highly correlated.Thus, it was crucial to initially examine the linearity between dependent and independent variables and then to verify the correlation between chosen independent variables.
Following the approach by Grant et al. (2012), we applied variance inflation factor (VIF) to remove combinations with increased correlation (Miles, 2014) and the Akaike information criterion (AIC) as a selection standard (Aho et al., 2014;Akaike, 1974;Sakamoto et al., 1986).See SI4 in Supporting Information S1 for detailed information.Figures 3a-3c show the fitness of the 3 best performing models: M1,M2 and M3.In an effort to develop a simple single-parameter model for all different bed and flow characteristics, D eff /D m versus every single of the 6 dimensionless numbers considered in our MLR analysis are plotted in Figures 3d-3i together with the corresponding R 2 values.It is shown that Re k and Re * show the best coefficient of determination, which will be explored as single parameters to develop empirical model for the prediction of the mass transfer coefficient.

A Unifying Model for the Hyporheic Mass Transfer
Using Re * and Re k , we were able to scale the mass transfer coefficient K L , forming a new relationship across diverse roughness and permeability conditions.To address the problem of defining the diffusive layer (δ DL ) in rough wall situations, we introduced a unifying transfer coefficient, KL = D eff /δ and K+ L = KL / u * , applicable to both hydrodynamically smooth and rough cases.This is done using the inflection point δ which can be estimated for all conditions using the equations in Figure 2.
The functions we proposed, represented as K+ L = αRe β i / (Re , where Re i = Re * or Re k , take inspiration from the sediment entertainment function by Garcia and Parker (1991).These functions describe how K+ L starts from zero and increases with Re * and Re k until it plateaus at higher Re * and Re k , aligning with self-similar plateau values proposed by Shaw and Hanratty (1977) and Steinberger and Hondzo (1999).Mass transfer and sediment entrainment are both induced by turbulence and momentum at the SWI (Higashino & Stefan, 2008;Lopez & Garcia, 1997;López & García, 1998).In the Garcia and Parker (1991) entrainment function, as the flow becomes more turbulent the near bed concentration of sediment reaches an asymptotic value (C b ∼ 0.3) which is the highest possible sediment concentration beyond which the high sediment concentration damps turbulence.Similarly, the normalized mass transfer coefficient, K + L , approaches a plateau value for high Re flows.This analogy is the main reason for adopting the same form for the single-parameter models proposed here.
We used data from O' Connor et al. (2009), Han et al. (2018) and Voermans et al. (2018b), Voermans et al. (2018bata set that include smooth to fully rough beds and low to high permeabilities (Elliott & Brooks, 1997;Marion et al., 2002;Nagaoka & Ohgaki, 1990;O'Connor et al., 2009;O'Connor & Hondzo, 2008;Packman et al., 2004;Steinberger & Hondzo, 1999;Tonina & Buffington, 2007;Voermans et al., 2017).All data can be found in Table S5 in SI5 in Supporting Information S1.Following regression analysis, we derived two equations, with results displayed in Figures 4a and 4b: For comparison, K+ L values are plotted using yellow cross symbols based on the effective diffusivity profile D eff estimated by Grant, Gomez-Velez, et al. (2020)) using experimental observations by Chandler et al. (2016) (Figures 4a and 4b).Numerical simulations applied to develop the unifying mass transfer models assumed a stationary/non-movable bed without considering effects of entrainment and resuspension on mass transfer (Teitelbaum et al., 2021;Wolke et al., 2019), but the unifying mass transfer models have been compared with data set from the literature that includes scenarios with movable beds/bedforms.
To validate the accuracy of our model, we tested it against the data set used for MLR analysis and compared with the performance of the zonal model (VM) by Voermans et al. (2018b), which is described in Equation 10.Two unified equations proposed in this study and VM are evaluated in Figure 4c, where δ values were estimated from closure equations depicted in Figure 2 with available D m , u * , k s , and K values in the data set.Note that the data set used for the regression in VM is the same data set used for MLR analysis and model evaluation.A similar zonal model has been introduced by Grant, Gomez-Velez, et al. (2020) for which the D eff /D m follows different linear equations that change slope and intercept in the dispersive (Re k < 1)and turbulent (Re k > 1) regimes respectively.
In Figure 4c, the roughness-based Reynolds number model demonstrates a higher R 2 value (0.71) than the permeability-based Reynolds number model (0.42).Voermans et al. (2018b)'s zonal model shows an R 2 value of 0.74, quite close to the Re * -based model, despite being derived from the validation data set.Lastly, the three best performing models derived using MLR (M1, M2 and M3) are also tested against the same data set (see Figure 4d).As expected, the multi-parameter models proposed here outperform both the single-parameter models and the zonal model by Voermans et al. (2018b), with R 2 values of 0.88, 0.80, and 0.81 respectively.
An advantage of the unifying model based on Re * is the fact that it is solely based on the Nikuradse roughness height rather than permeability, which can potentially be more challenging parameter to estimate, that is, the use of laboratory tests.In fact, even if compared with the prediction by Voermans et al. (2018b) zonal model (mean absolute error is 1.05), the MAE is 1.19 for the Re * -based model while the MAE for the Re k -based model is 1.66.Finally, the 3 MLR models have typically smaller errors (0.73, 0.97, 0.94 respectively); however, they introduce additional complexity and in some cases, that is, D eff /D m = Re 1.076 k Re 1.038 H b , parameters like the thickness of the bed permeable layer (H b ) in Re H b have more importance when we study the hyporheic exchange at a laboratory setting using shallow test flumes or require significant field work to estimate the elevation of any impervious bedrock layer.

Applications and Relevant Discussions of the Single-Parameter Model
The proposed Re * -based model outlined earlier requires minimal data inputs.These include some parameter estimations.
• Nikuradse roughness (k s ) can be determined based on characteristic bed diameter (D s ) as k s = α s D s for user selection.A comprehensive list of α s values can be found in Engelund (1970) and Garcia (2008)

Conclusion
A meta-analysis of pre-existing data sets, supplemented by high Reynolds number Computational Fluid Dynamics (CFD) results, was conducted to examine the scaling parameters that influence mass transfer in hyporheic zones.Using this enhanced data set, multiple linear regression was utilized to create multi-parameter predictive models for effective diffusivity in turbulent hyporheic flows.A novel unifying model was then introduced, aiming to estimate the mass transfer coefficient using the roughness height rather than bed permeability.This newly developed model underwent validation through comparisons with other models and existing literature data.It is designed to serve as a user-friendly tool that can provide essential data for estimating the mass transfer coefficient for solutes such as dissolved oxygen, particularly in scenarios where detailed bed characteristics might not be readily available.

Figure 1 .
Figure 1.(a) Schematic plot of the problem and definition of quantities, (b) instantaneous streamwise velocity (u x ) normalized by u * over the surrogate bed, (c) mean streamwise velocity profile (U x ) normalized by U b , (d) mean streamwise velocity gradient, (e) shear stress distribution.
Figures 2a-2d summarize the comparisons of one of the simulations performed herein in dashed lines for the mean velocity (U x ), Reynolds stresses ̅̅̅̅̅̅̅̅̅ u′w′ √ / u * ) and the variances (σ u /u * and σ w /u * ) with measurements performed by Manes et al. (2009) for five-layer low-Re case (Re k = 31.(b)and the LES results by Lian et al. (2021) for Re k = 24.2.Flows over highly permeable boundaries share the similarity relations introduced by Ghisalberti (2009): σ u /u * ∼ 1.8 and σ w /u * ∼ 1.1 at SWI are plotted for comparison.Figures 2e and 2f show the normalized δu * /ν thickness as it compares with the corresponding data in Voermans et al. (2018b) as functions of Re k and Re * , where Re * = k s u * /ν with k s = 2.5D based on roughness function for current flow and bed setup.δ is approximately ∼0.143 × 2.5D = 0.36D which is reasonably close to the 0.3 value reported by Voermans et al. (2017) (Note that δ = δ p* δ p ). Figures 2g and 2h show the mixing-length (<L m > = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ <u′w′ >/ (dU/ dz) 2 √ ) as functions of Re k and Re * .The predictions are reasonably close to those by Voermans et al. (2018b).SI2 in Supporting Information S1 summarizes the setup of all the cases in the present study.We aim to expand from laminar and transition regimes in existing data set (Re k ∼ O (10) and Re * ∼ O (1, 000)) to fully turbulent regime at higher Reynolds number (Re k ∼ O (1, 000) and Re * ∼ O (10, 000)).

Figure 2 .
Figure 2. (a), (b), (c) and (d) Typical CFD results and comparison against data from the literature for U x /U b , ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ <u′w′> √ / u * , σ u /u * and σ w /u * .(e), (f), (g) and (h) fitting relations of length scale δ and mixing length < L m > with Re k and Re * .

Figure 4 .
Figure 4. (a), (b) Non-dimensional K+ L versus Re * and Re k .(c) Evaluation of the accuracy of the two proposed models and comparison against zonal model in Voermans et al. (2018b).(d) Evaluation of the three best performing models based on the MLR analysis.
with different D s from D 35 to D 90 .The current study uses k s = 2.5D for uniform spheres in the rough bed.•Shear velocity (u * ) can be estimated by the friction slope (S f ) data as u * = ̅̅̅̅̅̅̅̅̅̅̅̅̅ gH w S f √ , or referred to nomographs and models for hyporheic friction factors, such as those inManes et al. (2012) andVoermans et al. (2018a).•Length scale δ can be calculated by the equations in Figure2as δu * /ν = 0.143Re 1.01 * (δ ∼ 0.143k s ).• Kinematic viscosity (ν) is a temperature relevant parameter.The ν of water can be estimated using the formula:ν = 1.79 × 10 6 / 1 + 0.3368T c + 0.00021T 2 c ) with T c in Celsius.With the parameters above, KL can be obtained from Equation 9 by Re * = k s u * /ν and u * .The effective diffusivity can be determined as D eff = KL δ.It is noted that the developed model from Equation 9 can yield D eff < D m for low Reynolds number flows.This can be adjusted by introducing a limiting parameter where D eff is the maximum of the predicted value and D m (also see SI5 in Supporting Information S1).

outer scales combined with multiple linear regression (MLR) over an extensive data set from the literature and found that D eff /D m has a strong relationship with permeability Reynolds number (Re k =
, δ p and δ p* on Re k .They also identified a critical Re k value (Re k ∼ 1-2) above which turbulence exchange effects dominate over dispersion at the SWI, corresponding to the findings of The financial support by Metropolitan Water Reclamation District of Greater Chicago (MWRDGC) is gratefully acknowledged.The computational support by National Science Foundation via Grant ACI-1548562 for Extreme Science and Engineering Discovery Environment (XSEDE), now migrated to Advanced Cyberinfrastructure Coordination Ecosystem: Services and Support (ACCESS), is appreciated.This work used XSEDE Stampede2 at the Texas Advanced Computing Center (TACC) The University of Texas at Austin through allocation TG-CTS190067 and Frontera at TACC through Pathways allocation CTS22005.