Richardson-Richards’ equation

\[ \require{physics} \]

1 Introduction

In this document, we present a mass conservative solution to Richardson-Richards’ equation. Richardson-Richards’ equation is known to be quite difficult to solve. Our aim is to use a derivation that is relatively simple, stable and fast, and that can be used as a backbone in our SciML implementation. The reason why we are not using existing implementations like SWAP or Hydrus-1D is the lack of modularity/granularity of these software tools. The aim is to replace part of this equation by machine learning models. We will start simple by replacing the top boundary condition (evapotranspiration based on Penman Monteith) by a simple artificial neural network (see here). But other variants seem also promising. For instance by replacing (parts of) Richardson-Richard’s equation by faster artificial neural networks.

2 Continuity equation

We will start with the continuity equation, which embodies the Law of conservation of mass:

\[ \frac{\partial \theta}{\partial t} = -\frac{\partial q}{\partial s} - S \tag{1}\]

where:

• \(\theta\): volume fraction of liquid (m³ m^-3);
• \(t\): time (s);
• \(q\): flux density, i.e., the volume of water passing through a plane, during a specific time interval (m³ m^-2 s^-1 = m s ^-1);
• \(S\): sink term (e.g., water uptake by roots m³ m^-3 s^-1 )
• \(s\): spatial coordinate in vertical direction (m, origin at the soil surface and positive upwards, following Koorevaar, Menelik, and Dirksen (1983)).

3 Darcy’s law

One dimensional vertical flow \(q\) in a porous medium is given by Darcy’s law:

\[ q = -K \frac{\partial H}{\partial s} = -K \frac{\partial (h + z)}{\partial s} = -K \left( \frac{\partial h}{\partial s} + 1\right) \tag{2}\]

where:

• \(H\): hydraulic head (m);
• \(h\): soil moisture pressure head (m);
• \(z\): gravitational head (m) and spatial coordinate in vertical direction (m, origin at soil surface and positive upwards);
• \(K\): hydraulic conductivity (m s^-1).

4 Richardson-Richards’ equation

Substituting Equation 2 into Equation 1 gives Richardson-Richards’ equation (Richardson 1922; Richards 1931): \[ \frac{\partial \theta}{\partial t} = \frac{\partial}{\partial s} \left[ K \left (\frac{\partial h}{\partial s} + 1 \right) \right] \tag{3}\]

This expression is also known as ‘Richards’ equation’. Its mixed form is given above as opposed to the \(h\)-based and \(\theta\)-based forms (Celia, Bouloutas, and Zarba 1990).

5 Numerical discretization

Rewriting Equation 3 gives:

\[ \frac{\partial \theta}{\partial t} - \frac{\partial}{\partial s} \left[K \frac{\partial h}{\partial s} \right] - \frac{\partial K}{\partial s} = 0 \tag{4}\]

5.1 Temporal discretization

Its backward (implicit) Euler approximation is:

\[ \frac{\theta^{n+1, m+1} - \theta^n}{\Delta_t} - \frac{\partial}{\partial s} \left[K^{n+1,m} \frac{\partial h^{n+1,m+1}}{\partial s} \right] - \frac{\partial K^{n+1,m}}{\partial s} = 0 \tag{5}\]

where superscripts \(n\) and \(m\) indicate the time step and the iteration step respectively. Note that these superscripts are indices and not exponents!

A first order Taylor series¹ approximation of \(\theta^{n+1,m+1}\) gives:

\[ \theta^{n+1,m+1} \approx \theta^{n+1,m} + \left. \frac{d\theta}{dh} \right|^{n+1,m} \left(h^{n+1,m+1} - h^{n+1,m} \right) \]

Substituting the differential (or specific) moisture capacity \(C = \frac{d\theta}{d h}\) into this approximation gives:

\[ \theta^{n+1,m+1} \approx \theta^{n+1,m} + \delta^m C^{n+1,m} \] where \(\delta^m = h^{n+1,m+1} - h^{n+1,m}\) is the iteration increment for timestep \(n+1\) between iterations \(m\) and \(m+1\). Hence,

\[ \frac{\partial \theta}{\partial t} \approx \frac{\theta^{n+1,m+1} - \theta^n}{\Delta_t} \approx \frac{C^{n+1,m}}{\Delta_t} \delta^m + \frac{\theta^{n+1,m} - \theta^n}{\Delta_t} \] Replacing the time-derivative of Equation 5 by this approximation gives:

