Hyporheic exchange in recirculating flumes under heterogeneous bacterial and morphological conditions

Hyporheic exchange (HE) contributes to the biogeochemical turnover of macro- and micro-pollutants in rivers. However, the spatiotemporal complexity and variability of HE hinder understanding of its role in the overall functioning of riverine ecosystems. The present study focuses on investigating the role of bacterial diversity and sediment morphology on HE using a multi-flume experiment. A fully coupled surface–subsurface numerical model was used to highlight complex exchange patterns between surface water and the underlying flow field in the sediments. Under the experimental conditions, the surface water flow induced by bedforms has a prominent effect on both local trajectories and residence time distributions of hyporheic flow paths, whereas mean hyporheic retention times are mainly modulated by average surface flowrates. In case of complex bedform morphologies, the numerical model successfully reproduces the HE estimated by means of salt dilution tests. However, the 2D numerical representation of the system falls short in predicting HE in absence of bedforms, highlighting the intrinsic complexity of water circulation patterns in real scenarios. Finally, results show that higher bacterial diversities in the stream sediments can significantly reduce hyporheic fluxes. This work provides a framework to interpret micropollutants turnover in light of the underlying physical transport processes in the hyporheic zone. The study emphasizes the importance of better understanding the tradeoff between physically driven transport processes and bacterial dynamics in the hyporheic zone to quantify the fate of pollutants in streams and rivers.


