We greatly appreciate your time and effort in reviewing our manuscript and providing constructive and detailed comments. Below is our point-by-point response addressing your concerns, with reviewer comments in italic and author responses in bold.

"Major comments. Model comparison:  The authors compare a stand-alone GRU model and their hybrid model approaches. I think the comparison between these two models is necessary, valuable, and the results are presented clearly. However, the fact that the hybrid performs better than the stand-alone process-based model is expected. With the hybrid approach, you have a model with more degrees of freedom, and the embedded NN can compensate for structural deficiencies in the process-based part, which will increase performance."

We appreciate your feedback regarding the model comparison. Indeed, the hybrid model, for example GRNN.U has more degrees of freedom compared to the stand-alone GR.U model. However, it is important to note that the inputs and outputs of the flux correction model are physically consistent and of the same dimension as the original model. This design allows the hybrid model to learn non linearities in the internal flux laws, which we analyze thoroughly in the flux correction analysis in both time and space throughout the paper.

While the hybrid model (GRNN.U) does not necessarily have more conceptual parameters (maintaining the same number of reservoirs and connections here), it does introduce more non linearity in the internal flux laws corrections with neural network Φ₁. This added complexity effectively increases the model's degrees of freedom, contributing to its enhanced performance. 

This more complex hybrid model outperforms the original model in calibration, which is expected given its added complexity. Remarkably, it maintains robustness in both spatial and temporal validations, as evidenced by the numerical results.  The hybrid approaches have been rigorously tested through temporal and spatial validation, demonstrating their robustness and performance not only on calibration but also on various validation scenarios. Moreover, what sets our hybrid approach apart is its physical interpretability, which demonstrates its strength compared to pure machine learning or deep learning  approaches.

We will further clarify, analyze, and justify these aspects in the revision to address your concerns.

"What I think is missing, to have a better idea of where the hybrid stands, is a comparison with a purely data-driven approach. For example, having a stand-alone LSTM, trained regionally with lumped meteorological inputs (e.g., catchment average values) would be a good benchmark. Or use as inputs, not only the basin-averaged values of precipitation, temperature, etc,… but also include other basin-averaged statistics (mean, std, max and min) that you can compute from the gridded products. This way we can see how the hybrid approach performs against purely data-driven methods, and if the extra effort of going distributed is worth it."

Thank you for your suggestion. This work focuses on a spatially distributed conceptual model based on physics and its physical hybridization. Our emphasis is on the rigorous presentation and analysis of this framework over a large sample, its performance in calibration and spatio-temporal extrapolation, and the interpretability of internal fluxes corrected with the hybrid approach, as you mentioned in the previous comment that the comparison between the GR and the hybrid GRNN models is "necessary, valuable, and the results are presented clearly".

First, we believe that analyzing the question of lumped versus spatialized models is not within the scope of this study. A spatially distributed approach is essential given the high spatial variabilities involved within the flash flood-prone catchments of this dataset.

Second, unlike traditional and hybrid process-based models, which rely on physical conservation equations of mass, momentum, energy, and empirical closures, pure LSTM or ML models do not inherently impose physical constraints, leading to reduced physical interpretability and generalizability, especially under extreme or unseen hydrological conditions (Beven, 2020; Sit et al., 2020; Shen, 2018). Given this context, we believe that a comparison with a stand-alone LSTM is beyond the scope of this study, which focuses on the hybridization and physical interpretability of spatially distributed process-based models. Our goal is not to benchmark the hybrid approach against purely data-driven methods but to demonstrate the improvements achieved through the hybridization of a well-established, spatially distributed, differentiable numerical hydrological model. A benchmark represents the scope of another full study.

Moreover, building a stand-alone LSTM, operating at an hourly time step, trained regionally, and capable of accounting for basin-average forcings/descriptors and spatial information, would be interesting but would represent a significant undertaking based on our previous experience with daily LSTM models (Hashemi et al., 2022). This is the scope of another full study. Such pure deep learning models, with even more degrees of freedom than our relatively parsimonious hybrid model, are hardly interpretable or extrapolable beyond the training set - as other neural networks. Therefore, we believe that these two approaches are not comparable within this context, where the present paper features a rigorous presentation and detailed analysis, especially of optimized quantities and internal fluxes, over a large sample and in regionalization. We believe this design and study are scientifically solid, this will be better explained in the revised manuscript to answer such questions.

