# Stability Analysis of the Explicit Difference Scheme for Richards Equation

^{1}

^{2}

^{3}

^{*}

Next Article in Journal

Next Article in Special Issue

Next Article in Special Issue

Previous Article in Journal

Previous Article in Special Issue

Previous Article in Special Issue

School of Mathematical and Physical Sciences, Dalian University of Technology, Panjin 124221, China

Institute of Mathematics for Industry, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan

College of Sciences, Northeastern University, Shenyang 110819, China

Author to whom correspondence should be addressed.

Received: 6 January 2020
/
Revised: 14 March 2020
/
Accepted: 15 March 2020
/
Published: 18 March 2020

(This article belongs to the Special Issue Applications of Nonlinear Diffusion Equations)

A stable explicit difference scheme, which is based on forward Euler format, is proposed for the Richards equation. To avoid the degeneracy of the Richards equation, we add a perturbation to the functional coefficient of the parabolic term. In addition, we introduce an extra term in the difference scheme which is used to relax the time step restriction for improving the stability condition. With the augmented terms, we prove the stability using the induction method. Numerical experiments show the validity and the accuracy of the scheme, along with its efficiency.

A knowledge of the way of infiltration of water into the ground is crucially important for predicting disasters, such as river floods and landslides, when heavy rain attacks. A classical mathematical model for describing the fluid motion in unsaturated zone in a porous medium is the Richards equation, a nonlinear degenerate advection-diffusion equation. Two main research directions arose recently for the Richards equation. One is to find exact solutions by finding a specific form of coefficient functions so as to make the equation completely integrable. An exact solution helps to clearly capture the physical mechanism of the phenomena and to pursue controllability [1,2,3]. The other research direction is to seek an approximate solution by numerical methods, for coefficient functions to match with a practical situation. Finite elements, finite difference or finite volumes methods are carried out in [4,5] and the reference therein. The above schemes commonly used fully implicit schemes based on a backward Euler format. Adaptive time stepping is studied in [6,7]. In these studies, some conservative schemes are devised and numerical tests show that they have good stability and some order accuracy, but no theoretical proof is given.

So far, we have only known one paper to prove the stability for the mixed finite element discretization of the Richards equation in [8,9]. In [8], they introduced an implicit mixed element scheme and applied the Kirchhoff transformation to deal with the degeneracy of the Richards equation. The Kirchhoff transformation could be used in the continuous inner product, but we have to take the discrete inner product in the analysis for the difference scheme, so we take a new way which is by adding a perturbation to the coefficient function of parabolic term to overcome the degeneracy.

In [8], the scheme they used is implicit, the stability condition is certainly superior to the explicit scheme. So they do not pay attention to the stability condition in the analysis for the implicit scheme. However, the explicit numerical schemes always are stable only in rigorous restriction for mesh ratio. To improve the stability condition, we introduced an extra term in the difference scheme to relax the time step restriction for improving the stability condition. Also, the theoretical analysis for the implicit numerical schemes and the explicit difference scheme is completely different.

As we know, these early attempts are all implicit schemes based on backward Euler difference scheme and the central difference scheme, although in certain case, the implicit scheme may have to be used to avoid instability. However a strongly nonlinear algebraic system must be solved at each time level, even though these iterative methods [10,11] are used, it still needs huge calculation. Explicit scheme is a good choice to improve the computation efficiency, but the classical explicit scheme cannot be used for the Richards equation due to its degeneracy and the severe time step length restriction.

The main purpose of this work is to provide an efficient explicit numerical scheme for the Richards equation and prove the stability. The key objectives of this work are threefold: First, we add a perturbation to the coefficient function of parabolic term to overcome the degeneracy. Secondly, we introduce a stabilization term with constant coefficient in the difference scheme to relax the stability restriction on the time step. Please note that a similar technique has been used in the simulation of the Cahn–Hilliard equation [12] and the MBE models [13]. The Cahn–Hilliard equation and the MBE models are fourth-order parabolic partial differential equations, so they introduced a second-order stabilization term in the Fourier spectral scheme and finite element scheme respectively. The Richards equation is a second-order equation, so the stabilization term which is added in the explicit difference scheme is completely different. Finally, we prove the stability by induction method and perform some numerical experiments.

