A simple, efficient, mass-conservative approach to solving Richards' equation (openRE, v1.0)

A simple numerical solution procedure – namely the method of lines combined with an off-the-shelf ordinary differential equation (ODE) solver – was shown in previous work to provide efficient, mass-conservative solutions to the pressure-head form of Richards' equation. We implement such a solution in our model openRE. We developed a novel method to quantify the boundary fluxes that reduce water balance errors without negative impacts on model runtimes – the solver flux output method (SFOM). We compare this solution with alternatives, including the classic modified Picard iteration method and the Hydrus 1D model. 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 off-the-shelf ODE solvers. We evaluated alternative SciPy ODE solvers that are available 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 concise, 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.

1 Introduction

Richards' equation (RE) describes the movement of water in variably saturated porous media. Almost any practical application of RE requires a numerical solution; yet RE remains challenging to solve reliably and accurately for a given set of boundary conditions and soil hydraulic properties (Farthing and Ogden, 2017). RE has been extensively reviewed (e.g., Vereecken et al., 2016; Farthing and Ogden, 2017, and references therein). RE has practical limitations in representing the flow processes in real soils containing macropores, especially with modeling infiltration and rapid percolation processes (Beven and Germann, 2013). The strength of RE is its ability to represent matrix drainage and capillary flows, which control evapotranspiration processes, and its ability to be coupled to heat and solute transport models. For this reason, RE remains a common approach to simulate soil moisture in many terrestrial system models (e.g., vadose zone models, ecohydrology models, and land-surface models; Vereecken et al., 2016; Clark et al., 2015, 2021).

Our objective in this paper is to (i) implement a simple and practical approach to solve RE that is efficient, is mass conservative, and uses open-source software (written in Python) that is readily available, including as a teaching tool; and (ii) present an improved mass balance calculation procedure for use with ordinary differential equation (ODE) solvers that apply adaptive time-stepping (ATS) schemes. We investigate how to maximize the efficiency and accuracy of ODE solvers and provide guidance on the subtle challenges that arise in evaluating the boundary fluxes and the water balance. For our purposes, we define the following five criteria for success for an RE solver: (i) the solver should successfully reproduce benchmark solutions for ψ(t,z) or θ(t,z); (ii) the solver should be mass conservative, with errors that are negligible based on the specific application; (iii) truncation errors in the simulated boundary fluxes should be negligible for the given time step; (iv) the solver should be computationally efficient; and (v) all else being equal on criteria (i)–(iv), the simplest code should be preferred. Simple code should be the most human readable/editable code, which means it should be clean, concise, modular, and free from redundancies.

The remainder of this paper is organized as follows. In Sect. 2, we describe a simple yet powerful approach to solve RE numerically using ODE solvers that we implement and test in the Python and MATLAB programming languages. In Sect. 2, we also discuss the complexities of closing the water balance in RE and present a novel solution that can be applied with any ODE solver: openRE. In Sect. 3, we benchmark the performance of our proposed solution against existing numerical models and solutions, including Hydrus 1D, and the solution method proposed by Celia at al. (1990). In Sect. 4, we summarize our recommendations.

2 Solving Richards' equation

RE is derived from the mass continuity equation applied to a control volume of soil, ΔxΔyΔz (L³), and for one-dimensional vertical flow passing through the area ΔxΔy (L²), we have

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial m}{\partial t}}=-{\displaystyle \frac{\partial \left(q\mathit{\rho }\right)}{\partial z}},\end{array}$$

where m (M L⁻³) is the mass of water per control volume of soil, ρ (M L⁻³) is the density of water, q (L T⁻¹) is the flux of water, and t (T) is time and z (L) is depth below some fixed datum. Assuming that density is constant, we may write

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}=-{\displaystyle \frac{\partial q}{\partial z}},\end{array}$$

where θ (L³ L⁻³) is the volumetric water content, defined as the volume of water per control volume of soil. The vertical flux is given by Darcy's law,

$$\begin{array}{}\text{(3)}& q_z=-K\left(\mathit{\psi }\right){\displaystyle \frac{\mathrm{d}h}{\mathrm{d}z}}=-K\left(\mathit{\psi }\right)\left({\displaystyle \frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}z}}-\mathrm{1}\right),\end{array}$$

where h (L) and ψ (L) are the hydraulic head and matric potential head, respectively, and K(ψ) (L T⁻¹) is the hydraulic conductivity. Combining Eqs. (2) and (3) we have

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K\left(\mathit{\psi }\right)\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}-\mathrm{1}\right)\right),\end{array}$$

which is the mixed form of RE (Celia et al., 1990). If we let $C\left(\mathit{\psi }\right)=\mathrm{d}\mathit{\theta }/\mathrm{d}\mathit{\psi }$ (L⁻¹), we can write

$$\begin{array}{}\text{(5)}& C\left(\mathit{\psi }\right){\displaystyle \frac{\partial \mathit{\psi }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K\left(\mathit{\psi }\right)\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}-\mathrm{1}\right)\right),\end{array}$$

which is the ψ-based form of RE (Celia et al., 1990). We can also express this as the θ form of RE, given by

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{K\left(\mathit{\theta }\right)}{C\left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}-K\left(\mathit{\theta }\right)\right),\end{array}$$

where the constitutive relationships C(θ) and K(θ) are expressed as functions of θ (Celia et al., 1990). Changes in storage in Eqs. (4)–(6) are associated with the filling and draining of soil pores. Because $\mathrm{d}\mathit{\theta }/\mathrm{d}\mathit{\psi }=\mathrm{0}$ in saturated, or close to saturated, soils, an elastic storage term is needed to solve RE in these conditions. This term represents the compression of the pore water and the expansion of the pore space as a function of increasing pore water pressure (though the latter is orders of magnitude larger than the former). The ψ form of RE with elastic storage is written

$$\begin{array}{}\text{(7)}& \left(S_{\mathrm{S}}{\displaystyle \frac{\mathit{\theta }\left(\mathit{\psi }\right)}{{\mathit{\theta }}_{\mathrm{s}}}}+{\displaystyle \frac{\partial \mathit{\theta }}{\partial \mathit{\psi }}}\right){\displaystyle \frac{\partial \mathit{\psi }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K\left(\mathit{\psi }\right)\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}-\mathrm{1}\right)\right),\end{array}$$

where S_S (L⁻¹) is the specific storage coefficient and θ_s is the saturated water content, equal to the porosity (Kavetski et al., 2001). If S_S is large, it may have a non-negligible impact on the water balance, which needs to be accounted for (e.g., as in Clark et al., 2021, Eq. 80). S_S is often treated as a numerical smoothing factor for RE when conditions are saturated or close to saturation, and its impact on the water balance is neglected (e.g., as in Tocci et al., 1997; Ireson et al., 2009). For convenience, hereafter we define $C\left(\mathit{\psi }\right)=S_{\mathrm{S}}\frac{\mathit{\theta }\left(\mathit{\psi }\right)}{{\mathit{\theta }}_{\mathrm{s}}}+\frac{\partial \mathit{\theta }}{\partial \mathit{\psi }}$, such that elastic storage can be ignored by setting S_S=0.

The ψ form or mixed form of RE with elastic storage included can be shown to work for saturated conditions (Miller et al., 1998) and so can be considered a general governing equation for flow in soils, aquitards, and confined and unconfined aquifers. Some numerical solutions to the ψ form of RE are reportedly subject to poor mass conservation (Milly, 1984, 1985; Celia et al., 1990; Farthing and Ogden, 2017; Clark et al., 2021; Tubini and Rigon, 2022), though it has been shown that mass balance errors can be effectively controlled using ATS (Rathfelder and Abriola, 1994; Tocci et al., 1997). We carefully assess the mass balance performance of our model in this paper.

