The authors build a decoupled hybrid architecture combining graph attention units, a differentiable Monod physics layer and a lightweight residual corrector, and integrate it into the mature CPqPy COMSOL-PHREEQC platform to resolve high computation costs and convergence defects of traditional geochemical solvers. Overall, this work presents an physically consistent surrogate workflow that reconciles computational speed with mechanistic accuracy for stiff biogeochemical reaction systems. However, several limitations merit further attention:

1. In the Methodology, P-GNN predicts dynamic k_max,j values, which are passed to the differentiable physical layer for Multi-Monod kinetic calculations. In field-scale long-term bioremediation, microbial biomass varies by multiple orders of magnitude: abundant substrates drive rapid growth, while toxic metabolites cause severe die-off. Since training data only covers limited concentration and biomass ranges, P-GNN may output unreliable k_max,j under out-of-distribution biomass conditions and induce numerical failures in the differentiable kinetic layer. Existing benchmarks only examine spatial scalability to 10⁵ grid nodes, lacking quantitative analysis of the surrogate’s long-term robustness against out-of-range biomass variations. It is recommended to add quantitative assessments of the framework’s extrapolation performance under extreme biomass fluctuations over long simulation horizons within the Discussion section.
2. Within the Methodology section, to mitigate long-term drift in predicted pH and pe values, the authors adopt a hybrid solving scheme that invokes full PHREEQC geochemical recalibration every ten timesteps throughout transport simulations. Fractured aquifers have highly heterogeneous flow and steep reaction fronts, making fixed recalibration intervals inefficient: smooth transport zones waste computing resources, while steep-gradient areas suffer delayed error accumulation. The authors are advised to discuss the feasibility of an adaptive recalibration strategy triggered by spatial residual variability or abrupt chemical gradients as a potential future improvement within the Discussion section.
3. The manuscript only reports results of the hybrid CPgPy workflow (PHYGNET + periodic PHREEQC recalibration), without long-term comparisons of three core benchmarks: original CPqPy, standalone PHYGNET, and the proposed hybrid scheme. No quantitative error curves of drift-sensitive variables including pH, pe and vinyl chloride over ≥500 h are provided, so the necessity and optimal value of the 10-step calibration cannot be verified. The authors are advised to supplement sensitivity tests with 5, 10 and 20-step recalibration frequencies, quantify the magnitude of long-term prediction drift under non-recalibrated conditions, plot comparative results of all three schemes in unified figures, and incorporate this analysis into Section 4.3 on computational efficiency.
4. The model is trained on 50000 single-step state transitions generated by PHREEQC (Δt = 1 h). However, when embedded in COMSOL for 240-hour RTM simulations, the surrogate performs autoregressive sequential predictions. The paper only reports spatial residuals at the final time step (Fig. 8), without quantifying how prediction errors evolve over time or whether they saturate or diverge.
5. The impressive 3524-fold speedup is a central claim, yet the manuscript lacks detailed specifications of the hardware (CPU/GPU models), software versions (PHREEQC, PyTorch, COMSOL), batch size for PHYGNET, and whether PHREEQC was run in serial or parallel mode. These details are essential for reproducibility.

This paper proposes a Physics-Constrained Graph Neural Network surrogate model (PHYGNET) to address the issues of low computational efficiency in microbial kinetic reactions, high computational costs of traditional geochemical solvers, and susceptibility to numerical convergence failures in existing Reactive Transport Models (RTMs). The method maps the microbial reaction network into a graph structure, introduces a Monod kinetics physical layer to ensure mass conservation, and incorporates a residual correction network to improve prediction accuracy. The model is further coupled with COMSOL to enable rapid simulation of complex reactive transport processes. Experimental results demonstrate that the proposed model significantly improves computational efficiency while maintaining prediction accuracy, and exhibits good interpretability. I have raised numerous comments below, including unscientific wording, unclear key steps that are not adequately explained, repetitive expressions, and contradictory statements throughout the manuscript.

