Anomalous transport through free-flow-porous media interface: Pore-scale simulation and predictive modeling

Journal Pre-proof

Anomalous Transport through Free-Flow-Porous Media Interface: Pore-Scale Simulation and Predictive Modeling Jun Song Kim , Peter K. Kang PII: DOI: Reference:

S0309-1708(19)30741-9 https://doi.org/10.1016/j.advwatres.2019.103467 ADWR 103467

To appear in:

Advances in Water Resources

Received date: Revised date: Accepted date:

22 August 2019 3 November 2019 17 November 2019

Please cite this article as: Jun Song Kim , Peter K. Kang , Anomalous Transport through Free-FlowPorous Media Interface: Pore-Scale Simulation and Predictive Modeling, Advances in Water Resources (2019), doi: https://doi.org/10.1016/j.advwatres.2019.103467

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Highlights 

We study pore-scale flow and transport through free-flow-porous media interface.



Porosity exerts dominant control over solute transport through the coupled system.



Pore-scale vortices and turbulence structures control anomalous transport.



Spatial Markov model successfully predicts effective particle transport.

-1-

Anomalous Transport through Free-Flow-Porous Media Interface: Pore-Scale Simulation and Predictive Modeling

Jun Song Kima, b, Peter K. Kanga, b, *

a

Department of Earth and Environmental Sciences, University of Minnesota, Minneapolis, MN 55455, USA

b

Saint Anthony Falls Laboratory, University of Minnesota, Minneapolis, MN 55414, USA

*corresponding author, e-mail: [email protected]

Abstract Pore-scale velocity and turbulence structures near streambeds may control solute transport and dispersion in streams. In this study, pore-scale flow and transport simulations are performed to investigate the effects of pore-scale processes on anomalous transport, which is often manifested by remarkably long residence times of solute particles, in coupled free-flow and porous media systems. By solving the 2D Reynolds-averaged Navier-Stokes equations integrated with a renormalization group k   turbulence model, we resolve complex porescale flows featured by vortices and preferential flows. Then, we incorporate the simulated velocity and turbulence fields into a Lagrangian particle tracking model that considers advection, turbulent diffusion, and molecular diffusion. Simulation results reveal that the interplay between pore-scale vortices and turbulence structures near the free-flow-permeable bed interface controls anomalous transport. High porosity induces strong turbulence penetration and preferential flow paths within permeable beds. The enhanced subsurface -2-

turbulence facilitates the escape of solute particles from recirculation zones via turbulent diffusion, causing steep power-law slopes in breakthrough curves (BTCs). In contrast, low porosity introduces heavy tailing in BTCs from particles that are trapped in the near-interface recirculation zones characterized by low velocities and limited turbulence. We upscale and predict particle transport via a Spatial Markov model (SMM) honoring the interplay between Lagrangian velocity distribution and velocity correlation. The SMM reproduces anomalous transport behaviors obtained from the numerical simulations. These results demonstrate that Lagrangian velocity statistics effectively encode anomalous transport mechanisms in the coupled systems.

Keywords: anomalous transport; free-flow-porous media interface; pore-scale simulation; recirculation zone; turbulence structures; spatial Markov model

1. Introduction