Introduction
In streams and rivers, hyporheic exchange fluxes between surface water and sediment porewater are driven by a variety of physical processes playing on different temporal and spatial scales (Vogt et al. 2010;Stonedahl et al. 2010;Stanford and Ward 1988;Boano et al. 2014;Lewandowski et al. 2011;Krause et al. 2012). The biogeochemical cycling of nutrients and pollutants mediated by microbial communities in the hyporheic zone is critically affected by the coupled transport processes between surface turbulent flows and laminar flows in the sediment ). Hyporheic flow processes control flow path distributions, exchange rates, exchange volumes, travel and residence times of water (and waterborne compounds) across the riverbed (Bencala and Walters 1983;Wörman et al. 2002;Marion et al. 2008;Mojarrad et al. 2019a). These variables bear a substantial signature on the biogeochemical cycling of nutrients and pollutants in the sediments and on the development of 1 3 234 Page 2 of 18 microbial and macroinvertebrate communities therein (Fischer et al. 2005;Boano et al. 2010;Marzadri et al. 2010; Bardini et al. 2012;Duff and Triska 2000;Battin et al. 2008;Peralta-Maraver et al. 2018).
In the hyporheic zone (HZ) strong redox and temperature gradients resulting from the interaction between surface water and groundwater provide a unique chemically and biologically active environment . Microbial activity, together with chemical and physical processes such as redox reactions, precipitation, sorption, settling and filtration, affect solute concentration and transport in complex feedback loops (Bolton et al. 1999;Brunke 1999;Packman and MacKay 2003;Caruso et al. 2017). As a result of these interwoven processes, hyporheic zones play an important role in improving water quality in streams and rivers (Boulton 2007;Boulton et al. 2010;Mulholland et al. 2008;Grimm and Fisher 1984). The ecosystem services associated with hyporheic exchanges have encouraged significant research efforts that, over the last few decades, have provided new mechanistic and phenomenological insights on the interactions between surface and subsurface waters .
A variety of approaches have been adopted to study the complex physical and biochemical functioning of the hyporheic zone. For example, field observations of hydrological and biochemical variables have been combined with analytical and conceptual lumped hydrological models to estimate hyporheic exchange (e.g. Bencala and Walters 1983;Wörman et al. 2002;Marion et al. 2003). The calibration of physically meaningful parameters with the aid of observed breakthrough curves of conservative and non-conservative tracers allowed models to provide valuable insights on reach-scale exchange processes and biochemical cycling (Chakraborty et al. 2009;Marion et al. 2008). However, the interpretation of exchange rates based on field data can be cumbersome as rivers integrate a variety of processes that are often characterized by high temporal and spatial variability (Jaeger et al. 2019b;Ren et al. 2019a, b). As a consequence, it is often difficult to decouple the contribution of different physical and biogeochemical processes on the turnover of nutrients and pollutants and to resolve the temporal and spatial heterogeneity of degradation rates along river reaches.
To address this problem, hyporheic exchange has been studied in laboratory or outdoor flumes where the environment can be controlled more easily (Fox et al. 2014;Endreny et al. 2011;Eylers et al. 1995;Kim et al. 1990). In such settings, it is possible to target specific research questions by reducing the degrees of freedom of the system.
Analytical flow models are also valuable tools to quantify exchange rates between surface and subsurface waters (Elliott and Brooks 1997;Kazezyilmaz-Alhan 2008;Marzadri et al. 2011;McCallum et al. 2012;Worman 1998). Their neat physical-based formulation enables the assessment of the sensitivity of important variables such as exchange fluxes and residence times to boundary conditions and geometric constraints of the system. However, closed form solutions of flow and transport in porous media are generally limited to simplified geometries and/or boundary conditions.
In recent years, the widespread improvement of computational power, combined with advanced fluid dynamics algorithms, have opened new avenues for a refined description of hyporheic flow processes also for non-stationary conditions in complex morphological configurations (Cardenas and Wilson 2007;Boano et al. 2009;Sheibley et al. 2003;Runkel and Chapra 1993;Saenger et al. 2005;Mojarrad et al. 2019b;Al-Mansori et al. 2020). Provided that they can be constrained by reliable boundary conditions and supplemented by suitable constitutive relations, numerical models are powerful aids in describing the physical processes that modulate the environmental cycling of pollutants (Ren et al. 2019a, b;Caruso et al. 2017).
Nowadays, also the use of mesocosms experiments employing flumes is becoming popular for studying the biogeochemical turnover of pollutant mediated by hyporheic processes. However, it is challenging to quantify HE without disturbing the sediment in an ongoing experiment and non-stationary physical and biological boundary conditions during long-term studies substantially challenge the interpretation of experimental results. In fact, microbial communities in the stream sediments can substantially alter surface-subsurface water exchanges (Battin and Sengschmitt 1999). Numerical models have shown that bioclogging can significantly affect the hydraulic properties of stream sediments and, in turn, modulate the fluxes and retention times in the hyporheic zone (Caruso et al. 2017). Nonetheless, no systematic assessment has yet been attempted to quantify the role of microbial diversity on HE fluxes in large experimental contexts.
This work combines numerical simulations and experimental observations to characterize HE in a set of 22 recirculating flumes (see Fig. 1 for a research outline). The study is part of a larger experiment which investigated-from a biochemical perspective-the effect of bacterial diversity and sediment morphology on the removal rates of micropollutants (Jaeger et al. 2019a;Posselt et al. 2020). Since an accurate characterization HE is a prerequisite to interpret the observed cycling of pollutants, the following analyses will be crucial to quantify the impact of flow processes on the biochemical turnover of waterborne compounds in forthcoming studies. In fact, despite the importance of surface-subsurface water interaction in enhancing biogeochemical cycling, hyporheic exchange processes are usually not assessed in mesocosms experiments investigating the turnover of micropollutants. The current study addresses this gap and generates knowledge on flow path and residence time distribution that can be integrated into microbiological and chemical experimental results, eventually impacting on the design of future mesocosm studies. Specifically, the current study primarily aims to: 1. Characterize hyporheic exchange processes in the experiment and quantify flow path and residence times with a specific focus on their spatial variability. 2. Constrain and evaluate the capacity of a fully coupled 2D numerical model to properly describe HE in experimental conditions. 3. Evaluate long-term changes in hyporheic fluxes and discuss the effect of bacterial diversity on HE.
Numerical simulations of hyporheic flows generally rely on pre-established pressure distributions assigned on the upper sediment interface (e.g. Caruso et al. 2017;Salehin et al. 2004). Pressure distributions are mostly obtained via the semi-analytic approach of Elliott and Brooks (1997) as a function of key morphological and surface flow features (e.g. bedform wavelength and amplitude, average surface flow rate, water depth). However, because of the complexity of the processes controlling the interaction between surface and subsurface exchange, a fully coupled approach was adopted in this study. In fact, directly coupling surface and subsurface flow fields allows feedback interactions between the two flow domains to be accounted for, providing a more accurate description of the forcing excreted by complex surface flows on the underlying transport in the porous media.
This study is organized as follows: the next section presents the experimental design and the methods adopted to acquire, analyze and interpret the data, followed by which the numerical representation of the experimental setup and its parametrization are introduced. The two subsequent sections present the experimental and the numerical results, respectively. The experimental and numerical results are then jointly discussed. The final section concludes the study.
In each flume, the sediments were arranged according to three alternative geometries: (i) 0 bedforms (i.e. flat sediment surface); (ii) 3 bedforms on a single side of the flume; (iii) 6 bedforms (3 on each side of the circular flume). Based on the increasing number of bedforms (BFs), the three  Table SI.1 for a summary of the notation used in this study). Three levels of bacterial diversity were tested by mixing clean sand (oven-dried at 120 °C for 24 h) with increasing amounts of sediment extracted from the riverbed of the Erpe, a creek located in Berlin, Germany (52°28′31.9″N 13°37′46.6″E). The Erpe is heavily impacted by treated wastewater as its streamflow can be up to 80% constituted by the urban effluents of the Münchehofe wastewater treatment plant. The plant has a capacity of 220,000 population equivalents and reaches the tertiary stage of wastewater treatment (Lewandowski et al. 2011;Jaeger et al. 2019b). The three levels of bacterial diversity (denoted as S1, S3 and S6) correspond to a dilution ratio of Erpe sediments to commercial sand of 1:10 1 , 1:10 3 and 1:10 6 , respectively.
The setup follows a "dilution-to-extinction" approach. The approach is based on experimental evidence proving that, after an initial incubation phase allowing re-growth of bacteria, the taxonomic and functional richness of bacterial communities directly relate to their initial concentration (Stadler et al. 2018;Jaeger et al. 2019a). The method is a standard in microbiology and it has been successfully applied in similar studies (Stadler et al. 2018;Ylla et al. 2013;Peter et al. 2011;Szabó et al. 2007). Real-time PCR and sequencing of bacterial DNA in sediment conducted during the final phase of the pre-incubation period confirmed the reduced bacterial diversity in flumes characterized by higher initial dilution rates of Erpe sediments (Jaeger et al. 2019a). Two flumes containing a 1:10 mix with sterilized Erpe sediments were additionally included as controls.  Table 1, whereas the specific characteristics of single flumes are reported in the table included in Fig. 3.
In each flume, surface water flow was induced by a 6-cm axial propeller (NWA 1.6 adj 2.6 W, Newa Wave Industria, Loreggia, Italy). The flumes were located inside a white tent to minimize the interaction between the experimental setups and the external environment (e.g. rainfall, debris, wild animals) and to maximize the scattering of direct solar radiation (i.e. flumes are expected to be homogenously irradiated). The average air temperature during the experimental period was 16.3 °C (standard deviation: 3.7 °C).
After a 12-day incubation phase (i.e. a period of bacterial acclimation), a set of 31 micropollutants was injected and monitored for 78 days on both the surface and in the porewater.
To evaluate biochemical cycling in the sediments, single porewater samplers (Standard Rhizons, Rhizosphere Research Products B.V., Wageningen, The Netherlands) were deployed in selected bedforms (Fig. 2a). Additionally, the first (i.e. upstream) bedform in flume 18 (B3, S1) was equipped with three equally spaced (4 cm apart) porewater samplers located 1.5 cm above the bottom of the flume (see Fig. 2b). Porewater was sampled at these locations for additional analysis on micropollutants concentrations . Further details on the experimental design can be found in Jaeger et al. (2019a).

