This study presents a Bayesian data-driven framework that incorporates multiple physical constraints to correct structural errors in groundwater contaminant transport models. The two case studies suggest that, compared with the unconstrained and single-constraint approaches, incorporating multiple constraints can improve parameter identifiability, reduce predictive uncertainty, and enhance physical consistency. However, several issues need to be clarified or further addressed before the manuscript can be considered for publication. Please see the comments below.

1. In Bayesian inference, MCMC sampling generally requires the target posterior distribution to remain fixed during the sampling process. In this study, however, the weights of the sub-likelihood functions are dynamically updated during sampling, which may cause the target posterior distribution to vary with iterations and thus affect the validity of the Bayesian inference. The authors are therefore encouraged to discuss whether the proposed dynamic weight-updating strategy affects MCMC convergence. Appropriate convergence diagnostics should also be provided.
2. The physical constraint strictness, sliding window length, and minimum tolerance variance may directly influence the posterior parameter distributions and predictive results. The authors should provide a clearer justification for the selected values.
3. The manuscript states that incorporating multiple physical constraints makes the structural error simpler and easier to characterize using GPR. However, this conclusion is currently supported mainly by qualitative interpretation. Quantitative statistics of the structural-error term should be provided to support this claim.
4. The posterior parameter distributions show that some parameters still exhibit pronounced multimodality under the multiple-constraint scenario. This suggests that parameter equifinality or limited parameter identifiability may still remain. The authors should discuss the possible causes of this phenomenon and its implications for the reliability of parameter estimation.
5. Increasing the number of physical constraints does not necessarily improve prediction, and additional constraints may be redundant or even conflict with existing constraints. The authors are encouraged to discuss the principles for selecting appropriate physical constraints.
6. The proposed method is tested only using a laboratory sand column experiment and an idealized synthetic model. The authors should clarify in the Discussion the conditions under which the proposed method can be extended to real contaminated sites or watershed-scale applications, as well as its potential limitations.
7. Since the proposed method introduces additional sub-likelihood functions and dynamic weight updating, the computational cost should be discussed.
8. The values of porosity and bulk density appear to be inconsistent between the text and the table in Case Study 2. Please check and clarify these values.

This manuscript presents a data-driven framework for correcting structural errors in groundwater contaminant transport models. Gaussian process regression is used to represent systematic model discrepancy, while multiple physical constraints are introduced through a combined likelihood function. An adaptive weighting strategy and a stopping criterion are proposed to balance the observational likelihood and the likelihood terms associated with different physical constraints. The framework is tested using a synthetic three-dimensional PCE reactive-transport model with simplified preferential pathways and a Cd–phosphate sand-column experiment in which the cotransport mechanism is omitted from the physical model. The results show that using multiple physical constraints can improve validation-period predictions, reduce mass-balance errors, and narrow posterior parameter and prediction intervals.  The topic is relevant to groundwater modeling and uncertainty analysis. Model structural error is an important source of predictive uncertainty but is often treated less explicitly than parameter uncertainty. The main contribution of the study is the extension of a GPR-based structural-error correction method from a single physical constraint to multiple constraints, together with an adaptive weighting procedure. The two case studies represent different types of model inadequacy and provide a useful first demonstration of the method. The manuscript is generally well organized, and the results are encouraging. However, several methodological details and interpretations should be clarified before publication. These revisions do not require changing the main framework, but they would improve reproducibility and strengthen the conclusions.

Specific comments

1.Please clarify how adaptive weighting is implemented within the MCMC procedure.

The dynamic weighting strategy is central to the manuscript, but it is not fully clear whether the weights are updated independently for each Markov chain or shared among all chains. Please also explain how samples generated before and after weight stabilization are used in the posterior analysis. A practical approach would be to treat weight adjustment as an adaptive or burn-in stage, fix the weights after the stopping criterion is reached, and use only the subsequent samples for posterior inference. It would also be helpful to report the iteration at which the weights stabilize, together with standard MCMC diagnostics such as acceptance rate, R^, and effective sample size. 

2.The formulation and scaling of the physical constraints need further explanation.

