The paper entitled "A simple, efficient, mass conservative approach to solving Richards’ Equation (openRE, v1.0)" outlines a straightforward implementation to solve the one-dimensional Richards' equation using off-the-shelf ODE solvers, with a novel amendment to effectively track the cumulative mass flux through the boundaries. The approach is rigorously compared to approaches and test cases from the literature. The paper is very well-written and structured. The contribution is somewhat novel (I have colleagues teaching solution of the advection dispersion equation using method of lines with basic ODE solvers at the undergraduate level; this is not a super-new idea), but the degree of rigour in assessment of the various libraries, tolerance and time step choices, and introduction of the SFOM flux tracking method puts this into the range of publishable contribution for a technical note in GMD. 

Some minor nitpicking comments that the authors may want to consider here

1) the use of Q_j->j+1 (introduced in eqn 18) seems like subscript overkill - why not just Q_j ?

2) it would be useful to report the domain extent and model simulation duration for Mathias' solution in section 3.1.3 (these are implicitly in the figure, but would provide a more complete problem statement in the text)

Some other things to consider in the future:

1) I envision the method of lines may perform even better in relation to other methods for cases with non-constant space steps and layering of different media. It would have been nice to see a case study in this vein, but I would by no means require it here. Just something worth toying around with.

2) The use of arithmetic mean for calcualting hydraulic conductivity for the 1-D problem struck me as strange - the effective resistance to flow is typically treated using the harmonic mean for such problems by default (and this is well-documented even in the source they provided).

3) It would be very interesting to see how this approach performs in the more relaxed domains simulated in land surface schemes, with inherently much larger space steps by default. That is, you are looking at the perfect limits against analytic solutions, but how does this approach do 'in the trenches' for practical problems where we can't afford the burden of 0.001s time steps and 0.0025 m space steps? Is it worth the effort of deploying for these types of problems?

I have been reviewing papers for 20 years and this is only the second initial submission where I have recommended acceptance 'as is'. Thanks for making my job as reviewer easy.

The manuscript focuses on numerical solutions to Richards equation, which is the de facto standard equation for simulating flows through variably saturated soils. 

Richards equation has long attracted the attention of hydrogeologists/modellers with interest in numerical solutions, hence there is a rich history of publications on numerical (as well as analytical and semi-analytical) solutions for this equation. The integration of Richards equation in time using "time stepping" schemes is one of such common research topics.

Overall I enjoyed reading the manuscript and its topic fits well within the scope of GMD. It is great to see open source tools becoming more and more common for geophysical models. The implementations in Python are generally instructive (but see comments below).

Unfortunately I have several major concerns with the current version of the manuscript:

1. The contribution of the study is stated in a confusing and potentially misleading way.

2. The benchmarking reported in the manuscript is carried out against common numerical schemes dating back to the 90s. However it is not acknowledged that many existing publications have already reported new (at the time) schemes that substantially outperform these common "old" schemes. Hence the claim of efficiency needs to be heavily qualified.

3. The theory section is unfortunately not of high quality and makes numerous inaccurate / confusing statements. The description of at least one of the benchmarked methods is clearly incorrect (referred to as "no iteration" when the cited reference uses an iterative algorithm)

These concerns are elaborated below (see "Detailed comments").

I see the real contributions of this work as follows:

- the new method to quantify the boundary fluxes that reduces water balance errors without negative impacts on model runtimes. I like the idea of adding an extra mass balance term to the vector of unknowns for the ODE solver to keep a more accurate track of the mass balance fluxes. This is a neat idea for off-the-shelf methjds.

- an open source implementation of common algorithms for (1D) Richards equation. This is indeed very nice for learning how these algorithms operate. However here please double-check/correct your implementation of the pressure-form solution against the description in the cited Celia et al 1990 publication

- benchmarking of these common algorithms against selected alternatives. This contribution is also useful but should be more carefully delineated. First the benchmarking using CPU time depends heavily on how much care is taken to optimise the code. Second the methods included in the comparison are indeed representative of common methods for Richards equation, but many methods (reported in the literature) that outperform these common methods are not included in the comparison. This is a pretty pertinent piece of information to be noted to avoid confusing or even misleading readers.

I recommend a major revision to address these concerns, and rebalance the manuscript towards the new findings.

-----

Detailed review comments

1. Contribution

The results presented upfront as the major findings are already well known and established in the literature on Richards equation.

For example, consider the first sentence of the abstract, which states "We show that a simple numerical solution procedure – namely the method of lines combined with an off-the-shelf ODE solver – can provide efficient, mass conservative solutions to the pressure-head form of Richards’ equation".

In the Richards eqn literature this result was established as far back as 1997 by Tocci et al who indeed used the method of lines, an off-the-shelf ODE solver (DASPK) - and applied it to the pressure form to achieve good mass balance and high numerical accuracy. 

From the Abstract of Tocci et al's paper: "We show how a differential algebraic equation implementation of the method of lines can give solutions to [NB: pressure-form] RE that are accurate, have good mass balance properties, explicitly control temporal truncation error, and are more economical than standard approaches for a wide range of solution accuracies".

