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.

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 q_ex (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 q_n = n·q: c = c_ext (x, t) (Dirichlet) if q_n < 0 (inflow) and n·(D ▽ c) = 0 (zero diffusive flux) if q_n ≥ 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, U_b = 0.6, and D₅₀ = 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.

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.

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.