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.
- Preprint
(1763 KB) - Metadata XML
-
Supplement
(237 KB) - BibTeX
- EndNote
Status: open (until 14 Sep 2025)
-
RC1: 'Comment on egusphere-2025-2631', Anonymous Referee #1, 20 Aug 2025
reply
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
reply
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
reply
-
RC2: 'Comment on egusphere-2025-2631', Anonymous Referee #2, 31 Aug 2025
reply
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
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 714 | 44 | 17 | 775 | 22 | 15 | 23 |
- HTML: 714
- PDF: 44
- XML: 17
- Total: 775
- Supplement: 22
- BibTeX: 15
- EndNote: 23
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1