General evaluation:

This paper presents “H2MV”, a Global Physically-Constrained Deep Learning Water Cycle Model with Vegetation. The proposed hybrid hydrological model extends a traditional physics-based hydrological model, combined with a static and a dynamic data-driven module. These modules learn temporal parameters (coefficients) as well as spatially-varying parameters. The novelty lies in explicitly representing vegetation constraints that play a relevant role in the water cycle, with the goal to achieve a more realistic and interpretable partitioning of evapotranspiration. The model is optimized against several observation products, and is tested for equifinality. The model code is publicly available.

Overall, the paper presents a valuable advance in hydrological modelling and is generally well written. However, the placement into existing literature could be more comprehensive. Further, the discussion of results is somewhat hard to follow, since the “analysis route” is not laid out beforehand, and many abbreviations destroy the flow.

With some minor edits, I believe this manuscript would qualify for publication with GMD. Please find my comments below.

Dear reviewer,

We appreciate your detailed review and your assistance in enhancing the clarity of our manuscript. Thank you for highlighting the general improvements needed, particularly in the literature review, methods, results, and the clarification of abbreviations. We will follow the reviewer’s suggestions to improve clarity in the revised version.

Best regards,

Zavud Baghirov on behalf of the authors

Specific comments:

1. l. 37: “Hybrid (or differentiable) modeling…” – I would see differentiable modeling as a specific case of hybrid modeling. Please be more specific. 

Thank you for pointing this out. We agree that the sentence is misleading and that hybrid modeling is not synonymous with differentiable modeling. We will clarify this and be more specific in the revised version.

2. l. 40 ff.: There are many more studies that demonstrate the use of LSTMs (or other networks) to “hybridize” physics-based hydrological models. Please provide a brief but broader overview of the literature before going into the specific model versions by Kraft et al. that you aim to improve here. (Part of Section 3.4 addresses this to some extent – consider moving this discussion of existing approaches to the introduction.)

We agree that conducting a broader literature review before discussing the specific study we build upon is indeed useful. We also agree that moving the discussion of existing approaches from Section 3.4 to the introduction will be beneficial. We will relocate the relevant parts of Section 3.4 to the introduction and broaden the review in the revised version of the manuscript.

3. Section 2.1: The introduction to the dataset is very short. Please mention at least which inputs are used, and provide a little context.

Thank you for pointing this out. We will extend the discussion of the datasets used in the revised version.

4. Eq. 2: Since fAPAR and all alphas are space- and time-dependent and are multiplied with each other, could you comment on the expected identifiability of your introduced parameters?

Thank you for this question. It highlights an important point---our introduced parameters are indeed susceptible to equifinality. Addressing equifinality is one of the biggest challenges in hybrid modeling, as it is in other approaches, including pure process-based modeling and pure machine learning. This is why we delve deeper into this aspect in the manuscript and quantify equifinality using a simple cross-validation (CV) based approach.

While we do not directly focus on the alphas in the manuscript, we do examine the parameters modeled using alphas. We quantify the variability of these estimated parameters across CV folds, which effectively assesses our model's sensitivity to both the initialization of neural networks' weights and the different training/validation sets for each fold. We find that some parameters exhibit high variability (indicating more equifinality), while others show relatively smaller variability. Please refer to Sections 2.4.1 and 3.3 for more details.

Regarding fAPAR predictions specifically, we do not expect equifinality, as it is directly constrained in the model against observations. We will add discussion about this in the revised version. 

5. Given the different alphas and their individual constraints, is there still a mass balance constraint over the entire hydrological model, or is this given up in the spirit of “local adjustments”?

Thank you for highlighting this important point. When we designed the conceptual hydrological model, we ensured that it fully adheres to the mass balance principle. The proposed hydrological model equations guarantee that the amount of water entering the system (e.g., rainfall) equals the amount of water leaving the system plus any change in storage within the system. This ensures that the neural networks (NNs) are also constrained to obey the mass balance.

Therefore, our predictions cannot be perfectly fitted to the observational data (e.g., terrestrial water storage), even when using highly data-adaptive NNs. We will clarify this in the revised version.

6. Section 2.4, “Model Evaluation”: It should be mentioned here that for model performance evaluation, RMSE, Pearson’s r and SDR are used.

Agreed. We will mention the performance evaluation metrics, including RMSE, Pearson's r, and SDR, in section 2.4 in the revised version as suggested by the reviewer.

7. Section 2.4.2: It is not clear why you are particularly interested in TWS, and a decomposition thereof. Please add a motivation here to guide the reader in your model evaluation efforts. In fact, after reading through the full manuscript, it seems this decomposition is never used. Omit this completely or include related results. Overall, the evaluation strategy should be much better explained here, see related comments below.

