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.

The description of the hydrological model is not clear. Could the authors provide a brie description of what processes are included in the hydrological model? If there is no space in the main text, the authors should consider adding it to the supplementary text.

Thank you for the suggestion. We will clarify the description of the different model structures in the revised manuscript and strengthen the reference to the detailed online documentation available at https://smash.recover.inrae.fr. While it would be possible to replicate the full model structure in the supplementary material, we believe this is unnecessary given that the online documentation provides comprehensive and openly accessible information on the hydrological operators, which are classical but employed in an original way in our framework.

What is the assumption of the routing model in smash? There exist different types of routing models. For example, large-scale 1D river routing models are commonly applied at relatively coarse spatial resolution with the assumption that each grid cell contains a representative channel. A 2-dimensional routing model that solves kinematic, diffusion, or dynamic wave routing at high spatial resolution, such as less than 1km, does not need such an assumption. However, it may not be appropriate to apply such 2-dimensional routing model at coarse spatial resolution. Which type of routing model was implanted in smash? It has been applied at 3km resolution for CONUS and 500m resolution for France application. It would be helpful if the authors could elaborate on the routing method.

We will clarify in the revised manuscript that the routing module implemented in smash is a classical 1D conceptual kinematic wave model, applied over a D8 flow direction network. This reduces the routing problem to a 1D structure along predefined flow paths. The model is solved numerically using a linearized implicit scheme (Chow et al., 1998). Despite its simplicity, this approach has shown good performance and robustness in large-sample applications, both at 3 km resolution over CONUS and 500 m resolution over France. In addition to clarifying this in the main text, we will also elaborate on the current model’s limitations and potential future developments, particularly regarding hydraulic flood routing and inundation modeling.

Could the author report the computation time for the regionalization simulations vs. an individual forward run? This can help the readers to better understand the benefits of differentiable models. For example, in typical hydrological model calibrations, one has to run the model many times, e.g., 1000, with perturbed parameters, to obtain the optimal parameters. I think the differentiable model requires much fewer simulations to get the optimal parameter values, making its application to large scales more feasible.

We will report the computation times for the regionalization simulation, particularly for the most computationally demanding configuration using the ANN, as well as the corresponding forward run times. For reference, the full regionalization procedures (including calibration, input, and output operations) already take approximately 95 hours for the CONUS case and 31 hours for the Aude case. We will additionally run and report the computation times of the associated forward simulations, isolating the time spent on I/O and the pure forward model run.

There is no description of the calibrated parameters and their possible range. In addition, how are the parameters perturbed to solve the adjoint of the forward model in Figure 2 is not clear. Are the parameters perturbed at watershed level or are they perturbed at grid level within each watershed?

Thank you for pointing this out. We will include the calibration parameter ranges for all tested model structures in an appendix, in the form of a table listing parameter names and their respective ranges. A link to the online documentation will also be provided to clarify the meaning and role of each parameter, so as to avoid overloading the main text.

Regarding the computation of cost function derivatives with respect to the parameters, these are not obtained through perturbation methods. Instead, they are computed accurately using the numerical adjoint model, derived via automatic code differentiation as described in the manuscript. This approach is more suitable for high-dimensional and nonlinear models, where perturbation methods tend to be less reliable. We will make this explanation more explicit in the revised version.

Finally, the calibration can be performed in various modes: (i) uniform calibration with a single parameter set per watershed, (ii) fully distributed calibration with one parameter set per grid cell, or (iii) constrained calibration using physical descriptor maps for example, using binary masks over sub-basins to enable semi-distributed parameter configurations.

Line 61: Please provide the full name of MPR.

We will provide the full name of MPR.

Line 60 – Line 65: There are more large sample studies in terms of calibration of spatially distributed models using efficient statistical approaches. For example,

Yang, Y., Pan, M., Beck, H. E., Fisher, C. K., Beighley, R. E., Kao, S. C., Hong, Y., and Wood, E. F.: In quest of calibration density and consistency in hydrologic modeling: distributed parameter calibration against streamflow characteristics, Water. Resour Res., 55, 7784–7803, 2019.

Xu, D., Bisht, G., Sargsyan, K., Liao, C., and Leung, L. R.: Using a surrogate-assisted Bayesian framework to calibrate the runoff-generation scheme in the Energy Exascale Earth System Model (E3SM) v1, Geosci. Model Dev., 15, 5021–5043, https://doi.org/10.5194/gmd-15-5021-2022, 2022.

Hirpa, F. A., Salamon, P., Beck, H. E., Lorini, V., Alfieri, L., Zsoter, E., and Dadson, S. J.: Calibration of the Global Flood Awareness System (GloFAS) using daily streamflow data, J. Hydrol., 566, 595–606, https://doi.org/10.1016/j.jhydrol.2018.09.052, 2018.