1. Lines 38–40 and Lines 49–51: The content expressed in these two sections is highly repetitive.
2. Line 43: The reference of "The cost" is unclear and ambiguous.
3. Lines 49–73: The introduction section lacks clarity regarding the development trajectory of AI-based methods and the research motivation. In the current text, terms such as AI, machine learning (ML), neural networks (NN), and graph neural networks (GNN) are used interchangeably, but a unified conceptual definition and clear hierarchical relationship among these terms are missing. The manuscript first points out that traditional neural networks are difficult to apply to the simulation of complex reaction networks, and then directly proposes GNN as a solution. However, GNN itself is also a type of neural network, which makes the logical transition problematic.
4. Line 85: The phrase "not only a theoretically appealing" and the subsequent "powerful choice for surrogate modelling of biogeochemical reaction systems" appear to convey largely the same meaning. It seems that the intended emphasis here should be on the advantage of embedding physical mechanisms into the network. The wording should be revised to clarify this point.
5. Line 99: The phrase "kinetic logic of given biogeochemical system" is missing an article. It should be revised to, for example, "the kinetic logic of a given biogeochemical system" or "the kinetic logic of the given biogeochemical system."
6. Line 106: The authors only used two reaction networks to validate the model, so it is inappropriate to describe this as a "dual-validation approach." Additionally, the expression "ready extensibility to" seems problematic and should be rephrased for clarity and accuracy.
7. Line 120: The subscripts *m* and *f* should be typeset in mathematical italicsrather than in Roman or default text style.
8. Line 130: The phrase "which are" needs to be checked for grammatical agreement. Also, regarding "reaction terms (S)": is the term singular or plural? If singular, it should be "term"; if plural, it should correspond to multiple (S) entries. The current usage of "terms" with a single (S) is inconsistent.
9. Line 144: The abbreviation "e.g.," is used to introduce an example, but only one example is provided. If multiple examples are not given, "e.g.," may not be appropriate; consider using "i.e.," or rephrasing accordingly.
10. Lines 145–146: The use of bold and non-bold "C" to represent different parameters is highly confusing and prone to misinterpretation. In addition, the use of conjunctions (e.g., "and") is problematic, and there are excessive commas throughout the sentence. The sentence following Equation (5) that defines the variables is overly long and structurally unclear, making it difficult for readers to follow.
11. Line 155: The phrase "network-level regulatory feedbacks" is mentioned here, but it is unclear where this concept was previously introduced or defined in the manuscript. The authors should clarify its origin or provide an appropriate reference.
12. Lines 144–165: The definitions of the formulas in this section are highly disorganized. Different parameters are defined using the same letter symbols, and several parameters are defined repeatedly throughout the text. This creates significant confusion and should be thoroughly revised for consistency and clarity.
13. Lines 179–181: The meanings expressed in these two sentences are redundant and overlap with each other. The authors should consolidate or remove the repetitive content.
14. Line 188: The term "latent embeddings" is used here without any prior definition or explanation. The authors should provide a clear definition of what latent embeddings refer to in the context of this work.
15. Line 193: The output of the MPL (Message Passing Layer) is described as defining "the outputs (concentrations, pH, and pe)," yet earlier the prediction target was stated to be "(kmax,j)." This raises a critical inconsistency: which one is actually the output layer of the neural network? The authors need to clarify the exact output of the model.
16. Line 194: The terms "decoupled residual block" and "tasks" are mentioned but have not been defined or introduced anywhere previously. The authors should provide clear explanations for these terms.
17. Line 196: The manuscript states that P-GNN is "for reaction rate estimation," but the main text says it predicts "the kinetic parameters (kmax,j)." These two descriptions are not equivalent—reaction rate estimation and kinetic parameter prediction are different objectives. The authors should reconcile this inconsistency.In addition, there is no full name of “AI RTM”
18. Line 188: The term "full-physics recalibrations" is ambiguous—what exactly does it refer to? Does it mean that the MLP refines the outputs (concentrations, pH, and pe) by correcting residual errors, or does it refer to the enforcement of fundamental mass conservation and stoichiometric consistency? Neither of these interpretations seems to qualify as "full-physics" in the strict sense. The authors should clarify what they mean by this term or consider using more precise wording.
19. Line 197: The phrase "To balance physical consistency with predictive accuracy" implies a trade-off between the two objectives. However, the manuscript does not previously mention or justify why such a trade-off exists in this context. The authors should explain the rationale behind this balance and clarify the potential conflict between physical consistency and predictive accuracy.
20. Line 200: PHYGNET is described as employing a multi-task learning approach that simultaneously captures chemical kinetic processes and system states. However, the subsequent description of "joint prediction" as a "unified loss function" is conceptually problematic: joint prediction is an output/result, whereas a unified loss function pertains to the training/optimization process of the network. These are distinct concepts and should not be conflated.
21. Lines 201–204: The explanation provided here does not clearly articulate how the multi-task learning is actually implemented, nor does it specify how the balance between physical consistency and predictive accuracy is achieved. The authors should provide a more detailed and rigorous description of the multi-task learning framework and the specific mechanisms used to balance these two objectives.
22. Line 218: The term "every node" is not clearly explained. Since the graph contains both species and reactions, which entities correspond to species nodes and which correspond to reaction nodes? Furthermore, the specific organization of the heterogeneous graph has not been clearly defined; the manuscript has not previously stated that this is a heterogeneous graph neural network. The authors should further clarify the node types, edge definitions, and connection relationships within the graph structure.
23. Line 211: The text states "A final readout or MLP head," which raises the question: which one is actually the output layer? Why is it presented as an "or" (i.e., a choice between the two)? The authors should specify the exact architecture of the output layer.
24. Line 222: The term "global environment outputs" has not appeared or been defined previously. Also, what are the "relevant nodes" mentioned here—are they species nodes or reaction nodes? How exactly is the aggregation performed over these nodes?
25. Line 223: What is meant by "learned embedding"? This term has not been clearly defined in the manuscript.
26. Line 227: The text mentions "Encoding P-GNN into PHYGNET." However, PHYGNET already incorporates P-GNN as a core component—so what does it mean to "encode" P-GNN into PHYGNET? This phrasing is confusing and should be clarified.
27. Line 233: The statement that data-driven models "will violate mass conservation" is overly absolute. Not all data-driven models necessarily violate mass conservation; the authors should qualify this claim or provide supporting evidence.
28. Line 234: What are "localized kinetic parameters"? And what are "dynamically predicted parameters"? These terms are introduced without clear definitions.
29. Line 237: Equations (6)–(8) do not appear to include any explicit mass conservation formula. The authors should check whether the mass conservation constraint is actually represented in these equations.
30. Line 241: The text states that "this code" can encode "elemental mass," but the preceding paragraph refers to "stoichiometric mass." These two terms seem to refer to different concepts. The authors should reconcile the terminology.
31. Line 245: What are "updated environmental features"? This term has not been mentioned previously in the manuscript.
32. Line 246: The claim that the MLP can maintain mass conservation because the residual is computed for concentrations, pH, and pe raises the question: where do the ground-truth values for these residuals come from? Residuals are calculated from predicted values versus measured/true values. If such ground-truth data are available, then in theory, one could directly build a data-driven model using these data. The authors should explain why this is not the case and justify the use of residual correction.
33. Line 254: The text says "see Fig. B," but no such figure is referenced or appears to be included. Please check and correct the reference.
34. Line 262: The manuscript states that PHREEQC is used for correction. Does this correction apply to the reaction rates at the current time step? The subsequent text mentions correction based on the "latest COMSOL output," but PHREEQC typically simulates reactions over a certain time period. Is this period equivalent to the ten time steps mentioned? This is not clearly explained. Moreover, if PHREEQC is used for correction, this implies that PHREEQC is assumed to be correct. If so, why not simply use PHREEQC directly instead of the surrogate model?
35. Line 265: Where do the "dynamic reaction rates" come from? Can COMSOL actually output dynamic rates, or does it only output a fixed numerical value? Thisshould be
36. Line 328: The text says "For Case A (BTEX), The sampled spac"—this sentence appears to be cut off or incomplete. Please revise.
37. Line 333: The phrase "a localized batch simulation" is used in reference to PHREEQC. It is unclear what is meant by "localized batch simulation" in this context. The authors should explain how PHREEQC performs such simulations.
38. Line 336: The description of data partitioning is insufficient. The authors only mention sampling at fixed time intervals but do not specify how the data are split into training, validation, and test sets.
39. Line 342: The term "Contaminant concentrations" is used here. However, in the context of this study, are these truly "contaminants" or simply chemical species of interest? The authors should verify whether this terminology is appropriate.
40. Line 351: The authors have not adequately explained why a dual-objective optimization approach is necessary, rather than simply using a single loss function to complete the training. The rationale for this design choice should be clarified.
41. Line 355: It is incorrect to state that MSE "ignores" low-concentration species, as the MSE is by definition computed over all samples. What the authors likely intend to express is that the error contribution from low-concentration species is imbalanced (i.e., overshadowed by high-concentration species). Furthermore, this imbalance is not due to the continuity of groundwater dynamics, but rather because different species concentrations span multiple orders of magnitude, and the successive reactions lead to uneven concentration distributions of intermediate products. The explanation should be revised accordingly.
42. Line 362: What are the specific values or settings of the weights used in the loss function? The authors should provide detailed information on how these weights are determined.
43. Line 368: The figure caption appears to mention the split between training and validation sets for the first time, but no details are provided regarding how this split was performed. The authors should describe the data partitioning strategy clearly.
44. Line 372: The description of the relationship between the surrogate model error and the reactive transport model error needs clarification. Improving the prediction accuracy of PHYGNET is intended to reduce the error propagation from the surrogate model when coupled with the RTM, rather than reducing the modeling error of the RTM itself. The authors should revise this statement to avoid misunderstanding.
45. Line 400: The phrase should read "the concentration of CH₄ is" (grammatical correction).
46. Line 406: The statement that PHREEQC data show "values near or below analytical detection thresholds" may not be appropriate or accurate in this context. The authors should reconsider this wording.
47. Line 414: The summation symbol Σ should encompass the entire expression within the parentheses. Please check and correct the mathematical notation.
48. Throughout the manuscript: The spelling and capitalization of "ethene" should be checked and unified throughout the text.
49. Line 433: The word "Specially" appears to be used incorrectly. Did the authors mean "Specifically"? Please correct if applicable.
50. Line 439: The format "Figure: 5." is incorrect. It should be written as "Figure 5" or "Fig. 5" following standard formatting conventions.
51. Line 446: The phrase "are treated to be" is grammatically awkward. It should be revised to "are treated as" or "are considered to be."
52. Line 448: The phrase "capture sharp hydraulic gradients" should be checked for clarity and accuracy in context.
53. Line 485: The phrase "surrogate for the traditional solve" is missing an article. It should be revised to, for example, "surrogate for the traditional solver" or "surrogate for the traditional solution."
54. Line 534: What is meant by "knowledge distillation"? This term has not been introduced or defined previously in the manuscript. The authors should clarify its meaning and relevance in this context.
55. Lines 568–578: These two paragraphs appear to convey the same or very similar meaning, resulting in redundancy. The authors should consolidate or remove the repetitive content.
56. Line 589: The term "krate" appears to be used for the first time here. Has this variable been defined earlier? If not, it should be introduced and defined properly.
57. Line 621: The phrase "could alleviates" is grammatically incorrect—"could" must be followed by the base form of the verb. It should be corrected to "could alleviate."

