Modeling surface water and groundwater mixing and mixing-dependent denitrification with bedform dynamics

The hyporheic zone (HZ), where surface water (SW) and groundwater (GW) interact and mix, acts as a critical interface that attenuates contaminants through enhanced biogeochemical cycling. While bedform migration significantly influences hyporheic exchange and non-mixing-driven reactions of solutes from upstream SW, the effects of bedform migration on SW-GW mixing dynamics and mixing-triggered biogeochemical reactions – particularly under gaining stream conditions – remain poorly understood. Establishing a coupled hydrodynamic and reactive transport model that incorporates bedform migration, this paper systematically examines nitrogen processing for scenarios of variable sediment grain size, stream velocities, and upwelling GW fluxes. Results of this study reveal that SW-GW mixing and mixing-triggered denitrification zones progressively transition from crescent shapes into uniform band-like configurations as bedforms migrate. Both hyporheic exchange flux and mixing flux increase with increasing stream velocity and associated bedform celerity. The mixing proportion and mixing zone size increase at the start of migration, while they reach approximately stable when turnover becomes the dominant hyporheic exchange mechanism. Slow to moderate migrated bedforms with enhanced mixing dynamics facilitate mixing-triggered denitrification, whereas fast stream flows and migrating bedforms shorten solute residence timescales and limit denitrification potential. Consequently, in fine to medium sandy sediments, groundwater-borne nitrate removal efficiency declines significantly with bedform migration. The self-purification capacity of the HZ, and particularly its functioning as a natural barrier against GW contamination, is hindered under such dynamic bedform conditions. These findings highlight the need to maintain stable bedform conditions in restoration projects to enhance the capacity of HZ contaminant attenuation.

1 Introduction

Anthropogenic activities such as the intensification of agricultural practices with its increased used of mineral and organic fertilizers, together with high livestock densities and emissions of inadequately treated domestic and industrial wastewater have significantly increased nitrogen loading to rivers and groundwater, which impacts water quality, causing eutrophication, hypoxic and related deterioration of ecosystem functions (Conley et al., 2009; Rouse et al., 1999). Long-term regulatory monitoring data (e.g., from the UK) indicate that nitrate levels have stabilized in many rivers, while nitrate concentrations in groundwater-fed rivers continue to increase (Burt et al., 2011; Howden and Burt, 2008). The persistence of nitrate contamination in groundwater and associated risks of a “nitrate time bomb” (Ascott et al., 2019) has highlighted the urgency of exploring the potential of natural microbial processes to mediate nitrate transformation and removal in riverbed sediments (Shelley et al., 2017; Lansdown et al., 2015; Rivett et al., 2008).

The hyporheic zone (HZ) has received significant attention for its potential to facilitate enhanced nitrate transformation and removal via denitrification that is a primary process permanently reducing nitrate, with hyporheic exchange flows (HEFs) acting as a critical mechanism for transporting nitrate-rich surface water to the riverbed sediments where microbial activities and biogeochemical reaction rates are enhanced (Boano et al., 2014; Boulton et al., 1998; Cardenas, 2015; Xian et al., 2022; Krause et al., 2022). It has for long been assumed that predominantly stream waters provide inputs of bioavailable (mainly dissolved) organic carbon (DOC), oxygen (O₂) and nitrate (NO${}_{\mathrm{3}}^{-}$) into the riverbed where the residence and reaction times determine the occurrence of aerobic respiration and the potential for shifts into anaerobic conditions that may facilitate denitrification along the HEF paths in the presence of enough remaining DOC (Zarnetske et al., 2011a, b). These hydrological and biogeochemical mechanisms are in this form mainly representative of headwater streams, where the HEF is induced by stream flow turbulence (Boano et al., 2011; Roche et al., 2018, 2019), local geomorphological setting (Cardenas et al., 2008; Marzadri et al., 2012; Tonina and Buffington, 2007), and flow obstacles such as woody debris, streambed engineering or restoration structures (Briggs et al., 2013; Wondzell et al., 2009), and hyporheic nitrate removal processes expected to mainly occur at the middle to end- hyporheic flow paths within HEF cell sediments.

In lowland systems, groundwater-fed streams and rivers in permeable catchments will be dominated by base flow of often nitrate enriched groundwaters. The subsurface hydrological conditions are driven by horizontal HEF as well as vertical fluxes resulting from regional groundwater flow toward (or from) the stream. Spatial variability in regional groundwater flow can significantly affect hyporheic exchange and biogeochemical cycling (Boano et al., 2013; Krause et al., 2013; Munz et al., 2011; Naranjo et al., 2015). It is important to note that the groundwaters of many piedmont plains and lowlands are often contaminated with high nitrate concentrations, but usually low in bio-available DOC (Krause et al., 2022). When nitrate-rich groundwater upwells through deeper sediments and reaches a region enriched in availability of DOC, nitrate reduction processes are significantly facilitated in the condition (Krause et al., 2009; Naranjo et al., 2015; Ping et al., 2023; Trauth et al., 2017). Stelzer and Bartsch (2012) developed such a conceptual framework of nitrate-rich gaining fluvial setting from 8 sites in the Waupaca River Watershed with three order magnitude in groundwater nitrate concentration. Lansdown et al. (2014) also measured high denitrification rate at deep sediment in the coarse-grained sediments typical of groundwater-fed system, located within the River Leith (Cumbria, UK) where diverse nitrogen concentration changes were confirmed earlier (Krause et al., 2009).

Turnover and removal of the large amounts of NO${}_{\mathrm{3}}^{-}$ from groundwater require DOC either from autochthonous streambed sources or from downwelling surface water to stimulate nitrate reduction (Krause et al., 2013, 2022; Ping et al., 2023; Sawyer, 2015; Trauth et al., 2017). For coarse grain or sandy riverbeds with low autochthonous organic matter content, stream-borne DOC dominates the supply of carbon sources for nitrogen transformation processes. Sandy sediments with less autochthonous organic carbon sources covering the majority of alluvial riverbeds are commonly characterized by topographical features such as ripples, dunes, and riffle-pool sequenced. The typical and multiple HEF cells induced by bedform topography are generally in crescent shapes (Fox et al., 2014; Wu et al., 2024). The downward advection of stream borne DOC provides electron donor and mixes with nitrate-rich and anoxic groundwater. It has been shown that the highest potential for mixing triggered denitrification is often found at the margin of the HEF cells, which represents the last natural protection before nitrate enter a stream (Hester et al., 2013, 2014; Gu et al., 2008; Nogueira et al., 2022). The effects of mixing triggered denitrification on groundwater borne nitrate transformation in HZs have been studied almost exclusively for the case of stationary, that is immobile bedforms (Hester et al., 2017, 2019; Trauth and Fleckenstein et al., 2017; Ping et al., 2023).

Bedforms are mobile in dynamic equilibrium or undergo constant changes during periods of moderate to high stream flow, and they are typically found in medium and larger waterways under realistic field conditions (Bartholdy et al., 2015; Risse-Buhl et al., 2023; Schindler et al., 2015). For example, Harvey et al. (2012) observed the migrating bedforms of dunes (with a median grain size D₅₀ of 380 µm for the riverbed sediments) at a velocity of 57.6 cm h⁻¹ during base flow in the “Clear Run” stream in eastern North Carolina, USA. Ahmerkamp et al. (2017) found that the ripple bedforms for sands (D₅₀=63 µm) ranged from 11 to 29 cm with a constant ratio of bedform height and length at $\mathrm{1}/\mathrm{9}$, and migrated at velocities of 0.7–6.5 cm h⁻¹ in the German Bight, Southeastern part of the North Sea. Bedforms migration complicates the development of hyporheic flow fields, facilitates solute exchange, alters redox conditions in riverbeds, and affects contaminant transport and transformation (Ahmerkamp et al., 2015; Schulz et al., 2023; Peleg et al., 2024). Previous studies focused on ripple-driven hyporheic exchange in headwater to midstream reaches under neutral conditions, and simulated bedform migration via a moving frame of reference (Ahmerkamp et al., 2015; Jiang et al., 2022; Kessler et al., 2015; Zheng et al., 2019). Ping et al. (2022) used a fixed reference frame and noted that bedform migration negatively impacts non-mixing-dependent denitrification rates and stream-borne nitrate removal efficiency. However, no studies have yet investigated the effects of bedform migration on surface water-groundwater mixing in the lowland reaches of streams under gaining conditions, nor have they explored its controls or implications for mixing-induced denitrification in groundwater-fed streams and rivers.

In this study, numerical modeling of hyporheic flow and multi-component solute (DOC, O₂, and NO${}_{\mathrm{3}}^{-}$) transport is used to evaluate the effect of bedform migration on mixing-dependent denitrification in the HZ of a gaining river, where the overlying stream water is induced into the sediment by ripple-type bedforms. The objectives of this study are to determine the effects of bedform migration on the overall extent and magnitude of mixing of surface water and upwelling groundwater, as well as its influences on groundwater-borne nitrate transport and transformation.

2 Methods

2.1 Model description

The stream flow over the bedform generates head gradients along the riverbed surface, drives hyporheic exchange, and induces the formation of periodic ripples. Ripples develop and migrate downstream above the riverbed surface as a consequence of bedload transport, with sediment grain eroded and deposited at the surface, whereas deeper streambed sediments remain relatively stable (Harvey et al., 2012; Precht et al., 2004). From the perspective of an observer who stands at the target single riverbed and watches the ripples passing by, the bedform shapes and associated pressure profiles migrate downstream over the immobile riverbed layer. Because ripple heights (H_d=2 cm) accounts for 20 % of the stream water depth (H=10 cm) and are far smaller than that of the underlying immobile domain (l=16 cm), the undulating riverbed can be reasonably approximated as a flat bed with time-varying pressure patterns during bedform migration. Besides, reactive solute transport also primarily takes places within the immobile riverbed sediment. Therefore, a two-dimensional (2D) rectangular domain is constructed with a length (λ) of 0.2 m and depth (l) of 0.16 m. A fixed reference frame is adopted, with its origin located at the bottom-left corner of the model. Here, the horizontal direction corresponds to the x-axis and the vertical direction to the y-axis. The model domain has four corners: A (0, 0), B (0, l), C (λ, l), and D (λ, 0). The lateral boundaries are AB (left, x=0) and CD (right, x=λ); the top boundary is BC (y=l), and the bottom boundary is AD (y=0).