We will expand the bibliography in this section.

Line 98: What type of mesh does smash support? Structured or unstructured?

smash only supports structured mesh.

Line 114 – Line 129: Are the snow, hydrological, and routing operators process-based functions or data-driven operators?

All operators in smash version 1.0 snow, hydrological, and routing are classical process-based conceptual models. We will clarify this point in the revised manuscript and include a reminder of the link to our detailed online documentation for further reference.

Line 159: Is equal weight used for w_g? Please specify how w_g was estimated for each gauge. In addition, how j_{reg} was estimated?

Thank you for highlighting this point. We acknowledge that this part was not clearly described in the manuscript and will clarify it in the revised version. In this study, equal weights were assigned to each watershed, which corresponds to minimizing the average of the individual cost functions. Additionally, no regularization term (j_{reg}) was applied in this work.

Line 215: Which year was chosen?

The evaluation period selected for this study was from July 31, 2010 to July 31, 2011. This information will be added to the revised manuscript for clarity.

Line 221: Please specify how the computation times were evaluated. According to my experience, it is unrealistic for a model to take only 0.1 seconds to run 365 time steps for 3846 grid cells. This relates to my previous comments: are the snow, hydrological, and routing operators process-based models?

The computation time reported in the manuscript corresponds to the average time for a pure forward or inverse run, excluding the time required for input and output processes. After re-evaluating the results, we confirm that this time reflects the computational cost of our relatively simple and efficient forward hydrological model, which involves only a small number of numerical operations per grid cell. Additionally, as mentioned earlier, the snow, hydrological, and routing operators are indeed process-based models.

Line 226: The better scaling in the inverse run is likely due to it being more computationally expensive than the forward run since the inverse run requires more computations.

We agree with your observation and this will be explicitly stated in the revised version.

Line 229: Why the routing scheme not be parallelized? There are a lot of global river routing models, and 2-dimensional hydrodynamic models are well parallelized.

You are correct, and we apologize for the confusion. The routing scheme is indeed parallelized, and we will correct this in the manuscript. What we intended to convey is that the routing scheme cannot be fully parallelized over the entire spatial domain, as it must be solved in a sequential manner from upstream to downstream. This distinction will be clarified in the revised version.

Figure 5: Please add longitude and latitude to the figure. Also, consider adding a colorbar to show what is plotted.

We will add both the longitude and latitude axes to the figure, as well as a colorbar to clearly indicate that the variable plotted is flow accumulation.

Line 268: Which 10 years are used for “warm-up”?

The "warm-up" period corresponds to the following:

• p1: from August 1, 1990, to July 31, 2000
• p2: from August 1, 1997, to July 31, 2007

Line 270: What parameters are calibrated?

A detailed description of the calibrated parameters and their respective optimization ranges has been added in the appendix for reference.

Line 275 – Line 280: Can you also report the default forward model performance?

Reporting the performance of the default forward model is not particularly meaningful in the context of this conceptual hydrological model. Our "default" model refers to a spatially uniform calibration of the parameters, which serves as a baseline for comparison.

Line 279: Please show the calibrated parameter values. Also, for the spatially uniform parameters, are they spatially uniform for the whole CONUS? Or are they only spatially uniform for each basin?

We will include statistics on the calibrated parameter values in the appendix. Regarding the spatially uniform parameters, they have been optimized independently for each watershed. As a result, we have a spatially uniform set of parameters within each watershed.

Line 279: I found the performance between spatially uniform and spatially distributed simulations is very close.

