Stability Analysis of the Explicit Difference Scheme for Richards Equation

Abstract

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.

1. Introduction

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 organization of the paper as follows. In Section 2, an explicit difference scheme is given for the Richards equation. In Section 3, we prove the stability of such scheme. In Section 4, some numerical experiments are given. We conclude the paper in Section 5.

2. Richards Equation and the Explicit Difference Scheme

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 ·K\left(h\right)\nabla h-\frac{\partial K}{\partial z}=0,$$

  (1)
• $\theta$-based

  $$\frac{\partial \theta }{\partial t}-\nabla ·D\left(\theta \right)\nabla \theta -\frac{\partial K}{\partial z}=0,$$

  (2)
• mixed

  $$\frac{\partial \theta }{\partial t}-\nabla ·K\left(h\right)\nabla h-\frac{\partial K}{\partial z}=0,$$

  (3)

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.

$$\theta \left(h\right)={\displaystyle \frac{\alpha ({\theta }_s-{\theta }_r)}{\alpha +{\left|h\right|}^{\beta }}}+{\theta }_r,K\left(h\right)=K_s{\displaystyle \frac{A}{A+{\left|h\right|}^{\gamma }}},$$

(4)

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

$$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}},K^{\prime }\left(h\right)={\displaystyle \frac{K_{s}A\gamma {\left|h\right|}^{\gamma -1}}{{(A+|h|}^{\gamma }{)}^2}}.$$

(5)

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,{\theta }_s=0.287,{\theta }_r=0.075,\beta =3.96,\hfill \\ K_s=0.00944\mathrm{cm}/\mathrm{s},A=1.175\times {10}^6,\gamma =4.74.\hfill \end{array}$$

$$h(40\mathrm{cm},t)=h_{top}=-20.7\mathrm{cm},h(0,t)=h_{bottom}=-61.5\mathrm{cm}.$$