We consider a finite-difference numerical solution for the ψ form of RE that applies the method of lines to reduce the partial differential equation (PDE) in Eq. (7) to a system of ODEs of the form

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}t}}=f(\mathit{\psi },t,\mathit{z}),\end{array}$$

where ψ and z represent vectors containing discrete values of ψ and z. A similar approach was taken by Tocci et al. (1997) and Farthing and Ogden (2017). There are various methods that can be applied to integrate Eq. (8) with respect to time – for our purposes, it is useful to consider three specific classes of methods.

• Fixed-step non-iterative methods (NIMs). These methods provide a solution for a fixed time interval and can be either explicit or implicit. Explicit methods are simple but impractical because they require very small time steps. Semi-implicit non-iterative methods take a single iteration based on the linearization of the spatially discretized governing equations. These methods are more stable than explicit methods but have limitations for solving RE over a fixed time interval because they do not iterate to improve the convergence of the problem, which is problematic due to the non-linear dependence of K and C on ψ (Celia et al., 1990). Celia et al. (1990) did not present results using a NIM, but their Eq. (4) can be solved directly (a semi-implicit non-iterative method, which we implemented and discuss in Sect. 3.1). We include NIM here because they may be instructive to new researchers learning about these problems.

• Fixed-step iterative methods (FIMs). Iterative implicit methods are very widely used (Celia et al., 1990; Rathfelder and Abriola, 1994). The iterations allow the solution to find more representative values of K and C over the time step. These methods can be applied to either the ψ form or mixed form of RE. For the mixed form of RE, the convergence criterion can be based directly on water balance closure, as in Celia's modified Picard iteration method. For the ψ form of RE, the convergence criterion is based on ψ values, and the solution is subject to larger water balance errors (Celia et al., 1990; Rathfelder and Abriola, 1994; Farthing and Ogden, 2017).

• Adaptive time-stepping (ATS) methods. An implicit or explicit method is used to find ψ_t+Δt over some calculation time step; the truncation error in the solution is assessed, for example, by comparing two different-order solutions, as in Kavetski et al. (2001); and depending on the size of the error, the time step is either reduced to improve the accuracy or increased to improve the efficiency. The solution marches forward until reaching what we term the reporting time step, where the state variables are output. The advantage of this approach is that the states and fluxes calculated at the intermediate calculation steps contain useful information that can be exploited in the output, as we demonstrate in this paper. Kavetski et al. (2001, 2002a, b) developed ATS solvers designed specifically to solve the different forms of RE, while other workers have applied readily available “black-box” (Kavetski, et al., 2001) ODE or DAE (differential algebraic equation) solvers to RE (Tocci et al., 1997; Ireson et al., 2009; Ireson and Butler, 2013; Mathias et al., 2015; Clark et al., 2021). A wide range of ODE solvers are available in many different programming languages. Solutions are simple to implement and, as we show in this study, can outperform other methods in terms of accuracy and efficiency.