This manuscript proposes PHYGNET, a GPU-accelerated surrogate model combining a physics-informed graph neural network and a residual MLP for nonlinear reactive transport simulations. The approach is interesting and potentially valuable. However, major concerns remain regarding benchmark fairness, temporal generalizability, long-term numerical stability, and manuscript organization. Detailed comments are provided below.

1. The manuscript compares GPU-based PHYGNET with CPU-based PHREEQC, but the benchmark lacks essential information, including the versions of COMSOL and PHREEQC, the coupling interface, CPU and GPU models, core counts, clock frequency, memory, and software environment. More importantly, PHYGNET appears to be compared with serial PHREEQC. This is not an equivalent benchmark because the USGS provides PhreeqcRM, which supports OpenMP and MPI parallelization for reactive transport simulations. A GPU-parallel versus CPU-serial comparison may substantially overestimate the reported speedup. The authors should provide complete software and hardware specifications and include comparisons with both serial PHREEQC and multicore PhreeqcRM. Wall-clock time, initialization, data-transfer cost, memory use, training cost, and inference cost should be reported separately.
2. All major experiments use a fixed time step of ( Δ t = 1 h). It remains unclear whether PHYGNET merely learns a one-hour state-transition mapping. Can the trained model be directly applied to other time steps, such as (0.5), (2), or (5 h), or would it need to be retrained?
3. The manuscript adopts a hybrid strategy in which PHREEQC is invoked every 10 time steps for physical correction. However, such periodic correction may constrain mass conservation and other physical requirements without eliminating the systematic trajectory drift caused by PHYGNET prediction errors. This external intervention may therefore mask the true accumulation of errors when PHYGNET is used as an independent surrogate. The authors should provide a comparison without PHREEQC correction, i.e., continuous prediction using pure PHYGNET, to evaluate its intrinsic numerical stability, long-term error accumulation, and the actual contribution of the periodic correction to prediction accuracy.
4. The sawtooth-like oscillations in (pe) during (t=10-60 h) in Fig. 8f suggest local overfitting and high-frequency numerical artifacts in this strongly nonlinear transition regime. However, the manuscript provides neither model performance across different datasets nor sufficient details of the training procedure. The authors should report the proportions and splitting strategy used for the training, validation, and independent test sets, together with quantitative metrics such as (R²) and RMSE for each dataset. Please also clarify whether overfitting-control strategies, such as smoothness regularization, early stopping, or cross-validation, were employed, and explain and correct the nonphysical oscillations observed in Fig. 8f.
5. In Fig. 8e, the prediction errors for major species, including (NO₃^-) and (O₂), increase noticeably after approximately 200 steps (200 h). However, the simulation is terminated at 240 steps, before the system reaches steady state. This early termination may conceal further error accumulation and long-term prediction instability. The authors should extend the simulation well beyond 240 steps until steady state or reaction completion is reached, explain the increase in errors after 200 steps, and evaluate the numerical stability of PHYGNET over longer prediction horizons.
6. The primary role of the final paragraph of the Introduction (Lines 85-96) should be to clearly delineate the study scope, experimental design, and evaluation objectives, namely, what the study does and aims to investigate, rather than prematurely presenting the results. I recommend that the authors restructure this concluding paragraph by clearly explaining the proposed model, the validation scenarios, and the evaluation criteria. Specifically, the paragraph should state: (1) what model is proposed, including PHYGNET and its core physics-constrained mechanisms; (2) the reactive transport scenarios in which it is validated, such as multispecies biogeochemical kinetic systems; and (3) the specific performance dimensions evaluated, such as long-term temporal stability, physical consistency, and computational speedup. This revision would allow the paragraph to serve as an effective roadmap for the remainder of the manuscript.
7. Lines 38-48 and 49-59 contain substantial overlap and redundancy. These passages should be merged and condensed to improve the clarity and conciseness of the manuscript.
8. The authors should move the methodological description of SHAP from Section 4.3 to Section 2 (Methodology), for example, around Line 204. Section 4.3 should retain only the SHAP results and their physical interpretation. This revision would improve the manuscript’s structural consistency and logical flow.