(6)

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},{\Delta }_h{H}_i^n=\frac{H_{i+1}^n-2H_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 ${|·|}_{1,h}$ and the discrete maximum norm ${|·|}_{\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}}{,\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.$$

(7)

For the degeneracy of the Richards equation, a special trick to handle the nonlinear parabolic term is devised. We add a positive constant ${ϵ}_1$ to $C\left(H_i^n\right)$ in (7) (in fact, ${ϵ}_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+{ϵ}_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.$$

(8)

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.

$$C\left(H_i^n\right){\delta }_t{H}_i^n+{ϵ}_1{\delta }_t{H}_i^n-{ϵ}_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,$$

(9)

where ${ϵ}_2$ is a positive constant to be determined so as to improve the stability condition.

Remark 1.

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.

3. Stability Analysis

Theorem 1.

The scheme (9) is stable with $L_{\infty }$ discrete norm, when the time-step length satisfies $\tau <(C_{ϵ}+{ϵ}_1)K_s/K_1^2$, where the $C_{ϵ}\ge 0$ is the lower bound of the $C\left(H_i^n\right)$.

Proof. 

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_{ϵ}>0$ such that $k_{ϵ}<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

$$I_1+I_2+I_3-I_4=0,$$

where $I_1$–$I_4$ satisfy

$$I_1=\langle (C\left(H^n\right)+{ϵ}_1)\frac{H^{n+1}-H^n}{\tau },H^{n+1}-H^n\rangle \ge (C_{ϵ}+\frac{{ϵ}_1}{\tau }){\parallel H^{n+1}-H^n\parallel }^2,$$

$$I_2=-\langle {ϵ}_2({\Delta }_h{H}^{n+1}-{\Delta }_h{H}^n),H^{n+1}-H^n\rangle ={ϵ}_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_{ϵ}}{2}(\parallel {\nabla }_h{H}^{n+1}{\parallel }^2-\parallel {\nabla }_h{H}^n{\parallel }^2),\hfill \end{array}$$

and

$$I_4=\langle K^{\prime }\left(H^n\right){\nabla }_h{H}^n,H^{n+1}-H^n\rangle \le \frac{{ϵ}_1}{2\tau }\parallel H^{n+1}-H^n{\parallel }^2+\frac{K_1^2\tau }{2{ϵ}_1}{\parallel {\nabla }_h{H}^n\parallel }^2.$$

Summing up, we obtain

$$\begin{array}{c}(C_{ϵ}+\frac{{ϵ}_1}{2\tau })\parallel H^{n+1}-H^n{\parallel }^2+({ϵ}_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{ϵ}_1}{\parallel {\nabla }_h{H}^n\parallel }^2,\hfill \end{array}$$

for ${ϵ}_2\ge K$. Simple calculation shows that

$$\parallel {\nabla }_h{H}^{n+1}{\parallel }^2\le (1+\frac{K_1^2\tau }{{ϵ}_1{K}_s}){\parallel {\nabla }_h{H}^n\parallel }^2,$$

(10)

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$.

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

$$E_1+E_2+F_1\le E_3+F_2+E_5,$$

where

$$E_1=\langle (C\left(H^n\right)+{ϵ}_1)\frac{H^{n+1}-H^n}{\tau },H^{n+1}\rangle \ge \frac{C_{ϵ}+{ϵ}_1}{2\tau }(\parallel H^{n+1}{\parallel }^2-\parallel H^n{\parallel }^2),$$

$$E_2=-\langle {ϵ}_2{\Delta }_h{H}^{n+1},H^{n+1}\rangle ={ϵ}_2{\parallel {\nabla }_h{H}^{n+1}\parallel }^2,$$

$$E_3=-\langle {ϵ}_2{\Delta }_h{H}^n,H^{n+1}\rangle \le \frac{{ϵ}_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}{\underbrace{\langle K\left(H^n\right){\nabla }_h{H}^n,{\nabla }_h{H}^n\rangle }}+\underset{F2}{\underbrace{\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_{ϵ}}{2}\parallel {\nabla }_h{H}^n{\parallel }^2+\frac{K_1^2}{2K_{ϵ}}{\parallel H^{n+1}\parallel }^2.$$

Applying Young’s inequality, we obtain

$$\begin{array}{c}{F}_1\ge K_{ϵ}\parallel {\nabla }_h{H}^n{\parallel }^2,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_{ϵ}+{ϵ}_1}{2\tau }-\frac{K_1^2}{2K_{ϵ}})\parallel H^{n+1}{\parallel }^2+(\frac{K_{ϵ}}{2}-\frac{{ϵ}_2}{2}+\frac{K}{2})\parallel {\nabla }_h{H}^n{\parallel }^2+(\frac{{ϵ}_2}{2}-\frac{K}{2})\parallel {\nabla }_h{H}^{n+1}{\parallel }^2\le \frac{C_{ϵ}+{ϵ}_1}{2\tau }{\parallel H^n\parallel }^2.$$

We conclude that for $\tau <(C_{ϵ}+{ϵ}_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.$$

(11)

Combining the above two results (11) and (10), using Sobolev’s embedding inequality, we can get $|H^{n+1}{|}_{\infty ,h}$ is bounded when $n\to \infty$. It implies the inductive hypothesis holds, completing the proof. □

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

Remark 2.

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

Remark 3.

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

4. 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,

$$\left\{\begin{array}{c}C\left(h\right)\frac{\partial h}{\partial t}-\nabla ·\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.$$

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

$$\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 & -\frac{5546t}{{\left(\frac{51z}{50}-\frac{tz(z-40)}{4T}+\frac{207}{10}\right)}^{237/50}+1175000}\hfill \\ \hfill & -\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 ${ϵ}_1=0$ and ${ϵ}_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 ${ϵ}_2$. The results are listed in

Table 1. Stability comparison with different ${ϵ}_2$ and $\tau$.

		$\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{ϵ}}_2$
${ϵ}_2=0$			Unstable	Unstable	Unstable	Unstable	$5.60\times {10}^{-4}$	$3.26\times {10}^{-5}$
${ϵ}_2=0.0001$			Unstable	Unstable	Unstable	Unstable	$1.38\times {10}^{-4}$	$3.75\times {10}^{-5}$
${ϵ}_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}$
${ϵ}_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}$

. It is evident that the improvement in the stability by use of the extra terms is significant. Moreover, in

Table 2. Accuracy.

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

, we set ${ϵ}_1=0,{ϵ}_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 ${ϵ}_1$ and errors.

5. Conclusions

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.