In this study, we implement each of the three possible methods for solving RE in scripts that are available from Ireson (2022, https://github.com/amireson/openRE, last access: 24 January 2023). All of the models in this study were coded in Python (version 3.8.11). We make use of the following libraries: NumPy (version 1.20.3); Matplotlib (version 3.4.2); SciPy (version 1.7.1), which contains various ODE solvers, described below; and Numba (version 0.53.1, Lam et al., 2015), which is a just-in-time (JIT) compiler that is optional but speeds up the model runs considerably. We organize and run the models using makefiles (Jackson, 2016).

We also provide a MATLAB version of our recommended solution (implemented in MATLAB R2017b installed on a Mac). The MATLAB implementation will not be described further, but compared with the optimal Python solution, the MATLAB model is somewhat simpler and achieves an equivalent performance in terms of simulated states and fluxes. We do not recommend Python over MATLAB (or vice versa) – both platforms work well, and the choice will come down to numerous factors, including, e.g., the availability of, and user familiarity with, either package or the need to use non-proprietary software to satisfy open-science requirements.

2.1 Fixed-step non-iterative and iterative solutions

For the NIM and FIM methods, we have coded up the numerical solutions from Celia et al. (1990). This provides a typical NIM method (first-order backward Euler implicit solution, Eq. 4 in Celia et al., 1990) and two alternative FIMs: Picard iteration that solves the ψ form of RE and their improved modified Picard iteration method (MPM) solution that solves the mixed form of RE. Fully reproducible details of these models were provided by Celia et al. (1990) and so will not be repeated here. These models were implemented in Python, and the code is provided at Ireson (2022, https://github.com/amireson/openRE, last access: 24 January 2023). The NIM model was coded up in 67 lines of code (note that blank lines and comment lines are not counted in the number of lines of code), with an additional 52 lines of code to configure the problem (define the grid, hydraulic properties, etc.). The FIM Picard iteration method was coded up in 79 lines of code and required the same 52 lines of code to configure the problem. The FIM modified Picard iteration was implemented using the Numba just-in-time compiler (Sect. A5) and was coded up in 90 lines of code, with an additional 59 lines of code to configure the problem. All solutions make use of the Thomas algorithm to solve the tridiagonal linear system arising from the implicit method.

2.2 Adaptive time-stepping solution

One major benefit of using a standard ODE solver to integrate RE is that the code is concise and easy to read and understand. Our objective is to write a function that will evaluate the right-hand side of Eq. (8) and feed this to the ODE solver. Within this function, the problem can be treated as instantaneous in time so that we only need to consider the spatial differences in our finite-difference solution scheme. We will use a cell-centered grid (Bear and Cheng, 2010, p. 533) in space; i.e., state variables are stored at nodes located at the center of grid cells, while fluxes are defined at the cell boundaries, as shown schematically in Fig. 1 (note that this is equivalent to what is called a “staggered grid” in fluid dynamics). For simplicity here, we consider a uniform grid (constant Δz), but it is straightforward to adapt these solutions to non-uniform grids. We introduce two spatial indices (Fig. 1): i represents the nodal values, and j represents the cell boundaries, both of which have an initial value of zero (because Python uses zero-based indexing). Hence, $j=i+\mathrm{1}/\mathrm{2}$. We will start with the ψ form of RE, for which the governing equation is given in the form

$$\begin{array}{}\text{(9)}& {\left.{\displaystyle \frac{\partial \mathit{\psi }}{\partial t}}\right|}_i=-{\displaystyle \frac{\mathrm{1}}{C\left({\mathit{\psi }}_i\right)}}{\left.{\displaystyle \frac{\partial q}{\partial z}}\right|}_i\end{array}$$

and

$$\begin{array}{}\text{(10)}& {\left.{\displaystyle \frac{\partial q}{\partial z}}\right|}_i={\displaystyle \frac{q_{i+\frac{\mathrm{1}}{\mathrm{2}}}-q_{i-\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{\Delta }z}}={\displaystyle \frac{q_{j+\mathrm{1}}-q_j}{\mathrm{\Delta }z}}.\end{array}$$

Figure 1 Schematic representation of the cell-centered finite-difference grid, for a soil column of depth L, showing the zero-based Python indices of the states (ψ_i) and fluxes (q_j) on the left and depths of the nodes on the right, assuming a regular grid. There are N states and N+1 fluxes. The upper- and lower-boundary conditions are q_j=0 and $q_{j_N}$, respectively.

Here, the fluxes are given by

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{q}_{i+\frac{\mathrm{1}}{\mathrm{2}}}=& q_{j+\mathrm{1}}=-{\displaystyle \frac{K\left({\mathit{\psi }}_{i+\mathrm{1}}\right)+K\left({\mathit{\psi }}_i\right)}{\mathrm{2}}}\\ & \times \left({\displaystyle \frac{{\mathit{\psi }}_{i+\mathrm{1}}-{\mathit{\psi }}_i}{\mathrm{\Delta }z}}-\mathrm{1}\right)\\ q_{i-\frac{\mathrm{1}}{\mathrm{2}}}=& q_j=-{\displaystyle \frac{K\left({\mathit{\psi }}_i\right)+K\left({\mathit{\psi }}_{i-\mathrm{1}}\right)}{\mathrm{2}}}\\ & \times \left({\displaystyle \frac{{\mathit{\psi }}_i-{\mathit{\psi }}_{i-\mathrm{1}}}{\mathrm{\Delta }z}}-\mathrm{1}\right).\end{array}\end{array}$$

In Eq. (11), we are using the arithmetic mean of K at the nodal points to estimate K at the cell boundaries, but other options are possible (see, e.g., Bear and Cheng, 2010, p. 535). It is possible to combine Eqs. (9)–(11), but keeping them separate keeps the code modular and simple.

There are three commonly used boundary conditions for this problem, namely

• i.

  a specified flux, q_T (L T⁻¹) (type-II boundary) is often used at the upper boundary to represent infiltration, where

  $$\begin{array}{}\text{(12)}& q_{j=\mathrm{0}}=q_{\mathrm{T}};\end{array}$$
• ii.

  a free-drainage boundary is often used at the lower boundary, where

  $$\begin{array}{}\text{(13)}& q_{j=N}=K\left({\mathit{\psi }}_{i=N-\mathrm{1}}\right);\end{array}$$
• iii.

  a fixed ψ (type-I) boundary maybe used at the upper boundary, ψ_T (L), typically to indicate a ponding depth, or at the lower boundary, ψ_B (L), typically to represent a fixed water table, where

  $$\begin{array}{}\text{(14)}& \begin{array}{rl}{q}_{j=\mathrm{0}}=& -{\displaystyle \frac{K\left({\mathit{\psi }}_{\mathrm{T}}\right)+K\left({\mathit{\psi }}_{i=\mathrm{0}}\right)}{\mathrm{2}}}\\ & \times \left({\displaystyle \frac{{\mathit{\psi }}_{\mathrm{T}}-{\mathit{\psi }}_{i=\mathrm{0}}}{\mathrm{\Delta }z/\mathrm{2}}}-\mathrm{1}\right)\\ q_{j=N}=& -{\displaystyle \frac{K\left({\mathit{\psi }}_{i=N-\mathrm{1}}\right)+K\left({\mathit{\psi }}_{\mathrm{B}}\right)}{\mathrm{2}}}\\ & \times \left({\displaystyle \frac{{\mathit{\psi }}_{i=N-\mathrm{1}}-{\mathit{\psi }}_{\mathrm{B}}}{\mathrm{\Delta }z/\mathrm{2}}}-\mathrm{1}\right).\end{array}\end{array}$$

Box 1 provides Python-based pseudo-code that implements this solution (Eqs. 9–14) in a function for $\mathrm{d}\mathit{\psi }/\mathrm{d}t$ with a type-II boundary at the upper boundary and a free-drainage boundary at the lower boundary and contains just seven lines of code. This function can be called by the ODE solver.

Box 1 Python-based pseudo-code implementation of RE with a cell-centered finite-difference approach in a function to be called by an ODE solver. Arrays are zero-indexed; q[0] and q[−1] refer to the first and last item in the q array, respectively, and K[1:] and $K[:-\mathrm{1}]$ refer to a slice of the array K from the second to the last node and from the first to the second-to-last node, respectively.

2.3 Mass balance closure

When RE is solved over some time interval, $t=t_{\mathrm{0}}\to t_M$ (where t_M−t₀ typically corresponds to multiple years in practical application) for a soil profile $\mathrm{0}\le z\le L$, the cumulative inflow minus outflow should equal the change in storage in the profile over the same interval; i.e., ϵ_B (mm), the bias error, defined in Eq. (15), should be zero.

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\mathit{ϵ}}_{\mathrm{B}}=& \underset{t=t_{\mathrm{0}}}{\overset{t_M}{\int }}\left(q\right(t,\mathrm{0})-q(t,z_N\left)\right)\mathrm{d}t\\ & -\underset{z=\mathrm{0}}{\overset{L}{\int }}\left(\mathit{\theta }\right(t_M,z)-\mathit{\theta }(t_{\mathrm{0}},z\left)\right)\mathrm{d}z\end{array}\end{array}$$

The bias error can be treated as a mass balance performance metric for the model, but this metric may underestimate the true errors in the water balance that occur within the time period simulated, which may cancel out over the entire run. The metric used by Celia et al. (1990), Rathfelder and Abriola (1994), and Tocci et al. (1997) has the same problem. A more rigorous mass balance performance metric is the root mean squared error of the daily (or some other time increment) cumulative net flux minus the change in storage, ϵ_R (mm), where

$$\begin{array}{}\text{(16)}& {\mathit{ϵ}}_{\mathrm{R}}=\sqrt{\sum _{j=\mathrm{1}}^M{\left(\begin{array}{l}\underset{t=t_j}{\overset{t_{j+\mathrm{1}}}{\int }}\left(q\right(t,\mathrm{0})-q(t,z_N\left)\right)\mathrm{d}t\\ -\underset{z=\mathrm{0}}{\overset{L}{\int }}\left(\mathit{\theta }\right(t_{j+\mathrm{1}},z)-\mathit{\theta }(t_j,z\left)\right)\mathrm{d}z\end{array}\right)}^{\mathrm{2}}/M},\end{array}$$

where j is an index in time, and M is the number of time steps considered. Reporting both metrics is informative – a high ϵ_R with a low ϵ_B indicates that daily errors are occurring but canceling one another out; a high ϵ_B with a low ϵ_R indicates that small daily errors are systematic in one direction and hence accumulate to give a high bias.

The fluxes in the mass balance calculation depend on the boundary conditions. For type-I (specified ψ) and free-drainage-type boundary conditions, the boundary flux depends on the simulated ψ value at the node closest to the boundary. Over a time increment t_j to t_j+1, ψ values will change continuously and hence so will the boundary flux. Due to the non-linearity of the K(ψ) relationship, the flux will change in a non-linear manner. The cumulative flux over the time increment, $Q_{j\to j+\mathrm{1}}$ (mm), is given by

$$\begin{array}{}\text{(17)}& Q_{j\to j+\mathrm{1}}=\underset{t=t_j}{\overset{t_{j+\mathrm{1}}}{\int }}q\left(t\right)\mathrm{d}t.\end{array}$$

$Q_{j\to j+\mathrm{1}}$ is estimated from discrete values of q (which in turn are approximated from discrete values of ψ, using either a forward difference approximation, where $Q_{j\to j+\mathrm{1}}\approx q_j\mathrm{\Delta }t$; a backward difference approximation (as in Celia et al., 1990), where $Q_{j\to j+\mathrm{1}}\approx q_{j+\mathrm{1}}\mathrm{\Delta }t$; or a central difference approximation (as in trapezoidal integration), where $Q_{j\to j+\mathrm{1}}\approx (q_j+q_{j+\mathrm{1}})\mathrm{\Delta }t/\mathrm{2}$. These discrete approximations for $Q_{j\to j+\mathrm{1}}$ can be poor if the time step is large – or, more precisely, if the changes in q over a time step are large and non-linear.

This leads to an important limitation with Celia's mixed form solution to RE (and other equivalent iterative solutions). This solution has excellent mass balance closure, but because it uses a fixed-step iterative solution procedure with a backward difference approximation for $Q_{j\to j+\mathrm{1}}$, the simulated boundary fluxes can be shown to be sensitive to the time step (see Sect. 3.2). Hence, even though for larger Δt the water balance is still perfectly closed, the actual terms within the water balance have changed, so there is less inflow and less change in storage.

ATS solvers can provide a practical solution to solving RE with good mass balance performance and boundary fluxes that do not depend on the user-specified time step. The basic idea behind ATS solvers is that when there are large changes in ψ in the model, small steps can be taken to capture the shape of C(ψ) and $Q_{j\to j+\mathrm{1}}$ and minimize the integration errors over a time step. When the changes in ψ are small, larger steps can be taken to maximize efficiency. We will refer to these adaptive time steps as calculation steps. The user specifies the time steps at which they wish the results to be saved, which we will refer to as the reporting time steps. In typical practical applications, the reporting step would be the resolution of the driving data (e.g., hourly or daily). The solver may take many calculation steps of varying lengths between the reporting steps, saving the outcomes internally each time that the error tolerance is satisfied and a successful step is taken.

To accurately calculate the boundary fluxes, it is necessary to use the ψ information from the calculation time steps because ψ may have evolved non-linearly over the reporting time step. However, this is not trivial. Two possible ways to do this (which are equivalent to one another) are to (i) enable dense output from the ODE solver (if this feature is supported by the ODE solver) or (ii) force the ODE solver to make the reporting steps equal to the calculation steps. However, both of these approaches are more computationally demanding in terms of memory and runtime – a significant disadvantage. We propose here a third method for calculating the boundary fluxes that can be used with any ODE solver. At an instance in time, the cumulative boundary flux, Q, is related to the instantaneous flux, q, by

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}=q.\end{array}$$

We can therefore use the ODE solver to integrate this expression and solve for Q. To do this, we add to the system of ODEs defined by Eq. (8) two new ODE expressions that represent the cumulative boundary fluxes. The dependent variable vector that is sent to the ODE solver now is F, defined

$$\begin{array}{}\text{(19)}& \mathit{F}=\left[\begin{array}{c}{Q}_{T:t_{\mathrm{0}}\to t}\\ {\mathit{\psi }}_{\mathrm{0}}\\ {\mathit{\psi }}_{\mathrm{1}}\\ \mathrm{⋮}\\ {\mathit{\psi }}_{N-\mathrm{2}}\\ {\mathit{\psi }}_{N-\mathrm{1}}\\ Q_{B:t_{\mathrm{0}}\to t}\end{array}\right],\end{array}$$

where $Q_{T:t_{\mathrm{0}}\to t}$ and $Q_{B:t_{\mathrm{0}}\to t}$ are the cumulative boundary fluxes at the current time, t, since the start of the simulation, t₀. The ODE solver will integrate the equation

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\mathrm{d}\mathit{F}}{\mathrm{d}t}}=f(\mathit{F},t,\mathit{z})\end{array}$$