The Richards equation could be written in three equivalent forms, with either pressure head $h\left[L\right]$ or moisture content $\theta [{L}^{3}/{L}^{3}]$ as the dependent variable. We recall that the hydraulic head $h+z$ is partitioned into the pressure head $h=p/\left(\rho g\right)$ and the gravity head z, the vertical coordinate increasing upwards, with the pressure p normalized by the gravity force. Here $\rho $ is the mass density of the fluid and g is the gravity acceleration. The constitutive relationship between $\theta =\theta (z,t)$ and $h=h(z,t)$ allows the conversion from one to another. The three forms can be identified as h-based, $\theta $-based, and mixed:

- h-based$$C\left(h\right)\frac{\partial h}{\partial t}-\nabla \xb7K\left(h\right)\nabla h-\frac{\partial K}{\partial z}=0,$$
- $\theta $-based$$\frac{\partial \theta}{\partial t}-\nabla \xb7D\left(\theta \right)\nabla \theta -\frac{\partial K}{\partial z}=0,$$
- mixed$$\frac{\partial \theta}{\partial t}-\nabla \xb7K\left(h\right)\nabla h-\frac{\partial K}{\partial z}=0,$$

where the real-valued functions $C\left(h\right)\equiv d\theta /dh$, $K\left(h\right)$, and $D\left(\theta \right)\equiv K\left(\theta \right)/C\left(\theta \right)$ respectively denote the specific moisture capacity function $[1/L]$, the unsaturated hydraulic conductivity $[L/T]$ and the unsaturated diffusivity $[{L}^{2}/T]$. The coefficient $K\left(h\right)$ describes the ease with which water can move through pore spaces, and depends on the intrinsic permeability of the material, degree of saturation, and the density and the viscosity of the fluid. The porous medium is assumed to be isotropic.

We consider the h-based form with the datum reported by Haverkamp et al. [14,15] which is used to solve an example of infiltration into soil column.
where $\theta \left(h\right)$ represents the moisture content. ${\theta}_{r}$ and ${\theta}_{s}$ are initial and saturated moisture content respectively. Moreover, $C\left(h\right)=d\theta /dh$, simple calculations show that

$$\theta \left(h\right)={\displaystyle \frac{\alpha ({\theta}_{s}-{\theta}_{r})}{\alpha +{\left|h\right|}^{\beta}}}+{\theta}_{r},\phantom{\rule{1.em}{0ex}}K\left(h\right)={K}_{s}{\displaystyle \frac{A}{A+{\left|h\right|}^{\gamma}}},$$

$$C\left(h\right)={\theta}^{\prime}\left(h\right)={\displaystyle \frac{\alpha ({\theta}_{s}-{\theta}_{r})\beta {\left|h\right|}^{\beta -1}}{{(\alpha +|h|}^{\beta}{)}^{2}}},\phantom{\rule{1.em}{0ex}}{K}^{\prime}\left(h\right)={\displaystyle \frac{{K}_{s}A\gamma {\left|h\right|}^{\gamma -1}}{{(A+|h|}^{\gamma}{)}^{2}}}.$$

For the given $\theta \left(h\right)$ and $K\left(h\right)$, we consider the following data [16].

$$\begin{array}{c}\alpha =1.611\times {10}^{6},\phantom{\rule{1.em}{0ex}}{\theta}_{s}=0.287,\phantom{\rule{1.em}{0ex}}{\theta}_{r}=0.075,\phantom{\rule{1.em}{0ex}}\beta =3.96,\hfill \\ {K}_{s}=0.00944\phantom{\rule{3.33333pt}{0ex}}\mathrm{cm}/\mathrm{s},\phantom{\rule{1.em}{0ex}}A=1.175\times {10}^{6},\phantom{\rule{1.em}{0ex}}\gamma =4.74.\hfill \end{array}$$

$$h(40\phantom{\rule{3.33333pt}{0ex}}\mathrm{cm},t)={h}_{top}=-20.7\phantom{\rule{3.33333pt}{0ex}}\mathrm{cm},\phantom{\rule{1.em}{0ex}}h(0,t)={h}_{bottom}=-61.5\phantom{\rule{3.33333pt}{0ex}}\mathrm{cm}.$$