The model schematic is sketched in Fig. 1. The model is relevant to the model used by Ping et al. (2022) but with significant differences. Specifically, Ping et al. (2022) focused on ripple-driven hyporheic exchange processes in the headwater to midstream sections of streams under neutral conditions. In contrast, the present study focuses on the downstream sections of rivers under gaining conditions, aiming to simulate surface water-groundwater interacting and mixing processes under dynamic bedforms. To our knowledge, previous migrating ripple models have extensively examined surface water-groundwater interaction processes (i.e., ripple-driven hyporheic exchange processes), yet the mixing process remains uninvestigated (Ahmerkamp et al., 2017; Jiang et al., 2022; Zheng et al., 2019; Kessler et al., 2015). In downstream sections, groundwater is typically nitrate-contaminated, and denitrification is often co-regulated by both hyporheic exchange and mixing regimes. The model is also developed with advancements of these previous models by including an important biogeochemical reaction in nitrate dynamics, that is mixing-triggered denitrification (see next section for details).

Figure 1 Schematic of the model domain with bed form geometry. Stream flow and bedform migration are from left to right. The dashed yellow lines represent the streambed surface and head profile at time t₀, and the solid black lines represent the streambed surface and head profile after migration at t₁. “SW in” is where surface water enters the riverbed layer, “SW out” is where surface water discharges to the stream, and “GW out” is along the upstream and downstream sides of bedforms where groundwater discharges to the stream. The yellow band represents the mixing zone of surface water and groundwater. NMD reaction = non-mixing-dependent reaction and MD reaction = mixing-dependent reaction. This figure is modified from Fig. 2 in Ping et al. (2022), with permission granted by Wiley.

In this study, hyporheic flow, multi-component solutes (DOC, O₂, and NO${}_{\mathrm{3}}^{-}$) transport, and biogeochemical reactions (both non-mixing-dependent NMD and mixing-dependent MD reactions) are modelled in saturated sediments beneath a riverbed layer. Stream flow over the riverbed is not simulated. The hyporheic exchange flow is driven by pumping (overlying flow over the bed form produces a pressure head $h{|}_y{=}_l$) and bedform migration (triangular shaped ripples migrate downstream by an average velocity u_c) and influenced by upwelling groundwater. Stream water transports DOC, O₂, and NO${}_{\mathrm{3}}^{-}$ into the riverbed, whereas groundwater polluted with NO${}_{\mathrm{3}}^{-}$ upwells toward the stream. In order to distinguish NMD and MD denitrification, we define stream-borne nitrate (s-NO${}_{\mathrm{3}}^{-}$) and groundwater-borne nitrate (g-NO${}_{\mathrm{3}}^{-}$) as two distinct reactants.

2.2 Model formulation

2.2.1 Governing equation

The pore water flow was calculated using Darcy's law and the groundwater flow equation:

$$\begin{array}{}\text{(1)}& \mathrm{\nabla }\cdot \left(-K\mathrm{\nabla }h\right)=\mathrm{0}\end{array}$$

where h [L] is the pressure head, and K [L T⁻¹] = $k\mathit{\rho }g/\mathit{\mu }$ is the hydraulic conductivity, k [L²] is the permeability of the riverbed, ρ [M L⁻³] are the density of water, g [L T⁻²] is the gravity acceleration, and μ [M L⁻¹ T⁻¹] is the dynamic viscosity of water.

The transport of reactive solutes within the streambed sediment was described by the advection-dispersion-reaction equation:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial c_i}{\partial t}}-\mathrm{\nabla }\cdot \left(D_{ij}\mathrm{\nabla }{c}_i\right)+\mathrm{\nabla }\cdot \left(\mathit{v}{c}_i\right)=R_i\end{array}$$

where c_i [M L⁻³] represents the concentration of reactive components, v [L T⁻¹] is the seepage or linear pore water velocity vector, D_ij [L² T⁻¹] is the hydrodynamic dispersion tensor, and R_i [M L⁻³ T⁻¹] denotes the biogeochemical rate of reactive components; the elements of the dispersion tensor D_ij are defined by Bear and Verruijt (1998):

$$\begin{array}{}\text{(3)}& D_{ij}=\left({\mathit{\alpha }}_{\mathrm{L}}-{\mathit{\alpha }}_{\mathrm{T}}\right)\cdot {\displaystyle \frac{v_i{v}_j}{\left|\mathit{v}\right|}}+{\mathit{\delta }}_{ij}\cdot \left({\mathit{\alpha }}_{\mathrm{T}}\left|\mathit{v}\right|+\mathit{\theta }\cdot \mathit{\iota }{D}_{\mathrm{m}}\right)\end{array}$$

where α_L [L] and α_T [L] are longitudinal and transverse dispersivities, respectively, D_m [L² T⁻¹] is molecular diffusion coefficient, v_i and v_j denote the seepage velocities in the i-direction and j-direction, respectively; and ι [–] is tortuosity.

The biogeochemical reactions: aerobic respiration (AR), non-mixing-dependent, and mixing-dependent denitrification (DN) are considered in the reactive transport model:

