In the paper entitled "A hybrid physics–ML framework for integrating groundwater dynamics into land surface modeling", the authors used a deep learning model to replace a physics-based solver in the ParFlow model. Overall this paper is well written. The structure is clear and description of the model is complete. However, there are still some issues to be clarified before the paper can be considered for a future publication in GMD.

1. In my understanding, because the training of the DL model is performed using the inputs and outputs of the simulations for the Pearl River Basin, China, the updated hybrid model is probably moslty suitable for this specific region. The authors may need to clarify it in the title and the abstract of the paper. Moreover, the name of the model  ParFlow may also need to appear in the title and the abstract as well.

2. L251, in my mind, the rest of WY2021 should be used for model testing instead of the full year, because the beginning of the year has been used for validating the model.

3. L360-361, the authors should use some numbers such as the computing time for one-month simulation, to indicate the speedup influence of the replacement more clearly.

General comment:

This manuscript develops a groundwater surrogate that combines static spatial information with dynamic hydrological inputs and is designed to replace ParFlow within a coupled land-groundwater modeling framework. The approach is demonstrated over a ∼34,000 km² basin in the Pearl River Basin, China, and achieves an approximately 20-fold speedup while maintaining generally good agreement with the physics-based reference model. The proposed design preserves the coupling interface with the land-surface model and therefore has potential value for efficient hybrid land-groundwater simulation.

However, the Results and Discussion section should be substantially reorganized. The authors should structure this section into clearly defined thematic subsections and provide more detailed discussion of the unresolved issues revealed by the results, particularly regarding spatial errors, long-term stability, and the advantages and limitations of the Hybrid model relative to the physics-based PF-CoLM system.

Detailed comments:

2 Methodology

2.1 Physics-based CoLM/ParFlow model

Line 118-120: The rooting-zone depth (3.43 m) and aquifer thickness (100 m) differ from those commonly used in previous ParFlow applications. Could the authors clarify the rationale for this vertical discretization and whether it was mainly designed to match the CoLM configuration?

2.2 ParFlow deep-learning surrogate model

Lines 176-177 and 234-235: The manuscript emphasizes that, unlike the original FSTR using meteorological forcing directly, the modified surrogate uses CoLM-derived root-zone net fluxes as external inputs. However, these fluxes are themselves responses to meteorological forcing processed through land-surface processes. The authors should clarify the advantage of using root-zone fluxes instead of meteorological forcing directly and how this flux-conditioned design differs from a conventional forcing-driven surrogate model. Does this modification mainly enable coupling with CoLM by preserving the ParFlow-CoLM interface, or does it also improve surrogate accuracy, stability, or physical consistency?

Lines 250-252 and 255-256: The surrogate is trained only on WY2019–2020 using 120 h temporal windows. The authors should justify the choice of these two years and clarify whether the training dataset sufficiently represents wet, dry, and extreme hydrological conditions.

Lines 250-256: The surrogate is trained exclusively using meteorological forcing and ParFlow simulations from the Pearl River Basin. It is therefore unclear whether the resulting surrogate represents a general approximation of ParFlow dynamics or is specific to the hydrological and spatial characteristics of the training basin. The authors should clarify this point.

Lines 267-268: Pressure heads are standardized using their mean and standard deviation for each vertical layer. Please discuss whether this standardization affects the representation of vertical hydraulic gradients, particularly under very dry or saturated conditions.

2.3 Additional architectural experiments

Lines 277-278 and 302-303: The purpose and added value of the temporal self-attention module are not sufficiently clear. The authors should clarify whether it is primarily intended to improve long-term dependency learning, prediction accuracy, autoregressive stability, or computational efficiency. If the improvement is limited, this section should be streamlined and greater emphasis placed on the core surrogate and coupling framework.

2.4 ParFlow surrogate-CoLM physics hybrid coupling

Lines 342-343: Given that the hybrid model is evaluated only on WY2021, please discuss whether a single evaluation year is sufficient to demonstrate long-term robustness across different hydroclimatic conditions.

Lines 358-360: Does the reported speedup apply only to the ParFlow-CoLM coupled simulation, excluding the ParFlow spin-up? If a full ParFlow spin-up is still required for application to a new basin, the overall computational benefit of the surrogate may be considerably smaller. Clarify this limitation and discuss the end-to-end computational cost.

3 Results and discussion

I do not recommend presenting the Results and Discussion primarily in the order of individual figures. In its current form, this section reads largely as a figure-by-figure description, making it difficult to identify the main scientific findings and distinguish the presentation of results from their interpretation.

The authors should reorganize this section into clearly defined subsections, with each subsection addressing a specific aspect of the results. The authors may determine the specific subsection structure according to the content. For example, Lines 373-391 could be consolidated into a subsection focusing on differences in spatial patterns between the hybrid and physics-based simulations. Such a structure would substantially improve the readability and scientific focus of this section.

The manuscript evaluates WTD, soil moisture, and latent heat flux, but runoff is not discussed. Please clarify why runoff is excluded and include an evaluation of runoff performance, given its importance in assessing the coupled hydrological system.

Figure S2 and Lines 347-348: The hybrid model is trained five times using different random initializations. Clarify what the labels in Fig. S2a represent and what specifically differs among these five realizations.

