Summary of changes to the paper:

* Fixed multiple typos and grammatical errors.

* Added several paragraphs to clarify the points raised by reviewer 2 related to

(1) the type of model solved by QuaSoARe,

(2) the implication of imposing Lipschitz continuity to flux functions,

(3) the different nature of variable S (state variable) and O (diagnostic variable),

(4) the possibility to estimate O using quadrature methods.

* Added supplementary material in response to reviewer 2 that illustrates the two following points

(1) Application of QuaSoARe to flux functions that are not Lipschitz continuous,

(2) Use of an approximate quadrature method to estimate O.

---------------------------------------

Comment: “This is a very useful article reminding modellers of the numerical issues that need to be considered. It should be of interest to all modellers. Many non-linear ODEs cannot be solved analytically, particularly when combined into a system of ODE equations. It is models that use such ODEs that will benefit greatly from using QuaSoARe.”

Response: We thank the reviewer for their positive comments on QuaSoARe.

---------------------------------------

Comment: “The issue with hydrological models is that typically, the time step employed is dictated by the available data, not what is needed to ensure a sufficiently accurate estimation of the solution of the ODEs included in the model. The use of a coarse time step in the input data means a loss of information about what is happening at finer time scales, leading to uncertainty in the model output.”

Response: In order to precise this point, we added the following comment above equation 1:

“More precisely, let us denote the reservoir volume by  and assume that the reservoir is submitted to forcing variables   that remain constant over the time step. This assumption corresponds to most practical hydrological modelling scenarios where the time distribution of the forcings during the time step is unknown (e.g. time average of rainfall or potential evapotranspiration).”

---------------------------------------

Comment: “In Appendix A, an alternative approach would be to take the Taylor series approximation about the centre point of the interval. This would not ensure a match at the ends of the interval, leading to a likely discontinuity between intervals that would not be desirable. The approach taken ensures a continuous approximation of the function (noting that it will be discontinuous in the first derivative) and is effective providing the function is sufficiently close to a quadratic form. This will depend on the width of the interval and the variability of the function within that interval. The illustrative example used is S3, so a quadratic approximation would work well providing that the interval is not so large that the value of S ranges from near zero to a large enough value within the interval (see discussion of the illustrative example on page 8). This then defines the acceptable interval width that should be used. In general applications, this will depend on the form of the function f.”

Response: We agree with the reviewer regarding the fact that QuaSoARe performance will depend on the form of the flux functions and on the quality of the interpolation. This is why we highlighted this point in Section 4 of the paper among other limitations of the method. Based on a comment from reviewer 2, we expanded this section in order to highlight the issue of flux functions that are not Lipschitz continuous.

We hope that this section will suffice to flag potential pitfalls and their remedies when using QuaSoARe.

---------------------------------------

Comment: “Minor comments”

Response: We thank the reviewer for their careful inspection of the figures in our manuscript. All minor comments identified by the reviewer were fixed in the revised manuscript.

---------------------------------------

Comment: “Typographical/grammatical errors”

Response: We thank the reviewer for their careful inspection of our manuscript. All errors identified by the reviewer were fixed in the revised manuscript.

Summary of changes to the paper:

* Fixed multiple typos and grammatical errors.

* Added several paragraphs to clarify the points raised by reviewer 2 related to

(1) the type of model solved by QuaSoARe,

(2) the implication of imposing Lipschitz continuity to flux functions,

(3) the different nature of variable S (state variable) and O (diagnostic variable),

(4) the possibility to estimate O using quadrature methods.

* Added supplementary material that illustrates the two following points

(1) Application of QuaSoARe to flux functions that are not Lipschitz continuous,

(2) Use of a quadrature method to estimate O.

---------------------------------------

Comment: The paper is written in a very straightforward way and is relatively easy to follow with some familiarity in analytical and numerical integration (apart from the very first part where continuity assumptions are stated, see comments below). I am overall supportive of the approach which has a great potential for speeding up current formulations of many lumped models, but I still have questions regarding the kind of solutions which should be put forward to cure the problems raised by Clark and Kavetski about the numerical soundness of many such models, in the case of multistate model structures.