We appreciate your suggestion and will keep it in mind for future benchmarking studies that compare a broader range of modeling approaches.

"Section 3: Here you present two datasets, with which you run two sets of experiments. The first dataset includes 235 non-nested catchments in France with 13 years of data. In this dataset, you test the effect of having a NN for process parameterization. In the second dataset, you have 21 catchments in the Mediterranean region, both nested and independent, with 7 years of data. This one you use to test the model regionalization. Is there a reason why this last test cannot be made in the first dataset? One can evaluate regionalization from catchment to catchment, and not only inside the same catchment. Moreover, having results in 235 catchments for the second experiment will give more robust tests. Also, you can mix everything in a single dataset with 256 catchments.  I was just wondering why did you make this division?"

Thank you for your valuable feedback. We will provide clarifications on the complementary goals of the experiments on the two datasets, which should address your points.

For the dataset with 235 Catchments:

• Objective: Test the performance of ϕ₁ (NN for flux correction) using local calibrations only at downstream gauges. This setup demonstrates the efficiency of the ϕ₁ NN in improving model performance.
• Reason for separate testing: conducting a multi-catchment setup across the entire mesh of France is computationally challenging, given the high spatio-temporal resolution of the data and model (see Huynh et al. (2024), for details on computational costs of the adjoint model). Moreover this would require to investigate cost functions adapted for meaningfull information selection/weighting over a set of catchments with contrasted area and physics (see issues for regionalization with downstream catchments in Huynh et al. (2024)).

For MedEst dataset:

• Objective: Evaluate regional calibration in a multi-catchment setup. This dataset tests the performance of both ϕ₁ (flux correction NN) and ϕ₂ (regionalization NN).
• Reason for focused regionalization: Analyzing physical interpretability in a national multi-catchment setup is complex. For this initial study, we focused on regionalization performance within a specific and known study zone. It is worth noting that regionalization performance over a larger zone, covering approximately 1/4 of France, has already been studied in Huynh et al. (2024). Future studies can certainly explore a national multi-catchment setup, as you suggested.

By separating the datasets, we aimed to provide a clear and focused analysis of both the flux correction and regionalization capabilities of our hybrid model. We appreciate your suggestion and will consider expanding the scope in future research.

"Line 282-294: The differences between the models are quite small. For example, the difference shown in Figure 4 between the median NSE for the GR.U and the GRNN.U  is 0.008 and between the GRD and the GRNN.D is 0.014. Are the differences between the reported distributions statistically significant? I think this point should be further discussed. Because the hybrid approach has a higher flexibility than the process-based model. The embedded neural networks produce flux-correction parameters for each pixel and timestep, and if the differences between the hybrid and the stand-alone process-based model are small, it would be interesting to find out why. Maybe the physical dissipation of the basins makes it unnecessary to have so much detail, if one is just interested in the simulated discharge at a specific point. Or maybe the meteorological data is restricting further increases in quality."

We appreciate your insightful comments regarding the differences between the models and the significance of the reported distributions.

First, while the median improvements may appear small in case of temporal validation as you pointed out (knowing that evident higher performance of hybrid models observed in calibration), it is important to consider the entire distribution. In addition to the median values, we observe notable enhancements in other statistical measures, such as the interquartile range (0.25 and 0.75 quantiles) and whiskers in the boxplots.

Another important point is that for catchments that already exhibit satisfactory performance, the effect of hybridization is relatively small (leading to similar median, 0.75 and 0.95 quantile values). However, for poorly performing basins, the hybrid models provide substantial improvements, as evidenced by enhanced performance in the lower quartiles. To further illustrate this, an additional catchment-by-catchment comparison graph may be helpful, and we will consider adding this in the revised version. Note that even for slight global performance NSE improvement, the proposed method enables to learn spatio-temporal corrections of model internal fluxes within a basin/region, as shown by the analysis of the obtained internal flux corrections for each of the 235 catchments accross France. This property of the proposed approach is very promising to exploit other data, for improving evapotranspiration, infiltration, groundwater exchanges and other processes.