There is also the follow-up paper by Miller et al (1998) which further developed these results. The studies by Tocci et al and Miller et al are well cited and amongst the classics for the time integration of Richards equation. So this is by no means some poorly recognised result being brought to light in the new manuscript.

It would be important to clarify how the submitted manuscript contributes beyond this existing knowledge.

The superior mass balance of the Modified Picard method (applied to the mixed form of Richards equation), which is stated as another key finding in the Abstract (line 17), is also not a new result. The Modified Picard method is by construction mass conservative, as designed by Celia et al in 1990 (and further explained by Rathfelder and Abriola 1994). This result is also common knowledge in the field and stated in the abstract of the 1990 paper which has 1000+ citations.

Line 38. The first objective is stated as "present a simple and practical approach to solve RE that is efficient, mass conservative," - please clarify if this is an existing approach or a new approach (and if so how it differs from existing approaches). See comments above on the use of existing off-the-shelf solvers for Richards eqn.

Abstract Line 18-19 "Our solution has comparable run-times to Hydrus and better water balance performance (though both models have excellent water balance closure for all the problems we considered)."

Please clarify what is meant by "our solution". Do you mean an off-the-shelf (i.e. pre-existing) ODE algorithm applied to Richard eqn? If so please detail which off-the-shelf algorithm is being used.

Please also clarify the term "better mass balance performance". Given that one of the contributions of the study is a new approach for estimating mass balance, does "mass balance performance" refer to actual mass balance errors, or to how are they estimated within the algorithm? Also please clarify if the same approach for estimating mass balance errors is used in the benchmarking of the algorithms - as that seems pertinent for benchmarking.

Line 97 "though it has been shown that mass conservative solutions are possible (Rathfelder and Abriola, 1994, Tocci et al., 1997). This issue will be explored in detail in this paper."

Given that this is one of the aims, please make it clearer across the paper, as well as in the Abstract/Conclusions, how this study extends what is already known on this topic from earlier work.

Section 2.3 Mass balance closure

This section presents one of the interesting/novel contributions of the study, which is the embedding of the mass balance calculation into the off-the-shelf ODE algorithm.

However this contribution is not so well delineated within this section. For clarity it would make sense to have this contribution in its own (sub)section, rather than presented alongside existing approaches. Then it would be clearer to compare the new approaches to the existing approaches.

In many ways this is the centerpiece of the work (and its contribution as mentioned in the abstract and objectives) - so worth making it stand out.

It is also not clear from the current presentation (at least the Abstract/Conclusions) what are the practical improvements in mass balance accuracy. Can you quantify these improvements in the key sections of the manuscript? This would help re-orient the presentation towards the new results.

2. The algorithms included in the benchmarking are representative of current practices (as in "practical work" rather than "research work"), but do not include algorithms found to be the most efficient in many of the studies cited in the manuscript.

The main method proposed in this work is an off-the-shelf DE solver applied to the pressure form of Richards equation. This solver uses truncation error control through adaptive time stepping.

The benchmarking is done against fixed-step implementations of methods dating back to 1990s, namely the Standard Picard (likely/possibly implemented inconsistently with the cited original paper - see later comments) and Modified Picard (which is arguably the de facto standard in practical modelling tools).

Note that the manuscript presents Hydrus 1D as a distinct method, however (from memory) this software also implements the Modified Picard iteration. Indeed, from the Hydrus 1D description in https://www.ars.usda.gov/pacific-west-area/riverside-ca/agricultural-water-efficiency-and-salinity-research-unit/docs/model/hydrus-1d-model/, "The water content term is evaluated using the mass conservative method [NB: Modified Picard] proposed by Celia et al. [1990]."

So it is hardly surprising that the modern off-the-shelf method is more efficient.

The same comparisons of adaptive methods have been carried out in 1997 by Tocci et al (and their subsequent work using high-order DASPK solvers), by Kavetski et al 2001 (and subsequent work using first/second order implicit methods), by Baron et al (2017), and others. In all these publications adaptive methods were found to outperform fixed-step implementations (and heuristic time stepping implementations). And at least some of these publications have already compared multiple adaptive time stepping methods.

So when the manuscript claims the proposed method is efficient, it should be noted that this is in comparison to methods that have already been found to be inefficient decades ago in the research literature. This seems a rather important point to make to avoid confusing or potentially misleading the readers.

3. Presentation of background theory / methods

The classification of methods on page 4 (lines 105-132) is quite confused and there are many substantial inaccuracies in the descriptions.

Adaptive time stepping (ATS) methods are not a different class of methods to "iterative methods that take one step" or to "one-step methods without iteration".

The name "methods that take one step" is non-standard and quite confusing because all the methods here take multiple steps (indeed thousands to millions of them).

I appreciate that classifying the wide range of numerical methods into neat self-contained categories is not always easy, but there are plenty of good textbook-level references that provide more rigorous classifications and standard naming conventions. For example the standard term is "single-step methods" in reference to methods that carry out a time step using information from that time step alone (in contrast to "multi-step" methods that carry information from preceding time steps).

