Surface Transient Storage Under Low‐Flow Conditions in Streams With Rough Bathymetry

In natural and restored rivers, protruding river‐bed elements can cause complex flow patterns. We hypothesize that, particularly at very low flows, solute spreading in such streams is controlled by different mechanisms than at high flows, questioning existing parameterizations. We investigate two‐dimensional flow and transport in a synthetic channel with rough bathymetry under low‐flow conditions with a numerical model solving the full shallow‐water equations. Recirculation cells form behind protruding bed elements and trap water, macroscopically acting as transient‐storage domain. We fit the one‐dimensional mobile‐immobile transport model to the concentration breakthrough curves of the two‐dimensional model and relate the fitted parameters to easy‐to‐measure properties of the stream. The velocity in the mobile domain undergoes a power‐law relation with the discharge and the channel slope, differing from the relationship expected for normal flow in streams with large water depths. Both the fitted dispersion coefficient and the mass‐transfer coefficient between the mobile and immobile domains scale linearly with the velocity in the mobile domain. The ratio of immobile to mobile storage volumes β first increases then decreases with increasing discharge. We suggest measuring the surface area of bed elements protruding the water surface and that of the open water. The product of the corresponding area ratio with the velocity exerts a power‐law control on β. Altogether, we can relate the coefficients of the mobile‐immobile transport model to a few easy‐to‐measure properties of the stream. The findings of the present work contribute to the understanding of surface transient storage in natural and restored streams.