$$\begin{array}{}\text{(R1)}& {\displaystyle }\mathrm{AR}:{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{CO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\text{(R2)}& {\displaystyle }\begin{array}{rl}\mathrm{DN}:& {\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+\mathrm{0.8}{\mathrm{NO}}_{\mathrm{3}}^{-}+\mathrm{0.8}{\mathrm{H}}^{+}\\ & \to {\mathrm{CO}}_{\mathrm{2}}+\mathrm{1.4}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}+\mathrm{0.4}{\mathrm{N}}_{\mathrm{2}}\end{array}\end{array}$$

Reactions (R1) and (R2), with DOC as the electron donor and O₂ and NO${}_{\mathrm{3}}^{-}$ as sequential electron acceptors, capture the primary mechanism of NO${}_{\mathrm{3}}^{-}$ cycling and are widely used in studies on hyporheic zone nitrogen dynamics (Bardini et al., 2012; Hester et al., 2019; Zheng et al., 2019). Denitrification is the well-recognized critical process for NO${}_{\mathrm{3}}^{-}$ transformation and reduction in riparian and hyporheic zones. The reaction stoichiometry and electron acceptor utilization order in Reactions (R1) and (R2) match well-established principles from lab incubations and field investigations (Hedin et al., 1998; Liu et al., 2017; Zarnetske et al., 2011a, b). Dissimilatory nitrate reduction to ammonium (DNRA) and anaerobic ammonium oxidation (ANAMMOX) are not included in the model, as these processes play secondary role on nitrogen cycling and require highly specific conditions to occur (Zarnetske et al., 2012). Ammonia (NH${}_{\mathrm{4}}^{+}$) was excluded, as it is unstable in the study environment: NH${}_{\mathrm{4}}^{+}$ in surface water or groundwater is prone to nitrification (converting NH${}_{\mathrm{4}}^{+}$ to NO${}_{\mathrm{3}}^{-}$) either within hyporheic flow cells or in upwelling groundwater upon O₂ exposure, leading to relatively low NH${}_{\mathrm{4}}^{+}$ concentrations (Hester et al., 2014).

Reaction rates were defined using Monod kinetics (Monod, 1949):

$$\begin{array}{}\text{(4)}& {\displaystyle }{R}_{\mathrm{AR}}=V_{\mathrm{AR}}\times X_{\mathrm{AR}}\times {\displaystyle \frac{c_{{\mathrm{O}}_{\mathrm{2}}}}{c_{{\mathrm{O}}_{\mathrm{2}}}+K_{{\mathrm{O}}_{\mathrm{2}}}}}\times {\displaystyle \frac{c_{\mathrm{DOC}}}{c_{\mathrm{DOC}}+K_{\mathrm{DOC}}}}\text{(5)}& {\displaystyle }\begin{array}{rl}{R}_{\mathrm{DN}}& =V_{\mathrm{DN}}\times X_{\mathrm{DN}}\times {\displaystyle \frac{c_{{\mathrm{NO}}_{\mathrm{3}}^{-}}}{c_{{\mathrm{NO}}_{\mathrm{3}}^{-}}+K_{{\mathrm{NO}}_{\mathrm{3}}^{-}}}}\\ & \times {\displaystyle \frac{c_{\mathrm{DOC}}}{c_{\mathrm{DOC}}+K_{\mathrm{DOC}}}}{\displaystyle \frac{K_{\mathrm{inh}}}{K_{\mathrm{inh}}+c_{{\mathrm{O}}_{\mathrm{2}}}}}\end{array}\end{array}$$

The reaction terms R_i was given by

$$\begin{array}{}\text{(6)}& {\displaystyle }{R}_{{\text{s-NO}}_{\mathrm{3}}^{-}}=-R_{\mathrm{DN}}\left(c_{{\text{s-NO}}_{\mathrm{3}}^{-}}\right)\text{(7)}& {\displaystyle }{R}_{{\text{g-NO}}_{\mathrm{3}}^{-}}=-R_{\mathrm{DN}}\left(c_{{\text{g-NO}}_{\mathrm{3}}^{-}}\right)\text{(8)}& {\displaystyle }{R}_{{\mathrm{O}}_{\mathrm{2}}}=-R_{\mathrm{AR}}\text{(9)}& {\displaystyle }{R}_{\mathrm{DOC}}=-R_{\mathrm{AR}}-\mathrm{0.8}\times \left(R_{\mathrm{DN}}\right(c_{{\text{s-NO}}_{\mathrm{3}}^{-}})+R_{\mathrm{DN}}(c_{{\text{g-NO}}_{\mathrm{3}}^{-}}\left)\right)\end{array}$$

where V_AR and V_DN [T⁻¹] are the maximum reaction rate of AR and DN, X_AR and X_DN [M L⁻³] are the biomass of functional microbial groups facilitating the reaction components of AR and DN. K_inh [M L⁻³] is a non-competitive inhibition factor used for representing inhibition of DN given oxygen availability.

2.2.2 Boundary condition

The pressure head imposed along the top boundary was described as a sinusoidal function that moves downstream by the ripple migration velocity (u_c [L T⁻¹]; Boano et al., 2013):

$$\begin{array}{}\text{(10)}& {\left.h\right|}_{y=l}=h_m\cdot \sin m(x-u_{\mathrm{c}}\mathrm{d}t)\end{array}$$

where m [–] is the wave number of the variation ($m=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$), dt [T] denotes a timestep in the simulation, and h_m [L] is the amplitude of the head variation (Elliott and Brooks, 1997):

$$\begin{array}{}\text{(11)}& h_m=a{\displaystyle \frac{U^{\mathrm{2}}}{\mathrm{2}g}}{\left({\displaystyle \frac{H_d/H}{\mathrm{0.34}}}\right)}^e\end{array}$$

where a=0.28 [–] is a dimensionless coefficient, U [L T⁻¹] is the average stream velocity, H_d [L] is the height of the ripple, H [L] is the water depth. The exponent e equals to $\mathrm{3}/\mathrm{8}$ if H_d<0.34H and $\mathrm{3}/\mathrm{2}$ otherwise. Given the periodic characteristics typically observed in riverbeds under field conditions (Ahmerkamp et al., 2017; Harvey et al., 2012), a single ripple was chosen in this study. Quasi-periodic pressure and periodic solute boundary conditions were applied to side boundaries AB and CD to replicate the repetitive geometric constraints of the streambed. This is a well-established practice in hyporheic exchange modeling, as it minimizes edge effects without the need for computationally intensive full-scale simulations of multiple consecutive bedforms (Fig. 2):

$$\begin{array}{}\text{(12)}& {\displaystyle }P(\mathrm{0},y,t)=P(\mathit{\lambda },y,t)+\mathrm{\Delta }P\text{(13)}& {\displaystyle }\begin{array}{rl}{c}_i(\mathrm{0},y,t)& =c_i(\mathit{\lambda },y,t){\left.\mathit{n}\cdot \left(-D_{ij}\mathrm{\nabla }{c}_i\right)\right|}_{x=\mathrm{0}}\\ & ={\left.\mathit{n}\cdot \left(D_{ij}\mathrm{\nabla }{c}_i\right)\right|}_{x=\mathit{\lambda }}\end{array}\end{array}$$

where P [M L⁻¹ T⁻²] is the pressure calculated using the hydrostatic head from Eq. (10) and elevation, $P=\mathit{\rho }g(h+y)$. The additional pressure drop ΔP [M L⁻¹ T⁻²] was derived from the streambed gradient and calculated using ΔP=Sρgλ. For the validation of model generalization, simulation of three consecutive ripples was conducted. The middle ripple was focused on to compare its vertical boundary-related pressure and solute concentration with outcomes from the single-ripple model. Simultaneously, comparisons were made between the two models (single-ripple and three-consecutive-ripple) regarding simulated riverbed flow fields, nitrate plumes, influxes, and reaction rates. Although moderate left/right boundary pressure differences were detected between their simulation results, such variations had negligible effects on reactive solute transport. This comparison is provided in Sect. S2 of the Supplement.

At the bottom boundary, an upward groundwater flux was defined to mimic the process of groundwater discharge. For groundwater-borne nitrate, the bottom boundary was specified as an open Dirichlet boundary with constant g-NO${}_{\mathrm{3}}^{-}$ concentration in groundwater:

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{ll}c=c_i& \mathit{n}\cdot \mathit{v}<\mathrm{0}\\ \mathit{n}\cdot \left(-D_{ij}\mathrm{\nabla }{c}_i\right)=\mathrm{0}& \mathit{n}\cdot \mathit{v}\ge \mathrm{0}\end{array}\right.\end{array}$$

where the outward unit normal is n. Besides, the top boundary for the g-NO${}_{\mathrm{3}}^{-}$ was assigned as outflow condition.

For stream-borne solutes includes DOC, O₂, and s-NO${}_{\mathrm{3}}^{-}$, an open boundary with a constant solute concentration was imposed on the top boundary, and outflow boundary was applied on the bottom of the domain. The concentrations of the all solutes were initially assumed to be zero. A quantitative mass-balance check for the model was included in Sect. S2.

Figure 2 Schematic diagram of the simulated domain and boundary conditions.

2.3 Bedform migration

In this study, ripples formation was qualitatively determined using the bedform stability diagram (Fig. 3). This diagram is a summarized reference for examining equilibrium bed configurations in unidirectional flow, derived from a series of flume experiments and field studies (Ashley, 1990). Besides, a series of quantitative criteria (Sect. S1) was examined to ensure ripple formation under the modeled scenarios and the achievement of dynamic equilibrium. For a specific grain size of streambed sediment (D₅₀), particular stream velocities (U) satisfying all these criteria are selected and displayed in Fig. 3. All simulation scenarios were listed in the Table S1 in the Supplement.

Figure 3 The bedform stability diagram (modified from Ashley, 1990) showing the bedform properties and hydraulic conditions considered in this study. Different colors indicate different median grain sizes of the streambed sediments. For each sediment grain-size scenario, various streamflow conditions were simulated, corresponding to the velocity values shown on the y-axis.

Ripple migration velocities were implemented using an empirical relation after Coleman and Melville (1994), which was derived from flume experiments:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{u_{\mathrm{c}}}{\left(u^{\ast }-u_{\mathrm{cr}}^{\ast }\right)\left({\mathit{\tau }}^{\ast }-{\mathit{\tau }}_{\mathrm{cr}}^{\ast }\right)}}{\left(H_d/D_{\mathrm{50}}-\mathrm{3.5}\right)}^{\mathrm{1.3}}=\mathrm{40}\end{array}$$

where u^∗ = (gHS)^0.5 [L T⁻¹] is the bed shear velocity, S [–] is the stream gradient and calculated by Chezy equation ($U=H^{\mathrm{2}/\mathrm{3}}{S}^{\mathrm{1}/\mathrm{2}}/n$, where n [–] is the Manning coefficient and assumed to be 0.02 for sand). $u_{\mathrm{cr}}^{\ast }$ [L T⁻¹] is the critical bed shear velocity and it can be calculated by the critical Shield parameter ${\mathit{\tau }}_{\mathrm{cr}}^{\ast }$ (${\mathit{\tau }}_{\mathrm{cr}}^{\ast }={\mathit{\tau }}_{\mathrm{cr}}/g({\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho })D_{\mathrm{50}}$, $u_{\mathrm{cr}}^{\ast }=({\mathit{\tau }}_{\mathrm{cr}}/\mathit{\rho }{)}^{\mathrm{0.5}}$, ρ_s [M L⁻³] and ρ [M L⁻³] are the density of sediment and water), and τ^∗ is the shield number related to the bed shear velocity (${\mathit{\tau }}^{\ast }=u^{\ast \mathrm{2}}/rg{D}_{\mathrm{50}}$; r [–] is the submerged specific gravity of sediment = (${\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho })/\mathit{\rho }$). The critical shields parameter ${\mathit{\tau }}_{\mathrm{cr}}^{\ast }$ defines the threshold for the initialization of motion, and it is determined by the dimensionless particle parameter D^∗ [–] (Soulsby, 1997; Zheng et al., 2019):

$$\begin{array}{}\text{(16)}& {\mathit{\tau }}_{\mathrm{cr}}^{\ast }={\displaystyle \frac{\mathrm{0.3}}{\mathrm{1}+\mathrm{1.2}{D}^{\ast }}}+\mathrm{0.055}\left[\mathrm{1}-\exp \left(-\mathrm{0.02}{D}^{\ast }\right)\right]\end{array}$$

where $D^{\ast }=D_{\mathrm{50}}\cdot (rg/{\mathit{\upsilon }}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{3}}$, υ is the kinematic viscosity of water [L² T⁻¹]. The derived celerity was substituted into Eq. (10) to determine the migration of sinusoidal pressure head profile with ripples moving.

2.4 Governing non-dimensional numbers

The characteristics of the modeled system were depicted by a series of non-dimensional numbers, which represent the relative dominance of various forces, transport, and reaction processes in this system. Firstly, the Reynolds number Re was used to characterize the flow condition of surface water (Cardenas and Wilson, 2006):

$$\begin{array}{}\text{(17)}& \mathit{Re}={\displaystyle \frac{U{H}_d}{\mathit{\upsilon }}}\end{array}$$

where υ [L² T⁻¹] represents kinematic viscosity of water.

The dimensionless parameter U_r [–] was introduced to represent the relative magnitude of bedform celerity and the pressure-induced pore water velocity driven by pressure variation over the ripple surface and upwelling groundwater:

$$\begin{array}{}\text{(18)}& U_{\mathrm{r}}={\displaystyle \frac{\mathit{\theta }\cdot u_{\mathrm{c}}-u_{\mathrm{s}}}{u_{\mathrm{p}}}}\end{array}$$

where u_s [L T⁻¹] is the seepage velocity of the underflow induced by stream gradient (u_s=KS), and thus the characteristic horizontal velocity is $u_{\mathrm{c}}-u_{\mathrm{s}}/\mathit{\theta }$. $u_{\mathrm{p}}/\mathit{\theta }$ [L T⁻¹] is the pore water velocity induced by pumping process and is calculated using the analytical solution after Boano et al. (2009) and Fox et al. (2014) accounting for vertical groundwater flux (u_q [L T⁻¹]):

$$\begin{array}{}\text{(19)}& {\displaystyle }\begin{array}{rl}{u}_{\mathrm{p}}& =u_{\mathrm{p},\mathrm{0}}\sqrt{\mathrm{1}-{\left(u_{\mathrm{q}}/\mathit{\pi }{u}_{\mathrm{p},\mathrm{0}}\right)}^{\mathrm{2}}}\\ & +\left(\left|u_{\mathrm{q}}\right|/\mathit{\pi }\right){\sin }^{-\mathrm{1}}\left(\left|u_{\mathrm{q}}\right|/\mathit{\pi }{u}_{\mathrm{p},\mathrm{0}}\right)-\left(\left|u_{\mathrm{q}}\right|/\mathrm{2}\right)\end{array}\text{(20)}& {\displaystyle }{u}_{\mathrm{p},\mathrm{0}}=a{\displaystyle \frac{K{U}^{\mathrm{2}}}{g\mathit{\lambda }}}{\left({\displaystyle \frac{H_d/H}{\mathrm{0.34}}}\right)}^e\end{array}$$

where u_p,0 [L T⁻¹] represents the hyporheic exchange solely driven by pressure variation over the sediment-water interface; if U_r>1, turnover process dominates and controls the hyporheic exchange, otherwise, the system is pumping process dominated (Jiang et al., 2022).

The relative magnitude of hyporheic exchange flux driven by pressure variation along the riverbed surface and upwelling groundwater flux was determined as:

$$\begin{array}{}\text{(21)}& U_{\mathrm{b}}={\displaystyle \frac{u_{\mathrm{q}}}{u_{\mathrm{p},\mathrm{0}}}}\end{array}$$

The relative dominance of hyporheic exchange and biogeochemical reaction in nitrate removal can be defined by the Damköhler number (Ocampo et al., 2006; Zarnetske et al., 2012; Zheng et al., 2019):

$$\begin{array}{}\text{(22)}& \mathit{Da}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{T}}}{{\mathit{\tau }}_{\mathrm{R}}}}\end{array}$$