Line 107 makes a very confusing statement that "methods with no iteration" are defined as methods "where a single evaluation of Eq. (9) is performed for each time step". 

Eq 9 is the continuous time ODE - how does one evaluate it once for each time step? If you are referring to the RHS of that equation - that is also not correct because many methods without iterations evaluate the RHS several times within each time step.

Line 111: "Semi-implicit one-step methods may employ the Thomas Algorithm, or similar, to solve the tri-diagonal matrix problem".

This statement is problematic for many reasons, e.g. it applies only to 1D problems with a particular spatial discretisation (which yields tridiagonal matrices). Further, iterative methods for such 1D problems will also typically use the Thomas algorithm within their iterations, otherwise they would be a rather inefficient implementation!

Line 113. "which is problematic due to the non-linear dependence of K and C on psi (Celia et al., 1990, reproduced below)".

This statement seems to confuse explicit methods (which indeed may not be a good practical choice for Richards eqn) and semi-implicit methods (which include some of the most efficient methods reported for Richards equation).

Line 115: "The iterations allow the solution to find the correct values of K and C at

the next time step"

Strange to refer to any of these values as "correct" - they are all approximations.

Line 128: What is a "black-box" solver in this context?

Overall the text/classification on 105-133 needs to be more or less rewritten from scratch. Here I would recommend contrasting fixed-step vs adaptive-step methods, and then contrasting the methods used within the time steps, etc.

Description of methods used in benchmarking

Line 149-151. "typical one-step method (first-order backward Euler implicit solution), and an excellent iterative method in their improved Modified Picard Iteration Method (MPM) solution"

This description is very confusing and stems from the earlier confusion on lines 105-133. Both the backward Euler method ("Standard Picard") and the Modified Picard method used by Celia et al (1990) are single-step *iterative* methods. It is strange - and very confusing! - to refer to one as "one-step" and to the other as "iterative".

Line 367. "the 'No iteration scheme' (Note: referenced earlier to Celia et al 1990) uses the psi form of RE and solves the problem with a single backward implicit step and no iteration (which we achieved using the Thomas Algorithm)."

As noted earlier, all schemes used by Celia et al 1990 are iterative - for the Standard Picard method (backward Euler solution of pressure form in Celia et al study), this should be clear from Celia et al's equations 4-5. So this sentence describing one of the schemes being benchmarked is incorrect in a very obvious way.

Further, as noted earlier the Thomas Algorithm is also used in the Modified Picard iteration for 1D problems - or at least *should be used*, as otherwise there would almost certainly be a major loss of efficiency. See equation 17 in Celia et al where the tridiagonal coupling between the unknowns can be seen.

Line 158 "All solutions make use of the Thomas Algorithm to solve the tridiagonal linear system arising from the implicit method."

Tridiagonal systems arise only in 1D (in space) problems with certain simple spatial discretisations (e.g. second order central finite differences). To avoid confusion, earlier text should make it clearer you are restricting your attention to 1D problems and the specific spatial discretisation.

Specifically, Line 59 and the rest of the presentation would be clearer if you stated "*1D* vertical flow", and then somewhere around Line 100 would be important to state that you are using second-order differencing in space, which yields tridiagonal matrices.

Line 394. "This is perhaps an under-appreciated limitation of Celia’s MPM solution, and solutions to the mixed form of RE generally – which is that mass balance is a necessary but insufficient criterion for model performance assessment, and truncation errors can still be present in the fluxes with perfect water balance closure"

Please clarify briefly why you consider this to be under-appreciated?

Section 2.4 "Improving efficiency" / subsection 2.4.1 "Defining the Jacobian pattern"

I think this section would make good sense when working with 2D/3D Richards equation where the (Jacobian) matrices are large and sparse. However as far as I understand (and note above) the authors are working with 1D Richards equation and second order spatial finite differences. If so, then the matrices appearing in the discretised (Picard iteration) equations are all tridiagonal and the Jacobian matrix is also tridiagonal.

For tridiagonal matrices, the Thomas algorithm is standard and will be more efficient than alternatives. So I don't get the point of experimenting with different Jacobian patterns given that the actual Jacobian pattern is known exactly and the best linear solver algorithm for it is also well known.

Other

The comparison of compiled languages (listed here as Fortran/C) vs interpreted languages (Python) in Section 2.4.2 is rather subjective and should be indicated as such. For example, line 332 essentially states that interpreted languages are "easier to learn, with excellent teaching resources widely and freely available".

The first of these statements (ease of learning) is clearly subjective - the authors should qualify this as either being their opinion, or (preferably) provide references to some kind of surveys on this topic. The latter two claims are puzzling because there are certainly plenty of excellent resources for Fortran and C, and free compilers are available for both languages. Indeed, the only non-freely available language that comes to mind is Matlab, which happens to be an interpreted language.

Worth mentioning C++ which is perhaps the most modern and widely used compiled language, and for which once again excellent resources and free compilers are commonplace.

This comment is not to deny the obvious appeal of intepreted languages, but a manuscript that focuses on computational software can surely use more accurate and defensible statements when discussing language choices.

Lines 338-344: A citation to a numba reference would seem appropriate here.