Finally, a particularly interesting result from this comparison is the performance of the hybrid model with spatially uniform control (GRNN.U) in multi-gauge regional calibration. Despite maintaining spatially uniform parameters, this model achieves performance comparable to the original model with spatially distributed parameters (GR.D). This result is remarkable, as it highlights the ability of ϕ₁ to leverage both the memory effects of ODEs and the NN’s capacity to learn spatially distributed forcings and learn spatio-temporal compensations of modeling rigidity/uncertainties, which is a very promising property in view of regionalization with more diverse information from multi-source data. 

We will include a more detailed comparison of these statistical indicators in the revised version to highlight the improvements more comprehensively.

"As an additional question, do the flux correction parameters allow the model to artificially increase/decrease the amount of water (violate the mass-conservation principle) in the control volume?"

Thank you for this valuable comment. The simulated water balance is influenced by the correction of k_exc, which represents the exchange flux and can result in either gains or losses relative to the original model (which is already non-conservative). The flux correction is illustrated in the graphs, showing its spatio-temporal variability, which makes it complex to analyze directly. We will add this clarification in the revised version.

"Line 295-300: You indicate that you are evaluating the performance of the model in flash floods, and then you evaluate it in 2700 events during the 6-year validation period. Are these 2700 events flash floods or just regular floods?  How did you classified them?"

Thank you for this comment. The selection of flood events is performed using an automatic segmentation algorithm based on Huynh et al. (2023), which identifies events where peak flows exceed a certain quantile threshold. Therefore, we acknowledge that these 2700 events are not necessarily flash floods but rather general flood events that include flash floods. To ensure clarity, we will remove the term “flash” from the text.

"Line 327-333: In these lines (and Figure 7) you compare the NSE for 143 flood events, indicating that the hybrid models perform better. Even if this is true, all the models performed quite badly. For the GR.U and GRNN.U the median NSEs are -0.48 and 0.09, which is a clear indication that the models do not work at all. Just taking the mean of the observed data would yield to a NSE of 0. For the other two models, the NSE did improve, but was still quite low (0.19 and 0.37). You should expand the discussion here and try to understand why all models are performing so badly."

We acknowledge that the NSE values computed for the 143 flood events are relatively low across all models. It is important to note that NSE for flood events, which are short time series with high values, is highly sensitive to small timing errors. Even slight discrepancies in peak timing can lead to substantial decreases in NSE. Additionally, data and modeling uncertainties may vary between events, making the accurate prediction of highly contrasted events particularly challenging.

The models are calibrated on the entire time series, and in this figure, we evaluate the validation results specifically for flood events, where classical approaches often struggle to accurately estimate water dynamics. This discrepancy highlights the difficulty in capturing the rapid and intense nature of flood events, even with advanced hybrid models.

We also recognize the need to investigate potential sources of error, including input data quality and model structural limitations, and the impact of using a calibration metric based solely on flood events, which could explain the overall challenges in flood event simulation. This will be clarified in the revised version.

"Minor comments."

We greatly appreciate your detailed minor comments and valuable suggestions. We will address each of these comments with a thorough explanation in the revision. Thank you again for your careful and constructive review of our manuscript.

Ngo Nghi Truyen Huynh and Pierre-André Garambois, on behalf of the authors.

References:

Beven, K., 2020. Deep learning, hydrological processes and the uniqueness of place. Hydrological Processes 34, 3608–3613.

Hashemi, R., Brigode, P., Garambois, P.A., Javelle, P., 2022. How can we benefit from regime information to make more effective use of long short-term memory (lstm) runoff models? Hydrology and Earth System Sciences 26, 5793–5816.

Huynh, N.N.T., Garambois, P.A., Colleoni, F., Javelle, P., 2023. Signatures-and-sensitivity-based multi-criteria variational calibration for distributed hydrological modeling applied to mediterranean floods. Journal of Hydrology 625, 129992.

Huynh, N.N.T., Garambois, P.A., Colleoni, F., Renard, B., Roux, H., Demargne, J., Jay-Allemand, M., Javelle, P., 2024. Learning regionalization using accurate spatial cost gradients within a differentiable high-resolution hydrological model: Application to the french mediterranean region. Water Resources Research 60, e2024WR037544.

Shen, C., 2018. A transdisciplinary review of deep learning research and its relevance for water resources scientists. Water Resources Research 54, 8558–8593.

Sit, M., Demiray, B.Z., Xiang, Z., Ewing, G.J., Sermet, Y., Demir, I., 2020. A comprehensive review of deep learning applications in hydrology and water resources. Water Science and Technology 82, 2635–2670.