\[ \frac{C^{n+1,m}}{\Delta_t} \delta^m + \frac{\theta^{n+1,m} - \theta^n}{\Delta_t} - \frac{\partial}{\partial s} \left[K^{n+1,m} \frac{\partial h^{n+1,m+1}}{\partial s} \right] - \frac{\partial K^{n+1,m}}{\partial s} = 0 \] Recalling that \(h^{n+1,m+1} = \delta^m + h^{n+1,m}\) gives: \[ \frac{C^{n+1,m}}{\Delta_t} \delta^m - \frac{\partial}{\partial s} \left[K^{n+1,m} \frac{\partial (\delta^m + h^{n+1,m})}{\partial s} \right] = \frac{\partial K^{n+1,m}}{\partial s} - \frac{\theta^{n+1,m} - \theta^n}{\Delta_t} \] It follows that \[ \begin{align} \frac{C^{n+1,m}}{\Delta_t} \delta^m - \frac{\partial}{\partial s} \left[K^{n+1,m} \frac{\partial \delta^m}{\partial s} + K^{n+1,m} \frac{\partial h^{n+1,m}}{\partial s}\right] =\\~\\ \frac{\partial K^{n+1,m}}{\partial s} - \frac{\theta^{n+1,m} - \theta^n}{\Delta_t} \end{align} \] and therefore \[ \begin{align} \frac{C^{n+1,m}}{\Delta_t} \delta^m - \frac{\partial}{\partial s} \left[K^{n+1,m} \frac{\partial \delta^m}{\partial s} \right] = \frac{\partial}{\partial s} \left[ K^{n+1,m} \frac{\partial h^{n+1,m}}{\partial s}\right] + \\~\\ \frac{\partial K^{n+1,m}}{\partial s} - \frac{\theta^{n+1,m} - \theta^n}{\Delta_t} \end{align} \tag{6}\]

5.2 Spatial discretization

Figure 1 schematically depicts the discetization of a soil profile in space and time.

\begin{tikzpicture}
    \draw[step=10mm,gray,ultra thin] (0,0) grid (7.5,9);
    \draw[thick,->] (0,9) -- (7.5,9) node[anchor=west] {$t$};
    \draw[thick,->] (0,9) -- (0,-0.5) node[anchor=north] {$s$};
    \draw[thick,<->] (1,-0.5) -- (2,-0.5) node[pos=0.5, anchor=north] {$\Delta_t$};
    \draw[thick,<->] (7.5,4) -- (7.5,5) node[pos=0.5, anchor=west] {$\Delta_s$};
    \draw[gray](0,9) -- (0,9.2) node[anchor=south] {$0$};
    \draw[gray](1,9) -- (1,9.2) node[anchor=south] {$1$};
    \draw[gray](2,9) -- (2,9.2) node[anchor=south] {$2$};
    \draw[gray](3,9) -- (3,9.2) node[anchor=south] {$\dots$};
    \draw[gray](4,9) -- (4,9.2) node[anchor=south] {$n-1$};
    \draw[gray](5,9) -- (5,9.2) node[anchor=south] {$n$};
    \draw[gray](6,9) -- (6,9.2) node[anchor=south] {$n+1$};
    \draw[gray](7,9) -- (7,9.2) node[anchor=south] {$\dots$};
    \draw[gray](0,9) -- (-0.2,9) node[anchor=east] {$0$};
    \draw[gray](0,8) -- (-0.2,8) node[anchor=east] {$1$};
    \draw[gray](0,7) -- (-0.2,7) node[anchor=east] {$2$};
    \draw[gray](0,6) -- (-0.2,6) node[anchor=east] {$\vdots$};
    \draw[gray](0,5) -- (-0.2,5) node[anchor=east] {$i-1$};
    \draw[gray](0,4) -- (-0.2,4) node[anchor=east] {$i$};
    \draw[gray](0,3) -- (-0.2,3) node[anchor=east] {$i+1$};
    \draw[gray](0,2) -- (-0.2,2) node[anchor=east] {$\vdots$};
    \draw[gray](0,1) -- (-0.2,1) node[anchor=east] {$n_i$};
    \draw[gray](0,0) -- (-0.2,0) node[anchor=east] {$n_i+1$};
    \foreach \s in {0,...,9}
        \filldraw[red](0,\s) circle (2pt);
    \foreach \t in {1,...,7}
        \filldraw[red](\t,9) circle (2pt);
    \foreach \t in {1,...,7}
        \filldraw[red](\t,0) circle (2pt);
    \foreach \s in {1,...,8}{
        \foreach \t in {1,...,7}
            \filldraw[blue] (\t,\s) circle (2pt);
    }