Estimating sediment hydraulic conductivity
Hyporheic fluxes-together with water retention times in the sediments-are critical controls for biochemical cycling   Summary of the variable levels and replicates in the experimental design. The three levels of the bedform morphology (B) and bacterial dilution (S) are denoted as: B0 (flat sediment); B3 (three bedforms); B6 (six bedforms); S1 (dilution 1:10 1 ); S3 (dilution 1:10 3 ); S6 (dilution 1:10 6 ).  in both experimental and in environmental settings (Arnon et al. 2013;Boano et al. 2014;Duff and Triska 2000;Battin et al. 2008;Runkel 2007;Krause et al. 2009). The hydraulic conductivity of the sediments modulates hyporheic fluxes and, in turn, it controls water retention times in the hyporheic zone. Homogeneous and isotropic hydraulic conductivity is expected in the experimental design as the sediments were thoroughly mixed and mainly constituted by commercial sand.

Number of flumes
The hydraulic conductivity of the sand used in the experiment was estimated through four falling head tests. The tests were performed using a vertical Plexiglas tube (5 cm diameter) partially filled with 30 cm of the experimental sediment mixture (the sediments were retained at the bottom of the device by a thin geotextile net). The remaining upper portion of the tube (about 200 cm) was initially filled with clean water. The hydraulic conductivity K s was evaluated by fitting the observed water level in the seepage device by means of Eq. 1: where h(t) is the water level at time t (measured from the bottom of the sediments), h 0 is the water level at t = 0 and L is the thickness of sediments.

Estimating average surface flow velocity
Estimating average flowrates in the experimental flumes is critical for establishing reliable boundary conditions in the numerical model. Given the experimental setup, it was not possible to accurately modulate the flowrate (flumes are closed loops and the pumps operated at a constant rate). Consequently, flow velocity is expected to be affected by hydrodynamic energy losses due to the different geometrical configurations. In particular, as the number of bedforms increase, flowrates are expected to decrease.
The following approach was adopted to account for the effect of the different morphological conditions on the average surface flowrate. Surface flow velocities were estimated in the B0 and B3 settings by measuring the time interval required by a floater to travel along one side of the flume. Measures were performed on the bedform-free sides of the B3 flumes to minimize the disturbing effect of bedforms on the flow field. Measured surface flow velocities were then used to constrain a logarithmic vertical flow profile. The logarithmic vertical flow profile is described by the law of the wall (Schlichting and Gersten 2017): where y is the vertical distance from the sediments, u is the shear velocity (defined as u = √ w ∕ , where w is the shear (2) u(y) = 1 u ln yu + C + u , stress at the sediment-water interface and is the density of water), is the kinematic viscosity of water, and and C + are constants (Schlichting and Gersten 2017). The kinematic viscosity in Eq. 2 is evaluated based on its dependence on the temperature (a reference temperature of 17 °C is assumed), whereas = 0.41 and C + = 5 are set to typical literature values (Schlichting and Gersten 2017). The shear velocity u is obtained by implicitly solving Eq. 2 after imposing u B0 y = y max = u B0 and u B3 y = y max = u B3 , where y max is the maximum water depth on the bedform-free side of the flumes and u B0 and u B3 are the measured surface flow velocities in the 0 and 3 bedforms settings. Finally, the depth-averaged flow U is computed by integrating u(y) along the vertical profile as U = 1∕y max ∫ y max 0 u(y)dy. In the B6 experimental setup, it was necessary to adopt an alternative approach to estimate the average flowrate (it was not possible to measure "undisturbed" surface flow velocities due to the presence of dunes on both sides of the flumes). In these cases, it was assumed that the hydraulic potential experienced two types of energy losses. The first (E 1 ) combines distributed energy losses (i.e. friction with the sediments and walls of the flumes) with localized energy losses in correspondence with the bends of the flume. The second (E 2 ) is due to localized energy dissipations caused by the expansion and contraction of the flow field in correspondence with the bedforms. The total energy loss in the B0, B3 and B6 setups ( E B0 , E B3 , E B6 ) can be expressed as a combination of the two components of the energy losses E 1 and E 2 as Since the power of the pumps is equal in all three morphological settings, it can be assumed that the total loss of energy per unit of time equals the power of the pump. As a consequence, Assuming ũ = 5 cm s as a first-order approximation of the flow velocity in the flumes, the corresponding Reynolds number for an open channel flow reads Re =ũ4D∕ ≈ 10 5 (Streeter 1962), where and are the density and the viscosity of the water at 17 °C and D is the hydraulic radius of the cross section of the flumes (5 cm). As the Reynolds number is larger than 2900, turbulent flow is established (Schlichting and Gersten 2017). Since turbulent energy losses scale with the squared flow velocity (Streeter 1962;White 1999), it is possible to express the two components of the energy losses E 1 B * and E 2 B * in Eq. 3 as where * refers to the considered bedform morphology (i.e. B0, B3 or B6) and and are constants. By combining and rearranging Eqs. (3), (4) and (5), it is possible to calculate the average flow velocity U 6 as a function of the average flows U 0 and U 3 as