Indeed, for this daily and 3 km resolution model, the performance between spatially uniform and spatially distributed simulations is very close. However, we do observe a slight improvement in the KGE when using a spatially distributed calibration compared to a uniform one. At higher spatio-temporal resolutions, the performance gap becomes more pronounced, with uniform parameters yielding lower performance than a distributed calibration. This was clearer in a similar study conducted with an hourly and kilometric model over France (https://hal.science/hal-04989183/document).

Line 341: For the non-inertial shallow water models, do you mean diffusion wave routing?

The shallow water model we refer to includes the convective inertia term in the momentum equation, which allows for a simple finite difference solution.

Line 342- 344: I wonder if the authors can give some comments on how to implement the differentiable model to more complex land surface models or hydrological models?

Any hydrological or land surface model can be integrated into the smash framework, provided it is compatible with the constraints of automatic differentiation imposed by the Tapenade tool. In particular, the model must be written in a differentiable form and comply with Tapenade's supported syntax and structures. We provide technical guidance on how automatic differentiation is implemented within smash in the online documentation, in the section: Automatic Differentiation – Development Process Details (https://smash.recover.inrae.fr/contributor_guide/development_process_details.html#automatic-differentiation).

Line 368: Why change the evaluation metric from KGE to NSE in the application of France?

The change in evaluation metric from KGE to NSE in the France application serves two purposes: (i) to demonstrate that our advanced calibration and regionalization approaches can be used with different cost functions, and (ii) to focus on a case study with a relatively fine spatial and temporal resolution, where greater emphasis is placed on the high-water aspects. Hence, the use of NSE as both a calibration and evaluation metric in this context.

Line 403: I agree with the authors that the smash framework represents an advancement in data assimilation in the hydrological model. However, I disagree that the smash framework represents an advancement in the hydrological model. Based on the method section, smash implements simplified hydrological processes.

We agree with your observation and will revise the sentence accordingly. We believe that the smash framework represents an advancement in modular, regionalizable, differentiable numerical modeling, as well as in hydrological data assimilation.

François Colleoni and Pierre-André Garambois, on behalf of the authors.

References:

Chow, V.T., Maidment, D.R., Mays, L.W., 1998. Applied Hydrology. McGraw-Hill Series in Water Resources and Environmental Engineering.

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.

Data assimilation is not implemented or validated in this study, so I suggest removing it from the title and abstract.

The proposed framework adopts a variational data assimilation (VDA) approach, utilizing a numerically differentiable solver to enable accurate estimation of cost function gradients for solving the inverse problem. In the study, the inverse problem consists in searching parameter maps from discharge data, without or with descriptors-to-parameters regionalization mapping in the forward model (the so called Hybrid Data Assimilation and Parameter Regionalization (HDA-PR) approach Huynh et al. (2024)), over long time windows with the VDA algorithm. The smash framework also enables to perform VDA at shorter time scales and to infer initial states (not illustrated). These aspects will be clarified in the revised version, including a distinction between data assimilation for reanalysis and near real-time applications, and upgrades required for background and observation covariances matrices, aimed at inferring uncertain parameters—including initial states.

Introduction: Differentiable modeling for large sample and high-resolution hydrologic studies have seen notable progress in recent years. I encourage the authors to mention some related works to fully reflect recent progress in this area.

We will expand the bibliography in this section.

Line 118-120: What does GR stand for? It needs to be explained above.

"GR" stands for "Génie Rural," which translates to "Rural Engineering" in English.

Line 137: Typo “[l_k, u_k]”

Thank you for pointing that out. We will correct the typo in the revised manuscript.

Line 170: I think here should be the optimization/training problem rather than a typical variational data assimilation method, which typically uses near-real-time observations to update model states and forecasts. As mentioned above, data assimilation is not examined in this study.

VDA is traditionally used to infer uncertain or unknown model parameters and initial states from available data using gradient based optimization, which is enable by the proposed smash framework that incorporates a VDA algorithm. You are right, the primary focus of this study is on model parameters training in reanalysis, but using the VDA algorithm, rather than "real-time" data assimilation. This will be clarified.

Line 221-220: Is the model parallelized across grids using OpenMP? How many threads per CPU were used for parallel execution? Please clarify.

Thank you for your comment. Yes, the model is parallelized across the spatial grid using OpenMP. We made an error in referring to CPUs; in this case, only one thread is used per CPU. In the revised manuscript, we will consistently refer to threads rather than CPUs to avoid any confusion.

Line 233-234: The inverse run includes an iteration loop from the optimization algorithm, and within each iteration, there is a time loop. Is this correct? Please clarify.

There is some confusion in the manuscript between the terms "inverse run" and "inverse routine." To clarify, an "inverse run" refers to a single call to the adjoint model during the optimization process. The number of adjoint model calls within each optimization iteration may vary depending on the optimization algorithm used, such as L-BFGS-B. On the other hand, the "inverse routine" described in Figure 2 refers to multiple calls to the adjoint model. To resolve this confusion, we will maintain the term "inverse routine" in Figure 2 and rename "inverse run" to "adjoint run" in the 'Computational performance' section.

Figure 4: Instead of using the number of CPUS, it is more reasonable to use the number of cores/threads as x axis.

You are correct. We will make this adjustment in the revised manuscript.

The difference between the experiments in Section 4.2 and Section 4.3 is not very clear. Please clarify.

The main distinction lies in the "precision" of the modeling approach. Both sections use the same regionalization method, with differences in the watersheds selected for calibration and validation, as well as the cost function. In Section 4.2, we use global data over a large spatio-temporal scale and apply the regionalization algorithm at a country-wide scale. In Section 4.3, we focus on national data at a finer spatio-temporal scale, applying the method at the watershed level, which is more relevant for operational flood forecasting.

François Colleoni and Pierre-André Garambois, on behalf of the authors.

References:

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. URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2024WR037544, doi:,arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2024WR037544.