Response: We thank the reviewer for their overall support and excellent comments that have helped us clarified several important aspects of the manuscript.

---------------------------------------

Comment: I think that the introduction somehow lacks a brief presentation of the context in which unbalanced, cumulative outputs from model reservoirs could happen.

Response: We agree with the reviewer and added the following two sentences in the first paragraph:

“Such reservoirs are extensively used in rainfall-runoff models such as GR4J (Perrin et al., 2003), HBV (Bergstrom and Forsman, 1973), IHACRES (Croke and Jakeman, 2004) and SAC-SMA (Burnash and Ferral, 1981) where the reservoir dynamic is described by a differential equation relating the change in storage with inputs and outputs fluxes. If the model is applied to a time step long enough for storage to change significantly, this equation must be integrated to obtain the storage level at the end of the time step and the flux totals. However, apart from a few simple cases, there is no analytical solution to this mathematical problem, and one has to revert to numerical approximations (Clark and Kavetski, 2010; Kavetski and Clark, 2010). Furthermore, flux functions in hydrological models are often highly non-linear, which magnifies numerical errors when using inappropriate numerical schemes (Kavetski and Clark, 2011). This, in turn, degrades model simulation and calibration due to the extra-parameterisation needed to compensate for these errors (Kavetski et al., 2006).”

---------------------------------------

Comment: Restricting the analysis to such autonomous ODEs (Clark and Kavetski, 2010) is very common and is a lesser loss of generality than assumption (i), but it is still worth mentioning.

Response: We agree with the reviewer and added the following clarification above Equation 1:

“More precisely, let us denote the reservoir volume by  and assume that the reservoir is submitted to forcing variables   that remain constant over the time step. This assumption corresponds to most practical hydrological modelling scenarios where the time distribution of the forcings during the time step is unknown (e.g. average rainfall or potential evapotranspiration).”

---------------------------------------

Comment: My first question would be: in validating model outputs against observations (e.g. runoff, actual evapotranspiration, recharge, etc.), why wouln’t Eq. 3 be suitable for flux comparison, and why do we absolutely need to compute integrated quantities such as the O_i^′s (runoff depth over δ, evapotranspiration depth, recharge depth, etc.) which have dimension L (length), rather than to their instantaneous counterparts  (runoff rate, evapotranspiration rate, etc. with dimension L · T⁻¹) at successive times t,t + δ,t + 2δ, etc.

Response: We agree with the reviewer that it is possible to use an approximation of Eq 2 to compute the flux terms Oi. We thank the reviewer for flagging this critical point which we have now clarified by adding the following sentence above Eq 3:

“In addition, aside of applying the right ODE integrator to Eq. 1, solving the reservoir equation requires a numerical integration of Eq. 2 using a potentially different algorithm. For example, a simple quadrature method can be used to estimate the integral in the right-hand-side of Eq. 2 based on the two values S0 and S(δ). However, this approach can be highly inaccurate if the solution s(t) does not vary linearly with time as demonstrated in the supplementary material. A more accurate method is to expand Eq1 into a system of differential equations by adding one scalar differential equation for each flux.”

In addition, we added a section in the supplementary material to show a comparison between computing the flux from a quadrature approach, as suggested by the reviewer, and by integrating the extended ODE as initially indicated in our paper. This comparison suggests that the quadrature approach is most of the time accurate but can generate large errors if the solution s(t) varies non-linearly over the time step (e.g., during high flow regime).

---------------------------------------

Comment: Clearly, solving flux equations individually is a sufficient condition to fulfill the conditions (B), but I rather disagree on the fact that it is a necessary one, as stated in lines 60–70 (“Finally, jointly solving Eq. 1 and Eq. 2 using a numerical solver requires transforming the scalar equation Eq. 1 to a system of differential equations by adding one scalar equation for each flux”). If we have a procedure to solve the single, state-space formulated ODE in S, attributing residual errors in cumulative balance to any of the outputs would yield the same result at lower computational cost. If we chose an output O_i0 for balance error compensation, there are several possible fixes allowing to fullfill conditions (B):