Line 356: The checkpoint with the lowest validation RMSE was selected for subsequent evaluation. Identify the selected checkpoint and report the corresponding random initialization and training epoch/iteration. Please also clarify how sensitive the reported evaluation is to the checkpoint selection.

Lines 378-380: The authors attribute the finer-scale heterogeneity to soil texture and land-surface properties. However, if the hybrid model exhibits stronger or finer-scale spatial variability than the physics-based simulation, the authors should clarify why this occurs and whether these additional spatial features are physically meaningful or artifacts introduced by the surrogate.

Line 385: Large errors appear to be concentrated in the southeastern part of the domain. The authors should provide a more detailed discussion of the possible causes of this spatial error pattern.

The Results and Discussion should focus more clearly on the performance gain of the hybrid model relative to the original PF-CoLM system. The central contribution of this manuscript is the development of a surrogate model to replace or accelerate ParFlow within the coupled framework. Therefore, the key comparisons should be hybrid vs. PF-CoLM and hybrid vs. observations/reference datasets.

Lines 412-419 and 443-449: Considerable emphasis is currently placed on comparisons between PF-CoLM and the PFFD/no-lateral-flow configuration. While PFFD can serve as a useful secondary benchmark for illustrating the effects of simplified groundwater treatment, these comparisons do not directly demonstrate the added value of the hybrid model and should not dominate the main Results and Discussion.

I therefore suggest substantially reducing the discussion of PFFD vs. PF-CoLM and retaining PFFD only as a secondary reference case.

Figure S4d: I do not recommend including the Spearman-correlation comparison between the free-drainage and lateral-flow ParFlow configurations. Since both are physics-based ParFlow configurations representing different subsurface treatments, this comparison provides limited information for evaluating the hybrid model.

Lines 421-423 and 438-440, Figure 7: The authors should provide a quantitative discussion of the latent heat flux errors shown in Figs. 7b and 7c, particularly for the comparison between hybrid model and ERA5-Land. The key performance metrics and their values shown in Fig. 7 should also be reported explicitly in the text to allow readers to assess model performance without repeatedly referring to the figure.

Lines 450-452: Although the hybrid model provides substantial computational speedup, its improvement in simulation performance appears relatively modest. For example, the physics-based model already shows a high KGE against ERA5-Land for latent heat flux, with only a small improvement from hybrid. The authors should clearly demonstrate and discuss the specific advantages of hybrid over the physics-based model beyond computational efficiency.

Figure 9: The river-network-related spatial patterns remain clearly visible in the physics-based simulation but are much weaker in the hybrid results, including the bias and correlation maps. Could this indicate a deficiency of the hybrid model in reproducing the river-network features captured by the physics-based model? The authors should explicitly discuss this issue and explain its possible causes.

Lines 481-485 and Figure 8: Figure 8 appears to provide limited additional information, as the upward propagation of subsurface-state errors into surface fluxes has already been demonstrated and discussed earlier. The authors should either clarify the distinct contribution of this figure or consider removing/streamlining it to avoid repetition.

Lines 510-511: Large hybrid-Physics discrepancies in the southeastern region are mentioned repeatedly, but the manuscript does not provide a detailed explanation for this persistent spatial error pattern. The authors should discuss the possible causes of the larger errors in this region.

Line 525: The progressive deterioration is attributed to heavy precipitation, but the supporting evidence is unclear. Please specify whether this interpretation is based on Fig. 6 or Fig. 10 subplot a and cite the corresponding figure. A spatial precipitation map should also be provided to further examine whether regions of large hybrid errors coincide with areas of high precipitation.

Figures 5, 9, and 11: In the full-year simulation shown in Fig. 9, the hybrid model does not reproduce the river-network features in soil moisture as clearly as the physics-based model. The same behavior is also evident in Fig. 5, whereas Fig. 11 shows much clearer river-related WTD patterns during a selected short-term heavy-rainfall event. The authors should discuss the possible reasons for this difference between long-term and short-term predictions.

Figures 6 and 10: For the same period (hours 5761-6480), the median WTD bias is approximately 0.3-0.4 m in the continuous full-year simulation (Fig. 6) but remains below about 0.1 m in the segmented 720-h prediction (Fig. 10). Please explain why the same period produces such different error magnitudes under the two prediction methods.

If Fig. 9 represents a long-term simulation and Fig. 11 represents a short-term heavy-rainfall event, the results seem to suggest that the hybrid model loses clear river-network features during long-term simulation but can reproduce them during short-term high-rainfall conditions. The authors should discuss whether this indicates that the hybrid model is more suitable for short-term heavy-rainfall simulations than for long-term hydrological simulations.

Lines 573-574: I agree that the hybrid model appears promising for short-term flood simulations, as it provides substantial computational speedup while reproducing river-network features and overall model performance comparable to the physics-based simulation. However, the current results do not convincingly demonstrate its suitability for long-term or decadal simulations. Long-term application requires substantial training data, and the hybrid model appears to lose clear river-network features, whereas the physics-based model maintains good performance. Under these conditions, the main demonstrated advantage of the hybrid model appears to be computational efficiency. The authors should discuss this limitation more explicitly.

Conclusion

The hybrid model is developed and evaluated only for the Pearl River Basin. The authors should clarify whether the trained surrogate model can be directly applied to other basins, or whether new ParFlow simulations, training data, and retraining would be required for each new domain. The manuscript should also distinguish between the general applicability of the proposed framework and the transferability of the trained surrogate model.