In Eqs. (9)–(10), the constrained variable includes the physical-model prediction, GPR structural error, and measurement error. Since measurement error is not expected to satisfy mass conservation, the authors may consider applying the physical constraints to the latent corrected state, f(θ)+b(x,ϕ), rather than to the observation noise. Alternatively, the current implementation should be justified more clearly. Please also explain why the same constraint variance of 0.2 is used for constraints with different units and numerical scales. A short sensitivity test using one or two alternative values, or a description of how the constraint residuals were normalized, would be sufficient. 

3.Please provide a more complete assessment of predictive uncertainty.

The manuscript mainly interprets narrower 95% prediction intervals as reduced uncertainty. Narrower intervals are useful only when they retain adequate coverage. Please report, at least for the validation period, the coverage rate and average width of the 95% prediction intervals. An interval score could also be included if readily available. For the synthetic case, it would be useful to indicate whether the known true parameter values fall within the corresponding posterior credible intervals. These additions would distinguish improved uncertainty quantification from simple interval contraction. 

4.Several case-study details should be checked and clarified.

The measured bulk density of the sand column is reported as 1.69 g cm^{-3}, whereas Table 2 gives 1.24 g cm^{-3}. Please confirm which value was used in the model. The upper Cd concentration bound of 0.75 mg cm−3 also appears inconsistent with the stated influent concentration of 0.1 mM and may represent a normalized concentration; the unit and definition should therefore be verified.  In addition, the exponential and gamma priors in Tables 3 and 4 should be defined by their distribution parameters, not only by ranges. 

5.A simple comparison would help isolate the value of the adaptive weighting strategy.

The current scenarios compare no constraint, one constraint, and multiple constraints. It is therefore difficult to determine how much of the improvement comes from adding physical constraints and how much comes from dynamically updating their weights. Please consider adding one simple baseline using multiple constraints with fixed equal weights or fixed weights after residual normalization. This comparison could be limited to one representative case if computational cost is a concern. The discussion should also slightly moderate statements such as “general framework” and “physical constraints function analogously to additional observations.” More precise wording would be that the framework has been demonstrated in two cases and that physical constraints provide additional prior or regularization information. Because the remaining mass-balance errors are not negligible, “reduces violations of mass conservation” would also be more accurate than implying that conservation is fully satisfied. 

Overall, the manuscript presents a useful and promising approach. Addressing the points above would make the statistical formulation clearer, improve reproducibility, and provide stronger support for the conclusions.

This manuscript presents a physics-informed data-driven framework for correcting structural errors in groundwater reactive transport models. The topic is timely and relevant, and the proposed framework has the potential to improve predictive accuracy in groundwater models. I recommend publication after the following points are addressed.

- (Section 2.3) The authors introduced a dynamic weight-updating strategy to improve parameter identification and prediction accuracy when coupling multiple physical constraints. However, this is somewhat expected because additional physical constraints introduce more information into the inverse problem. Would similar improvements be obtained if fixed weights were assigned to multiple constraints? A comparison between constant-weight likelihood formulation may better demonstrate the necessity of the proposed adaptive weighting scheme.

- (Section 4.1) The physical constraints considered in the synthetic case are based on carbon and chlorine mass conservation, which can be accurately quantified in the example. In practical applications, however, physical constraints are often subject to uncertainty, e.g., boundary fluxes or source terms. Can the authors discuss how the proposed framework would handle situations when the physical constraints themselves are uncertain?

- (Section 5.1) It is expected that the posterior distributions contract compared to the case with single constraint or without physical constraint. I suggest plotting the reference distributions in Figures 4 and 5. This would allow readers to assess whether they converge toward the correct parameter values.

- (Section 5.4.1) The proposed framework iteratively updates the weights of the sub-likelihood functions. I have two questions regarding this procedure: (1) In case 1, the sub-likelihood values differ by several orders of magnitude. Does this imply that certain physical constraints are inherently more informative than the others? (2) Since one of the motivations of the framework is its application to practical problems, the computational cost of the adaptive weighting procedure needs to be addressed.

- (Section 5.4.3) This study does not compare the developed framework against other physics-informed machine learning approaches. As such, it is difficult to evaluate how much improvement comes from this model compared with existing approaches, such as PINNs. I understand implementing these methods is beyond the scope of this work; can the authors include a discussion comparing the proposed framework with existing physics-informed approaches?