Breakthrough of acesulfame
After an initial incubation phase of 12 days, a mixture of micropollutants-including acesulfame-was injected in the surface water of each flume aiming at a concentration of 10 µg L −1 per compound. The concentrations of micropollutants were repeatedly measured using three porewater samplers located in the first bedform of flume 18 (B3, S1; see "Recirculating flumes"). Acesulfame is a widespread artificial sweetener that is ubiquitously found in treated wastewater. Since acesulfame has been shown to be unaffected by transport retardation (due for example to sorption processes on the porous matrix (Schaper et al. 2019;Jaeger et al. 2019a)) it can act as a proxy to quantify advective flow in the HZ. As such, the breakthrough of acesulfame across the sampling locations in the bedforms provided a first-order estimate of the timescales of advective hyporheic transport and, in turn, it provided useful insights to benchmark the numerical simulations.

Salt dilution tests
During the final phase of the experiment, salt dilution tests were conducted to quantify the average water exchange fluxes between surface water and sediments (Arnon et al. 2007;Packman et al. 2004;Mutz et al. 2007;Fox et al. 2014;Galloway et al. 2019). Specifically, 60 g of sodium chloride (NaCl) were diluted in the surface water of the flumes. Electrical conductivity EC (i.e. a proxy for chloride concentration) was automatically sampled from the surface water (CTD-Diver, van Essen Instruments, Delft, the Netherlands) at 5 min temporal resolution in all flumes except for flumes 11-20 where EC was manually sampled. All measurements were corrected by 5-step calibrations of EC loggers and meters following the termination of the test. NaCl is commonly considered as a non-reactive tracer for solute transport. Specifically, background EC, temperature and interaction with other chemicals is assumed to be negligible given the large increment in EC during the test (EC ranged from about 1.1 mS cm −1 to more than 3 mS cm −1 upon NaCl addition). The exchange rate between surface water and porewater can be evaluated by measuring the temporal decrease of the electrical conductivity in the surface water. In fact, through hyporheic exchanges, surface water is forced into the sediments while water with lower salinity (i.e. lower EC) is pushed out from the sediments. As a consequence, the EC in the surface water tends to decrease until equilibrium. Equilibrium is reached when a homogeneous concentration of NaCl is established across both the surface water and the porewater affected by HE. The decrease in the surface water concentration of NaCl can be analytically described by an exponential function of the type (Mutz et al. 2007;Salehin et al. 2004): where C(t) represents the EC of the surface water at time t, C 0 and C eq are, respectively, the initial EC ( C 0 = C(t = 0) ) and the equilibrium EC ( C eq = C(t → ∞)) , and k [1/time] is the exchange rate constant between surface and sediment water. The parameters k and C eq can be obtained by fitting Eq. 7 to the decays of EC observed through time.
From a mass balance accounting for the equilibrium concentration C eq , it is possible to estimate the volume of the porewater V s affected by hyporheic exchange as where V w is the surface water volume. Furthermore, assuming that during the inception phase of the salt dilution test the amount of NaCl exiting the sediments is negligible compare to the amount of NaCl that is advected into the sediments, the following mass balance equation can be established: where Q in is the flux of water entering the sediments. Substituting Eq. 7 and its derivative into Eq. 9, the latter reads: By taking the limit of Eq. 10 for t → 0 (initial phase of the dilution test), it is possible to quantify the exchange flux between surface water and sediment water as Under the assumption of perfect mixing of water in the sediments or, alternatively, in case of homogeneous residence time of the advective porewater flow through the sediments, a first-order estimate of the hyporheic retention time can be obtained as