\end{tikzpicture}

Figure 1: Numerical discretization of Richardson-Richards’ equation in space (vertically) and time (horizontally). The vertical dimension is postive upward. Boundary conditions are given as red dots, unknown pressure heads as blue dots.

The spatial discretization of individual components of Equation 6 is:

\[ \begin{align} \frac{\partial}{\partial s} \left[K^{n+1,m} \frac{\partial \delta^m}{\partial s} \right] \approx \frac{K_{i+\textonehalf}^{n+1,m} \frac{\delta_{i+1}^m - \delta_i^m}{\Delta_s} - K_{i-\textonehalf}^{n+1,m} \frac{\delta_i^m - \delta_{i-1}^m}{\Delta_s}}{\Delta_s} = \\~\\ \frac{K_{i+\textonehalf}^{n+1,m} \left(\delta_{i+1}^m - \delta_i^m \right) - K_{i-\textonehalf}^{n+1,m} \left(\delta_i^m - \delta_{i-1}^m \right)}{\Delta_s^2} \end{align} \]

and

\[ \frac{\partial}{\partial s} \left[K^{n+1,m} \frac{\partial h^{n+1,m}}{\partial s} \right] \approx \frac{K_{i+\textonehalf}^{n+1,m} \left(h_{i+1}^m - h_i^m \right) - K_{i-\textonehalf}^{n+1,m} \left(h_i^m - h_{i-1}^m \right)}{\Delta_s^2} \]

where \(i = \frac{s_i}{\Delta_s}\) is the vertical space index, and \(s_i\) the vertical space coordinate \(s_i = i \Delta_s\).

Substituting these components in Equation 6 gives: \[ \begin{align} \frac{\delta_i^m C_i^{n+1,m}}{\Delta_t} - \frac{K_{i+\textonehalf}^{n+1,m} \left(\delta_{i+1}^m - \delta_i^m \right) - K_{i-\textonehalf}^{n+1,m} \left(\delta_i^m - \delta_{i-1}^m \right)}{\Delta_s^2} = \\ ~\\~ \frac{K_{i+\textonehalf}^{n+1,m} \left(h_{i+1}^{n+1,m} - h_i^{n+1,m} \right) - K_{i-\textonehalf}^{n+1,m} \left(h_i^{n+1,m} - h_{i-1}^{n+1,m} \right)}{\Delta_s^2} + \\ ~\\~ \frac{K_{i+\textonehalf}^{n+1,m} - K_{i-\textonehalf}^{n+1,m}}{\Delta_s} - \frac{\theta_i^{n+1,m} - \theta_i^n}{\Delta_t} \equiv R_i^{n+1,m} \end{align} \tag{7}\]

Upon convergence, iteration increments \(\delta^m\) approach zero, and so is the residual associated with the Picard iteration \(R_i^{n+1,m}\).

This solution to Richardson-Richards’ equation is more complicated than the one presented in Feddes, Kowalik, and Zaradny (1978), but is more stable as global mass conservation is guaranteed (Celia, Bouloutas, and Zarba 1990).

6 Implementation

6.1 Matrix notation

The left-hand side in Equation 7 can be rewritten as: \[ \begin{align} \frac{\delta_i^m C_i^{n+1,m}}{\Delta_t} - \frac{K_{i+\textonehalf}^{n+1,m} \left(\delta_{i+1}^m - \delta_i^m \right) - K_{i-\textonehalf}^{n+1,m} \left(\delta_i^m - \delta_{i-1}^m \right)}{\Delta_s^2} = \\ ~\\~ \frac{C_i^{n+1,m}}{\Delta_t}\delta_i^m - \frac{K_{i+\textonehalf}^{n+1,m}}{\Delta_s^2} \delta_{i+1}^m + \frac{K_{i+\textonehalf}^{n+1,m}}{\Delta_s^2} \delta_i^m + \frac{K_{i-\textonehalf}^{n+1,m}}{\Delta_s^2} \delta_i^m - \frac{K_{i-\textonehalf}^{n+1,m}}{\Delta_s^2} \delta_{i-1}^m = \\ ~\\~ -\frac{K_{i-\textonehalf}^{n+1,m}}{\Delta_s^2} \delta_{i-1}^m + \left( \frac{C_i^{n+1,m}}{\Delta_t} + \frac{K_{i-\textonehalf}^{n+1,m}}{\Delta_s^2} + \frac{K_{i+\textonehalf}^{n+1,m}}{\Delta_s^2} \right) \delta_i^m - \frac{K_{i+\textonehalf}^{n+1,m}}{\Delta_s^2} \delta_{i+1}^m\\ \end{align} \]