to solve for F. The function that is called by the ODE solver will evaluate the vector

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\mathrm{d}\mathit{F}}{\mathrm{d}t}}=\left[\begin{array}{c}\frac{\mathrm{d}{Q}_{T:t_{\mathrm{0}}\to t_j}}{\mathrm{d}t}\\ \frac{\mathrm{d}{\mathit{\psi }}_{\mathrm{0}}}{\mathrm{d}t}\\ \frac{\mathrm{d}{\mathit{\psi }}_{\mathrm{1}}}{\mathrm{d}t}\\ \mathrm{⋮}\\ \frac{\mathrm{d}{\mathit{\psi }}_{N-\mathrm{2}}}{\mathrm{d}t}\\ \frac{\mathrm{d}{\mathit{\psi }}_{N-\mathrm{1}}}{\mathrm{d}t}\\ \frac{\mathrm{d}{Q}_{B:t_{\mathrm{0}}\to t_j}}{\mathrm{d}t}\end{array}\right]=\left[\begin{array}{c}{q}_{\mathrm{T}}\\ -\frac{\mathrm{1}}{C\left({\mathit{\psi }}_{\mathrm{0}}\right)}{\left.\frac{\partial q}{\partial z}\right|}_{\mathrm{0}}\\ -\frac{\mathrm{1}}{C\left({\mathit{\psi }}_{\mathrm{1}}\right)}{\left.\frac{\partial q}{\partial z}\right|}_{\mathrm{1}}\\ \mathrm{⋮}\\ -\frac{\mathrm{1}}{C\left({\mathit{\psi }}_{N-\mathrm{2}}\right)}{\left.\frac{\partial q}{\partial z}\right|}_{N-\mathrm{2}}\\ -\frac{\mathrm{1}}{C\left({\mathit{\psi }}_{N-\mathrm{1}}\right)}{\left.\frac{\partial q}{\partial z}\right|}_{N-\mathrm{1}}\\ q_{\mathrm{B}}\end{array}\right].\end{array}$$

We note that each term in Eq. (21) is expressed at a single instant in time, t, and subscripts $\mathrm{0},\mathrm{1},\mathrm{\dots },N-\mathrm{2},N-\mathrm{1}$ refer to the indices of the finite-difference discretization points, and we use zero-based indexing to be consistent with the Python language and Fig. 1. After solving for F, the first and last rows of F correspond to the cumulative boundary fluxes, which at time t are $Q_{T:t_{\mathrm{0}}\to t}$ and $Q_{B:t_{\mathrm{0}}\to t}$. The cumulative fluxes over each time step are obtained from

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& Q_{T:t_j\to t_{j+\mathrm{1}}}=Q_{T:t_{\mathrm{0}}\to t_{j+\mathrm{1}}}-Q_{T:t_{\mathrm{0}}\to t_j}\\ & Q_{B:t_j\to t_{j+\mathrm{1}}}=Q_{B:t_{\mathrm{0}}\to t_{j+\mathrm{1}}}-Q_{B:t_{\mathrm{0}}\to t_j}.\end{array}\end{array}$$

Solving systems of ODEs in this way is straightforward – requiring the user to pack and unpack the dependent variable vector. Python-based pseudo-code showing how this can be implemented is given in Box 2. Note that here we multiply the fluxes by 1000 within the solver (equivalent to converting the fluxes from m d⁻¹ to mm d⁻¹), so that the magnitude of the fluxes is comparable with the magnitude of the changes in ψ (typically in m), which can improve the water balance estimate. Any arbitrary scaling factor can be applied here, as long as the output fluxes are re-scaled to the correct units. We shall refer to this method for calculating the boundary fluxes as the solver flux output method (SFOM).

Box 2 Python-based pseudo-code implementation of an ODE solver solution. The ODE function returns the time derivative of the dependent variable array, given in Eq. (21). Note that $F[:,\mathrm{1}:-\mathrm{1}]$ is a slice through the array F that takes all the items in the first (time) dimension and the second to second-to-last items in the second dimension.

2.4 Improving efficiency

Here, we provide details of two techniques that can be used to improve the computational runtime of the model. These methods have no impact on the accuracy of the solution, so they are optional, but combined they can result in better than a factor of 10 reduction in the runtime at the cost of only a few additional lines of code. In Appendix A, we investigate the impact of a range of different possible model decisions/assumptions on the accuracy, mass balance, efficiency, and simplicity of the model.

2.4.1 Defining the Jacobian pattern

For the form of RE given by Eq. (9), the Jacobian matrix, J, is an n×n matrix, where the cell in each row, i, and column, j, is defined by the derivative