These data provide some real numbers from an example of infiltration into soil column. From the data, we can verify that there are upper bounds for $K\left(h\right)$ and ${K}^{\prime}\left(h\right)$ easily, i.e., $K\left(h\right)\le K<{K}_{s},{K}^{\prime}\left(h\right)\le {K}_{1}$ for $h\in \mathbb{R}$. This is to be used in the stability proof.

Let $\Delta z=L/M$ be the uniform step length, where M is a positive integer. We divide the domain of time T with N segments, let $\tau =T/N,{t}^{n}=n\tau $ be the uniform time length. Then for a function $h(t,z)$, denote ${H}_{i}^{n}=h({z}_{i},{t}_{n})$, where ${z}_{i}=m\Delta z,m=0,1,\dots ,M$, and $\overline{\mathsf{\Omega}}=\left\{{z}_{i}\right|i=0,1,\dots M\}$, ${t}_{n}=n\tau ,n=0,1,\dots ,N$. Let $\lambda =\tau /\Delta {z}^{2}$ be the mesh ratios.

Define the following difference operators

$${\delta}_{t}{H}_{i}^{n}=\frac{{H}_{i}^{n+1}-{H}_{i}^{n}}{\tau},0\le n\le N-1,$$

$${\nabla}_{h}{H}_{i}^{n}=\frac{{H}_{i+1}^{n}-{H}_{i-1}^{n}}{2\Delta z},\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}{\Delta}_{h}{H}_{i}^{n}=\frac{{H}_{i+1}^{n}-2{H}_{i}^{n}+{H}_{i-1}^{n}}{{\Delta z}^{2}},1\le i\le M-1.$$

Now, we introduce the discrete ${L}^{2}$ inner product as $\langle u,v\rangle =\Delta z\{\frac{1}{2}({u}_{0}{v}_{0}+{u}_{M}{v}_{M})+{\displaystyle \sum _{i=1}^{M-1}}{u}_{i}{v}_{i}\},$ and the coresponding discrete ${L}^{2}$ norm is ${\parallel v\parallel}_{h}={\langle v,v\rangle}^{\frac{1}{2}}$. Moreover, the discrete ${H}^{1}$ seminorm ${|\xb7|}_{1,h}$ and the discrete maximum norm ${|\xb7|}_{\infty ,h}$ are defined as

$${\left|v\right|}_{1,h}={\left[\Delta z\sum _{i=1}^{M-1}{\left({\nabla}_{h}{v}_{i}\right)}^{2}\right]}^{\frac{1}{2}}{,\phantom{\rule{1.em}{0ex}}\left|v\right|}_{\infty ,h}=\underset{i}{sup}\left|{v}_{i}\right|,i=1,\dots M.$$

A classical first-order explicit difference scheme is

$$C\left({H}_{i}^{n}\right){\delta}_{t}{H}_{i}^{n}-{\nabla}_{h}\left(K\left({H}_{i}^{n}\right){\nabla}_{h}{H}_{i}^{n}\right)-{K}^{\prime}\left({H}_{i}^{n}\right){\nabla}_{h}{H}_{i}^{n}=0.$$

For the degeneracy of the Richards equation, a special trick to handle the nonlinear parabolic term is devised. We add a positive constant ${\u03f5}_{1}$ to $C\left({H}_{i}^{n}\right)$ in (7) (in fact, ${\u03f5}_{1}$ can be seen as a small positive bound of $C\left({H}_{i}^{n}\right)$). Then modified first-order explicit scheme is of the form

$$C\left({H}_{i}^{n}\right){\delta}_{t}{H}_{i}^{n}+{\u03f5}_{1}{\delta}_{t}{H}_{i}^{n}-{\nabla}_{h}\left(K\left({H}_{i}^{n}\right){\nabla}_{h}{H}_{i}^{n}\right)-{K}^{\prime}\left({H}_{i}^{n}\right){\nabla}_{h}{H}_{i}^{n}=0.$$

With the Richards equation featured by a convection dominated diffusion problems, numerical experiments show that the stability of (8) is restricted by the mesh ratio, meaning that the scheme is stable only in very tiny time step, as is expected. So we add extra diffusion terms in (8) to improve the stability condition so that relax the restriction of the time step.
where ${\u03f5}_{2}$ is a positive constant to be determined so as to improve the stability condition.

