# Implicit and semi-implicit second-order time stepping methods for the
Richards equation

Sana Keita<EMAIL_ADDRESS>Abdelaziz Beljadid<EMAIL_ADDRESS>Yves Bourgault<EMAIL_ADDRESS>International Water Research Institute,
Mohammed VI Polytechnic University, Morocco  Department of Mathematics and
Statistics, University of Ottawa, Canada

###### Abstract

This study concerns numerical methods for efficiently solving the Richards
equation where different weak formulations and computational techniques are
analyzed. The spatial discretizations are based on standard or mixed finite
element methods. Different implicit and semi-implicit temporal discretization
techniques of second-order accuracy are studied. To obtain a linear system for
the semi-implicit schemes, we propose second-order techniques using
extrapolation formulas and/or semi-implicit Taylor approximations for the
temporal discretization of nonlinear terms. A numerical convergence study and
a series of numerical tests are performed to analyze efficiency and robustness
of the different schemes. The developed scheme, based on the proposed temporal
extrapolation techniques and the mixed formulation involving the saturation
and pressure head and using the standard linear Lagrange element, performs
better than other schemes based on the saturation and the flux and using the
Raviart-Thomas elements. The proposed semi-implicit scheme is a good
alternative when implicit schemes meet convergence issues.

###### keywords:

Richards equation , Linearization schemes , Numerical convergence analysis ,
Second-order time-accuracy , Galerkin finite elements , Mixed finite elements

††journal:

[table]capposition=below

##  1 Introduction

The prediction of flows through unsaturated porous media is of great
importance in a wide variety of applications in science and engineering. There
are plenty of relevant applications in many fields such as hydrology,
agriculture, water resources management, risk assessment of environmental
contaminants, enhanced oil recovery, etc. The Richards equation [63] describes
fluid flows in unsaturated porous medium due to the actions of gravity and
capillarity. This equation was derived based on a multiphase flow extension of
Darcy’s law and the principle of mass conservation [10, 63]. Some analytical
solutions of the Richards equation are available, but they are limited to very
simple geometries and specific initial and boundary conditions. Most practical
situations require the use of numerical techniques to produce approximate
solutions. The design and analysis of appropriate numerical schemes for the
Richards equation is a very challenging task due to the highly nonlinear
nature of the equation. Despite intensive research, there is still a strong
need for more robust numerical schemes for computing multiphase flows in
unsaturated porous media.

Much research has been done on the numerical treatment of the Richards
equation. Various classes of approaches are available for the spatial
discretization. Standard finite element methods [3, 32, 69] and mixed finite
element methods [4, 8, 9, 12, 13, 30, 58, 60, 61, 67] are commonly used, often
combined with mass lumping to ensure nonoscillatory solutions [12, 21].
Locally mass-conservative finite volume methods for the Richards equation were
proposed in [29, 31, 48, 49, 33], while others studies have used finite
difference methods [3, 21, 23, 36, 50] and a control volume finite element
method [53]. Another well-known numerical approach used to solve the Richards
equation is the Kirchhoff transformation technique [4, 14, 15, 58, 60, 71].
While this technique is very promising, it is limited to specific functions of
the capillary pressure [71]. In terms of temporal discretization techniques, a
proper treatment of the temporal derivative is required for reliable numerical
solutions due to the highly nonlinear and stiff nature of the equation. Most
of the numerical models employ an implicit or semi-implicit scheme for the
temporal discretization. The main reason is the need of a stable
discretization allowing reasonable time steps since, due to the presence of a
stiff nonlinear diffusion term, an explicit scheme requires very small time
step which makes this scheme impractical to use in terms of efficiency. The
backward Euler method is the most commonly used [3, 4, 8, 9, 12, 29, 32, 53,
60, 61, 67] and it results in a nonlinear system. The conventional approach
consists in solving this nonlinear system using an iterative procedure. There
are several iterative methods available for the linearization of the resulting
nonlinear system. The Newton method, which is quadratically convergent, is
often used. However, its convergence occurs only when one starts with an
initial guess which is very close to the solution. A method for Richards’
equation has been proposed in [20] that exploits the typical behavior of the
soil hydraulic functions to define nested iterations that guarantee mass
conservation and quadratic convergence of a Newton-type method. The Picard
method, which is easy to implement, is also widely used, although it is only
linearly convergent. A modified Picard method is used to develop a mass-
conservative numerical scheme in [21]. In [47], an acceleration technique is
applied to the modified Picard method to provide faster convergence while
maintaining the advantages, i.e., the mass conservation and lower memory
requirements with respect to Newton’s method. A comparison between Picard and
Newton iteration approaches and a combination of the modified Picard and
Newton methods, as well as $L$-scheme and its combination with Newton method,
can be found in [43, 46, 55].

Newton and Picard methods can produce good results but these methods entail
computational costs associated with having to solve a nonlinear system at each
time step. Noniterative methods require a single matrix assembly and
factorization per time step and thus are an attractive alternative to
traditional iterative methods for the Richards equation [40, 54, 64, 79, 80].
Second order time discretization has shown benefits over traditional first
order temporal approximation for modeling unsaturated flows through porous
media [30, 40, 54, 80]. Furthermore, in our last studies [42, 41] we showed
the advantages of semi-implicit techniques in solving highly nonlinear
parabolic equations. This motivates an investigation of second order
noniterative approaches to solving the nonlinear Richards equation. This study
focuses on finite element methods for the Richards equation where three
different weak formulations of the equation are considered. The first one is
expressed in terms of the saturation as the only unknown. The second is a
mixed formulation based on the saturation and the flux. This mixed formulation
is mostly studied and the flux is approximated using the Raviart-Thomas
elements [9, 12, 13, 30, 58, 61] or the Brezzi-Douglas-Marini element [8]. The
third one is also a mixed formulation, but it consists in approximating
simultaneously the saturation and the pressure head by using the standard
linear Lagrange element. All the formulations include a regularization
technique to handle degeneracies that could occur in the gradient of the
pressure head for some constitutive relationships of the capillary pressure
when the saturation reaches the value of one. For each weak formulation,
iterative and noniterative second-order methods are studied. We use the
second-order backward differentiation formula (BDF2) for the time integration
with fixed time step. We note that there exists adaptive BDF2 methods which
allow the use of varying time steps [11, 25, 27]. In [7], the fully implicit
adaptive BDF2 method is applied to Richards’ equation. Our schemes can be
modified to allow varying time steps by using the techniques developed in [7,
11, 25, 27], but this would require modifications in the treatment of the
nonlinear terms since some of our schemes are based on semi-implicit BDF2
methods. To obtain a linear system for the noniterative schemes, a second-
order extrapolation formula and/or a semi-implicit Taylor expansion of second-
order truncation error are used for the nonlinear terms. To the best of our
knowledge, the linear mixed scheme based on the saturation and the pressure
head, using the standard Lagrange linear element for the spatial
discretization and the linearization techniques based on Taylor expansion, has
not been explored in previous studies. A main goal is to determine which weak
formulation using such temporal discretization is the most competitive among a
relatively large set of methods. Here, competitive is to be understood in
terms of accuracy versus computational requirements, but also ease of
implementation. Numerical convergence analysis is performed using manufactured
and reference solutions to investigate the spatial and temporal order of
convergence of each scheme.

The outline of the paper is as follows. In section 2, the Richards equation is
introduced. In section 3, three different weak formulations of the Richards
equation including a regularization technique are presented. In section 4,
several numerical schemes are proposed. In section 5, convergence rates,
accuracy and robustness of the numerical schemes are investigated and compared
by means of manufactured and exact solutions and selected test cases. Section
6 provides concluding remarks.

##  2 Model equations

Let $\Omega\subset\mathbb{R}^{2}$ be an open bounded set with sufficiently
smooth boundary $\Gamma$ and $T>0$ a fixed time. The Richards equation is a
standard, commonly-used continuum model for describing flow in unsaturated
porous media. It can be written in three different forms, namely the “pressure
head form”, “water content form” and “mixed form” [21]. The mixed form is
given by

$\dfrac{\partial\theta}{\partial
t}-\nabla\cdot[K_{s}(\bm{x})K_{r}(\psi)\nabla(\psi+z)]=0\quad\text{in}\quad(0,T)\times\Omega,$
(2.1)

where $\theta$ is the water content, $\psi$ is the pressure head, $K_{s}$ is
the saturated hydraulic conductivity, $K_{r}$ is the water relative
permeability, and $z$ is the vertical coordinate and is assumed to be positive
upward. The water relative permeability $K_{r}$ accounts for the effect of
partial saturation. While $K_{r}$ and $\psi$ include physical parameters which
depend on the spatial heterogeneity of the medium, for simplicity we write
$K_{r}(\psi)$. To complete the model formulation of the Richards equation
(2.1), one must specify the constitutive relationship to describe the
interdependence between the relative permeability, pressure head and water
content. There are several empirical relationships that are used in modeling,
for instance, van Genuchten [76], Brooks-Corey [18] and Gardner [35] models.

By introducing the effective saturation

$S=\dfrac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}},$ (2.2)

where $\theta_{s}$ and $\theta_{r}$ are the saturated and residual volumetric
water contents, respectively, the Richards equation (2.1) can be written as

$\phi(\bm{x})\dfrac{\partial S}{\partial
t}-\nabla\cdot[K_{s}(\bm{x})K_{r}(\psi)\nabla(\psi+z)]=0\quad\text{ in
}(0,T)\times\Omega,$ (2.3)

where $\phi$ is given by

$\phi=\theta_{s}-\theta_{r}.$ (2.4)

The water flux $\bm{q}$ is expressed by the extended form of Darcy’s law [57]:

$\bm{q}=-K_{s}K_{r}\nabla(\psi+z).$ (2.5)

Equation (2.3) is to be supplemented by an initial condition

$S(\bm{x},0)=S_{0}(\bm{x})$ (2.6)

and boundary conditions

$\bm{q}\cdot\bm{n}=q_{\Gamma}\quad\text{on}\quad\Gamma,$ (2.7)

where $\bm{n}$ is the outward unit normal vector to the boundary $\Gamma$ of
the domain $\Omega$ and $q_{\Gamma}\in\mathbb{R}$.

##  3 Weak formulations

For most constitutive relationships used in practice, the pressure head can be
expressed as

$\psi(S)=h_{cap}(\bm{x}){J}(S),$ (3.1)

where $h_{cap}:\Omega\rightarrow\mathbb{R}^{+}$ is the capillary rise function
which depends on space for heterogeneous medium and
$J:(0,1]\rightarrow(-\infty,0]$ is the Leverett $J$-function [45] which is of
class $\mathcal{C}^{1}$ on $(0,1)$. The Leverett $J$-function tends to
infinity as the effective saturation $S$ goes to zero, and for some
constitutive relationships, $J$ is not differentiable at $S=1$, for instance,
the constitutive relationship used by Haverkamp et al. [36] (see (5.4) below
for the explicit form of the $J$-function) and the van Genuchten model [76]
(see (5.15) below). The computations of solutions may encounter numerical
issues when the saturation $S$ approaches the value of zero or one due to the
lack of regularity of the constitutive functions. Some $J$-functions are
differentiable at $S=1$ such as for the Brooks-Corey model [18]. In case of
lack of differentiability, mathematical regularization techniques have been
proposed in [59, 68] to overcome the degeneracy in the pressure head. In this
work, we consider the regularization consisting in replacing the derivative
$J^{\prime}$ by $J_{\delta}^{\prime}$ which satisfies

