Benthic Uptake Rate Due to Hyporheic Exchange: the Effects of Streambed Morphology for Constant and Sinusoidally Varying Nutrient Loads

Hyporheic exchange carries reactive solutes, which may include biological oxygen demand (BOD), dissolved oxygen (DO) and reactive dissolved inorganic nitrogen (Nr), into the sediment, where biochemical reactions consume DO. Here, we study the impact of streambed morphology, stream-reactive solute loads and their diel oscillations on the DO benthic uptake rate (BUR) due to hyporheic processes. Our model solves the hyporheic flow field and the solute transport equations analytically, within a Lagrangian framework, considering advection, longitudinal diffusion and reactions modeled as first order kinetics. The application of the model to DO field measurements over a gravel bar-pool sequence shows a good match with measured DO concentrations with an overall agreement of 58% and a kappa index of 0.46. We apply the model to investigate the effects of daily constant and sinusoidally time varying stream BOD, DO and Nr loads and of the morphodynamic parameters on BUR. Our modeling results show that BUR varies as a function of bedform size and of nutrient loads and that the hyporheic zone may consume up to 0.06% of the stream DO at the pool-riffle bedform scale. Daily oscillations of stream BOD and DO loads have small effects on BUR, but may have an important influence on local hyporheic processes and organisms' distribution.