the characteristic timescale for the transport of solutes through the ripple was estimated as (Azizian et al., 2015):

$$\begin{array}{}\text{(23)}& {\mathit{\tau }}_{\mathrm{T}}={\displaystyle \frac{\mathit{\lambda }\mathit{\theta }}{{\mathit{\pi }}^{\mathrm{2}}{u}_{\mathrm{p}}}}\end{array}$$

and the reaction timescale (τ_R) represents the time needed to consume dissolved oxygen of hyporheic water to a prescribed anoxic environment threshold (0.2 mg L⁻¹). The reaction timescale was described as:

$$\begin{array}{}\text{(24)}& {\mathit{\tau }}_{\mathrm{R}}={\displaystyle \frac{\ln \left(c_{{\mathrm{O}}_{\mathrm{2}}}/c_{{\mathrm{O}}_{\mathrm{2}},lim}\right)}{V_{\mathrm{AR}}}}\end{array}$$

biogeochemical reactions are transport-limited when Da < 1. The biogeochemical reactions depended on reaction kinetics due to the brevity of the time that reactants spend within the HZ. Under these low Da conditions, the HZ remains oxic conditions, resulting in a minimal or no denitrification to occur. Conversely, when Da > 1, the residence time of reactants exceeds the reaction time, and thus oxygen is consumed and favors for the occurrence of denitrification in anoxic conditions (Jiang et al., 2022; Zarnetske et al., 2011a).

2.5 Model setup and parametrization

All parameter values in this study were shown in Table 1. The bedform geometry of Ping et al. (2022) was adopted for this study (the ripples located at λ_c=0.15 m with a height of H_d=0.02 m). Here grain sizes D₅₀ of 0.08, 0.15, 0.22 and 0.36 mm were considered, typically falling within characteristic grain diameters on sandy riverbeds (Ahmerkamp et al., 2017; Harvey et al., 2012). The permeability of riverbed was calculated using the empirical relation $k={\mathit{Da}}^{\prime }\times \mathrm{735}\times {\mathrm{10}}^{\mathrm{6}}\times D_{\mathrm{50}}^{\mathrm{2}}$ (where Da^′ = $\mathrm{9.869}\times {\mathrm{10}}^{-\mathrm{13}}$ is the conversion factor for unit Darcy to m²; Gangi, 1985).

The concentrations of DOC, O₂, and s-NO${}_{\mathrm{3}}^{-}$ in stream were specified as 30, 8, and 5 mg L⁻¹. This configuration represents a pristine stream characterized by moderate nutrient levels (Ocampo et al., 2006). The g-NO${}_{\mathrm{3}}^{-}$ in groundwater was set as 15 mg L⁻¹, representing the chemical signature of nitrate-contaminated groundwater that lacks both oxygen and organic matter (Hester et al., 2014). The maximum reaction rate and corresponding functional microbial concentration for AR and DN were listed in Table 1, the choose biogeochemical values are consistent with the parameter setting of nutrient cycling in hyporheic zones and riparian zones (Gu et al., 2008; Nogueira et al., 2021; Zarnetske et al., 2012).

Table 1 Model parameters used in numerical simulations.

^a Janssen et al. (2012). ^b Ping et al. (2022). ^c Zarnetske et al. (2012). ^d Sawyer (2015). ^e Ahmerkamp et al. (2015). ^f Bardini et al. (2012).

The following distinct model experiments were carried out: The Reynolds number, i.e., mean stream velocity, was varied for Re = 2000–6000 in intervals of 500, with corresponding stream water velocities of U = 0.1–0.3 m s⁻¹. The range of U_b was set from 0.3 to 0.7 in intervals of 0.1; correspondingly, the upwelling groundwater flux u_q was ranged between 0.3 × u_p,0 and 0.7 × u_p,0. A larger upward flux than 0.9 × u_p,0 would eliminate the entire hyporheic flow cell, thus the maximum boundary flux was set at a value slightly below this threshold. Meanwhile a minimum of 0.2 × u_p,0 ensures that upwelling groundwater is still mixing with surface water with minor influences on hyporheic flow cell. Longitudinal and transverse dispersivities (α_L and α_T) are important values for obtaining accurate mixing. Here sensitivity analysis of local dispersion was conducted by varying α_L from 0.001 to 0.015 m, with α_T set to $\mathrm{1}/\mathrm{10}$ of α_L. For the base case model, we set α_L equal to 0.001 m and α_T equal to 0.0001 m for groundwater borne conservation tracer. Molecular diffusion was fixed at 10⁻⁹ m² s⁻¹, which has a secondary effect on solute dispersion process. This sensitivity analysis was provided in Sect. S4.

The finite element software, COMSOL Multiphysics (version 6.1) was used to solve the Darcy flow and multi-component solute transport model. Three computational grids were evaluated to ensure simulated results independence of grid resolution, confirming that the grid is sufficiently refined to capture mixing dynamics and minimize numerical dispersion. The fine, base and coarse grid sizes were 1.5, 2, and 2.5 mm, respectively. The grid independence analysis was provided in Sect. S2. For the simulations, the base grid size was selected, with the domain discretized into 19 940 cells. To maintain a constant bedform displacement (Δx) per timestep, the simulation was conducted with Δx=2 mm, while dt was adjusted in inverse proportion to the migration celerity u_c. The total duration of the simulation was set to be equal to the time needed for hundreds of ripples to travel across the modeled domain until the hyporheic exchange and biogeochemical processes reached quasi-steady states.

2.6 Model Metrics

2.6.1 Mixing of surface water and groundwater

Here, macroscopic mixing is defined as the colocation of surface water (SW) and groundwater (GW) within a specific aquifer volume (e.g., a numerical model element or cell), a process that causes solutes to be present simultaneously in an overlapping region (Nogueira et al., 2022). There is no established standard for the threshold governing SW and GW fractions during the mixing process; nonetheless, the general recognition is that mixing occurs when the GW proportion falls within the range of 10 % to 90 % (Hester et al., 2013; Woessner, 2000). This interval effectively distinguishes the occurrence of mixing, where SW and GW interact dynamically, from the two endmembers: pure GW (>90 %) and pure SW (<10 %). To evaluate the influence of threshold variations on mixing metrics, sensitivity analysis was conducted to examine three alternative GW fraction ranges: 10 %–90 % (wider interval), 16 %–84 %, and 20 %–80 % (narrower range). Thus, mixing flux and mixing zone were calculated based on these three thresholds.

A constant concentration boundary condition of conservative tracer (c_gw=1 mg L⁻¹) in GW was specified to represent the tracer signature of groundwater, following the methods outlined by Hester et al. (2013, 2014). The groundwater-borne tracer that undergoes mixing as it transitions from flow paths originating at the bottom boundary to those emerging at the streambed surface was quantified to determine the net effect of mixing. The streambed surface was divided into three zones (Fig. 1): “SW IN”, where surface water enters the riverbed layer; “SW OUT”, where surface water discharges back into the overlying water column; and “GW OUT”, where upwelling groundwater discharges into the stream. The conservation tracer flux for the SW OUT zone was used to describe mixing. If no mixing occurred, all the conservative tracer entering the model at the bottom boundary would exit through the GW OUT zone on the riverbed surface. The arrival of conservative solute at the riverbed surface is defined as the start of mixing, and the complete leave of all the conservation solute from the riverbed is defined as the end of mixing. Due to dispersion, the groundwater-borne tracer gradually spreads through the streambed and eventually crosses the sediment-water interface. The terms “mixing start and end” denote the period over which mixing is calculated.

The mixing flux across the riverbed surface was determined by integrating the outward volumetric flux containing groundwater-borne conservation solute within the specified concentration range (i.e., c_gw = 0.1–0.9 mg L⁻¹, c_gw = 0.16–0.84 mg L⁻¹, c_gw = 0.2–0.8 mg L⁻¹). This flux corresponds to groundwater flowing through “SW OUT” zone. Concurrently, the mixing zone was calculated by integrating the riverbed area that contains groundwater-borne conservation solute within the specified concentration range. The average mixing flux and mixing zone over the whole mixing period were computed and termed the mixing flux (MF) and mixing zone, respectively. The mixing fraction (F_mix) was computed as the ratio of the mixing flux to the total hyporheic flux, which was derived from integrating the volumetric flux along the sediment-water interface (including both SW OUT and GW OUT). F_mix serves as a metric to quantify the proportion of SW-GW mixing flux relative to the total flux of SW and GW interaction across the sediment-water interface. A_mix is defined as the ratio of the mixing zone to the entire modeled domain, representing the proportion of the riverbed occupied by the mixing zone.

