A simple, efficient, mass conservative approach to solving Richards’ Equation (openRE, v1.0)
 ^{1}Global Institute for Water Security, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
 ^{2}School for Environment and Sustainability, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
 ^{3}Department of Computer Science, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
 ^{4}Department of Geography and Planning, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
 ^{5}Department of Engineering, Durham University, Durham, UK
 ^{1}Global Institute for Water Security, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
 ^{2}School for Environment and Sustainability, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
 ^{3}Department of Computer Science, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
 ^{4}Department of Geography and Planning, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
 ^{5}Department of Engineering, Durham University, Durham, UK
Abstract. We show that a simple numerical solution procedure – namely the method of lines combined with an offtheshelf ODE solver – can provide efficient, mass conservative solutions to the pressurehead form of Richards’ equation. We present a novel method to quantify the boundary fluxes that reduces water balance errors without negative impacts on model runtimes. We compare this solution with alternatives, including the classic Modified Picard Iteration method of Celia et al. (1990) and the Hydrus 1D model (Šimůnek et al., 2005, 2016). We reproduce a set of benchmark solutions with all models. We find that Celia’s solution has the best water balance, but it can incur significant truncation errors in the simulated boundary fluxes, depending on the time steps used. Our solution has comparable runtimes to Hydrus and better water balance performance (though both models have excellent water balance closure for all the problems we considered). Our solution can be implemented in an interpreted language, such as MATLAB or Python, making use of offtheshelf ODE solvers. We investigated alternative scipy ODE solvers in Python and make practical recommendations about the best way to implement them for Richards’ equation. There are two advantages of our approach: i) the code is short and simple, making it ideal for teaching purposes; and ii) the method can be easily extended to represent alternative properties (e.g., novel ways to parameterize the K(ψ) relationship) and processes (e.g., it is straightforward to couple heat or solute transport), making it ideal for testing alternative hypotheses.
Andrew M. Ireson et al.
Status: final response (author comments only)

RC1: 'Comment on gmd2022185', James Craig, 12 Sep 2022
The paper entitled "A simple, efficient, mass conservative approach to solving Richards’ Equation (openRE, v1.0)" outlines a straightforward implementation to solve the onedimensional Richards' equation using offtheshelf 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 wellwritten 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 supernew 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 nonconstant 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 1D 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 welldocumented 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.
 AC1: 'Reply on RC1', Andrew Ireson, 15 Nov 2022

RC2: 'Comment on gmd2022185', Anonymous Referee #2, 17 Oct 2022
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 semianalytical) 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 offtheshelf 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 doublecheck/correct your implementation of the pressureform 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. ContributionThe 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 offtheshelf ODE solver – can provide efficient, mass conservative solutions to the pressurehead 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 offtheshelf 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: pressureform] 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 followup 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 offtheshelf solvers for Richards eqn.
Abstract Line 1819 "Our solution has comparable runtimes 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 offtheshelf (i.e. preexisting) ODE algorithm applied to Richard eqn? If so please detail which offtheshelf 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 closureThis section presents one of the interesting/novel contributions of the study, which is the embedding of the mass balance calculation into the offtheshelf 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 reorient 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 offtheshelf 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 fixedstep 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/pacificwestarea/riversideca/agriculturalwaterefficiencyandsalinityresearchunit/docs/model/hydrus1dmodel/, "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 offtheshelf 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 highorder 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 fixedstep 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 105132) 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 "onestep methods without iteration".
The name "methods that take one step" is nonstandard 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 selfcontained categories is not always easy, but there are plenty of good textbooklevel references that provide more rigorous classifications and standard naming conventions. For example the standard term is "singlestep methods" in reference to methods that carry out a time step using information from that time step alone (in contrast to "multistep" 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: "Semiimplicit onestep methods may employ the Thomas Algorithm, or similar, to solve the tridiagonal 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 nonlinear 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 semiimplicit 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 "blackbox" solver in this context?
Overall the text/classification on 105133 needs to be more or less rewritten from scratch. Here I would recommend contrasting fixedstep vs adaptivestep methods, and then contrasting the methods used within the time steps, etc.
Description of methods used in benchmarking
Line 149151. "typical onestep method (firstorder 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 105133. Both the backward Euler method ("Standard Picard") and the Modified Picard method used by Celia et al (1990) are singlestep *iterative* methods. It is strange  and very confusing!  to refer to one as "onestep" 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 45. 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 secondorder differencing in space, which yields tridiagonal matrices.
Line 394. "This is perhaps an underappreciated 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 underappreciated?
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.
OtherThe 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 nonfreely 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 338344: A citation to a numba reference would seem appropriate here.
 AC2: 'Reply on RC2', Andrew Ireson, 15 Nov 2022
Andrew M. Ireson et al.
Model code and software
openRE Ireson, A. M. https://doi.org/10.5281/zenodo.6939855
Andrew M. Ireson et al.
Viewed
HTML  XML  Total  BibTeX  EndNote  

270  104  15  389  4  4 
 HTML: 270
 PDF: 104
 XML: 15
 Total: 389
 BibTeX: 4
 EndNote: 4
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1