Introduction
Human activities have recently increased nutrient loads in aquatic environments, which would typically be nutrient-limited [1].This increase has led to water body eutrophication, thereby threatening the integrity of aquatic ecosystems [1][2][3][4][5].In freshwater, dissolved inorganic reactive nitrogen, Nr, primarily under the form of nitrates, NO3 − and ammonium NH4 + and dissolved organic matter have cascading effects on dissolved oxygen (DO) concentrations [6,7].For example, nitrogenous biological oxygen demand (NBOD = L N ) is associated with the transformation of NH4 + to NO3 − , whereas respiration and metabolic activities are characterized by the carbonaceous biological oxygen demand (CBOD = L C ). High CBOD is considered the most increasing water pollution problem in the world [8].
These quantities show intertwined complex dynamics over a hierarchy of scales, from daily to seasonal [9,10].Sinusoidal oscillations have been observed to approximate the daily variations well for DO [10], temperature [11,12] and BOD [13].Daily oscillations of DO are due to photosynthesis/respiration plant cycles, those of temperature to the day/night cycle, while BOD variations are chiefly correlated with discharge from wastewater treatment plants.
Following the previous works of Rutherford et al. [16], we propose that this reduction in DO concentrations between downwelling and upwelling waters could be used to define the oxygen consumption due to biogeochemical activities occurring within the sediment.This oxygen consumption is also referred to as the benthic uptake rate, BUR.
Although field experiments show the importance of hyporheic residence time in riverine-hyporheic biogeochemistry [22,[31][32][33], models coupling surface-subsurface exchange fluxes with biogeochemical reactions to quantify BUR are limited [16], with more publications addressing nitrate cycling (e.g., [31,[33][34][35][36][37]). Rutherford et al. [16] analyzed BUR for hyporheic exchange induced by two-dimensional dune-like bedforms, while neglecting solute diffusion within the streambed and nitrification processes.Boano et al. [38] numerically analyzed biogeochemical zonation of intra-meander bars due to differences in the residence time of hyporheic water, but they did not evaluate stream oxygen consumption due to this process, nor the influence of stream bedforms on BUR.In a similar fashion, Marzadri et al. [31,35] proposed a simplified, yet effective, model of the nitrogen cycle within the hyporheic zone, but neglecting BOD loads.Similarly, other recent investigations focused on nitrate cycling and modeled dissolved oxygen transport to quantify hyporheic zone redox conditions (e.g., [34,36,37]).More recently, we investigated the effect of streambed morphology on the hyporheic thermal regime and studied the effects of daily temperature variations of stream water on dissolved oxygen distribution within the hyporheic zone [36].In particular, we showed that mean daily temperatures suffice to characterize the hyporheic zone biogeochemistry at the bedform scale, while daily temperature oscillations may affect local biogeochemistry in sediment volumes with short hyporheic residence times.However, we did not quantify the effects of streambed size nor of BOD loads on BUR.
In the present work, we close this gap by extending the Rutherford et al. [16] model to gravel-bed rivers with a three-dimensional alternate bar morphology, which is a ubiquitous and ecologically important bedform [39,40].We develop a semi-analytical model for BUR, coupling surface and subsurface fluxes with aerobic reactions (respiration and nitrification) within the streambeds.Our model considers solute advection, longitudinal diffusion and first-order kinetics, under the assumption of low nutrient concentrations, with temperature-dependent reaction rate coefficients.
Our first objective is to test the capability of our model to predict DO concentration patterns in a natural setting by comparing predicted and measured DO concentrations.Here, we use the field data reported by Malard et al. [41] and by Rouch [42] for a 15 m-long pool-riffle reach of the Lachein River (Pyrennees, France).Our second objective is to show how the effects of the daily time varying stream BOD, DO and Nr loads (constant or sinusoidally) change the global hyporheic response in terms of BUR and how the local oxic hyporheic pattern changes according to the DO distribution.

Hyporheic Hydraulic Model
We model the hyporheic flow as a Darcian flow [18,43], with the seepage velocity given by Darcy's law: u = −Kh∇h/n, where Kh is the hydraulic conductivity, n is the sediment porosity and h is the energy head.We assume homogenous and isotropic hydraulic properties of the streambed sediment with stationary flow conditions, such that the governing equation is the following: (1) We solve Equation (1) for a three-dimensional alternate bar morphology under the assumptions that the lateral sides (stream banks) and the bottom of the alluvium are impervious and that the hydraulic head gradient is equal to the streambed slope.The stream water surface elevation, which approximates the energy head distribution in the case of submerged bars, is decomposed through the fast Fourier transform.Notice that water surface elevation can be either measured [35,36], obtained through numerical models [40,44] or obtained from analytical hydrodynamic solutions [18,45].Owing to the linearity of Equation (1), the main components of the Fourier decomposition are applied separately as Dirichlet boundary conditions at the upper planar surface that delimitates the subsurface domain.Previous studies showed that replacing the streambed surface with the plane coinciding with the average streambed elevation, as done here, has little effects on the flow and residence time distributions [18,45].The solutions obtained with the single FFT components as the boundary condition are then superimposed to obtain the final head distribution within the subsurface domain.Here, we represent the head distribution from two different sources: field measured heads reported by Malard et al. [41] over a 0 2 = ∇ h pool-riffle unit and analytically predicted by the model of Colombini et al. [46] for the case of a fully-submerged three-dimensional alternate bar in equilibrium with the stream discharge (Figure 1).Once the flow field is known, the transport equation is solved by means of the particle tracking technique [47].The works of Marzadri et al. [45,48] on pool-riffle-induced hyporheic flow show that residence time distributions and biogeochemical transformations depend on bedform size through the ratio Y * BM = HBM/Y0 between the bedform amplitude HBM and the water hydraulic depth, Y0.

Biogeochemical Model
Following the approach of Marzadri, et al. [31], the governing equation for DO transport within the hyporheic zone is written in a Lagrangian framework by applying mass conservation along a stream tube [49]: where α = D/u 2 is the ratio between the longitudinal diffusion coefficient (D) and the module square of the mean velocity (u 2 ) (u = |u|).In addition, KL, KN and KR are the (linear) reaction rates of carbonaceous BOD, nitrogenous BOD and respiration, respectively (Figure 1c).The linear reaction rates of the consumption of DO are used, owing to its relatively low concentration.The last term on the right-hand side of Equation ( 2) accounts for dissolved oxygen consumption due to other biological sources besides the downwelling carbonaceous and nitrogenous BODs.The model neglects the lateral exchange of mass between adjacent stream tubes, and the reaction rates are assumed to be temperature dependent.Because the effects of daily temperature oscillations are small on the biogeochemical response of the hyporheic zone at the bedform scale [50], reaction rates are assumed to change with the mean daily temperature, according to the Arrhenius law [51].
We assume full mixing within the stream, such that the solute concentrations are spatially constant within the downwelling surface, while they vary in time.The initial zero concentration of all solutes is assumed within the hyporheic zone.With these hypotheses, we solve Equation (2) for constant and sinusoidally varying DO, L C and L N stream water concentrations.The latter assumption derives from experimental results showing that daily variations of nutrient concentrations are well approximated by a sinusoidal function [13].Solutions of Equation ( 2) under these initial and boundary conditions are reported in Appendix A.
We introduce two dimensionless numbers: the first is the ratio between the dissolved oxygen concentration in the stream water, DO0, and the total oxygen demand (L0 C + L0 N ), RDO = DO0/(L0 C + L0 N ), and the second is the ratio between carbonaceous and nitrogenous BOD concentrations, RBN = L0 C /L0 N , with the objective to facilitate the analysis of the effects of nutrient loads on BUR.

Benthic Uptake Rate
BUR is defined as the difference between mass fluxes of DO entering and leaving the streambed sediments, which are given by the first and the second term on the right-hand side of Equation (3), respectively:

BUR H x y q x y DO t dxdy dt H x y q x y DO x y t dxdydt
where A = 2 b λ is the streambed area, τ(x,y) is the hyporheic residence time for a particle downwelling at x, y, z = 0 and upwelling after the (residence) time τ since injection, λ is the bedform wavelength (the distance between two pools) and 2b the width of the channel.The time τmax is the maximum hyporheic residence time for the given bedform and defines the time interval over which BUR is quantified.
In addition, t0 is the initial time, which is set to a value larger than maximum residence time τmax in order to eliminate in Equation (3) the effect of the initial conditions.The step function H is unitary within the downwelling areas and is zero outside.Notice that mass balance along a single streamline has been used in the second right-hand term of Equation (3) to transform the integral over the upwelling area in an equivalent integral over the downwelling area [31,35], under the assumption that hyporheic flow is at the steady state.The steady-state flow condition allows us to rewrite Equation (3) in the following form: where q is the mean downwelling flux averaged over the entire streambed area, and time averaged dissolved oxygen concentrations are defined as: (5) and: ( at the downwelling and upwelling areas, respectively. Because the objective of the present analysis is to quantify the contribution of the hyporheic zone to the consumption of stream DO, we make the BUR dimensionless with the mean DO mass flux in the stream: (7) where Qstream is the stream discharge.Similarly, the BUR can be made dimensionless with respect to the hyporheic flux by replacing in Equation ( 7) Qstream with q A: (8) Hereafter, this dimensionless expression of BUR will be indicated with BUR * HZ.
The following two end members can be defined for BUR * when all of the upwelling fluxes are anoxic, such that all of the oxygen carried by the stream is consumed within the hyporheic zone: (9) which have been obtained by considering that carbonaceous and nitrogen oxygen consumption is limited by the DO of stream water, such that = when RD0 < 1 and BUR = L0 C + L0 N when RD0 > 1. Equation ( 8) provides a first approximation analysis for BUR * , which depends only on the hyporheic (q 2bλ) and stream (Qstream) discharges, the hyporheic residence time distribution and the nutrient loads (RDO = DO0/(L0 C + L0 N )) in the stream.The hyporheic values, τmax and q, can be estimated with models presented in the literature [18,19,45,[52][53][54][55].

Field Data
The works of Malard et al. [41] and Rouch [42] report both the near-bed head and the dissolved oxygen concentration distributions at the sediment-water interface, supplemented by streambed sediment hydraulic conductivity and porosity along a 15 m-long pool-riffle bar reach of the Lachein stream (Pyrennees, France).In the present work, the head distribution derived from the contour lines depicted in Figure 3c of Malard et al. [41] by superimposing a regular two-dimensional grid is first decomposed via FFT (Figure 2).Owing to the linearity of the flow Equation ( 1) and Darcy's law, the most relevant FFT components are used as the upper Dirichlet boundary condition for modeling the 3D hyporheic flow field, through the solution of Laplace and Darcy equations, separately for each component, under the assumption of homogeneous and isotropic alluvium, as described in Section 2.1.The resulting velocity fields are then superimposed to obtain the final three-dimensional velocity field distribution.We set the alluvium thickness to zd = 2 m, which is deep enough to make negligible its effect on the hyporheic flow field, the hydraulic conductivity to Kh = 8.39 × 10 −4 m/s, which is the mean of the hydraulic conductivity values reported in Figure 3b of Malard et al. [41], and the porosity to n = 0.2664, which is the mean of the porosity values reported by Rouch [42] for this study site.The residence time of 8192 particles injected within the downwelling area was obtained by particle tracking with the analytical velocity obtained as explained above.We do not account for particles whose trajectories do not upwell in the stream.For the biogeochemical model, we assumed a constant dissolved oxygen concentration in the stream water, DO0 = 10 mg/L, which is the mean of the measured values over the downwelling areas, and a saturated DO concentration of DOsat = 11.83mg/L, which corresponds to the observed water temperature of T = 8 °C.In addition, we set the reaction rate of respiration equal to the mean value calculated from the field data reported by Beaulieu et al. [56] for streams with a water temperature near 20 °C: KR(20 °C) = 5.01 d −1 .We neglect all of the other parameters, DOlim, L C , L N and α, whose values were not reported in the work of Malard et al. [41] and Rouch [42].The coefficient KR implicitly accounts for the consumption of DO due to both carbonaceous and nitrogenous BOD, whose values were not reported, and thus, we set KL and KN equal to 0. These working assumptions allow us to use only measured values, avoiding the need to calibrate the model.We used the error matrix technique, which is a common GIS (Geographic Information System) method for map comparison [57][58][59], to assess the accuracy of the modeled spatially-distributed DO concentrations against field observations.The error matrix provides two indexes of agreement between predicted and observed spatially-distributed quantities on a cell-by-cell basis: an overall agreement, OA, and a kappa index.The OA is calculated as the proportion of correctly-predicted cells versus the total number of cells.Differently from OA, the kappa index removes the effect of by-chance agreement.In fact, the probability that two successive random extractions of a categorical variable belong to the same category is larger than zero, by definition.Since the kappa index filters out the match of predicted and observed DO cells due to random coincidence, its value is typically lower than OA [60].Both indexes vary from 0 to 1 with 0 being poor and 1 perfect agreement.For the K statistic, Landis and Koch [59] suggested the following classification: <0 = poor agreement, 0-0.19 slight, 0.20-0.39fair, 0.4-0.59moderate, 0.60-0.79substantial and 0.80-1.00almost perfect agreement.We calculated OA and K with the DO concentrations divided into the 3 categories reported by Malard et al. [41].The predicted and observed DO category needs to coincide exactly to be counted as a match and not an error.This is a strict method, because it does not take into account if the match is a little bit offset spatially; for instance, a nearby cell could have the correct value [60].

Numerical Simulations
We considered 18 streams characterized by an alternate bar morphology, whose bedforms are in equilibrium with the flow discharge.The bedform dimensions, whose morphodynamic parameters are reported in Table 1, are within those typically encountered in natural streams [39,[61][62][63].For each stream, we ran 10 simulations, changing the nutrient loads (Table 2).We ran a set of simulations with constant loads and a set with sinusoidal varying CBOD, NBOD and DO concentrations with the same mean value as the constant case to explore the effects of time-dependent nutrient releases.Table 1.Morphodynamic parameters of the pool-riffle bedform.β is the aspect ratio (β = b/Y0); θ is the Shields' number; ds is the relative submergence of the particles; d50 is the median grain size; Y * BM is the dimensionless mean flow depth; s is the streambed slope; ∆ is the bar amplitude; λ is the bar length; and Qstream is the stream water discharge.Table 2. Biogeochemical parameters used in the simulations.RDO is the ratio between dissolved oxygen concentration (DO) and the total oxygen demand; RBN is the ratio between carbonaceous (L C ) and nitrogenous (L N ) BOD; L0 C and L0 N are the constant stream concentrations of carbonaceous (CBOD) and nitrogenous (NBOD) BOD, respectively.LA C and LA N are the stream daily fluctuations of CBOD and NBOD, respectively.KL, KN and KR are the decay rates of CBOD, NBOD and DO, respectively.In all of the simulations, the concentration of dissolved oxygen concentration DO0 = 10 mg/L, while the stream daily fluctuation of oxygen is DOA = 4 mg/L.These simulations are tailored to investigate the relationship between BUR * , RBN and RDO, which identify the nutrient loads, and Y * BM, an index of stream size.Large streams are typically characterized by large Y * BM and long residence times, whereas small streams by small Y * BM and short residence times [31].

Comparison with Field Data
Figure 3 shows the contour map of the computed hyporheic residence times (Figure 3a), the categorical map of measured dissolved oxygen concentration (red contours) (Figure 3b, which corresponds to Figure 4d in Malard et al. [41]) and the map of model-predicted dissolved oxygen concentrations (Figure 3c) at the sediment-water interface.The overall agreement value between the predicted and observed DO concentration distribution is 0.58, and the kappa index is 0.46, which correspond to a good agreement and are comparable to the values reported in the GIS literature [57,59,64].Because, to our knowledge, this is the first comparison between observed and modeled spatially-distributed hyporheic quantities, we cannot compare these values with others reported in the hyporheic literature.This result indicates that hyporheic processes are those that chiefly control the DO concentration distribution at and below the water-sediment interface and, consequently, affect the habitat of benthic and hyporheos organisms.The predicted patterns of τ and DO concentrations also match the contours of the categorical map of DO concentrations measured in the field (Figure 4d in [41]).This further confirms that residence time distribution is a key factor in the fate of reactive solutes within the hyporheic zone.In our case, the transformations depend entirely on the reaction time, because the biogeochemical reaction rates are spatially and temporally constant.Malard et al. [41] observed that organisms colonized preferential areas of the streambed depending on DO concentrations and described the behavior in terms of morphological units: bar, pool and bar fringe.Patches with a high dissolved oxygen concentration (>74% of saturation) coincide with depositional areas, such as bars, which are mostly predicted as downwelling zones, and bar fringes, which are predicted as upwelling areas with a short residence time.These zones show dense communities of harpacticoid and hypogean species; whereas streambed areas with low DO concentrations, which are identified as erosional areas by Malard et al. [41], like the pool, and for which our model predicts upwelling areas with long residence times, have fewer harpacticoids and hypogean organisms.All of these results suggest that the near-bed DO concentration is chiefly controlled by stream bedforms, which therefore control the habitat within and at the streambed sediment.

Effect of Stream Morphology
Figure 4a shows BUR * , i.e., the BUR made dimensionless with respect to the stream flux, as a function of the dimensionless depth, Y * BM, for hyporheic zones induced by three-dimensional pool-riffle morphology for RDO = 1.2 and RBN = 0.5.In addition, Figure 4b shows how the BUR (BUR * HZ) behaves when it is normalized with respect to the hyporheic flux.BUR * initially increases with river size: it reaches a peak and then decreases for large streams.The peak for BUR * is caused by the initial faster increase of downwelling fluxes with long hyporheic residence times, which potentially allow more DO consumption, than stream discharge, which also increases with stream size [54].Eventually, the increase in stream discharge is large enough to offset the increase in hyporheic exchange, such that BUR * decreases.
Differently from BUR * , BUR * HZ increases monotonically with river size, because of the longer residence time that water spends within the hyporheic zone as the size of the river increases.Whereas the BUR * HZ increase is fast around Y * BM = 1, it slows down significantly at larger values of Y * BM.This suggests that short flow paths may be always present, reducing the efficiency of the hyporheic zone to consume DO to a value smaller than 100%, even for large streams.
Figure 5 shows a snapshot of hyporheic DO concentrations along the streamlines of the 3D hyporheic flow field provided by our model.Within the alluvium, the concentration of DO decreases with depth according to the time that water spends traveling along a given streamline.
The longer the path is, the greater is the probability to find streamlines with DO < 2 mg/L (red path of the trajectory) and, thus, volumes of sediment in anaerobic conditions.The complex DO distribution depicted in Figure 5 demonstrates the importance of considering the entire three-dimensional flow patterns when studying alternate bar morphologies.As the stream dimension increases, the volume of streambed sediment saturated with stream water by the hyporheic flow increases and, consequently, so do the residence times.Thus, in large streams, the red portion of the flow lines would dominate.However, the presence of short residence times at the transitional zone between upwelling and downwelling areas would still have a blue tone.This supports the results of Figure 4, which suggest that BUR * HZ may be limited to a value below 100%.However, this upper threshold depends also on biogeochemical rates, because faster rates would require shorter times to reach anoxic conditions [32,65].  1.

Effects of Constant BOD and DO Loads
Figure 6 shows BUR * as a function of Y * BM in the case of three-dimensional alternate bars for several RDO values with RBN = 0.5 (Figure 6a) and RBN = 2 (Figure 6b).Keeping in mind the discussion in Section 3.2 concerning the mechanism leading to the peak in the BUR * (Figure 4a), here, we can observe the effect of nutrient concentrations on the location of this peak.In particular, as RDO decreases (because stream DO decreases or L0 C + L0 N increases), the total oxygen demand (L C + L N ) is larger than the stream DO concentrations; the BUR * peak increases, and it occurs in smaller streams (smaller Y * BM) for low rather than high RDO.For instance, the BUR * peak occurs at Y * BM = 0.75 for RDO = 1.2 and at Y * BM = 0.48 for RDO = 0.2, regardless of RBN.Both cases with a different RBN show similar patterns with slightly higher BUR * associated with larger RBN.The faster kinetics of L C than L N justifies this behavior.However, the effect of nutrient loads on BUR * is important only for small streams.As the stream size increases, both the RDO and RBN effects become negligible, and all curves tend to the same value.Our simulation suggests that the threshold is near Y * BM = 1.

Effects of Sinusoidally Varying and Constant Stream Solute Concentrations
Figure 7 shows BUR * for sinusoidally varying and constant L C and L N loads as a function of Y * BM for RBN = 0.5 with RDO = 0.2 (Figure 7a) and RDO = 1.2 (Figure 7b).The most important result shown in Figure 7a,b is that including daily variations of the stream DO provides little benefit with respect to the use of a constant stream DO for large streams.Some effects of the time-varying in-stream concentration are observed for small streams.This behavior can be explained by considering the hyporheic residence time as the main controlling factor of DO consumption for the same biological activity.Path lines with long residence times show larger portions with low DO concentrations, which are less affected by the daily DO fluctuations of stream water.Conversely, large portions of path lines with short residence times show DO concentrations that vary temporally, depending on stream nutrient and DO concentrations.Thus, systems with predominantly long residence times, like streams with large Y * BM, show a similar global response for both boundary conditions; whereas the BUR * response of streams with small Y * BM shows a dependence on temporal variations of stream nutrient and DO concentrations.However, as shown in Figure 7a,b, the dependence on DO fluctuations in stream water is small, yet appreciable, thereby suggesting that the mean daily concentrations of DO are sufficient information for studying the global response of the hyporheic zone in most cases of practical interest.
BUR * shows the overall system response, which is certainly important for understanding the mean behavior of the system for different types of stream nutrient loads (Figure 7a,b).However, the lifetime of many aquatic organisms is strongly related to the local behavior of the system, which could be quite different from the overall response.In this context, Figure 8 shows a three-dimensional snapshot of simulated DO concentrations within the hyporheic zone in the case of sinusoidally-periodic (Figure 8a) and constant (Figure 8b) DO, L C and L N concentrations for RDO = 1.2 and RBN = 0.5.The variation in time of the two nutrients and DO in the stream water induce differentiated spatial and temporal distributions of DO concentrations within the hyporheic zone.The temporal changes of DO concentrations are not systematic within the domain, but they are relevant for short residence times and become negligible along the streamlines as the residence time increases.Thus, the hyporheic zone may show a dynamic volume of sediment, whose DO concentrations fluctuate following those of the stream water, and a stable volume, where DO concentrations are low and stationary.We should expect a spatial specialization of organisms dwelling in the hyporheic zone following the DO profile as those observed at the sediment-water interface by Malard et al. [41].

Conclusions
We presented an analytical model, which accounts for advection, diffusion and both carbonaceous and nitrogenous BODs, for predicting DO concentrations within the hyporheic zone.The application of the model on a gravel-bed pool-riffle reach shows that the proposed semi-analytical tool is able to predict the spatial distribution of measured DO concentrations.The overall agreement between predicted and observed DO concentrations is 58%.Our results show that hyporheic DO consumption varies as a function of bedform morphology, size and nutrient loads.Our model uses spatially-constant reaction rates, which probably may change spatially, and homogenous hydraulic conductivities.These two limitations may be important to detect micro-biogeochemical hotspots, which may be linked to localized slow velocities [65].
Our simulations suggest that river morphology plays an important role in controlling the BOD removal capacity of the hyporheic zone.However, the system responds differently at local and global scales, depending on the temporal variations of stream DO concentrations.The bedform scale response of the hyporheic zone expressed by BUR * is only minimally influenced by daily oscillations of stream water nutrient and DO concentrations.Thus, daily mean concentrations of these solutes can be used as an index for quantifying the hyporheic response at the bedform scale.Conversely, daily oscillations of stream nutrient concentrations may influence the spatial distribution and density of aquatic organisms locally.Finally, we show that for a given nutrient concentration, there is a stream size that maximizes the removal of L C and L N and, consequently, the consumption of stream dissolved oxygen.
for constant and sinusoidally time varying DO concentration in stream water, respectively.In Equation (A1), L0 C and L0 N are, respectively, the constant concentrations of L C and L N in the stream water.On the other hand, in Equation (A2), LA C and LA N and DOA are the amplitude of the in-stream oscillation signal for L C , L N and DO, respectively, and ωL, ωN and ω are their periods of oscillation.
The solution of Equation ( 2) with initial and boundary conditions expressed in Equation (A1), which provides the contribution of one stream tube in DO transport, assumes the following form: ( In addition, the solution of Equation ( 2) with initial and boundary conditions expressed in Equation (A2), which, similarly to the previous case, provides the contribution of one stream tube in DO transport, assumes the following form: ( where Fi (i = 1 ÷ 6) and Gi (i = 1 ÷ 6) are the functions reported in Marzadri et al. [45] When the stationary conditions are reached, the solution of Equation (A4) further simplifies to: Carbonaceous (L C ) and nitrogenous (L N ) BOD, transport and transformation within the hyporheic zone are modeled solving the transport Equations (A7) and (A8) for both constant (A1) and sinusoidally (A2) varying in-stream solute concentrations: The solution of Equations (A7) and (A8) with initial and boundary conditions expressed in Equation (A1), which provides the contribution of one stream tube in carbonaceous and nitrogenous BOD transport, respectively, assumes the following forms:  In addition, the solution of Equations (A7) and (A8) with initial and boundary conditions expressed in Equation (A2), which, similarly to the previous case, provides the contribution of one stream tube in carbonaceous and nitrogenous BOD transport, respectively, assumes the following forms:  ( ) where Fi (i = 1 ÷ 6) and Gi (i = 1 ÷ 6) are the functions reported in Marzadri et al. [45] and in Equation (A5).Under steady-state conditions, the solution of Equations (A11) and (A12) simplify to: ( )

Figure 1 .
Figure 1.Sketch of (a) the planar view and (b) the cross-section of the three-dimensional alternate bar morphology; and (c) diagram of the biogeochemical model.DO, dissolved oxygen; CBOD, carbonaceous biological oxygen demand; NBOD, nitrogenous biological oxygen demand.

Figure 2 .
Figure 2. Head distribution in the 15 m-long reach of the Lachein stream measured by Malard et al.'s [41] contour lines and reconstructed by FFT (color plot).

Figure 3 .
Figure 3. Maps of: (a) the computed residence time distribution of water within the hyporheic zone, τup; (b) the measured dissolved oxygen concentrations, (DO), in % of saturation reproduced from Figure 4d in Malard et al. [41]; and (c) the model-predicted DO along the 15-m reach of the Lachein stream.

Figure 4 .
Figure 4. (a) Stream, BUR * , and (b) hyporheic, BUR * HZ, dimensionless BUR, as a function of dimensionless mean flow depth, Y * BM , for three-dimensional alternate bars with RDO = 1.2 and RBN = 0.5, with RDO being the ratio between in-stream DO (DO0) and total oxygen demand (L0 C + L0 N ) and RBN the ratio between in-stream carbonaceous and nitrogenous BOD concentrations, RBN = L0 C /L0 N .BUR, benthic uptake rate.

Figure 5 .
Figure 5. Snapshot of dissolved oxygen concentrations (DO) simulated within the hyporheic zone of a river with the morphological parameters of Test 5 in Table1.

Figure 6 .
Figure 6.Dimensionless stream BUR, BUR * , as a function of the dimensionless mean flow depth, Y * BM, and RDO for (a) RBN = 0.5 and (b) RBN = 2.0, with RDO the ratio between in-stream DO (DO0) and total oxygen demand (L0 C + L0 N ) and RBN the ratio between in-stream carbonaceous and nitrogenous BOD concentrations, RBN = L0 C /L0 N .

Figure 7 .
Figure 7. Dimensionless stream BUR, BUR * , as a function of dimensionless mean flow depth, Y * BM, for periodic and constant in-stream boundary conditions and RBN = 0.5 with (a) RDO = 0.2 and (b) RDO = 1.2, with RDO the ratio between in-stream DO (DO0) and total oxygen demand (L0 C + L0 N ) and RBN the ratio between in-stream carbonaceous and nitrogenous BOD concentrations, RBN = L0 C /L0 N .

Figure 8 .
Figure 8. Maps of dissolved oxygen concentrations, DO, within the hyporheic zone for (a) periodic and (b) constant in-stream concentrations of DO, CBOD and NBOD.In both cases, RDO = 1.2 and RBN = 0.5 (Test 5-A), with RDO the ratio between in-stream DO (DO0) and total oxygen demand (L0 C + L0 N ) and RBN the ratio between in-stream carbonaceous and nitrogenous BOD concentrations, RBN = L0 C /L0 N .