2.6.2 Nitrate reaction rate and removal efficiency

When the hydro-physical and biogeochemical conditions reach a quasi-steady state, the average flux of stream- or groundwater- borne nitrate entering the riverbed layer over the final serval periods of ripple migration was calculated:

$$\begin{array}{}\text{(25)}& {\displaystyle }\begin{array}{rl}{M}_{{\text{in,s-NO}}_{\mathrm{3}}^{-}}& ={\displaystyle \frac{\mathrm{1}}{w}}{\sum }_{T}B\int \mathit{v}{c}_{{\text{s-NO}}_{\mathrm{3}}^{-}}\\ & -\mathit{\theta }{D}_{ij}\cdot \mathrm{\nabla }{c}_{{\text{s-NO}}_{\mathrm{3}}^{-}}\mathrm{d}{L}_{\mathrm{top}}\end{array}\text{(26)}& {\displaystyle }\begin{array}{rl}{M}_{{\text{in,g-NO}}_{\mathrm{3}}^{-}}& ={\displaystyle \frac{\mathrm{1}}{w}}{\sum }_{T}B\int \mathit{v}{c}_{{\text{g-NO}}_{\mathrm{3}}^{-}}\\ & -\mathit{\theta }{D}_{ij}\cdot \mathrm{\nabla }{c}_{{\text{g-NO}}_{\mathrm{3}}^{-}}\mathrm{d}{L}_{\mathrm{bottom}}\end{array}\end{array}$$

where B [L] is the per unit width, θ [–] is the sediment porosity, and $w=T/\mathrm{d}t$ [–]; L_top and L_bottom refer to the top boundary BC and bottom boundary AD, respectively.

The total amount of nitrate removed by NMD denitrification and MD denitrification were calculated as follows during the same time interval:

$$\begin{array}{}\text{(27)}& {\displaystyle }{M}_{\mathrm{NMD}}={\displaystyle \frac{\mathrm{1}}{w}}{\sum }_{T}B\int R_{\mathrm{DN}}\left(c_{{\text{s-NO}}_{\mathrm{3}}^{-}}\right)\mathrm{d}A\text{(28)}& {\displaystyle }{M}_{\mathrm{MD}}={\displaystyle \frac{\mathrm{1}}{w}}{\sum }_{T}B\int R_{\mathrm{DN}}\left(c_{{\text{g-NO}}_{\mathrm{3}}^{-}}\right)\mathrm{d}A\end{array}$$

where A [L²] is the area of the streambed.

Nitrate removal efficiency was quantified as the rate of nitrate removed from the riverbed layer divided by the nitrate flux entering the riverbed layer:

$$\begin{array}{}\text{(29)}& N_{\text{RE-NMD}}={\displaystyle \frac{M_{\mathrm{NMD}}}{M_{{\text{in,s-NO}}_{\mathrm{3}}^{-}}}},N_{\text{RE-MD}}={\displaystyle \frac{M_{\mathrm{MD}}}{M_{{\text{in,g-NO}}_{\mathrm{3}}^{-}}}}\end{array}$$

3 Results

3.1 Model validation

3.1.1 Validation of flat bed model with triangular ripple model

Model validation was conducted in two steps, with two numerical methods employed to represent moving bedforms. Beyond simulating porewater flow and reactive solute transport in a fixed frame of reference with a flat riverbed surface, a reference frame moving with triangular ripples was also adopted to capture bedform migration, following the methods of Bottacin-Busolin and Marion (2010), Ahmerkamp et al. (2015), and Kessler et al. (2015). The latter method captures bedform migration while preserving the ripple shape, and comparisons between the two approaches were conducted to examine the impact of neglecting ripple geometry on the simulated results.

A typical case (Re = 3000, U_b=0.6 and D₅₀=0.15 mm) was chosen here for comparison. The relative differences between these two modeling approaches were compared after normalizing the triangular ripple shape to a flat riverbed by projecting the values on the sloped edges onto the horizontal straight boundary. The overall patterns of pressure head, Darcy velocity, and stream-borne solute plume distributions were mostly similar (Fig. 4). The difference in groundwater-borne nitrate plumes is attributed to the distinct advection velocities for reactive solute transport. Specifically, the moving reference frame involves a transformation of the horizontal coordinate system, based on the assumption that the frame of reference moves in the direction of bedform migration and at the same migration celerity; in this frame, the modeled advection velocity of reactive solutes is the seepage velocity minus the celerity. Thus, in the triangular ripple model, g-NO${}_{\mathrm{3}}^{-}$ flows upstream with larger advection velocities, resulting in a more uniform solute plume; in contrast, it exits the riverbed across the entire riverbed surface in the flat bed model (see the white streamlines). Correspondingly, discrepancies emerge in the crest area of the MD denitrification zone.

Figure 4 Comparison of hyporheic flow and reactive solute transport modeling results from the flat bed model and the triangular ripple model. Streamlines and arrows in the Darcy flow field indicate porewater flow, while streamlines in nitrate solute plumes represent the advection velocity of reactive solute transport. Relative differences in hyporheic and biogeochemical parameters between the two models are also shown.

Compared with the triangular ripple model, the flat bed model exhibits larger hyporheic exchange flux, mixing flux, and stream-borne nitrate influx. The relative difference is approximately 5 %. Meanwhile, both the total MD and NMD denitrification rates are higher in the triangular ripple model than in the flat bed model, with relative difference around 10 %. Additionally, the spatial relative differences between these two models are depicted in Fig. 4. For any given parameter, the relative difference is defined as the absolute value of the difference between the two models' calculated results at the identical coordinate location, divided by the maximum value derived from the flat bed model. Note that the largest differences occur in the MD denitrification hotspots. Spatially, the flat bed model captures the hotspots of the denitrifying zone but underestimates the rates, potentially leading to an underestimation of groundwater borne nitrate removal efficiency. To further assess the differences between the two models and their reliability, two methods were utilized to perform validation against the experiments by Wolke et al. (2020).

3.1.2 Validation of numerical models with laboratory experiments

The developed models were validated through comparison with Wolke et al.'s (2020) flume experiments, which investigated the evolution of oxygen in the riverbed under varying conditions of mean stream velocity (0.16–0.32 m s⁻¹) and bedform migration celerity (0–0.394 cm h⁻¹). The experiment was designed with a total of 5 operating conditions, each of which was repeated twice and labeled as Set 1 and Set 2. Oxygen distribution within sediments was measured via planar optodes, and the oxygenated zone was defined as regions where oxygen saturation exceeds 15 %. Oxygen consumption was observed within riverbed layer and aerobic respiration is widely recognized as the dominant process driving oxygen consumption within riverbed sediments (Ahmerkamp et al., 2017; Reeder et al., 2018). In their experiments, a fixed amount of NaCl was added to surface water; this NaCl was then used to compute hyporheic exchange flux (HEF) from the decline in its concentration, following the method by Fox et al. (2018). The oxygen influx (F_in) and outflux (F_out) were calculated by multiplying HEF by the surface water oxygen concentration and the mean oxygen concentration in the oxygenated zone, respectively.

Two numerical models: the flat bed model and the triangular ripple model, were developed under the same hydraulic conditions, with their parameters adjusted to fit HEF, F_in, F_out, and A_O to experimental data. The model parameters were summarized in Sect. S3. All physical and hydraulic parameters (e.g., ripple geometry, streamflow, bedform migration velocity, porosity, and permeability) were used exactly as reported in Wolke et al. (2020) experiments without any adjustment. Only the biogeochemical reaction parameters (e.g., maximum oxygen consumption rate, and half-saturation constant for oxygen) were calibrated to match the observed oxygen plume, as these parameters are difficult to measure directly in dynamic streambed environments and are commonly calibrated in reactive transport models.

Figure 5 Comparison of numerical modeled oxygen distribution, (e) hyporheic exchange flux, (f) summed oxic zone and (g) oxygen influx and (h) oxygen outflux and experimental measurements by Wolke et al. (2020) under various conditions. The units of RMSE are (e) cm d⁻¹, (f) cm², (g) mmol m⁻² d⁻¹, and (h) mmol m⁻² d⁻¹. RMSE values are computed by comparing simulated results with averaged experimental data under identical hydraulic conditions.

In the numerical model, HEF was calculated by integrating the volumetric flux over the inflow zones along the sediment-water interface. Oxygen transport was governed by the advection-dispersion-reaction equation, where oxygen consumption occurred through aerobic respiration following the Monod kinetics described by Eq. (8). The F_in and F_out were then computed by multiplying the HEF by the oxygen concentration in the surface water and the oxygen concentration over the outflow zones along the riverbed surface, respectively. The oxygenated area was calculated by integrating the riverbed portions where oxygen saturation exceeds 15 %. Based on the criteria for ripple migration, it was determined that under the hydraulic conditions of Run 5, ripples could not migrate while maintaining their shapes due to increased stream flow velocity. Therefore, the model validation simulations considered four hydraulic conditions of Run 1 to 4.

For stationary and slow-migrating beds, the spatial distribution of oxygen creates a conchoidally shaped plume beneath the bedform; in contrast, beneath fast-migrating beds, the oxygen plume shifts to a more uniform front. The HEF, F_in, F_out, and A_O simulated using both the ripple bedform model and the flat bed model were presented in Fig. 5; these results are in good qualitative and quantitive agreement with the measurements of Wolke et al. (2020). The root means square error (RMSE) between the experimentally measured values and those simulated by the flat bed model were further calculated. All RMSE values fall within reasonable ranges, indicating the critical processes of bedform migration and oxygen dynamic are captured by the numerical models. This comparison with Wolke et al. (2020) validates the model's capacity to reproduce bedform migration. Conversely, for mixing dynamics and mixing-triggered reactions, most prior studies depend solely on numerical simulations, with no experimental data available to validate these mixing processes under gaining conditions.

