the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 05 Nov 2024
Technical note: Quadratic solution of the approximate reservoir equation (QuaSoARe)
Abstract. This paper presents a method to solve the reservoir equation, a special type of scalar ordinary differential equation controlling the dynamic of conceptual reservoirs found in most hydrological models. The method called “Quadratic Solution of the Approximate Reservoir Equation” (QuaSoARe) is applicable to any reservoir equation regardless of its non-linearity or the number of fluxes entering and leaving the reservoir. The method is based on a piecewise quadratic interpolation of the flux functions, which lead to an analytical and mass conservative solution. It is applied to two routing models and two rainfall-runoff stores that are representatives of hydrological model components and evaluated on six catchments located in Eastern Australia that experienced one of the most extreme floods in recent Australian history. A comparison of the method against two standard numerical schemes, the Radau fifth order implicit and Runge-Kutta of order 5(4) explicit schemes suggests that it can reach similar accuracy while reducing runtime by a factor of 10 to 50 depending on the model considered. At the same time, the model code is simple enough to be presented as a short pseudo-code included in our paper. Beyond solving a given reservoir equation, the method constitutes a promising avenue to define flexible models where flux functions are defined as piecewise quadratic functions, which can be solved exactly with QuaSoARe.
- Preprint
(1682 KB) - Metadata XML
- BibTeX
- EndNote
Status: closed
-
RC1: 'Comment on egusphere-2024-3184', Anonymous Referee #1, 22 Dec 2024
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.
As noted in the paper, numerical methods are needed in situations where an analytical solution is not possible. For example, linear ODEs that are commonly used in unit hydrograph components of hydrological models (e.g. IHACRES) can be solved analytically, so a numerical solution is not necessary. Some non-linear ODEs can also be solved analytically (e.g. the drainage equation in the IHACRES CMD module). 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.
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.
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.
Overall, the paper is scientifically very sound and relevant to the general modelling community. It is suitable for publication once the minor errors noted below are fixed.
Minor comments
- Line 142: should be “m-1 intervals”
- Figure 5: readers might not notice the scale factor on the top right of each panel. May be better to have this in the axis label on the right side – e.g. (1e-3 mm)
- Figure 6: the use of a shaded white font for the median values is difficult to read at times.
Typographical/grammatical errors
- Line 26: “a history”
- Line 33: “to bridge”
- Line 45: comma after “example”
- Line 58: “requires significant”
- Line 77: “, for example, the Saint-Venant …”
- Line 97: “large systems” or “a large system “
- Line 100: “that” instead of “which”
- Line 112: “Sections 2.2 and 2.3, which an accompanying Python …”
- Line 123: delete second “given”
- Line 149: “functions”
- Line 154: comma after “analytically”
- Line 165: “requirement”
- Line 207: comma before “such”
- Line 210: “as blue lines” – delete “in”
- Line 240: “If this is the case”
- Line 263: “that do not exist”
- Line 270: “conditions”
- Line 287: “outflow”
- Line 351: “For each catchment”
- Line 372: “quadratic functions”
- Line 411: “(500) that lead to”
- Line 421: “reach an error”
- Line 430: “polynomials”
- Line 445: “is satisfactory”
- Line 450: comma before “which” and another after “work” on the next line
- Line 454: “that” rather than “which”
- Line 470: “context, like most rainfall-runoff models, remain arbitrary"
- Line 473: “currently being explored”
- Line 483: insert a comma before “which”
- Line 484: “method’s” and “reservoirs”
- Line 486: “... is an order of magnitude smaller than the typical ..."
Citation: https://doi.org/10.5194/egusphere-2024-3184-RC1 -
AC2: 'Reply on RC1', Julien Lerat, 01 Feb 2025
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.
Citation: https://doi.org/10.5194/egusphere-2024-3184-AC2
-
RC2: 'Comment on egusphere-2024-3184', Anonymous Referee #2, 08 Jan 2025
The technical note proposed by J. Lerat introduces a new method for the integration of a scalar Ordinary Differential Equation (ODE), describing the evolution of a conceptual component such as a reservoir in a lumped hydrological model. The rationale for the development of the method is that, given this ODE formulated in state-space, it is difficult to compute together both the evolution of the single state variable (denoted S) over a timestep, as well as the different outputs from the system integrated over the timestep. The QuaSoARe method uses analytical (quadratic) substitutes for the instantaneous fluxes entering/leaving the reservoir.
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.
Please see PDF in attachment for detailed comments.
-
AC1: 'Reply on RC2', Julien Lerat, 01 Feb 2025
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 Oi′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−1) 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 Oi0 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 Oi’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 ∈ Rm 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.
-
AC1: 'Reply on RC2', Julien Lerat, 01 Feb 2025
Status: closed
-
RC1: 'Comment on egusphere-2024-3184', Anonymous Referee #1, 22 Dec 2024
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.
As noted in the paper, numerical methods are needed in situations where an analytical solution is not possible. For example, linear ODEs that are commonly used in unit hydrograph components of hydrological models (e.g. IHACRES) can be solved analytically, so a numerical solution is not necessary. Some non-linear ODEs can also be solved analytically (e.g. the drainage equation in the IHACRES CMD module). 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.
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.
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.
Overall, the paper is scientifically very sound and relevant to the general modelling community. It is suitable for publication once the minor errors noted below are fixed.
Minor comments
- Line 142: should be “m-1 intervals”
- Figure 5: readers might not notice the scale factor on the top right of each panel. May be better to have this in the axis label on the right side – e.g. (1e-3 mm)
- Figure 6: the use of a shaded white font for the median values is difficult to read at times.
Typographical/grammatical errors
- Line 26: “a history”
- Line 33: “to bridge”
- Line 45: comma after “example”
- Line 58: “requires significant”
- Line 77: “, for example, the Saint-Venant …”
- Line 97: “large systems” or “a large system “
- Line 100: “that” instead of “which”
- Line 112: “Sections 2.2 and 2.3, which an accompanying Python …”
- Line 123: delete second “given”
- Line 149: “functions”
- Line 154: comma after “analytically”
- Line 165: “requirement”
- Line 207: comma before “such”
- Line 210: “as blue lines” – delete “in”
- Line 240: “If this is the case”
- Line 263: “that do not exist”
- Line 270: “conditions”
- Line 287: “outflow”
- Line 351: “For each catchment”
- Line 372: “quadratic functions”
- Line 411: “(500) that lead to”
- Line 421: “reach an error”
- Line 430: “polynomials”
- Line 445: “is satisfactory”
- Line 450: comma before “which” and another after “work” on the next line
- Line 454: “that” rather than “which”
- Line 470: “context, like most rainfall-runoff models, remain arbitrary"
- Line 473: “currently being explored”
- Line 483: insert a comma before “which”
- Line 484: “method’s” and “reservoirs”
- Line 486: “... is an order of magnitude smaller than the typical ..."
Citation: https://doi.org/10.5194/egusphere-2024-3184-RC1 -
AC2: 'Reply on RC1', Julien Lerat, 01 Feb 2025
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.
Citation: https://doi.org/10.5194/egusphere-2024-3184-AC2
-
RC2: 'Comment on egusphere-2024-3184', Anonymous Referee #2, 08 Jan 2025
The technical note proposed by J. Lerat introduces a new method for the integration of a scalar Ordinary Differential Equation (ODE), describing the evolution of a conceptual component such as a reservoir in a lumped hydrological model. The rationale for the development of the method is that, given this ODE formulated in state-space, it is difficult to compute together both the evolution of the single state variable (denoted S) over a timestep, as well as the different outputs from the system integrated over the timestep. The QuaSoARe method uses analytical (quadratic) substitutes for the instantaneous fluxes entering/leaving the reservoir.
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.
Please see PDF in attachment for detailed comments.
-
AC1: 'Reply on RC2', Julien Lerat, 01 Feb 2025
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 Oi′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−1) 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 Oi0 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 Oi’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 ∈ Rm 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.
-
AC1: 'Reply on RC2', Julien Lerat, 01 Feb 2025
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 216 | 59 | 15 | 290 | 10 | 8 |
- HTML: 216
- PDF: 59
- XML: 15
- Total: 290
- BibTeX: 10
- EndNote: 8
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1