$$C\left({H}_{i}^{n}\right){\delta}_{t}{H}_{i}^{n}+{\u03f5}_{1}{\delta}_{t}{H}_{i}^{n}-{\u03f5}_{2}({\Delta}_{h}{H}_{i}^{n+1}-{\Delta}_{h}{H}_{i}^{n})-{\nabla}_{h}\left(K\left({H}_{i}^{n}\right){\nabla}_{h}{H}_{i}^{n}\right)-{K}^{\prime}\left({H}_{i}^{n}\right){\nabla}_{h}{H}_{i}^{n}=0,$$

In [17], for the one-dimensional Richards equation, we also established a linearized semi-implicit finite difference scheme and analyzed the stability. Compared to the scheme of [17], the explicit difference scheme (9) is to avoid solving a linear algebraic equations at every time step. If we divide the domain of time with N segments, the explicit difference scheme could reduce the computational cost to $1/N$ of that. So the explicit difference scheme (9) is more concise and the speed of its numerical simulation is faster.

The scheme (9) is stable with ${L}_{\infty}$ discrete norm, when the time-step length satisfies $\tau <({C}_{\u03f5}+{\u03f5}_{1}){K}_{s}/{K}_{1}^{2}$, where the ${C}_{\u03f5}\ge 0$ is the lower bound of the $C\left({H}_{i}^{n}\right)$.

From (4), the lower bound of the $K\left({H}_{i}^{n}\right)$ is non-zero for the bounded $\left|h\right|$. Now we assume that there is a constant ${k}_{\u03f5}>0$ such that ${k}_{\u03f5}<K\left({H}_{i}^{n}\right),i=0,1,\dots ,M$ for fixed n. Taking the inner product of (9) with ${H}^{n+1}-{H}^{n}$ gives
where ${I}_{1}$–${I}_{4}$ satisfy
and

$${I}_{1}+{I}_{2}+{I}_{3}-{I}_{4}=0,$$

$${I}_{1}=\langle (C\left({H}^{n}\right)+{\u03f5}_{1})\frac{{H}^{n+1}-{H}^{n}}{\tau},{H}^{n+1}-{H}^{n}\rangle \ge ({C}_{\u03f5}+\frac{{\u03f5}_{1}}{\tau}){\parallel {H}^{n+1}-{H}^{n}\parallel}^{2},$$

$${I}_{2}=-\langle {\u03f5}_{2}({\Delta}_{h}{H}^{n+1}-{\Delta}_{h}{H}^{n}),{H}^{n+1}-{H}^{n}\rangle ={\u03f5}_{2}{\parallel {\nabla}_{h}{H}^{n+1}-{\nabla}_{h}{H}^{n}\parallel}^{2},$$

$$\begin{array}{cc}\hfill {I}_{3}=& \langle K\left({H}^{n}\right){\nabla}_{h}{H}^{n},{\nabla}_{h}{H}^{n+1}-{\nabla}_{h}{H}^{n}\rangle \hfill \\ \hfill =& \langle K\left({H}^{n}\right)({\nabla}_{h}{H}^{n+1}-{\nabla}_{h}{H}^{n}),{\nabla}_{h}{H}^{n}-{\nabla}_{h}{H}^{n+1}+{\nabla}_{h}{H}^{n+1}\rangle \hfill \\ \hfill \ge & -K\parallel {\nabla}_{h}{H}^{n+1}-{\nabla}_{h}{H}^{n}{\parallel}^{2}+\frac{{K}_{\u03f5}}{2}(\parallel {\nabla}_{h}{H}^{n+1}{\parallel}^{2}-\parallel {\nabla}_{h}{H}^{n}{\parallel}^{2}),\hfill \end{array}$$

$${I}_{4}=\langle {K}^{\prime}\left({H}^{n}\right){\nabla}_{h}{H}^{n},{H}^{n+1}-{H}^{n}\rangle \le \frac{{\u03f5}_{1}}{2\tau}\parallel {H}^{n+1}-{H}^{n}{\parallel}^{2}+\frac{{K}_{1}^{2}\tau}{2{\u03f5}_{1}}{\parallel {\nabla}_{h}{H}^{n}\parallel}^{2}.$$