$J_{\delta}^{\prime}(S)=\left\\{\begin{aligned} &J^{\prime}(S),\quad\forall
S\in(0,1-\delta),\\\ &J^{\prime}(1-\delta),\quad\forall
S\in[1-\delta,1],\end{aligned}\right.$ (3.2)

where $0<\delta\ll 1$. The regularized $J_{\delta}^{\prime}$ is a continuous
function on the interval $(0,1]$. This is very useful for the convergence of
numerical methods involving the computation of the derivative of $J$. The
$J$-function in (3.1) may include parameters which depend on the heterogeneity
of the medium. Here, the impact of the variability due to heterogeneity is
considered dominant by the variation of the intrinsic permeability included in
the capillary rise $h_{cap}$ [45]. By expanding (3.1), we get

$\nabla\psi=J(S)\nabla h_{cap}+h_{cap}J_{\delta}^{\prime}(S)\nabla
S,\quad\forall S\in(0,1].$ (3.3)

Substituting (3.3) in the Richards equation (2.3) leads to

$\displaystyle\phi\dfrac{\partial S}{\partial
t}-\nabla\cdot[K_{s}K_{r}(\psi(S))\big{(}J(S)\nabla
h_{cap}+h_{cap}J_{\delta}^{\prime}(S)\nabla S+\nabla z\big{)}]=0.$ (3.4)

Multiplying (3.4) by a test function $v$, integrating over the domain $\Omega$
and using (2.7) give rise to

$\displaystyle\int_{\Omega}\phi\dfrac{\partial S}{\partial t}v\,d\bm{x}$
$\displaystyle+\int_{\Omega}\big{[}K_{s}K_{r}(\psi(S))\big{(}J(S)\nabla
h_{cap}+h_{cap}J_{\delta}^{\prime}(S)\nabla S+\nabla
z\big{)}\big{]}\cdot\nabla v\,d\bm{x}$ (3.5)
$\displaystyle+\int_{\Gamma}q_{\Gamma}v\,d\bm{x}=0,\quad\forall v\in
H^{1}(\Omega).$

A weak formulation for (3.4) is: Find $S\in H^{1}(\Omega)$ such that (3.5) is
satisfied.

######  Remark 3.1.

Dirichlet conditions could be imposed on the boundary $\Gamma$. In this case,
the solution $S$ is imposed on the boundary $\Gamma$ and the test function $v$
is set to zero on $\Gamma$. The same idea applies to all the variational
formulations below.

By substituting (2.5) and (3.3), the Richards equation (2.3) can be written in
the mixed form

$\displaystyle\phi\dfrac{\partial S}{\partial t}+\nabla\cdot\bm{q}=0,$ (3.6)
$\displaystyle\bm{q}+K_{s}K_{r}(\psi)\big{(}J(S)\nabla
h_{cap}+h_{cap}J_{\delta}^{\prime}(S)\nabla S+\nabla z\big{)}=0.$

Multiplying each equation in (3.6) by a test function and integrating over the
domain lead to

$\displaystyle\int_{\Omega}\phi\dfrac{\partial S}{\partial
t}v\,d\bm{x}-\int_{\Omega}\bm{q}\cdot\nabla
v\,d\bm{x}+\int_{\Gamma}q_{\Gamma}v\,d\bm{x}=0,\quad\forall v\in
H^{1}(\Omega)$ (3.7)
$\displaystyle\int_{\Omega}\bm{q}\cdot\bm{w}\,d\bm{x}+\int_{\Omega}\big{[}K_{s}K_{r}(\psi(S))\big{(}J(S)\nabla
h_{cap}+h_{cap}J_{\delta}^{\prime}(S)\nabla S+\nabla
z\big{)}\big{]}\cdot\bm{w}\,d\bm{x}=0,\quad\forall\bm{w}\in(L^{2}(\Omega))^{2}.$

A weak formulation for (3.6) is: Find $(S,\bm{q})\in H^{1}(\Omega)\times
H(div;\Omega)$ such that (3.7) is satisfied, where

$H(div;\Omega)=\Big{\\{}\bm{w}\in L^{2}(\Omega)\times
L^{2}(\Omega):\nabla\cdot\bm{w}\in L^{2}(\Omega)\Big{\\}}.$ (3.8)

By using the pressure head relation (3.1), the Richards equation (2.3) can
also be written in the alternative mixed formulation

$\displaystyle\phi\dfrac{\partial S}{\partial
t}-\nabla\cdot[K_{s}K_{r}(\psi)\nabla(\psi+z)]=0,$ (3.9) $\displaystyle\psi-
h_{cap}J(S)=0.$

Multiplying each equation in (3.9) by a test function and integrating over the
domain give

$\displaystyle\int_{\Omega}\phi\dfrac{\partial S}{\partial
t}v\,d\bm{x}+\int_{\Omega}[K_{s}K_{r}(\psi)\big{(}\nabla\psi+\nabla
z\big{)}]\cdot\nabla
v\,d\bm{x}+\int_{\Gamma}q_{\Gamma}v\,d\bm{x}=0,\quad\forall v\in
H^{1}(\Omega),$ (3.10) $\displaystyle\int_{\Omega}\psi
w\,d\bm{x}-\int_{\Omega}h_{cap}J(S)w\,d\bm{x}=0,\quad\forall w\in
H^{1}(\Omega).$

A weak formulation for (3.9) is: Find $(S,\psi)\in H^{1}(\Omega)\times
H^{1}(\Omega)$ such that (3.10) is satisfied. The mixed formulation (3.10) is
quite simple and leads to a smaller linear system to solve than the mixed
formulation (3.7) as we will see in the next section. It is also advantageous
when it comes to the choice of the finite element spaces for the spatial
discretization since the standard linear Lagrange element is used instead of
the Brezzi-Douglas-Marini element or the Raviart-Thomas elements, which are
often used for the mixed formulation (3.7) and are more difficult to
implement.

##  4 Numerical methods

###  4.1 Temporal discretization

The time interval $[0,T]$ is discretized as

$\Delta t=\dfrac{T}{N},\quad t_{n}=n\Delta t,\quad n=0,1,2,...,N,$ (4.1)

where $\Delta t$ is the time step used. We set $S^{n}\simeq S(t_{n})$,
$\psi^{n}\simeq\psi(t_{n})$ and $\bm{q}^{n}\simeq\bm{q}(t_{n})$ for
$n=0,...,N$.
We propose implicit and semi-implicit time stepping discretizations for each
of the formulations (3.4), (3.6) and (3.9). Few second order accurate time-
stepping schemes were proposed for Richards equation based on the Crank-
Nicolson method [19, 79, 80] and BDF methods [7, 26] in their fully implicit
form requering Newton or Picard iterations. The Crank-Nicolson method is
A-stable but lacks L-stability, and may lead to non-monotone solutions for
larger time steps. Among the simplest methods for solving stiff ODEs, the
second-order backward differentiation formula (BDF2) is L-stable. A concise
discussion on these concepts and two time-stepping methods are presented in
[44]. Consequently, all the methods proposed here approximate the time
derivative using the BDF2 formula. We first consider the fully implicit BDF2
method, which is the most stable and accurate, but requires the solution of
stiff nonlinear systems at each time step. To avoid the solution of nonlinear
systems, semi-implicit (also called IMEX) variants of the BDF2 method are
available. These variants are usually based on the extrapolation of nonlinear
terms [1, 5, 28, 42], at the expense of a slightly more restrictive time step
condition for stability, although with time steps that are still independent
from the space step. The extension of these semi-implicit methods to Richards’
equation requires some care, since those were so far applied to PDEs with
other forms of nonlinear terms. The semi-implicit BDF2 method developed in
[42] for the Cahn-Hilliard equation can be extended for solving Richards’
equation. To obtain a linear scheme of second-order time-accuracy for the
semi-implicit temporal discretizations, a second-order extrapolation formula
is employed for the relative permeability function $K_{r}(\psi)$ and the
regularized derivative $J_{\delta}^{\prime}(S)$, while the Leverett
$J$-function $J(S)$ is evaluated using a Taylor approximation with second-
order truncation error.

####  4.1.1 Formulation (3.4)

A second-order implicit temporal discretization for (3.4) is given by

$\displaystyle\phi\dfrac{3S^{n+1}-4S^{n}+S^{n-1}}{2\Delta
t}-\nabla\cdot[K_{s}K_{r}(\psi(S^{n+1}))\big{(}J(S^{n+1})\nabla h_{cap}$
($D_{1i}$)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+h_{cap}J_{\delta}^{\prime}(S^{n+1})\nabla
S^{n+1}+\nabla z\big{)}]=0.$

For a second order semi-implicit temporal discretization of (3.4), we propose

$\displaystyle\phi\dfrac{3S^{n+1}-4S^{n}+S^{n-1}}{2\Delta
t}-\nabla\cdot\Big{(}\big{[}K_{s}\big{(}2K_{r}(\psi(S^{n}))-K_{r}(\psi(S^{n-1}))\big{)}\big{]}$
($D_{1s}$)
$\displaystyle\quad\big{[}\big{(}J(S^{n})+J_{\delta}^{\prime}(S^{n})(S^{n+1}-S^{n})\big{)}\nabla
h_{cap}+\big{(}2J_{\delta}^{\prime}(S^{n})-J_{\delta}^{\prime}(S^{n-1})\big{)}h_{cap}\nabla
S^{n+1}+\nabla z\big{]}\Big{)}=0.$

####  4.1.2 Formulation (3.6)

An implicit discretization of second-order accuracy for (3.6) reads

$\displaystyle\phi\dfrac{3S^{n+1}-4S^{n}+S^{n-1}}{2\Delta
t}-\nabla\cdot\bm{q}^{n+1}=0,$ ($D_{2i}$)
$\displaystyle\bm{q}^{n+1}+K_{s}K_{r}(\psi(S^{n+1}))\big{[}J(S^{n+1})\nabla
h_{cap}+h_{cap}J_{\delta}^{\prime}(S^{n+1})\nabla S^{n+1}+\nabla z\big{]}=0$

and a second-order semi-implicit discretization is given by

$\displaystyle\phi\dfrac{3S^{n+1}-4S^{n}+S^{n-1}}{2\Delta
t}-\nabla\cdot\bm{q}^{n+1}=0,$ ($D_{2s}$)
$\displaystyle\bm{q}^{n+1}+\big{[}K_{s}\big{(}2K_{r}(\psi(S^{n}))-K_{r}(\psi(S^{n-1}))\big{)}\big{]}\big{[}\big{(}J(S^{n})+J_{\delta}^{\prime}(S^{n})(S^{n+1}-S^{n})\big{)}\nabla
h_{cap}$
$\displaystyle\qquad\qquad\qquad\qquad+\big{(}2J_{\delta}^{\prime}(S^{n})-J_{\delta}^{\prime}(S^{n-1})\big{)}h_{cap}\nabla
S^{n+1}+\nabla z\big{]}=0.$

####  4.1.3 Formulation (3.9)

For the formulation (3.9), we consider the second order implicit time
integration

$\displaystyle\phi\dfrac{3S^{n+1}-4S^{n}+S^{n-1}}{2\Delta
t}-\nabla\cdot[K_{s}K_{r}(\psi^{n+1})\nabla(\psi^{n+1}+z)]=0,$ ($D_{3i}$)
$\displaystyle\psi^{n+1}-h_{cap}J(S^{n+1})=0,$

and the semi-implicit discretization of second-order accuracy

$\displaystyle\phi\dfrac{3S^{n+1}-4S^{n}+S^{n-1}}{2\Delta
t}-\nabla\cdot[K_{s}\big{(}2K_{r}(\psi^{n})-K_{r}(\psi^{n-1})\big{)}\nabla(\psi^{n+1}+z)]=0,$
($D_{3s}$)
$\displaystyle\psi^{n+1}-h_{cap}\big{(}J(S^{n})+J_{\delta}^{\prime}(S^{n})(S^{n+1}-S^{n})\big{)}=0.$

###  4.2 Spatial discretization

Let $\mathcal{T}_{h}$ be a triangulation of the domain $\Omega$ into disjoint
triangles $\kappa$, with $h_{\kappa}:=diam(\kappa)$ and
$h:=\max_{\kappa\in\mathcal{T}_{h}}h_{\kappa}$, so that

$\Omega=\bigcup_{\kappa\in\mathcal{T}_{h}}\kappa.$ (4.2)

Let $\mathbf{J}$ be the set of nodes of $\mathcal{T}_{h}$, with coordinates
$\\{\bm{x}_{j}\\}_{j\in\mathbf{J}}$. We consider the following finite element
spaces

$\displaystyle\mathcal{V}_{h}$
$\displaystyle=\Big{\\{}v_{h}\in\mathcal{C}^{0}(\Omega,\mathbb{R}):{v_{h}}_{|_{\kappa}}\text{
is linear,}\quad\forall\kappa\in\mathcal{T}_{h}\Big{\\}},$ (4.3)
$\displaystyle\mathcal{S}_{h}$
$\displaystyle=\Big{\\{}v_{h}\in{\mathcal{V}_{h}}:0<v_{h}(\bm{x}_{j})\leqslant
1,\quad\forall\bm{x}_{j}\in\mathcal{T}_{h}\Big{\\}},$ (4.4)
$\displaystyle\mathcal{R}_{h}$ $\displaystyle=\Big{\\{}\bm{v}_{h}\in
H(div;\Omega):{\bm{v}_{h}}_{|_{\kappa}}\in
RT_{0},\quad\forall\kappa\in\mathcal{T}_{h}\Big{\\}},$ (4.5)

where $RT_{0}$ is the Raviart-Thomas finite elements of degree $0$ [16]. In
what follows, we give a finite element method for each of the discretizations
in Section 4.1. The following equations are obtained from the ones in Section
4.1 by integrating over the domain, using the boundary condition (2.7) and
replacing the Sobolev spaces by the corresponding finite element spaces.

A finite element approximation for ($D_{1i}$) is: Given a suitable
approximation of the initial solution $S_{h}^{0}=\Pi_{h}S_{0}$ and a proper
initialization for $S_{h}^{1}\in\mathcal{S}_{h}$, find
$\widehat{S}_{h}^{n+1}\in{\mathcal{V}}_{h}$ such that

$\displaystyle\int_{\Omega}\phi_{h}\dfrac{3\widehat{S}_{h}^{n+1}-4S_{h}^{n}+S_{h}^{n-1}}{2\Delta
t}v_{h}\,d\bm{x}$ (4.6)
$\displaystyle\qquad+\int_{\Omega}\big{[}K_{s}K_{r}(\psi(\widehat{S}_{h}^{n+1}))\big{(}J(\widehat{S}_{h}^{n+1})\nabla
h_{cap}+h_{cap}J_{\delta}^{\prime}(\widehat{S}_{h}^{n+1})\nabla\widehat{S}_{h}^{n+1}+\nabla
z\big{)}\big{]}\cdot\nabla v_{h}\,d\bm{x}$
$\displaystyle\qquad+\int_{\Gamma}q_{\Gamma}v_{h}\,d\bm{x}=0,\quad\forall
v_{h}\in{\mathcal{V}}_{h}.$

Then, find $S_{h}^{n+1}\in\mathcal{S}_{h}$, solution of

$\left\\{\begin{aligned} &S_{h}^{n+1}=\text{arg
}\min_{S_{h}\in{\mathcal{V}_{h}}}G(S_{h}),\text{ where
}G(S_{h}):=\dfrac{1}{2}\|S_{h}-\widehat{S}_{h}^{n+1}\|^{2}_{L^{2}(\Omega)},\\\
&\text{ subject to }\quad c(S_{h}):=S_{h}-1\leqslant 0.\\\
\end{aligned}\right.$ ($\mathbf{P}$)

Our approach is based on the techniques we developed for the Cahn-Hilliard
equation [41], consisting in solving an optimization problem after the
variational problem to force the discrete solution $S_{h}$ to satisfy the
desired physical property $S\leqslant 1$. The optimization problem (P) can be
solved by using the Uzawa algorithm [2, 22, 75] or a projection technique
[41], and it takes about 3-4 iterations for the convergence. We refer to [41]
for more discussions on these techniques, their implementation, and
computational cost. A finite element approximation for ($D_{1s}$) is: Given a
suitable approximation of the initial solution $S_{h}^{0}=\Pi_{h}S_{0}$ and a
proper initialization for $S_{h}^{1}\in\mathcal{S}_{h}$, find
$\widehat{S}_{h}^{n+1}\in\mathcal{V}_{h}$ such that

$\displaystyle\int_{\Omega}\phi_{h}\dfrac{3\widehat{S}_{h}^{n+1}-4S_{h}^{n}+S_{h}^{n-1}}{2\Delta
t}v_{h}\,d\bm{x}+\int_{\Omega}\Big{(}\big{[}K_{s}\big{(}2K_{r}(\psi(S_{h}^{n}))-K_{r}(\psi(S_{h}^{n-1}))\big{)}\big{]}$
(4.7)
$\displaystyle\quad\big{[}\big{(}J(S_{h}^{n})+J_{\delta}^{\prime}(S_{h}^{n})(\widehat{S}_{h}^{n+1}-S_{h}^{n})\big{)}\nabla
h_{cap}+\big{(}2J_{\delta}^{\prime}(S_{h}^{n})-J_{\delta}^{\prime}(S_{h}^{n-1})\big{)}h_{cap}\nabla\widehat{S}_{h}^{n+1}+\nabla
z\big{]}\Big{)}\cdot\nabla v_{h}\,d\bm{x}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\int_{\Gamma}q_{\Gamma}v_{h}\,d\bm{x}=0,\quad\forall
v_{h}\in{\mathcal{V}}_{h},$

Then, find $S_{h}^{n+1}\in\mathcal{S}_{h}$, solution of $(\mathbf{P})$. In the
following, the schemes (4.6) and $(\mathbf{P})$, and (4.7) and $(\mathbf{P})$,
derived from the formulation (3.5) that uses $S$ as the sole unknown, will be
referred to as the implicit $S$-scheme and the semi-implicit $S$-scheme,
respectively.

The weak formulation (3.7)-(3.8) naturally leads to mixed finite element
methods. These methods became popular for modeling flow in porous media due to
their capability of handling general irregular grids and allowing simultaneous
approximation of water content (or pressure head or saturation) and flux as
part of the formulation [9, 12, 13, 58, 61]. The lowest order Raviart-Thomas
elements $(RT_{0})$ [16] is used in these studies to approximate the flux. We
consider the mixed finite element method for ($D_{2i}$): Given a suitable
approximation of the initial solution $S_{h}^{0}=\Pi_{h}S_{0}$ and a proper
initialization for $S_{h}^{1}\in\mathcal{S}_{h}$, find
$(\widehat{S}_{h}^{n+1},\bm{q}_{h}^{n+1})\in{\mathcal{V}}_{h}\times\mathcal{R}_{h}$
such that

$\displaystyle\int_{\Omega}\phi_{h}\dfrac{3\widehat{S}_{h}^{n+1}-4S_{h}^{n}+S_{h}^{n-1}}{2\Delta
t}v_{h}\,d\bm{x}-\int_{\Omega}\bm{q}_{h}^{n+1}\cdot\nabla
v_{h}\,d\bm{x}+\int_{\Gamma}q_{\Gamma}v_{h}\,d\bm{x}=0,\quad\forall
v_{h}\in\mathcal{V}_{h},$ (4.8)
$\displaystyle\int_{\Omega}\bm{q}_{h}^{n+1}\bm{w}_{h}\,d\bm{x}+\int_{\Omega}K_{s}K_{r}(\psi(\widehat{S}_{h}^{n+1}))\big{[}J(\widehat{S}_{h}^{n+1})\nabla
h_{cap}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+h_{cap}J_{\delta}^{\prime}(\widehat{S}_{h}^{n+1})\nabla\widehat{S}_{h}^{n+1}+\nabla
z\big{]}\cdot\bm{w}_{h}\,d\bm{x}=0,\quad\forall\bm{w}_{h}\in\mathcal{R}_{h}.$

Then, find $S_{h}^{n+1}\in\mathcal{S}_{h}$ solution of $(\mathbf{P})$.
A mixed finite element approximation for ($D_{2s}$) is: Given a suitable
approximation of the initial solution $S_{h}^{0}=\Pi_{h}S_{0}$ and a proper
initialization for $S_{h}^{1}\in\mathcal{S}_{h}$, find
$(\widehat{S}_{h}^{n+1},\bm{q}_{h}^{n+1})\in\mathcal{V}_{h}\times\mathcal{R}_{h}$
such that

$\displaystyle\int_{\Omega}\phi_{h}\dfrac{3\widehat{S}_{h}^{n+1}-4S_{h}^{n}+S_{h}^{n-1}}{2\Delta
t}v_{h}\,d\bm{x}-\int_{\Omega}\bm{q}_{h}^{n+1}\cdot\nabla
v_{h}\,d\bm{x}+\int_{\Gamma}q_{\Gamma}v_{h}\,d\bm{x}=0,\quad\forall
v_{h}\in\mathcal{V}_{h},$ (4.9)
$\displaystyle\int_{\Omega}\bm{q}_{h}^{n+1}\cdot\bm{w}_{h}\,d\bm{x}+\int_{\Omega}\Big{(}\big{[}K_{s}\big{(}2K_{r}(\psi(S_{h}^{n}))-K_{r}(\psi(S_{h}^{n-1}))\big{)}\big{]}\big{[}\big{(}J(S_{h}^{n})+J_{\delta}^{\prime}(S_{h}^{n})(\widehat{S}_{h}^{n+1}-S_{h}^{n})\big{)}\nabla
h_{cap}$
$\displaystyle\qquad+(2J_{\delta}^{\prime}(S_{h}^{n})-J_{\delta}^{\prime}(S_{h}^{n-1}))h_{cap}\nabla\widehat{S}_{h}^{n+1}+\nabla
z\big{]}\Big{)}\cdot\bm{w}_{h}\,d\bm{x}=0,\quad\forall\bm{w}_{h}\in\mathcal{R}_{h}.$

Then, find $S_{h}^{n+1}\in\mathcal{S}_{h}$, solution of $(\textbf{P})$. The
schemes (4.8) and $(\textbf{P})$, and (4.9) and $(\textbf{P})$, derived from
the mixed formulation (3.7) dealing with $S$ and $\bm{q}$ as unknowns, will be
referred to as the implicit $(S,\bm{q})$-scheme and the semi-implicit
$(S,\bm{q})$-scheme, respectively.

The weak formulation (3.10) also leads to mixed finite element methods. But
unlike the previous mixed formulation, this one is based on the saturation and
the pressure head. This mixed formulation is simpler and more advantageous for
the spatial discretization since standard Lagrange linear elements are used
instead of the Raviart-Thomas finite elements or the Brezzi-Douglas-Marini
finite elements which lead to larger systems to solve, and are therefore more
demanding in terms of computational cost. For instance, even with the lowest
order Raviart-Thomas elements, the number of degrees of freedom for $\bm{q}$
is equal to the number of edges in the mesh, while with linear Lagrange
element the number of unknowns for $\psi$ is equal to the number of triangle
vertices and there are roughly half many vertices as edges in a 2-D mesh (if
not less). We propose a mixed finite element method for ($D_{3i}$): Given a
suitable approximation of the initial solution $S_{h}^{0}=\Pi_{h}S_{0}$ and a
proper initialization for $S_{h}^{1}\in\mathcal{S}_{h}$, find
$(\widehat{S}_{h}^{n+1},\psi_{h}^{n+1})\in\mathcal{V}_{h}\times\mathcal{V}_{h}$
such that

$\displaystyle\int_{\Omega}\phi_{h}\dfrac{3\widehat{S}_{h}^{n+1}-4S_{h}^{n}+S_{h}^{n-1}}{2\Delta
t}v_{h}\,d\bm{x}+\int_{\Omega}\big{[}K_{s}K_{r}(\psi_{h}^{n+1})(\nabla\psi_{h}^{n+1}+\nabla
z)\big{]}\cdot\nabla v_{h}\,d\bm{x}$ (4.10)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\int_{\Gamma}q_{\Gamma}v_{h}\,d\bm{x}=0,\quad\forall
v_{h}\in{\mathcal{V}}_{h},$
$\displaystyle\int_{\Omega}\psi_{h}^{n+1}w_{h}\,d\bm{x}-\int_{\Omega}h_{cap}J(\widehat{S}^{n+1})w_{h}\,d\bm{x}=0,\quad\forall
w_{h}\in\mathcal{V}_{h}.$

Then, find $S_{h}^{n+1}\in\mathcal{S}_{h}$, solution of $(\textbf{P})$.
A mixed finite element approximation for ($D_{3s}$) is: Given a suitable
approximation of the initial solution $S_{h}^{0}=\Pi_{h}S_{0}$ and a proper
initialization for $S_{h}^{1}\in\mathcal{S}_{h}$, find
$(\widehat{S}_{h}^{n+1},\psi_{h}^{n+1})\in\mathcal{V}_{h}\times\mathcal{V}_{h}$
such that

$\displaystyle\int_{\Omega}\phi_{h}\dfrac{3\widehat{S}_{h}^{n+1}-4S_{h}^{n}+S_{h}^{n-1}}{2\Delta
t}v_{h}\,d\bm{x}+\int_{\Omega}\big{[}K_{s}\big{(}2K_{r}(\psi_{h}^{n})-K_{r}(\psi_{h}^{n-1})\big{)}\big{(}\nabla\psi_{h}^{n+1}+\nabla
z\big{)}\big{]}\cdot\nabla v_{h}\,d\bm{x}$ (4.11)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\int_{\Gamma}q_{\Gamma}v_{h}\,d\bm{x}=0,\quad\forall
v_{h}\in{\mathcal{V}}_{h},$
$\displaystyle\int_{\Omega}\psi_{h}^{n+1}w_{h}\,d\bm{x}-\int_{\Omega}h_{cap}\big{(}J(S_{h}^{n})+J_{\delta}^{\prime}(S_{h}^{n})(\widehat{S}_{h}^{n+1}-S_{h}^{n})\big{)}w_{h}\,d\bm{x}=0,\quad\forall
w_{h}\in\mathcal{V}_{h}.$

Then, find $S_{h}^{n+1}\in\mathcal{S}_{h}$, solution of $(\textbf{P})$. To the
best of our knowledge, the linear second-order time stepping mixed finite
element method (4.11) for the Richards equation (2.3) has not been explored
before. The schemes (4.10) and $(\mathbf{P})$, and (4.11) and $(\mathbf{P})$
will be referred to as the implicit $(S,\psi)$-scheme and the semi-implicit
$(S,\psi)$-scheme, respectively.

######  Remark 4.1.

All the numerical algorithms are two steps schemes and therefore, they require
the use of a starting procedure to obtain an approximation of the solution at
the first step.

######  Remark 4.2.

We will use the Newton method for the linearization of the implicit
$(S,\psi)$-scheme. The use of this method for the implicit $S$\- and
$(S,\bm{q})$-schemes requires the computation of the second derivative of the
Leverett $J$-function. The resulting $J$-functions from many of the
constitutive relationships, in particular equations (5.4) and (5.15) used in
this study, are not sufficiently differentiable (e.g., $J^{\prime\prime}$ is
not defined at $S=0$ and $S=1$). This leads to convergence issues when the
Newton method is applied. Therefore, we will use the Picard method for the
linearization of the implicit $S$\- and $(S,\bm{q})$-schemes. The stopping
criteria for the iterative methods is
$\|u_{h}^{n,k+1}-u_{h}^{n,k}\|_{L^{2}(\Omega)}\leqslant\epsilon$, where
$0<\epsilon\ll 1$ is the chosen tolerance, $u_{h}^{n,k+1}$ and $u_{h}^{n,k}$
are the solutions at the iterations $k+1$ and $k$, respectively.

######  Remark 4.3.

For each of the above schemes, one can easily show the conservation of the
variable $\widehat{S}_{h}$ over time by taking $v_{h}$ to be constant on the
domain $\Omega$ and considering no-flux boundary conditions, i.e., setting
$q_{\Gamma}=0$. However, solving directly the additional optimization problem
(P) by an Uzawa algorithm or a simple projection technique could lead to a
defect of conservation for the variable $S_{h}$. We have quantified this error
in a previous work on the Cahn-Hilliard equation (see Figure 1 in [41]). We
also proposed in [41] new conservative projection techniques for solving the
optimization problem while maintaining the conservation. These techniques can
be easily applied to the Richards equation to satisfy the conservation. We
also compare the time evolution of the total mass of the exact and numerical
solutions in the test case in section 5.2. Our tests show a good agreement for
the total mass.

######  Remark 4.4.

The discrete time-derivative introduces a “distributed” mass matrix (the term
in $\hat{S}_{h}^{n+1}$) in the global algebraic system for all our finite
element methods. As noted in [21, 39], this distributed mass matrix (without
lumping) leads to a numerical method that no longer satisfies a discrete
maximum principle when the time step is taken smaller than a prescribed
threshold value. Without discrete maximum principle, the solution may suddenly
become non-monotone with unphysical oscillations from one spatial node to the
next. This loss of discrete maximum principle for time steps below a threshold
was proven for linear parabolic equations, such as the heat equation [62]. For
Richards equation, the same threshold on the time step holds in numerical
tests for the loss of monotonicity of the numerical solutions, though no proof
is available even for simple time-stepping schemes such as backward Euler
method [39]. A common way employed in finite element methods for the Richards
equation to recover a discrete maximum principle is the mass lumping
technique. An early study on the subject [21] illustrated that mass lumping
can effectively prevent nonoscillatory solutions, by comparing solutions with
and without mass lumping on the time derivative. Several forms of mass lumping
have been proposed in the standard and mixed finite element methods for
unsaturated flow problems [12, 24, 51, 78]. A simple mass lumping technique,
based on the row-sum of all the terms in the time derivative, is used in our
study and was effective at removing oscillations that may occur near the
drainage front. However, proving a discrete maximum principle may turn out to
be challenging for our two-step methods.

##  5 Numerical tests

A series of numerical problems is given in this section to test and compare
the performance of the numerical methods discussed in section 4. All the test
cases are in two dimensions. We use the sparse linear solver UMFPACK which
provides a relatively efficient solution. Computations and results could thus
carry over to 3D in a straightforward manner. However, this would require
modifications in the implementation of the schemes. All the algorithms for the
test cases are implemented using the FreeFem++ software [37, 38].

###  5.1 Manufactured solutions

In this subsection, we numerically investigate the rates of convergence in
space and time of the different schemes. Following the framework presented in
[7, 70], we employ a manufactured solution to test the methods. The idea is to
use an arbitrary sufficiently differentiable function as exact solution with a
source term for the computation of numerical errors. We perform numerical
simulations using the computational domain
$\Omega=[0\,cm,4\,cm]\times[0\,cm,20\,cm]$ and the final computational time
$T=120\,s$. We consider the manufactured solution given by

$\psi(x,z,t)=20.4\tanh\big{(}0.5(z+\frac{t}{12}-15)\big{)}+c,\quad
c\in\mathbb{R},$ (5.1)

where $c$ is a constant satisfying $c\leqslant-20.4$ to ensure the non-
positivity of $\psi$. The Richards equation (2.3) is then solved with
Dirichlet boundary conditions and a source term such that (5.1) is a solution,
i.e., the closed-form expression (5.1) is used to calculate the adequate
source term as well as the initial and boundary conditions. We use the
following relationships between the relative permeability, pressure head and
water content [36]:

$\theta=\theta_{r}+\dfrac{\theta_{s}-\theta_{r}}{1+|\tilde{\alpha}\psi|^{\beta}},\quad
K_{r}=\dfrac{1}{1+|\tilde{A}\psi|^{\gamma}},$ (5.2)

where

$\displaystyle\theta_{s}=0.287,\qquad\tilde{\alpha}=0.0271\,cm^{-1},\qquad
K_{s}=9.44\times 10^{-3}\,cm.s^{-1},\qquad\gamma=4.74,$ (5.3)
$\displaystyle\theta_{r}=0.075,\qquad\beta=3.96,\qquad\qquad\qquad\tilde{A}=0.0524\,cm^{-1}.$

Using (5.2), the pressure head $\psi$ can be expressed as in (3.1), where

$h_{cap}(\bm{x})=\dfrac{1}{\tilde{\alpha}}\quad\text{and}\quad
J(S)=-(S^{-1}-1)^{1/\beta}.$ (5.4)

The saturation is given by the relation

$S=\dfrac{1}{1+|\tilde{\alpha}\psi|^{\beta}},$ (5.5)

that one uses to calculate the initial and boundary conditions for $S$ from
the conditions on $\psi$.

####  5.1.1 Numerical example without regularization of the Leverett
$J$-function

For the first test case, we take $c=-41.1$ in (5.1). The profile of the
manufactured solutions is shown in Figure 1.

|
---|---
Figure 1: Test case without regularization of the Leverett $J$-function:
Profile of the manufactured solution (5.1) at different times for $c=-41.1$.

The saturation $S$ is far from $1$ for all $t\in[0,T]$. Therefore,
$J_{\delta}^{\prime}(S)$ in (3.2) is exactly given by $J^{\prime}(S)$. Note
also that no specific treatment is needed in the numerical discretization of
the source term. We use a constant time step $\Delta t=0.2$ for all the
schemes. The tolerance for the stopping criteria for the iterative schemes is
$\epsilon=10^{-5}$. For each scheme, the $L^{2}$\- and $H^{1}$-error on the
effective saturation $S_{h}$ and the pressure head $\psi_{h}$ are computed on
different meshes. The numerical results are presented in Figure 2. We can see
that all the schemes have the same order of convergence around $k=1$ for the
$H^{1}$-error on the variables $S_{h}$ and $\psi_{h}$. For the $L^{2}$-error,
an order of convergence varying around $k=2$ is obtained on both variables
$S_{h}$ and $\psi_{h}$ with each scheme. This order is slightly larger with
the semi-implicit schemes. There is no significant difference on the errors
produced by the different schemes. However, the error on $\psi_{h}$ is larger
than the one on $S_{h}$. The computational times for each scheme with respect
to the mesh refinement are plotted in Figure 3. Here, the semi-implicit
schemes are the winner in terms of computational cost since we use the same
time step for all the schemes. Larger time step could be used for the implicit
schemes as we will see in the next test cases. With respect to the
formulations, our numerical simulations show that the $(S,\psi)$-schemes are
faster than the $S$-schemes which are faster than the $(S,\bm{q})$-schemes.
This shows the efficiency of the mixed formulation based on the saturation and
the pressure head with respect to the other formulations.

|
---|---
|
Figure 2: Manufactured solution — Test case without regularization of the
Leverett $J$-function: $L^{2}$\- and $H^{1}$-error on $S_{h}$ and $\psi_{h}$
as a function of the mesh size $h$ for each scheme with respect to the
manufactured solution (5.1) for $c=-41.1$. Figure 3: Manufactured solution —
Test case without regularization of the Leverett $J$-function: Computational
time in seconde (s) as a function of the mesh size $h$ for each scheme.

For the error and convergence analysis in time, we compute a reference
solution $(S_{h}^{\text{ref}}$,$\psi_{h}^{\text{ref}})$ with the semi-implicit
$(S,\psi)$-scheme on a computational grid of size $h=0.176$ ($32\times 160$
elements) by using a small time step $\Delta t=0.001\,s$ to get a good
approximation of the semi-discretized solution of the problem. We then analyze
the numerical convergence in time by varying the time step. For each time
step, the $L^{2}$-error on $S_{h}$ and $\psi_{h}$ with respect to the
reference solution are recorded at the final time $T=120\,s$. Figure 4 shows
the errors as function of the time step for each scheme. An order of
convergence in time varying slightly around $p=2$ is obtained on the variables
$S_{h}$ and $\psi_{h}$ with each scheme. The schemes produce almost the same
error for each variable. However, the error on $\psi_{h}$ is larger than the
error on $S_{h}$.

|
---|---
Semi-implicit schemes | Implicit schemes
Figure 4: Manufactured solution — Test case without regularization of the
Leverett $J$-function: $L^{2}$-error on $S_{h}$ and $\psi_{h}$ as a function
of the time step $\Delta t$ for each scheme.

Most numerical methods for the Richards equation (2.1) are based on the
classical backward Euler method [3, 4, 8, 9, 12, 29, 32, 53, 60, 61, 67]. For
the purpose of comparison of this standard technique with the proposed
schemes, we consider the finite element methods for the formulations (3.4) and
(3.6), respectively, using the classical backward Euler method for the
temporal discretization. These schemes are named backward Euler $S$-scheme for
the formulation (3.4) and backward Euler $(S,\bm{q})$-scheme for the
formulation (3.6), and are compared against the closest among the proposed
schemes, namely the $S$-schemes. Table 1 shows the $L^{2}$-error on $S_{h}$
and $\psi_{h}$, and the computational times for the solutions computed with
the $S$-schemes, backward Euler $S$-scheme and backward Euler
$(S,\bm{q})$-scheme, respectively, on the mesh $32\times 160$ elements with
the time step $\Delta t=4\,s$. Our results show that the errors with classical
methods are at least seven times larger than the errors with the proposed
implicit $S$-scheme. The proposed semi-implicit $S$-scheme produces almost the
same errors as the standard methods but with a computational time at least six
times smaller. By reducing the time step, the semi-implicit $S$-scheme becomes
more accurate and efficient than the standard methods. This confirms the
efficiency of the proposed techniques with respect to the standard methods.

Methods |  $L^{2}$-error on $S_{h}$ |  $L^{2}$-error on $\psi_{h}$ | CPU(s)
---|---|---|---
Implicit $S$-scheme | $0.0075831$ | $0.53518$ | $237$
Semi-implicit $S$-scheme | $0.0775579$ | $3.37313$ | $040$
Backward Euler $S$-scheme | $0.0603825$ | $3.84694$ | $262$
Backward Euler $(S,\bm{q})$-scheme | $0.0595938$ | $3.80333$ | $496$
Table 1: Manufactured solution — Test case without regularization of the
Leverett $J$-function: Comparison of errors and computational times between
standard and proposed methods.

####  5.1.2 Numerical example with regularization of the Leverett $J$-function

For the second test case, we set $c=-20.4$ in (5.1). The profile of the
manufactured solution is shown in Figure 5.

|
---|---
Figure 5: Test case with regularization of the Leverett $J$-function: Profile
of the manufactured solution (5.1) at different times for $c=-20.4$.

The purpose of this test case is to verify the convergence behavior of the
different schemes with the regularization (3.2) when the saturation solution
$S$ reaches the value of one in some parts of the domain. We have noticed in
practice that the $(S,\psi)$-schemes allow smaller values of $\delta$ than the
$S$\- and $(S,\bm{q})$-schemes. With respect to the temporal discretization,
the implicit $S$\- and $(S,\bm{q})$-schemes allow smaller values of $\delta$
than their corresponding semi-implicit versions. We take $\delta=10^{-3}$ for
the semi-implicit $S$\- and $(S,\bm{q})$-schemes, $\delta=10^{-6}$ for the
implicit $S$\- and $(S,\bm{q})$-schemes, and $\delta=10^{-10}$ for the
$(S,\psi)$-schemes. The tolerance for the iterative methods used for the
linearization of the implicit schemes is set to $5\times 10^{-4}$. For each
scheme, the $L^{2}$\- and $H^{1}$-error on $S_{h}$ and $\psi_{h}$ on different
meshes are presented in Figure 6 and the computational time for each scheme
with respect to the mesh refinement is plotted in Figure 7.

|
---|---
|
Figure 6: Manufactured solution — Test case with regularization of the
Leverett $J$-function: $L^{2}$\- and $H^{1}$-error on $S_{h}$ and $\psi_{h}$
as a function of the mesh size $h$ for each scheme with respect to the
manufactured solution (5.1) for $c=-20.4$. Figure 7: Manufactured solution —
Test case with regularization of the Leverett $J$-function: Computational time
in seconde (s) as a function of the mesh size $h$ for each scheme.

Numerical results show that the $(S,\psi)$-schemes produce almost the same
error for a fixed mesh size, and converge with an order around $k=1$ for the
$H^{1}$-error and an order about $k=2$ for the $L^{2}$-error on both variables
$S_{h}$ and $\psi_{h}$. The implicit $S$\- and $(S,\bm{q})$-schemes converge
slowly with an order less than $1$ on fine meshes, while a lack of convergence
is noticed with the semi-implicit $S$\- and $(S,\bm{q})$-schemes. This shows
again that the linear scheme using the mixed formulation consisting in
simultaneously constructing approximations of the saturation and the pressure
head with Lagrange linear element performs better in terms of accuracy and
efficiency than the schemes using the mixed formulation based on the
saturation and the flux using Raviart-Thomas elements. The implicit schemes
perform better than the semi-implicit schemes in terms of efficiency. However,
the iterative methods used for the linearization of the implicit schemes could
not converge for some less regular solution as we will see with the test case
in section 5.3. In such cases, semi-implicit schemes could be an alternative.

We perform numerical tests to obtain the temporal order of convergence of the
schemes. To this end, we compute a reference solution
$(S_{h}^{\text{ref}}$,$\psi_{h}^{\text{ref}})$ with the implicit
$(S,\psi)$-scheme using a computational grid of size $h=0.088$ ($64\times 320$
elements) and a time step $\Delta t=0.1\,s$. We then analyze the numerical
convergence in time by varying the time step for each scheme. For each time
step, the $L^{2}$-error on $S_{h}$ and $\psi_{h}$ with respect to the
reference solution is recorded at the final time $T=120\,s$. Figure 8 shows
the $L^{2}$-error on $S_{h}$ and $\psi_{h}$ as a function of the time step for
each scheme. An order of convergence in time around $p=2$ is obtained with the
semi-implicit $(S,\psi)$-scheme and an order around $p=1.8$ is obtained with
the implicit $(S,\psi)$-scheme on the variables $S_{h}$ and $\psi_{h}$. As for
the spatial convergence, a lack of convergence in time is noticed with the
$S$\- and $(S,\bm{q})$-schemes.

A second-order time-accuracy and optimal spatial convergence rates are
obtained with the $(S,\psi)$-schemes for solutions with regularization,
attesting the robustness of the $(S,\psi)$-schemes for solutions in which the
saturation reaches the value of one in some parts of the domain.

|
---|---
Semi-implicit schemes | Implicit schemes
Figure 8: Manufactured solution — Test case with regularization of the
Leverett $J$-function: $L^{2}$-error on $S_{h}$ and $\psi_{h}$ as a function
of the time step $\Delta t$ for each scheme.

###  5.2 2-D analytical solution of the Green-Ampt problem

In this section, we use an analytical solution to (2.3) for checking the
accuracy of our numerical schemes. Closed-form one-dimensional solutions for
the Richards equation are available [17, 34, 56, 65, 66]. Two- and three-
dimensional analytical solutions have been derived in [72, 73]. Here, we use a
two-dimensional analytical solution proposed in [73].

Figure 9(i) shows a 2-D soil sample of dimensions $a\times L$.

LaRainfall(i)

(ii)

Figure 9: (i) Two-dimensional soil sample $[0,a]\times[0,L]$. (ii) Analytical
solution for effective saturation $S$ at time $t=10$ days.

The column of soil is initially dry ($\psi=\psi_{d}$) and the infiltration is
considered at the top boundary. This problem is referred to as the Green-Ampt
problem in [73]. An analytical solution to this problem using Gardner
formulation [35] for the relative hydraulic conductivity $K_{r}$ and a simple
formulation from [77] for the water content $\theta$,

$K_{r}=e^{(\alpha\psi)},\qquad\theta=\theta_{r}+(\theta_{s}-\theta_{r})e^{(\alpha\psi)},$
(5.6)

is derived in [73], by imposing a specific pressure head at the top that is
zero in the middle and tapers rapidly to $\psi_{d}$ at $x=0$ and $x=a$, and is
maintained at $\psi=\psi_{d}$ along the bottom and vertical sides of the soil.
The initial and boundary conditions are [73]:

$\displaystyle\psi(x,z,0)=\psi_{d},\quad(x,z)\in[0,a]\times[0,L],$ (5.7)
$\displaystyle\psi(x,0,t)=\psi(a,z,t)=\psi(0,z,t)=\psi_{d},\quad
x\in[0,a],\text{ }z\in[0,L],\text{ }t\geqslant 0,$
$\displaystyle\psi(x,L,t)=\dfrac{1}{\alpha}\log\Big{(}\varepsilon+\big{(}1-\varepsilon\big{)}\Big{[}\frac{3}{4}\sin\big{(}\frac{\pi
x}{a}\big{)}-\frac{1}{4}\sin\big{(}\frac{3\pi x}{a}\big{)}\Big{]}\Big{)},\quad
x\in[0,a],\text{ }z\in[0,L],\text{ }t\geqslant 0$

and the analytical solution of this problem is given by [73]

$\psi(x,z,t)=\dfrac{1}{\alpha}\log(\varepsilon+\tilde{\psi}(x,z,t)),$ (5.8)

where

$\displaystyle\tilde{\psi}=(1-\varepsilon)\exp\big{(}\frac{\alpha}{2}(L-z)\big{)}\Big{\\{}\frac{3}{4}\sin\big{(}\frac{\pi
x}{a}\big{)}\Big{[}\dfrac{\sinh(\beta_{1}z)}{\sinh(\beta_{1}L)}+\frac{2}{Lb}\sum_{k=1}^{\infty}(-1)^{k}\frac{\lambda_{k}}{\gamma_{1k}}\sin(\lambda_{k}z)\exp(-\gamma_{1k}t)\Big{]}$
(5.9)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-\frac{1}{4}\sin\big{(}\frac{3\pi
x}{a}\big{)}\Big{[}\dfrac{\sinh(\beta_{3}z)}{\sinh(\beta_{3}L)}+\frac{2}{Lb}\sum_{k=1}^{\infty}(-1)^{k}\frac{\lambda_{k}}{\gamma_{3k}}\sin(\lambda_{k}z)\exp(-\gamma_{3k}t)\Big{]}\Big{\\}}$

with

$\lambda_{k}=\dfrac{k\pi}{L},\text{ }\varepsilon=e^{(\alpha\psi_{d})},\text{
}b=\dfrac{\alpha(\theta_{s}-\theta_{r})}{K_{s}},\text{
}\tilde{\lambda}_{i}=\dfrac{i\pi}{a},\text{
}\beta_{i}=\sqrt{\frac{\alpha^{2}}{4}+\tilde{\lambda}_{i}^{2}},\text{
}\gamma_{ik}=\dfrac{\beta_{i}^{2}+\lambda_{k}^{2}}{b}.$ (5.10)

Using (5.6), the pressure head $\psi$ can be expressed in the form (3.1),
where

$\quad h_{cap}(\bm{x})=\dfrac{1}{\alpha},\quad J(S)=\log(S).$ (5.11)

In practice, one should truncate the series in (5.9) to calculate the exact
solution. Here, we consider the first $200$ terms of the series.

We solve the equation (2.3) with the semi-implicit $(S,\psi)$-scheme in the
domain represented in Figure 9(i) by using the following values for the
parameters

$L=a=50m,\text{ }K_{s}=0.2m/day,\text{ }\psi_{d}=-50m,\text{
}\theta_{s}=0.45,\text{ }\theta_{r}=0.15,\text{ }\alpha=0.1m^{-1}.$ (5.12)

The analytical solution for the effective saturation at time $t=10$ days is
represented in Figure 9(ii). The $L^{2}$\- and $H^{1}$-error norms on the
effective saturation $S_{h}$ and the pressure head $\psi_{h}$ are shown in
Table 2. The contours of the analytical and numerical solutions for the
effective saturation at time $t=10$ days are plotted in Figure 10. The results
show that the numerical solution converges to the analytical solution as the
mesh is refined. An order of convergence varying around $k=1.8$ for the
$L^{2}$-error and an order varying around $k=0.9$ for the $H^{1}$-error are
obtained on both variables $S_{h}$ and $\psi_{h}$. However, the error on the
pressure head is larger.

|  |
---|---|---
Mesh $25\times 25$ | Mesh $50\times 50$ | Mesh $100\times 100$
Figure 10: Contours of the analytical (black curves) and numerical (red curves) solutions on different meshes at time $t=10$ days. |  $L^{2}$-error |  $H^{1}$-error
---|---|---
$\Delta t$ | Mesh | $S_{h}$ | $\psi_{h}$ | $S_{h}$ | $\psi_{h}$
$0.010$ | $25\times 25$ | $0.055429$ | $26.3803$ | $0.125187$ | $41.3671$
$0.005$ | $50\times 50$ | $0.016745$ | $8.72881$ | $0.057976$ | $22.2810$
$0.0025$ | $100\times 100$ | $0.004397$ | $2.45371$ | $0.027922$ | $11.9616$
$0.00125$ | $200\times 200$ | $0.001182$ | $0.54719$ | $0.013805$ | $6.20522$
Table 2: Analytical solution: $L^{2}$\- and $H^{1}$-error on $S_{h}$ and
$\psi_{h}$ for the solution computed with the semi-implicit $(S,\psi)$-scheme
on different meshes.

To test the conservation property of the scheme, we compare the total mass in
the computational domain of the exact and numerical solutions as time evolves.
The time evolution of the total mass of the exact and numerical solutions for
the saturation and the relative error are displayed in Figure 11. The results
show a good agreement for the total mass of the exact and numerical solutions.

|
---|---
(a) | (b)
Figure 11:  (a) Time-evolution of the total mass in the computational domain
of the exact and numerical solutions and (b) relative error on the mesh
$100\times 100$ with $\Delta t=0.0025$ day.

###  5.3 Infiltration in heterogeneous medium

The objective of the next test cases is to investigate the performance of the
algorithms for infiltration in heterogeneous medium.

####  5.3.1 Curvilinearly layered soil

This test case was used by Manzini et al. [48], where two layers of soil are
considered. The capillary pressure is modeled by the van Genuchten
constitutive relationship [76]:

$\theta=\dfrac{\theta_{s}-\theta_{r}}{[1+(\alpha|\psi|)^{n}]^{m}}+\theta_{r}$
(5.13)

and the relative hydraulic conductivity is expressed as follows [52, 76]:

$K_{r}=\dfrac{[1-(\alpha|\psi|)^{n-1}(1+(\alpha|\psi|)^{n})^{-m}]^{2}}{[1+(\alpha|\psi|)^{n}]^{m/2}},\quad
m=1-1/n.$ (5.14)

Using (5.13), the pressure head $\psi$ can be expressed in the form (3.1),
where

$\quad h_{cap}(\bm{x})=\dfrac{1}{\alpha},\quad J(S)=-(S^{-1/m}-1)^{1/n}$
(5.15)

and the saturation can then be expressed as

$S=[1+(\alpha|\psi|)^{n}]^{-m}.$ (5.16)

The computational domain (see Figure 12) is the square
$[0\,cm,100\,cm]\times[0\,cm,100\,cm]$.

$z\geqslant\xi(x)$$z\leqslant\xi(x)$0100100 Figure 12: Geometry of two layers
of soil separated by a curved interface.

The two different regions of the layered soil are separated by the curved
interface

$\xi(x)=100\big{(}0.1(1-\cos(\pi x/100))+0.45\big{)}.$ (5.17)

The parameter values for (5.13)-(5.14) are given by

$\displaystyle
z\geqslant\xi(x):\quad\theta_{s}=0.50,\quad\theta_{r}=0.120,\quad\alpha=0.028\,cm^{-1},\quad
n=3.00,\quad K_{s}=0.25\,cm.h^{-1},$ (5.18) $\displaystyle
z\leqslant\xi(x):\quad\theta_{s}=0.46,\quad\theta_{r}=0.034,\quad\alpha=0.016\,cm^{-1},\quad
n=1.37,\quad K_{s}=2.00\,cm.h^{-1}.$

Homogeneous Dirichlet boundary conditions are set on the top and bottom sides
of the domain, while homogeneous Neumann conditions are set on the two
vertical sides. We consider the initial condition $\psi(x,z,t=0)=-z$, the mesh
size $h=1\,cm$ and the time step $\Delta t=\frac{1}{60}$ hour. The Newton
method applied to the implicit $(S,\psi)$-scheme does not converge with the
parameter values (5.18) because the relative permeability function $K_{r}$
expressed by (5.14) is not differentiable at point $\psi=0$ for $n<2$. For
physical values of $n$ with $0<n<2$, the implicit $(S,\psi)$-scheme may not
converge. In this case, we use the semi-implicit $(S,\psi)$-scheme for the
numerical simulations. The time-evolution of the water content is represented
in Figure 13. Our numerical results agree very well with the ones presented by
Manzini et al. [48]. This shows that linear semi-implicit $(S,\psi)$-scheme is
a good alternative when implicit schemes encounter converging difficulties.

|  |  |
---|---|---|---
$t=0$ |  $t=0.25$ day |  $t=0.5$ day |  $t=0.75$ day
|  |  |
$t=1$ day |  $t=1.25$ day |  $t=1.5$ day |  $t=1.75$ day
Figure 13: Curvilinearly layered soil — Test case 1: Time-evolution of the
water content $\theta_{h}$ computed with the semi-implicit $(S,\psi)$-scheme.

To test the performance of our scheme for large degree of heterogeneity, we
reconduct the test case using large difference between the saturated hydraulic
conductivities of the layers and we keep the same values for the other
parameters. We perform computations for $K_{s}=2.5\,cm.h^{-1}$ and
$K_{s}=25\,cm.h^{-1}$ for $z\leqslant\xi(x)$, respectively, using different
time steps. Numerical results (not shown here) show that the semi-implicit
$(S,\psi)$-scheme remains stable in modeling unsaturated flow through highly
heterogeneous soils for a wide range of time steps.

To investigate the order of convergence of the semi-implicit $(S,\psi)$-scheme
with this non regular solution, we compute a reference solution on a grid of
mesh size $h=0.28\,cm$ ($500\times 500$ elements) with a time step $\Delta
t=\frac{1}{240}$ hour at time $t=1$ day. The $L^{2}$\- and $H^{1}$-error norms
are computed with respect to this reference solution for different values of
the mesh size to evaluate the spatial order of convergence. Figure 14 presents
the errors as function of the mesh size. An order of convergence around
$k=1.5$ for the $L^{2}$-error and around $k=1.2$ for the $H^{1}$-error are
observed on the pressure head $\psi_{h}$ while an order of convergence around
$k=0.6$ for the $L^{2}$-error and around $k=0.1$ for the $H^{1}$-error are
observed on $S_{h}$. The lack of regularity of this solution deteriorates
significantly the order of convergence on the saturation variable as observed
in [74]. The errors on the pressure head $\psi_{h}$ are larger than the errors
on the saturation $S_{h}$. However, the magnitude of the difference decreases
with the mesh refinement since the orders for $\psi_{h}$ are larger than the
orders for $S_{h}$.

|
---|---
Figure 14: Curvilinearly layered soil — Test case 1: $L^{2}$\- and
$H^{1}$-error on $S_{h}$ and $\psi_{h}$ as a function of the mesh size $h$
computed with the semi-implicit $(S,\psi)$-scheme with respect to the
reference solution.

Another test case with this geometry was used in [6] but for two layers of
soil having the same characteristics except for the saturated hydraulic
conductivity which is $K_{s}$ at top ($z\geqslant\xi(x)$) and $2K_{s}$ at
bottom ($z\leqslant\xi(x)$). The parameter values for capillary pressure
function (5.13) and the relative permeability function (5.14) are

$\theta_{s}=0.5,\quad\theta_{r}=0.12,\quad\alpha=0.02\,cm^{-1},\quad n=3,\quad
K_{s}=6.944\times 10^{-5}cm.s^{-1}.$ (5.19)

The initial condition is $\psi(x,z,t=0)=-z$. A homogeneous Neumann condition
is enforced on the lateral parts of the domain, while a homogeneous Dirichlet
one is imposed on top and bottom sides. We reproduce this test case using the
implicit $(S,\psi)$-scheme. The mesh size is $h=1\,cm$ and the time step
$\Delta t=\frac{1}{36}$ hour. The time-evolution of the relative pressure
$\psi_{h}^{t}-\psi_{h}^{0}$ is represented in Figure 15. A good agreement of
our numerical results with the ones presented in [6] is observed.

|  |
---|---|---
$t=0$ | $t=36h$ | $t=72h$
Figure 15: Curvilinearly layered soil — Test case 2: Relative pressure
$\psi_{h}^{t}-\psi_{h}^{0}$ at different times $t$ computed with the implicit
$(S,\psi)$-scheme.

####  5.3.2 Layers of soil of L-shape form

Here, we perform numerical simulations using the test case presented in [7]
for two layers of soil with same characteristics given by the relations
(5.13)-(5.14) for the capillary pressure and the relative hydraulic
conductivity with the parameter values (5.19), except for the saturated
hydraulic conductivity $K_{s}$ which changes in the medium. The computational
domain is the square $[0\,cm,100\,cm]\times[0\,cm,100\,cm]$ with variable
saturated hydraulic conductivity as shown in Figure 16. The initial condition
is $\psi(x,z,t=0)=-z$. A homogeneous Neumann condition is enforced on the
lateral parts of the domain, while a homogeneous Dirichlet one is imposed on
top and bottom sides. We discretize the domain using the mesh size $h=1\,cm$.
The time step is $\Delta t=\frac{1}{180}$ hour. The time-evolution of the
effective saturation computed with the implicit $(S,\psi)$-scheme is
represented in Figure 17. Good agreement with the results presented in [7] is
observed.

$K_{s}$$3K_{s}$0100100 Figure 16: Geometry of two layers of soil of L-shape form. |  |
---|---|---
$t=12h$ | $t=24h$ | $t=48h$
Figure 17: L-shape geometry test case: Saturation $S_{h}$ at different times
computed with the implicit $(S,\psi)$-scheme.

To evaluate the order of convergence of the implicit $(S,\psi)$-scheme with
this solution, we compute a reference solution on a grid of mesh size
$h=0.28\,cm$ ($500\times 500$ elements) at time $t=24h$. Figure 18 shows the
errors as function of the mesh size. A convergence rate around $k=1.8$ for the
$L^{2}$-error and around $k=1.3$ for the $H^{1}$-error are observed on each
variable. The errors on the pressure head $\psi_{h}$ are about $70$ times
larger than the errors on the saturation $S_{h}$.

|
---|---
Figure 18: L-shape geometry test case: $L^{2}$\- and $H^{1}$-error on $S_{h}$
and $\psi_{h}$ as a function of the mesh size $h$ computed with the implicit
$(S,\psi)$-scheme with respect to the reference solution.

######  Remark 5.1.

All the test cases in section 5.3 have two saturated boundary conditions on
both top and bottom sides. The environments tend to full saturation.
Therefore, they serve also as an additional check of convergence behavior of
the $(S,\psi)$-schemes with the regularization (3.2) when the saturation
reaches the value of one.

##  6 Conclusion

The present study focuses on several iterative and noniterative second-order
time stepping standard/or mixed finite element methods to approximate the
mixed form of the Richards equation written in three different weak
formulations. The first weak formulation uses the saturation as the sole
unknown and is discretized using the standard linear Lagrange elements. The
second weak formulation uses the saturation and the vector flux as unknowns
and is discretized using the linear Lagrange elements for the saturation and
the lowest order Raviart-Thomas elements for the vector flux. The third weak
formulation is based on the saturation and the pressure head and is
discretized using the linear Lagrange element for both variables. The
different methods are studied and compared with respect to accuracy and
computational requirements using exact solutions and commonly used test cases.
Some numerical results are presented involving solutions that require
regularization techniques to avoid degeneracy of the capillary pressure
function when the saturation reaches the value of one in some regions of the
domain.

Our numerical investigations show that the $(S,\psi)$-schemes (derived from
the weak formulation using the saturation and the pressure head as unknowns)
perform better in terms of accuracy and efficiency than the
$(S,\bm{q})$-schemes (derived from the weak formulation using the saturation
and the flux head as unknowns) and the $S$-schemes (derived from the weak
formulation using the saturation as the sole unknown). The results of our
study also indicate that second order accurate time stepping methods are more
efficient than the traditional first order accurate schemes.

The implicit $(S,\psi)$-scheme allows the use of larger time steps. However,
the iterative procedure may not converge because of the nonlinearity of the
soil constitutive laws. The semi-implicit $(S,\psi)$-scheme is linear, as
accurate as the implicit $(S,\psi)$-scheme, easy to implement and requires
less memory. In summary, our results show that the proposed semi-implicit
$(S,\psi)$-scheme is the best among the methods studied and is a good
alternative to traditional iterative methods for the Richards equation.

## Acknowledgments

The first author was supported through a Postdoctoral Fellowship of the
UM6P/OCP group of Morocco and a Postdoctoral Fellowship of the Fields
Institute. The second author acknowledges funding from UM6P/OCP group of
Morocco. The research of the third author is funded through a Discovery Grant
of the Natural Sciences and Engineering Research Council of Canada.

## References

  * [1] G. Akrivis.  Stability of implicit-explicit backward difference formulas for nonlinear parabolic equations.  SIAM J. Numer. Anal., 53:464–484, 2015.
  * [2] G. Allaire.  Numerical Analysis and Optimization.  Oxford University Press, New York, 2007.
  * [3] T. Arbogast, M. Obeyesekere, and M. F. Wheeler.  Numerical methods for the simulation of flow in root-soil systems.  SIAM J. Numer. Anal., 30:1677–1702, 1993.
  * [4] T. Arbogast and M. F. Wheeler.  A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media.  SIAM J. Numer. Anal., 33:1669–1687, 1996.
  * [5] U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton.  Implicit-explicit methods for time-dependent partial differential equations.  SIAM J. Numer. Anal., 32:797–823, 1995.
  * [6] V. Baron.  Méthodes numériques pour les écoulements en milieu poreux : estimations à posteriori et stratégie d’adaptation. PhD thesis, LMJL-Laboratoire de Mathématiques Jean Leray, 2015.
  * [7] V. Baron, Y. Coudière, and P. Sochala.  Adaptive multistep time discretization and linearization based on a posteriori error estimates for the Richards equation.  Appl. Numer. Math., 112:104–125, 2017.
  * [8] M. Bause.  Higher and lowest order mixed finite element approximation of subsurface flow problems with solutions of low regularity.  Adv. Water Resour., 31:370–382, 2008.
  * [9] M. Bause and P. Knabner.  Computation of variably saturated subsurface flow by adaptive mixed hybrid finite element methods.  Adv. Water Resour., 27:565–581, 2004.
  * [10] J. Bear.  Dynamics of Fluids in Porous Media.  American Elsevier Publishing Company, New York, 1972.
  * [11] J. Becker.  A second order backward difference method with variable steps for a parabolic problem.  BIT Numer. Math., 38:644–662, 1998.
  * [12] B. Belfort, F. Ramasomanana, A. Younes, and F. Lehmann.  An efficient lumped mixed hybrid finite element formulation for variably saturated groundwater flow.  Vadose Zone J., 8:352–362, 2009.
  * [13] L. Bergamaschi and M. Putti.  Mixed finite elements and Newton-type linearizations for the solution of Richards’ equation.  Int. J. Numer. Meth. Engng., 45:1025–1046, 1999.
  * [14] H. Berninger, R. Kornhuber, and O. Sander.  Fast and robust numerical solution of the Richards equation in homogeneous soil.  SIAM J. Numer. Anal., 49:2576–2597, 2011.
  * [15] H. Berninger, R. Kornhuber, and O. Sander.  A multidomain discretization of the Richards equation in layered soil.  Comput. Geosci., 19:213–232, 2015.
  * [16] F. Brezzi and M. Fortin.  Mixed and Hybrid Finite Element Method.  Springer, New York, 1991.
  * [17] P. Broadbridge and I. White.  Constant rate rainfall infiltration: A versatile nonlinear model: I. Analytic solution.  Water Resour. Res., 24:145–154, 1988.
  * [18] R. H. Brooks and A. T. Corey.  Properties of porous media affecting fluid flow.  J. Irrig. Drain. Div. Am. Soc. Civ. Eng., 92:61–88, 1966.
  * [19] W. Cai, J. Wang, and K. Wang.  Convergence analysis of Crank-Nicolson Galerkin-Galerkin FEMs for miscible displacement in porous media.  J. Sci. Comput., 83:25, 2020.
  * [20] V. Casulli and P. Zanolli.  A nested Newton-type algorithm for finite volume methods solving Richards’ equation in mixed form.  SIAM J. Sci. Comput., 32:2255–2273, 2010.
  * [21] M. A. Celia, E. T. Bouloutas, and R. L. Zarba.  A general mass-conservative numerical solution for the unsaturated flow equation.  Water Resour. Res., 26:1483–1496, 1990.
  * [22] P. G. Ciarlet.  Introduction to Numerical Linear Algebra and Optimisation.  Cambridge University Press, Cambridge, 1989.
  * [23] T. P. Clement, W. R. Wise, and F. J. Molz.  A physically based, two-dimensional, finite-difference algorithm for modeling variably saturated flow.  J. Hydrol., 161:71–90, 1994.
  * [24] R. L. Cooley.  Some new procedures for numerical solution of variably saturated flow problems.  Water Resour. Res., 19:1271–1285, 1983.
  * [25] M. Crouzeix and F. J. Lisbona.  The convergence of variable-stepsize, variable-formula, multistep methods.  SIAM J. Numer. Anal., 21:512–534, 1984.
  * [26] B. Cumming, T. Moroney, and I. Turner.  A mass-conservative control volume-finite element method for solving Richards’ equation in heterogeneous porous media.  BIT Numer. Math., 51:845–864, 2011.
  * [27] E. Emmrich.  Convergence of the variable two-step BDF time discretisation of nonlinear evolution problems governed by a monotone potential operator.  BIT Numer. Math., 49:297–323, 2009.
  * [28] M. Ethier and Y. Bourgault.  Semi-implicit time-discretization schemes for the bidomain model.  SIAM J. Numer. Anal., 46:2443–2468, 2008.
  * [29] R. Eymard, M. Gutnic, and D. Hilhorst.  The finite volume method for Richards equation.  Comput. Geosci., 3:259–294, 1999.
  * [30] M. W. Farthing, C. E. Kees, and C. T. Miller.  Mixed finite element methods and higher order temporal approximations for variably saturated groundwater flow.  Adv. Water Resour., 26:373–394, 2003.
  * [31] M. W. Farthing and C. T. Miller.  A comparison of high-resolution, finite-volume, adaptive–stencil schemes for simulating advective-dispersive transport.  Adv. Water Resour., 24:29–48, 2000.
  * [32] P. A. Forsyth and M. C. Kropinski.  Monotonicity considerations for saturated-unsaturated subsurface flow.  SIAM J. Sci. Comput., 18:1328–1354, 1997.
  * [33] P. A. Forsyth, Y. S. Wu, and K. Pruess.  Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media.  Adv. Water Resour., 18:25–38, 1995.
  * [34] H. Fujita.  The exact pattern of a concentration-dependent diffusion in a semi-infinite medium, Part ii.  Text. Res. J., 22:823–827, 1952.
  * [35] W. R. Gardner.  Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table.  Soil Sci., 85:228–232, 1958.
  * [36] R. Haverkamp, M. Vauclin, J. Touma, P. J. Wierenga, and G. Vachaud.  A comparison of numerical simulation models for one-dimensional infiltration.  Soil Sci. Soc. Am. J., 41:285–294, 1977.
  * [37] F. Hecht.  New development in FreeFem++.  J. Numer. Math., 20:251–265, 2012.
  * [38] F. Hecht.  FreeFem++: Version 4.1.  https://doc.freefem.org/documentation/index.html, 2019.
  * [39] M. Karthikeyan, T.-S. Tan, and K.-K. Phoon.  Numerical oscillation in seepage analysis of unsaturated soils.  Can. Geotech. J., 38:639–651, 2001.
  * [40] D. Kavetski, P. Binning, and S. W. Sloan.  Noniterative time stepping schemes with adaptive truncation error control for the solution of Richards equation.  Water Resour. Res., 38:29–1–29–10, 2002.
  * [41] S. Keita, A. Beljadid, and Y. Bourgault.  Mass-conservative and positivity preserving second-order semi-implicit methods for high-order parabolic equations.  J. Comput. Phys., Revision submitted, arXiv:2010.11913 [math.NA], 2020.
  * [42] S. Keita, A. Beljadid, and Y. Bourgault.  Efficient second-order semi-implicit finite element method for fourth-order nonlinear diffusion equations.  Comput. Phys. Commun., 258:107588, 2021.
  * [43] F. Lehmann and P. Ackerer.  Comparison of iterative methods for improved solutions of the fluid flow equation in partially saturated porous media.  Transp. Porous Media, 31:275–292, 1998.
  * [44] R. J. LeVeque.  Finite Difference Methods for Ordinary and Partial Differential Equations.  Society for Industrial and Applied Mathematics, 2007.
  * [45] M. C. Leverett.  Capillary behavior of porous solids.  Transaction of AIME, 142:152–169, 1941.
  * [46] F. List and F. A. Radu.  A study on iterative methods for solving Richards’ equation.  Comput. Geosci., 20:341–353, 2016.
  * [47] P. A. Lott, H. F. Walker, C. S. Woodward, and U. M. Yang.  An accelerated Picard method for nonlinear systems related to variably saturated flow.  Adv. Water Resour., 38:92–101, 2012.
  * [48] G. Manzini and S. Ferraris.  Mass-conservative finite volume methods on 2-D unstructured grids for the Richards’ equation.  Adv. Water Resour., 27:1199–1215, 2004.
  * [49] D. McBride, M. Cross, N. Croft, C. Bennett, and J. Gebhardt.  Computational modelling of variably saturated flow in porous media with complex three-dimensional geometries.  Int. J. Numer. Meth. Fluids, 50:1085–1117, 2006.
  * [50] C. T. Miller, G. A. Williams, C. T. Kelley, and M. D. Tocci.  Robust solution of Richards’ equation for nonuniform porous media.  Water Resour. Res., 34:2599–2610, 1998.
  * [51] P. C. D. Milly.  A mass-conservative procedure for time-stepping in models of unsaturated flow.  Adv. Water Resour., 8:32–36, 1985.
  * [52] Y. Mualem.  A new model for predicting the hydraulic conductivity of unsaturated porous media.  Water Resour. Res., 12:513–522, 1976.
  * [53] A. A. H. Oulhaj, C. Cancés, and C. Chainais-Hillairet.  Numerical analysis of a nonlinearly stable and positive control volume finite element scheme for Richards equation with anisotropy.  ESAIM: M2AN, 52:1533–1567, 2018.
  * [54] C. Paniconi, A. A. Aldama, and E. F. Wood.  Numerical evaluation of iterative and noniterative methods for the solution of the nonlinear Richards equation.  Water Resour. Res., 27:1147–1163, 1991.
  * [55] C. Paniconi and M. Putti.  A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems.  Water Resour. Res., 30:3357–3374, 1994.
  * [56] G. W. Parkin, A. W. Warrick, D. E. Elrick, and R. G. Kachanoski.  Analytical solution for one-dimensional drainage: Water stored in a fixed depth.  Water Resour. Res., 31:1267–1271, 1995.
  * [57] J. R. Philip.  Theory of infiltration.  Adv. Hydrosci., 5:215–296, 1969.
  * [58] I. S. Pop, F. Radu, and P. Knabner.  Mixed finite elements for the Richards’ equation: linearization procedure.  J. Comput. Appl. Math., 168:365–373, 2004.
  * [59] I. S. Pop and B. Schweizer.  Regularization schemes for degenerate Richards equations and outflow conditions.  Math. Models Methods Appl. Sci., 21:1685–1712, 2011.
  * [60] F. Radu, I. S. Pop, and P. Knabner.  Order of convergence estimates for an Euler implicit, mixed finite element discretization of Richards’ equation.  SIAM J. Numer. Anal., 42:1452–1478, 2004.
  * [61] F. A. Radu and W. Wang.  Convergence analysis for a mixed finite element scheme for flow in strictly unsaturated porous media.  Nonlinear Anal. Real World Appl., 15:266–275, 2014.
  * [62] E. Rank, C. Katz, and H. Werner.  On the importance of the discrete maximum principle in transient analysis using finite element methods.  Int. J. Numer. Meth. Engng., 19:1771–1782, 1983.
  * [63] L. A. Richards.  Capillary conduction of liquids through porous mediums.  Physics, 1:318–333, 1931.
  * [64] P. J. Ross.  Modeling soil water and solute transport–fast, simplified numerical solutions.  Agron. J., 95:1352–1361, 2003.
  * [65] G. C. Sander, I. F. Cunning, W. L. Hogarth, and J.-Y. Parlange.  Exact solution for nonlinear, nonhysteretic redistribution in vertical soil of finite depth.  Water Resour. Res., 27:1529–1536, 1991.
  * [66] G. C. Sander, J.-Y. Parlange, V. Kühnel, W. L. Hogarth, D. Lockington, and J. P. J. O’Kane.  Exact nonlinear solution for constant flux infiltration.  J. Hydrol., 97:341–346, 1988.
  * [67] E. Schneid, P Knabner, and F. A. Radu.  A priori error estimates for a mixed finite element discretization of the Richards’ equation.  Numer. Math., 98:353–370, 2004.
  * [68] B. Schweizer.  Regularization of outflow problems in unsaturated porous media with dry regions.  J. Differ. Equ., 237:278–306, 2007.
  * [69] M. Slodicka.  A robust and efficient linearization scheme for doubly nonlinear and degenerate parabolic problems arising in flow in porous media.  SIAM J. Sci. Comput., 23:1593–1614, 2002.
  * [70] P. Sochala.  Méthodes numériques pour les écoulements souterrains et couplage avec le ruissellement.  PhD thesis, École Nationale des Ponts et Chaussées, 2008.
  * [71] H. Suk and E. Park.  Numerical solution of the Kirchhoff-transformed Richards equation for simulating variably saturated flow in heterogeneous layered porous media.  J. Hydrol., 579:124213, 2019.
  * [72] F. T. Tracy.  Clean two- and three-dimensional analytical solutions of Richards’ equation for testing numerical solvers.  Water Resour. Res., 42:W08503, 2006.
  * [73] F. T. Tracy.  Analytical and numerical solutions of Richards’ equation with discussions on relative hydraulic conductivity.  In Prof. Lakshmanan Elango, editor, Hydraulic Conductivity. IntechOpen, 2011.
  * [74] L. Traverso, T. N. Phillips, and Y. Yang.  Mixed finite element methods for groundwater flow in heterogeneous aquifers.  Comput. Fluids, 88:60–80, 2013.
  * [75] H. Uzawa.  Iterative methods for concave programming.  In K. J. Arrow, L. Hurwicz, H. Uzawa (eds.), Studies in Linear and Non-Linear Programming, pages 154–165. Stanford University Press, Stanford, 1958.
  * [76] M. Th. van Genuchten.  A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.  Soil Sci. Soc. Am. J., 44:892–898, 1980.
  * [77] A. Warrick.  Soil Water Dynamics.  Oxford University Press, New York, 2003.
  * [78] A. Younes, P. Ackerer, and F. Lehmann.  A new mass lumping scheme for the mixed hybrid finite element method.  Int. J. Numer. Meth. Engng., 67:89–107, 2006.
  * [79] Y. Zha, L. Shi, M. Ye, and J. Yang.  A generalized Ross method for two- and three-dimensional variably saturated flow.  Adv. Water Resour., 54:67–77, 2013.
  * [80] Y. Zha, M. Tso, L. Shi, and J. Yang.  Comparison of noniterative algorithms based on different forms of Richards’ equation.  Environ. Model. Assess., 21:357–370, 2016.