Thank you for the feedback on the model evaluation section of the manuscript. We will add more information and motivation to clarify our motivations, specifically regarding TWS decomposition. TWS anomaly observations from GRACE are one of our main constraints because these are direct satellite observations of an integrated hydrological state available globally and time-resolved. In our model, TWS is the sum of three main water storages that lack direct observations (except for snow water equivalent, which is derived from observational data). Therefore, it is important and interesting to assess which component of TWS is more dominant and where (spatially), in particular because previous studies highlighted large modelling uncertainties associated to this.

8. Also, in the results part, IAV plays a dominant role, but has never been introduced.

Thank you for pointing this out. We will introduce this and explain how we define interannual variability (IAV) in the methods section in the revised version of the manuscript.

9. It would make sense to bring Fig. B1 to the main body of the manuscript.

We agree that Figure B1 is an informative plot and deserves to be included in the main body of the manuscript. We will bring the figure to the main body in the revised manuscript. Thank you for the suggestion.

10. Fig. 5 is referenced before Fig. 3, so reorder the Figures, or reorder the storyline of the results Section (Fig. 5 is easier to understand after the discussion of Figs. 3 and 4).

Thank you for the suggestion. We will reorder the storyline of the results section to first reference and discuss Figures 3 and 4 before discussing Figure 5.

11. Fig. 3: How is the target obtained? Is it the mean value over all spatial domains contained in the testing set? Why is the variability of the target not shown? That would help judge the credibility of the model. And/or show individual spatial regions.

Indeed, the target in Figure 3 represents the observed mean fAPAR value over the testing set (spatial domain). The main purpose of this figure is to evaluate the general performance of our model on the testing set. We actually display our model's performance for fAPAR across the major regions in Figure C4.

We will clarify this in the figure caption in the revised version of the manuscript.

12. Fig. 8: The term “Equifinality index” is used here for the first time. Introduce in Section 2.4.1 or replace the term in the results part.

Thank you for pointing this out. We will define the term "equifinality index" in Section 2.4.1.

13. I would have expected a discussion of equifinality (Section 3.3) before Section 3.2 (interpretation of results with respect to emerging global patterns), and a judgement of these interpretations based on the findings about equifinality (how robust are the conclusions you draw in Section 3.2, given that some states are not perfectly constrained?).

We agree that it makes more sense to discuss equifinality before addressing the emerging global patterns. In the revised version of the manuscript, we will move the discussion of equifinality (Section 3.3) before Section 3.2. Thank you for the suggestion.

Technical comments:

1. l. 32: maybe “physics-based” instead of “physical” (they are still computer models…)
2. l. 41: “a … network” instead of “a … networks”
3. l. 230 “interannual” instead of “interranual” (several instances throughout the manuscript, and other misspellings such as “interannaul”)
4. l. 238: “observed patterns … are” instead of “is”

Thank you for the technical comments. We will follow the reviewer's suggestions in the revised version of the manuscript.

1. Scope

The article is within the scope of GMD.

2. Summary

The authors present the H2MV model, a further development of the H2M model (Kraft et al., 2022). H2MV is a hybrid model for the global terrestrial water cycle. It consists of a conceptual process based hydrological model including the main terrestrial stocks and fluxes of water, and connected neural networks (partly static, partly dynamic with memory) processing static catchment attributes and dynamic forcing to deliver space-variable, time-static catchment attributes (maximum vegetation-reachable soil water SMmax), space-time-variable vegetation states (fraction of absorbed photosynthetically active radiation fAPAR), and various space-time-variable parameters of the process-based hydrological model. Model forcing includes precipitation (P), radiation (Rad) and air temperature (T). Further observations used for model training are terrestrial water storage (TWS), fraction of absorbed photosynthetically active radiation (fAPAR), snow water equivalent (SWE) and runoff (Q). H2MV is trained in a Cross-validation (C-V) approach on 10 spatially mutually exclusive datasets, and validated on an additional spatial holdout set. Model performance is discussed for all predictive variables TWS, fAPAR, SWE and Q on various temporal aggregations (monthly, seasonal, interannual). The authors conclude that generally, model performance is acceptable and shows space-time patterns in agreement with expert expectations and the literature. Further, the authors discuss model equifinality, here expressed as the predictive variability of the target variables among the 10 C-V models. Here, the authors conclude that mainly soil-related parameters are uncertain, andthat model errors are dominated by phase shifts.

3. Evaluation

Overall, the work presented by the authors is an interesting and relevant contribution to global and surface modeling. The presentation style is mainly clear and complete, and the conclusion are supported by the results. So there are only minor revisions required to increase clarity and completeness before publication.

Dear Uwe Ehret,

We would like to thank you for your detailed review of our manuscript and your encouraging words about our work. We also appreciate your assistance in making our manuscript clearer and more readable.

Best regards,

Zavud Baghirov on behalf of the authors