Summing up, we obtain
for ${\u03f5}_{2}\ge K$. Simple calculation shows that
which implies that there exists a positive constant ${c}_{1}$, being independent of $\Delta z$ and $\tau $, such that $\parallel {\nabla}_{h}{H}^{n+1}{\parallel}^{2}\le {c}_{1}$.

$$\begin{array}{c}({C}_{\u03f5}+\frac{{\u03f5}_{1}}{2\tau})\parallel {H}^{n+1}-{H}^{n}{\parallel}^{2}+({\u03f5}_{2}-K){\parallel {\nabla}_{h}{H}^{n+1}-{\nabla}_{h}{H}^{n}\parallel}^{2}\hfill \\ +\frac{{K}_{s}}{2}(\parallel {\nabla}_{h}{H}^{n+1}{\parallel}^{2}-\parallel {\nabla}_{h}{H}^{n}{\parallel}^{2})\le \frac{{K}_{1}^{2}\tau}{2{\u03f5}_{1}}{\parallel {\nabla}_{h}{H}^{n}\parallel}^{2},\hfill \end{array}$$

$$\parallel {\nabla}_{h}{H}^{n+1}{\parallel}^{2}\le (1+\frac{{K}_{1}^{2}\tau}{{\u03f5}_{1}{K}_{s}}){\parallel {\nabla}_{h}{H}^{n}\parallel}^{2},$$

We are now prepared to prove the stability with respect to the discrete ${L}^{2}$ norm. Taking the inner product of (9) with ${H}^{n+1}$, we have
where

$${E}_{1}+{E}_{2}+{F}_{1}\le {E}_{3}+{F}_{2}+{E}_{5},$$

$${E}_{1}=\langle (C\left({H}^{n}\right)+{\u03f5}_{1})\frac{{H}^{n+1}-{H}^{n}}{\tau},{H}^{n+1}\rangle \ge \frac{{C}_{\u03f5}+{\u03f5}_{1}}{2\tau}(\parallel {H}^{n+1}{\parallel}^{2}-\parallel {H}^{n}{\parallel}^{2}),$$

$${E}_{2}=-\langle {\u03f5}_{2}{\Delta}_{h}{H}^{n+1},{H}^{n+1}\rangle ={\u03f5}_{2}{\parallel {\nabla}_{h}{H}^{n+1}\parallel}^{2},$$

$${E}_{3}=-\langle {\u03f5}_{2}{\Delta}_{h}{H}^{n},{H}^{n+1}\rangle \le \frac{{\u03f5}_{2}}{2}(\parallel {\nabla}_{h}{H}^{n+1}{\parallel}^{2}+\parallel {\nabla}_{h}{H}^{n}{\parallel}^{2}),$$

$$\begin{array}{cc}\hfill {E}_{4}=& \langle K\left({H}^{n}\right){\nabla}_{h}{H}^{n},{\nabla}_{h}{H}^{n+1}\rangle =\langle K\left({H}^{n}\right){\nabla}_{h}{H}^{n},{\nabla}_{h}{H}^{n+1}-{\nabla}_{h}{H}^{n}+{\nabla}_{h}{H}^{n}\rangle ,\hfill \\ \hfill =& \underset{{F}_{1}}{\underset{\u23df}{\langle K\left({H}^{n}\right){\nabla}_{h}{H}^{n},{\nabla}_{h}{H}^{n}\rangle}}+\underset{F2}{\underset{\u23df}{\langle K\left({H}^{n}\right){\nabla}_{h}{H}^{n},{\nabla}_{h}{H}^{n+1}-{\nabla}_{h}{H}^{n}\rangle}},\hfill \end{array}$$

$${E}_{5}=\langle {K}^{\prime}\left({H}^{n}\right){\nabla}_{h}{H}^{n},{H}^{n+1}\rangle \le \frac{{K}_{\u03f5}}{2}\parallel {\nabla}_{h}{H}^{n}{\parallel}^{2}+\frac{{K}_{1}^{2}}{2{K}_{\u03f5}}{\parallel {H}^{n+1}\parallel}^{2}.$$

Applying Young’s inequality, we obtain