of contaminants and nutrients, being of importance for environmental, biogeochemical, and ecological studies Liu et al., 2020).
The 1-D transient storage model (TSM) proposed by Bencala and Walters (1983) has been extensively applied to describe solute transport in both large rivers and small streams (Anderson & Phanikumar, 2011;Boano et al., 2014;Cheong et al., 2007;Gooseff et al., 2007;Liu et al., 2020;Trevisan & Periáñez, 2016). However, the nature of transient-storage mechanisms is not always clear. While the hyporheic-zone community tends to attribute transient storage mainly to the exchange between the stream and the pore water in sediments, it has been recognized that transient-storage zones also exist in the surface water itself, and separating hyporheic from surface transient storage has been identified as key problem for reactive transport studies, because the type and intensity of biogeochemical reactions in these zones differ significantly (Guillet et al., 2019;Liu et al., 2020;Riml et al., 2013;Robinson et al., 2007). Due to different mechanisms, the residence time distributions in surface and hyporheic transient-storage zones typically differ significantly with surface transient storage zones undergoing more rapid exchange, with the main flow compared to hyporheic transient-storage zones. However, there is a significant overlap, particularly when considering short hyporheic flow-paths through ripples. The current study aims at relating surface-transient-storage parameters to easy-to-measure properties of the stream so that additional tailing in breakthrough curves of stream tracers can be assigned to hyporheic exchange.
Surface transient storage is sensitive to the physical conditions and properties of the stream, including bed topography, discharge, channel structure, and channel roughness, which are largely affected by boulders, vegetation, precipitation, among others (Zarnetske et al., 2007). Previous investigations have extensively explored the mass exchange between surface-transient-storage zones and the main flow from a hydrodynamic perspective, in which the surface-transient-storage zones are characterized by flow fields of water recirculation . Flow structures such as vortex interactions and large-scale velocity fluctuations generally regulate the exchange process to a large extent. It has been postulated that the residence time within surface-transient-storage zones can be approximated either by an exponential distribution in small streams (Gooseff et al., 2005), which is mathematically equivalent to assuming a well-mixed transient-storage zone undergoing first-order exchange with the main stream, or by a power-law distribution in large rivers (Anderson & Phanikumar, 2011), which can be achieved, e.g., by particular variants of the Multi-Rate Mass Transfer (MRMT) parameterization (Haggerty et al., 2000;Haggerty & Reeves, 2002). The mean residence time within surface transient-storage zones has been related to the structure of the surface-transient-storage zone, discharge, gyre distribution, and others.
Classical one-dimensional river-hydraulic models generally assume a shear-flow profile over the entire width of the stream (Fischer, 1973;Fischer et al., 1979). This assumption breaks down in streams containing many rock blocks or other macro-roughness elements that are not fully submerged at low flows (Boyero & Bosch, 2004;Gooseff et al., 2008;Liu et al., 2020;Niayifar et al., 2018). Under such conditions, surface transient storage can be explained by persistent recirculation cells behind protruding roughness elements. At ultralow flow, the stream starts to resemble a porous medium, in which the major variability in flow velocity is caused by the arrangement of smaller and large openings between protruding roughness elements. Under such conditions, solute particles can sample the autocorrelated velocity variations by advection only, resulting in hydrodynamic dispersion with asymptotic dispersion coefficients scaling linearly with the mean velocity, instead of quadratically as in pure shear flows (Bear, 2013). With the increasing probability of low-flow events due to climate change (Dai, 2013;Sheffield & Wood, 2008), rough bathymetry exerts greater influence on river transport and transient storage than before. At low flows, many natural or restored rivers with rough bathymetry consist more or less of connected pools and riffles, leading to some potential surface-transient-storage zones surrounded by high bed elements that inherently affect nutrient cycling and stream connectivity (Sandoval et al., 2019). At sufficiently high discharge, the bed topography will be totally submerged under the water surface, diminishing the role of surface transient storage in solute transport (Ward et al., 2013;Zarnetske et al., 2007), even though preasymptotic solute spreading in shear flow may also be parameterized by transient-storage models (e.g., Reichert & Wanner, 1991;Schmalle & Rehmann, 2014).
Several studies have analyzed the effects of flow recirculation and still-water zones on macroscopic stream transport using the one-dimensional mobile-immobile transport model (Anderson & Phanikumar, 2011;Gooseff et al., 2008;Jackson, Haggerty, Apte & O'Connor, 2013McCoy et al., 2008;Sandoval et al., 2019), typically fitting the model to tracer-test data. Tracer-derived estimates of these coefficients have the disadvantage of not being transferable to the same reach at different discharge, not to mention other stream reaches. As a consequence, it has been difficult to predict how solute transport coefficients change with changing discharge or topography in the stream of interest. Some studies tried to elucidate the dependence of solute transport coefficients on riverbed topography. For example, Jackson, Haggerty, Apte, & O'Connor (2013) proposed an empirical relationship between hydraulic and geomorphic parameters to the mean residence time distribution of a single surface-transient-storage zone. Femeena et al. (2019) performed a meta-analysis of previous tracer studies proposing concise laws to predict transient-storage parameters to river characteristics. Through laboratory-scale experiments and numerical simulations, Aubeneau et al. (2015) related the residence-time distribution to the fractality of the rough-bed topography obtained from laser scanning at a very high resolution, followed by Lee et al. (2020) in a numerical work exploring the dependence of both interfacial fluxes and hyporheic travel times on the fractal properties. Some studies have explored the surface transient storage caused by structured cavities or groynes (McCoy et al., 2008;Jackson, Haggerty, Apte, & O'Connor, 2013. To the best of our knowledge, however, no work has specifically related surface-transient-storage parameters to stream-bed topography or discharge in channels with rough bed under low-flow conditions. In this study, we focus on solute transport in a synthetic channel with rough bathymetry under low flow conditions. The main goal is to quantitatively relate surface-transient-storage parameters to hydraulic and morphological characteristics. In particular, we ask how dispersion scales with velocity at very low flows, and which easy-to-determine properties of a stream can explain the transport parameters. The study is motivated by previous research performed at the Steinlach River in Tübingen, Germany (Guillet et al., 2019;Liu et al., 2020;Schwientek & Selle, 2016). Upon river restoration, the vertical hard basin drop structures of the mostly straight river channel were replaced by grouted sloping boulders, and individual rock blocks were placed between these drops. The natural sediments are dominated by medium sized gravel. It is difficult to precisely measure the rough bathymetry of this stream including details of riffles and pools at the scale of a river reach. As an alternative, we generate a synthetic channel mimicking a river in which the bed is predominantly made of macro-roughness elements and study flow and transport at different flow rates and for different bed slopes in this synthetic river. Toward this end, we use a two-dimensional numerical flow model, capable of simulating recirculation zones, and an associated two-dimensional solute transport model, and characterize the resulting flow field and transport characteristics. Subsequently, we fit the coefficients of a standard one-dimensional mobile-immobile transport model to the numerical breakthrough curves of the two-dimensional model. We then search for relationships between the coefficients of the mobile-immobile transport model and easy-to-measure properties of the river, namely, the discharge, the channel slope, and the ratio of the surface area of protruding rocks to that of open water. This analysis provides an easy-to-use means to evaluate the transport coefficients in streams with rough bathymetry under low-flow conditions. While the analyses are limited to an idealized synthetic channel that only mimics a real stream with similar characteristics, the results give insight into the mechanism of surface transient storage and the dependence of one-dimensional mobile-immobile transport models on river characteristics.

Two-Dimensional Flow and Transport Equations
Depth-integrated flow without volumetric source/sink terms is simulated by the two-dimensional shallow-water equations, consisting of the continuity and momentum equations: Accounting for turbulent momentum transfer is crucial, particularly in the direction transverse to the local flow direction. Without eddy viscosity, recirculation zones could not form. In particular, the popular diffusive-wave approximation of the two-dimensional shallow-water equations leads to flow fields without closed streamlines. For the eddy viscosity ν t , we use the following approximation: in which D T is a dimensionless empirical constant.
Solute transport in the two-dimensional model is simulated by the advection-dispersion equation: in which c [ML −3 ] is the solute concentration and D t = ν t the turbulent diffusion coefficient, being identical to the turbulent viscosity. Equation 4 is subject to no-flux boundary conditions along the banks and at boundaries to protruding bed elements, to a zero-dispersive-flux condition at the downstream end of the reach, a known flux condition at the upstream end of the reach, and an initial condition of zero throughout the domain.

One-Dimensional Mobile-Immobile Transport Model
We fit the one-dimensional mobile-immobile transport model to the flux-weighted cross-sectional average of the concentration time series:  Table A1).
The fitted model assumes the instantaneous injection of a solute mass m [M] injected into the flux at location x = 0: where A m is the cross-sectional area of the mobile domain, and the Dirac delta-function δ has the inverse dimension of its argument. As second auxiliary condition, we assume that the concentration drops at least exponentially in the far-distance limit: In the initial state, both domains are assumed to be solute-free: ( 0) = 0 (9) ( 0) = 0. (10)

Numerical Implementation and Description of the Test Case
We generated a synthetic channel with random, spatially periodic, auto-correlated bed topography using the spectral random-field generating method of Dietrich and Newsam (1997), amended by a longitudinal slope and slopes at the river banks. The synthetic-channel generation is implemented in Matlab. Figure 1 shows the resulting bathymetry of a unit cell with 100 m length and 5 m width with a spatial resolution of Δx = Δy = 0.1 m. The auto-covariance function of the bed fluctuations is an isotropic squared-exponential function with a correlation length of 0.5 m and a standard deviation of 0.2 m. The two banks have a slope of 45°, and the three mean channel slopes tested are 1, 2, and 3‰, respectively. Because the unit cell is periodic, longer continuous fields can be generated by adding several unit cells with seamless transition. The synthetic channel generated in Matlab was saved as a geo-referenced tif-file. The corresponding Matlab code is made available in a public repository (see the acknowledgments).
We simulated steady-state two-dimensional flow for different combinations of channel slope and discharge with the 2-D flow model of HEC-RAS 5.0.7 (Brunner, 2016), which uses a cell-centered Finite-Volume scheme. The Matlab-generated tif-file was imported into HEC-RAS. The structured grid used in HEC-RAS coincided with the Matlab-generated bathymetry field. To achieve steady state, we simulated transient flow until the head field did not change anymore. The Manning's roughness coefficient of the bed and bank surface is set to 0.06 s/m 1/3 . The dimensionless coefficient D T to compute the eddy viscosity by Equation 3 is set to 0.11. We considered 11 constant discharges from 0.005 to 1.0 m 3 /s (0.005, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, and 1.0 m 3 /s). To achieve a periodic flow field, we simulated the shallow-water equations for a domain made of eight unit cells, applying a constant discharge boundary condition at the inlet and local normal-flow depth at the outlet of the extended domain. We chose the full shallow-water equations rather than the diffusive-wave approximation as the computational method in the transient-flow analysis, because the latter cannot simulate recirculation zones. Since the inlet and outlet boundary conditions affect the flow field up to a substantial distance to the boundaries, we extracted the flow field of a single unit cell from the interior of the domain, where the water depth and velocity field over the inlet and outlet boundaries are almost the same. Under these conditions, the transition of the velocity field from the outlet to the inlet boundary of the unit cell is seamless. Within the unit cell, however, the rough bathymetry causes strong spatial variations in flow velocity and water depth.
The water depths in all cells and the velocities at the interfaces of the 2-D HEC-RAS simulations were stored as HDF5-file and imported into Matlab for postprocessing and transport simulations. The corresponding Matlab codes are made available in the same reporsitory as the other codes. As the steady-state depth-integrated velocities are non-divergent, we could compute a field of the stream-function. Contour lines of the stream function are streamlines. Streamtubes bounded by streamlines connecting the inlet and outlet boundaries belong to the main channel flow, whereas recirculation zones exhibit closed streamlines. We identified all recirculation cells by extracting all closed streamlines and excluding closed streamlines that lay within other closed streamlines. The remaining outer closed streamline of each recirculation cell defines its position, shape, horizontal surface area and water volume, which will be displayed in the following section.
In the transport simulations, we added eight unit cells of the periodic velocity field, resulting in a transport domain of 800 × 5 m. We chose this value to ensure that solute dispersion reaches its asymptotic stage. We simulated two-dimensional advective-dispersive transport using the cell-centered Finite Volume Method with upwind differentiation of advection for spatial discretization and implicit Euler integration in time with a step-size of 0.1Δ ∕̄ , in which is the mean velocity in the longitudinal direction. The constant inlet concentration was set to unity. We then computed flux-weighted breakthrough curves at cross-sections with different distances from the inlet boundary. With an inflow concentration of unity, the resulting breakthrough curves can be interpreted as The one-dimensional mobile-immobile transport model was solved in the Laplace domain and numerically back-transformed into the time domain by the method of de Hoog et al. (1982), implemented in Matlab. This model was fit to the breakthrough curves at all observation planes by minimizing the sum of squared residuals using the Trust-Region Reflective Least-Square method of the function lsqnonlin of the optimization toolbox of Matlab (Coleman & Li, 1996). Subsequently, we applied linearized uncertainty propagation to obtain the covariance matrix of the parameter uncertainties: in which θ is the n θ × 1 vector of parameters, θ LS is the associated least-square fit, p 2D (τ i ) and p MIM (τ i , θ) are the probability densities of travel time τ i obtained by the two-dimensional model and the one-dimensional mobile-immobile model, respectively, J is the n τ × n θ Jacobian matrix of partial derivatives derived about the least-square fit, n θ is the number of parameters, and n τ the number of discrete travel times at which the comparison of the two models is made. The resulting covariance matrix C θθ contains the estimation variances on the main diagonal. The square-root of the estimation variances are the corresponding standard deviations of estimation : Figure 2 provides the visual impression that with increasing discharge the surface area A r of recirculation cells decreases, whereas the area A w of the mobile water increases. Figure 3a shows this dependence in a quantitative way, revealing that at ultralow flow the surface area of the recirculation cells is smaller than at very low flows. Figure 3a also includes the surface-area fraction A r /A w of recirculation cells with respect to the main flow as function of discharge. Figure 3b illustrates the corresponding dependencies on discharge for the water volumes rather than surface areas. Again, we see a nonmonotonic behavior of V r /V w with discharge at ultralow flows. The size of the recirculation cells depends on the intricate interplay of lateral shear, causing rotation, and bottom friction, suppressing the latter, and the placement of flow-blocking bathymetric high elements.

Two-Dimensional Flow Fields and Concentration Distributions
Previous studies on surface and hyporheic transient storage postulated that the cross-sectional area ratio of the immobile to the mobile domain β (which is identical to the volume ratio) monotonically decreases with increasing discharge (Zarnetske et al., 2007). In our setup, this is confirmed only once a certain threshold discharge, and thus mean water depth, has been exceeded. Indeed, at relatively larger flow, a considerable part of recirculation cells merge into the main flow, reducing both the absolute and relative area and volume of the recirculation cells. At sufficiently large flows, no recirculation cells can be observed any more.
With the given bathymetry, the maximum volume fraction of recirculation cells ranges between 20% and 30%. Larger channel slopes result in both larger areas and volumes of recirculation cells and in smaller surface areas and volumes of total open water. The smaller water depth in steeper channels leads to a reduction of the surface area and volume of the open water. However, the higher velocity corresponds to larger kinetic energy, which strengthens the blockage effect of protruding bed elements, which explains the larger surface areas of the recirculation cells overcompensating the smaller water depths when considering the volume of the recirculation cells. velocities near the banks are reduced by the no-slip condition, resulting in a macroscopic shear-flow profile causing Taylor-Aris dispersion (Fischer et al., 1979).

Parameter Estimation for the One-Dimensional Mobile-Immobile Transport Model
There are two different approaches of obtaining the coefficients of the mobile-immobile transport model. In the first approach, we assign the delineated recirculation cells of the flow simulations to the immobile zones, directly giving us the value of the immobile-to-mobile volume ratio β. In this approach, we can compute the mean advective velocity v in the mobile domains from the discharge, the length of the reach, and the mobile water volume: in which L [L] is the length of the channel, and V m [L 3 ] is the total volume of mobile water in the 2-D numerical model. The remaining coefficients D and k need to be fitted from the breakthrough curves of the two-dimensional model. While this approach associates the immobile domain with a physically meaningful entity, it requires knowing the size of the recirculation zones, which is not realistic for natural streams. In Section 3.3.3, we will discuss how to estimate β from easier-to-measure properties of the river.
In the second approach, which is adopted in the present work, all coefficients are fitted. We used the Trust-Region Reflective Least-Square method implemented in the function lsqnonlin of Matlab's optimization toolbox (Coleman & Li, 1996) to jointly fit several breakthrough curves and estimate uniform values of D, k, β and v. The mobile-immobile model could be fitted for discharge values in the range of 0.005 m 3 /s to 0.7 m 3 /s. At discharge values larger than 0.7 m 3 /s, the size of recirculation cells becomes negligible, and solute transport simplifies to advection and dispersion without transient storage effects. Figure 5 compares the travel-time distributions derived from the breakthrough curves of the two-dimensional model to the travel-time distributions of the fitted one-dimensional mobile-immobile transport model. The agreement is excellent with some differences at the shortest travel distances, where the solute plume has not yet sampled the velocity field in a representative way. The standard deviations of estimation for all parameters, computed by Equations 11-13 are all orders of magnitude smaller than the corresponding fitted parameters, indicating that not only the agreement between the fitted one-dimensional mobile-immobile model is good, but the parameters are also identifiable, when jointly fitting the travel-time distributions observed at several travel distances. Table   Figure 3. Surface area and volume of recirculation cells, mobile water and the ratio between the two for the entire channel length of 800 m: blue circle: S = 1‰; red square: S = 2‰; green diamond: S = 3‰. S1 in Supporting Information S1 lists the fitted parameters of the one-dimensional mobile-immobile transport model and some parameters of the two-dimensional numerical model for all slopes and discharge values. In the following, we relate the fitted parameters to characteristics of the rivers.
As mentioned above, four coefficients of D, k, β and v are fitted through the 1-D mobile-immobile transport model, of which β and v are compared with those that can be obtained by postprocessing the flow field of the 2-D numerical model. Figure 6 shows good agreement with RMSE-values of 0.018 for the β and 0.002 m/s for the velocity, especially the excellent agreement of velocities. The fitted coefficients are physically representative of the corresponding real values. It implies that the four coefficients obtained from the two different approaches should be almost identical.

Dependence of Coefficients on Easy-to-Measure Quantities
A key difficulty in transient-storage modeling is that the coefficients are typically obtained by fitting the model to measured breakthrough curves, thus requiring the results of a solute-transport experiment to predict solute transport. In the following, we relate the fitted parameters of the one-dimensional mobile-immobile transport to easy-to-measure quantities of the stream, namely the discharge, slope, and the ratio of the surface area associated to protruding bed elements to that of the open water body.

Dependence of Water Depth and Velocity on Discharge and Channel Slope
The mean water depth h of the open water generally increases with discharge while it decreases with channel slope. The h(Q) relationship shown in Figure 7a can be fitted by the following number-value equation: requiring that the water depth is given in meters, and the discharge in cubic meters per second. With an RMSE value of 0.004 m, the fitted dependence of water depth h on discharge Q and channel slope S is robust. Considering the rough bed topography, the water remains in unconnected pools allowing no flow as long as the water depth is below the percolation threshold of 0.206 m. Field studies on Q(h) relationships in streams with rough bed (e.g., Ferguson et al., 2019) often lack such percolation thresholds. If the streambed is permeable, there can be flow through water pools that are not connected on the surface. For the water depth above that threshold, we could fit a power-law dependence. We assume that the form of Equation 15 should hold in most natural streams with rough beds, even though the exact coefficients will be different. The open markers in Figure 7a represent the average water depth h m measured only in the main flow. Most recirculation cells appear in pools behind protruding bed elements, so that the mean water depth in the recirculation cells is somewhat larger than that of the main flow, leading to h m being slightly smaller than h. This effect is enhanced at low discharge at which the pools store a considerable fraction of the water.
Note that applying Manning's friction law under normal-flow conditions with the simplifying assumption that the hydraulic radius equals the mean water depth would lead to the dependence of h ∝ S 0.3 Q 0.4 , which differs considerably from the fitted power-law exponents of 0.226 for the slope and 0.58 for discharge. This deviation may be explained by the mean hydraulic radius depending nonlinearly on water depth at low flow for the given complex bathymetry. In fact, the water depth in the present case is on the same order as the height of the rough bed elements, making the effective Manning's coefficient n increase rapidly with decreasing water depth. Hence Manning's formula with a discharge-and stage-independent n-value is not applicable to predict the velocity in such streams. The open markers in Figure 7b represent the results of an empirical equation derived by Froehlich (2012) for the velocity v emp in channels with a rough bed:  [-] is the Darcy-Weisbach friction factor, f b [-] is the bed friction factor, f w [-] is the bank friction factor (Cheng & Chua, 2005), w [L] is the mean width of the river at the surface, Re is the Reynolds number of the main flow, and ν = 1.004 × 10 −6 m 2 /s is the kinematic viscosity of water. Similar to Froehlich (2012), we approximate the roughness height k s as 1.7 times half the bed topography length scale of 0.25, resulting in k s = 1.7 × 0.25 m = 0.425 m. By fitting Equation 17 to the velocity of the one-dimensional mobile-immobile transport model, the dimensionless coefficient B is found to be 1.432. The fit is excellent with an RMSE-value of 0.0057 m.

Dependencies of the Dispersion Coefficient and the Mobile-Immobile Exchange Rate Coefficient on Velocity, Discharge, and Channel Slope
For solute transport in shear flow, the asymptotic dispersion coefficient D is proportional to v 2 (Aris, 1956;Taylor, 1953). It is also widely accepted that D scales approximately quadratically with v in natural rivers (Deng et al., 2001;Fischer, 1973;Fischer et al., 1979;Seo & Cheong, 1998). However, the empirical data to confirm the latter relationship were collected in streams in which the water depth is considerably larger than the height of the roughness elements. This is different in our setting. Figure 8 which implies a "geometric" dispersion process that has been well documented for Darcy flow in porous media (Bear, 2013). In fact, the visual impression of the flow field given by Figures 2 and 4 makes the river look like a high-porosity porous medium, in which the main flow at low discharge is quite tortuous, split up into different flow paths that circumvent protruding bed elements and recirculation cells. These flow paths differ in width, depth, and velocity, and show no persistent connectivity of high-and low-velocity sections. Solute particles introduced into the main flow can thus sample high and low velocities by advection alone with some temporal auto-correlation of velocity along the flow path. Such a transport regime results in a longitudinal dispersion coefficient D within the main flow that scales linearly with the mean velocity in the main flow. For the given geometry, the corresponding dispersivity α = D/v [L] was 0.869 m, regardless of the channel slope.
Besides advection and dispersion in the main flow, the solutes undergo exchange with the recirculation cells by turbulent mixing, where they are intermittently trapped before being released back to the main flow. Figure 9 Higher velocities at larger discharge enhance turbulence and thus accelerate the mass transfer between recirculation zones and the main flow. With higher exchange rates, the mean residence time in the transient-storage zone decreases accordingly. Comparing the two relationships shown in Figures 8 and 9, the exchange rate coefficient is k ≈ 0.096 [m −2 ] × D.
Combining the empirical relationships v(Q), D(v), and k(v), it is straightforward to derive the dependence of the dispersion coefficient D in the mobile domain and the exchange-rate coefficient k of the mobile-immobile transport model on the discharge Q and the channel slope S. Using the units as before, the corresponding number-value equations are Figure 10 shows that the resulting relationships are acceptable, with RMSE values of 0.013 m 2 /s for the dispersion coefficient and 0.003/s for the exchange-rate coefficient. At this point, we have related the mean velocity v and dispersion coefficient D in the mobile domain, and the exchange-rate coefficient k to discharge and channel slope, in which the power-law exponent of discharge is considerably higher than that of the slope. The only remaining coefficient of the mobile-immobile transport model is the ratio β of the cross-sectional areas in the immobile and mobile domain.

Estimating the Immobile-to-Mobile Storage Ratio β From Easy-to-Measure Properties
In this section, we propose to relate the storage ratio β (, i.e., the volume in the immobile domain over the volume in the mobile domain) to the ratio of the horizontal area of all bed elements ("stones") in the stream that protrude the water surface divided by the horizontal area of the remaining water surface. This metric could be easily evaluated at a river by image analysis of aerial photos. We denote this dimensionless metric the area ratio of stone to water, i.e., A sw [-]. Figure 11 shows the relationship of the area ratio A sw on the discharge Q and the channel slope S for the given bathymetry with an RMSE value of 0.032. Using the same units as before, the fitted number-value equation is We see an exponential decrease of A sw with increasing discharge and thus water depth. At sufficiently high discharge, all bed elements are submerged, and A sw becomes zero. At ultralow discharge, some areas with low bed surface surrounded by high-bed elements may be completely disconnected from the main flow, so that they remain dry in the simulation and contribute to the protruding elements in the area calculations of high bed elements and water surface. This may explain the deviation from the exponential A sw (Q)-relationship at ultralow discharge. If transient-storage zones can be associated to recirculation zones, and if recirculation zones occur in the wake of protruding macro-roughness elements, the stone-to-water area ratio A sw should be a good and comparably easyto-assess indicator of the immobile-to-mobile storage ratio β.
In the discussion of the two-dimensional numerical model, we have observed that both the area and volume ratio of the recirculation cells to open water, and thus the real β-value in the 2-D numerical model, first increases and then decreases with increasing discharge (see Figure 3). This is different to the other coefficients of the mobile-immobile transport model, i.e., v, D, and k, which all increase monotonically with discharge. Most recirculation cells appear behind the protruding stones; they are formed under the effect of lateral shear stress induced by the flow velocity difference between the slowly rotating recirculation cell and the fast moving main flow. Larger portions of the rough bed protrude the water surface at lower discharge, thus enhancing the blocking effect which should lead to larger recirculation cells, the faster moving main flow at larger discharge exerts larger inertia which elongates the recirculation cells. The resulting nonmonotonic dependence on discharge inspired us to make β dependent on the product A sw × v rather than the stone-to-water area ratio A sw alone. Figure 12a shows that a power-law relationship between the fitted β-values of the mobile-immobile transport model and A sw × v yields a good agreement with an RMSE value of 0.035. Assuming the same units as before, the fitted number-value equation is Combining the relationships given in Equations 16, 28 and 29, we can now derive an empirical dependence of the immobile-to-mobile storage ratio β on discharge Q and channel slope S for the given bathymetry, again expressed as number-value equations assuming SI units:  Figure 12a is somewhat more reliable (RMSE = 0.035). Previous field studies indicated a monotonous decrease of β with increasing discharge, most likely because their measurements did not include the very-low-discharge regime (Zarnetske et al., 2007).

Discussion
With global warming, droughts will become more frequent (Dai, 2013;Sheffield & Wood, 2008). Under the resulting low-flow conditions, the bed topography of the rivers will become more significant for surface flow and longitudinal solute spreading in streams. At very low discharge, the main flow splits into various paths circumventing the high-bed elements and water recirculation zones, violating the assumption of regular shear-flow underlying the Taylor-Aris dispersion model. The numerous water recirculation zones trap solutes and result in long tails of the solute breakthrough curves. We consider the one-dimensional mobile-immobile transport model adequate to describe macroscopic solute transport at low flows, because the standard advection-dispersion equation misses the tailing whereas more complex parameterizations are difficult to calibrate, not to speak about predicting their parameters. The present study aimed at relating the coefficients of the transient-storage model to easy-accessible metrics of rivers with rough bathymetry. Working with a specific synthetic bathymetry, the exact numbers of these relationships are less important than the functional shapes.
Numerous preceding studies have developed relationships among mean flow velocity, discharge, channel slope, water surface width and other stream parameters (e.g., Carrillo et al., 2021;Ferguson et al., 2019). For channel with rough bed, like pool-riffle channels, some field works indicated that the exponent describing the mean velocity with discharge is between 0.48 and 0.6 (Comiti et al., 2007), of which our numerical work was also within the range. While it has been established that discharge acts as the first-order control on hydraulic geometry (Comiti et al., 2007), only general trends of the dependence on bed topography were given for channels with rough bed because data on the bathymetry was lacking. Kang and Sotiropoulos (2011) numerically explored the detailed flow hydrodynamics in a short meandering channel with rough bed and determined the recirculation zones induced by the boundary, but not at low discharge at which the bed macroroughness exerts considerable influence on flow. As expected, downstream of most protruding high-bed elements, two parallel reversely rotating recirculation zones develop (Figure 2). Our random bed topography results in recirculation zones of various sizes and shapes, which are less predictable than recirculation zones in structured cavities or after regular obstacles (McCoy et al., 2008;Jackson, Haggerty, Apte & O'Connor, 2013;Jackson et al., 2015).
Our work specifically focused on surface transient storage in streams at very low flow at the reach scale, considering the average effects of many protruding elements and recirculation zones. A key physical insight was that we could explain coefficients of the one-dimensional model by the product of the area ratio of protruding elements with the velocity. While Gooseff et al. (2005) applied an exponential surface-transient-storage residence time distribution and related it to the averaged pool geometry in a bedrock stream, we physically integrated the hydraulic and morphological metrics of the stream into transient-storage metrics. Zarnetske et al. (2007) reported that the immobile-to-mobile storage ratio β linearly decreases with increasing discharge Q, because the dominant storage zones were assimilated into the main flow with larger discharge. We see the same, but not at ultralow flows, at which an increase in the water level causes a more rapid expansion of dead-end recirculation zones than that in the main channel. The maximum value of β in our simulations was about 30% (Figure 12b), under very low-flow conditions. Field surveys revealing higher β-values must be affected by hyporheic exchange. In addition, we could relate the water depth, mean velocity, dispersion coefficient and mass-transfer coefficient to the discharge and slope, which will be useful to predict the solute transport over different flow regimes and thus constrain the associated biogeochemical processes such as photodegradation. While the exact coefficients of the relationships depend on the details of the virtual bathymetry, we believe that we can point river engineers and scientists to easy-to-measure properties of the stream that control surface transient storage.
The dominant mobile-immobile exchange mechanism is the turbulent vortex between the main flow and the recirculation zones. Higher velocities at larger discharge enhance turbulence and thus accelerate the mass transfer between recirculation zones and the main flow. In general, the strength of the turbulent vortex should be positively correlated with the mean velocity. However, the flow in channels with different discharge and slope results in recirculation zones of different shapes and sizes, the velocity of the main flow dragging the recirculation zones varies considerably, so that the strength of the turbulent vortex does not depend that linearly on the mean velocity, and the exchange rate coefficient k in Figures 9 and 10b is a bit scattered.
As with any synthetic study, one may question how realistic the virtual reality is and thus how solid the determined relationships are, but the question of generality would have also arisen if we had worked with a detailed survey of the bathymetry in a specific river -without the computational convenience of periodicity. The synthetic rough channel fundamentally captured the high and low bed elements that are seemingly randomly distributed in real streams (Gooseff et al., 2008;Liu et al., 2020;Thompson, 2013). In real streams, fluvial geomorphology is affected by flow, sediments, obstacles, and other influences. For example, the flow shear may erode the bed and form scour holes in the vicinity of obstacles, while sediments may accumulate in low-flow zones caused by the same or other obstacles. Therefore, the real rough bed topography of a river will be somewhat different from our synthetic test case, especially in streams with a high density of sediments (Kleinhans & van den Berg, 2011;Inoue et al., 2016). Another difference between the synthetic channel and natural streams is that we deliberately excluded hyporheic exchange by considering an impermeable bed. This was needed to separate surface transient storage from hyporheic transient storage, and is justified by the small influence of hyporheic fluxes on surface hydrodynamics. However, in a follow-up study we would like to analyze the individual versus combined effects of surface and subsurface transient storage using the same bathymetry.
Our parameterization has been developed for very low flows. At higher discharge, when the water level is significantly above the largest bed elements, solute transport will revert to Taylor-Aris dispersion with a dispersion coefficient scaling with the mean velocity squared. We have not yet explored at which point the transition between the quasi-porous-media and classical shear-flow behavior occurs. Most likely, this would require a larger model domain to achieve the asymptotic behavior.

Conclusions
Rough river beds are common in natural streams and restored rivers. Climate change will lead to a reduction of stream flows, thus making the flow scenario discussed in this study more likely. Because of the large bed-surface-area to water-volume ratio, biogeochmical transformations and pollutant turnover are expected to be accelerated in streams with rough bathymetry at low flows. Enhanced hyporheic exchange is also expected, but as we show surface transient storage cannot be neglected and should be estimated separately from hyporheic exchange. In order to explore the surface transient storage emerging under such conditions, we studied two-dimensional flow and transport in a synthetic channel with rough bed topography under low-flow conditions and fitted the one-dimensional mobile-immobile transport model to the breakthrough curves for different combinations of discharge and channel slope. The key goal was to establish relationships between easy-to-measure properties of the stream and the fitted coefficients of the mobile-immobile transport model. The principal findings of this research point to the following: 1. Streams with many macroroughness elements resemble porous media when the water depth is below the height of these elements. The water flows through gaps between the protruding bed elements experiencing large velocity contrasts in seemingly random auto-correlated patterns. Recirculation cells generally form behind the protruding bed elements caused by the shear stress between nearly stagnant water within the recirculation cells and the fast-moving flow in the main stream. The recirculation zones act as surface transient-storage zones that can intermediately trap solutes. Hence, the macroscopic description of solute transport by the one-dimensional mobile-immobile transport equation is adequate. The simplest version of this model, assuming first-order exchange between a well-mixed mobile domain and an equally well-mixed immobile domain with constant coefficients appears adequate if the spatial statistics of the bathymetry variations is stationary 2. Considering an impermeable riverbed, there is a percolation threshold of the water depth in our model. Below a minimum value, no flow occurs because the open water is not connected. With a permeable bed, hyporheic flow could connect otherwise disconnected pools thus allowing non-zero discharge even below the percolation threshold. Above this threshold, the water depth undergoes a power-law dependence on discharge and slope. Also the velocity in the mobile domain undergoes a power-law dependence on discharge and slope, but the corresponding power-law exponents differ from what is expected for normal flow in rivers with large water depths 3. The dispersion coefficient in the mobile domain of the mobile-immobile transport model scales linearly with velocity. This is the classical behavior of hydromechanic dispersion in porous media and caused by the same mechanism: Solute particles can experience the variability of the auto-correlated velocity field by advection alone, establishing a "geometric" dispersion process 4. The rate coefficient for mass transfer between the mobile and immobile domains also scales linearly with velocity, though shows a bit scatter 5. At ultralow flow, the immobile-to-mobile storage ratio β increases with increasing discharge. The smaller the discharge, the more bed elements protrude the water surface giving the chance to produce recirculation cells. Conversely, low flow velocities imply less drag by lateral transfer of momentum. At ultralow flow, the latter effect dominates. From a still very low value of discharge on, β decreases with increasing flow, which is first mainly attributed to an increase in the volume of the main-flow domain. The larger the water depth, the fewer bed elements protrude the water surface until all are submerged. At sufficiently high discharge, no stagnant recirculation zones remain 6. We propose to measure the surface area of the protruding bed elements and of open water. The ratio of these two surface areas, times the velocity, exerts a power-law control on the immobile-to-mobile storage ratio β 7. With the given bed topography, we could relate all coefficients of the mobile-immobile transport equation to the slope, discharge, and stone-to-water area ratio In this work, we have analyzed a synthetic river with rough bed topography, in which the variability of the bathymetry is a stationary random space variable. Setups in which grouted sloping boulder drops and boulder-free stream stretches alternatively might behave slightly different. Thus, further field confirmations of our qualitative findings are necessary. Also it is worthwhile to study the dependence of the empirical coefficients on geometric properties of the bathymetry.

Symbol Dimension Description
Mobile-immobile transport model coefficients