Inputs and parametrization of the numerical model
Surface flow velocity and hydraulic conductivity of porous media are critical variables modulating hyporheic exchange (Ren et al. 2019a, b). Flow velocity controls the hydrodynamic potential at the interface between surface water and sediments, in turn imposing the internal boundary condition across the two flow domains (Elliott and Brooks 1997;Bottacin-Busolin and Marion 2010). In fact, hydrostatic and hydrodynamic pressure gradients, jointly with turbulent momentum transfer between surface water and porewater, are the main drivers of the flow field in the sediments (Boano et al. 2007;Elliott and Brooks 1997;O'Connor et al. 2012;Taylor 1954;Nagaoka and Ohgaki 1990;Dade 2001). Hyporheic exchanges have nonlinear dependence on surface flow velocity as hyporheic exchange driven by both advective pumping and turbulent momentum transfer tends to scale as a quadratic function of surface flow velocity (Elliott and Brooks 1997;Packman et al. 2004;O'Connor and Harvey 2008;Arnon et al. 2007;Salehin et al. 2004). The hydraulic conductivity of sediments (K s ) is also key in controlling hyporheic exchanges since flow velocities in porous media linearly depends on K s (O'Connor et al. 2012;Elliott and Brooks 1997;Al-Rasool Ali and Shakir 2019;Fazelabdolabadi and Golestan 2020). Hence, the estimates of K s and of the average surface flow rate obtained as described in the previous sections are crucial to constrain the numerical representation of the system implemented as described hereafter.

Coupled surface water and pore-water flow model
A 2D numerical model was implemented to couple the surface flow field (described by the Navier-Stokes equations with L-VEL turbulence closing model) with the underlying flow in the porous media (described by the Brinkman equations) (Ren et al. 2019a;O'Connor et al. 2008). The Brinkman equations account for the laminar momentum transfer in terms of kinetic and viscous forces at the interface between the two flow domains. They represent an extension of Darcy's theory, which assumes that flow in a fully saturated porous medium is simply driven by pressure gradients (Amiri and Vafai 1998). In the following analyses, the flow equations were further extended to account for turbulent momentum transfer between the surface flow and the porous domain by accounting for the Forchheimer drag in the Brinkman equation (Amiri and Vafai 1998). The model was implemented in COMSOL ® Multiphysics, a finite element platform designed to solve and couple complex systems of partial differential equation describing different physical processes across multiple domains (www. comsol. com). The fully coupled flow field between surface water and porewater was simulated on a single side of the recirculating flumes (in the B3 case the side with bedforms). The geometry of the system was constructed to accurately reproduce the experimental layout (see Fig. 2 and Fig. SI.1), whereas the hydraulic properties of the sediments include the hydraulic conductivity (estimated as described in "Estimating sediment hydraulic conductivity") and porosity (n = 0.35). A dynamic viscosity of 1.08 mPa · s and a density of 1.0 g/cm 3 were assigned to water. The Forchheimer term in the Brinkman equation 150nk is evaluated based on the (known) values of water density ( ) , sediment porosity (n), and sediment permeability (k).
The boundary conditions adopted in the numerical model include a no-slip constraint on the lower boundary (i.e. in correspondence of the interface between the sediment and the bottom of the flume) and a fixed inlet and outflow. Inlet and outlet flows were set to the corresponding average flows U 0 , U 3 and U 6 obtained as described in "Estimating average surface flow velocity". To avoid convergence instability related to the high degree of freedom on the problem, a fixed upper boundary condition (no stress, slip condition) was applied in correspondence with the atmospheric boundary. The assumption is justified by the reduced longitudinal variability of the free surface elevation that was observed during the experiment.
Computations have been performed on a flow domain discretized using a free triangular mesh having a minimum (maximum) element size of approximately 0.3 (2.0) mm.

Analysis of the hyporheic flow paths
To reproduce and evaluate the breakthrough curves of acesulfame observed as described in "Breakthrough of acesulfame", hyporheic flow paths were simulated by combining the results of the numerical flow model with a specifically developed inverse tracking algorithm. The aim of the analysis was to evaluate the ability of the numerical model to reproduce the time-to-peak of the breakthrough of acesulfame at the different sampling locations and to estimate the uncertainty in the predicted flow path length and travel times in different bedforms throughout a stochastic approach. The agreement between model estimates and experimental observations is a requisite for using the numerical model to interpret the biochemical cycling of micropollutants observed in Jaeger et al. (2019a). Note that the implemented approach aims to capture the advective component of the transport process, which generally dominates over dispersion in case of bedform-induced HE (Elliott and Brooks 1997;Salehin et al. 2004;Savant et al. 1987).
The experimental sampling points in the first bedform of flume 18 (see sections. "Recirculating flumes" and "Breakthrough of acesulfame") were recreated in the numerical model. Moreover, to compare the trajectories of hyporheic flow paths in downstream bedforms, three additional sampling locations were simulated in the second bedform.
The inverse tracking algorithm adopted in the analysis backtracks the trajectories of the water parcels reaching the sampling locations by accounting for uncertainties associated with the knowledge of the flow field and/or of the precise sampling locations. The framework allows assessment of how sampling different flow paths in the surrounding of the expected sampling points can affect flow path lengths and travel times (or, in other words, how the variance in sampling locations propagates to the variance of flow path length and travel times).
Specifically, the method tracks 30,000 particles seeded at the sampling points (10,000 for each sampling location) back until they enter the sediments. The particles are initially seeded around the sampling location according to a bivariate normal distribution with horizontal and vertical variance, respectively, of 2 x = 5 mm 2 , 2 y = 2.5 mm 2 . The variances 2 x and 2 y are assigned based on a first-order estimate of the uncertainties expected on the precise location of the porewater samplers. Trajectories are delineated by identifying the upstream nodes on the flow field mesh (obtained as described in "Coupled surface water and pore-water flow fields") whose seepage velocity leads to the initial seeding point. The algorithm is recursively repeated until the surface of the sediments is reached. The trajectory of each particle is computed together with the corresponding travel time and flow path length.

Experimental results
The hydraulic conductivity (K s ) and permeability (k) obtained from the falling head experiments were 2.35 · 10 -4 m/s and k = 2.40 · 10 -11 m 2 , respectively (see Sect. SI.2 for details), while the average surface flow velocities in the three morphological configurations B0, B3 and B6 was 8.7, 6.6 and 5.5 cm/s, respectively (see Sect. SI.3).
The concentration of acesulfame peaked after about one day in the first (A) and second (B) sampling location (see Fig. 2). At sampling location C the breakthrough of acesulfame is more spread and it is peaking after about two days (see Sect. SI.4).
The decay of electrical conductivity during the initial phase of the salt dilution test is shown in Fig. 3. Overall, the analytical model captures the initial phase of the salt dilution test except for flumes 21 and 22 where the curves tend to overestimate the very initial dilution. Figure SI.5 (supporting information) summarizes and highlights an enhanced inter-flume variability in terms of hyporheic exchange (flux and volume) and average water retention time across the different bacterial and morphological configurations.
The exchange variables (mean fluxes, exchange volumes and residence times) for each combination of bedform morphology and bacterial diversity are aggregated in Table 2. Notably, the mean exchange fluxes tend to be lower in S1 settings (lower sediment dilution, i.e. larger bacterial diversity) regardless of streambed morphology (Table 2a). On the other hand, the mean residence times are shorter in flumes with lower sediment dilution (Table 2c). Furthermore, data suggest a moderate increase in exchange volumes and mean residence times with increasing morphological diversity (as well as a decrease in mean exchange fluxes).
Analysis of variance confirms a significant difference in the hyporheic fluxes in S1 compared to S3 and S6 settings (P ANOVA = 0.025), whereas no statistically significant differences are found in HE metrics for different morphological conditions (P ANOVA > 0.05). However, despite the relatively large number of replicates in the experiment (i.e. flumes with different combinations of morphological and bacterial conditions), statistical significance (as quantified by ANOVA) might be affected by the available sample sizes. Figure 4 aggregates the theoretical fits of the initial phase of the salt dilution tests as a function of bacterial diversity in the flumes (i.e. decreasing sediment dilution). Figure 4 highlights how bacterial diversity modulates the flow between surface and porewater. Flumes with higher bacterial diversity in fact display slower decays of EC through time (i.e. reduced hyporheic exchange) compared to flumes featuring lower bacterial diversity (see also Table 2). However, the relationship between HE and bacterial diversity appears to be nonlinear as a threshold seems to exist between the larger bacterial diversity (S1) and the lower bacterial diversities (S3 and S6). In fact, no significant difference appears between HE fluxes in S3 and S6 conditions (P ANOVA = 0.97).

Numerical model: results
Coupled surface water and pore-water flow fields Figure 5 shows the surface water flow field in B0, B3 and B6 settings (note that in B0 the entire flume was simulated to quantify the effects on hyporheic exchange induced by the initial discontinuity of the sediment cover). In case of B3 and B6, the stretch of the flume featuring the bedforms was simulated.
Results on the flat sediment bed were as expected (Fig. 5a), featuring a homogeneous longitudinal flow field and a vertical velocity gradient consistent with the law of the wall (see "Estimating average surface flow velocity" and Fig. SI.3).
On the other hand, bedforms enhanced spatial gradients in both surface and subsurface flow fields, especially in B3 settings (i.e. larger surface flow rate compared to B6 settings, see Fig. 5b, c). exchange pattern on both sides in the B6 setting, the total exchange fluxs in the two case results: Q model B3 = 11.1l∕d and Q model B6 = 2 ⋅ 6.8 = 13.6l∕d. The first assumption is justified due to the negligible exchange fluxes predicted by the numerical model in case of a flat sediment surface (i.e. B0 setting, see Fig. 5a), while the second assumption is motivated by the symmetrical configuration of the experimental layout in the B6 case. Note that turbulence and pressuredriven hyporheic fluxes in correspondence with the curves of the flumes are not accounted by the two-dimensional flow model. The longitudinal exchange velocity profile is not shown in case of flat sediment as exchange was small ( Q model B0 = 0.4 l∕d) and limited to the initial (and final) discontinuity in the sediment cover.

Analysis of hyporheic flow paths
The simulation of flow paths for the 30,000 particles backtracked from the modelled sampling locations in the first two bedforms revealed different trajectories in the first compared to the second bedform (Fig. 7). In contrast, no considerable difference can be noted in their trajectories as the surface flow velocity varies (i.e. in B3 vs B6 settings). In fact, the distribution of the lengths of the flow paths reaching the three sampling locations in the B3 and B6 cases displays small variability (see Sect. SI.6 in the supporting material). Figure 8 shows the distribution of the advective travel times in the sediments. Unlike for flow path length, the morphological configuration of the flumes (B3 vs B6) and the bedforms (first vs second) strongly modulate both the average and the variance of the residence time distributions (see Fig. 8). In the B6 setting-i.e. where surface flow velocity is lower-travel times are on average longer and more heterogeneous (they are featured by a larger variability) compared to the B3 case (Table 3). Furthermore, it can be noted how travel time variability scales with the average travel time and travel times are shorter in the second bedform (for both the B3 and B6 settings). Interestingly-unlike average travel times (which are consistently shorter in the second bedform  Fig. 7 Simulated trajectories of the water parcel reaching the three sampling points (denoted by crosses) located in the first and in the second bedform, for both the B3 and B6 case. Flow paths were estimated by backtracking 10,000 particles normally distributed around each sampling location for both B3 and B6 settings)-the standard deviations of the travel times distributions tend to be smaller (larger) in the first bedform of B3 (B6) configurations (see Table 3). Larger variabilities of traveltime distributions in downstream bedforms in B3 settings can be a consequence of the enhanced turbulence in the surface water induced by upstream bedforms combined with a larger surface flow rate. The modelled average travel time to the sampling locations A, B and C in the first bedform of B3 settings are 15.3, 24.3 and 43.3 h, respectively (see Fig. 8 and Table 3). The estimates correspond to observed time to time-to-peak in the breakthrough of acesulfame of about 1 day for sampling locations A and B and 2 days for sampling location C (see Sect. SI.4 in the supporting materials). Although the temporal resolution of the observed breakthrough of acesulfame was not sufficient to clearly differentiate the times-to-peaks at the first two sampling locations, the numerical simulations are compatible with the empirical data and suggest reduced travel times in case of sampling location A.

Discussion
The analytical relation used to estimate the hydraulic conductivity of the sediments accurately represents the fallinghead tests (Fig. SI.2). Interestingly, the observed decay rates (i.e. K s /L, see Eq. 1) tend to decrease slightly as the tests are repeated, possibly indicating a moderate reduction of K s due to sediment compaction in the first phase of the seepage tests (clogging is less likely as clean water was used).
Surface flow velocity in the experimental flumes tends to decreases as the number of bedforms increases (Fig. 5 and Sect. SI.3). This is an expected consequence of the energy losses localized in correspondence of the bedforms. Flow estimate obtained as described in section "Estimating average surface flow velocity" for B6 settings is appropriate. In fact, it corresponds to the average flow rate that would be required to match the surface flow velocities observed over the crests of the bedforms (about 10 cm/s). Figure 5 does not highlight any considerable flow field in the sediments in the B0 settings. On the other hand, Fig. 5 shows that in presence of bedforms the porewater flow fields become more intense as surface flow velocity increases (cf. B3 and B6 configurations). The numerical simulations highlight how surface flow velocity sharply rises over the bedforms and it mildly reduces in between bedforms. This indicates that upstream bedforms affect the flow field impacting the subsequent bedforms and  hydrodynamic pressure hits mainly the crest of downstream bedforms. As a result, the flow field in the first bedform is better established across the entire bedform reaching larger depths (i.e. longer flow paths), whereas hyporheic fluxes affect mainly the crest of the downstream bedforms.
Although the first bedform is more exposed to the flow, Fig. 6 shows that the maximum exchange flux remains stable across bedforms in flumes sharing the same morphological configuration (i.e. B3 or B6) and it is located at the bedform's crests (upstream: downwelling; downstream: upwelling). It can be noted how hyporheic exchange scales with flow velocity (while exchange patterns are preserved). Figure 6 additionally highlights a complex longitudinal hyporheic exchange pattern. The entire upstream facing slope of the first bedform is affected by downwelling. On the contrary, the upstream face of the following bedforms are affected by more heterogeneous hyporheic circulation patterns (alternation of upwelling and downwelling). This is a consequence of the turbulence in the surface flow induced by upstream bedforms that in turn affect the pressure field on the boundary of downstream bedforms.
Although the stronger (and more coherent) pressure gradient enhances the flow field in the first dune, modelled travel times to the three sampling locations are shorter in downstream dunes (Fig. 8). The unexpected result derives from the different trajectories of the flow paths leading to the sampling location as it can be noted in Fig. 7. In fact, Fig. 7 and Fig. SI.6 show that flow paths are longer in the first bedform. Hence, the position of the bedform has a strong impact on the length (and trajectory) of the flow paths. In particular, the upstream bedform being directly impacted by the flow has deeper flow paths whereas shorter and shallower flow paths are typical of the downstream bedforms. Nonetheless, Fig. 7 and Fig. SI.6 show that differences in surface flow velocity between B3 and B6 have a limited influence on the length of the flow paths when considering the same bedform (i.e. the water parcels tend to follows the same trajectories). Figure SI.6 shows that surface flow rate, bedform ordering and sampling position do not modulate the variability of the flow path lengths (i.e. variability is similar across different sampling locations and different bedforms for different flow rates). On the other hand, Fig. 8 highlights a potentially large variability in the travel times of sampled water, especially for longer flow paths and for reduced hyporheic exchange rates (i.e. lower surface flow rate). Despite the fact that variability of flow path lengths remains comparable between the different settings (i.e. B3 vs B6 settings, first vs second bedform, A vs B vs C sampling locations), the variability of the travel times can vary substantially and generally increases with increasing travel time. Figure 8 shows an increase in both the average and the variability of the travel times as surface flow velocity decreases (i.e. going from B3 to B6). Moreover, travel times are generally shorter in the second bedform (for both the B3 and B6 cases) because of the disturbance in the surface flow induced by the first bedform. This highlights the importance of coupling the surface and the subsurface flow fields when studying specific HE settings. In fact, peculiar patterns in the surface flow field can locally play as critical as the average surface flow rate in constraining hyporheic exchange and the ensuing condition for biochemical cycling in stream sediments.
The agreement between modelled travel times and observed time-to-peaks of acesulfame in the sediments (see "Analysis of hyporheic flow paths") support the assumptions underlying the conceptual and the numerical modelling framework adopted in this study (e.g. boundary conditions, numerical algorithms). The inverse particle tracking analysis suggests a broad range of conditions for biogeochemical cycling in settings with lower surface flow rates because of the large variability of possible travel times to the sampling locations. In fact, a wider set of compounds can find suitable biochemical conditions for being degraded when a wider range of residence times is available. On the other hand, in cases of higher surface flow rate (B3 settings), larger exchange fluxes can enhance the overall yield of targeted processes within a narrow range of residence times.
The mean exchange fluxes for the B3 and B6 configurations derived numerically ( Q model B3 = 11.1l∕d and Q model B6 = 13.6l∕d ) agree with experimental estimates ( Q exp B3 = 7.9 l∕d and Q exp B6 = 7.5l∕d , see Table 2). Experimental data confirm that doubling the number of bedforms does not double the hyporheic exchange (exchange fluxes are comparable in both morphological settings), emphasizing the strong and non-linear sensitivity of HE to surface flow rates.
The reduced fluxes estimated by means of the salt dilution tests (compared to numerical simulations) is possibly a consequence of the clogging of the sediment pores caused by microbial communities (Battin and Sengschmitt 1999;Orr et al. 2009). In fact, the presence of microbial communities has been shown to significantly affect hyporheic exchanges in the experimental settings (see Fig. 4). Additionally, reduced hyporheic exchanges estimated by means of the salt dilution can be due to alterations in bedform topography. The initial experimental geometry was used to design the numerical model. However, smoothed bedforms might have reduced the pumping effects induced by hydrodynamic pressure gradients and, in turn, decreased hyporheic fluxes during the final phase of the experiment. Additional numerical simulations performed to account for the 20% reduction in bedform amplitude observed during the final experimental conditions in B6 settings confirm the decrease of hyporheic exchange in case of smoothed sediment morphology (see Sect. SI.7 in the supporting materials). In this case, the modeled hyporheic exchange ( Q model B6 final = 8.8l∕d ) accurately matches the estimates obtained by means of the salt dilution tests in case of low bacterial diversity ( Q exp. B6,S6 = 9.0l∕d ). The lower HE observed in S1 and S3 settings can be the consequence of a reduced hydraulic conductivities triggered by a higher bacterial diversity therein (see Table 2a).
Interestingly, for flat sediment beds (B0), the salt dilution tests provide exchange fluxes one order of magnitude larger than those estimated by means of the numerical model. This can be due to (i) small riffles that developed during the experiment over the originally flat sediment surface; (ii) the use of a 2-D model which is unable to account for the full complexity of the flow field in the flumes. Unlike the numerical model, salt dilution tests could capture HE along short flow paths across small bedforms and/or in correspondence with the bends. Indeed, residence times estimates in the B0 case are the shortest (compared to the B3 and B6 settings, see Table 2). As a result, comparable hyporheic exchanges can be experimentally observed in case of larger flow velocities over small and dense bedforms and in case of slower flows on larger dunes.
The differences between modelled and experimental estimates of hyporheic fluxes in absence of bedforms also suggest that the simplified boundary conditions and the idealized flow domains that are commonly employed in numerical models can bias the estimates of exchange flows, especially in simple geometrical settings. This can hold also in cases where a coupled approach is adopted and additional terms are included in the flow equations to account for laminar and turbulent momentum transfer across the sediment-water interface. The results confirm that mixing between surface water and groundwater is often underestimated by numerical simulation (Cardenas and Wilson 2007;Hester et al. 2013). Although numerical models are a fundamental tool for process understanding, it should be acknowledged that modelling efforts might fall short in capturing the full complexity of natural processes.
Morphological variability has been shown to bear a consistent signature on the half-lives of acesulfame in the experimental flumes (Jaeger et al. 2019a). In general, more bedforms resulted in enhanced degradation of acesulfame. However, the effect of bedforms was more evident in the flumes featuring low bacterial diversity (Jaeger et al. 2019a, b). Specifically, the effect of bedform-induced degradation decreased for decreasing sediment dilution rates (i.e. for increasing bacterial diversity). Thus, higher bacterial diversity appeared to reduce the effect of morphologically driven removal of acesulfame. The effect is consistent with the observed decrease in hyporheic exchange associated with reduced hydraulic conductivities of the sediments in flumes featuring larger bacterial diversity (Table 2 and Fig. 4).
Interestingly, in S3 settings (i.e. the medium bacterial diversity settings) the B3 flumes displayed longer half-lives compared to the B0 flumes. In these cases, differences in travel time distributions might have played jointly with the presence of different microbial communities (Arnon et al. 2007;Ren et al. 2018). As a result, the short (and fast) exchange patterns across the small riffles that formed over time on the flat sediments might have favored the degradation of acesulfame compared to the longer residence times in the larger bedforms. Field studies have indeed shown that shallow hyporheic exchange patterns characterized by short residence times can be more effective in biodegrading a series of common micropollutants (Schaper et al. 2019). However, the shorter residence times modelled in B3 settings (compared to B6 settings) do not correspond to higher acesulfame degradation rates. In fact, for all bacterial diversities (S1, S3 and S6), the B6 treatments resulted in the shortest half-lives (Jaeger et al. 2019a). In B6 cases, longer (and more heterogeneous travel time distribution, see Fig. 8) might have prevailed over the shorter (and less variable) residence times in the B3 treatment, providing enhanced degradation rates. Indeed, the biochemical turnover of pollutants through the hyporheic zone is the result of the complex interplay between surface-subsurface exchange rates mediated by the probability distribution of water retention times in the sediments (Ensign and Doyle 2006;Marcé and Armengol 2009;Poole et al. 2008;Boano et al. 2010).
The results suggest that the fate of micropollutants can be extremely sensitive to surface-subsurface water exchange patterns. Hence, even moderate variability in critical descriptors of hyporheic exchange (e.g. residence time, flow path length, exchange fluxes) should be accurately quantified because they can substantially impact the turnover of micropollutants, both in experimental and in natural conditions.

Conclusions
This study characterized the physical processes controlling the coupled surface and porewater flow field in a large experiment featuring recirculating flumes by combining experimental results with numerical modelling. The modelling results reproduce the observed hyporheic exchange except for the case of flumes with a flat sediment bed. The discrepancies between modelling and experimental results highlight: (i) the limits of numerical models to reproduce the full range of complex processes featuring surfaceporewater interaction; (ii) the challenges in describing a slowly evolving biological system over a long period of time by means of a stationary description of the underlying hyporheic exchange mechanisms.
The limitations of the analysis are mainly related to the experimental setup (and especially to the recirculating nature of the flumes). The geometrical setup challenged the quantification of the boundary conditions (such as surface flow rates), and the establishment of complex 3D flow fields in the bends may not have been adequately described by the 2D flow model. Sources of uncertainties related to the subsurface flow field can be assessed by means of the backtracking algorithm that was introduced in this study. The tool proves useful to assess the sensitivity of residence time and flow path length to uncertainties in hyporheic trajectories, finding profitable applications in a broad set of contexts. Specifically, for the current experiment, the degradation rates of micropollutant observed at different sampling locations will be interpreted in forthcoming studies in light of the underlying variability that characterizes travel times at each sampling point.
The simulations show how the turbulent flow of surface water plays a critical role on the trajectories of the subsurface flow field and on the variance of the residence time distributions in the hyporheic zone, while average retention times in the sediments are mainly controlled by the average surface water flow. Consequently, the turnover of micropollutants through hyporheic exchanges can display considerable variability between different streambed morphologies. Additionally, the study shows how HE is subject to consistent dynamics during long-term experiments and by comparing experimental and numerical results it is possible to decouple the bacterial from the morphological contribution on the decrease of HE over time.
Finally, the study shows that bacterial diversity can significantly modulate HE. This effect can becloud the contribution of streambed morphology in controlling the interaction between surface water and sediment water. Nonetheless, the complex non-linear dependency between bacterial diversity and hyporheic exchanges calls for additional research. Moreover, the tradeoff between physical transport processes and bacterial dynamics on the ensuing biochemical cycling in the hyporheic zone requires further investigations to better understand the fate of pollutants in streams and rivers.