$$\begin{array}{c}{F}_{1}\ge {K}_{\u03f5}\parallel {\nabla}_{h}{H}^{n}{\parallel}^{2},\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}{F}_{2}\le \frac{K}{2}(\parallel {\nabla}_{h}{H}^{n+1}{\parallel}^{2}-\parallel {\nabla}_{h}{H}^{n}{\parallel}^{2}).\hfill \end{array}$$

Eventually, we obtain

$$(\frac{{C}_{\u03f5}+{\u03f5}_{1}}{2\tau}-\frac{{K}_{1}^{2}}{2{K}_{\u03f5}})\parallel {H}^{n+1}{\parallel}^{2}+(\frac{{K}_{\u03f5}}{2}-\frac{{\u03f5}_{2}}{2}+\frac{K}{2})\parallel {\nabla}_{h}{H}^{n}{\parallel}^{2}+(\frac{{\u03f5}_{2}}{2}-\frac{K}{2})\parallel {\nabla}_{h}{H}^{n+1}{\parallel}^{2}\le \frac{{C}_{\u03f5}+{\u03f5}_{1}}{2\tau}{\parallel {H}^{n}\parallel}^{2}.$$

We conclude that for $\tau <({C}_{\u03f5}+{\u03f5}_{1}){K}_{s}/{K}_{1}^{2}$, there is a positive constant ${c}_{2}$ which is independent of $\Delta z$ and $\tau $ such that

$$\parallel {H}^{n+1}{\parallel}^{2}\le {c}_{2}.$$

Similar estimation techniques were used in other useful applications [18].

We exploit the constant ${\u03f5}_{2}$ to improve the stability. If ${\u03f5}_{2}>K$, the scheme (9) is stable. Please note that $K<{K}_{s}$ and ${K}_{s}=0.00944cm/s$ makes it possible to take ${\u03f5}_{2}$ sufficiently small in case of excessive errors.

To make sure the numerical scheme (9) is convergent, the time step should be chosen to satisfy $\tau <({C}_{\u03f5}+{\u03f5}_{1}){K}_{s}/{K}_{1}^{2}$. If ${C}_{\u03f5}=0$, we have to take a non-zero perturbation ${\u03f5}_{1}$ to ensure the stability of the scheme. Because of ${C}_{\u03f5}>0$, we are allowed to take ${\u03f5}_{1}=0$ in numerical experiments.

In this section, we illustrate the numerical stability by a numerical experiment of the infiltration process based on a generalized h-based Richards equation. Since it is difficult to obtain the exact solution of this model, to verify the theoretical results, we take the following non-homogeneous model,
where the boundary conditions remain unchanged as (6). If we suppose an exact solution $h=-1.02-20.7+t(z-40)/\left(4T\right)$, a simple calculation shows that