$$\begin{array}{}\text{(23)}& J_{i,j}={\displaystyle \frac{\mathrm{d}}{\mathrm{d}{\mathit{\psi }}_j}}\left({\left.{\displaystyle \frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}t}}\right|}_i\right)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}{\mathit{\psi }}_j}}\left(-{\displaystyle \frac{\mathrm{1}}{C\left({\mathit{\psi }}_i\right)}}{\left.{\displaystyle \frac{\partial q}{\partial z}}\right|}_i\right).\end{array}$$

Each entry in J_i,j can be evaluated in a function that is passed to the ODE solver (assuming that the particular ODE solver being used has this functionality), in order to speed up the solution process. For the spatial discretization scheme described above, all the terms of the Jacobian are zero, except for where $j=i-\mathrm{1}$, j=i, and $j=i+\mathrm{1}$ (ignoring when $i-\mathrm{1}<\mathrm{1}$ and $i+\mathrm{1}>n$, which would be zero terms anyway, for any boundary condition). A simpler alternative to defining the full Jacobian matrix is to define the Jacobian pattern. The Jacobian pattern is a matrix of ones and zeros that defines the location of the structurally non-zero elements of the Jacobian – that is, where the terms are not identically zero. For the spatially discretized RE as given in Eq. (9), the Jacobian pattern is a simple tridiagonal matrix, with ones on the three main diagonals and zeros everywhere else. To implement this requires an ODE solver capable of using the Jacobian pattern (also referred to as the Jacobian sparsity matrix). The SciPy ODE solvers ode and solve_ivp have the ability to define a banded Jacobian pattern: setting uband and lband arguments to 1 tells the ODE solver that the Jacobian is a tridiagonal matrix. The SciPy ODE solver solve_ivp can also handle a general n×n Jacobian pattern, which is more adaptable for multi-dependent variable coupled problems (e.g., Goudarzi et al., 2016). The MATLAB ODE solvers can read the Jacobian sparsity pattern matrix from the JPattern argument. We report on the relative performance/complexity of each of these methods in Sect. A4.

2.4.2 Just-in-time compilation

In this paper, we are providing guidance for the use of interpreted programming languages (e.g., Python or MATLAB) to solve RE. Interpreted languages have a number of advantages over compiled languages (such as FORTRAN, C, and C$++)$, including that they are, at least in our opinion, easier to learn, with excellent teaching resources widely and freely available; they tend to have higher level abstractions, so that the same task can be completed in fewer lines of code; and they are cross platform and typically easier to install. Interpreted languages are not pre-compiled and are hence slower to execute than compiled languages. A nice compromise between the simplicity of interpreted languages and efficiency of compiled languages is to use a just-in-time compiler. In Python, the Numba library is such a just-in-time compiler (Lam et al., 2015). Numba compiles selected Python functions once at the start of the runtime, and then all subsequent calls to the code run much faster. We find that using Numba in conjunction with our preferred ODE solver solution described above, results in up to 10× faster code execution (see Sect. A5). The drawback to using Numba is that some re-factoring of the code may be necessary to make a script that previously ran without Numba work using Numba – in particular, there are complications around how variables are allocated into NumPy arrays. We include code in Ireson (2022, https://github.com/amireson/openRE, last access: 24 January 2023) that demonstrates how to successfully implement Numba.

3 Benchmarking the model performance

In this section, we run our RE solver, openRE, on a number of benchmark problems, comparing the different solution procedures and assessing the performance of all solutions against the five success criteria identified in the introduction, namely (i) accuracy of θ(t,z) and ψ(t,z), (ii) mass balance performance, (iii) consistent boundary fluxes with Δt, (iv) computational efficiency (i.e., runtime), and (v) simplicity of the code. For the purposes of comparing efficiency (iv), all simulations were run on the same laptop computer, and we report the runtimes as a measure of relative performance. For the purposes of comparing simplicity of the code (v) we use a very simple metric of lines of code, which is reported above in Sect. 2.

3.1 Published model benchmarks

3.1.1 Celia's problem

Celia's test case (Celia et al., 1990) is used to compare our ATS solution with the different solutions schemes previously proposed by Celia et al. (1990). The test problem uses a 40 cm deep vertical soil profile (z_N=40) with a uniform 1 cm space step (dz=1), a 360 s duration ($t_{\mathrm{0}}=\mathrm{0};t_M=\mathrm{360}$) with a 1 s time step (Δt=1), and the following initial and type-I boundary conditions:

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& \mathit{\psi }(t=t_{\mathrm{0}},\mathrm{0}\le z\le z_N)=-\mathrm{61.5}\mathrm{cm}\\ & \mathit{\psi }(t_{\mathrm{0}}\le t\le t_M,z=\mathrm{0})=-\mathrm{20.7}\mathrm{cm}\\ & \mathit{\psi }(t_{\mathrm{0}}\le t\le t_M,z=z_N)=-\mathrm{61.5}\mathrm{cm}.\end{array}\end{array}$$

Figure 2 Reproduction of Celia's model benchmark: the ψ breakthrough curve using a backward difference implicit Euler solution with the no-iteration scheme, backward difference implicit Euler with Picard iteration, Celia's modified Picard iteration method (MPM), and the ODE solver with adaptive time stepping. The ODE solver produces consistent ψ(tz) results independently of the mass balance calculation procedure. The time steps reported are calculation time steps for Celia's solutions but are reporting time steps for the ODE solver, which uses an adaptive time step for calculation steps.

The soil hydraulic properties are given by

$$\begin{array}{}\text{(25)}& \begin{array}{rl}& \mathit{\theta }={\displaystyle \frac{\mathit{\alpha }({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})}{\mathit{\alpha }+{\left|\mathit{\psi }\right|}^{\mathit{\beta }}}}+{\mathit{\theta }}_{\mathrm{r}}\\ & K=K_{\mathrm{s}}{\displaystyle \frac{A}{A+{\left|\mathit{\psi }\right|}^{\mathit{\gamma }}}},\end{array}\end{array}$$

where the parameter values are α = 1.611 × 10⁶, θ_s = 0.287, θ_r = 0.075, β = =3.96, K_s = 00944 cm s⁻¹, A = 1.175 × 10⁶, and γ=4.74. Celia et al. (1990) presented three solution schemes: the “no-iteration scheme” uses the ψ form of RE and solves the problem with a single backward implicit step and no iteration (which we achieved using the Thomas algorithm), the “Picard iteration scheme” also solves the ψ form of RE but uses the Picard iteration method to improve the solution, with errors in ψ used as a convergence criterion, and finally the “modified Picard iteration method” (MPM) uses the mixed-form of RE and uses errors in θ as a convergence criterion. The MPM is mass conservative because the iteration ensures that the cumulative change in fluxes (the right-hand side of RE) balances the changes in storage (the left-hand side of the mixed form of RE). Celia's three solutions were reproduced in Python scripts (https://github.com/amireson/openRE, last access: 24 January 2023) and compared with our ATS/SFOM solution. Results from all three solutions are shown in Fig. 2. All methods are consistent for very small time steps. The fixed-step method with no iteration performs poorly, with delayed breakthrough of the wetting front when the time step is large. The solution is improved by using the Picard iteration, but there are still some delays. The MPM has a much better performance but, as with all implicit Euler time-stepping schemes, is still subject to some numerical dispersion for larger time steps (van Genuchten and Gray, 1978). The ATS solution reproduces the ψ breakthrough curve but with no dispersion and no differences associated with the time step. We note that the time steps for plotting the ATS solution only represent the reporting time step – the underlying calculation time steps are likely much smaller (Sect. 2.2).

Figure 3 Mass balance of Celia's three models and our adaptive time-stepping model. Panel (a) shows a reproduction of Celia's mass balance calculation. Panel (b) shows the cumulative inflow that is simulated.

In Fig. 3, we show the cumulative inflow simulated by each of these models for Celia's benchmark problem, along with the mass balance bias error, for different reporting time steps. The fixed-step solution with no iteration and the Picard iteration solution both have poor mass balance performance unless the time step is very small – on this basis, we do not consider these solutions further. We see that the MPM method is perfectly mass conservative for any Δt used in the model, as we should expect. However, we can see that the cumulative inflow is sensitive to Δt. Hence, even though for larger Δt the water balance is still perfectly closed, the actual terms within the water balance have changed, so there is less inflow and less change in storage. 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 even with perfect water balance closure.

In Fig. 3, we also show the water balance performance of our ATS solution, using either reporting-time-step information for the water balance calculation or using calculation-time-step information (i.e., using the SFOM described in Sect. 2.3). Using reporting-time-step information is the easiest and most intuitive approach to take – you numerically integrate (1) discrete θ values over depth to get storage and (2) discrete q values during reporting time steps to get Q (e.g., using trapezoidal integration). However, this approach fails to capture non-linear changes in q over a reporting time step and results in large water balance errors and errors in the cumulative fluxes, as is clear in Fig. 3. Using the SFOM, we see that the water balance is almost exactly closed, and the boundary fluxes are independent of the reporting time step. It is also important to note that the discrete ψ(tz) values simulated by both ATS solution procedures here are identical (see Fig. 2) – the only difference is how the boundary fluxes are calculated.

3.1.2 Miller's saturated infiltration pulse problem

Miller et al. (1998) investigated solutions to RE that aimed to address numerical convergence problems associated with challenging boundary conditions, and they specifically looked at the problem of infiltration from a ponded upper boundary into a hydrostatic soil profile with a fixed water table at the lower boundary. This is a good benchmark because it requires the model to deal with perched saturated conditions over unsaturated conditions and involves highly non-linear changes in properties over short distances and time steps. The problem uses van Genuchten (1980) hydraulic properties, given by

$$\begin{array}{}\text{(26)}& {\displaystyle }{S}_{\mathrm{e}}=(\mathrm{1}+(\mathit{\alpha }\mathit{\psi }{)}^n{)}^{-m},\text{(27)}& {\displaystyle }\mathit{\theta }={\mathit{\theta }}_{\mathrm{r}}+({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})S_{\mathrm{e}},\text{(28)}& {\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}\mathit{\psi }}}={\displaystyle \frac{-\mathit{\alpha }m({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})}{\mathrm{1}-m}}{S}_{\mathrm{e}}^{\mathrm{1}/m}(\mathrm{1}-S_{\mathrm{e}}^{\mathrm{1}/m}{)}^m,\text{(29)}& {\displaystyle }K=K_{\mathrm{s}}{S}_{\mathrm{e}}^{\mathrm{1}/\mathrm{2}}(\mathrm{1}-(\mathrm{1}-S_{\mathrm{e}}^{\mathrm{1}/m}{)}^m),\end{array}$$