Let

\[ G_i^{n+1,m} = -\frac{\Delta_s^2}{\Delta_t} C_i^{n+1,m} - K_{i-\textonehalf}^{n+1,m} - K_{i+\textonehalf}^{n+1,m} \]

then it follows that

\[ \mqty[ G_1^{n+1,m} & K_{1\textonehalf}^{n+1,m} & \cdot & \cdot & \cdot & \cdot & \cdot \\ K_{1 \textonehalf}^{n+1,m}& G_2^{n+1,m} & K_{2 \textonehalf}^{n+1,m} & \cdot & \cdot & \cdot & \cdot \\ \cdot & K_{i-\textonehalf}^{n+1,m} & G_i^{n+1,m} & K_{i+\textonehalf}^{n+1,m} & \cdot & \cdot & \cdot \\ \cdot & \cdot & \ddots & \ddots & \ddots & \vdots & \cdot \\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \cdot\\ \cdot & \cdot & \cdot & \cdot & K_{n_i-1\textonehalf}^{n+1,m} & G_{n_i-1}^{n+1,m} & K_{{n_i}-\textonehalf}^{n+1,m}\\ \cdot & \cdot & \cdot & \cdot & \cdot & K_{n_i-\textonehalf}^{n+1,m} & G_{n_i}^{n+1,m} ] \mqty[ \delta_1^m \\ \delta_2^m \\ \vdots \\ \delta_{n_i}^m] = -\Delta_s^2 \mqty[ r_1 \\ r_2 \\ \vdots \\ r_{n_i}] \]

6.2 Boundary conditions

6.2.1 Initial boundary condition

Initial boundary conditions are given as red dots along the vertical axis in Figure 1. It is recommended to use a burn-in period to eliminate the effect of the initial boundary condtion on the final results. As long as the burn-in period is sufficiently long, we may select a constant pressure head or a moisture profile that is in equilibrium with the groundwater table depth.

6.2.2 Top and bottom boundary conditions

Since we are solving for iteration increments it follows that \(\delta_0^m = \delta_{n_i + 1}^m = 0\). Indeed we don’t want to iteratively modify these (known) boundary condition. We only need to specify \(h_0^{n+1, m}\) and \(h_{n_i + 1}^{n+1, m}\) in vector \(\vb{r}\).

7 References

Braun, M. 1993. Differential Equations and Their Applications. An Introduction to Applied Mathematics. Springer.

Celia, M. A., E. T. Bouloutas, and R. L. Zarba. 1990. “A General Mass-Conservative Numerical Solution for the Unsaturated Flow Equation.” Water Resources Research 26 (7): 1483–96. https://doi.org/https://doi.org/10.1029/WR026i007p01483.

Feddes, R. A., P. J. Kowalik, and H. Zaradny. 1978. Simulation of Field Water Use and Crop Yield. Simulation Monographs. Pudoc.

Koorevaar, P., G. Menelik, and C. Dirksen. 1983. Elements of Soil Physics. Elsevier Science Publishers.

Richards, L. A. 1931. “Capillary Conduction of Liquids Through Porous Mediums.” Physics 1 (5): 318–33. https://doi.org/10.1063/1.1745010.

Richardson, L. F. 1922. Weather Prediction by Numerical Process. Cambridge, The University press.

Footnotes

1. A first order Tailor series approximation of a function \(f(x)\) with respect to \(a\) is given by \(f(x) \approx f(a) + \frac{f'(a)}{1!}(x - a)\) or expressed as a time derivative \(y(t_k + h) \approx y(t_k) + h \frac{dy(t_k)}{dt}\) (Braun 1993). The backward or implicit Euler method is based on this approximation and is given by \(y(t_{k + 1}) = y(t_k) + hf(t_{k+1}, y_{k+1})\) where \(f(t, y)=\frac{dy}{dt}\) or \(\frac{y(t_{k + 1}) - y(t_k)}{h} = f(t_{k+1}, y_{k+1})\)