the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 29 Jul 2025
Modeling surface water and groundwater mixing and mixing-dependent denitrification with bedform dynamics
Abstract. 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. Pioneering 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 remain approximately constant when turnover becomes the dominant water exchange mechanism for fine-medium sandy riverbed. Fast stream flows and migrating bedforms reduce solute residence timescales and limits denitrification opportunities. Consequently, nitrate removal efficiency from both stream- and groundwater-borne sources decreases significantly with bedform migration in fine-medium sandy sediments. 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.
-
Notice on discussion status
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
-
Preprint
(1763 KB)
-
Supplement
(237 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(1763 KB) - Metadata XML
-
Supplement
(237 KB) - BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-2631', Anonymous Referee #1, 20 Aug 2025
-
AC1: 'Reply on RC1', Zhang Wen, 29 Aug 2025
Thank you for your valuable and insightful comments, which are greatly helpful in improving the clarity and overall quality of our manuscript. We will carefully address all the comments if we get a chance to revise the manuscript. Please find our briefly response to your comments in the following:
Regarding the parameter explanations: We will conduct a systematic review of all parameters in the manuscript, supplementing and expanding their definitions, including physical meanings and units, to improve readability. In addition, for the writing style: All instances of subjective phrasing (e.g., “we”) will be revised into objective expressions in accordance with standard academic writing conventions.
For the specific comments raised by this reviewer, we totally agree and we will accept all of your suggestions during the revision process of the paper.
Once again, thank you very much for your time spending on review our manuscript. If you have any other questions, please feel free to contact us.
Citation: https://doi.org/10.5194/egusphere-2025-2631-AC1
-
AC1: 'Reply on RC1', Zhang Wen, 29 Aug 2025
-
RC2: 'Comment on egusphere-2025-2631', Anonymous Referee #2, 31 Aug 2025
Ping X. et al reported a modeling study of the mixing process and mixing-dependent denitrification beneath a moving bed form. The topic is very important, and results are valuable. However, there are several critical issues in conceptual model description, ripple geometry simplification, nitrate concentration portioning, shear stress estimation, and model validation. Major issues were summarized below for the authors for consideration. Detailed comments are attached at the end.
Major comments:
- Simplifying ripples as a flat surface need justification. For non-mobile bedform, it is well known that non-uniform pressure distribution applied on the ripple surface cause pressure gradients horizontally and vertically, which further control the complex pressure, velocity, and concentration distributions inside and beneath the bedform. Due to its driving force, it is critical to include the ripple geometry in simulations, which is a common practice in most modeling research for stationary hyporheic zone (HZ) (see Cardenas2007, WRR). For mobile bedform, one of the few experiments on hyporheic zone with mobile bedforms show that the distribution of concentrations is mostly complex inside and near moving ripple (Wolke2020, Water). This also suggests the most important role of ripple geometry in reproducing realistic flow and transport distributions in HZ. Without representing ripple geometry in the model, the exchange flux, concentration distributions, and shear stress on, inside, and beneath the ripple will likely have high uncertainty if they are estimated from existing empirical equations. I suggest the authors justify such assumption by adding a new simulation with ripple geometry and then evaluate how much difference could be in exchange flux, concentration flux, and shear stress from models with and without ripple.
- The conceptual model presentation and Figure 1 needs improvements to better align with text descriptions, Figure 1 and texts in Section 2.2-2.6. In the current modeling approach for HZ, popular approaches are two-layer and three-layer conceptual model. Figure 1 suggests a two-layer model, i.e., splitting the domain as a surface water layer (regions above ripple top boundary) and groundwater layer (regions below ripple top boundary) with riverbed surface (ripple top boundary) as a permeable interface. However, when reading the details in Section 2.1-2.6, the model is a three-layer structure: surface water layer, riverbed layer (a thick 3D layer instead of a riverbed surface), and groundwater layer. Different from traditional 3-layer structure, this paper does not explicitly solve Navier-Stockes equations for surface water layer and Darcy-equation for groundwater layer. Insteady, this work only solves the pressure head and concentration for the riverbed layer with simplified boundary conditions for pressure head/concentrations at the top (boundary BC in Figure 1) and bottom (boundary AD in Figure 1) boundary of the riverbed layer. Further, this paper split the riverbed layer into two pools, i.e., non-mixing dependent reaction pool (blue region in Figure 1) and mixing dependent reaction pool (red region in Figure 1). Note that such “two” pools share the same computational domain (bounded by ABCD). They do not have an internal interface (critical for traditional two-layer or three-layer models) to conserve water and concentration mass. Insteady, this work defines two variables (c_{s-NO3-} and c_{g-NO3-}) and apply the advection-diffusion-reaction equation (equation 4) twice for these two variables. This means the model is essentially a one-domain multiple-variable model. However, most of this information was not clearly presented in Section 2.1, which is supposed to serve for this purpose. I suggest the authors modify Section 2.1 and relevant texts in Sections 2.2-2.7 to reflect the key conceptual model setups in Section 2.1. Please be aware that further readers may not be familiar with different conceptual models for modeling HZ. Your precise introduction is critical for them to correctly understand your work.
- It is not clear if mass conversation for pressure head and concentration can be achieved based on the current model descriptions. The boundary BC and DC is fine due to use simple “periodic” boundary. For boundary BC and AD, it is not clear to me how to guarantee consistent at boundary conditions at BC and AD for pressure head, head gradient (used to further compute exchange flux), and concentration gradient due to unclear boundary condition descriptions.
- Partitioning the total Nitrate into surface water layer contribution (c_{s-NO3-}) and groundwater layer contribution (c_{g-NO3-}) maybe problematic due to the involvement of non-linear reaction rate in the Monod kinetics. From equation 7, the reaction rate of denitrification R_{DN} is a function of Nitrate concentration, i.e., R_{DN} = f(NO3). If simulating the total NO3- with portioning, there will be two reaction rate R_{DN}(c_{s-NO3-) (equation 8) and R_{DN}(c_{g-NO3-) (equation 9). Due to non-linearity of Monod equation f, R_{DN} does not equal the summation of the reaction rate from two partition components, i.e., f(NO3) ≠ f(sNO3) + f(gNO3). This means partition approach is not mass-conservation. Additionally, introducing two concentration variables means need 4 more boundary conditions at top and boundary conditions for c_{s-NO3-} and c_{g-NO3-}. This paper does not clearly explain how such boundary conditions are assigned. In short, the current Nitrate concentration approach causes mass conservation issues and needs 4 more boundary conditions data which are not well described. I suggest the authors clarify these issues in the revision.
- Estimation of shear stress that drives sand transport is likely not reliable. The estimation of migration velocity uc depends on shear velocity and shear stress on the ripple surface. Due to the simplification of ripple into a flat surface, calculation of shear stress/shear velocity is not possible. This paper estimates the shear velocity using the Cheezy equation. The Cheezy equation usually works well for flat surface with homogenous sediments. In this work, the sand is very small (0.08-0.36 mm) and homogenous. However, the ripple is much bigger (height 0.02 m). The correct shear stress is the one on the surface with both large ripple and small sand. The use of Cheezy is likely not adequate to estimate such variables. The most accurate way will be using numerical model to directly compute shear stress with the ripple geometry as the top boundary. If the authors want to use the Cheezy equation, it needs to be verified against at least one CFD simulated or observed shear stress to quantify potential uncertainty.
- What were compared in the model validation is not clear. Figure 5 shows a comparison between model and observation. However, how HEF, AO, Fin, Fout were computed is not explained clearly. The model results are distributed data. The authors need to explain how to convert distributed values into single values like HEF, AO, Fin, Fout. Second, the oxygen concentration from the model is affected by reaction R1. It is not clear if the experiment also implemented a similar reaction. If the model has a reaction R1 and the experiment does not, then it is not reasonable to compare model predicted RO with experiment observed RO. Please add more information about the details of model validation. Especially, please clarify if the model setup is consistent with the experiment setup, especially the reaction part.
- The validation result from Figure 5 suggests model prediction for HEF and concentration flux may be reliable (if potential issues in Comment 5 get clarified). However, simulation for concentration distribution (which could be used to compute AO in Figure 5b) shows much larger uncertainty. Because this work focuses on analyzing denitrification that depends on accurate prediction of nitrate concentration. I am concerned that the prediction for Nitrate may have similar uncertainty level as the RO shown in Figure 5b. I suggest the authors provide additional information to verify if the prediction uncertainty for nitrate is negligible for your analysis. If lab experiment is not available. A potential solution is to run your model with the ripple geometry explicitly represented by mesh and compute the difference of Nitration distribution with and without ripple geometry. According to Line 349, the authors seem to acknowledge that the larger error in RO shown in Figure 5b could be caused by the exclusion of the ripple geometry. Such addition modeling work will help to explain Figure 5b and enhance the prediction reliability for nitrate distribution.
- This work has several figures (e.g., figures 1-4) and most equations similar to the Ping X. 2022 (10.1029/2022WR033258), please highlight major novelty in this paper compared to the existing paper.
Other minor comments are listed below for reference:
Line 22: missing a comma between “migration” and “this”.
Line 123: Riverbed and riverbed surface mean different things. If the “Riverbed” at this sentence means “the top surface of the riverbed”, it is better to use “riverbed surface”.
Figure 1: Figure 1 looks very similar to Figure 2 in Ping2022 (10.1029/2022WR033258). The ripple shapes are the same in two papers. I guess this figure is modified from Figure 2 in Ping 2022. In this scenario, please add the following information in figure caption: [Figure is modified from Figure 2 in Ping 2022, permission has granted from Wiely]
Figure 1: Because HZ is a complex system, it involves concepts such as domain, layer, surface, interface, flux on surface, in from domain/surface, out from domain/surface. I suggest the authors modifying the Figure 1 and captions to carefully choose vocabularies to “concisely” describe such concepts. For example, riverbed could mean riverbed surface or riverbed layer. Clearly selecting words will help readability of the paper. This comment is also related to Major comment 2.
Line 144: I am not quite agree with this. From Elliot1997 paper, the size of ripple heights is 1-2.5 cm, and surface water depth is of 3-7 cm, which means ripple heights are 15-20% of the water depth. In Cardenas2007 CFD simulation paper (10.1016/j.advwatres.2006.06.009), dune height is 5 cm and water depth is 45 cm. Their ratio is 10% (Figure 2 Cardenas2007). At line 144, I don't see the actual size for ripple and depth, so it is reasonable to say ripple heights are negligible. Such negligible assumption will lead to the simplification of "triangular ripple” into a "flat" surface. From both Elliot1997 and Cardenas 2007 paper, they do not simplify the ripple/dune into flat surface. Instead, Elliot1997 shows the pressure varies sinuously along the surface of the ripple. In Cardenas 2007, it shows that pressure head variations are affected by Reynolds number. The sinusoidal function reported in Elliot1997 is valid only at low Reynolds number and such match is applied only half of the region along the ripple surface (see Figure 9 in Cardenas 2007 10.1016/j.advwatres.2006.06.009). Therefore, whether you can simplify the ripple as flat surface (which will affect how you design the computational domain) and using a sinusoidal function (depending on your Reynolds number) requires more justification. This comment is also related to Major comment 1.
Line 142: show lambda on Figure 1. Also please visualize ripple height in Figure 1.
Line 146: A periodic morphological feature does not mean periodic conditions for pressure, velocity, and concentration. Previous work shows that the flow in upstream of ripple influences flow downstream (see Sinha2017 WRR, 10.1002/2016WR019662). Due to this issue, numerical simulation will need to use at least two ripples to minimize impacts of upstream (see Cardenas2007, 10.1016/j.advwatres.2006.06.009). In this work, a period boundary is applied at the location B and C (see Figure 1) where pressure and velocity are not periodic (see Figure 2 in Cardenas2007). However, all the information mentioned above is for stationary ripple/dune. For the moving ripple, whether you should/should not use a periodic boundary may be different but requires further justification.
Line 153: Is this "BC" the top edge of the blue region in Figure 1? In this work, the 'BC' is a flat surface. However, equation 2 is usually applied to the top surface of the ripple which is not a flat surface. Please justify if such simplification affects results (e.g., exchange flux and distribution of nitrate in the computational domain).
Line 154: equation 2 was first introduced in Boano2013. It is better to cite correct references.
Line 163: For equation 4/5, you use bold (v), non-bold v, vi, vj, and |v|. In general, bold (v) means vector. For non-bold v, there is no definition. Also, no definition for vi, vj. Please correct such variables to make sure they are consistent and mathematically correct.
Line 164 equation 4: gradient v.ci should be gradient dot (v*ci).
Line 173: Does this sentence mean you solved both c_{sNO3} and c_{gNO3}? Note that most previous work does not distinguish which portion of nitrate from surface water and which portion is from groundwater, which means they only solve for one concentration for NO3. From this sentence, it seems you are solving both concentrations of Nitrate. Please confirm if this is the case. Also, for models that solve for one concentration for NO3, you only need to assign two boundary condition for NO3, one is at the top boundary (BC) and one is for the lower boundary (AD). But if you are solving for two concentrations for NO3, you need to assign 4 boundary conditions. You need boundary conditions at both BC and AD for both c_{sNO3} and c_{gNO3}. Please clarify how each of these boundary conditions is determined.
Line 176 equation R1/R2： it is better to add a citation for equations R1 and R2. It is also important to add a justification for why you choose such a reaction. Note that in lab or field experiments, the reaction network is quite complex, which is much different than the equation R1/R2 used here. R1/R2 is a kind of idealized or conceptualized reaction, which could be used as qualitative understating. However, to compare the results to experiment or field observation, you may need to add more information regarding how the assumed reaction R1/R2 is consistent with lab/field reaction network.
Line 180 equation 7: I am not quite sure the meaning of C_{s/g-NO3-}. Does this mean the summation of C_{s-NO3-} and C_{s-NO3-}? Or it means C_{s-NO3-} or C_{s-NO3-}? I also have doubts about the mass balance issue caused by this equation. Please see my major comment 3 for details.
Line 190: the period boundary ab AB and DC needs more justification. From previous work (Cardenas2007 Figure 2 10.1016/j.advwatres.2006.06.009), the distribution of pressure at AB and DC could be quite different depending on Reynolds number. A better way is to use two ripples and apply period at the beginning of the first ripple and the end of the second ripple. Please justify the current setup.
Line 190: If P(0) = P(lamba) + delta P, this is not called period boundary. You can just use the mathematical formula to describe it, but do not call it "periodic". Also, equation 1 means you are solving for pressure head h, but no pressure P. You need to assign the boundary for h, but not P.
Line 193: How do you determine the concentrations at boundary BC? For pressure head h, it is determined by the sinusoidal function from previous work, but for concentration, how do you get it.
Line 195: What is the boundary condition for h at AD boundary?
Line 213: using the Chezy equation to estimate shear stress seems not reliable. Please see my major comment 4 for details.
Line 215: How do you calculate the tau_{cr}?
Line 217: what is r in the equation used for computing shields number at line 217?
Line 231: Is the S at line 231 the same S at line 212? Do you define K somewhere? Do you provide the value for K?
Line 241: Is uq upwelling groundwater flux? How do you decide this value?
Line 264: Is the Da at line 264 the same as the Da in equation 18? If not, please use a different variable.
Line 266: How are these values estimated? Are these assumed values? Please provide a justification for why you choose such values.
Table 1: Please visualize l, lc, lambda, Hd, H in Figure 1. What is the difference between ripple crest and height of ripple crest. In Figure 1, the ripple height is clearly visualized, but the ripple crest is not visualized. Please describe their difference.
Tabel 1: Provide citation for porosity, longitudinal and transverse dispersity or any justification to use such values.
Line 283: Does COMSOL solve the equations 1-11? Do you develop your own COMSOL solver using equations 1-11 or you use the default equations in COMSOL? Is the COMSOL setup and data available in public data repository?
Line 291: I am confused by the term "surface water". From Figure 1, your computational domain is the region bounded by ABCD. The whole region is beneath the surface (or stream) water. This means the model framework does not solve anything related to surface water that is not governed by Darcy flow. Please clarify this.
Line 301: Is bottom boundary AD and streambed surface BC?
Line 302: Is “streambed surface” the "surface water layer" (a layer with pure water)?
Line 303: more accurate description: surface water enters riverbed layer.
Line 304-305: From Figure 1, the groundwater is separated from the surface water layer by "the sediment layer", it seems groundwater cannot be directly transported to the surface water layer. When you use "groundwater discharged to stream", are you trying to identify which portion of water is from groundwater when such groundwater is first "up-welled" to the sediment layer and further "up-welled" to the surface water layer? If this is the case, you may need to explicitly describe this purpose. And you may add a sentence that "such identification could be achieved using tracer or different color/size of tracer".
Line 317: Where is the ripple surface? Is ripple surface the boundary BC? How do you compute Fswi?
Line 322: Is this streambed "boundary BC"?
Figure 3b: From Wolke2019, the oxic zone occurs mainly inside the dune (the triangular region as shown in Figure 1). However, in this paper, the computational domain is a rectangular shape, and the triangular shape is "neglected". How do you compare the oxic zone as reported in the Wolke2019 paper? From Figure 3b, the uncertainty for AO between model and experiment is much larger than the rest 3 variables, could this larger error be caused by the over-simplification of the ripple to a flat surface? Also see major comment 6.
Figure 3: Is this the average upwelling flux at the boundary BC? If not, please clarify where and how the HEF is computed. From the numerical, the flux is a distribution, it is necessary to clarify how you get a single value from a distribution. See also major comment 5.
Line 345: From Figure 3b, the relative error between model and experiments could be unto 100%. The statement of "slightly lower" is not consistent with the data. Please compute the relative error between model results and observation and then make a more accurate conclusion for the accuracy using the computed data.
Line 349: This is related to my major comment 1 and 6. As suggested by the author, the missing of the mobile (triangular ripple) can cause significant error for estimating the total area of the oxygen zone. As the rest of the paper is to study the denitrification which will use data of concentrations of nitrate. My guess is that the missing of the triangular ripple zone (mobile zone) will also significantly affect the accuracy of nitration concentration prediction (could be similar uncertainty of oxygen as shown in Figure 3b). I suggest the authors clarify how much impact will be by simplifying the mobile ripple as flat surface.
Line 375 and Figure 4: This is something uncleaned at lower left corner of Figure 4.
Line 474: This paper has a lot of similarity to Ping 2022 (10.1029/2022WR033258). It may be over statement to use “first time” due to the existing of Xue2022 paper.
Line 852: citation year is not correct.
Citation: https://doi.org/10.5194/egusphere-2025-2631-RC2 -
AC2: 'Reply on RC2', Zhang Wen, 17 Sep 2025
Thank you for your valuable and insightful comments, which are really helpful in improving the clarity and overall quality of our manuscript. We agree with all your concerns, and we will carefully address all the comments if we get a chance to revise the manuscript. Please find our briefly response to your major comments in the following:
1. Simplifying ripples as a flat surface need justification. For non-mobile bedform, it is well known that non-uniform pressure distribution applied on the ripple surface cause pressure gradients horizontally and vertically, which further control the complex pressure, velocity, and concentration distributions inside and beneath the bedform. Due to its driving force, it is critical to include the ripple geometry in simulations, which is a common practice in most modeling research for stationary hyporheic zone (HZ) (see Cardenas2007, WRR). For mobile bedform, one of the few experiments on hyporheic zone with mobile bedforms show that the distribution of concentrations is mostly complex inside and near moving ripple (Wolke2020, Water). This also suggests the most important role of ripple geometry in reproducing realistic flow and transport distributions in HZ. Without representing ripple geometry in the model, the exchange flux, concentration distributions, and shear stress on, inside, and beneath the ripple will likely have high uncertainty if they are estimated from existing empirical equations. I suggest the authors justify such assumption by adding a new simulation with ripple geometry and then evaluate how much difference could be in exchange flux, concentration flux, and shear stress from models with and without ripple.
Response to comment #1: We acknowledge the importance of ripple geometry in HZ modeling and address your concern by supplementing with additional simulations and literature support. Actually, simplifying undulating riverbeds to a horizontal configuration and applying a sinusoidal head distribution at the water-sediment interface exerts minimal impacts on both the hyporheic exchange flux (HEF) driven by the riverbed and solute transport within the riverbed, which has been proven by many researchers (Elliott & Brooks, 1997a, b; Eylers, 1994; Rutherford et al., 1995). Indeed, we have conducted an additional simulation incorporating ripple geometry and compared its results with those of the flat-bed model as the reviewer suggested; the two sets of results exhibit similar characteristics in hyporheic exchange flux and solute plume distributions (This simulation will be added in the revised version). Furthermore, to reproduce the experiment by Wolke (2020), we developed two models: a flat-bed model and a triangular bedform model. Based on the parameter adjustments in the previous manuscript text, we refined the numerical model to achieve a better fit to the experimental data. Its adequacy is further supported by our successful reproduction of Wolke’s (2020) results, thereby confirming its credibility for addressing the study’s objectives. Therefore, in the revised manuscript, we are planning to add two critical contexts: first, a comparison between the triangular bedform model and the flat-bed model; second, model validation against Wolke (2020) using both of these models.
2. The conceptual model presentation and Figure 1 needs improvements to better align with text descriptions, Figure 1 and texts in Section 2.2-2.6. In the current modeling approach for HZ, popular approaches are two-layer and three-layer conceptual model. Figure 1 suggests a two-layer model, i.e., splitting the domain as a surface water layer (regions above ripple top boundary) and groundwater layer (regions below ripple top boundary) with riverbed surface (ripple top boundary) as a permeable interface. However, when reading the details in Section 2.1-2.6, the model is a three-layer structure: surface water layer, riverbed layer (a thick 3D layer instead of a riverbed surface), and groundwater layer. Different from traditional 3-layer structure, this paper does not explicitly solve Navier-Stockes equations for surface water layer and Darcy-equation for groundwater layer. Insteady, this work only solves the pressure head and concentration for the riverbed layer with simplified boundary conditions for pressure head/concentrations at the top (boundary BC in Figure 1) and bottom (boundary AD in Figure 1) boundary of the riverbed layer. Further, this paper split the riverbed layer into two pools, i.e., non-mixing dependent reaction pool (blue region in Figure 1) and mixing dependent reaction pool (red region in Figure 1). Note that such “two” pools share the same computational domain (bounded by ABCD). They do not have an internal interface (critical for traditional two-layer or three-layer models) to conserve water and concentration mass. Insteady, this work defines two variables (c_{s-NO3-} and c_{g-NO3-}) and apply the advection-diffusion-reaction equation (equation 4) twice for these two variables. This means the model is essentially a one-domain multiple-variable model. However, most of this information was not clearly presented in Section 2.1, which is supposed to serve for this purpose. I suggest the authors modify Section 2.1 and relevant texts in Sections 2.2-2.7 to reflect the key conceptual model setups in Section 2.1. Please be aware that further readers may not be familiar with different conceptual models for modeling HZ. Your precise introduction is critical for them to correctly understand your work.
Response to comment #2: We apologize for the confusion and acknowledge that the conceptual model was not clearly explained in the submitted version. It is notable that this model is a one-domain, multiple-variable model which focus on porewater dynamics and multicomponent solute transport within riverbed sediments; its domain corresponds to the actual rectangular riverbed geometry (bounded by ABCD) with no further subdivision of the riverbed interior into two pools. Instead of simulating turbulent flow over bedforms, we used an idealized scenario where sinusoidal pressure variations are applied to a flat-bed (boundary BC in Figure 1). We regret any misunderstanding caused by Figure 1 and will revise the domain schematic to accurately reflect the conceptual model’s configurations. We will overhaul Section 2 to serve as the “conceptual model foundation”: it will explicitly state the model type, define the modeled domain, clarify key processes and assumptions (e.g., no surface water flow simulation, idealized top boundary pressure), and detail the two nitrate variables (origins and boundary conditions). Additionally, we will systematically cross-check and revise Sections 2.2-2.7 to ensure fully consistency with Section 2.1.
3. It is not clear if mass conversation for pressure head and concentration can be achieved based on the current model descriptions. The boundary BC and DC is fine due to use simple “periodic” boundary. For boundary BC and AD, it is not clear to me how to guarantee consistent at boundary conditions at BC and AD for pressure head, head gradient (used to further compute exchange flux), and concentration gradient due to unclear boundary condition descriptions.
Response to comment #3: First, we would like to explain the boundaries in detail. The top boundary BC is prescribed with a Dirichlet head condition (exchange flux is not prescribed, but computed via qex (x, t) = -n·K▽h|BC). The bottom boundary AD is set as Neumann boundary with specified upward flux: -n·K▽h|AD = q (x) (where q > 0 denotes upwelling groundwater into the modeled domain). These two boundaries are mutually compatible, one is prescribing head and the other one is normal flux, ensuring a well-posed mixed boundary setup. For reactive solute transport, the top boundary serves as an open boundary for stream-borne solutes, while the bottom boundary is set as an outflow condition. In contrast, these two boundary conditions are swapped for groundwater-borne nitrate. Open boundary where the flow direction switches steadily between inflow and outflow in accordance with the computed normal Darcy flux qn = n·q: c = cext (x, t) (Dirichlet) if qn < 0 (inflow) and n·(D ▽ c) = 0 (zero diffusive flux) if qn ≥ 0 (outflow). This standard inflow-outflow treatment prevents artificial sources/sinks at outflow and is conservative by construction. In addition, we verified fluid mass balance by calculating the total influx and outflow fluxes across all boundaries over simulated period for Re = 3000, Ub = 0.6, and D50 = 0.15 mm as an example. We have quantified the fluid influx and outflux, as well as the solute input and output, along with the internal consumption. Notably, the relative error remains less than the acceptable threshold (5%) for hydrological models, which effectively confirms the conservation of fluid and solute mass. We will clarify the boundary naming, state the governing equations, and give the exact boundary conditions used for both porewater flow and reactive solute transport, and add the results of mass-balance checks in the revised manuscript.
4. Partitioning the total Nitrate into surface water layer contribution (c_{s-NO3-}) and groundwater layer contribution (c_{g-NO3-}) maybe problematic due to the involvement of non-linear reaction rate in the Monod kinetics. From equation 7, the reaction rate of denitrification R_{DN} is a function of Nitrate concentration, i.e., R_{DN} = f(NO3). If simulating the total NO3- with portioning, there will be two reaction rate R_{DN}(c_{s-NO3-) (equation 8) and R_{DN}(c_{g-NO3-) (equation 9). Due to non-linearity of Monod equation f, R_{DN} does not equal the summation of the reaction rate from two partition components, i.e., f(NO3) ≠ f(sNO3) + f(gNO3). This means partition approach is not mass-conservation. Additionally, introducing two concentration variables means need 4 more boundary conditions at top and boundary conditions for c_{s-NO3-} and c_{g-NO3-}. This paper does not clearly explain how such boundary conditions are assigned. In short, the current Nitrate concentration approach causes mass conservation issues and needs 4 more boundary conditions data which are not well described. I suggest the authors clarify these issues in the revision.
Response to comment #4: We defined stream-borne nitrate and groundwater-borne nitrate as two distinct reactants. Rather than partitioning the nitrate pool into two separate fractions, we solve for the transport of these two nitrate types in the porous medium separately using the advection-dispersion-reaction equation. As outlined in our response to Comment #3, each nitrate type (differentiated by origin) has its own set of boundary conditions. We quantified the input and output of stream-borne and groundwater-borne nitrate, as well as their internal consumption via denitrification. Notably, the relative error remains below the acceptable threshold (less than 5%), confirming the conservation of solute mass. We will provide the exact boundary conditions used for both porewater flow and reactive solute transport, and add the results of mass-balance verification in the revised manuscript.
5. Estimation of shear stress that drives sand transport is likely not reliable. The estimation of migration velocity uc depends on shear velocity and shear stress on the ripple surface. Due to the simplification of ripple into a flat surface, calculation of shear stress/shear velocity is not possible. This paper estimates the shear velocity using the Cheezy equation. The Cheezy equation usually works well for flat surface with homogenous sediments. In this work, the sand is very small (0.08-0.36 mm) and homogenous. However, the ripple is much bigger (height 0.02 m). The correct shear stress is the one on the surface with both large ripple and small sand. The use of Cheezy is likely not adequate to estimate such variables. The most accurate way will be using numerical model to directly compute shear stress with the ripple geometry as the top boundary. If the authors want to use the Cheezy equation, it needs to be verified against at least one CFD simulated or observed shear stress to quantify potential uncertainty.
Response to comment #5: We acknowledge the limitations of the classical Chezy equation in capturing ripple-induced shear stress. In our study, it was only used to estimate baseline shear velocity for incipient sediment motion. To account for bedform development, we incorporated Coleman and Melville’s (1994) empirical relationship for ripple migration celerity, indirectly considering ripple-associated form drag. This two-step approach (grain-scale resistance for particle initiation and empirical ripple-scale form drag) is widely adopted and validated, and successfully reproduced Wolke’s (2020) results, supporting its credibility for our objectives. Besides, our primary focus was on the effects of bedform migration on solute transport within the riverbed. Rather than simulating turbulent flow over triangular bedforms and calculating shear stress via CFD, we employed two approaches: an idealized model in which sinusoidal pressure variations are applied to the sediment layer surface, and empirical equations to determine bedform migration celerity. We fully agree that CFD simulations would yield more precise estimates of pressure profile and shear stress, but this lies beyond the scope of our current study. In the revised manuscript. We will explicitly state the limitations of our study, propose CFD validation as a valuable direction for future work, and confirm that the current empirical approach is appropriate for the study’s context and objectives.
6. What were compared in the model validation is not clear. Figure 5 shows a comparison between model and observation. However, how HEF, AO, Fin, Fout were computed is not explained clearly. The model results are distributed data. The authors need to explain how to convert distributed values into single values like HEF, AO, Fin, Fout. Second, the oxygen concentration from the model is affected by reaction R1. It is not clear if the experiment also implemented a similar reaction. If the model has a reaction R1 and the experiment does not, then it is not reasonable to compare model predicted RO with experiment observed RO. Please add more information about the details of model validation. Especially, please clarify if the model setup is consistent with the experiment setup, especially the reaction part.
Response to comment #6: We regret the lack of clarity in the model validation description and will clarify this in both this response and the revised manuscript. Specifically, we will detail the experimental measurement methods and simulation-based calculation approaches for the key parameters: hyporheic exchange flux (HEF), oxygen influx/outflux (Fin/Fout), and oxygenated area (AO). The experimentally observed RO is calculated from oxygen influx/outflux and mean oxygen penetration depth, and represents the volumetric oxygen consumption rate within the whole domain. The good fit between simulated and measured Fin/Fout, together with the spatial distribution of oxygen plumes, will ensure a minor discrepancy in this RO-derived rate. Additionally, RO differs in meaning from R1 within the model: the latter specifically refers to the reaction rate per unit pore volume. In the modeling, oxygen transport was governed by the advection-dispersion-reaction equation, where oxygen consumption occurred through aerobic respiration following the Monod kinetic R1. Therefore, we did not compare the model-predicted R1 with the experimentally observed RO. We will provide a clear introduction to and description of the model setup and the corresponding parameter calculation methods in the revised manuscript.
7. The validation result from Figure 5 suggests model prediction for HEF and concentration flux may be reliable (if potential issues in Comment 5 get clarified). However, simulation for concentration distribution (which could be used to compute AO in Figure 5b) shows much larger uncertainty. Because this work focuses on analyzing denitrification that depends on accurate prediction of nitrate concentration. I am concerned that the prediction for Nitrate may have similar uncertainty level as the RO shown in Figure 5b. I suggest the authors provide additional information to verify if the prediction uncertainty for nitrate is negligible for your analysis. If lab experiment is not available. A potential solution is to run your model with the ripple geometry explicitly represented by mesh and compute the difference of Nitration distribution with and without ripple geometry. According to Line 349, the authors seem to acknowledge that the larger error in RO shown in Figure 5b could be caused by the exclusion of the ripple geometry. Such addition modeling work will help to explain Figure 5b and enhance the prediction reliability for nitrate distribution.
Response to comment #7: Two issues that previously reduced the model’s credibility have been addressed: first, the unclear description of model validation (including experimental measurement methods and simulation-based calculation approaches); second, the poor fit between simulated results and experimental data. Specifically, we have refined the numerical parameters to improve the fit of AO, Fin, and Fout, and further optimized the numerical model based on the parameter adjustments outlined in the original manuscript, achieving better alignment with the experimental data. As suggested by the reviewer, we additionally conducted a simulation incorporating ripple geometry (intended to reproduce both Wolke’s (2020) experiment and our model settings) and compared its results with those of the flat-bed model. Therefore, in the revised manuscript, we are planning to add two critical sections of context: first, a detailed comparison between the triangular bedform model and the flat-bed model; second, a comprehensive model validation against Wolke’s (2020) data using both of these models.
8. This work has several figures (e.g., figures 1-4) and most equations similar to the Ping X. 2022 (10.1029/2022WR033258), please highlight major novelty in this paper compared to the existing paper.
Response to comment #8: This is an extension of Ping et al. (2022), but it is quite different from the previous publication. The main novelties of this paper can be summarized as: 1) we focus on the downstream “gaining conditions” and first simulate surface water-groundwater mixing under dynamic bedforms; 2) beyond analyzing the variations in hyporheic exchange flux and cell size, we also assess how the surface water-groundwater mixing flux and mixing zone size change with bedform migration; 3) unlike the previous study, which only considers stream-borne nitrate (a limited scope), we quantify both the reduction rate and efficiency of stream-borne and groundwater-borne nitrate. While basic figures/equations overlap (for fundamental hyporheic processes), these innovations fill mixing-reaction coupling gaps and improve understanding of hyporheic nitrate biogeochemistry. We will add this context to the revised manuscript and clearly articulate the novel contributions of this study in the revised Introduction.
For the specific comments raised by this reviewer, we totally agree and we will accept all of your suggestions during the revision process of the paper.
Once again, thank you very much for your time spending on review our manuscript. If you have any other questions, please feel free to contact us.
Citation: https://doi.org/10.5194/egusphere-2025-2631-AC2
-
RC3: 'Comment on egusphere-2025-2631', Anonymous Referee #3, 19 Sep 2025
This study investigates the effects of migrating bedforms on surface-groundwater mixing and denitrification in a gaining stream reach. The authors employ a numerical model that couples Darcy flow with a translating sinusoidal head boundary and a multi-species reactive transport model. While the paper introduces new metrics such as “mixing flux,” “mixing fraction,” and “mixing-zone area,” my primary concern is the substantial overlap with the authors' previous work (2022). The modeling framework, numerical implementation, and validation appear nearly identical.
The claimed novelty rests on applying the existing model to a different scenario and introducing new diagnostic metrics. However, the justification for these metrics is insufficient, and their introduction does not lead to the discovery of new physical mechanisms. Furthermore, the main conclusion: bedform migration shortens residence times and thereby reduces denitrification—reiterates, is a key finding from the 2022 paper (https://doi.org/10.1029/2022WR033258). The manuscript reads more as an incremental application of prior work rather than a standalone scientific article with a distinct and significant contribution.
1.The literature review should be expanded to include a more thorough discussion of numerical model development in this field. Critically, this section must clearly and explicitly state how this study's model and objectives differ from and advance upon the 2022 paper. The work appears to be a case study applying a previously established model. It does not introduce a new mechanistic framework, a dimensionless parameter space map, or any analytical/scaling relationships that would constitute an independent and generalizable scientific contribution.2.Thresholds used to define the mixing zone (e.g., 16%–84% concentration range) are presented without a clear physical or statistical basis. The robustness of the results is questionable without a sensitivity analysis of these thresholds, as well as the numerical dispersion and grid resolution. Furthermore, the newly proposed metrics are not validated against any established, well-accepted measures of mixing from the fluid dynamics or hydrology literature, making their utility and interpretation difficult to assess.
3.The model's idealizations limit the generalizability of its findings. Key simplifications include: (a) representing bedform migration as a simple translating sinusoidal head, which neglects morphodynamic feedbacks; (b) using a 2-D domain with a single bedform wavelength and amplitude; and (c) omitting sediment heterogeneity, such as grain size variations or layering, which are known to strongly influence hyporheic exchange.
Citation: https://doi.org/10.5194/egusphere-2025-2631-RC3 -
AC3: 'Reply on RC3', Zhang Wen, 25 Sep 2025
General comments: This study investigates the effects of migrating bedforms on surface-groundwater mixing and denitrification in a gaining stream reach. The authors employ a numerical model that couples Darcy flow with a translating sinusoidal head boundary and a multi-species reactive transport model. While the paper introduces new metrics such as “mixing flux,” “mixing fraction,” and “mixing-zone area,” my primary concern is the substantial overlap with the authors' previous work (2022). The modeling framework, numerical implementation, and validation appear nearly identical. The claimed novelty rests on applying the existing model to a different scenario and introducing new diagnostic metrics. However, the justification for these metrics is insufficient, and their introduction does not lead to the discovery of new physical mechanisms. Furthermore, the main conclusion: bedform migration shortens residence times and thereby reduces denitrification—reiterates, is a key finding from the 2022 paper (https://doi.org/10.1029/2022WR033258). The manuscript reads more as an incremental application of prior work rather than a standalone scientific article with a distinct and significant contribution. Response to the general comments: We thank this reviewer for the constructive comments on our manuscript. We think the main concern of this reviewer is the difference between this study and our previous publication (Ping et al., 2022). We will classify the difference and the novelty of this study clearly in the revised version. Actually, this study is quite different from Ping et al. (2022) because: 1) Firstly, in this study, we focus on downstream section of rivers under gaining conditions firstly to simulate surface water-groundwater mixing processes under dynamic bedforms, while Ping et al. (2022) focused on hyporheic exchange processes in the headwater to midstream sections of rivers. therefore, the research focus is quite different for these two studies. 2) Secondly, the mixing process has been considered in this study. As far as we know, this process has not been thoroughly investigated in the migrating ripple model (Jiang et al., 2022; Zheng et al., 2019; Kessler et al., 2015). It is notable that the bedform migration affect mixing process significantly. 3) Thirdly, Ping et al. (2022) only considered stream-borne solutes and non-mixing dependent reactions, while in this study we evaluate and quantify the self-purification ability of the hyporheic zone for groundwater-borne contaminant. Therefore, we believe that this study is quite different from Ping et al. (2022). In addition, we have obtained several interesting findings in this study, which can be summarized as: 1) bedform migration reshape the surface water-groundwater mixing zone and mixing-trigged biogeochemical zone patterns; 2) at slow to moderate migration rates, bedform migration expands the surface water-groundwater mixing zone and increases the fraction of mixing in hyporheic flux, with both parameters stabilizing at a plateau thereafter; 3) a large bedform migration celerity limits mixing-dependent denitrification rate, weakening the hyporheic zone’s role as a natural barrier against groundwater contaminant. In summary, we appreciate the concern from this reviewer on the novelty of this manuscript. We will revise the paper carefully to make it clear on the novelty of this paper, and clearly state the difference between this study and Ping et al. (2022). We are planning to revise literature review, the schematic diagram (figure 1), and some other parts, if necessary, to make it clearer. We are pretty sure that this concern can be addressed and this paper can be considered for potential publication in HESS after revision. 1.The literature review should be expanded to include a more thorough discussion of numerical model development in this field. Critically, this section must clearly and explicitly state how this study's model and objectives differ from and advance upon the 2022 paper. The work appears to be a case study applying a previously established model. It does not introduce a new mechanistic framework, a dimensionless parameter space map, or any analytical/scaling relationships that would constitute an independent and generalizable scientific contribution. Response to comment #1: We agree with this reviewer on this point. In the revised manuscript, we will make targeted revisions in Introduction section to thoroughly discuss the development of numerical models in this field, explicitly state the difference from the 2022 paper, and highlight the novelty of this research. 2.Thresholds used to define the mixing zone (e.g., 16%–84% concentration range) are presented without a clear physical or statistical basis. The robustness of the results is questionable without a sensitivity analysis of these thresholds, as well as the numerical dispersion and grid resolution. Furthermore, the newly proposed metrics are not validated against any established, well-accepted measures of mixing from the fluid dynamics or hydrology literature, making their utility and interpretation difficult to assess. Response to comment #2: We will supplement targeted analyses and validation to address these concerns, and we elaborate on the revisions and scientific rationale as follows: 1) Actually, there is no standard for the threshold of the mixing zone; however, it generally accepted that the mixing zone ranges from 10% to 90% in terms of groundwater proportion (Hester et al., 2013; Woessner et al., 2000). This interval effectively distinguishes the mixing zone, where groundwater and surface water interact dynamically, from the two endmembers: pure groundwater (>90%) and pure surface water (<10%). In this study, we selected the 16–84% range, which corresponds to the ±1σ interval of a normal distribution, as proposed by Santizo et al. (2020). To assess how threshold variations influence results, we will perform a sensitivity analysis that tests three alternative concentration ranges: 10%–90% (wider interval), 16%–84%, and 20%–80% (narrower range) in the revised paper. 2) We agree that the numerical dispersion and grid resolution can influence solute transport and mixing results. We will test three grid sizes to assess spatial discretization impacts: baseline grid (the original resolution), fine grid and coarse grid. Critical parameters would be compared for each grid size, confirming that the baseline grid is sufficiently resolved to capture mixing dynamics and minimize numerical dispersion. 3) The calculations of mixing metrics in this study follow the methods outlined by Hester et al. (2013, 2014). These metrics are not newly proposed by this study; instead, they have been widely applied and validated in previous studies (Santizo et al., 2020, 2022; Nogueira et al., 2022; Woessner et al., 2000), the evidence that supports their rationality. We will state that clearly in the revised manuscript. 3.The model's idealizations limit the generalizability of its findings. Key simplifications include: (a) representing bedform migration as a simple translating sinusoidal head, which neglects morphodynamic feedbacks; (b) using a 2-D domain with a single bedform wavelength and amplitude; and (c) omitting sediment heterogeneity, such as grain size variations or layering, which are known to strongly influence hyporheic exchange. Response to comment #3: (a) Actually, simplifying undulating riverbeds to a horizontal configuration and applying a sinusoidal head distribution at the sediment-water interface exerts minimal impacts on both the hyporheic exchange flux (HEF) driven by the riverbed and solute transport within the riverbed, which has been proven by many researchers (Elliott & Brooks, 1997a, b; Eylers, 1994; Rutherford et al., 1995). Indeed, we have conducted an additional simulation incorporating ripple geometry and compared it to the flat-bed model: both yield similar HEF and solute plume characteristics (this simulation will be added to the revised manuscript). (b) Given the periodicity of bedform geometry in our model, we defined the model domain using a single bedform wavelength and amplitude. This is a well-established practice in modeling, as it minimizes edge effects without the computational intensity of full-scale simulations of multiple consecutive bedforms. To validate this simplification, we simulated three consecutive ripples, focusing on the middle one: comparisons of its pressure, solute concentration, and vertical boundary gradients with the single ripple model showed only minor differences. We will state that the model applies most directly to straight, low-curvature streams with uniform bedform spacing (typical of agricultural/urban downstream reaches, consistent with our focus; Hester et al., 2014; Krause et al., 2022). Here, hyporheic exchange is dominated by streamwise-vertical flow cells, with lateral hyporheic flux contributing little to total exchange. (c) We omitted this factor in the current study to isolating the effects of bedform migration on surface water-groundwater mixing and mixing-triggered reaction. This work specifically targets downstream gaining reaches characterized by fine sand beds, which are a relatively homogeneous sediment type. To address this limitation, we will add relevant context to the Discussion section of the revised manuscript. We will propose that future studies integrate stochastic K fields to explore how sediment heterogeneity interacts with bedform migration, for instance, whether high-K hotspots amplify or dampen migration-driven mixing.
Citation: https://doi.org/10.5194/egusphere-2025-2631-AC3 -
AC4: 'Reply on RC3', Zhang Wen, 25 Sep 2025
General comments: This study investigates the effects of migrating bedforms on surface-groundwater mixing and denitrification in a gaining stream reach. The authors employ a numerical model that couples Darcy flow with a translating sinusoidal head boundary and a multi-species reactive transport model. While the paper introduces new metrics such as “mixing flux,” “mixing fraction,” and “mixing-zone area,” my primary concern is the substantial overlap with the authors' previous work (2022). The modeling framework, numerical implementation, and validation appear nearly identical. The claimed novelty rests on applying the existing model to a different scenario and introducing new diagnostic metrics. However, the justification for these metrics is insufficient, and their introduction does not lead to the discovery of new physical mechanisms. Furthermore, the main conclusion: bedform migration shortens residence times and thereby reduces denitrification—reiterates, is a key finding from the 2022 paper (https://doi.org/10.1029/2022WR033258). The manuscript reads more as an incremental application of prior work rather than a standalone scientific article with a distinct and significant contribution. Response to the general comments: We thank this reviewer for the constructive comments on our manuscript. We think the main concern of this reviewer is the difference between this study and our previous publication (Ping et al., 2022). We will classify the difference and the novelty of this study clearly in the revised version. Actually, this study is quite different from Ping et al. (2022) because: 1) Firstly, in this study, we focus on downstream section of rivers under gaining conditions firstly to simulate surface water-groundwater mixing processes under dynamic bedforms, while Ping et al. (2022) focused on hyporheic exchange processes in the headwater to midstream sections of rivers. therefore, the research focus is quite different for these two studies. 2) Secondly, the mixing process has been considered in this study. As far as we know, this process has not been thoroughly investigated in the migrating ripple model (Jiang et al., 2022; Zheng et al., 2019; Kessler et al., 2015). It is notable that the bedform migration affect mixing process significantly. 3) Thirdly, Ping et al. (2022) only considered stream-borne solutes and non-mixing dependent reactions, while in this study we evaluate and quantify the self-purification ability of the hyporheic zone for groundwater-borne contaminant. Therefore, we believe that this study is quite different from Ping et al. (2022). In addition, we have obtained several interesting findings in this study, which can be summarized as: 1) bedform migration reshape the surface water-groundwater mixing zone and mixing-trigged biogeochemical zone patterns; 2) at slow to moderate migration rates, bedform migration expands the surface water-groundwater mixing zone and increases the fraction of mixing in hyporheic flux, with both parameters stabilizing at a plateau thereafter; 3) a large bedform migration celerity limits mixing-dependent denitrification rate, weakening the hyporheic zone’s role as a natural barrier against groundwater contaminant. In summary, we appreciate the concern from this reviewer on the novelty of this manuscript. We will revise the paper carefully to make it clear on the novelty of this paper, and clearly state the difference between this study and Ping et al. (2022). We are planning to revise literature review, the schematic diagram (figure 1), and some other parts, if necessary, to make it clearer. We are pretty sure that this concern can be addressed and this paper can be considered for potential publication in HESS after revision. 1.The literature review should be expanded to include a more thorough discussion of numerical model development in this field. Critically, this section must clearly and explicitly state how this study's model and objectives differ from and advance upon the 2022 paper. The work appears to be a case study applying a previously established model. It does not introduce a new mechanistic framework, a dimensionless parameter space map, or any analytical/scaling relationships that would constitute an independent and generalizable scientific contribution. Response to comment #1: We agree with this reviewer on this point. In the revised manuscript, we will make targeted revisions in Introduction section to thoroughly discuss the development of numerical models in this field, explicitly state the difference from the 2022 paper, and highlight the novelty of this research. 2.Thresholds used to define the mixing zone (e.g., 16%–84% concentration range) are presented without a clear physical or statistical basis. The robustness of the results is questionable without a sensitivity analysis of these thresholds, as well as the numerical dispersion and grid resolution. Furthermore, the newly proposed metrics are not validated against any established, well-accepted measures of mixing from the fluid dynamics or hydrology literature, making their utility and interpretation difficult to assess. Response to comment #2: We will supplement targeted analyses and validation to address these concerns, and we elaborate on the revisions and scientific rationale as follows: 1) Actually, there is no standard for the threshold of the mixing zone; however, it generally accepted that the mixing zone ranges from 10% to 90% in terms of groundwater proportion (Hester et al., 2013; Woessner et al., 2000). This interval effectively distinguishes the mixing zone, where groundwater and surface water interact dynamically, from the two endmembers: pure groundwater (>90%) and pure surface water (<10%). In this study, we selected the 16–84% range, which corresponds to the ±1σ interval of a normal distribution, as proposed by Santizo et al. (2020). To assess how threshold variations influence results, we will perform a sensitivity analysis that tests three alternative concentration ranges: 10%–90% (wider interval), 16%–84%, and 20%–80% (narrower range) in the revised paper. 2) We agree that the numerical dispersion and grid resolution can influence solute transport and mixing results. We will test three grid sizes to assess spatial discretization impacts: baseline grid (the original resolution), fine grid and coarse grid. Critical parameters would be compared for each grid size, confirming that the baseline grid is sufficiently resolved to capture mixing dynamics and minimize numerical dispersion. 3) The calculations of mixing metrics in this study follow the methods outlined by Hester et al. (2013, 2014). These metrics are not newly proposed by this study; instead, they have been widely applied and validated in previous studies (Santizo et al., 2020, 2022; Nogueira et al., 2022; Woessner et al., 2000), the evidence that supports their rationality. We will state that clearly in the revised manuscript. 3.The model's idealizations limit the generalizability of its findings. Key simplifications include: (a) representing bedform migration as a simple translating sinusoidal head, which neglects morphodynamic feedbacks; (b) using a 2-D domain with a single bedform wavelength and amplitude; and (c) omitting sediment heterogeneity, such as grain size variations or layering, which are known to strongly influence hyporheic exchange. Response to comment #3: (a) Actually, simplifying undulating riverbeds to a horizontal configuration and applying a sinusoidal head distribution at the sediment-water interface exerts minimal impacts on both the hyporheic exchange flux (HEF) driven by the riverbed and solute transport within the riverbed, which has been proven by many researchers (Elliott & Brooks, 1997a, b; Eylers, 1994; Rutherford et al., 1995). Indeed, we have conducted an additional simulation incorporating ripple geometry and compared it to the flat-bed model: both yield similar HEF and solute plume characteristics (this simulation will be added to the revised manuscript). (b) Given the periodicity of bedform geometry in our model, we defined the model domain using a single bedform wavelength and amplitude. This is a well-established practice in modeling, as it minimizes edge effects without the computational intensity of full-scale simulations of multiple consecutive bedforms. To validate this simplification, we simulated three consecutive ripples, focusing on the middle one: comparisons of its pressure, solute concentration, and vertical boundary gradients with the single ripple model showed only minor differences. We will state that the model applies most directly to straight, low-curvature streams with uniform bedform spacing (typical of agricultural/urban downstream reaches, consistent with our focus; Hester et al., 2014; Krause et al., 2022). Here, hyporheic exchange is dominated by streamwise-vertical flow cells, with lateral hyporheic flux contributing little to total exchange. (c) We omitted this factor in the current study to isolating the effects of bedform migration on surface water-groundwater mixing and mixing-triggered reaction. This work specifically targets downstream gaining reaches characterized by fine sand beds, which are a relatively homogeneous sediment type. To address this limitation, we will add relevant context to the Discussion section of the revised manuscript. We will propose that future studies integrate stochastic K fields to explore how sediment heterogeneity interacts with bedform migration, for instance, whether high-K hotspots amplify or dampen migration-driven mixing.
-
AC5: 'Reply on AC4', Zhang Wen, 25 Sep 2025
The directly uploaded text does not show paragraph breaks; however, the complete response is also accessible via this attachment.
Citation: https://doi.org/10.5194/egusphere-2025-2631-AC5
-
AC5: 'Reply on AC4', Zhang Wen, 25 Sep 2025
-
AC3: 'Reply on RC3', Zhang Wen, 25 Sep 2025
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-2631', Anonymous Referee #1, 20 Aug 2025
The hyporheic zone plays a vital role in mitigating contaminants through interactions between surface water and groundwater. This study investigates the effects of bedform migration on the mixing dynamics of these waters and the accompanying biogeochemical reactions, particularly under gaining stream conditions, by utilizing a coupled hydrodynamic and reactive transport model. The research focus is meaningful, and the numerical model is well validated. However, I believe the overall clarity and presentation of the manuscript can be further improved. My suggestions are as follows:
Parameter explanations: I strongly recommend that the authors carefully review and expand the explanations of all parameters. Many parameters are currently not well defined, which may reduce accessibility for future readers.
Writing style: I suggest avoiding the use of "we" throughout the manuscript, as this expression is less common in formal academic writing. Rephrasing such sentences into more objective forms would enhance the manuscript’s style.
Specific comments:
Line 22: Please consider adding a comma between "migration" and "this" ("incorporates bedform migration this paper" → "incorporates bedform migration, this paper").
Notation consistency: I recommend keeping nitrate consistently as NO3^- throughout the manuscript (e.g., Line 81 uses NO3^-1, while Line 55 uses NO3).
Line 100: Please explain D50 here.
Line 114: Since the manuscript refers to "multi-component," I recommend using brackets to quickly introduce these components.
Line 125: Providing the dimensions of parameters (as done in Line 152) would be helpful.
Figure 1: The figure is somewhat confusing. Using a different color scheme for the mixing zone would improve clarity, as the current colors resemble a geological layer. Adding a note to indicate that the blue curve represents river elevation would also help. In the caption, I suggest rephrasing as: "NMD reaction = non-mixing-dependent reaction and MD reaction = mixing-dependent reaction."
Line 133: Please revise "'and SW out is …'" to "'SW out' is …."
Line 150: Since the bedform has a nonuniform distribution, would it be possible to consider anisotropic hydraulic conductivity in this study?
Line 153: As the abbreviation BC first appears here, it should be explained. (Since it also appears in Fig. 1 as line BC, I suggest not using BC as the abbreviation for boundary condition to avoid confusion.) Similarly, uc and dt in Eq. (2) require explanation. Given the large number of parameters, I recommend creating a notation table summarizing all variables with their definitions.
Line 161: Please clarify whether n and a have physical meaning or relationships. This would improve the organization of the manuscript.
Line 165: Explanation of Ri is missing.
Line 168: Explanations of vi and vj are missing. Also, note that v is a vector in Eq. (4) but a scalar in Eq. (5). Please clarify if they represent the same variable.
Line 171: Emphasizing why the three chemical species were selected would strengthen the manuscript (e.g., due to dominant reactions or research focus).
Line 178: Since Monod kinetics is referenced, I recommend directly citing Monod’s original work.
Line 190: In Line 150, you refer to hydraulic head, while here pressure is mentioned. Are these the same variable? Consistency is important. Additionally, please clarify the physical meaning of applying a periodic boundary condition.
Line 197: The phrase "a set of criteria" is unclear—are these criteria listed in the manuscript? Figure 2 is also difficult to interpret; more explanation would help.
Line 260: Just a positive note here—I appreciate how the manuscript justifies parameter values with references or physical reasoning. This strengthens the credibility of the work.
Figure 3: Including RMSE values (simulation vs. experiment) in each subplot would enhance the validation of the model.
Conclusion: The conclusion could be further strengthened by including quantitative results that highlight the main findings.
Citation: https://doi.org/10.5194/egusphere-2025-2631-RC1 -
AC1: 'Reply on RC1', Zhang Wen, 29 Aug 2025
Thank you for your valuable and insightful comments, which are greatly helpful in improving the clarity and overall quality of our manuscript. We will carefully address all the comments if we get a chance to revise the manuscript. Please find our briefly response to your comments in the following:
Regarding the parameter explanations: We will conduct a systematic review of all parameters in the manuscript, supplementing and expanding their definitions, including physical meanings and units, to improve readability. In addition, for the writing style: All instances of subjective phrasing (e.g., “we”) will be revised into objective expressions in accordance with standard academic writing conventions.
For the specific comments raised by this reviewer, we totally agree and we will accept all of your suggestions during the revision process of the paper.
Once again, thank you very much for your time spending on review our manuscript. If you have any other questions, please feel free to contact us.
Citation: https://doi.org/10.5194/egusphere-2025-2631-AC1
-
AC1: 'Reply on RC1', Zhang Wen, 29 Aug 2025
-
RC2: 'Comment on egusphere-2025-2631', Anonymous Referee #2, 31 Aug 2025
Ping X. et al reported a modeling study of the mixing process and mixing-dependent denitrification beneath a moving bed form. The topic is very important, and results are valuable. However, there are several critical issues in conceptual model description, ripple geometry simplification, nitrate concentration portioning, shear stress estimation, and model validation. Major issues were summarized below for the authors for consideration. Detailed comments are attached at the end.
Major comments:
- Simplifying ripples as a flat surface need justification. For non-mobile bedform, it is well known that non-uniform pressure distribution applied on the ripple surface cause pressure gradients horizontally and vertically, which further control the complex pressure, velocity, and concentration distributions inside and beneath the bedform. Due to its driving force, it is critical to include the ripple geometry in simulations, which is a common practice in most modeling research for stationary hyporheic zone (HZ) (see Cardenas2007, WRR). For mobile bedform, one of the few experiments on hyporheic zone with mobile bedforms show that the distribution of concentrations is mostly complex inside and near moving ripple (Wolke2020, Water). This also suggests the most important role of ripple geometry in reproducing realistic flow and transport distributions in HZ. Without representing ripple geometry in the model, the exchange flux, concentration distributions, and shear stress on, inside, and beneath the ripple will likely have high uncertainty if they are estimated from existing empirical equations. I suggest the authors justify such assumption by adding a new simulation with ripple geometry and then evaluate how much difference could be in exchange flux, concentration flux, and shear stress from models with and without ripple.
- The conceptual model presentation and Figure 1 needs improvements to better align with text descriptions, Figure 1 and texts in Section 2.2-2.6. In the current modeling approach for HZ, popular approaches are two-layer and three-layer conceptual model. Figure 1 suggests a two-layer model, i.e., splitting the domain as a surface water layer (regions above ripple top boundary) and groundwater layer (regions below ripple top boundary) with riverbed surface (ripple top boundary) as a permeable interface. However, when reading the details in Section 2.1-2.6, the model is a three-layer structure: surface water layer, riverbed layer (a thick 3D layer instead of a riverbed surface), and groundwater layer. Different from traditional 3-layer structure, this paper does not explicitly solve Navier-Stockes equations for surface water layer and Darcy-equation for groundwater layer. Insteady, this work only solves the pressure head and concentration for the riverbed layer with simplified boundary conditions for pressure head/concentrations at the top (boundary BC in Figure 1) and bottom (boundary AD in Figure 1) boundary of the riverbed layer. Further, this paper split the riverbed layer into two pools, i.e., non-mixing dependent reaction pool (blue region in Figure 1) and mixing dependent reaction pool (red region in Figure 1). Note that such “two” pools share the same computational domain (bounded by ABCD). They do not have an internal interface (critical for traditional two-layer or three-layer models) to conserve water and concentration mass. Insteady, this work defines two variables (c_{s-NO3-} and c_{g-NO3-}) and apply the advection-diffusion-reaction equation (equation 4) twice for these two variables. This means the model is essentially a one-domain multiple-variable model. However, most of this information was not clearly presented in Section 2.1, which is supposed to serve for this purpose. I suggest the authors modify Section 2.1 and relevant texts in Sections 2.2-2.7 to reflect the key conceptual model setups in Section 2.1. Please be aware that further readers may not be familiar with different conceptual models for modeling HZ. Your precise introduction is critical for them to correctly understand your work.
- It is not clear if mass conversation for pressure head and concentration can be achieved based on the current model descriptions. The boundary BC and DC is fine due to use simple “periodic” boundary. For boundary BC and AD, it is not clear to me how to guarantee consistent at boundary conditions at BC and AD for pressure head, head gradient (used to further compute exchange flux), and concentration gradient due to unclear boundary condition descriptions.
- Partitioning the total Nitrate into surface water layer contribution (c_{s-NO3-}) and groundwater layer contribution (c_{g-NO3-}) maybe problematic due to the involvement of non-linear reaction rate in the Monod kinetics. From equation 7, the reaction rate of denitrification R_{DN} is a function of Nitrate concentration, i.e., R_{DN} = f(NO3). If simulating the total NO3- with portioning, there will be two reaction rate R_{DN}(c_{s-NO3-) (equation 8) and R_{DN}(c_{g-NO3-) (equation 9). Due to non-linearity of Monod equation f, R_{DN} does not equal the summation of the reaction rate from two partition components, i.e., f(NO3) ≠ f(sNO3) + f(gNO3). This means partition approach is not mass-conservation. Additionally, introducing two concentration variables means need 4 more boundary conditions at top and boundary conditions for c_{s-NO3-} and c_{g-NO3-}. This paper does not clearly explain how such boundary conditions are assigned. In short, the current Nitrate concentration approach causes mass conservation issues and needs 4 more boundary conditions data which are not well described. I suggest the authors clarify these issues in the revision.
- Estimation of shear stress that drives sand transport is likely not reliable. The estimation of migration velocity uc depends on shear velocity and shear stress on the ripple surface. Due to the simplification of ripple into a flat surface, calculation of shear stress/shear velocity is not possible. This paper estimates the shear velocity using the Cheezy equation. The Cheezy equation usually works well for flat surface with homogenous sediments. In this work, the sand is very small (0.08-0.36 mm) and homogenous. However, the ripple is much bigger (height 0.02 m). The correct shear stress is the one on the surface with both large ripple and small sand. The use of Cheezy is likely not adequate to estimate such variables. The most accurate way will be using numerical model to directly compute shear stress with the ripple geometry as the top boundary. If the authors want to use the Cheezy equation, it needs to be verified against at least one CFD simulated or observed shear stress to quantify potential uncertainty.
- What were compared in the model validation is not clear. Figure 5 shows a comparison between model and observation. However, how HEF, AO, Fin, Fout were computed is not explained clearly. The model results are distributed data. The authors need to explain how to convert distributed values into single values like HEF, AO, Fin, Fout. Second, the oxygen concentration from the model is affected by reaction R1. It is not clear if the experiment also implemented a similar reaction. If the model has a reaction R1 and the experiment does not, then it is not reasonable to compare model predicted RO with experiment observed RO. Please add more information about the details of model validation. Especially, please clarify if the model setup is consistent with the experiment setup, especially the reaction part.
- The validation result from Figure 5 suggests model prediction for HEF and concentration flux may be reliable (if potential issues in Comment 5 get clarified). However, simulation for concentration distribution (which could be used to compute AO in Figure 5b) shows much larger uncertainty. Because this work focuses on analyzing denitrification that depends on accurate prediction of nitrate concentration. I am concerned that the prediction for Nitrate may have similar uncertainty level as the RO shown in Figure 5b. I suggest the authors provide additional information to verify if the prediction uncertainty for nitrate is negligible for your analysis. If lab experiment is not available. A potential solution is to run your model with the ripple geometry explicitly represented by mesh and compute the difference of Nitration distribution with and without ripple geometry. According to Line 349, the authors seem to acknowledge that the larger error in RO shown in Figure 5b could be caused by the exclusion of the ripple geometry. Such addition modeling work will help to explain Figure 5b and enhance the prediction reliability for nitrate distribution.
- This work has several figures (e.g., figures 1-4) and most equations similar to the Ping X. 2022 (10.1029/2022WR033258), please highlight major novelty in this paper compared to the existing paper.
Other minor comments are listed below for reference:
Line 22: missing a comma between “migration” and “this”.
Line 123: Riverbed and riverbed surface mean different things. If the “Riverbed” at this sentence means “the top surface of the riverbed”, it is better to use “riverbed surface”.
Figure 1: Figure 1 looks very similar to Figure 2 in Ping2022 (10.1029/2022WR033258). The ripple shapes are the same in two papers. I guess this figure is modified from Figure 2 in Ping 2022. In this scenario, please add the following information in figure caption: [Figure is modified from Figure 2 in Ping 2022, permission has granted from Wiely]
Figure 1: Because HZ is a complex system, it involves concepts such as domain, layer, surface, interface, flux on surface, in from domain/surface, out from domain/surface. I suggest the authors modifying the Figure 1 and captions to carefully choose vocabularies to “concisely” describe such concepts. For example, riverbed could mean riverbed surface or riverbed layer. Clearly selecting words will help readability of the paper. This comment is also related to Major comment 2.
Line 144: I am not quite agree with this. From Elliot1997 paper, the size of ripple heights is 1-2.5 cm, and surface water depth is of 3-7 cm, which means ripple heights are 15-20% of the water depth. In Cardenas2007 CFD simulation paper (10.1016/j.advwatres.2006.06.009), dune height is 5 cm and water depth is 45 cm. Their ratio is 10% (Figure 2 Cardenas2007). At line 144, I don't see the actual size for ripple and depth, so it is reasonable to say ripple heights are negligible. Such negligible assumption will lead to the simplification of "triangular ripple” into a "flat" surface. From both Elliot1997 and Cardenas 2007 paper, they do not simplify the ripple/dune into flat surface. Instead, Elliot1997 shows the pressure varies sinuously along the surface of the ripple. In Cardenas 2007, it shows that pressure head variations are affected by Reynolds number. The sinusoidal function reported in Elliot1997 is valid only at low Reynolds number and such match is applied only half of the region along the ripple surface (see Figure 9 in Cardenas 2007 10.1016/j.advwatres.2006.06.009). Therefore, whether you can simplify the ripple as flat surface (which will affect how you design the computational domain) and using a sinusoidal function (depending on your Reynolds number) requires more justification. This comment is also related to Major comment 1.
Line 142: show lambda on Figure 1. Also please visualize ripple height in Figure 1.
Line 146: A periodic morphological feature does not mean periodic conditions for pressure, velocity, and concentration. Previous work shows that the flow in upstream of ripple influences flow downstream (see Sinha2017 WRR, 10.1002/2016WR019662). Due to this issue, numerical simulation will need to use at least two ripples to minimize impacts of upstream (see Cardenas2007, 10.1016/j.advwatres.2006.06.009). In this work, a period boundary is applied at the location B and C (see Figure 1) where pressure and velocity are not periodic (see Figure 2 in Cardenas2007). However, all the information mentioned above is for stationary ripple/dune. For the moving ripple, whether you should/should not use a periodic boundary may be different but requires further justification.
Line 153: Is this "BC" the top edge of the blue region in Figure 1? In this work, the 'BC' is a flat surface. However, equation 2 is usually applied to the top surface of the ripple which is not a flat surface. Please justify if such simplification affects results (e.g., exchange flux and distribution of nitrate in the computational domain).
Line 154: equation 2 was first introduced in Boano2013. It is better to cite correct references.
Line 163: For equation 4/5, you use bold (v), non-bold v, vi, vj, and |v|. In general, bold (v) means vector. For non-bold v, there is no definition. Also, no definition for vi, vj. Please correct such variables to make sure they are consistent and mathematically correct.
Line 164 equation 4: gradient v.ci should be gradient dot (v*ci).
Line 173: Does this sentence mean you solved both c_{sNO3} and c_{gNO3}? Note that most previous work does not distinguish which portion of nitrate from surface water and which portion is from groundwater, which means they only solve for one concentration for NO3. From this sentence, it seems you are solving both concentrations of Nitrate. Please confirm if this is the case. Also, for models that solve for one concentration for NO3, you only need to assign two boundary condition for NO3, one is at the top boundary (BC) and one is for the lower boundary (AD). But if you are solving for two concentrations for NO3, you need to assign 4 boundary conditions. You need boundary conditions at both BC and AD for both c_{sNO3} and c_{gNO3}. Please clarify how each of these boundary conditions is determined.
Line 176 equation R1/R2： it is better to add a citation for equations R1 and R2. It is also important to add a justification for why you choose such a reaction. Note that in lab or field experiments, the reaction network is quite complex, which is much different than the equation R1/R2 used here. R1/R2 is a kind of idealized or conceptualized reaction, which could be used as qualitative understating. However, to compare the results to experiment or field observation, you may need to add more information regarding how the assumed reaction R1/R2 is consistent with lab/field reaction network.
Line 180 equation 7: I am not quite sure the meaning of C_{s/g-NO3-}. Does this mean the summation of C_{s-NO3-} and C_{s-NO3-}? Or it means C_{s-NO3-} or C_{s-NO3-}? I also have doubts about the mass balance issue caused by this equation. Please see my major comment 3 for details.
Line 190: the period boundary ab AB and DC needs more justification. From previous work (Cardenas2007 Figure 2 10.1016/j.advwatres.2006.06.009), the distribution of pressure at AB and DC could be quite different depending on Reynolds number. A better way is to use two ripples and apply period at the beginning of the first ripple and the end of the second ripple. Please justify the current setup.
Line 190: If P(0) = P(lamba) + delta P, this is not called period boundary. You can just use the mathematical formula to describe it, but do not call it "periodic". Also, equation 1 means you are solving for pressure head h, but no pressure P. You need to assign the boundary for h, but not P.
Line 193: How do you determine the concentrations at boundary BC? For pressure head h, it is determined by the sinusoidal function from previous work, but for concentration, how do you get it.
Line 195: What is the boundary condition for h at AD boundary?
Line 213: using the Chezy equation to estimate shear stress seems not reliable. Please see my major comment 4 for details.
Line 215: How do you calculate the tau_{cr}?
Line 217: what is r in the equation used for computing shields number at line 217?
Line 231: Is the S at line 231 the same S at line 212? Do you define K somewhere? Do you provide the value for K?
Line 241: Is uq upwelling groundwater flux? How do you decide this value?
Line 264: Is the Da at line 264 the same as the Da in equation 18? If not, please use a different variable.
Line 266: How are these values estimated? Are these assumed values? Please provide a justification for why you choose such values.
Table 1: Please visualize l, lc, lambda, Hd, H in Figure 1. What is the difference between ripple crest and height of ripple crest. In Figure 1, the ripple height is clearly visualized, but the ripple crest is not visualized. Please describe their difference.
Tabel 1: Provide citation for porosity, longitudinal and transverse dispersity or any justification to use such values.
Line 283: Does COMSOL solve the equations 1-11? Do you develop your own COMSOL solver using equations 1-11 or you use the default equations in COMSOL? Is the COMSOL setup and data available in public data repository?
Line 291: I am confused by the term "surface water". From Figure 1, your computational domain is the region bounded by ABCD. The whole region is beneath the surface (or stream) water. This means the model framework does not solve anything related to surface water that is not governed by Darcy flow. Please clarify this.
Line 301: Is bottom boundary AD and streambed surface BC?
Line 302: Is “streambed surface” the "surface water layer" (a layer with pure water)?
Line 303: more accurate description: surface water enters riverbed layer.
Line 304-305: From Figure 1, the groundwater is separated from the surface water layer by "the sediment layer", it seems groundwater cannot be directly transported to the surface water layer. When you use "groundwater discharged to stream", are you trying to identify which portion of water is from groundwater when such groundwater is first "up-welled" to the sediment layer and further "up-welled" to the surface water layer? If this is the case, you may need to explicitly describe this purpose. And you may add a sentence that "such identification could be achieved using tracer or different color/size of tracer".
Line 317: Where is the ripple surface? Is ripple surface the boundary BC? How do you compute Fswi?
Line 322: Is this streambed "boundary BC"?
Figure 3b: From Wolke2019, the oxic zone occurs mainly inside the dune (the triangular region as shown in Figure 1). However, in this paper, the computational domain is a rectangular shape, and the triangular shape is "neglected". How do you compare the oxic zone as reported in the Wolke2019 paper? From Figure 3b, the uncertainty for AO between model and experiment is much larger than the rest 3 variables, could this larger error be caused by the over-simplification of the ripple to a flat surface? Also see major comment 6.
Figure 3: Is this the average upwelling flux at the boundary BC? If not, please clarify where and how the HEF is computed. From the numerical, the flux is a distribution, it is necessary to clarify how you get a single value from a distribution. See also major comment 5.
Line 345: From Figure 3b, the relative error between model and experiments could be unto 100%. The statement of "slightly lower" is not consistent with the data. Please compute the relative error between model results and observation and then make a more accurate conclusion for the accuracy using the computed data.
Line 349: This is related to my major comment 1 and 6. As suggested by the author, the missing of the mobile (triangular ripple) can cause significant error for estimating the total area of the oxygen zone. As the rest of the paper is to study the denitrification which will use data of concentrations of nitrate. My guess is that the missing of the triangular ripple zone (mobile zone) will also significantly affect the accuracy of nitration concentration prediction (could be similar uncertainty of oxygen as shown in Figure 3b). I suggest the authors clarify how much impact will be by simplifying the mobile ripple as flat surface.
Line 375 and Figure 4: This is something uncleaned at lower left corner of Figure 4.
Line 474: This paper has a lot of similarity to Ping 2022 (10.1029/2022WR033258). It may be over statement to use “first time” due to the existing of Xue2022 paper.
Line 852: citation year is not correct.
Citation: https://doi.org/10.5194/egusphere-2025-2631-RC2 -
AC2: 'Reply on RC2', Zhang Wen, 17 Sep 2025
Thank you for your valuable and insightful comments, which are really helpful in improving the clarity and overall quality of our manuscript. We agree with all your concerns, and we will carefully address all the comments if we get a chance to revise the manuscript. Please find our briefly response to your major comments in the following:
1. Simplifying ripples as a flat surface need justification. For non-mobile bedform, it is well known that non-uniform pressure distribution applied on the ripple surface cause pressure gradients horizontally and vertically, which further control the complex pressure, velocity, and concentration distributions inside and beneath the bedform. Due to its driving force, it is critical to include the ripple geometry in simulations, which is a common practice in most modeling research for stationary hyporheic zone (HZ) (see Cardenas2007, WRR). For mobile bedform, one of the few experiments on hyporheic zone with mobile bedforms show that the distribution of concentrations is mostly complex inside and near moving ripple (Wolke2020, Water). This also suggests the most important role of ripple geometry in reproducing realistic flow and transport distributions in HZ. Without representing ripple geometry in the model, the exchange flux, concentration distributions, and shear stress on, inside, and beneath the ripple will likely have high uncertainty if they are estimated from existing empirical equations. I suggest the authors justify such assumption by adding a new simulation with ripple geometry and then evaluate how much difference could be in exchange flux, concentration flux, and shear stress from models with and without ripple.
Response to comment #1: We acknowledge the importance of ripple geometry in HZ modeling and address your concern by supplementing with additional simulations and literature support. Actually, simplifying undulating riverbeds to a horizontal configuration and applying a sinusoidal head distribution at the water-sediment interface exerts minimal impacts on both the hyporheic exchange flux (HEF) driven by the riverbed and solute transport within the riverbed, which has been proven by many researchers (Elliott & Brooks, 1997a, b; Eylers, 1994; Rutherford et al., 1995). Indeed, we have conducted an additional simulation incorporating ripple geometry and compared its results with those of the flat-bed model as the reviewer suggested; the two sets of results exhibit similar characteristics in hyporheic exchange flux and solute plume distributions (This simulation will be added in the revised version). Furthermore, to reproduce the experiment by Wolke (2020), we developed two models: a flat-bed model and a triangular bedform model. Based on the parameter adjustments in the previous manuscript text, we refined the numerical model to achieve a better fit to the experimental data. Its adequacy is further supported by our successful reproduction of Wolke’s (2020) results, thereby confirming its credibility for addressing the study’s objectives. Therefore, in the revised manuscript, we are planning to add two critical contexts: first, a comparison between the triangular bedform model and the flat-bed model; second, model validation against Wolke (2020) using both of these models.
2. The conceptual model presentation and Figure 1 needs improvements to better align with text descriptions, Figure 1 and texts in Section 2.2-2.6. In the current modeling approach for HZ, popular approaches are two-layer and three-layer conceptual model. Figure 1 suggests a two-layer model, i.e., splitting the domain as a surface water layer (regions above ripple top boundary) and groundwater layer (regions below ripple top boundary) with riverbed surface (ripple top boundary) as a permeable interface. However, when reading the details in Section 2.1-2.6, the model is a three-layer structure: surface water layer, riverbed layer (a thick 3D layer instead of a riverbed surface), and groundwater layer. Different from traditional 3-layer structure, this paper does not explicitly solve Navier-Stockes equations for surface water layer and Darcy-equation for groundwater layer. Insteady, this work only solves the pressure head and concentration for the riverbed layer with simplified boundary conditions for pressure head/concentrations at the top (boundary BC in Figure 1) and bottom (boundary AD in Figure 1) boundary of the riverbed layer. Further, this paper split the riverbed layer into two pools, i.e., non-mixing dependent reaction pool (blue region in Figure 1) and mixing dependent reaction pool (red region in Figure 1). Note that such “two” pools share the same computational domain (bounded by ABCD). They do not have an internal interface (critical for traditional two-layer or three-layer models) to conserve water and concentration mass. Insteady, this work defines two variables (c_{s-NO3-} and c_{g-NO3-}) and apply the advection-diffusion-reaction equation (equation 4) twice for these two variables. This means the model is essentially a one-domain multiple-variable model. However, most of this information was not clearly presented in Section 2.1, which is supposed to serve for this purpose. I suggest the authors modify Section 2.1 and relevant texts in Sections 2.2-2.7 to reflect the key conceptual model setups in Section 2.1. Please be aware that further readers may not be familiar with different conceptual models for modeling HZ. Your precise introduction is critical for them to correctly understand your work.
Response to comment #2: We apologize for the confusion and acknowledge that the conceptual model was not clearly explained in the submitted version. It is notable that this model is a one-domain, multiple-variable model which focus on porewater dynamics and multicomponent solute transport within riverbed sediments; its domain corresponds to the actual rectangular riverbed geometry (bounded by ABCD) with no further subdivision of the riverbed interior into two pools. Instead of simulating turbulent flow over bedforms, we used an idealized scenario where sinusoidal pressure variations are applied to a flat-bed (boundary BC in Figure 1). We regret any misunderstanding caused by Figure 1 and will revise the domain schematic to accurately reflect the conceptual model’s configurations. We will overhaul Section 2 to serve as the “conceptual model foundation”: it will explicitly state the model type, define the modeled domain, clarify key processes and assumptions (e.g., no surface water flow simulation, idealized top boundary pressure), and detail the two nitrate variables (origins and boundary conditions). Additionally, we will systematically cross-check and revise Sections 2.2-2.7 to ensure fully consistency with Section 2.1.
3. It is not clear if mass conversation for pressure head and concentration can be achieved based on the current model descriptions. The boundary BC and DC is fine due to use simple “periodic” boundary. For boundary BC and AD, it is not clear to me how to guarantee consistent at boundary conditions at BC and AD for pressure head, head gradient (used to further compute exchange flux), and concentration gradient due to unclear boundary condition descriptions.
Response to comment #3: First, we would like to explain the boundaries in detail. The top boundary BC is prescribed with a Dirichlet head condition (exchange flux is not prescribed, but computed via qex (x, t) = -n·K▽h|BC). The bottom boundary AD is set as Neumann boundary with specified upward flux: -n·K▽h|AD = q (x) (where q > 0 denotes upwelling groundwater into the modeled domain). These two boundaries are mutually compatible, one is prescribing head and the other one is normal flux, ensuring a well-posed mixed boundary setup. For reactive solute transport, the top boundary serves as an open boundary for stream-borne solutes, while the bottom boundary is set as an outflow condition. In contrast, these two boundary conditions are swapped for groundwater-borne nitrate. Open boundary where the flow direction switches steadily between inflow and outflow in accordance with the computed normal Darcy flux qn = n·q: c = cext (x, t) (Dirichlet) if qn < 0 (inflow) and n·(D ▽ c) = 0 (zero diffusive flux) if qn ≥ 0 (outflow). This standard inflow-outflow treatment prevents artificial sources/sinks at outflow and is conservative by construction. In addition, we verified fluid mass balance by calculating the total influx and outflow fluxes across all boundaries over simulated period for Re = 3000, Ub = 0.6, and D50 = 0.15 mm as an example. We have quantified the fluid influx and outflux, as well as the solute input and output, along with the internal consumption. Notably, the relative error remains less than the acceptable threshold (5%) for hydrological models, which effectively confirms the conservation of fluid and solute mass. We will clarify the boundary naming, state the governing equations, and give the exact boundary conditions used for both porewater flow and reactive solute transport, and add the results of mass-balance checks in the revised manuscript.
4. Partitioning the total Nitrate into surface water layer contribution (c_{s-NO3-}) and groundwater layer contribution (c_{g-NO3-}) maybe problematic due to the involvement of non-linear reaction rate in the Monod kinetics. From equation 7, the reaction rate of denitrification R_{DN} is a function of Nitrate concentration, i.e., R_{DN} = f(NO3). If simulating the total NO3- with portioning, there will be two reaction rate R_{DN}(c_{s-NO3-) (equation 8) and R_{DN}(c_{g-NO3-) (equation 9). Due to non-linearity of Monod equation f, R_{DN} does not equal the summation of the reaction rate from two partition components, i.e., f(NO3) ≠ f(sNO3) + f(gNO3). This means partition approach is not mass-conservation. Additionally, introducing two concentration variables means need 4 more boundary conditions at top and boundary conditions for c_{s-NO3-} and c_{g-NO3-}. This paper does not clearly explain how such boundary conditions are assigned. In short, the current Nitrate concentration approach causes mass conservation issues and needs 4 more boundary conditions data which are not well described. I suggest the authors clarify these issues in the revision.
Response to comment #4: We defined stream-borne nitrate and groundwater-borne nitrate as two distinct reactants. Rather than partitioning the nitrate pool into two separate fractions, we solve for the transport of these two nitrate types in the porous medium separately using the advection-dispersion-reaction equation. As outlined in our response to Comment #3, each nitrate type (differentiated by origin) has its own set of boundary conditions. We quantified the input and output of stream-borne and groundwater-borne nitrate, as well as their internal consumption via denitrification. Notably, the relative error remains below the acceptable threshold (less than 5%), confirming the conservation of solute mass. We will provide the exact boundary conditions used for both porewater flow and reactive solute transport, and add the results of mass-balance verification in the revised manuscript.
5. Estimation of shear stress that drives sand transport is likely not reliable. The estimation of migration velocity uc depends on shear velocity and shear stress on the ripple surface. Due to the simplification of ripple into a flat surface, calculation of shear stress/shear velocity is not possible. This paper estimates the shear velocity using the Cheezy equation. The Cheezy equation usually works well for flat surface with homogenous sediments. In this work, the sand is very small (0.08-0.36 mm) and homogenous. However, the ripple is much bigger (height 0.02 m). The correct shear stress is the one on the surface with both large ripple and small sand. The use of Cheezy is likely not adequate to estimate such variables. The most accurate way will be using numerical model to directly compute shear stress with the ripple geometry as the top boundary. If the authors want to use the Cheezy equation, it needs to be verified against at least one CFD simulated or observed shear stress to quantify potential uncertainty.
Response to comment #5: We acknowledge the limitations of the classical Chezy equation in capturing ripple-induced shear stress. In our study, it was only used to estimate baseline shear velocity for incipient sediment motion. To account for bedform development, we incorporated Coleman and Melville’s (1994) empirical relationship for ripple migration celerity, indirectly considering ripple-associated form drag. This two-step approach (grain-scale resistance for particle initiation and empirical ripple-scale form drag) is widely adopted and validated, and successfully reproduced Wolke’s (2020) results, supporting its credibility for our objectives. Besides, our primary focus was on the effects of bedform migration on solute transport within the riverbed. Rather than simulating turbulent flow over triangular bedforms and calculating shear stress via CFD, we employed two approaches: an idealized model in which sinusoidal pressure variations are applied to the sediment layer surface, and empirical equations to determine bedform migration celerity. We fully agree that CFD simulations would yield more precise estimates of pressure profile and shear stress, but this lies beyond the scope of our current study. In the revised manuscript. We will explicitly state the limitations of our study, propose CFD validation as a valuable direction for future work, and confirm that the current empirical approach is appropriate for the study’s context and objectives.
6. What were compared in the model validation is not clear. Figure 5 shows a comparison between model and observation. However, how HEF, AO, Fin, Fout were computed is not explained clearly. The model results are distributed data. The authors need to explain how to convert distributed values into single values like HEF, AO, Fin, Fout. Second, the oxygen concentration from the model is affected by reaction R1. It is not clear if the experiment also implemented a similar reaction. If the model has a reaction R1 and the experiment does not, then it is not reasonable to compare model predicted RO with experiment observed RO. Please add more information about the details of model validation. Especially, please clarify if the model setup is consistent with the experiment setup, especially the reaction part.
Response to comment #6: We regret the lack of clarity in the model validation description and will clarify this in both this response and the revised manuscript. Specifically, we will detail the experimental measurement methods and simulation-based calculation approaches for the key parameters: hyporheic exchange flux (HEF), oxygen influx/outflux (Fin/Fout), and oxygenated area (AO). The experimentally observed RO is calculated from oxygen influx/outflux and mean oxygen penetration depth, and represents the volumetric oxygen consumption rate within the whole domain. The good fit between simulated and measured Fin/Fout, together with the spatial distribution of oxygen plumes, will ensure a minor discrepancy in this RO-derived rate. Additionally, RO differs in meaning from R1 within the model: the latter specifically refers to the reaction rate per unit pore volume. In the modeling, oxygen transport was governed by the advection-dispersion-reaction equation, where oxygen consumption occurred through aerobic respiration following the Monod kinetic R1. Therefore, we did not compare the model-predicted R1 with the experimentally observed RO. We will provide a clear introduction to and description of the model setup and the corresponding parameter calculation methods in the revised manuscript.
7. The validation result from Figure 5 suggests model prediction for HEF and concentration flux may be reliable (if potential issues in Comment 5 get clarified). However, simulation for concentration distribution (which could be used to compute AO in Figure 5b) shows much larger uncertainty. Because this work focuses on analyzing denitrification that depends on accurate prediction of nitrate concentration. I am concerned that the prediction for Nitrate may have similar uncertainty level as the RO shown in Figure 5b. I suggest the authors provide additional information to verify if the prediction uncertainty for nitrate is negligible for your analysis. If lab experiment is not available. A potential solution is to run your model with the ripple geometry explicitly represented by mesh and compute the difference of Nitration distribution with and without ripple geometry. According to Line 349, the authors seem to acknowledge that the larger error in RO shown in Figure 5b could be caused by the exclusion of the ripple geometry. Such addition modeling work will help to explain Figure 5b and enhance the prediction reliability for nitrate distribution.
Response to comment #7: Two issues that previously reduced the model’s credibility have been addressed: first, the unclear description of model validation (including experimental measurement methods and simulation-based calculation approaches); second, the poor fit between simulated results and experimental data. Specifically, we have refined the numerical parameters to improve the fit of AO, Fin, and Fout, and further optimized the numerical model based on the parameter adjustments outlined in the original manuscript, achieving better alignment with the experimental data. As suggested by the reviewer, we additionally conducted a simulation incorporating ripple geometry (intended to reproduce both Wolke’s (2020) experiment and our model settings) and compared its results with those of the flat-bed model. Therefore, in the revised manuscript, we are planning to add two critical sections of context: first, a detailed comparison between the triangular bedform model and the flat-bed model; second, a comprehensive model validation against Wolke’s (2020) data using both of these models.
8. This work has several figures (e.g., figures 1-4) and most equations similar to the Ping X. 2022 (10.1029/2022WR033258), please highlight major novelty in this paper compared to the existing paper.
Response to comment #8: This is an extension of Ping et al. (2022), but it is quite different from the previous publication. The main novelties of this paper can be summarized as: 1) we focus on the downstream “gaining conditions” and first simulate surface water-groundwater mixing under dynamic bedforms; 2) beyond analyzing the variations in hyporheic exchange flux and cell size, we also assess how the surface water-groundwater mixing flux and mixing zone size change with bedform migration; 3) unlike the previous study, which only considers stream-borne nitrate (a limited scope), we quantify both the reduction rate and efficiency of stream-borne and groundwater-borne nitrate. While basic figures/equations overlap (for fundamental hyporheic processes), these innovations fill mixing-reaction coupling gaps and improve understanding of hyporheic nitrate biogeochemistry. We will add this context to the revised manuscript and clearly articulate the novel contributions of this study in the revised Introduction.
For the specific comments raised by this reviewer, we totally agree and we will accept all of your suggestions during the revision process of the paper.
Once again, thank you very much for your time spending on review our manuscript. If you have any other questions, please feel free to contact us.
Citation: https://doi.org/10.5194/egusphere-2025-2631-AC2
-
RC3: 'Comment on egusphere-2025-2631', Anonymous Referee #3, 19 Sep 2025
This study investigates the effects of migrating bedforms on surface-groundwater mixing and denitrification in a gaining stream reach. The authors employ a numerical model that couples Darcy flow with a translating sinusoidal head boundary and a multi-species reactive transport model. While the paper introduces new metrics such as “mixing flux,” “mixing fraction,” and “mixing-zone area,” my primary concern is the substantial overlap with the authors' previous work (2022). The modeling framework, numerical implementation, and validation appear nearly identical.
The claimed novelty rests on applying the existing model to a different scenario and introducing new diagnostic metrics. However, the justification for these metrics is insufficient, and their introduction does not lead to the discovery of new physical mechanisms. Furthermore, the main conclusion: bedform migration shortens residence times and thereby reduces denitrification—reiterates, is a key finding from the 2022 paper (https://doi.org/10.1029/2022WR033258). The manuscript reads more as an incremental application of prior work rather than a standalone scientific article with a distinct and significant contribution.
1.The literature review should be expanded to include a more thorough discussion of numerical model development in this field. Critically, this section must clearly and explicitly state how this study's model and objectives differ from and advance upon the 2022 paper. The work appears to be a case study applying a previously established model. It does not introduce a new mechanistic framework, a dimensionless parameter space map, or any analytical/scaling relationships that would constitute an independent and generalizable scientific contribution.2.Thresholds used to define the mixing zone (e.g., 16%–84% concentration range) are presented without a clear physical or statistical basis. The robustness of the results is questionable without a sensitivity analysis of these thresholds, as well as the numerical dispersion and grid resolution. Furthermore, the newly proposed metrics are not validated against any established, well-accepted measures of mixing from the fluid dynamics or hydrology literature, making their utility and interpretation difficult to assess.
3.The model's idealizations limit the generalizability of its findings. Key simplifications include: (a) representing bedform migration as a simple translating sinusoidal head, which neglects morphodynamic feedbacks; (b) using a 2-D domain with a single bedform wavelength and amplitude; and (c) omitting sediment heterogeneity, such as grain size variations or layering, which are known to strongly influence hyporheic exchange.
Citation: https://doi.org/10.5194/egusphere-2025-2631-RC3 -
AC3: 'Reply on RC3', Zhang Wen, 25 Sep 2025
General comments: This study investigates the effects of migrating bedforms on surface-groundwater mixing and denitrification in a gaining stream reach. The authors employ a numerical model that couples Darcy flow with a translating sinusoidal head boundary and a multi-species reactive transport model. While the paper introduces new metrics such as “mixing flux,” “mixing fraction,” and “mixing-zone area,” my primary concern is the substantial overlap with the authors' previous work (2022). The modeling framework, numerical implementation, and validation appear nearly identical. The claimed novelty rests on applying the existing model to a different scenario and introducing new diagnostic metrics. However, the justification for these metrics is insufficient, and their introduction does not lead to the discovery of new physical mechanisms. Furthermore, the main conclusion: bedform migration shortens residence times and thereby reduces denitrification—reiterates, is a key finding from the 2022 paper (https://doi.org/10.1029/2022WR033258). The manuscript reads more as an incremental application of prior work rather than a standalone scientific article with a distinct and significant contribution. Response to the general comments: We thank this reviewer for the constructive comments on our manuscript. We think the main concern of this reviewer is the difference between this study and our previous publication (Ping et al., 2022). We will classify the difference and the novelty of this study clearly in the revised version. Actually, this study is quite different from Ping et al. (2022) because: 1) Firstly, in this study, we focus on downstream section of rivers under gaining conditions firstly to simulate surface water-groundwater mixing processes under dynamic bedforms, while Ping et al. (2022) focused on hyporheic exchange processes in the headwater to midstream sections of rivers. therefore, the research focus is quite different for these two studies. 2) Secondly, the mixing process has been considered in this study. As far as we know, this process has not been thoroughly investigated in the migrating ripple model (Jiang et al., 2022; Zheng et al., 2019; Kessler et al., 2015). It is notable that the bedform migration affect mixing process significantly. 3) Thirdly, Ping et al. (2022) only considered stream-borne solutes and non-mixing dependent reactions, while in this study we evaluate and quantify the self-purification ability of the hyporheic zone for groundwater-borne contaminant. Therefore, we believe that this study is quite different from Ping et al. (2022). In addition, we have obtained several interesting findings in this study, which can be summarized as: 1) bedform migration reshape the surface water-groundwater mixing zone and mixing-trigged biogeochemical zone patterns; 2) at slow to moderate migration rates, bedform migration expands the surface water-groundwater mixing zone and increases the fraction of mixing in hyporheic flux, with both parameters stabilizing at a plateau thereafter; 3) a large bedform migration celerity limits mixing-dependent denitrification rate, weakening the hyporheic zone’s role as a natural barrier against groundwater contaminant. In summary, we appreciate the concern from this reviewer on the novelty of this manuscript. We will revise the paper carefully to make it clear on the novelty of this paper, and clearly state the difference between this study and Ping et al. (2022). We are planning to revise literature review, the schematic diagram (figure 1), and some other parts, if necessary, to make it clearer. We are pretty sure that this concern can be addressed and this paper can be considered for potential publication in HESS after revision. 1.The literature review should be expanded to include a more thorough discussion of numerical model development in this field. Critically, this section must clearly and explicitly state how this study's model and objectives differ from and advance upon the 2022 paper. The work appears to be a case study applying a previously established model. It does not introduce a new mechanistic framework, a dimensionless parameter space map, or any analytical/scaling relationships that would constitute an independent and generalizable scientific contribution. Response to comment #1: We agree with this reviewer on this point. In the revised manuscript, we will make targeted revisions in Introduction section to thoroughly discuss the development of numerical models in this field, explicitly state the difference from the 2022 paper, and highlight the novelty of this research. 2.Thresholds used to define the mixing zone (e.g., 16%–84% concentration range) are presented without a clear physical or statistical basis. The robustness of the results is questionable without a sensitivity analysis of these thresholds, as well as the numerical dispersion and grid resolution. Furthermore, the newly proposed metrics are not validated against any established, well-accepted measures of mixing from the fluid dynamics or hydrology literature, making their utility and interpretation difficult to assess. Response to comment #2: We will supplement targeted analyses and validation to address these concerns, and we elaborate on the revisions and scientific rationale as follows: 1) Actually, there is no standard for the threshold of the mixing zone; however, it generally accepted that the mixing zone ranges from 10% to 90% in terms of groundwater proportion (Hester et al., 2013; Woessner et al., 2000). This interval effectively distinguishes the mixing zone, where groundwater and surface water interact dynamically, from the two endmembers: pure groundwater (>90%) and pure surface water (<10%). In this study, we selected the 16–84% range, which corresponds to the ±1σ interval of a normal distribution, as proposed by Santizo et al. (2020). To assess how threshold variations influence results, we will perform a sensitivity analysis that tests three alternative concentration ranges: 10%–90% (wider interval), 16%–84%, and 20%–80% (narrower range) in the revised paper. 2) We agree that the numerical dispersion and grid resolution can influence solute transport and mixing results. We will test three grid sizes to assess spatial discretization impacts: baseline grid (the original resolution), fine grid and coarse grid. Critical parameters would be compared for each grid size, confirming that the baseline grid is sufficiently resolved to capture mixing dynamics and minimize numerical dispersion. 3) The calculations of mixing metrics in this study follow the methods outlined by Hester et al. (2013, 2014). These metrics are not newly proposed by this study; instead, they have been widely applied and validated in previous studies (Santizo et al., 2020, 2022; Nogueira et al., 2022; Woessner et al., 2000), the evidence that supports their rationality. We will state that clearly in the revised manuscript. 3.The model's idealizations limit the generalizability of its findings. Key simplifications include: (a) representing bedform migration as a simple translating sinusoidal head, which neglects morphodynamic feedbacks; (b) using a 2-D domain with a single bedform wavelength and amplitude; and (c) omitting sediment heterogeneity, such as grain size variations or layering, which are known to strongly influence hyporheic exchange. Response to comment #3: (a) Actually, simplifying undulating riverbeds to a horizontal configuration and applying a sinusoidal head distribution at the sediment-water interface exerts minimal impacts on both the hyporheic exchange flux (HEF) driven by the riverbed and solute transport within the riverbed, which has been proven by many researchers (Elliott & Brooks, 1997a, b; Eylers, 1994; Rutherford et al., 1995). Indeed, we have conducted an additional simulation incorporating ripple geometry and compared it to the flat-bed model: both yield similar HEF and solute plume characteristics (this simulation will be added to the revised manuscript). (b) Given the periodicity of bedform geometry in our model, we defined the model domain using a single bedform wavelength and amplitude. This is a well-established practice in modeling, as it minimizes edge effects without the computational intensity of full-scale simulations of multiple consecutive bedforms. To validate this simplification, we simulated three consecutive ripples, focusing on the middle one: comparisons of its pressure, solute concentration, and vertical boundary gradients with the single ripple model showed only minor differences. We will state that the model applies most directly to straight, low-curvature streams with uniform bedform spacing (typical of agricultural/urban downstream reaches, consistent with our focus; Hester et al., 2014; Krause et al., 2022). Here, hyporheic exchange is dominated by streamwise-vertical flow cells, with lateral hyporheic flux contributing little to total exchange. (c) We omitted this factor in the current study to isolating the effects of bedform migration on surface water-groundwater mixing and mixing-triggered reaction. This work specifically targets downstream gaining reaches characterized by fine sand beds, which are a relatively homogeneous sediment type. To address this limitation, we will add relevant context to the Discussion section of the revised manuscript. We will propose that future studies integrate stochastic K fields to explore how sediment heterogeneity interacts with bedform migration, for instance, whether high-K hotspots amplify or dampen migration-driven mixing.
Citation: https://doi.org/10.5194/egusphere-2025-2631-AC3 -
AC4: 'Reply on RC3', Zhang Wen, 25 Sep 2025
General comments: This study investigates the effects of migrating bedforms on surface-groundwater mixing and denitrification in a gaining stream reach. The authors employ a numerical model that couples Darcy flow with a translating sinusoidal head boundary and a multi-species reactive transport model. While the paper introduces new metrics such as “mixing flux,” “mixing fraction,” and “mixing-zone area,” my primary concern is the substantial overlap with the authors' previous work (2022). The modeling framework, numerical implementation, and validation appear nearly identical. The claimed novelty rests on applying the existing model to a different scenario and introducing new diagnostic metrics. However, the justification for these metrics is insufficient, and their introduction does not lead to the discovery of new physical mechanisms. Furthermore, the main conclusion: bedform migration shortens residence times and thereby reduces denitrification—reiterates, is a key finding from the 2022 paper (https://doi.org/10.1029/2022WR033258). The manuscript reads more as an incremental application of prior work rather than a standalone scientific article with a distinct and significant contribution. Response to the general comments: We thank this reviewer for the constructive comments on our manuscript. We think the main concern of this reviewer is the difference between this study and our previous publication (Ping et al., 2022). We will classify the difference and the novelty of this study clearly in the revised version. Actually, this study is quite different from Ping et al. (2022) because: 1) Firstly, in this study, we focus on downstream section of rivers under gaining conditions firstly to simulate surface water-groundwater mixing processes under dynamic bedforms, while Ping et al. (2022) focused on hyporheic exchange processes in the headwater to midstream sections of rivers. therefore, the research focus is quite different for these two studies. 2) Secondly, the mixing process has been considered in this study. As far as we know, this process has not been thoroughly investigated in the migrating ripple model (Jiang et al., 2022; Zheng et al., 2019; Kessler et al., 2015). It is notable that the bedform migration affect mixing process significantly. 3) Thirdly, Ping et al. (2022) only considered stream-borne solutes and non-mixing dependent reactions, while in this study we evaluate and quantify the self-purification ability of the hyporheic zone for groundwater-borne contaminant. Therefore, we believe that this study is quite different from Ping et al. (2022). In addition, we have obtained several interesting findings in this study, which can be summarized as: 1) bedform migration reshape the surface water-groundwater mixing zone and mixing-trigged biogeochemical zone patterns; 2) at slow to moderate migration rates, bedform migration expands the surface water-groundwater mixing zone and increases the fraction of mixing in hyporheic flux, with both parameters stabilizing at a plateau thereafter; 3) a large bedform migration celerity limits mixing-dependent denitrification rate, weakening the hyporheic zone’s role as a natural barrier against groundwater contaminant. In summary, we appreciate the concern from this reviewer on the novelty of this manuscript. We will revise the paper carefully to make it clear on the novelty of this paper, and clearly state the difference between this study and Ping et al. (2022). We are planning to revise literature review, the schematic diagram (figure 1), and some other parts, if necessary, to make it clearer. We are pretty sure that this concern can be addressed and this paper can be considered for potential publication in HESS after revision. 1.The literature review should be expanded to include a more thorough discussion of numerical model development in this field. Critically, this section must clearly and explicitly state how this study's model and objectives differ from and advance upon the 2022 paper. The work appears to be a case study applying a previously established model. It does not introduce a new mechanistic framework, a dimensionless parameter space map, or any analytical/scaling relationships that would constitute an independent and generalizable scientific contribution. Response to comment #1: We agree with this reviewer on this point. In the revised manuscript, we will make targeted revisions in Introduction section to thoroughly discuss the development of numerical models in this field, explicitly state the difference from the 2022 paper, and highlight the novelty of this research. 2.Thresholds used to define the mixing zone (e.g., 16%–84% concentration range) are presented without a clear physical or statistical basis. The robustness of the results is questionable without a sensitivity analysis of these thresholds, as well as the numerical dispersion and grid resolution. Furthermore, the newly proposed metrics are not validated against any established, well-accepted measures of mixing from the fluid dynamics or hydrology literature, making their utility and interpretation difficult to assess. Response to comment #2: We will supplement targeted analyses and validation to address these concerns, and we elaborate on the revisions and scientific rationale as follows: 1) Actually, there is no standard for the threshold of the mixing zone; however, it generally accepted that the mixing zone ranges from 10% to 90% in terms of groundwater proportion (Hester et al., 2013; Woessner et al., 2000). This interval effectively distinguishes the mixing zone, where groundwater and surface water interact dynamically, from the two endmembers: pure groundwater (>90%) and pure surface water (<10%). In this study, we selected the 16–84% range, which corresponds to the ±1σ interval of a normal distribution, as proposed by Santizo et al. (2020). To assess how threshold variations influence results, we will perform a sensitivity analysis that tests three alternative concentration ranges: 10%–90% (wider interval), 16%–84%, and 20%–80% (narrower range) in the revised paper. 2) We agree that the numerical dispersion and grid resolution can influence solute transport and mixing results. We will test three grid sizes to assess spatial discretization impacts: baseline grid (the original resolution), fine grid and coarse grid. Critical parameters would be compared for each grid size, confirming that the baseline grid is sufficiently resolved to capture mixing dynamics and minimize numerical dispersion. 3) The calculations of mixing metrics in this study follow the methods outlined by Hester et al. (2013, 2014). These metrics are not newly proposed by this study; instead, they have been widely applied and validated in previous studies (Santizo et al., 2020, 2022; Nogueira et al., 2022; Woessner et al., 2000), the evidence that supports their rationality. We will state that clearly in the revised manuscript. 3.The model's idealizations limit the generalizability of its findings. Key simplifications include: (a) representing bedform migration as a simple translating sinusoidal head, which neglects morphodynamic feedbacks; (b) using a 2-D domain with a single bedform wavelength and amplitude; and (c) omitting sediment heterogeneity, such as grain size variations or layering, which are known to strongly influence hyporheic exchange. Response to comment #3: (a) Actually, simplifying undulating riverbeds to a horizontal configuration and applying a sinusoidal head distribution at the sediment-water interface exerts minimal impacts on both the hyporheic exchange flux (HEF) driven by the riverbed and solute transport within the riverbed, which has been proven by many researchers (Elliott & Brooks, 1997a, b; Eylers, 1994; Rutherford et al., 1995). Indeed, we have conducted an additional simulation incorporating ripple geometry and compared it to the flat-bed model: both yield similar HEF and solute plume characteristics (this simulation will be added to the revised manuscript). (b) Given the periodicity of bedform geometry in our model, we defined the model domain using a single bedform wavelength and amplitude. This is a well-established practice in modeling, as it minimizes edge effects without the computational intensity of full-scale simulations of multiple consecutive bedforms. To validate this simplification, we simulated three consecutive ripples, focusing on the middle one: comparisons of its pressure, solute concentration, and vertical boundary gradients with the single ripple model showed only minor differences. We will state that the model applies most directly to straight, low-curvature streams with uniform bedform spacing (typical of agricultural/urban downstream reaches, consistent with our focus; Hester et al., 2014; Krause et al., 2022). Here, hyporheic exchange is dominated by streamwise-vertical flow cells, with lateral hyporheic flux contributing little to total exchange. (c) We omitted this factor in the current study to isolating the effects of bedform migration on surface water-groundwater mixing and mixing-triggered reaction. This work specifically targets downstream gaining reaches characterized by fine sand beds, which are a relatively homogeneous sediment type. To address this limitation, we will add relevant context to the Discussion section of the revised manuscript. We will propose that future studies integrate stochastic K fields to explore how sediment heterogeneity interacts with bedform migration, for instance, whether high-K hotspots amplify or dampen migration-driven mixing.
-
AC5: 'Reply on AC4', Zhang Wen, 25 Sep 2025
The directly uploaded text does not show paragraph breaks; however, the complete response is also accessible via this attachment.
Citation: https://doi.org/10.5194/egusphere-2025-2631-AC5
-
AC5: 'Reply on AC4', Zhang Wen, 25 Sep 2025
-
AC3: 'Reply on RC3', Zhang Wen, 25 Sep 2025
Peer review completion
Journal article(s) based on this preprint
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 1,278 | 194 | 52 | 1,524 | 73 | 43 | 64 |
- HTML: 1,278
- PDF: 194
- XML: 52
- Total: 1,524
- Supplement: 73
- BibTeX: 43
- EndNote: 64
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Yang Xian
Menggui Jin
Stefan Krause
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(1763 KB) - Metadata XML
-
Supplement
(237 KB) - BibTeX
- EndNote
- Final revised paper
The hyporheic zone plays a vital role in mitigating contaminants through interactions between surface water and groundwater. This study investigates the effects of bedform migration on the mixing dynamics of these waters and the accompanying biogeochemical reactions, particularly under gaining stream conditions, by utilizing a coupled hydrodynamic and reactive transport model. The research focus is meaningful, and the numerical model is well validated. However, I believe the overall clarity and presentation of the manuscript can be further improved. My suggestions are as follows:
Parameter explanations: I strongly recommend that the authors carefully review and expand the explanations of all parameters. Many parameters are currently not well defined, which may reduce accessibility for future readers.
Writing style: I suggest avoiding the use of "we" throughout the manuscript, as this expression is less common in formal academic writing. Rephrasing such sentences into more objective forms would enhance the manuscript’s style.
Specific comments:
Line 22: Please consider adding a comma between "migration" and "this" ("incorporates bedform migration this paper" → "incorporates bedform migration, this paper").
Notation consistency: I recommend keeping nitrate consistently as NO3^- throughout the manuscript (e.g., Line 81 uses NO3^-1, while Line 55 uses NO3).
Line 100: Please explain D50 here.
Line 114: Since the manuscript refers to "multi-component," I recommend using brackets to quickly introduce these components.
Line 125: Providing the dimensions of parameters (as done in Line 152) would be helpful.
Figure 1: The figure is somewhat confusing. Using a different color scheme for the mixing zone would improve clarity, as the current colors resemble a geological layer. Adding a note to indicate that the blue curve represents river elevation would also help. In the caption, I suggest rephrasing as: "NMD reaction = non-mixing-dependent reaction and MD reaction = mixing-dependent reaction."
Line 133: Please revise "'and SW out is …'" to "'SW out' is …."
Line 150: Since the bedform has a nonuniform distribution, would it be possible to consider anisotropic hydraulic conductivity in this study?
Line 153: As the abbreviation BC first appears here, it should be explained. (Since it also appears in Fig. 1 as line BC, I suggest not using BC as the abbreviation for boundary condition to avoid confusion.) Similarly, uc and dt in Eq. (2) require explanation. Given the large number of parameters, I recommend creating a notation table summarizing all variables with their definitions.
Line 161: Please clarify whether n and a have physical meaning or relationships. This would improve the organization of the manuscript.
Line 165: Explanation of Ri is missing.
Line 168: Explanations of vi and vj are missing. Also, note that v is a vector in Eq. (4) but a scalar in Eq. (5). Please clarify if they represent the same variable.
Line 171: Emphasizing why the three chemical species were selected would strengthen the manuscript (e.g., due to dominant reactions or research focus).
Line 178: Since Monod kinetics is referenced, I recommend directly citing Monod’s original work.
Line 190: In Line 150, you refer to hydraulic head, while here pressure is mentioned. Are these the same variable? Consistency is important. Additionally, please clarify the physical meaning of applying a periodic boundary condition.
Line 197: The phrase "a set of criteria" is unclear—are these criteria listed in the manuscript? Figure 2 is also difficult to interpret; more explanation would help.
Line 260: Just a positive note here—I appreciate how the manuscript justifies parameter values with references or physical reasoning. This strengthens the credibility of the work.
Figure 3: Including RMSE values (simulation vs. experiment) in each subplot would enhance the validation of the model.
Conclusion: The conclusion could be further strengthened by including quantitative results that highlight the main findings.