where S_e (–) is the effective saturation; α (L⁻¹), n (–), and m (–) are parameters that determine the shape of the θ(ψ) curve; θ_r (–) and θ_s (–) are the residual and saturated volumetric water contents; and K_s (L T⁻¹) is the saturated hydraulic conductivity. Miller's problem uses the parameters in Table 1 for three different soil types.

Table 1 Soil hydraulic properties used the Miller et al. (1998) problem.

Table 2 Runtime and water balance performance for the Celia MPM and ATS SFOM solutions, applied to the problem in Miller et al. (1998). We also show here the solver settings. For the ATS solutions, we had to reduce the relative and absolute tolerances (atol/rtol) of the integrator from the default settings and increase the maximum number of calculations steps allowed within a reporting step (nsteps). Note that time steps (dt) in ATS solutions are reporting time steps.

A hydrostatic initial condition is used, with a fixed water table at depth of 10, 5, and 2 m below ground surface for sand, loam, and clay loam, respectively. At the upper boundary, 0.1 m of ponding is applied throughout the simulation runtime of 0.18, 2.25, and 1.0 d for sand, loam, and clay loam, respectively. We simulated this problem with our ATS solution and with Celia's MPM model for comparison purposes. Both models faced challenges with this problem. For Celia's MPM, we had to use a small time step to get the solver to produce accurate ψ(z) profiles (Fig. 4). For the ATS solutions using the default ODE solver settings, the models failed to propagate the wetting front into the soil correctly. It was necessary to increase the maximum number of calculation steps allowed per reporting time step (we increased this from the default 500 to 10 000) so that very small time steps could be taken (Table 2). It was also necessary (for loam and clay loam) or beneficial (for sand) to reduce the solver error tolerances – see values in Fig. 4 and Table 2. The results from these simulations are shown in Fig. 4 and are consistent with those reported in Fig. 1 of Miller et al. (1998), showing that both models are able to successfully reproduce this benchmark. The runtimes and water balance for each solution are tabulated in Table 2. Celia's MPM has consistently better water balance performance, though we think the water balance errors in both models are acceptably low. The ATS solution is slower for sand, faster for loam, and about the same for clay loam. We note that the runtimes and water balances of the ATS solution are sensitive to the reporting time step Δt and the solver settings nsteps, atol, and rtol – an improved solution might be attainable by optimizing these settings. On the other hand, for Celia's MPM solution, we only needed to optimize Δt.

Figure 4 Reproduction of the Miller infiltration pulse result, using Celia's MPM model and our ATS SFOM solution. Both models are satisfactorily consistent with the output reported by Miller et al. (1998), in their Fig. 1.

3.1.3 Mathias' solution for horizontal infiltration

Mathias and Sander (2021) developed a pseudospectral similarity solution for horizontal infiltration (i.e., solving RE without gravity) that is fast and accurate. This solution assumes a semi-infinite horizontal soil column ($\mathrm{0}\le x<\mathrm{\infty }$) with a uniform initial condition ($\mathit{\psi }(t=\mathrm{0})={\mathit{\psi }}_{\mathrm{0}}$) and a type-I boundary condition on the left boundary ($\mathit{\psi }(x=\mathrm{0})={\mathit{\psi }}_{\mathrm{L}}$). The model was run for 100 min ($\mathrm{0}\le t\le$ 100 min). The solution can resolve very large gradients in saturation and ψ at the boundary that propagate into the soil rapidly – and as such this is another challenging problem for a numerical RE model to reproduce. We solved this problem for three soil types, namely, Hygiene sandstone, silt loam GE 3, and Beit Netofa clay, with properties from van Genuchten (1980) as listed in Table 3. S_S was set to 0, consistent with Mathias' solution. To configure our model for horizontal flow, it is necessary to remove the gravity term from the flux calculations (i.e., Eq. 11). We solved this problem for a left-hand boundary effective saturation of 0.99 and an initial saturation of 0.01. The grid is configured such that the wetting pulse does not reach the right-hand boundary over the simulation runtime.

Table 3 Soil hydraulic properties used the Miller et al. (1998) problem.

Figure 5 Comparison of our ATS solver flux output model for horizontal infiltration with the Mathias and Sander (2021) pseudospectral similarity solution (denoted analytical in the legend).

Table 4 Runtime and water balance performance for the Celia MPM and ATS SFOM solutions, applied to the problem described by Mathias and Sander (2021).