Response: See previous response.

---------------------------------------

Comment: Using the language of atmospheric modeling we might say that the O_i’s are some kind of diagnostic variables that could be computed at the end of the simulated time period, S being the only variable truly belonging to the category of state variables. Adding this distinction in the formulation of the problem could improve the paper and help modelers think about the way their models are written (either discrete time, or state-space).

Response: We agree with the reviewer and thank them for this very useful comment. We incorporated the term “diagnostic variable” in a sentence at the end of section 1.1:

“It is important to note that the computation of Oi does not affect the solution S(t). To follow the terminology of Atmospheric Modelling,  Oi is a diagnostic variable (American Meteorological Society, 2024) whereas S is a prognostic variable.”

---------------------------------------

Comment: This would be my main question about a possible extension of QuaSoARe to the case of multistate models: can we avoid operator splitting, and if no, shouldn’t we still favor direct integration of the full, coupled system of ODEs describing the evolution of the state vector S ∈ R^m using the semi-implicit Euler (SIE) scheme for example?

Response:  The starting point of our paper, stated in the second paragraph of Section 1.1, is that direct integration of coupled reservoir systems using modern ODE integrators has not been adopted widely in hydrology. In this paragraph and in Section 1.2, we suggest that this is due to the complexity of the existing numerical methods and the additional burden they add to a model code. Motivated by this unsatisfactory situation, we developed QuaSoARe while making clear (in the abstract, at the end of the first and second paragraphs of section 1.1, at the beginning of section 4 and in the conclusion) that our method is limited to solving a scalar ODE.

We also acknowledged in the first sentence of Section 4, that it is preferrable to solve multiple stores jointly, hence following the reviewer suggestion. To make this point clear, we reworded the second sentence of Section 4 as follows:

“This is an important limitation of QuaSoARe as most hydrological models contain multiple stores that could benefit from a joint solution using a Runge−Kutta method, as presented by Clark and Kavetski (2010).”

---------------------------------------

Comment: The introduction of the paper quickly gets us into the swing of things with the mathematical assumptions stated p.2 (l. 36–41). It is a bit difficult to understand the practical implications of the Lipschitz-continuity assumption: does it mean that some functions are not admissible as flux functions?

Response: The reviewer is correct in stating that certain flux functions are not admissible as per our Lipschitz-continuity assumption. To clarify this point, we added the following paragraph in Section 4:

“Another limitation of QuaSoARe comes from the mathematical assumptions introduced in Section 1, and more specifically the need for flux functions to be Lipschitz continuous, which is equivalent to having a bounded derivative if the function is absolutely continuous. This assumption is required to ensure the unicity of the solutions of Eq. 1, and hence its monotonous nature as demonstrated at the beginning of Section 2.2 but eliminates many common flux functions encountered in hydrological systems (for example power functions of S with an exponent lower than 1). A first solution to this problem is to alter the flux functions to obtain smoother functions with bounded derivatives following Kavetski et al. (2006, see Section 5 of this paper) and, hence, revert to the domain of applicability of QuaSoARe. If this is not an option, QuaSoARe may still generate reasonable simulations if the solution of Eq. 1 does not come close to any discontinuity in the flux functions derivatives. This is, of course, case specific and subtle because it goes beyond the theoretical validity of QuaSoARe. An example of such a case is presented in supplementary material.”

---------------------------------------

Comment: If we consider the classical equation of a draining tank, (…) This function is not Lipschitz-continuous on the possible range of storage values (…). Does it mean that we have to check each flux function?

Response: The reviewer is correct is saying that all functions need to be checked for Lipschitz-continuity. However, as said in the response above, we believe that QuaSoARe could still generate reasonable predictions if the solution s(t) is not approaching the points where the derivative of the flux functions becomes unbounded. To illustrate this point we added an example in the supplementary material.