3.2 Effect of bedform migration on mixing and solute dynamic

To simulate the range of natural environmental conditions, the reactive transport equations were solved for different stream velocities, grain sizes and groundwater upwelling fluxes, which include a corresponding range of ripple migration celerites and sediment permeabilities derived from the empirical relations. As an example, the patterns of pore water transport as well as SW and GW mixing are shown for a grain size of 0.15 mm, a constant ratio between pumping driven hyporheic exchange flux and upwelling GW flux U_b=0.6, and five different stream velocities, that is five different Re numbers.

For low surface water flow velocity (Re = 2500), no migration of bedform is predicted by the model. SW enters the riverbed layer in the high-pressure region on the stoss side, flows through the porous medium, and exits the riverbed layer in the low-pressure region on the lee side, forming a typical conchoidally shaped hyporheic flow cell. Upwelling GW is diverted around the hyporheic flow cell, mixes with SW, and then exits into the overlying water from both sides, in patterns similar to those shown previously by Fox et al. (2014) and Hester et al. (2019). SW and GW mixing zone (i.e. the mixing area where the fraction of GW ranges between 16 % and 84 %) emerges as a thin band along the hyporheic flow cell, and covers over 2.97 % of the modeled domain. Meanwhile, NMD denitrification occurs below the oxygen plume with the reactive zone in a conchoidal shaped distribution, while MD denitrification reactive zone develops along the edge of the mixing zone, where DOC from SW meets nitrate from GW (row 1 in Fig. 6).

Figure 6 Effect of bedform migration on riverbed biogeochemistry for U_b=0.6 and D₅₀=0.15 mm. Shown are profiles of (column 1) conservation solute representing groundwater fraction (c_gw), (column 2) DOC concentration (c_DOC), (column 3) O₂ concentration (c_DO), (column 4) s-NO${}_{\mathrm{3}}^{-}$ concentration ($c_{{\text{s-NO}}_{\mathrm{3}}^{-}}$), (column 5) g-NO${}_{\mathrm{3}}^{-}$ concentration ($c_{{\text{g-NO}}_{\mathrm{3}}^{-}}$), (column 6) non-mixing-dependent (NMD) denitrification rate ($R_{{\text{s-NO}}_{\mathrm{3}}^{-}}$) and (column 7) mixing-dependent (MD) denitrification rate ($R_{{\text{g-NO}}_{\mathrm{3}}^{-}}$).

At a surface water Re of 3000, bedforms migrate at a low celerity of 1.04 cm h⁻¹. This low migration rate has minimal impact on the plumes of stream-borne O₂ and s-NO${}_{\mathrm{3}}^{-}$. The outflow zone of groundwater-borne solutes (c_gw and g-NO${}_{\mathrm{3}}^{-}$) exhibits a slight upstream shift, approaching the middle of the riverbed surface. The size of SW and GW mixing zone (GW faction: 16 %–84 %) increases to 10.06 % of the modeled domain. Furthermore, both the NMD and MD denitrifying rates are enhanced at the start of bedform migration (row 2 in Fig. 6).

With increasing stream velocity (Re = 3500), bedform migration drives continuous changes in pressure distribution. The simulated hyporheic flow cells move downstream and decrease in size at the same time. The shape of the SW and GW mixing zone changes distinctly, forming a horizontal band with a wider range 10.71 % of the whole domain (GW faction: 16 %–84 %). The penetration depths of stream-derived solutes into the streambed are reduced, with a more gradual concentration gradient in both horizontal and vertical directions, whereas the g-NO${}_{\mathrm{3}}^{-}$ plume is uniformly distributed horizontally. Both NMD and MD denitrification hot spots form in the central position near the sediment-water interface as the bedform moving (row 3 in Fig. 6).

When bedform migration is further increased (Re = 4500 and 5000), bedform migration fully dominates over the pore water flow, and hence, continuous solute layers are found in the subsurface (as depicted in row 4 and 5 of Fig. 6). The penetration depths of stream-borne solutes are decreased in comparison to those in slow- to medium- fast migrating bedforms. The NMD and MD denitrification zones become thin and move upward with decreased reaction rates. Similar to the conclusions obtained in previous studies (Kessler et al., 2015; Zheng et al., 2019), bedform migration reduces the penetration depths of stream-borne solutes and the scope of hyporheic exchange cell. We also found that bedform migration increases the size of mixing zone between SW and GW. The SW-GW mixing zone accounts for 9.30 % and 9.14 % of the domain area, respectively (take the GW fraction 16 %–84 % for example).

3.3 Effect of migration celerity on mixing regimes

The mixing intensity across the riverbed surface and the size of the mixing area within the subsurface were estimated by simulating four different grain sizes and five upwelling GW fluxes, under varying stream velocities associated with the corresponding bedform celerity. The evolutions of the hyporheic exchange flux, net mixing flux and size of mixing zone with increasing Re number were summarized in Fig. 7. The 16 %–84 % GW fraction range was displayed here to evaluate how migration celerity affects mixing regimes.

As shown in Fig. 7, across various grain sizes, both HEF and MF increase with increasing stream velocity and bedform celerity. Nevertheless, the mixing fraction exhibits a significant increase at the onset of bedform migration, with a tendency to stabilize or decrease moderately as celerity further increases. This phenomenon is attributed to enhanced mixing between SW and discharging GW along the entire horizontal riverbed in moving bedforms – contrasting with stationary bedforms, where mixing occurs solely along the typical crescent-shaped hyporheic exchange flow cells (Fig. 6). Additionally, when bedforms start moving and move at moderate velocities, the SW-GW mixing zone shows a significant increase compared to stationary ones. As bedform migration celerity accelerates, the extent of the mixing zone gradually approaches a plateau, with slight reductions observed at higher celerities.

Figure 7 The variation of hyporheic exchange flux (HEF, row 1) and mixing flux (MF, row 2) across the sediment-water interface, the proportion of mixing flux to hyporheic exchange flux (F_mix, row 3) and the ratio of mixing zone to whole domain (A_mix, row 4) with Re across different medium grain size when U_b = 0.3–0.7.

Although the GW fraction range determines both mixing flux and scope, it does not alter the effect of migration celerity on these mixing patterns (Figs. S4 and S5 in the Supplement). A broader groundwater fraction range induces a larger mixing flux, an increased F_mix, and a wider mixing zone; nonetheless, A_mix persists as a relatively small scope (about 10 %–15 % of the full domain) within the riverbed. Relatively moderate differences in computed MF and F_mix are observed between the narrower range (20 %–80 %) and the baseline range (16 %–84 %). For bedforms of different sediment grain sizes, MF represents around 60 % of the total HEF when bedforms are moderately to fast-moving, compared with approximately 20 % when they are stationary or slow-moving.

Taking a grain size of 0.15 mm and a constant ratio between pumping-driven hyporheic exchange flux and upwelling groundwater flux (U_b=0.6) as an example, we varied α_L from 0.001 to 0.015 m. While increased longitudinal and transverse dispersivities significantly enhance F_mix and A_mix, they do not modify the influence of bedform migration on these mixing metrics. This confirms that greater dispersivities facilitate SW and GW mixing, leading to a larger proportion of hyporheic exchange flux involved in mixing and an expanded mixing zone. As such, reliable dispersion parameter estimates are essential for accurate mixing magnitude assessments. Additionally, F_mix and A_mix values obtained under various GW fraction ranges and longitudinal and transverse dispersivities were presented in Sect. S4.

3.4 Impact of ripples migration on nitrate removal

To assess the impact of ripple migration on the removal of s-NO${}_{\mathrm{3}}^{-}$ and g-NO${}_{\mathrm{3}}^{-}$ within domains of varying medium grain sizes, the influx of nitrate into the riverbed layer and the total reaction rate within the riverbed layer were determined.

Figure 8 The non-mixing-dependent denitrification rates (R_NMD) as functions of U_b and Re for different medium grain sizes.

For s-NO${}_{\mathrm{3}}^{-}$, the NMD denitrification rate increases in moving bedforms composed of very fine sand (D₅₀=0.08 mm; Fig. 8). This is likely caused by higher flow velocities driving more reactants into the sediment, thereby enhancing denitrification in reaction-limited systems (Da > 5.04). In contrast, in riverbeds of fine to medium sand (D₅₀ = 0.15–0.36 mm) with higher permeabilities, the rate of s-NO${}_{\mathrm{3}}^{-}$ is negatively correlated with the mean stream velocity when the system becomes transport-limited. This is likely because s-NO${}_{\mathrm{3}}^{-}$ travels fast along shorten flow paths and does not undergo denitrification within the moving bedforms (Fig. 6). Additionally, the migrating bedforms enhance the delivery of s-NO${}_{\mathrm{3}}^{-}$ into the sediment due to increased HEF. Consequently, the removal efficiency of s-NO${}_{\mathrm{3}}^{-}$ decreases monotonically across various medium grain sizes (Fig. 9).

Figure 9 The removal efficiencies of stream borne nitrate (N_RE-NMD) as functions of U_b and Re for different medium grain sizes.

For g-NO${}_{\mathrm{3}}^{-}$, the increase in MD denitrification is also seen for fine sandy sediments (D₅₀ = 0.08–0.22 mm) at low to medium Re (Fig. 10). Compared to s-NO${}_{\mathrm{3}}^{-}$, the advective flow paths and residence times of g-NO${}_{\mathrm{3}}^{-}$ in groundwater are longer. Additionally, the enhanced mixing dynamics between SW and GW further facilitates the occurrence of MD denitrification as the bedforms start moving. Consequently, the reduction rate of g-NO${}_{\mathrm{3}}^{-}$ decreases only in riverbeds consisting of medium sand with fast bedform migration celerity (Re > 4000), as the solute residence time is significantly reduced. Interestingly, the rise in the MD denitrification rate compensates for the increased g-NO${}_{\mathrm{3}}^{-}$ influx in very fine sand (D₅₀ = 0.08 and 0.15 mm) at moderate stream velocity (Re < 4000). Most of g-NO${}_{\mathrm{3}}^{-}$ that enters the sediment is consumed before entering the overlying water column. For fine to medium sand riverbed, the g-NO${}_{\mathrm{3}}^{-}$ removal efficiency decreases strongly with increasing Re. The natural protective role of the SW and GW mixing zone in preventing nitrate-contaminated GW from entering rivers is being hindered in fast moving bedforms (Fig. 11).