We solved this problem with our ATS solution, with Celia's MPM solution, and with the pseudospectral similarity solution (Mathias and Sander, 2021, implemented in MATLAB). The results in Fig. 5 show that both the ATS solution and the Celia solution do an excellent job of reproducing this solution for θ(tx) (where x, m, is horizontal distance). The runtimes and water balance for each solution are tabulated in Table 4. Here, we see that the ATS and Celia solutions have the same performance in terms of the water balance and the cumulative fluxes simulated. Runtimes vary between models: both are the same for sandstone, Celia's solution is faster for silt loam, and the ATS solution is faster for clay.

3.2 Comparison with Hydrus 1D

Hydrus 1D (Šimůnek et al., 2005, 2016) is a widely used one-dimensional RE solver. The calculations within Hydrus are undertaken using openly available FORTRAN source code, and the software runs through a (closed-source) graphical user interface on Microsoft Windows. The FORTRAN code can be compiled using gfortran on the macOS operating system and run from the command line, which we did here, so that the runtime comparisons with our model are fair. Within Hydrus, the user interface provides somewhat limited control over the error tolerances. We were unable to modify any settings to improve the water balance, and so we present here model runs that use the default iteration criteria (maximum number of iterations = 100; water content tolerance = 0.001; pressure head tolerance = 10 mm; lower/upper optimal iteration range = 0.7/1.3; lower/upper time step multiplication factor = 1.3/0.7; lower/upper limit of tension interval = 10⁻⁶/10³ cm).

Table 5 Comparison of model runtimes and mass balance performance.

Figure 6 Comparison of water balance performance from Hydrus 1D, the ATS solution, and our implementation of Celia's MPM solution. The water balance error is $q_{\mathrm{T}}-q_{\mathrm{B}}-\mathrm{\Delta }S$, and we show the balance for each time step (a) and cumulative balances since the start of the simulation (b). In panels (c) and (d) the time step error is plotted against the infiltration rate.

Figure 7 Simulated storage (a) and drainage (b) using Hydrus 1D, the ATS solution, and our implementation of Celia's MPM solution.

We configured Hydrus 1D, our ATS solution, and our implementation of Celia's MPM solution for a simple numerical experiment, where we simulate the infiltration of a 10-year time series of daily precipitation into a 1.5 m deep soil column, with a free-drainage lower-boundary condition. The minimum, mean, and maximum annual precipitation was 265, 484, and 680 mm yr⁻¹, and the maximum daily precipitation was 55 mm d⁻¹. We used silt loam GE 3 soil hydraulic properties from van Genuchten (1980), where θ_r = 0.131, θ_s = 0.396, α = 0.423 m⁻¹, K_s = 0.0496 m d⁻¹, and n=2.06. We set S_S = 10⁻⁶ m⁻¹ and used a uniform ψ initial condition of −3.59 m. The results of the simulations with each model are given in Table 5, showing the runtime and water balance performance; Fig. 6, showing the detailed water balance performance; and Fig. 7, showing the simulated storage and drainage.

In our ATS solution, we can trade off between water balance error and the runtime by modifying the rtol argument for the ODE solver. We found that the default rtol of 10⁻⁶ had the fastest runtime, but the water balance performance, whilst good enough for all practical purposes, was the worst overall (Table 5, Fig. 6). Therefore, we reduced rtol to 10⁻⁷, which improved the water balance performance but increased the runtime. Even though the water balance errors reported here are all very small, it is still interesting to look closely at how these compare for the different models, as shown in Fig. 6. The first thing to note is that the Celia MPM solution has water balance errors of essentially zero, which we expect, because this solution enforces water balance closure. Celia's solution did have the longest runtime – approximately 40 % slower than the other solutions. In the ATS solutions, on a daily basis, the water balance errors are much smaller than in Hydrus. However, in Hydrus, the water balance errors appear random, with a mean of zero, and hence when looking at the cumulative errors in Hydrus, there is no systematic accumulation in the errors. In the ATS solutions, in the lower-right panel of Fig. 6, we can see that the water balance errors are strongly correlated to the infiltration flux at the upper boundary – larger fluxes result in larger errors. Hence, for the ATS solution with rtol = 10⁻⁶, we see that the errors accumulate, and after 4 years, the cumulative errors in the ATS solution exceed those in Hydrus. For the ATS solution with rtol = 10⁻⁷, the errors do not accumulate monotonically, and the long-term cumulative errors tend to oscillate about zero. The water balance performance of the ATS solution with rtol = 10⁻⁷ is therefore better than the performance in Hydrus (Table 5, Fig. 6), while the runtimes of these models are essentially the same (Hydrus is slightly faster, with a runtime of 2.21 vs. 2.30 s, Table 5).

Looking at the simulated storage and discharge in Fig. 7, the two ATS solutions are visually indistinguishable and are both broadly consistent with the Hydrus 1D model outputs. The Celia MPM solution has non-negligible differences with all other solutions. This is because the MPM solution applies an iterative solution procedure to solve the model at a daily time step, and the boundary fluxes are therefore subject to errors, as discussed above. The solution scheme imposes mass balance on the problem but does not track the truncation errors in the fluxes. The avoidance of this issue represents a significant advantage of adaptive time-stepping solutions.

4 Summary and recommendations

