A Physically Integrated GNN Surrogate for Microbially-Mediated Kinetic Reactive Transport: PHYGNET

Wang, Jinbo; Zhang, Kunfeng; Illman, Walter; Chen, Shuai; Liu, Mingzhu

doi:10.5194/egusphere-2026-2509

Preprints

https://doi.org/10.5194/egusphere-2026-2509

Preprints

13 May 2026

| 13 May 2026

Status: this preprint is open for discussion and under review for Hydrology and Earth System Sciences (HESS).

A Physically Integrated GNN Surrogate for Microbially-Mediated Kinetic Reactive Transport: PHYGNET

Jinbo Wang, Kunfeng Zhang, Walter Illman, Shuai Chen, and Mingzhu Liu

Abstract. Reactive transport modelling (RTM) is essential for subsurface environmental management but is fundamentally constrained by traditional geochemical solvers. These solvers incur prohibitive computational costs and frequently suffer from numerical instabilities such as convergence failures, particularly in microbially-mediated kinetic reaction systems. While machine learning surrogates offer acceleration, they often lack physical consistency when dealing with stiff biogeochemical dynamics. Here we propose PHYGNET (Physically GNN Network), which maps microbial reaction networks into a directed graph by representing species and reactions as nodes and edges. It embeds Monod-type kinetics within a physics layer to enforce mass conservation and thermodynamic hierarchies, and incorporates a residual corrector for refinement. By successfully coupling with COMSOL, PHYGNET demonstrates the capability to execute full reactive transport simulations. Benchmark tests reveal that, in contrast to the severe super-linear time penalties faced by traditional solvers at engineering scales, PHYGNET maintains stable sub-linear scaling via tensor parallelism. At a scale of 10⁵ grid nodes, PHYGNET achieved an acceleration (up to 3524-fold) without numerical crashes. Furthermore, its escalating speedup ratio establishes a "small-sample training, ultra-large-scale inference" paradigm that effectively offsets initial data generation costs. Overall, PHYGNET provides an efficient and physically consistent framework for accelerating Monod-type microbial reactive transport simulations, offering a practical pathway for large-scale environmental applications.

Received: 30 Apr 2026 – Discussion started: 13 May 2026

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Download & links

Jinbo Wang, Kunfeng Zhang, Walter Illman, Shuai Chen, and Mingzhu Liu

Status: open (until 06 Aug 2026)

Post a comment Subscribe to comment alert

RC1:
'Comment on egusphere-2026-2509', Anonymous Referee #1, 25 Jun 2026 reply
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:
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.

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.

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.

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.

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.

Reply
Citation: https://doi.org/10.5194/egusphere-2026-2509-RC1
RC2:
'Comment on egusphere-2026-2509', Anonymous Referee #2, 10 Jul 2026 reply
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.

Lines 38–40 and Lines 49–51: The content expressed in these two sections is highly repetitive.

Line 43: The reference of "The cost" is unclear and ambiguous.

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.

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.

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."

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.

Line 120: The subscripts *m* and *f* should be typeset in mathematical italicsrather than in Roman or default text style.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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”

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.

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.

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.

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.

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.

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.

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?

Line 223: What is meant by "learned embedding"? This term has not been clearly defined in the manuscript.

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.

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.

Line 234: What are "localized kinetic parameters"? And what are "dynamically predicted parameters"? These terms are introduced without clear definitions.

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.

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.

Line 245: What are "updated environmental features"? This term has not been mentioned previously in the manuscript.

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.

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.

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?

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

Line 328: The text says "For Case A (BTEX), The sampled spac"—this sentence appears to be cut off or incomplete. Please revise.

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.

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.

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.

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.

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.

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.

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.

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.

Line 400: The phrase should read "the concentration of CH₄ is" (grammatical correction).

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.

Line 414: The summation symbol Σ should encompass the entire expression within the parentheses. Please check and correct the mathematical notation.

Throughout the manuscript: The spelling and capitalization of "ethene" should be checked and unified throughout the text.

Line 433: The word "Specially" appears to be used incorrectly. Did the authors mean "Specifically"? Please correct if applicable.

Line 439: The format "Figure: 5." is incorrect. It should be written as "Figure 5" or "Fig. 5" following standard formatting conventions.

Line 446: The phrase "are treated to be" is grammatically awkward. It should be revised to "are treated as" or "are considered to be."

Line 448: The phrase "capture sharp hydraulic gradients" should be checked for clarity and accuracy in context.

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."

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.

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.

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.

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."

Reply
Citation: https://doi.org/10.5194/egusphere-2026-2509-RC2

Jinbo Wang, Kunfeng Zhang, Walter Illman, Shuai Chen, and Mingzhu Liu

Viewed

Total article views: 275 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
211	46	18	275	14	16

HTML: 211
PDF: 46
XML: 18
Total: 275
BibTeX: 14
EndNote: 16

Views and downloads (calculated since 13 May 2026)

Month	HTML	PDF	XML	Total
May 2026	185	40	13	238
Jun 2026	16	1	2	19
Jul 2026	10	5	3	18

Cumulative views and downloads (calculated since 13 May 2026)

Month	HTML	PDF	XML	Total
May 2026	185	40	13	238
Jun 2026	16	1	2	19
Jul 2026	10	5	3	18

Viewed (geographical distribution)

Total article views: 259 (including HTML, PDF, and XML) Thereof 259 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 16 Jul 2026

Short summary

PHYGNET, a physical AI model, designed to accelerate microbial reactive transport simulations. This new approach is over 3,500 times faster than traditional methods and handles complex, large-scale simulations without crashing. By allowing "small-sample training, ultra-large-scale inference," it cuts weeks of environmental modeling work down to hours, making it practical for real-world cleanup and risk assessment projects.


Total:	0
HTML:	0
PDF:	0
XML:	0