Figure 10 The mixing-dependent denitrification rates (R_MD) as functions of U_b and Re for different medium grain sizes.

Figure 11 The removal efficiencies of groundwater borne nitrate (N_RE-MD) as functions of U_b and Re for different medium grain sizes.

4 Discussion

This study quantified the effect of bedform migration on SW and GW mixing process as well as mixing triggered denitrification. Previous research has primarily focused on the potential impacts of bedform migration on hyporheic exchange driven by streambed morphological features, as well as non-mixing-dependent biogeochemical processes where reactants are assumed to be predominantly in SW. However, such studies represent only a small subset of possible streambed environmental conditions, focusing exclusively on specific headwater stream conditions (Jiang et al., 2022; Kessler et al., 2015; Ping et al., 2022; Zheng et al., 2019). The impact of bedform migration on the conceptual model of bedform-induced hyporheic exchange, which is influenced by GW upwelling and/or ambient lateral GW flow in the mid-stream section of lowland rivers, has received relatively less attention and examination.

In streams and rivers that are fed by regional GW and possess undulating bedforms, SW gets mixed with GW throughout the local hyporheic exchange process. The mixing zone exhibits a typical crescent shape along the periphery of typical hyporheic exchange cells within a stationary streambed (as reported by Fox et al., 2014; Hester et al., 2019; Nogueira et al., 2022). The sizes of the SW-GW mixing zone (e.g., thickness and area) occupy a small proportion of the whole HZ. In the immobile bedform, the thin mixing zone occupying ∼ 5 % is consistent with prior work (Hester et al., 2013; Santizo et al., 2020). However, at the onset of bedform migration, the mixing pattern, extent, and intensity of SW-GW interactions are altered. A continuous SW-GW mixing zone is formed within the ripples of the medium- to fast-moving bedform (Fig. 6), and the area of mixing zone increases to approximately 10 %–15 % at this time. Besides, the net flux of SW and GW mixing increases with stream velocity and bedform migration celerity significantly (Fig. 7). As a result, bedform migration controls and determines the hotspots and magnitudes of the SW and GW mixing. The bedforms are typically assumed to be immobile potentially making underestimations of SW-GW mixing flux and mixing zone in a HZ.

Instead of the typical crescent-shaped MD denitrifying zone observed in stationary bedforms (Naranjo et al., 2015; Hester et al., 2014, 2019), the MD reaction zonation changes to the layer shape distributed at the fringe of the HEF cells, where mixing between SW and GW develop to a largest degree in fast-moving bedforms (Fig. 6). Such a situation exists where the stream flows into with a relative homogeneous sandy riverbed with low autochthonous organic carbon content and encounters with nitrate enrich GW. The heterogeneous streambeds including buried autochthonous organic matter (Sawyer, 2015; Ping et al., 2022), deposited particulate organic particles (Drummond et al., 2017; Ping et al., 2023), and biological aggregate (Xian et al., 2022), would complicate the hyporheic exchange process, induce the rough and irregular shapes and boundaries of HEF cells and mixing zone, and therefore affect redox microenvironments and biogeochemical zonations. MD denitrifying hotspot would also develop around available sources of DOC.

Figure 12 Variation of the dimensionless Damköhler number as a function of the dimensionless parameter U_r.

Previous studies have demonstrated that migrating bedforms constrains the penetration depths of stream-borne solutes, shorten their transport timescales, and reduces the removal efficiencies of stream-borne nitrate (Jiang et al., 2022; Kessler et al., 2015; Ping et al., 2022; Zheng et al., 2019). Our findings indicate that bedform migration exerts more complex influences on MD denitrification and groundwater-borne nitrate removal. Bedform migration facilitates increases in MD denitrification rate, driven by an enhanced mixing flux, a higher mixing proportion, and an expanded mixing zone. At the same time, it shortens the residence timescales of water and solutes, which is detrimental to the occurrence of denitrification. When Re < 4000 (slow to moderate stream velocities), the MD denitrification rate of moving bedforms is greater than or comparable to that of stationary bedforms; in contrast, for fast-moving bedforms (Re ≥ 4000), the MD denitrification rate decreases significantly. This is likely because the reduction in MD rate driven by shortened transport timescales dominates over the increase caused by enhanced mixing intensity; meanwhile, at this stage (∼ Re 4000), the mixing fraction of the total exchange flux and the extent of the mixing zone both approach plateaus. Fast-moving bedforms also reduce the HZ attenuation capacity for groundwater-borne nitrate. Notably, the g-NO${}_{\mathrm{3}}^{-}$ removal efficiency in HZs peaks for stationary bedforms and those at the onset of migration.

The causes of fast-moving bedforms' negative impact on g-NO${}_{\mathrm{3}}^{-}$ removal differ depending on the type of riverbed sediment (Fig. 12). For very fine sand (D₅₀=0.08 mm), a larger celerity results in a decline in transport timescale with less impact on removal efficiency as the system is rate-limited (Da ≫ 1). For fine to medium sand riverbed (D₅₀ = 0.15–0.36 mm), the transport-limited situation leads to a low denitrification rate with increasing U_r but constantly decreasing Da (Da < 1). It is important to note that the calculated reaction timescale is based on the consumption period of oxygen to a prescribed anoxic threshold. Under these conditions, the exhaustion of labile DOC would also lead to the cessation of denitrification (Zarnetske et al., 2011a, b). These results demonstrate that in order to evaluate the self-purification capacity of the HZ and its function as a natural barrier mitigating groundwater contamination, riverbed sediment transport dynamics and grain size distributions need to be considered. Stabilizing bedform configurations in restoration projects would enhance the natural attenuation capacity of HZs. Additionally, enhancing mixing intensity, along with extending the water residence timescale, would also facilitate the self-purification of HZs.

Different from previous studies that demonstrated the magnitude of NMD denitrification is often greater than that of MD denitrification (Hester et al., 2014; Trauth and Fleckenstein, 2017). The results in our study show that the total reaction rate of s-NO${}_{\mathrm{3}}^{-}$ is smaller than that of g-NO${}_{\mathrm{3}}^{-}$ in mobile bedforms. This phenomenon can be attributed to the following two reasons: first, the concentration of s-NO${}_{\mathrm{3}}^{-}$ is one-third that of g-NO${}_{\mathrm{3}}^{-}$; second, the reaction zone is reduced by migration celerity for NMD denitrification, while the reaction zone for MD denitrification is increased in moving bedforms. Mixing intensity increases with bedform migration, a process that facilitates MD denitrification more effectively (Hester et al., 2019; Nogueira et al., 2024; Trauth and Fleckenstein, 2017). More attention should be paid to the mixing dynamics and mixing triggered biogeochemical reactions, which is helpful to put forward appropriate stream restoration plans so as to enhance the health of the aquatic ecosystem (Hester et al., 2017; Lawrence et al., 2013).

The findings of this study, derived from 2D numerical models, are most directly generalizable to straight, low-curvature streams with periodic bedform distributions in riverbed sediments. In such stream systems, hyporheic exchange is dominated by streamwise-vertical flow cells, with lateral (cross-stream) hyporheic flux accounting for a small fraction of total exchange (Hu et al., 2014; Naranjo et al., 2015). Importantly, straight, low-curvature streams are prevalent in agricultural and urban downstream gaining reaches, which are characterized by relatively homogeneous sediment types (e.g., fine sands). This is consistent with the focus of our study, as these reaches often face nitrate pollution challenges and depend on natural attenuation to maintain water quality (Hester et al., 2014; Trauth et al., 2017). In practice, riverbeds may exhibit far more complex three-dimensional (3D) bedform morphologies coupled with sediment heterogeneity, for example, in highly meandering streams (Gomez-Velez et al., 2017; Nogueira et al., 2024; Pescimoro et al., 2019). Future studies should establish 3D numerical model and incorporate stochastic hydraulic conductivity (K) fields to explore how sediment heterogeneity interacts with bedform migration, evaluating whether high-K hotspots enhance or reduce migration-driven mixing.

A critical limitation of this study is the reliance on empirical equations for bedform migration celerity. The Chezy equation (suited for flat, homogeneous sediment beds) was used to estimate the shear velocity for incipient sediment motion, with the experimentally derived empirical relationship by Coleman and Melville (1994) further adopted to account for the additional shear stress induced by ripples. The two-step approach (shear stress for particle initiation and ripple-scale form drag via empirical formulations) reproduced Wolke et al.'s (2020) results and proved useful for investigating the effects of bedform migration on hyporheic solute transport (Ahmerkamp et al., 2017; Zheng et al., 2019). We did not simulate turbulent flow over triangular bedforms nor compute shear stress via computational fluid dynamics (CFD; a numerical method that would yield more precise pressure profiles and shear stress distributions) given that the core focus of this study was to investigate how bedform migration affects hyporheic solute transport and mixing-dependent denitrification within the riverbed layer. As future work, the empirical approach should be validated using CFD simulations to enhance the accuracy of bedform-related stress and migration calculations. This validation will facilitate the quantification of potential uncertainties and refine the mechanistic understanding of bedform dynamic-hyporheic zone interactions.