L3 “. . . we explicitly represent vegetation states by the fraction of absorbed photosynthetically active radiation (fAPAR), and by the maximum soil moisture capacity (SMmax), . . . ”. This is misleading, as it suggests SMmax represents a vegetation state. From the rest of the text, I take it that SMmax is a spatially variable but temporally invariant representation of the maximum (vegetation-reachable) soil water content, i.e. a purely soil-related property. Therefore I suggest rephrasing with a better distinction between and explanation of the abiotic and biotic controls of soil water capacity.

Thank you for pointing this out. We agree that the sentence may cause confusion. Indeed, we define SM_max as spatially varying and temporally invariant. It represents the maximum storage of water that can be accessed and transpired by vegetation. Conceptually, it is dependent on both, physical soil properties, and effective rooting depth. We will clarify this in the revised version.

L8 The authors use the term ‘constrain’ throughout the manuscript to refer to observables used in an objective function during model training/calibration. As not all readers will be familiar with this use of terms, I suggest adding a related clarification, e.g. in L44.

Thank you for this suggestion. We agree that defining what we mean by “constraints” before using the term will help readers better understand our methodology. We will add a brief introduction of the term in the revised version of the manuscript.

L37 “Hybrid (or differentiable) modeling aims to address this challenge”. The sentence suggests that hybrid modeling is synonymous to differentiable modeling. This is not the case. There are hybrid modeling approaches that do not require differentiability of the process-based part, and not all differentiable model are hybrids. Therefore I suggest rephrasing.

We agree that hybrid modeling should not be used synonymously with differentiable modeling. We will be more precise here in the revised version.

L81 Please add a short information about the length of the available data.

Thank you for the suggestion. We will add brief information about the length of the datasets used in the revised version.

L165 C-V approach: The authors use 10 validation sets (and one common test data set common) that are mutually exclusive in terms of space, but not time. I understand that a time-exclusive CV approach for the validation sets may not be possible due to limited data, but, as high spatial correlation may exist between validation and testing sets for the same time, the testing set should be differing from the validation sets in terms of both space and time. This will help to better assess the models space-time generalization capabilities. My suggestion: Train the model on all but the last two years. Use all but the last two years for space-only CV testing in the same way as done now. Use the last two years for space-time independent testing.

Thank you for your suggestion. We will conduct the suggested experiment of holding out the last few years of data to test the model’s space-time generalizability and include the results in the revised manuscript.

L 174 Loss function: Eq. (9) calculates the total loss over all observed targets (TWS, fAPAR, SWE, Q). The targets come with different units, so their influence on total L might be different. How is equal weighting of each target in L assured? Is the loss calculated from the Z-transformed data as in Kraft et al. (2022)?

Indeed, different targets have different units, so we standardize them using Z transformation before calculating each individual loss and then summing all the loss terms. We will make this clear in the revised version of the manuscript.

L201 It is unclear to me what the authors mean by “each estimated process”. Please add this information to the text. Also, in L208, the authors refer to “parameters” rather than “processes”. Please clarify.

We apologize for the confusion. Throughout the manuscript, we had used “estimated processes” and “estimated parameters” interchangeably, and we agree that this needs clarification. In the revised version, we will consistently use “estimated parameters” for clarity.

Figs 3, 4, C1, C2, C3. Please add a x-axis label to plots a)

We will add the corresponding x-axis labels to the upper left plots in Figures 3, 4, C1, C2, and C3. Thank you for the suggestion.

L261 As H2MV is a further development of H2M, a performance comparison between the two is important. The authors provide this comparison in the Appendix in Fig. C5, but do not discuss it. In Fig. C5, it becomes apparent that H2MV performs worse than H2M in terms of at least

• RMSE for TWS, SWE, ET
• SDR for TWS, SWE

While I do not think that a new model generation needs to outperform a previous one for all metrics, the reader will benefit from a more detailed discussion of the performance differences between H2MV and H2M.

Thank you for highlighting this important point. Indeed, H2MV performs slightly worse than H2M on these parameters. However, we believe this is expected, as H2MV includes more process formulations than H2M, resulting in more 'physical' constraints (indirect regularization). In contrast, H2M is more flexible in adapting to the data, which likely accounts for the slight performance drop of H2MV compared to H2M.

We will add this discussion to Section C2, where Figure C5 is shown.

Parameter stability: I wonder how time-stable or time-variant the LSTM-predicted parameters of the hydrological model are. Ideally, if the hydrological model would fully contain all relevant processes, the parameters should be static. Time-variations would point at functional deficiencies of the hydrological model, and the time-patterns could point at the nature of these functional deficiencies. See e.g. Fig. 8 in Acuna Espinoza et al. (2024). I do not require that such a discussion is added to the current paper, rather it is a suggestion for further work.

Yours sincerely,

Uwe Ehret

We agree on this general consideration and that it deserves further attention. We are planning to address this in our future work.