Fluid flow and transport in coupled free-flow and porous media systems is of great importance due to its wide applicability in natural and industrial processes such as nuclear waste storage and fuel cells (You and Liu, 2002; Bourgeat et al., 2009), environmental engineering including flow and pollutant dispersion over aquatic canopies and porous beds in rivers as well as buildings in urban areas (Zhou and Mendoza, 1993; Britter and Hanna, 2003; Nepf, 2012), and biological systems of biofilm formation (Landa-Marbán et al., 2019). Flow and mixing of surface and subsurface waters in a permeable bed, also termed hyporheic flow, is of significant interest in hydrology. Hyporheic flow is an essential hydrodynamic process controlling fluid flow, mass transport, and mixing in rivers and streams (Harvey et al., 1996; -3-

Elliott and Brooks, 1997; Packman et al., 2004). Streamflows are typified by fast turbulent flows in a free-flow layer and slow laminar flows in a porous layer. However, turbulence in surface water can extend to a shallow portion of the porous media, defined as a turbulent transition layer corresponding to a non-Darcy flow regime, as described in Fig. 1 (Nagaoka and Ohgaki, 1990; Pokrajac et al., 2007; Manes et al., 2009). Nevertheless, conventional models often assume that hyporheic flow is laminar, modeling it as Darcy flow (Wörman et al., 2002; Cardenas and Wilson, 2007; Bottacin‐ Busolin and Marion, 2010; Engdahl et al., 2010; Hester et al., 2013; Hester et al., 2019). These approaches often underestimate the hyporheic exchange rate by orders of magnitude, particularly in gravel beds with high permeability (O’Connor et al., 2008; O'Connor et al., 2012). This is mainly because the extended turbulence into the transition layer, commonly ignored in typical models, significantly contributes to mass transfer across the free-flow-porous media interface (Higashino and Stefan, 2011; Roche et al., 2018; Voermans et al., 2018). To accurately model the effects of hyporheic exchange on solute transport in hydrologic systems, it is critical to improve our mechanistic understanding of turbulent flow and transport through the free-flowporous media interface. Velocimetry techniques like particle image velocimetry (PIV) and laser Doppler velocimetry (LDV) are increasingly being employed in exploring flow interactions between a free stream and a permeable bed (Goharzadeh et al., 2005; Liu et al., 2008; Pokrajac and Manes, 2009; Blois et al., 2012; Yang et al., 2015; Kim et al., 2018; Wu and Mirbod, 2018). A common observation from these experiments is that turbulent shear stress is maximized at the free-flow-porous media interface, with the strong interfacial turbulence causing extension of the turbulence structure into permeable beds. From the experimental findings, we can infer that solutes may significantly diffuse into and out of the streambed by interfacial turbulent -4-

mixing of surface and subsurface waters. During such an exchange process, retarded subsurface velocities increase residence times of the solutes, thereby producing anomalous (non-Fickian) transport in streams (Fernald et al., 2001; Cardenas, 2008a; Aubeneau et al., 2014). The key features of the anomalous transport are the unusually early arrival of tracers downstream and remarkably long residence time of tracers. In field tracer tests, breakthrough curves (BTCs) commonly exhibit anomalous late-time tailings attributed to the turbulent hyporheic exchange, which cannot be described with Fickian advection-dispersion equation (ADE) models (Harvey and Bencala, 1993; Haggerty et al., 2002; Jonsson et al., 2003; Aubeneau et al., 2014). Recent studies investigated interfacial solute transport across the interface between the free-flow region and the porous media region at Darcy scale. Aquino et al. (2015), Li et al. (2017) and Roche et al. (2019) implemented particle tracking methods on theoretically and experimentally derived profiles of streamwise velocity and turbulent diffusion by assuming that solute particles are only advected with positive streamwise velocities. Sherman et al. (2019) simulated particle trajectories subject to instantaneous flow fields from direct numerical simulation (DNS), modeling subsurface flows with the drag coefficient representing flow resistance, like a Darcy model. These studies, however, modeled subsurface flows using simple analytical solutions or with a drag force term without directly considering or resolving the interfacial mass transfer phenomenon at the pore scale. These macroscale approaches may underestimate the inertial effects of complex pore-scale flow fields including preferential flows and vortices on particle transport. Several studies based on PIV measurements revealed that vortices (recirculation zones) emerge at pore spaces between coarse grains beneath open channel flows (Blois et al., 2012; Kim et al., 2018). The recirculation zones often act as transient storage zones, where solutes -5-

are trapped and released slowly into channel flows, highlighting a possible source of anomalous transport (Cardenas et al., 2007; Cardenas, 2008b; Crevacore et al., 2016; Zhou et al., 2019). Therefore, explicitly resolving subsurface flow structures at the pore-scale is necessary, to fully understand the non-Fickian nature of transport in hyporheic flows. Some studies involved fully resolved pore-scale flow simulations of coupled turbulent free-flow and porous media flow, where the porous media is idealized with an array of cubic or spherical grains. The flow simulation were performed using Reynolds-averaged NavierStokes (RANS) equations with turbulence models like a standard k   model and a Reynolds stress model (Prinos et al., 2003; Yang et al., 2018), as well as large eddy simulation (Leonardi et al., 2018; Fang et al., 2018; Lian et al., 2019) and DNS (Breugem and Boersma, 2005; Kuwata and Suga, 2017). These studies suggest that recirculating flow structures dominate flow in some portions of porous media, and moreover, geometric characteristics like subsurface porosity strongly correlate to the size of these pore-scale vortices. Yet, these studies only involved flow analysis without studying transport phenomena. Thus, the effects of the pore-scale flow patterns on solute transport through the free-flowporous media interface remain largely unknown. The goal of this study is to evaluate the effects of the complex interplay between pore-scale velocity and turbulence structures near the free-flow-porous media interfaces, on flow and mass transport through coupled free-flow and porous media systems. More specifically, we investigate how the pore-scale flow structures such as vortices and preferential flows influence anomalous transport in a coupled free-flow and porous media system. In this study, we fully resolve pore-scale flow fields directly coupled with overlying turbulent free-flows by solving 2D RANS equations incorporating a turbulence model. This hydrodynamic model is first validated with the experimental data (Prinos et al., 2003), and then applied to -6-

investigate cases with different porosities and a wide range of Reynolds numbers. Based on simulated velocity and turbulence fields, we perform particle tracking simulations involving advection, turbulent

diffusion, and molecular diffusion, to advance fundamental

understanding of flow and transport through free-flow-porous media interfaces. We then perform detailed analysis to elucidate the influence of pore-scale vortices on anomalous transport. Finally, we upscale solute transport using a stochastic random-walk model parametrized with particle trajectory information from the numerical simulations, and then conduct predictions of effective transport using the upscaled model.

2. Methods

2.1 Hydrodynamic model

We analyze mean and turbulent flow characteristics in coupled free-flow and porous media systems using the open-source fluid dynamics library, OpenFOAM which adopts a finite volume method for solving governing equations. We use 2D RANS equations as the governing equations for fluid flow. 2D turbulence models can provide the general picture of the turbulent flows when one spatial direction is considerably constrained by geometry or applied body forces (Boffetta and Ecke, 2012; Xia and Francois, 2017). 2D turbulence model is relevant in this study because the mass exchange between the free-flow zone and the porous media zone is mostly triggered by the vertical flows (Boano et al., 2014; Gomez‐ Velez and Harvey, 2014). Unsteady forms of continuity and momentum equations in Einstein notation in Cartesian coordinates can be written as:

-7-

U i 0 xi

U i U i 1 P  U j   t x j  xi x j

(1a)

  U U j   i    x j xi

  ui u j  Si     x  j 

(1b)

where U i and U j are the time-averaged velocity vectors in the xi - and x j -directions of Cartesian coordinates, respectively. Herein, xi and x j represent x -axis and y -axis properties when i , j = 1, while xi and x j denote y -axis and x -axis properties when i ,

j = 2. P is the fluid pressure;  is the density of the fluid;  is the kinematic viscosity of the fluid; ui and u j are the fluctuating velocity vectors in the xi - and x j -directions, respectively; Si is the source term representing the pressure gradient to drive the flow. In addition, ui u j indicates the Reynolds stress, a product of the Reynolds decomposition, and the term should be modeled with turbulence closure models. The RANS turbulence models employ an eddy viscosity concept derived from Boussinesq approximation which relates Reynolds stresses to mean properties of turbulent flows. A standard k   model is one of the most widely used turbulence models in environmental hydraulics applications because of its good convergence with a minor computational overhead (Mößner and Radespiel, 2015). To close the governing equations, this turbulence model includes two additional transport equations to solve for turbulent quantities of flows such as turbulent kinetic energy, k and its dissipation rate,  . Nevertheless, this conventional method has a limitation in accurately reproducing complex wall-bounded flows -8-

at low Reynolds number, owing to the overestimation of an eddy viscosity (Shih et al., 1995). Hence, we adopt a renormalization group (RNG) k   model (Yakhot et al., 1992) to enhance the accuracy of flow simulation in the free-flow and porous media systems. The RNG k   model better resolves the near-wall turbulent flows, and many studies have reported that the model adequately reproduces flow characteristics near-wall regions such as recirculation zones (Speziale and Thangam et al., 1992; Dargahi, 2006; Cassan and Belaud, 2011; Sinha et al., 2017). The corresponding transport equations for k and  are as follows:

     k     kU j   t xi x j

 t    k 

 k    Gk    Sk   x j 

(2a)

        U j   t xi x j

 t     

    *    C1 Gk  C 2    S   x k  j 

(2b)

where C1 ,  k and   are the empirical constants set to 1.42, 0.7194 and 0.7194, respectively (Yakhot et al., 1992);  and t are the molecular viscosity and the dynamic turbulent viscosity, respectively; Gk is the generation of turbulent kinetic energy arising from the mean velocity gradients. S k and S are the moduli of the mean rate of strain tensors of k and  , respectively. The main difference from the standard k   model is that the RNG k   model treats C *2 as a variable parameter to consider the different scales of motion by changing the production term which can be expressed as:

-9-

C *2  C2 

C (1   / 0 )

(3)

1   3

where C2 , C ,  and  0 are constants set to 1.68, 0.09, 0.012 and 4.38, respectively (Yakhot et al., 1992). The variable parameter,  , is obtained as:



k



2Sij Sij

(4)

where Sij is the strain rate tensor.

2.2 Particle tracking model

We simulate transport and spreading of solute particles using a 2D Lagrangian particle tracking (LPT) model. Particle motion is governed by advection, molecular diffusion, and turbulent diffusion across the free-flow-porous media interface (O'Connor and Harvey, 2008; Voermans et al., 2018). In consequence, the particle trajectories are determined by the stochastic Langevin equation described as follows in Einstein notation.

xi , p (t  t )  xi , p (t )  U i ( xi , p , t )  ui ( xi , p , t )  t  i (t ) 2Dm t

(5)

where xi , p is the particle position; i (t ) is the unit-Gaussian distributed random number; and Dm is the molecular diffusion coefficient. The fluctuating velocity, ui in Eq. (5) - 10 -

characterizes turbulent diffusion whose length scale is usually much larger than that of molecular diffusion so that the former typically governs solute mixing in rivers and streams (Rutherford, 1994). However, the effects of molecular diffusion can be important in porous media, and we consider both turbulent diffusion and molecular diffusion in this study. Since ui is not directly available from the RANS simulation, we apply a discrete random walk (DRW) technique (Gosman and Loannides, 1983) to consider the effects of turbulent mixing on particle trajectories. We model turbulent diffusion of particles as a succession of interactions between particles and eddies which have the finite length and time scales (Dehbi, 2008). The DRW method incorporates the mean turbulent flow properties, and has been shown to successfully reproduce particle dispersion in environmental flows (Shams et al., 2002; Tian and Ahmadi, 2007; Gao et al., 2012; Jayaraju et al., 2015). According to the DRW method, ui can be coupled with k obtained from the flow simulation, where k is the space-dependent variable, as follows:

ui ( xi , p , t )   i (t ) 2k / 3

(6)

where  i (t ) is the unit-Gaussian distributed random number for the velocity perturbation which is temporally varying so that  i (t ) is updated after each eddy lifetime,  e , defined as:

k

 e  C 3/4 ln(r ) 

(7)

- 11 -

where r is the random number between zero and unity. Thus, particles are assumed to be interacting with an eddy which has a constant value of  i during  e . When the time from the origin time of the eddy exceeds  e , another particle-eddy interaction is generated by introducing new values of  i and  e (Rybalko et al. 2012). This process is repeated to consider spatio-temporal evolution of turbulent eddies.

2.3 Upscaled Transport Model: Spatial Markov model

We upscale interfacial solute transport phenomena with an upscaled stochastic transport model, known as Spatial Markov Model (SMM) which honors spatial Markov properties of Lagrangian velocity transitions (Le Borgne et al., 2008). This effective random walk model, which is parameterized with Lagrangian velocity statistics, has been successfully applied to predict anomalous transport across a broad range of flows (Le Borgne et al., 2011; de Anna et al., 2013; Kang et al., 2014; Dentz et al., 2016; Sund et al., 2017). This modeling framework incorporates Lagrangian velocity statistics by honoring the interplay between velocity distribution and velocity correlation. Velocity distribution captures velocity heterogeneity with one-point statistics, and velocity correlation captures spatial memory of velocity transitions with two-point statistics (Dentz and Bolster, 2010; Kang et al., 2015b). The series of particle positions and travel times in streamwise direction are described as:

x( n1)  x( n) 

(8)

t ( n1)  t ( n )   ( n )

(9) - 12 -

where  ( n )  / v( n) is the successive transition time of particles over the fixed spatial jump, , at the finite step which varies from zero to ntotal where ntotal is the channel length, L , (n) . v is the corresponding Lagrangian velocity so that we can characterize the

divided by

velocity distribution using a probability density function (PDF) of  , p( ) . In the SMM framework,  =  ( n ) is correlated to a particle transition time at a previous step,   =  ( n 1) as the transition times are shown to exhibit spatial Markovianity if sampled in space (Le Borgne et al., 2008). Hence,  can be modeled using a transition matrix that characterizes one-step correlation. To assess the transition probabilities from the particle trajectories simulated with the 2D LPT model, we discretize p( ) into N classes,

Ci 1i N ,

where i increases with the increase in  , and further define the conditional

probability density, r ( |  ) . In this study, the transition times are discretized equidistance in logarithmic scale to provide a better discretization for small transition times that control long residence time of particles (Kang et al., 2015a). With these Lagrangian properties, the transition matrix is subsequently defined as (Le Borgne et al., 2008):

 i1

Tij 

 j 1

  i

r ( |  ) p( )d d

j

 j 1



(10)

p( )d 

j

where Tij describes the transition probability from the current transition time class j to the next transition time class i .

- 13 -

2.4 Configuration of numerical experiments

2.4.1 Model validation cases

The RANS simulation incorporating the RNG k   turbulence model yields timeaveraged properties of turbulent flows. For the model validation, we compare our simulation results with experimental data available from Prinos et al. (2003). They investigated the 2D mean and turbulent flow characteristics in a 10 m long and 0.25 m wide open channel with an underlying permeable bed consisting of regular rods with a diameter, d = 0.01 m in the non-staggered pattern, as shown in Fig. 2. Table 1 displays the geometric and hydrodynamic characteristics of the flow. The free-flow depth, h f , is set to either 0.03 m or 0.05 m referred to as Case 250-30 and Case 250-50, respectively. The porous layer depth, hp , is 0.055 m and porosity is 0.8286. With these experimental setups, Prinos et al. (2003) also conducted the numerical simulation using the 2D RANS equations based on the standard k   turbulence model.

2.4.2 Study cases

After validating the hydrodynamic model with the experimental data, we further study two different subsurface-porosity scenarios where porosity is either 0.5 or 0.8, as shown in Fig. 3. The diameter of grains is set to 0.01 m so that this setup can be considered as coarse-gravel beds. Here, h f and hp are fixed to 0.045 m and 0.075 m, respectively. In addition, we vary Reynolds number, Re f = U f h f / , by varying the inlet velocity, U in , where U f denotes - 14 -

the mean velocity of the surface layer over 0 ≤ y /h f ≤ 1. The values of Re f vary approximately from 2,000 to 12,000 which lead to Froude number, Fr = U f / gh f , to vary between 0.11 and 0.59. The investigated Fr values are observable values in natural stream environments. Most of natural streams flow at Fr below 0.3, and the flow at Fr above 0.6 are infrequently observed (Davidian and Cahal, 1963). The geometric and hydraulic conditions as well as the computational mesh resolution for each study case are summarized in Table 2. By studying these cases, we can assess the sensitivity of flow and transport dynamics to porosity and Re f . For the 2D LPT simulation, the total particle number of 5 x 104 is uniformly injected at the inlet of the free-flow layer. In this study, we consider tracer injection at streams and no particles are injected at the porous layer mimicking tracer injection at a stream. The time step is set to values which retain the Courant number, Cr = U i t /x , less than unity, and Dm is fixed to 2  10-9 m2/s representing the typical

molecular diffusivity of liquid phases (Li and Gregory, 1974).

2.4.3 Simulation conditions

The computational meshes are composed of quadrilateral cells. The dense grids satisfying the dimensionless wall distance, y  = yU * / , to be less than unity near the fluid-grain interfaces are generated. Note that U * =

gh f S is the shear velocity. This grid refinement

is necessary to ensure the application of the enhanced wall treatment for accurately simulating the flows in the boundary layers (Salim and Cheah, 2009; Escue and Cui, 2010). The resolutions of generated meshes for the model validation and study cases are listed in Table 1 and Table 2, respectively. We apply the periodic boundary conditions to simulate the - 15 -

mean and turbulent flows for representative computational domains, as illustrated in Fig. 2, because the validation and study cases have periodically repeating geometries in the streamwise direction. At the free surface (zero-shear slip wall), the symmetry boundary condition is applied under the condition of no surface oscillation, which is valid for Fr less than unity (subcritical flow). At the solid-fluid boundaries (at the grain surfaces), the no-slip boundary conditions are used. The convection, diffusion and gradient terms of the momentum equation are discretized using the Gauss linear scheme. This second-order and unbounded scheme is also used for discretizing the turbulence transport equations. The backward scheme of the second order and implicit method handles the discretization of the time derivative term, for which the time step, t is set to the value satisfies Cr to be less than 0.5. The Pressure-Implicit Split-Operator

(PISO) algorithm is employed for coupling pressure and velocity to ensure mass conservation (Prasad et al., 2015). We simulate until U i , P , k and  in Eqs. (1) and (2) converge to the steady-state values corresponding to the fully developed flow condition. These values at the steady state are used for further flow analysis and particle tracking simulation.

3. Results and discussion

3.1 Results of model validation

Fig. 4 shows vertical profiles of the mean velocity magnitude and the turbulent kinetic energy for Case 250-30 and Case 250-50. Both the mean velocity and turbulent kinetic energy are normalized with U * to directly compare our simulation results from the RNG - 16 -

k   turbulence model with the experimental and numerical data from Prinos et al. (2003).

To calculate U * , we derive the hydraulic slope, S , from the streamwise pressure gradient associated with Si in Eq. (1b). As shown in Figs. 4a and b, the presence of a permeable bed reduces the mean velocity in the free-flow layer and produces the slip velocity at the freeflow-porous media interface due to the interfacial momentum exchange. The RNG k   turbulence model accurately reproduces these effects of the permeable bed on the mean velocity of the surface water, thereby showing better agreement with the experimental data, compared to the standard k   turbulence model. Yet, relatively large discrepancies between the simulated velocity and the measured velocity are found adjacent to the interface. These differences may be attributed to the fact that the porous layer of the experimental channel is not long enough to reach the fully developed flow as reported by Prinos et al. (2003). As aforementioned, in the 2D LPT model, the turbulent kinetic energy, k in Eq. (2a) is a critical variable that controls effective particle dispersion via turbulent diffusion. Thus, it is critically important to accurately simulate turbulent kinetic energy profiles because the turbulent diffusion term in the 2D LPT model estimated by the turbulent kinetic energy, as shown in Eq. (6), controls the particle diffusion into the porous media. As shown in Fig. 4, the RNG k   turbulence model reproduces spatially varying behaviors of turbulence, and estimates peak values of turbulent kinetic energy and their locations more accurately than the standard k   turbulence model. Furthermore, the interaction of the turbulent free-flows with the permeable bed results in the enhanced interfacial turbulence propagating through the porous regions.

- 17 -

3.2 Analysis of hydrodynamics

3.2.1 Mean and turbulent flow characteristics

We apply the validated hydrodynamic model to explore the influence of porosity and Re f on the turbulent flow characteristics over and within permeable beds using the study cases shown in Table 2. Fig. 5 presents vertical profiles of the simulated mean velocity magnitude and turbulent kinetic energy for each study case. The mean velocity of the permeable bed is significantly higher in the porosity of 0.8 than in the porosity of 0.5 (Figs. 5a and b). This is because the higher porosity corresponds to more flows through the permeable bed, and accordingly relatively low surface velocity is obtained to honor mass conservation (Table 2). Moreover, Fig. 5b inset shows that porosity determines U p /U f , where U p is the mean velocity of the subsurface layer over -1 ≤ y /hp ≤ 0, because of the strong preferential flow paths observed between grains for the high-porosity cases. The porosity of 0.8 yields U p /U f ≈ 0.1 which is one order of magnitude greater than a value of U p /U f ≈ 0.01 for the porosity of 0.5. On the other hand, this velocity ratio tends to remain relatively constant with respect to U in because the mean velocity structures in the limited water depth are more dependent on porosity than Re f (Goharzadeh et al., 2005; Wu and Mirbod, 2018). The porosity further plays a vital role in turbulence structures. Figs. 5c and d show that the higher porosity induces the greater penetration of turbulent kinetic energy into the permeable bed. In order to quantify the penetration extent of turbulent flows, we calculate the

- 18 -

penetration length scale of turbulent kinetic energy,  k , normalized by pore diameter, d , as follows:

k  

1 (kint / U *2 )d



y 0 m y 0.075 m

f ( y )dy

(11)

*2 where f ( y) indicates the vertical profile of k /U *2 inside the porous regions; and kint /U

is k /U *2 at the free-flow-porous media interface of y = 0 m. Fig. 5d inset shows that the porosity of 0.8 introduces  k ranging 0.420-0.667, which is about three times the values of

 k ranging 0.172-0.187 for the porosity of 0.5. In the high-porosity cases,  k slightly increases with increasing Re f , which is originated from the fluctuation of the small turbulent kinetic energy below y /hp ≈ -0.2 in the turbulent kinetic energy profile (Fig. 5d). However, the turbulent kinetic energy near the interface of y /hp between -0.2 and 0 is less sensitive to Re f . This implies that slight increase in  k may have negligible effects on solute transport

because solute exchange between the free-flow and porous layers mostly takes place in the vicinity of the interface. In summary, the flow analysis results indicate that velocity and turbulence structures are significantly more sensitive to porosity than Re f .

3.2.2 Recirculation zones

The exchange of momentum between free-flow and permeable layers, which is associated with the interfacial slip velocities and turbulence, leads to the occurrence of recirculation - 19 -

zones beneath the free-flow-porous media interface (Breugem and Boersma, 2005; Kuwata and Suga, 2017), as depicted in Fig. 6. This figure also shows the spatial maps of mean flow velocity and turbulent kinetic energy for Case 5010 and Case 8010. The structures of recirculation zones rely strongly on porosity as the high-porosity bed produces a large single vortex, while the low-porosity bed produces a pair of relatively small vortices directly below the interface. The vertical stretch of the first vortex from the interface (squares with red dashed lines in Fig. 6) enlarges proportional to  k and U p /U f , extending to the porous depth of y /hp ≈ -0.1 and y /hp ≈ -0.05 for the porosity of 0.8 and 0.5, respectively. This tendency is consistent with simulation results from Yang et al. (2018) who reported that the vortex size increases monotonously with an increase in U p /U f . The length scale of the vortices may significantly affect interfacial mass-transfer dynamics because the greater extent of the near-interface recirculating flows can cause solutes to penetrate deeper into the porous regions. As solute exchange between the free-flow and the porous media mainly occurs near the interface due to the steep decay of interfacial turbulence towards the channel bed (Figs. 5c and d), the mean and turbulent flow characteristics of the near-interface recirculation zones can be considered as one of the key factors governing mass transfer between the surface and subsurface layers. As shown in Fig. 6, the higher mean velocity and the stronger turbulence are observed when the porosity increases. For the vortex regions located within the first grain layer (-0.133 ≤ y /hp ≤ 0), the spatial-averaged values of U /U * and k /U *2 for the porosity 0.8 case are U /U * = 0.58 and k /U *2 = 0.91, respectively, which are about 1.5 times the values of U /U * = 0.34 and k /U *2 = 0.65 for the porosity 0.5 case. Such difference in the properties of the recirculation zones would be originated from the interfacial momentum - 20 -

exchange intrinsic to porosity. The slow recirculating flows with the diminished turbulence for the low porosity cases may result in the strong trapping effects of the pore-scale vortices, which can invoke the heavy tailing in BTCs by not only increasing residence times of solutes but also restricting the solutes to escape from these recirculation zones. The interaction between the recirculation zones and solute transport will be further discussed in the section 3.3.3.

3.3 Analysis of transport dynamics

3.3.1 Tracer breakthrough curves

We now analyze transport behaviors using BTCs (residence time distributions). Fig. 7 shows the BTCs at the downstream station located at L  3 m obtained from 2D LPT simulations. Note that particles can breakthrough through both the free-flow and porous layers. As discussed in the section 3.2.1, the mean and turbulent flow structures governing the particle transport patterns are not sensitive to the change in Re f . This is further demonstrated by showing that the BTCs normalized by U in can be collapsed for the different values of Re f (Figs. 7a and b insets). The dotted lines in Figs. 7a and b show the BTCs of particles

that did not travel through porous media (below y /hp = 0) before reaching the downstream stations. The residence times of these particles are much shorter than those of the particles that interact with the porous layer (solid lines in Figs. 7a and b). This consequently reveals that solute exchange at the free-flow-porous media interface is the crucial dynamics

- 21 -

responsible for anomalous transport patterns distinguished by the tailing behaviors in the BTCs. We investigate the effects of porosity on mass transfer of particles to the permeable bed by measuring the percent mass of particles accumulated in the porous regions, %M sub , which is defined as the number of particles in the porous layer divided by the number of injected particles, for each time step during the simulation periods. Fig. 8 shows that %M sub initially increases as time elapses, and then converges to specific values. For the porosity of 0.8, particles accumulate up to %M sub ≈ 16% which is almost twice the value of the porosity 0.5 cases. This infers that the enhanced turbulence inside the high-porosity bed stimulates the particles to diffuse more easily into the porous regions. As expected, whereas, we cannot observe the significant difference in %M sub with the variation of Re f . Interestingly, despite the smaller number of the particles captured in the low-porosity bed owing to the limited extent of turbulence penetration into the permeable bed, the porosity of 0.5 introduces stronger tails in BTCs compared to those of the porosity of 0.8, as shown in Fig. 9a. This seemingly contradictory behavior in the BTCs can be explained by understanding the interplay between vortices and turbulence structures within the permeable bed. Once particles diffuse into the subsurface regions, the particles are significantly more stagnant in the low-porosity bed than in the high-porosity bed due to the stronger trapping effects by the recirculation zones. The strong trapping effects originate from the low velocities and limited turbulence, which lead to the strong tailing in the BTCs. On the other hand, the high-porosity bed allows particles to easily escape from the recirculation zones with the enhanced subsurface turbulence and transport through the preferential flow paths as shown in Fig. 6b. This key finding from the 2D LPT simulation is consistent with the result - 22 -

of field tracer tests from Aubeneau et al. (2014) who reported that pea-gravel beds (low porosity), which produces low subsurface velocities, resulted in longer retention times of solute tracers than coarse-gravel beds (high porosity).

3.3.2 Lagrangian velocity statistics

We characterize the Lagrangian velocity statistics using the Lagrangian velocity distribution and correlation. The motion of the solute particles is discretized equidistantly with the fixed spatial step of

in the main flow direction. Accordingly, the velocity

distribution can be represented by the distribution of the transition times,  ( n ) = previously defined. The transition times are sampled at every

/ v( n ) , as

= 0.05 m from the particle

trajectories. Fig. 9b shows the PDF of the transition times which represents the transition time distribution. From this figure, the probability of having the long transition times (low velocities) increases with the decreasing porosity, leading to the heavy-tailed transition time distribution. Furthermore, Fig. 9 demonstrates that the transition time PDF has the power-law slope identical to that of the corresponding BTCs. Therefore, this implies that the effects of porosity on the anomalous tailing behavior in the BTCs are encoded adequately in the transition time distributions. However, the transition time distribution alone does not include the property of the spatial velocity correlation structure which is another key factor that determines the anomalous transport behavior (de Anna et al., 2013; Kang et al., 2015b). Thus, we first evaluate the velocity autocorrelation functions for the given lag,  s (s) , and then estimate the correlation length scale,

c

by integrating  s (s) over the distance (Dentz et al., 2016; Kang et al., - 23 -

2017). Fig. 9b inset shows that

c

is larger for larger porosity cases. This tendency would be

the result of the fast-preferential flows developed in the high-porosity bed (Fig. 6b), which contributes to the increase in the velocity correlation. The existence of the finite correlation length along the particle trajectories implies that the Lagrangian velocities (transition times) can be effectively describes as a spatial Markov process (Le Borgne et al., 2008; Dentz et al., 2016; Kang et al., 2017). To consider the velocity correlation effects, the series of the transition times sampled in space is characterized using the transition matrix in Eq. (10). We construct the 30  30 transition matrix to adequately capture the entire ranges of correlation effects (Le Borgne et al., 2011), as shown in Fig. 10. Herein, each transition time class is spaced equidistantly in the log scale, and we choose length scale,

c

= 0.05 m which is smaller than the characteristic correlation

, to properly resolve the velocity correlation structure. The constructed

transition matrix is used to evaluate the one-step correlation between  ( n ) and  ( n1) . Fig. 10 shows that the correlation is stronger for the high velocities (short transition times) than for the low velocities (long transition times). The porosity exerts a significant impact on the transition matrices. The transition matrix for the porosity of 0.8 shows the emergence of the strong velocity correlation in the transition time classes between 9 and 10 (dashed box in Fig. 10b) corresponding to the preferential flow paths in porous media (Fig. 6b), where advection is dominant over molecular diffusion and turbulent diffusion. Particles can transport with high subsurface velocities through these preferential flow paths, as illustrated with the snapshot of velocity field overlain by particle trajectories in Fig. 10b. The dashed box in Fig 10a is the result of the recirculation zones attributed to the porosity of 0.5, which will be discussed in detail in the next section 3.3.3. - 24 -

3.3.3 Effects of recirculation zones on transport

In order to assess the contributions of pore-scale vortices on effective particle transport, we quantify particle residence times in these local recirculation zones by measuring travel times of particles within the flow separation regions enclosed by dividing streamlines, as indicated in Fig. 6. The residence time distribution of the particles captured in the recirculation zones defined above is shown in Fig. 11, where the tail in the residence time PDF is heavier in the porosity of 0.5 than in the porosity of 0.8. This is because, for the low porosity, the particles stay longer in the recirculation zones owing to the low velocity in the recirculation zones and the limited turbulence. From this result, one can notice that power-law slopes of the residence time distributions for the recirculation zones resemble those of the BTCs (Fig. 9a). This implies that the evolution of the late-time tailing in the BTCs is controlled by the mean and turbulent flow properties of the recirculation zones. Also, we consistently observe that the BTC tails depend strongly on porosity rather than Re f as the mean residence times of the particles in the recirculation zones remain relatively constant across different flow conditions, as shown in Fig. 11 inset. The transition matrix (Fig. 10) also reveals that the importance of the pore-scale vortices on particle transport. For the case with porosity 0.5, the particles belonging to the current transition time classes j = 6-29 have high probability of traveling with the transition times of the classes i ≈ 6 in the next

(dashed box in Fig. 10a). These next transition time

classes correspond to the flow velocities immediately above the free-flow-porous media interface. Once the particles enter the low-porosity bed, the limited vertical extent of the near- 25 -

interface vortices and turbulence penetration (Fig. 6a) makes particles to stay near the interface, as described with the snapshot of particle trajectories in Fig. 10a. Hence, the particles in the porous region either stay in recirculation zones or diffuse back to the freeflow regions through interfacial turbulent diffusion. Such transport mechanisms make the particles to mostly jump to the next

with the interfacial flow velocities represented by the

transition time classes i ≈ 6. The consequence of this intrinsic particle motion inherent to the low porosity is well illustrated in Fig. 9. This figure shows that, for the porosity of 0.5, the BTC measured at both the free-flow and porous layers (solid line) is almost identical with that measured at the free-flow layer only (dotted line), which proves that passage routes through the downstream station are mostly limited to the free-flow layer. For the porosity of 0.8, in contrast, the particles exit the study reach through both free-flow and porous layers as they can jump actively into the fast-preferential flow paths due to the extended stretch of the near-interface vortices and turbulence penetration (Fig. 6b).

3.4 Predictions with Spatial Markov Model

Fig. 12 shows BTCs estimated at three different downstream locations for Case 5010 and Case 8010. The comparison between the SMM predictions and the 2D LPT results shows that the SMM successfully reproduces the first passage time distributions at all downstream stations, and accurately captures the anomalous tailing behaviors in the BTCs, caused by the permeable bed-induced velocity and turbulence structures. Such good agreement between the SMM results and the 2D LPT model simulations validates that the Lagrangian velocity statistics such as the velocity distribution and correlation effectively encode the particle transport and dispersion dynamics in the coupled free-flow and porous media systems. Note - 26 -

that we applied the SMM to all study cases, but only presented two representative cases because no significant difference in the predicted BTCs is found as a function of the flow condition. In contrast, an uncorrelated continuous time random walk (CTRW) model, which does not consider the spatial velocity correlation between successive jumps, fails to predict the tailing behavior, as shown in Fig. 12(b) inset. This demonstrates that the spatial velocity correlation plays a critical role in determining the late-time tailing in BTCs. In addition, to further validate the predictability of the SMM, we parameterize the SMM with the transition times obtained from the first 200 s of simulation period. The SMM parametrized with the limited time window still predicts the BTCs successfully, as shown in Fig. 12(b) inset. Recent studies also showed effective ways of parameterizing transition matrices of transition times from experimental data without running transport simulations (Kang et al., 2015b; Sherman et al., 2017). However, the focus of this study is not in model parameter estimation but in demonstrating the applicability of the SMM in free-flow porous media interfaces. The tail truncation time of BTCs is one of the important features to evaluate the effects of hyporheic exchange on anomalous transport in streams (Drummond et al., 2012; Runkel 2015). Nevertheless, the truncation time cannot typically be observed in the field-based tracer experiments due to limited experimental periods, thereby being merely approximated by extrapolating the observation data (Aubeneau et al., 2014). To complement this limitation, we use results from SMM to predict the BTC truncation times. According to the BTCs predicted at x /L = 1 using the SMM (Fig. 12), the porosity of 0.5 introduces the truncation time of t ≈ 480 s which is about twice the value of t ≈ 240 s in the porosity of 0.8. Therefore, the truncation time also highlights the importance of porosity in determining the anomalous transport in the coupled free-flow-porous media system. - 27 -

4. Conclusions

This study highlights the importance of interfacial turbulence structures and recirculation flows in controlling solute transport through coupled free-flow-porous media systems by performing the pore-scale simulations. The mean and the turbulent flow profiles simulated via 2D RANS equations involving the RNG k   turbulence model are consistent with those of experimental data. The velocity and the turbulence structures within permeable beds are very sensitive to changes in porosity because the depth of turbulence penetration into the porous layer and the subsurface flow velocity increase with increasing porosity, and also extends the near-interface recirculation zones. The 2D LPT simulations revealed that particle transport is significantly influenced by the complex interplay between the pore-scale vortices and turbulence structures in permeable beds. The extended turbulence allows more active diffusion of particles from the recirculation zones into the free-flow regions or the preferential flow paths, where the high velocities cause rapid breakthrough of solute particles. Consequently, the high porosities cause steeper powerlaw slopes than low porosities. For the low porosity cases, the prominent tailings in BTCs are attributed to particles trapped in the near-interface recirculation zones characterized by low velocity and low turbulent kinetic energy compared to the higher porosity cases. Furthermore, the power-law slopes of the residence time distributions for particles captured in recirculation zones resemble the scaling of BTCs. These results reveal that pore-scale vortices exert dominant control over anomalous transport through free-flow-porous media interfaces. We applied the SMM for predictive modeling and successfully predicted solute transport. We parametrized the upscaled model using the Markovian transition matrix involving the - 28 -

Lagrangian velocity (transition time) distribution and correlation information retrieved from the LPT-simulated particle trajectories. The effects of the pore-scale vortices on transport are well captured in the transition matrix. The transition matrix demonstrates that the limited extent of the near-interface vortices and turbulence penetration inside the low-porosity beds cause retention of particles in the vicinity of the interface, with exit restricted to the free-flow regions. In contrast, for high-porosity beds, the transition matrix shows that particles can also readily enter the preferential flow zones via strong turbulent diffusion. The SMM accurately reproduced the BTCs, indicating that the effects of the pore-scale vortices on transport are effectively encoded in Lagrangian velocity statistics. This study elucidates the mechanisms of anomalous transport in the coupled system and suggests an effective way to incorporate the mechanisms into an upscaled model. However, to further generalize our findings, particle transport in more realistic pore geometries should be investigated. The mean and turbulent flow properties obtained in this study with the RANS equations may be sensitive to the choice of a turbulence model. LES or DNS methods can directly resolve eddies without a turbulence model and can more explicitly demonstrate the effects of recirculation zones on solute transport. However, we expect the main findings of this study will be independent of the turbulence models because the average velocity and turbulence fields should exert first-order control on effective solute transport. Finally, while this study is based on numerical analysis, complementary experiments using PIV and laser-induced fluorescence (LIF) will allow us to further advance our understanding of flow and transport mechanisms in the coupled free-flow and porous media systems. For example, Ling et al. (2018) combined pore-scale numerical simulations with microfluidics experiments to study saturated flow through channel-porous media coupled systems and showed promising results by comparing with upscaled models. - 29 -

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The authors gratefully acknowledge support from the Institute of Engineering Research and Institute of Construction and Environmental Engineering in Seoul National University, Seoul, South Korea. PKK also acknowledges the College of Science & Engineering at the University of Minnesota and the George and Orpha Gibson Endowment for its generous support of Hydrogeology, and a grant from Korea Environment Industry and Technology Institute (KEITI) through Subsurface Environmental Management (SEM) Project, funded by the Korea Ministry of Environment (MOE) (2018002440003). We thank the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for computational resources and support.

References

Aquino, T., Aubeneau, A., Bolster, D., 2015. Peak and tail scaling of breakthrough curves in hydrologic tracer tests. Advances in water resources, 78, 1-8. Aubeneau, A. F., Hanrahan, B., Bolster, D., Tank, J. L., 2014. Substrate size and heterogeneity control anomalous transport in small streams. Geophysical Research Letters, 41(23), 8335-8341.

- 30 -

Blois, G., Smith, G. S., Best, J. L., Hardy, R. J., Lead, J. R., 2012. Quantifying the dynamics of flow within a permeable bed using time-resolved endoscopic particle imaging velocimetry (EPIV). Experiments in Fluids, 53(1), 51-76. Boano, F., Harvey, J.W., Marion, A., Packman, A.I., Revelli, R., Ridolfi, L. and Wörman, A., 2014.

Hyporheic

flow

and

transport

processes:

Mechanisms,

models,

and

biogeochemical implications. Reviews of Geophysics, 52(4), 603-679. Boffetta, G. and Ecke, R.E., 2012. Two-dimensional turbulence. Annual Review of Fluid Mechanics, 44, 427-451. Bottacin‐ Busolin, A., Marion, A., 2010. Combined role of advective pumping and mechanical dispersion on time scales of bed form–induced hyporheic exchange. Water Resources Research, 46(8). Bourgeat, A., Jurak, M., Smaï, F., 2009. Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository. Computational Geosciences, 13(1), 29. Breugem, W. P., Boersma, B. J., 2005. Direct numerical simulations of turbulent flow over a permeable wall using a direct and a continuum approach. Physics of Fluids, 17(2), 025103. Britter, R. E., Hanna, S. R., 2003. Flow and dispersion in urban areas. Annual Review of Fluid Mechanics, 35(1), 469-496. Cardenas, M. B., Slottke, D. T., Ketcham, R. A., Sharp Jr, J. M., 2007. Navier‐ Stokes flow and transport simulations using real fractures shows heavy tailing due to eddies. Geophysical Research Letters, 34(14). Cardenas, M. B., 2008a. Surface water‐ groundwater interface geomorphology leads to scaling of residence times. Geophysical Research Letters, 35(8). - 31 -

Cardenas, M.B., 2008b. Three‐ dimensional vortices in single pores and their effects on transport. Geophysical Research Letters, 35(18). Cardenas, M. B., Wilson, J. L., 2007. Dunes, turbulent eddies, and interfacial exchange with permeable sediments. Water Resources Research, 43(8). Cassan, L., Belaud, G., 2011. Experimental and numerical investigation of flow under sluice gates. Journal of Hydraulic Engineering, 138(4), 367-373. Crevacore, E., Tosco, T., Sethi, R., Boccardo, G., Marchisio, D. L., 2016. Recirculation zones induce non-Fickian transport in three-dimensional periodic porous media. Physical Review E, 94(5), 053118. Dargahi, B., 2006. Experimental study and 3D numerical simulations for a free-overflow spillway. Journal of Hydraulic Engineering, 132(9), 899-907. Davidian, J., Cahal, D. I., 1963. Distribution of shear in rectangular channels. US Geological Service, Professional Paper, (475-C). de Anna, P., Le Borgne, T., Dentz, M., Tartakovsky, A. M., Bolster, D., Davy, P., 2013. Flow intermittency, dispersion, and correlated continuous time random walks in porous media. Physical Review Letters, 110(18), 184502. Dehbi, A., 2008. Turbulent particle dispersion in arbitrary wall-bounded geometries: A coupled CFD-Langevin-equation based approach. International Journal of Multiphase Flow, 34(9), 819-828. Dentz, M., Bolster, D., 2010. Distribution-versus correlation-induced anomalous transport in quenched random velocity fields. Physical Review Letters, 105(24), 244301. Dentz, M., Kang, P. K., Comolli, A., Le Borgne, T., Lester, D. R., 2016. Continuous time random walks for the evolution of Lagrangian velocities. Physical Review Fluids, 1(7), 074004. - 32 -

Drummond, J. D., Covino, T. P., Aubeneau, A. F., Leong, D., Patil, S., Schumer, R., Packman, A. I., 2012. Effects of solute breakthrough curve tail truncation on residence time estimates: A synthesis of solute tracer injection studies. Journal of Geophysical Research: Biogeosciences, 117(G3). Elliott, A. H., Brooks, N. H., 1997. Transfer of nonsorbing solutes to a streambed with bed forms: Theory. Water Resources Research, 33(1), 123-136. Engdahl, N.B., Vogler, E.T., Weissmann, G.S., 2010. Evaluation of aquifer heterogeneity effects on river flow loss using a transition probability framework. Water Resources Research, 46(1). Escue, A., Cui, J., 2010. Comparison of turbulence models in simulating swirling pipe flows. Applied Mathematical Modelling, 34(10), 2840-2849. Fang, H., Han, X., He, G., Dey, S., 2018. Influence of permeable beds on hydraulically macro-rough flow. Journal of Fluid Mechanics, 847, 552-590. Fernald, A. G., Wigington Jr, P. J., Landers, D. H., 2001. Transient storage and hyporheic flow along the Willamette River, Oregon: Field measurements and model estimates. Water Resources Research, 37(6), 1681-1694. Gao, N., Niu, J., He, Q., Zhu, T., Wu, J., 2012. Using RANS turbulence models and Lagrangian approach to predict particle deposition in turbulent channel flows. Building and Environment, 48, 206-214. Goharzadeh, A., Khalili, A., Jørgensen, B. B., 2005. Transition layer thickness at a fluidporous interface. Physics of Fluids, 17(5), 057102. Gomez‐ Velez, J.D. and Harvey, J.W., 2014. A hydrogeomorphic river network model predicts where and why hyporheic exchange is important in large basins. Geophysical Research Letters, 41(18), 6403-6412. - 33 -

Gosman, A. D., Loannides, E., 1983. Aspects of computer simulation of liquid-fueled combustors. Journal of Energy, 7(6), 482-490. Haggerty, R., Wondzell, S. M., Johnson, M. A., 2002. Power‐ law residence time distribution in the hyporheic zone of a 2nd‐ order mountain stream. Geophysical Research Letters, 29(13), 18-1. Harvey, J. W., Bencala, K. E., 1993. The effect of streambed topography on surface‐ subsurface water exchange in mountain catchments. Water Resources Research, 29(1), 89-98. Harvey, J. W., Wagner, B. J., Bencala, K. E., 1996. Evaluating the reliability of the stream tracer approach to characterize stream‐ subsurface water exchange. Water Resources Research, 32(8), 2441-2451. Hester, E. T., Eastes, L. A., Widdowson, M. A., 2019. Effect of Surface Water Stage Fluctuation on Mixing‐ Dependent Hyporheic Denitrification in Riverbed Dunes. Water Resources Research. Hester, E. T., Young, K. I., Widdowson, M. A., 2013. Mixing of surface and groundwater induced by riverbed dunes: Implications for hyporheic zone definitions and pollutant reactions. Water Resources Research, 49(9), 5221-5237. Higashino, M., Stefan, H. G., 2011. Non-linear effects on solute transfer between flowing water and a sediment bed. Water Research, 45(18), 6074-6086. Jayaraju, S. T., Sathiah, P., Roelofs, F., Dehbi, A., 2015. RANS modeling for particle transport and deposition in turbulent duct flows: Near wall model uncertainties. Nuclear Engineering and Design, 289, 60-72.

- 34 -

Jonsson, K., Johansson, H., Wörman, A., 2003. Hyporheic exchange of reactive and conservative solutes in streams—Tracer methodology and model interpretation. Journal of Hydrology, 278(1-4), 153-171. Kang, P. K., de Anna, P., Nunes, J. P., Bijeljic, B., Blunt, M. J., Juanes, R., 2014. Pore‐ scale intermittent velocity structure underpinning anomalous transport through 3‐ D porous media. Geophysical Research Letters, 41(17), 6184-6190. Kang, P. K., Dentz, M., Le Borgne, T., Juanes, R., 2015a. Anomalous transport on regular fracture networks: Impact of conductivity heterogeneity and mixing at fracture intersections. Physical Review E, 92(2), 022148. Kang, P. K., Dentz, M., Le Borgne, T., Lee, S., Juanes, R., 2017. Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes. Advances in Water Resources, 106, 80-94. Kang, P. K., Le Borgne, T., Dentz, M., Bour, O., Juanes, R., 2015b. Impact of velocity correlation and distribution on transport in fractured media: Field evidence and theoretical model. Water Resources Research, 51(2), 940-959. Kaufman, M. H., Cardenas, M. B., Buttles, J., Kessler, A. J., Cook, P. L., 2017. Hyporheic hot moments: Dissolved oxygen dynamics in the hyporheic zone in response to surface flow perturbations. Water Resources Research, 53(8), 6642-6662. Kim, T., Blois, G., Best, J. L., Christensen, K. T., 2018. Experimental study of turbulent flow over and within cubically packed walls of spheres: Effects of topography, permeability and wall thickness. International Journal of Heat and Fluid Flow, 73, 16-29. Kuwata, Y., Suga, K., 2017. Direct numerical simulation of turbulence over anisotropic porous media. Journal of Fluid Mechanics, 831, 41-71.

- 35 -

Landa-Marbán, D., Liu, N., Pop, I.S., Kumar, K., Pettersson, P., Bødtker, G., Skauge, T., Radu, F.A., 2019. A pore-scale model for permeable biofilm: Numerical simulations and laboratory experiments. Transport in Porous Media, 127(3), 643-660. Le Borgne, T., Bolster, D., Dentz, M., de Anna, P., Tartakovsky, A., 2011. Effective pore‐ scale dispersion upscaling with a correlated continuous time random walk approach. Water Resources Research, 47(12). Le Borgne, T., Dentz, M., Carrera, J., 2008. Lagrangian statistical model for transport in highly heterogeneous velocity fields. Physical Review Letters, 101(9), 090601. Leonardi, A., Pokrajac, D., Roman, F., Zanello, F., Armenio, V., 2018. Surface and subsurface contributions to the build-up of forces on bed particles. Journal of Fluid Mechanics, 851, 558-572. Levy, M., Berkowitz, B., 2003. Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. Journal of contaminant hydrology, 64(3-4), 203-226. Lian, Y. P., Dallmann, J., Sonin, B., Roche, K. R., Liu, W. K., Packman, A. I., Wagner, G. J., 2019. Large eddy simulation of turbulent flow over and through a rough permeable bed. Computers & Fluids, 180, 128-138. Ling, B., Oostrom, M., Tartakovsky, A.M., Battiato, I., 2018. Hydrodynamic dispersion in thin channels with micro-structured porous walls. Physics of Fluids, 30(7), p.076601. Liu, D., Diplas, P., Fairbanks, J. D., Hodges, C. C., 2008. An experimental study of flow through rigid vegetation. Journal of Geophysical Research: Earth Surface, 113(F4). Lowell, J. L., Gordon, N., Engstrom, D., Stanford, J. A., Holben, W. E., Gannon, J. E., 2009. Habitat heterogeneity and associated microbial community structure in a small-scale floodplain hyporheic flow path. Microbial Ecology, 58(3), 611-620.

- 36 -

Manes, C., Pokrajac, D., McEwan, I., Nikora, V., 2009. Turbulence structure of open channel flows over permeable and impermeable beds: A comparative study. Physics of Fluids, 21(12), 125109. Marzadri, A., Tonina, D., Bellin, A., 2012. Morphodynamic controls on redox conditions and on nitrogen dynamics within the hyporheic zone: Application to gravel bed rivers with alternate‐ bar morphology. Journal of Geophysical Research: Biogeosciences, 117(G3). Nagaoka, H., Ohgaki, S., 1990. Mass transfer mechanism in a porous riverbed. Water Research, 24(4), 417-425. Nepf, H. M., 2012. Flow and transport in regions with aquatic vegetation. Annual Review of Fluid Mechanics, 44, 123-142. O'Connor, B. L., Harvey, J. W., 2008. Scaling hyporheic exchange and its influence on biogeochemical reactions in aquatic ecosystems. Water Resources Research, 44(12). O'Connor, B. L., Harvey, J. W., McPhillips, L. E., 2012. Thresholds of flow‐ induced bed disturbances and their effects on stream metabolism in an agricultural river. Water Resources Research, 48(8). Packman, A. I., Salehin, M., Zaramella, M., 2004. Hyporheic exchange with gravel beds: basic hydrodynamic interactions and bedform-induced advective flows. Journal of Hydraulic Engineering, 130(7), 647-656. Pokrajac, D., Manes, C., 2009. Velocity measurements of a free-surface turbulent flow penetrating a porous medium composed of uniform-size spheres. Transport in Porous Media, 78(3), 367. Prasad, B., Hino, T., Suzuki, K., 2015. Numerical simulation of free surface flows around shallowly submerged hydrofoil by OpenFOAM. Ocean Engineering, 102, 87-94.

- 37 -

Prinos, P., Sofialidis, D., Keramaris, E., 2003. Turbulent flow over and within a porous bed. Journal of Hydraulic Engineering, 129(9), 720-733. Roche, K. R., Blois, G., Best, J. L., Christensen, K. T., Aubeneau, A. F., Packman, A. I., 2018. Turbulence links momentum and solute exchange in coarse‐ grained streambeds. Water Resources Research, 54(5), 3225-3242. Roche, K. R., Li, A., Bolster, D., Wagner, G. J., Packman, A. I., 2019. Effects of turbulent hyporheic mixing on reach‐ scale transport. Water Resources Research. Runkel, R. L., 2015. On the use of rhodamine WT for the characterization of stream hydrodynamics and transient storage. Water Resources Research, 51(8), 6125-6142. Rybalko, M., Loth, E., Lankford, D., 2012. A Lagrangian particle random walk model for hybrid RANS/LES turbulent flows. Powder Technology, 221, 105-113. Salim, S. M., Cheah, S., 2009. Wall Y strategy for dealing with wall-bounded turbulent flows. In Proceedings of the international multiconference of engineers and computer scientists (Vol. 2, pp. 2165-2170). Hong Kong. Shams, M., Ahmadi, G., Smith, D. H., 2002. Computational modeling of flow and sediment transport and deposition in meandering rivers. Advances in Water Resources, 25(6), 689699. Sherman, T., Fakhari, A., Miller, S., Singha, K. and Bolster, D., 2017. Parameterizing the spatial Markov model from breakthrough curve data alone. Water Resources Research, 53(12), 10888-10898. Sherman, T., Roche, K. R., Richter, D. H., Packman, A. I., Bolster, D., 2019. A Dual Domain stochastic lagrangian model for predicting transport in open channels with hyporheic exchange. Advances in Water Resources, 125, 57-67.

- 38 -

Shih, T. H., Liou, W. W., Shabbir, A., Yang, Z., Zhu, J., 1995. A new k-ϵ eddy viscosity model for high reynolds number turbulent flows. Computers & Fluids, 24(3), 227-238. Sinha, S., Hardy, R. J., Blois, G., Best, J. L., Sambrook Smith, G. H., 2017. A numerical investigation into the importance of bed permeability on determining flow structures over river dunes. Water Resources Research, 53(4), 3067-3086. Speziale, C. G., Thangam, S., 1992. Analysis of an RNG based turbulence model for separated flows. International Journal of Engineering Science, 30(10), 1379-IN4. Storey, R. G., Fulthorpe, R. R., Williams, D. D., 1999. Perspectives and predictions on the microbial ecology of the hyporheic zone. Freshwater Biology, 41(1), 119-130. Sund, N. L., Porta, G. M., Bolster, D., 2017. Upscaling of dilution and mixing using a trajectory based spatial Markov random walk model in a periodic flow domain. Advances in Water Resources, 103, 76-85. Tian, L., Ahmadi, G., 2007. Particle deposition in turbulent duct flows—comparisons of different model predictions. Journal of Aerosol Science, 38(4), 377-397. Voermans, J. J., Ghisalberti, M., Ivey, G. N., 2018. A model for mass transport across the sediment‐ water interface. Water Resources Research, 54(4), 2799-2812. Wörman, A., Packman, A. I., Johansson, H., Jonsson, K., 2002. Effect of flow‐ induced exchange in hyporheic zones on longitudinal transport of solutes in streams and rivers. Water Resources Research, 38(1), 2-1. Wu, Z., Mirbod, P., 2018. Experimental analysis of the flow near the boundary of random porous media. Physics of Fluids, 30(4), 047103. Xia, H. and Francois, N., 2017. Two-dimensional turbulence in three-dimensional flows. Physics of Fluids, 29(11), 111107.

- 39 -

Yakhot, V., Orszag, S. A., Thangam, S., Gatski, T. B., Speziale, C. G., 1992. Development of turbulence models for shear flows by a double expansion technique. Physics of Fluids A: Fluid Dynamics, 4(7), 1510-1520. Yang, J. Q., Kerger, F., Nepf, H. M., 2015. Estimation of the bed shear stress in vegetated and bare channels with smooth beds. Water Resources Research, 51(5), 3647-3663. Yang, G., Weigand, B., Terzis, A., Weishaupt, K., Helmig, R., 2018. Numerical simulation of turbulent flow and heat transfer in a three-dimensional channel coupled with flow through porous structures. Transport in Porous Media, 122(1), 145-167. You, L., Liu, H., 2002. A two-phase flow and transport model for the cathode of PEM fuel cells. International Journal of Heat and Mass Transfer, 45(11), 2277-2287. Yuan-Hui, L., Gregory, S., 1974. Diffusion of ions in sea water and in deep-sea sediments. Geochimica et Cosmochimica Acta, 38(5), 703-714. Zhou, D., Mendoza, C., 1993. Flow through porous bed of turbulent stream. Journal of Engineering Mechanics, 119(2), 365-383. Zhou, J. Q., Wang, L., Chen, Y. F., Cardenas, M. B., 2019. Mass Transfer Between Recirculation and Main Flow Zones: Is Physically Based Parameterization Possible?. Water Resources Research, 55(1), 345-362.

- 40 -

Table 1. The experimental and numerical setups for model validation cases from Prinos et al. (2003). Parameter

Description

Case 250-30

Case 250-50



Porosity

0.8286

0.8286

h f (m)

Surface flow depth

0.03

0.05

h f /H

Relative depth

0.353

0.476

U * (m/s)

Shear velocity

0.0243

0.0313

Re f

Reynolds number

7,581

14,370

Fr

Froude number

0.858

0.586

Ng

Number of cells

119,641

125,245

- 41 -

Table 2. The key hydraulic and geometric characteristic for the study cases. Parameter

Description

Case 5002

Case 5005

Case 5010

Case 8002

Case 8005

Case 8010



Porosity

0.50

0.50

0.50

0.80

0.80

0.80

U in (m/s)

Inlet velocity

0.02

0.05

0.10

0.02

0.05

0.10

h f (m)

Surface flow depth

0.045

0.045

0.045

0.045

0.045

0.045

h f /H

Relative depth

0.375

0.375

0.375

0.375

0.375

0.375

U f (m/s)

Mean surface velocity

0.053

0.131

0.261

0.047

0.116

0.230

U * (m/s)

Shear velocity

0.0058

0.0130

0.0241

0.0056

0.0126

0.0238

Re f

Reynolds number

2,388

5,915

11,750

2,202

5,230

10,367

Fr

Froude number

0.12

0.30

0.59

0.11

0.26

0.52

Ng

Number of cells

80,912

80,912

80,912

123,475

123,475

123,475

- 42 -

U

Surface flow

Free-flow layer

Flow direction Mass exchange

Interface Turbulent transition layer

Subsurface flow Porous layer

Fig. 1. A schematic showing flow structure and mass exchange in a coupled free-flow and porous media system.

- 43 -

Fig. 2. The geometric configuration of model validation cases from Prinos et al. (2003) and a computational domain with boundary conditions.

- 44 -

Fig. 3. The geometric setups for the study cases with two different porosity values.

- 45 -

(b)

(d)

(c)

y/hp

y/hf

(a)

k  k 

U/U*

U/U*

k/U*2

k/U*2

Fig. 4. The comparison of simulation results with experimental and numerical data from Prinos et al. (2003). The vertical profile of mean velocity magnitude for (a) Case 250-30 and (b) Case 250-50. The vertical profile of turbulent kinetic energy for (c) Case 250-30 and (d) Case 250-50.

- 46 -

(a)

(c)

(d)

δk

y/hp

Up/Uf

y/hf

(b)

x 103

x 103

Ref

U/U*

Ref

U/U*

k/U*2

k/U*2

Fig. 5. The vertical profile of mean velocity magnitude for porosity of (a) 0.5 and (b) 0.8. The vertical profile of turbulent kinetic energy for porosity of (c) 0.5 and (d) 0.8. The inset in (b) and (d) indicate ratios of mean subsurface velocity to mean surface velocity and penetration depths of turbulent kinetic energy, respectively as a function of Re f .

- 47 -

(a)

-0.07

y/hp

-0.14

TKE penetration

2.00 1.67

Vortex size 0

1.33 1.00

-0.21

U/U*

-0.14 -0.21

-0.28

2.00 1.67

-0.28

-0.35

1.33 1.00

-0.35

-0.42 -0.50

TKE penetration

-0.07

0.67 0.33 0

y/hp

Vortex size 0

(b)

k/U*2

0.67 0.33 0

Preferential flow path

-0.42 -0.50

Fig. 6. The velocity map with the recirculation zones and the map of turbulent kinetic energy (TKE) near the free-flow-porous media interface for (a) Case 5010 and (b) Case 8010. The red dashed boxes indicate extent of the first vortex developed immediately below the freeflow-porous media interface.

- 48 -

(a)

PDF

PDF

PDF

(b)

t/(L/Uin)

t/(L/Uin)

t (s)

t (s)

Fig. 7. The breakthrough curves at L  3 m for porosity of (a) 0.5 and (b) 0.8 with the variation of Re f . The dotted lines indicate breakthrough curves of particles traveling only through a free-flow layer. The insets indicate breakthrough curves normalized with inlet velocity.

- 49 -

%Msub (%)

Converging to 16% Converging to 9%

t/(L/Uin) Fig. 8. The temporal evolution of the percent mass of injected particles accumulated in subsurface regions during the simulation period.

- 50 -

(b) c

(m)

(a)

PDF

x 103

Ref

~t-1.8

~t-1.8

~t-3.2

~t-3.2

t/(L/Uin)

t/(L/Uin)

Fig. 9. The influence of porosity on (a) breakthrough curves at L  3 m and (b) corresponding transition time distributions of particles. The dotted lines in (a) indicates breakthrough curves of particles reaching the downstream location through a free-flow zone. The inset in (b) presents correlation length scales as a function of Re f .

- 51 -

next transition time class (Ci)

log10(PDF) t = 60 s

U/U*

(a)

3.0 2.5 2.0 1.5 1.0 0.5 0

next transition time class (Ci)

t = 60 s

(b)

current transition time class (Cj) Fig. 10. The transition matrices of particle transition times with snapshots of particle trajectories at t = 60 s for (a) Case5010 and (b) Case8010. The dashed boxes in (a) and (b) indicate regions of next transition time classes corresponding to the flow region directly above the free-flow-porous media interface, and the flow region of preferential flow paths in porous media, respectively.

- 52 -

t / ( L / U in )

PDF

~t-1.8

~t-3.2 x 103 Ref

t/(L/Uin) Fig. 11. Residence time distributions of particles captured in recirculation zones for Case 5010 and Case 8010. The inset indicates mean particle residence times inside recirculation zones as a function of Re f .

- 53 -

(b)

PDF

(a)

t (s)

Truncation time ≈ 480 s

t (s)

Truncation time ≈ 240 s

Fig. 12. The comparison between breakthrough curves obtained from 2D Lagrangian particle tracking (LPT) simulations and from a Spatial Markov model (SMM) for (a) Case5010 and (b) Case8010. The breakthrough curves are estimated at x /L = 0.17, 0.33 and 1.00. The inset in (b) shows breakthrough curves predicted at x /L = 1.00 by the SMM using full data (simulation period of 300 s) and partial data (simulation period of 200 s) of transition times (  ) as well as by an uncorrelated continuous time random walk (CTRW) model.

- 54 -