This study focuses on ripples and, more broadly, shorter-wavelength topographic roughness elements formed under low subcritical flow conditions in sandy riverbeds (Ashley, 1990; Gomez-Velez et al., 2015; Raudkivi, 1997). The undulating bedforms maintain dynamic equilibrium through geometric adjustments, with their geometry remaining unchanged as the stream velocity fluctuates within a specific range (10–30 cm s⁻¹). When stream velocities exceed this upper threshold, a condition commonly observed in fast-flowing rivers, bedform geometries can be altered, ultimately leading to bedform erosion (Boano et al., 2013; Harvey et al., 2012). This process is not accounted for in the current model. If small-scale ripples develop and merge with larger-scale ripples and dunes under high stream velocities, the removal efficiency of s-NO${}_{\mathrm{3}}^{-}$ and g-NO${}_{\mathrm{3}}^{-}$ may be enhanced due to the extended hyporheic flow paths and increased residence timescales (Harvey et al., 2012; Zomer and Hoitink, 2024). Otherwise, the removal of s-NO${}_{\mathrm{3}}^{-}$ would be highly hindered because of shorter residence time and fully oxic condition in fast moving bedforms and fast flowing rivers, while the removal of g-NO${}_{\mathrm{3}}^{-}$ would likely be less affected within the immobile streambeds.

The numerical models assume isotropic sediment K. In fact, natural riverbeds with small-scale bedforms (e.g., dunes and ripples) commonly exhibit strong anisotropy. For bedforms with flow-transverse crests, pore pathways parallel to the crests are relatively continuous and straight yielding high transverse K. In contrast, pathways perpendicular to the crests are highly tortuous: water here must move upward through erosional, high-porosity troughs and downward over depositional, low-porosity crests. This undulating flow path, combined with flow separation and local energy losses at morphological transitions, substantially increases hydraulic resistance and lowers longitudinal K (Dallmann et al., 2020; Salehin et al., 2004). Additionally, spatial variability in bedform dimensions (e.g., increasing dune wavelength downstream) enhances this anisotropy by introducing zones of differing flow resistance along the streamwise direction (Venditti et al., 2005). Given that the model is constructed along a longitudinal section with spatially consistent, periodic bedforms, the use of an isotropic hydraulic conductivity assumption is justified. Incorporating anisotropic K is critical for accurately modeling hyporheic exchange and related biogeochemical processes in 3D heterogeneous models.

DNRA was not incorporated in the model, given that denitrification is typically regarded as the predominant pathway for nitrate removal, whereas DNRA plays a secondary role in nitrate transformation (Zarnetske et al., 2012). Lansdown et al. (2012) and Quick et al. (2016) have demonstrated that approximately 5 % of ¹⁵NO${}_{\mathrm{3}}^{-}$ tracer in river sediment incubations underwent DNRA, while 85 % underwent denitrification. Nevertheless, DNRA competes with denitrification for NO${}_{\mathrm{3}}^{-}$ and DOC as electron acceptor and donors within HZs. When an oligotrophic and/or a pristine stream infiltrate into the streambed and subsequently interact and mix with nitrate-enriched GW, the MD DNRA would not occur due to the low $\mathrm{C}/\mathrm{N}$ ratio. When a eutrophication stream with higher DOC concentration, DNRA would have a greater influence in nitrate transformation because DNRA is prone to occur in NO${}_{\mathrm{3}}^{-}$-limited (that is DOC sufficient) conditions compared to denitrification (Zhu et al., 2023). The ammonia produced by MD DNRA would be further nitrified within the aerobic HEF cell, thereby potentially elevating the risk of nitrate pollution in SW.

In this model, stream velocity and upward GW flux are considered constant in the present model, yet they may change in time due to storm events, tidal pumping, snowmelt, or reservoir hydro-peaking (Liu et al., 2024; Nogueira et al., 2022; Song et al., 2018). Hester et al. (2019) demonstrated that increasing surface water stage would enhance both NMD and MD denitrification. Nogueira et al. (2024) and Trauth and Fleckenstein (2017) pointed out that GW discharge events increase the magnitude of SW-GW mixing, therefore affecting the prevalence of MD denitrification. The interactions among morphological dynamics, hyporheic exchange, and biogeochemical processes under transient conditions are key areas for future research.

5 Conclusion

The numerical model developed in this study was used to simulate the interaction and mixing of upwelling groundwater with bedform-induced hyporheic flow, examining how bedform migration influences SW and GW mixing and the processing of groundwater-borne nitrate within the HZ. Our analysis quantified the mixing flux and the size of mixing zone, as well as the mixing-dependent denitrification rates and removal efficiencies across riverbed sediments characterized by varying grain sizes, stream flow velocities, and groundwater discharge fluxes. These model simulations revealed that as bedforms migrate, the SW-GW mixing zone and the associated mixing-dependent denitrification zone progressively evolve into uniform, band-like structures. When turnover dominates the hyporheic exchange process, the mixing flux increases significantly (an increase of one order of magnitude). The mixing proportion of total exchange flux and mixing zone scope also experience significant increases at the onset of migration; as stream velocity further rises, the mixing proportion (∼ 60 %) and mixing size (∼ 10 %–15 % of the riverbed) gradually reach plateaus. Enhanced SW-GW mixing dynamics facilitates MD denitrification as bedforms begins migration; however, shortened solute residence timescales resulting from fast bedform movement limits MD denitrification occurrence. Under dynamic bedform conditions, the self-purification capacity of the HZ is reduced for fine to medium sand (D₅₀ = 0.15–0.36 mm), compromising its role as a natural barrier against groundwater contamination. Incorporating the identified and analyzed factors can enhance the management of riverbed sediment-associated aquatic systems, particularly for management objectives focused on removing groundwater-borne nitrate.

Appendix A

S	Stream slope [–]
H	Water depth [L]
U	Stream velocity [L T−1]
λ	Ripple wavelength [L]
uc	Bedform migration celerity [L T−1]
l	Streambed height [L]
Hd	Ripple height [L]
B	Riverbed per unit width [L]
D50	Median sediment size [L]
x	Horizontal coordinate, rightward positive [–]
y	Vertical coordinate, upward positive [–]
k	Sediment permeability [L2]
K	Sediment hydraulic conductivity [L T−1]
h	Hydraulic head [L]
hm	Amplitude of the sinusoidal head variation [L]
m	Wavenumber of the variation [–]
g	Gravity acceleration [L T−2]
ci	Concentration of reactive components [M L−3]
cgw	Groundwater tracer [M L−3]
v	Seepage velocity [L T−1]
vi	Seepage velocity in the i-direction [L T−1]
vj	Seepage velocity in the j-direction [L T−1]
θ	Sediment porosity [–]
αL	Longitudinal dispersivity [L]
αT	Transverse dispersivity [L]
Dij	Hydrodynamic dispersion [L2 T−1]
Dm	Molecular diffusion coefficient [L2 T−1]
ι	Tortuosity factor [–]
VAR	Maximum reaction rate of aerobic respiration [T−1]
VDN	Maximum reaction rate of denitrification [T−1]
Kinh	Non-competitive inhibition factor [M L−3]
KDOC	Half-saturation for dissolved organic carbon [M L−3]
$K_{{\mathrm{NO}}_{\mathrm{3}}^{-}}$	Half-saturation for nitrate [M L−3]
$K_{{\mathrm{O}}_{\mathrm{2}}}$	Half-saturation for oxygen [M L−3]
XAR	Microbial concentration facilitating aerobic respiration [M L−3]
XDN	Microbial concentration facilitating denitrification [M L−3]
ρs	Sediment density [M L−3]
ρ	Water density [M L−3]
P	Hydraulic pressure [M L−1 T−2]
D∗	Dimensionless particle parameter [–]
u∗	Bed shear velocity [L T−1]
$u_{\mathrm{cr}}^{\ast }$	Critical bed shear velocity [L T−1]
n	Manning coefficient [–]
τ∗	Shield parameter [–]
τcr	Critical shear stress [M L−1 T−2]
${\mathit{\tau }}_{\mathrm{cr}}^{\ast }$	Critical Shield parameter [–]
r	Submerged specific gravity of sediment [–]
$R_{{\mathrm{O}}_{\mathrm{2}}}$	Aerobic respiration rate [M L−3 T−1]
$R_{{\text{s-NO}}_{\mathrm{3}}^{-}}$	Non-mixing-dependent denitrification rate [M L−3 T−1]
$R_{{\text{g-NO}}_{\mathrm{3}}^{-}}$	Mixing-dependent denitrification rate [M L−3 T−1]
RDOC	Dissolved oxygen carbon consumption rate [M L−3 T−1]

υ	Kinematic viscosity of water [L2 T−1]
μ	Dynamic viscosity of water [M L−1 T−1]
us	Underflow seepage velocity induced by stream gradient [L T−1]
Re	Reynolds number [–]
Ur	ratio of bedform celerity to pore water velocity [–]
up	Darcy velocity induced by pumping process [L T−1]
uq	Vertical groundwater flux [L T−1]
Ub	Ratio of vertical groundwater flux to hyporheic exchange flux [–]
τR	Biogeochemical reaction timescale [T]
τT	Water transport timescale [T]
Da	Damköhler number [–]
Δx	Bedform migrating displacement per timestep [L]
dt	The length of per timestep [T]
Da′	Conversion factor for unit Darcy to m2 [–]
Fmix	Proportion of mixing flux to hyporheic exchange flux [–]
Amix	Proportion of mixing zone to the whole domain [–]
A	Streambed area [L2]
MNDN	Nitrate removed by non-mixing-dependent denitrification [M T−1]
MDN	Nitrate removed by mixing-dependent denitrification [M T−1]
NRE	Nitrate removal efficiency [–]
Min	Nitrate being introduced into streambed [M T−1]
Fin	Oxygen influx [M L−2 T−1]
Fout	Oxygen outflux [M L−2 T−1]
AO	Oxygenated area [L2]
Abbreviations
HZ	Hyporheic zone
SW	Surface water
GW	Groundwater
HEF	Hyporheic exchange flow
MF	Mixing flux
NMD	Non-mixing-dependent
MD	Mixing-dependent
AR	Aerobic respiration
DN	Denitrification
RMSE	Root means square error
DNRA	Dissimilatory nitrate reduction to ammonium
ANAMMOX	Anaerobic ammonium oxidation