$$\left\{\begin{array}{c}C\left(h\right)\frac{\partial h}{\partial t}-\nabla \xb7\left(K\left(h\right)\nabla h\right)-\frac{\partial K}{\partial z}=g(z,t),\hfill \\ h(z,0)=-1.02z-20.7,\hfill \end{array}\right.$$

$$\begin{array}{cc}\hfill g(z,t)=& \frac{5808800331328389z(z-40){\left(\frac{51z}{50}-\frac{tz(z-40)}{4T}+\frac{207}{10}\right)}^{74/25}}{17179869184T{\left[{\left(\frac{51z}{50}-\frac{tz(z-40)}{4T}+\frac{207}{10}\right)}^{99/25}+1611000\right]}^{2}}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& -\frac{5546t}{{\left(\frac{51z}{50}-\frac{tz(z-40)}{4T}+\frac{207}{10}\right)}^{237/50}+1175000}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& -\frac{3613000706430075\left(\frac{t(z-20)}{2T}-\frac{51}{50}\right){\left(\frac{51z}{50}-\frac{tz(z-40)}{4T}+\frac{207}{10}\right)}^{187/50}}{68719476736{\left[{\left(\frac{51z}{50}-\frac{tz(z-40)}{4T}+\frac{207}{10}\right)}^{237/50}+1175000\right]}^{2}}.\hfill \end{array}$$

By using the scheme (9), we show, in Figure 1, the variation trend of the pressure head with depth for the time interval from 10 s to 360 s, with the choice of ${\u03f5}_{1}=0$ and ${\u03f5}_{2}=0$.

In [14], the authors used an implicit numerical scheme, and took a large time step, for instance, 10 s, 30 s and 120 s, to save the computational workload. In this paper, we take the time step as small as ${10}^{-3}$ s. Correspondingly the space steps that they used are much larger than ours. The different scheme and the large grid gap may bring some but tolerable discrepancy.

We take $T=100$ s and $M=200$ to test the stability of the scheme with different time steps and different choices of ${\u03f5}_{2}$. The results are listed in Table 1. It is evident that the improvement in the stability by use of the extra terms is significant. Moreover, in Table 2, we set ${\u03f5}_{1}=0,{\u03f5}_{2}=0,M=200$ and $T=1$ s, and confirm that the expected order of convergence is obtained.

Figure 2 shows the linear relationship between ${\u03f5}_{1}$ and errors.

In this work, a stable explicit scheme for the Richards equation was developed and analyzed. We proposed techniques to avoid the degeneracy of the Richards equation and improve the stability condition of the finite difference scheme. A numerical example is provided to verify our theoretical analysis. Demonstration of the numerical stability over a long time, along with the error estimate as shown by Figure 2, is indicative of the physical stability of a typical solution of the Richards equation; infinitesimal perturbations to the solution do not grow. A rigorous mathematical analysis of the stability of the traveling-wave solution and its relevance to the numerical stability call for an independent investigation.

Compared to implicit numerical schemes and linearized numerical schemes, stable explicit numerical schemes improve the calculation efficiency. This paper is a first step toward the explicit difference schemes for the Richards equation, we only analysis the stability of such a scheme. It is our ongoing work to extend other high order accuracy explicit difference schemes and estimate the errors.

Conceptualization, Y.F., F.L. and X.Z.; Methodology, F.L. and X.Z.; Calculation, F.L.; Validation, X.Z. and Y.F.; Formal Analysis, F.L.; Data Curation, F.L.; Writing—Original Draft Preparation, F.L.; Writing—Review Editing, Y.F., F.L. and X.Z. All authors have read and agreed to the published version of the manuscript.

F.L. was supported by Fundamental Research Funds for Central Universities (no. DUT19RC(4)038), Y.F. was supported by a Grant-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (grant no.19K03672), and X.Z. was supported by China Postdoctoral Science Foundation (no. 2015M581689).

The authors declare no conflict of interest.

- Broadbridge, P.; Daly, E.; Goard, J. Exact Solutions of the Richards Equation With Nonlinear Plant-Root Extraction. Water Resour. Res.
**2017**, 53, 9679–9691. [Google Scholar] [CrossRef] - Broadbridge, P.; Edwards, M.; Kearton, J. Closed-form solutions for unsaturated flow under variable flux boundary conditions. Adv. Water Res.
**1996**, 19, 207–213. [Google Scholar] [CrossRef] - Broadbridge, P.; Triadis, D.; Hill, J.M. Infiltration from supply at constant water content: An integrable model. J. Eng. Math.
**2009**, 64, 193–206. [Google Scholar] [CrossRef] - Kumar, C. A numerical simulation model for one-dimensional infiltration. ISH J. Hyddraul. Eng.
**1998**, 4, 5–15. [Google Scholar] [CrossRef] - Lopez, F.; Vacca, G. Spectral properties and conservation laws in mimetic finite difference methods for PDEs. J. Comput. Appl. Math.
**2016**, 292, 760–784. [Google Scholar] [CrossRef] - Kavetski, D.; Binning, P.; Sloan, S.W. Adaptive time stepping and error control in a mass conservative numerical solution of the mixed form of Richards equation. Adv. Water Res.
**2001**, 24, 595–605. [Google Scholar] [CrossRef] - Williams, G.A.; Miller, C.T. An evaluation of temporally adaptive transformation approaches for solving Richards equation. Adv. Water Res.
**1999**, 22, 831–840. [Google Scholar] [CrossRef] - Radu, F.; Pop, I.S.; Knabner, P. Order of convergence estimates for an Euler implicit, mixed finite element discretization of Richards equation. SIAM J. Numer. Anal.
**2004**, 42, 1452–1478. [Google Scholar] [CrossRef] - Both, J.W.; Kumar, K.; Nordbotten, J.M.; Pop, I.S.; Radu, F.A. Iterative Linearisation Schemes for Doubly Degenerate Parabolic Equations. In Numerical Mathematics and Advanced Applications; Springer: Cham, Switzerland, 2017; Volume 126. [Google Scholar]
- Kuraz, M.; Mayer, P.; Pech, P. Solving the nonlinear Richards equation model with adaptive domain decomposition. J. Comput. Appl. Math.
**2014**, 270, 2–11. [Google Scholar] [CrossRef] - List, F.; Radu, F.A. A study on iterative methods for solving Richards’ equation. Comput. Geosci.
**2016**, 20, 341–353. [Google Scholar] [CrossRef] - Zhu, J.; Chen, L.-Q.; Shen, J.; Tikare, V. Coarsening kinetics from a variable-mobility Cahn-Hilliard equation: Application of a semi-implicit Fourier spectral method. Phys. Rev. E
**1999**, 60, 3564–3572. [Google Scholar] [CrossRef] [PubMed] - Xu, C.; Tang, T. Stability analysis of large time-stepping methods for epitaxial growth models. SIAM J. Numer. Anal.
**2006**, 44, 1759–1779. [Google Scholar] [CrossRef] - Haverkamp, R.; Vauclin, M. A note on estimating finite difference interblock hydraulic conductivity values for transient unsaturated flow problems. Water Resour. Res.
**1979**, 15, 181–187. [Google Scholar] [CrossRef] - Pop, I.S.; Radu, F.; Knabner, P. Mixed finite elements for the Richards’ equation: Linearization procedure. J. Comput. Appl. Math.
**2004**, 168, 365–373. [Google Scholar] [CrossRef] - Celia, M.A.; Bouloutas, E.T. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res.
**1990**, 26, 1483–1496. [Google Scholar] [CrossRef] - Liu, F.; Fukumoto, Y.; Zhao, X. A linearized finite difference scheme for the Richards equation under variable-flux boundary conditions. J. Sci. Comput.
**2020**, in press. [Google Scholar] [CrossRef] - Shang, Y. Fixed-time group consensus for multi-agent systems with non-linear dynamics and uncertainties. IET control Theory Appl.
**2018**, 12, 395–404. [Google Scholar] [CrossRef]

$\mathit{\tau}$ | $\mathit{\tau}=0.4$ | $\mathit{\tau}=0.2$ | $\mathit{\tau}=0.1$ | $\mathit{\tau}=0.05$ | $\mathit{\tau}=0.025$ | $\mathit{\tau}=0.0125$ | ||
---|---|---|---|---|---|---|---|---|

Accuracy | ||||||||

${\mathit{\u03f5}}_{2}$ | ||||||||

${\u03f5}_{2}=0$ | Unstable | Unstable | Unstable | Unstable | $5.60\times {10}^{-4}$ | $3.26\times {10}^{-5}$ | ||

${\u03f5}_{2}=0.0001$ | Unstable | Unstable | Unstable | Unstable | $1.38\times {10}^{-4}$ | $3.75\times {10}^{-5}$ | ||

${\u03f5}_{2}=0.0005$ | Unstable | $1.88\times {10}^{-1}$ | $5.47\times {10}^{-2}$ | $5.00\times {10}^{-3}$ | $3.78\times {10}^{-4}$ | $1.90\times {10}^{-4}$ | ||

${\u03f5}_{2}=0.001$ | $9.10\times {10}^{-2}$ | $1.44\times {10}^{-2}$ | $4.00\times {10}^{-3}$ | $1.50\times {10}^{-3}$ | $7.54\times {10}^{-4}$ | $3.79\times {10}^{-4}$ |

N | 1000 | 2000 | 4000 | 8000 | 16,000 |
---|---|---|---|---|---|

${|h-H|}_{\infty ,h}$ | $1.90\times {10}^{-3}$ | $9.65\times {10}^{-4}$ | $4.82\times {10}^{-4}$ | $2.41\times {10}^{-4}$ | $1.21\times {10}^{-4}$ |

Ratio | Non | 0.98 | 1.00 | 1.00 | 0.99 |

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).