We developed a simple adaptive time-stepping scheme (ATS) for RE using the interpreted language Python and making use of the SciPy ODE solver ode. We also developed a new solver flux output method (SFOM) whereby cumulative boundary fluxes can be included within the dependent variable vector, allowing the determination of highly accurate integrated fluxes over designated time periods. The SFOM is particularly useful for providing reliable assessment of mass balance closure. In principle, SFOM can be implemented in any ODE solver because it does not require any special output (such as dense output) to be available. Our model was coded up in Python and released with the name openRE (Ireson, 2022). Our model performed well against our five success criteria: (i) we successfully reproduced benchmark solutions for ψ(t,z) and θ(tz) from Celia et al. (1990), Miller et al. (1998), and Mathias and Sander (2021); (ii) we report negligibly low mass balance errors; (iii) we simulate boundary fluxes that are independent of the reporting time step (unlike Celia's solution, as demonstrated in Fig. 7); (iv) we have low runtimes (as good as Hydrus 1D); and (v) our code is very simple, concise (92 lines of code for the solver plus 68 lines of code for model configuration for the numerical experiment in Sect. 3.2), and easily adaptable to new problems. Our solution had the best balance of efficiency, accuracy, and simplicity as compared to alternative established solution procedures.

Appendix A: Navigating pitfalls in ODE solver solutions

There are several subtle decisions that must be made when solving RE using a generic ODE solver. Here, we test a number of alternative model configurations and report the impact of these decisions using the following metrics: for model accuracy (criteria i), we report the RMSE of ψ at all grid points in t and z between the current model run and a reference model run; for the mass balance (criteria ii), we report both the bias error (Eq. 15) and the more rigorous daily water balance RMSE (Eq. 16); for the model efficiency (criteria iv), we simply report the runtime, where all runs were undertaken on the same laptop computer. For these numerical experiments, we used the 10-year infiltration numerical experiment described in Sect. 3.2.

Figure A1 Water balance performance plot for 10-year infiltration experiment, with the best model configuration, showing the cumulative change in storage in the profile and the cumulative inflow minus outflow. The water balance bias was −0.018 mm, and the RMSE of the daily water balance errors was 8.06 × 10⁻⁵ mm.

The best model configuration, against which all other model configurations are compared, was as follows: use the SciPy ODE solver ode with the method BDF (backward differentiation formula, Brown et al., 1989), use our SFOM solution (Sect. 2.3/A2), use the analytical expression for C(ψ), use a banded Jacobian sparsity pattern matrix (Sect. 2.4.1/A4), and use the Numba JIT compiler (Sect. 2.4.2/A5). The water balance performance of this model, showing the cumulative change in storage against cumulative inflow (as infiltration at the surface) minus outflow (as drainage at the base), is plotted in Fig. A1.

A1 Alternative SciPy ODE solvers

Here, we compare the alternative ODE solvers that were available in SciPy at the time of writing, which includes ode, odeint, and solve_ivp. These functions are alternative wrappers to classic ODE solvers written in Fortran, of which we consider here VODE with the method BDF (Brown et al., 1989, available within ode and solve_ivp) and LSODA (available with all three functions, Petzold, 1983). Note that for all solutions reported here we used the banded Jacobian sparsity pattern, with the exception of solve_ivp BDF, which only allows for the full Jacobian sparsity pattern to be defined and which we found slowed the solution down – hence the results for the solve_ivp BDF model do not use any information about the Jacobian matrix.

Table A1 Model performance for the different ODE solvers/methods available in SciPy.

We see that solve_ivp underperforms in accuracy, water balance, and efficiency. The odeint solver has the best performance in terms of accuracy and water balance but is slower by a non-negligible amount. The ode BDF method is the most efficient but has slightly worse water balance performance – however, the water balance performance of all methods is extremely good, and errors are negligible for practical purposes. We therefore chose ode BDF as our preferred solution – but ode LSODA is also a good option. It is also possible to increase the error tolerances in the ODE solver, reduce the maximum number of time steps, and increase the minimum time step – all of which could result in a faster runtime at the cost of lower accuracy/water balance closure.

A2 Alternative boundary flux calculation methods

As detailed in Sect. 2.3, there are alternative ways to calculate the boundary fluxes for use in the water balance calculation. In Sect. 3.3, we developed a novel approach to calculating the boundary fluxes – the SFOM. In addition to this method, we consider methods that calculate the boundary fluxes based on the output model states at either reporting-step or calculation-step information. We also consider using forward, backward, or central difference approximations to integrate the flux over a time step (Eq. 18). The results of this analysis are provided in Table A2.

Table A2 Model performance using different calculation methods for the boundary fluxes.

Model state variables are unaffected by the different boundary flux calculation methods. The SFOM has the best water balance performance, both in terms of bias and RMSE. Using calculation-step-level information results in good water balance closure, with the central difference approximation giving the lowest errors. However, the efficiency of this is poor, with runtimes increased by more than a factor of 3. This is because many additional calculations need to be performed outside the ODE solver for each calculation time step. By default, the ODE solver allows up to 500 calculation time steps for every reporting time step – so this is very inefficient. Calculating the boundary fluxes using reporting-time-step information is very efficient and slightly faster than our method, but the water balance errors are significantly larger. These reporting-step errors will increase with an increased reporting time step, as is shown in Fig. 3. Overall then, the SFOM provides the performance of using calculation-step information without the loss of computational efficiency.

The key take-home point here is that the easiest and most obvious approach to calculating the boundary fluxes is to use reporting-step information. This is a bad idea – the mass balance errors are large, and if this is combined with other bad decisions (such as using discrete approximations for C(ψ) as discussed in the next section), the results can be catastrophic (water balance errors > 100 mm).

A3 Alternative estimation methods for $\mathrm{d}\mathit{\theta }/\mathrm{d}\mathit{\psi }$

When we use a parametric expression for θ(ψ), such as the van Genuchten equations (Eqs. 26–29), we can obtain an analytical expression for $\mathrm{d}\mathit{\theta }/\mathrm{d}\mathit{\psi }$, as in Eq. (28), and this can be used to calculate C(ψ) as implemented in RE in Eq. (5). However, depending on the numerical solution procedure that is adopted, this can lead to errors with mass conservation, and it is recommended by some researchers (Rathfelder and Abriola, 1994; Clark et al., 2021) that a discrete approximation is used for $\mathrm{d}\mathit{\theta }/\mathrm{d}\mathit{\psi }$, whereby

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}\mathit{\psi }}}={\displaystyle \frac{{\mathit{\theta }}_n-{\mathit{\theta }}_{n-\mathrm{1}}}{{\mathit{\psi }}_n-{\mathit{\psi }}_{n-\mathrm{1}}}},\end{array}$$

where here n is a time index. This approach could be seen as equivalent to solving the mixed form of RE and can minimize water balance errors in the model that arise because the changes in $\mathrm{d}\mathit{\theta }/\mathrm{d}\mathit{\psi }$ over a time step are non-linear, as shown by Celia et al., (1990). However, it is necessary to apply this very carefully in the context of ATS methods. The values of θ_n−1 and ψ_n−1 must be available from the previous calculation time step and not the previous reporting time step. If reporting-time-step information is used, the model will fail badly because as the calculation steps move forward in time over a reporting step, C(ψ) is constantly referenced back to the beginning of the reporting step. This is clearly an erroneous approach, resulting in mass balance errors of more than 100 mm for our problem. For the solver flux output method of calculating the boundary fluxes, it is necessary to output states at reporting steps, and therefore it is not possible to use the discrete approximation for C(ψ). The results in Table A3 all use calculation-step information. A more subtle issue is the order of the temporal integrator used by the ODE solver, which can be specified by the user. Here, we use either first-order or (variable) higher-order (as determined by the ODE solver) temporal integration methods. For the solver flux output method, we use higher-order temporal integration. The results are given in Table A3.

Table A3 Model performance using different approaches to calculate C(ψ).

Table A4 Model performance using different approaches to define the Jacobian matrix.

Table A5 Model performance with and without Numba JIT compilation.

We see in Table 3 that the discrete C(ψ) approach works quite well for first-order integration methods but is very slow. When higher-order integration methods are used, the model is faster, but the mass balance is chronically degraded. We think that this happens because with higher-order methods the model states evolve in a more complex manner (i.e., non-linear manner) over a calculation time step, so the linear approximation in Eq. (A1) is not good. It is noteworthy that the modeled ψ values were slightly modified using the discrete high-order approach. For comparison purposes, we looked at using analytical representations of C(ψ) with first-order and higher-order methods, and this time the higher-order methods performed better. Overall, we recommend against using discrete C(ψ) approximations, unless using a tailor-made ODE solver (such as Kavetski et al., 2001, 2002a, b).

A4 Alternative approaches to defining the Jacobian

As described in Sect. 4.1, providing the ODE solver with information about the Jacobian matrix is reported to improve the solution efficiency. Here we compare three approaches: no information provided about the Jacobian, defining the Jacobian pattern, and defining the full Jacobian matrix. For the last case, this was complex to define for our method, and therefore it was implemented for the high-order reporting-step solution procedure described in Sect. 4.3.2. The results are reported in Table A4.

We see that for both model configurations, defining the banded Jacobian sparsity pattern matrix led to improvements in performance of around 20 %. This is modest, but because it is trivial to define the banded matrix, this is worthwhile. For the reporting-step model, when we defined the full Jacobian matrix, this led to a very slight improvement in performance over the banded solution (1.58 s vs. 1.63 s). Defining the full Jacobian is challenging and requires an additional function/function call in the code – we therefore recommend against using the full Jacobian matrix and recommend instead defining the banded matrix.

A5 Running the model with and without Numba JIT compilation

The best model configuration was also run with and without the Numba JIT compiler, and the result is shown in Table A5. It can be seen that Numba has no impact on the model output (accuracy and mass balance are identical for each run) as expected, but using Numba improves the runtime by a factor of ∼ 15. All other model runs reported in this paper use Numba.