# An Explicit Stabilized Runge–Kutta–Legendre Super Time-Stepping Scheme for the Richards Equation

Authors: Ramesh Chandra Timsina; Harihar Khanal; Kedar Nath Uprety
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

We solve one-dimensional Kirchhof transformed Richards equation numerically using finite difference method with various time-stepping schemes, forward in time central in space (FTCS), backward in time central in space (BTCS), Crank–Nicolson (CN), and a stabilized Runge–Kutta–Legendre super time-stepping (RKL), and compare their performances.


---

## Body

## 1. Introduction

The prediction of fluid movement in unsaturated porous media is important in many branches of science and engineering such as soil sciences, agricultural engineering, environmental engineering, and ground water hydrology. There are so many problems related to the movement of water in the unsaturated soil. One of the most important problems is the transfer of the different contaminants from ground surface to groundwater through the unsaturated zone. Flow of water through unsaturated porous media is described by Richards equation which is derived by combining the Darcy’s law and mass conservation equation.The Richards equation can be written in different forms withψ (the pressure head), or θ (the volumetric moisture content), as the dependent variable. We consider the Richards equation in the following (mixed) form [1]:(1)∂θ∂t−∇·Kψ∇ψ+∂K∂z=0,where Kψ is the unsaturated hydraulic conductivity and z is the vertical dimension (downward positive).The constitutive relationship betweenθ=θz,t and ψ=ψz,t allows to express the Richards equation (1) with a single dependent variable ψ or θ.The hydraulic conductivityKψ describes the ease with which water can move through pore space and depends on the intrinsic permeability of the material and the properties of fluid such as degree of saturation, density, and viscosity [2]. There are many empirical formulations for the hydraulic conductivity Kψ and the moisture content θψ functions. We use the following popular model from groundwater hydrology due to Haverkamp et al. found at Celia’s paper [3] which describes these constitutive relations as the continuous functions of ψ.(2)Kψ=KsAα+ψγ,θψ=θr+αθs−θrα+ψβ,where θs and θr represent the saturated and residual moisture content, respectively, Ks corresponds to the saturated hydraulic conductivity, and A, α, β, and γ are dimensionless soil parameters.Both the hydraulic conductivity and the moisture content functionsK and θ in (2) are highly nonlinear since they can dramatically change over a small range of ψ. As such, the Richards equation (1) becomes a highly nonlinear partial differential equation and its analytic solutions are limited to some simplified cases with no practical importance. Because of the lack of analytic solution to the Richards equation, study on the computational methods for the unsaturated flow problem has grown significantly. Several numerical schemes based on finite difference, finite element, and finite volume spatial discretization techniques to the partial differential equation have been developed to approximate the solution of the Richards equation [4–9]. Different numerical approaches yield different accuracy. Adaptive time-stepping strategies are studied in [10].Due to the stability restriction of the explicit schemes for the parabolic partial differential equations, it seems there has not been much attention given to develop explicit schemes for the Richards equation and extensive amount of research work is devoted to developing implicit schemes. Due to the highly nonlinear nature of the problem, to implement implicit schemes, no mater which numerical approaches we use, the problem has to be linearized somehow at some stage. As such, some iterative approaches need to be applied to tackle this highly nonlinear problem [3, 11]. The implicit schemes are unconditionally stable, but with the iterative process involved, they all turn out to be computationally expensive and in certain circumstances, unreliable.Recent trends in computational approach are to develop efficient parallel algorithm for high-performing computers. To be able to parallelize the numerical procedures, the algorithm should be iteration free. For this, an explicit scheme is the best option. In [12], a linearized Richards equation model is studied. In [8], stability analysis of an explicit scheme for the Richards equation is studied.In this paper, we use Kirchhoff integral transform to reduce the highly nonlinear equation to a functional linear parabolic equation and solve it numerically.The main purpose of this paper is to implement a relatively new numerical approach (super time-stepping to stabilize an explicit scheme [13] to transformed Richards equation). The work here presented describes and verifies the employment and accuracy of a stabilized Runge–Kutta–Legendre super time-stepping strategy (RKL) to an explicit finite difference scheme to simulate flow in unsaturated porous media. The advantage of using our approach is that it is unconditionally stable, easy to implement, can be easily extended to problems in higher dimensions, and can be easily parallelized.This paper is organized as follows: in Section2, we present Kirchhoff transformation to (spatially) linearize the Richards equation. In Section 3, we present the numerical methods based on finite difference method with different time-stepping schemes. In Section 4, we solve multiple test cases and compare the results. Finally, in Section 5, we present our conclusions.


## 2. Simplified One-Dimensional Model

We consider Richards equation (1) in one space dimension with no source term.(3)∂θ∂t=∂∂zKψ∂ψ∂z−∂K∂z.Richards equation(3) is typically used to simulate infiltration experiments (in both lab and field scale). These experiments begin with a dry soil and then water is poured on top of the ground surface, showing a clear connection with the Darcy’s law. We assume that the infiltration with known surface flux does not exceed the infiltration intensity and does not generate runoff. That is, we use the following initial and boundary conditions:(4)ψz,0=ψ0z,0<z<L,Kψ−Kψ∂ψ∂z=qt,z=0,t>0,ψL,t=βt,t>0.


### 2.1. Kirchhof Integral Transform

To apply Kirchhof integral transformation to equation (3), we let h=ψ−z and define(5)ϕh=∫0hK¯λdλ.SinceKh>0 from (2), the function ϕh is strictly increasing with K¯h=Kψ. Taking derivative of both sides of the transformation, we obtain(6)∂ϕ∂z=∂ϕ∂h∂h∂z=K¯h∂ψ−z∂z=Kψ∂ψ∂z−1=Kψ∂ψ∂z−Kψ.Again, taking derivative of equation (6),(7)∂2ϕ∂z2=∂Kψ∂ψ/∂z∂z−∂∂zKψ.Using equation (7), the Richards equation (3) takes the form(8)∂θ¯∂t=∂2ϕ∂z2,with θ¯ϕ=θh. The corresponding initial and boundary conditions for the transformed equation (8) take the following form:(9)ϕz,0=ϕ0z,0<z<L,∂ϕ∂z=q¯t,z=0,t>0,ϕL,t=β¯t,t>0.The Kirchhoff transformation transformed the doubly nonlinear equation (3) to a nonlinear parabolic problem (8). Also, we note that the Kirchhoff transformation preserves the uniqueness result for the transformed problem.


## 2.1. Kirchhof Integral Transform

To apply Kirchhof integral transformation to equation (3), we let h=ψ−z and define(5)ϕh=∫0hK¯λdλ.SinceKh>0 from (2), the function ϕh is strictly increasing with K¯h=Kψ. Taking derivative of both sides of the transformation, we obtain(6)∂ϕ∂z=∂ϕ∂h∂h∂z=K¯h∂ψ−z∂z=Kψ∂ψ∂z−1=Kψ∂ψ∂z−Kψ.Again, taking derivative of equation (6),(7)∂2ϕ∂z2=∂Kψ∂ψ/∂z∂z−∂∂zKψ.Using equation (7), the Richards equation (3) takes the form(8)∂θ¯∂t=∂2ϕ∂z2,with θ¯ϕ=θh. The corresponding initial and boundary conditions for the transformed equation (8) take the following form:(9)ϕz,0=ϕ0z,0<z<L,∂ϕ∂z=q¯t,z=0,t>0,ϕL,t=β¯t,t>0.The Kirchhoff transformation transformed the doubly nonlinear equation (3) to a nonlinear parabolic problem (8). Also, we note that the Kirchhoff transformation preserves the uniqueness result for the transformed problem.


## 3. Numerical Method

To solve the transformed equation (8) numerically with the prescribed initial and boundary conditions (9), it is more convenient to have a single state variable. For this, we assume θ and ϕ are single valued continuous functions of one another and rearrange to obtain(10)∂θ∂t=∂θ∂ϕ∂ϕ∂t=1∂ϕ/∂θ∂ϕ∂t,∂ϕ∂θ=∂ϕ∂h∂h∂θ.Differentiating (2) and (9) with respect to h, we get(11)∂θ∂h=αθs−θrα+hβ−2⋅βhβ−1,∂ϕ∂h=K¯h=Kψ.Using (10) and (11), the transformed Richards equation (8) takes the form(12)cϕ∂ϕ∂t=∂2ϕ∂z2,where the functional coefficient c depends on ϕ through h as(13)cϕh=αβθs−θrhβ−1K¯hα+hβ2.


### 3.1. Finite Difference Discretization

LetΔz=L/M and Δt=T/N. We construct a grid zj,tn, with zj=jΔz,j=0,1,2,…,M and tn=nΔt,n=0,1,2,…,N. Let ϕjn denote ϕzj,tn. The partial differential equation (12) can be approximated using forward difference in time and central difference in space as(14)∂ϕ∂t|zj,tn≈ϕjn+1−ϕjnΔt,∂2ϕ∂z2|zj,tn≈ϕj−1n−2ϕjn+ϕj+1nΔz2.Let0≤λ≤1. Using a weighted average of the derivative ∂2ϕ/∂z2 at two time levels, tn and tn+1, equation (12) can be discretized as(15)ϕjn+1=ϕjn+σjnλϕj−1n+1−2ϕjn+1+ϕj+1n+1+1−λϕj−1n−2ϕjn+ϕj+1n,where σjn=Δt/cjnΔz2.Equation (15) is used to update the values of ϕjn+1 for the internal nodes. Using the second order central difference with ghost node approach at the upper boundary, we get(16)1+2λσjnϕ0n+1−2λσjnϕ1n+1=1−21−λσ0nϕ0n+21−λσ0nϕ1n−2Δzλσ0n+1qtn+1+1−λσ0nqtn,ϕMn+1=β¯tn+1.The numerical schemes (15) and (16) represent a forward in time central in space (FTCS), backward in time central in space (BTCS), and Crank-Nicolson (CN) schemes for λ=0, λ=1, and λ=1/2, respectively [14]. The error associated with this approximation is OΔz2+Δt for all λ≠1/2. In the case of Crank–Nicolson, it is OΔz2+Δt2, second order accurate in both space and time.We can express the above numerical procedure in a tridiagonal matrix system as(17)1+2λσ0n−2λσ0n0…000−λσ1n1+2λσ1n−λσ1n…000…………….….….…………….….….000…−λσM−1n1+2λσM−1n−λσM−1n000…001ϕ0n+1ϕ1n+1……ϕM−1n+1ϕMn+1=1−21−λσ0nϕ0n+21−λσ0nϕ1n−2Δzλσ0n+1qtn+1+1−λσ0nqtn1−λσ1nϕ0n+1−21−λσ1nϕ1n+1−λσ1nϕ2n………1−λσM−1nϕM−2n+1−21−λσM−1nϕM−1n+1−λσM−1nϕMnβ¯tn+1.The numerical scheme (17) can be used to update the transformed variable ϕjn to its value in the next time level ϕjn+1. But we cannot advance the algorithm to the next time level ϕjn+2 without evaluating the function cϕjn+1 which requires computing the intermediate variable hjn+1. For this, we employ equation (11) which can be approximated as(18)hjn+1=hjn+ϕjn+1−ϕjnK¯hjn.


### 3.2. Von Neumann Stability Analysis

Substituting Fourier componentϕjn=ϕ˜neji2πhξ in equation (15) and rearranging the terms, we get(19)ϕ˜jn+11−λσjne−2iπhξ−2+e2iπhξ=ϕ˜jn1+1−λσjne−2iπhξ−2+e2iπhξ.Utilizing the relations2hπξ=u and cosu=eiu+e−iu/2, equation (19) becomes(20)ϕ˜jn+11+2λσjn1−cosu=ϕ˜jn1−21−λσjn1−cosu.Defining the amplification factorG by(21)Gu=1−21−λσjn1−cosu1+2λσjn1−cosu=1−41−λσjnsin2u/21+4λσjnsin2u/2,we write equation (20) as(22)ϕ˜jn+1=Guϕ˜jn.For a stable solution, the absolute value ofG must be less or equal to 1 for all values of u. Since G′u=0 at u=0,±π, and G0=1, the stability requirement can be expressed as(23)1−41−λσjn1+4λσjn≤1.Since1−4σjn1−λ/1+4λσjn is always ≤1, stability condition (24) reduces to(24)4σjn1−2λ≤2.It is easy to see that ifλ≥1/2, inequality (24) is always satisfied and if λ<1/2, inequality (24) is satisfied only if(25)σjn≤121−2λ.Hence, we see that ifλ≥1/2, our scheme (15) is unconditionally stable and for λ<1/2 scheme is stable only if maxσjn≤1/2 which puts a severe restriction on the time-step size known as Courant–Friedrichs–Lewy (CFL) condition, as(26)Δt≤12mincϕjnΔz2.Note that whencϕ=0 in (12) becomes purely a diffusion equation. Thus, to avoid degeneracy, we need to adjust a small parameter in the calculation of cϕjn.


### 3.3. Super Time-Stepping Scheme

Whenλ=0, equation (17) takes the form(27)ϕ0n+1ϕ1n+1……ϕM−1n+1ϕMn+1=1−2σ0n2σ0n00000σ1n1−2σ1nσ1n…000…………….….….…………….….….000…σM−1n1−2σM−1nσM−1n000…000ϕ0nϕ1n……ϕM−1nϕMn+2σ0nΔzqtn0……0β¯tn+1,which can be expressed in the following form:(28)Φn+1=RΔtAΦn+ΔtFn,where the amplification factor R=I−ΔtA with the tridiagonal matrix A appropriately defined from the above linear system.The explicit scheme (26) is easy to implement as the right hand side of this equation is explicitly known at the time level tn and there is no need to invert the matrix or need to apply any iterative approach to solve it. However, there is a serious drawback of this scheme. It is restricted by the stability criterion, also known as CFL condition (18).As is well known, abovementioned algorithm (26) is subject to the restrictive stability condition which is equivalent to the fact that the spectral radius of matrix R is less than or equal to unity.Super time-stepping scheme relaxes restriction of the CFL condition by requiring stability at the end of one super time-stepΔτ consisting of a cycle of s substeps, rather than at the end of each time-step Δt, thus leading to a Runge–Kutta-like method with s stages. The inner values have, in general, no approximation properties and should only be considered as intermediate calculations.In Runge–Kutta–Chebychev (RKC) scheme [15, 16], the amplification factor R is defined as(29)Rs=∏k=1sI−τkA,and the intermediate substeps τ1,τ2,…,τs in one cycle of the super time-step are chosen as the roots of a modified Chebychev polynomial which maximizes the duration of the super time-step Δτ=∑k=1sτk and ensures stability condition ρRs<1 at the end of one super time-step. Similar to the Chebyshev polynomials, Legendre polynomials are also bounded in magnitude by unity and are useful to develop a stable scheme. Here, we use a stabilized Runge–Kutta–Legendre super time-stepping scheme where the amplification factor RsΔτA is defined in terms of the Legendre polynomials [13, 16].(30)Rsz=as+bsPswo+w1z.The parametersas=0, bs=1, w0=1, and w1=2/s2+s are chosen to satisfy the consistency conditions at first order Rs0=1 and Rs′0=1. Thus, the s-stage RK scheme takes the form(31)Φn+1=PsI−w1ΔτAΦn+ΔτFn.Using the three-point recursion property of the Legendre polynomials,(32)kPkx=2k−1xPk−1x−k−1Pk−2x.RKL scheme (30) can be written as(33)U0=Φn,U1=U0−μ1¯ΔτAU0+ΔτFn,Uk=μkUk−1+νkUk−2−μk¯ΔτAUk−1,2≤k≤s,Φn+1=Us,where(34)μk=2k−1k,μ¯k=2k−1kw1,νk=1−kk.The maximum stable super time-step for the aboves-stage RKL scheme is given by(35)Δτmax=Δtw1=Δts2+s2.Note that unlike RKC, the RKL method described by the recursive relation is consistent and stable at each of the intermediate stages which is more flexible for output.


## 3.1. Finite Difference Discretization

LetΔz=L/M and Δt=T/N. We construct a grid zj,tn, with zj=jΔz,j=0,1,2,…,M and tn=nΔt,n=0,1,2,…,N. Let ϕjn denote ϕzj,tn. The partial differential equation (12) can be approximated using forward difference in time and central difference in space as(14)∂ϕ∂t|zj,tn≈ϕjn+1−ϕjnΔt,∂2ϕ∂z2|zj,tn≈ϕj−1n−2ϕjn+ϕj+1nΔz2.Let0≤λ≤1. Using a weighted average of the derivative ∂2ϕ/∂z2 at two time levels, tn and tn+1, equation (12) can be discretized as(15)ϕjn+1=ϕjn+σjnλϕj−1n+1−2ϕjn+1+ϕj+1n+1+1−λϕj−1n−2ϕjn+ϕj+1n,where σjn=Δt/cjnΔz2.Equation (15) is used to update the values of ϕjn+1 for the internal nodes. Using the second order central difference with ghost node approach at the upper boundary, we get(16)1+2λσjnϕ0n+1−2λσjnϕ1n+1=1−21−λσ0nϕ0n+21−λσ0nϕ1n−2Δzλσ0n+1qtn+1+1−λσ0nqtn,ϕMn+1=β¯tn+1.The numerical schemes (15) and (16) represent a forward in time central in space (FTCS), backward in time central in space (BTCS), and Crank-Nicolson (CN) schemes for λ=0, λ=1, and λ=1/2, respectively [14]. The error associated with this approximation is OΔz2+Δt for all λ≠1/2. In the case of Crank–Nicolson, it is OΔz2+Δt2, second order accurate in both space and time.We can express the above numerical procedure in a tridiagonal matrix system as(17)1+2λσ0n−2λσ0n0…000−λσ1n1+2λσ1n−λσ1n…000…………….….….…………….….….000…−λσM−1n1+2λσM−1n−λσM−1n000…001ϕ0n+1ϕ1n+1……ϕM−1n+1ϕMn+1=1−21−λσ0nϕ0n+21−λσ0nϕ1n−2Δzλσ0n+1qtn+1+1−λσ0nqtn1−λσ1nϕ0n+1−21−λσ1nϕ1n+1−λσ1nϕ2n………1−λσM−1nϕM−2n+1−21−λσM−1nϕM−1n+1−λσM−1nϕMnβ¯tn+1.The numerical scheme (17) can be used to update the transformed variable ϕjn to its value in the next time level ϕjn+1. But we cannot advance the algorithm to the next time level ϕjn+2 without evaluating the function cϕjn+1 which requires computing the intermediate variable hjn+1. For this, we employ equation (11) which can be approximated as(18)hjn+1=hjn+ϕjn+1−ϕjnK¯hjn.


## 3.2. Von Neumann Stability Analysis

Substituting Fourier componentϕjn=ϕ˜neji2πhξ in equation (15) and rearranging the terms, we get(19)ϕ˜jn+11−λσjne−2iπhξ−2+e2iπhξ=ϕ˜jn1+1−λσjne−2iπhξ−2+e2iπhξ.Utilizing the relations2hπξ=u and cosu=eiu+e−iu/2, equation (19) becomes(20)ϕ˜jn+11+2λσjn1−cosu=ϕ˜jn1−21−λσjn1−cosu.Defining the amplification factorG by(21)Gu=1−21−λσjn1−cosu1+2λσjn1−cosu=1−41−λσjnsin2u/21+4λσjnsin2u/2,we write equation (20) as(22)ϕ˜jn+1=Guϕ˜jn.For a stable solution, the absolute value ofG must be less or equal to 1 for all values of u. Since G′u=0 at u=0,±π, and G0=1, the stability requirement can be expressed as(23)1−41−λσjn1+4λσjn≤1.Since1−4σjn1−λ/1+4λσjn is always ≤1, stability condition (24) reduces to(24)4σjn1−2λ≤2.It is easy to see that ifλ≥1/2, inequality (24) is always satisfied and if λ<1/2, inequality (24) is satisfied only if(25)σjn≤121−2λ.Hence, we see that ifλ≥1/2, our scheme (15) is unconditionally stable and for λ<1/2 scheme is stable only if maxσjn≤1/2 which puts a severe restriction on the time-step size known as Courant–Friedrichs–Lewy (CFL) condition, as(26)Δt≤12mincϕjnΔz2.Note that whencϕ=0 in (12) becomes purely a diffusion equation. Thus, to avoid degeneracy, we need to adjust a small parameter in the calculation of cϕjn.


## 3.3. Super Time-Stepping Scheme

Whenλ=0, equation (17) takes the form(27)ϕ0n+1ϕ1n+1……ϕM−1n+1ϕMn+1=1−2σ0n2σ0n00000σ1n1−2σ1nσ1n…000…………….….….…………….….….000…σM−1n1−2σM−1nσM−1n000…000ϕ0nϕ1n……ϕM−1nϕMn+2σ0nΔzqtn0……0β¯tn+1,which can be expressed in the following form:(28)Φn+1=RΔtAΦn+ΔtFn,where the amplification factor R=I−ΔtA with the tridiagonal matrix A appropriately defined from the above linear system.The explicit scheme (26) is easy to implement as the right hand side of this equation is explicitly known at the time level tn and there is no need to invert the matrix or need to apply any iterative approach to solve it. However, there is a serious drawback of this scheme. It is restricted by the stability criterion, also known as CFL condition (18).As is well known, abovementioned algorithm (26) is subject to the restrictive stability condition which is equivalent to the fact that the spectral radius of matrix R is less than or equal to unity.Super time-stepping scheme relaxes restriction of the CFL condition by requiring stability at the end of one super time-stepΔτ consisting of a cycle of s substeps, rather than at the end of each time-step Δt, thus leading to a Runge–Kutta-like method with s stages. The inner values have, in general, no approximation properties and should only be considered as intermediate calculations.In Runge–Kutta–Chebychev (RKC) scheme [15, 16], the amplification factor R is defined as(29)Rs=∏k=1sI−τkA,and the intermediate substeps τ1,τ2,…,τs in one cycle of the super time-step are chosen as the roots of a modified Chebychev polynomial which maximizes the duration of the super time-step Δτ=∑k=1sτk and ensures stability condition ρRs<1 at the end of one super time-step. Similar to the Chebyshev polynomials, Legendre polynomials are also bounded in magnitude by unity and are useful to develop a stable scheme. Here, we use a stabilized Runge–Kutta–Legendre super time-stepping scheme where the amplification factor RsΔτA is defined in terms of the Legendre polynomials [13, 16].(30)Rsz=as+bsPswo+w1z.The parametersas=0, bs=1, w0=1, and w1=2/s2+s are chosen to satisfy the consistency conditions at first order Rs0=1 and Rs′0=1. Thus, the s-stage RK scheme takes the form(31)Φn+1=PsI−w1ΔτAΦn+ΔτFn.Using the three-point recursion property of the Legendre polynomials,(32)kPkx=2k−1xPk−1x−k−1Pk−2x.RKL scheme (30) can be written as(33)U0=Φn,U1=U0−μ1¯ΔτAU0+ΔτFn,Uk=μkUk−1+νkUk−2−μk¯ΔτAUk−1,2≤k≤s,Φn+1=Us,where(34)μk=2k−1k,μ¯k=2k−1kw1,νk=1−kk.The maximum stable super time-step for the aboves-stage RKL scheme is given by(35)Δτmax=Δtw1=Δts2+s2.Note that unlike RKC, the RKL method described by the recursive relation is consistent and stable at each of the intermediate stages which is more flexible for output.


## 4. Simulation Results

### 4.1. Numerical Setup

The numerical procedure developed in the previous section is written in Python and ran on a laptop with 2.8 GHz Quad-Core Intel Core i7 processor. We inspect the behavior of the four numerical schemes presented in the previous section in a specific infiltration experiment. In this simulation, we consider a vertical soil column of depthL=70 cm in a time period of tmax=1 hr.Following Haverkamp et al. [7], we use the soil parameters and the characteristics relationship between the soil moisture content θψ and the hydraulic conductivity Kψ as follows:(36)θψ=θr+αθs−θrα+ψβ=0.075+1.616×106.287−.0751.616×106+ψ3.96,Kψ=KsAA+ψγ=34×1.175×1061.175×106+ψ4.74.The simulation starts with a uniform saturationθ=0.1 cm3/cm3 and a constant water head ψ=−61.5 cm is maintained at the bottom boundary z=L. At the upper boundary z=0 (the soil surface), a constant flux qt=13.69 cm/hr for t<0.3 hr and zero normal flux condition for t>0.3 hr.To compute the approximate solution using the explicit scheme FTCS, we have used a uniform spatial step sizeΔz=1 cm and the time-step size Δt=0.005 sec which guarantees CFL condition (26). This being highly resolved in temporal direction, we use it as the representative solution to compare with all other methods. The implicit schemes BTCS and Crank–Nicolson (CN), being unconditionally stable, have no restriction on the time-step size. The step size is limited only by the accuracy of approximation. By construction, the stabilized Runge–Kutta–Legendre super time-stepping (RKL) scheme is also unconditionally stable. The step size in RKL is determined by the choice of number of substeps in one super time-stepping cycle and, of course, this can be chosen by compromising the accuracy of approximation. As such, to analyze the performance of RKL compared to the implicit schemes BTCS and CN, we run BTCS and CN with Δt=∑k=1sΔτ for different values of s. To deal with the nonlinear dependence of the functional coefficient cϕ in the implicit schemes, we use fix point iteration with maximum allowable iteration MAXIT=10 and relative error tolerance TOL=10−6.


### 4.2. Results and Discussion

Richards equation is a highly nonlinear degenerate partial differential equation. The nonlinear behavior appears from the use of constitutive relationship betweenθ and ψ and K and ψ.The relative curves to the experimental data with appropriate parameter values in the moisture retention function and the hydraulic conductivity model as found in the study by Haverkamp et al. [7] are shown in Figure 1.Figure 1
The shape of functionθψ (a) and Kψ (b).
(a)(b)The numerical solutions of the Richards equation are computed using four different schemes, namely, FTCS, BTCS, CN, and RKL. First, we observed that FTCS scheme is conditionally stable, whereas BTCS, CN, and RKL schemes are unconditionally stable. Thus, to make a long simulation with the explicit scheme, we need to fix the CFL-satisfying time-step size. Since we do not have exact solution to estimate the errors, we use the numerical solution obtained from the fully explicit numerical scheme with a fine mesh (Δz=1 cm and Δt=0.002 sec) as the reference solution.For the performance comparison, we make simulations with a fairly fine mesh ofΔz=2 cm and find Δtexpl using CFL criteria and define Δτ=Δtexpls2+s/2 as the time-step size of a substep in a super time-step. Note that, in general, comparison of explicit and implicit scheme is not realistic. Explicit scheme requires too many time steps, and implicit scheme needs too many iterations. Since we are using linear system solver in the implicit schemes, it is impossible to compare the efficiency using the same time-step size. To make the comparison more realistic, we have used the time-step size for implicit schemes equal to the duration of one super-step in RKL, Δt=sΔτ, the duration of one super time-step size. The simulation results (speedup and the accuracy) are displayed in Tables 1–3 and Figures 2–4. The CPU time (in sec) of all the numerical schemes implemented in this work is shown in Table 1. The CPU timings of the RKL scheme for different values of s (the number of substeps in one super time-step) show the acceleration in the speed of the explicit scheme FTCS (s=1). In Figure 2, the comparison of the CPU timings of the RKL for various values of s with the implicit schemes BTCS and CN is plotted. Tables 2 and 3 show comparison of errors in the final values of moisture content θ at t = 1 hr obtained from RKL and the other three schemes in L∞-norm and L2-norm. Table 2 also shows the stability of the explicit scheme RKL as the overall error in L∞-norm with respect to the reference solution which is not increasing much with the increase in the values of s. Similar behavior can also be observed from Table 3 with the errors measured in L2-norm.Table 1
Timing of four methods.
sEXPLRKLIMPLCN533.1415.40101.23101.251032.668.2628.1828.342032.694.798.418.292533.044.315.926.053032.803.694.444.553532.763.273.603.624032.803.173.083.124532.752.992.672.855032.602.842.412.535532.682.642.222.346032.832.442.082.006532.872.481.892.058032.812.211.681.5310032.772.451.571.71Table 2
Test of RKL scheme with EXPL, IMPL, and CN inL∞.
sCPU time (sec)RKL−EE∞RKL−IE∞RKL−CN∞546.22.25192E−022.24732E−022.24649E−021024.781.37265E−022.24732E−021.34624E−022014.377.33902E−036.61120E−036.48773E−032512.935.80247E−034.63486E−034.43717E−033011.074.20458E−031.93189E−031.63057E−03359.813.27050E−035.05055E−041.19139E−04409.512.53303E−031.19531E−031.78467E−03458.971.61003E−034.21368E−034.69720E−03508.523.80157E−056.90071E−039.08979E−03557.925.36361E−047.11451E−037.43036E−03607.322.17285E−031.13744E−021.53756E−02657.443.23959E−031.15905E−021.69259E−02806.339.77096E−032.24249E−025.21857E−021007.357.96584E−032.15554E−024.52289E−02Table 3
Test of RKL scheme with EXPL, IMPL, and CN inL2.
sCPU time (sec)RKL−EE2RKL−IE2RKL−CN2546.23.23587E−023.22906E−023.22784E−021024.782.02935E−021.99448E−021.98998E−022014.371.10585E−029.97846E−039.79713E−032512.938.77606E−037.04349E−036.75439E−033011.076.38021E−032.99619E−032.55873E−03359.814.97194E−038.79174E−042.55407E−04409.513.84175E−031.69154E−032.53512E−03458.972.43071E−036.18330E−036.84142E−03508.524.63038E−051.01977E−021.32881E−02557.928.80764E−041.05394E−021.08956E−02607.323.32750E−031.67543E−022.22226E−02657.445.05941E−031.72343E−022.45424E−02806.331.35737E−023.17820E−027.01856E−021007.351.27567E−023.22055E−026.30587E−02Figure 2
Convergence of the methods (EXPL, RKL, IMPL, and CN):s = 20 in (a), s = 100 in (b), s = 50 in (c), and s = 35 in (d).
(a)(b)(c)(d)Figure 3
Test comparison of the four methods’ error in infinity norm (a), in 2 norms (b), and speedup with RKL (c).
(a)(b)(c)Figure 4
Variational trend of moisture content in depth (a) and pressure head (b).
(a)(b)Comparison of the solutionθsRKLz,t at t=0.3hr with different values of s to fully explicit θFTCSz,.3, implicit θBTCSz,.3, and Crank–Nicolson θCNz,.3 methods is shown in Figure 3. These values of s are selected randomly to show that as s becomes large, there will be less accurate approximation and as s⟶1, the method becomes the most accurate. From Figure 3, we can see that s=35 gives the best performance. It should be noted that RKL can be run with higher s by compromising with the low accuracy. Finally, we present some numerical simulation results in a very fine spatial mesh Δz=0.5 cm obtained from the RKL method s=100 and CN Δt=sΔτ . In Figures 4 and 5, we describe the profile of moisture content θz,t and pressure head ψz,t along the soil depth at various times t=0.1,0.2,…,1 hr.Figure 5
Variational trend of moisture content in depth and time (a) and pressure head (b).
(a)(b)


## 4.1. Numerical Setup

The numerical procedure developed in the previous section is written in Python and ran on a laptop with 2.8 GHz Quad-Core Intel Core i7 processor. We inspect the behavior of the four numerical schemes presented in the previous section in a specific infiltration experiment. In this simulation, we consider a vertical soil column of depthL=70 cm in a time period of tmax=1 hr.Following Haverkamp et al. [7], we use the soil parameters and the characteristics relationship between the soil moisture content θψ and the hydraulic conductivity Kψ as follows:(36)θψ=θr+αθs−θrα+ψβ=0.075+1.616×106.287−.0751.616×106+ψ3.96,Kψ=KsAA+ψγ=34×1.175×1061.175×106+ψ4.74.The simulation starts with a uniform saturationθ=0.1 cm3/cm3 and a constant water head ψ=−61.5 cm is maintained at the bottom boundary z=L. At the upper boundary z=0 (the soil surface), a constant flux qt=13.69 cm/hr for t<0.3 hr and zero normal flux condition for t>0.3 hr.To compute the approximate solution using the explicit scheme FTCS, we have used a uniform spatial step sizeΔz=1 cm and the time-step size Δt=0.005 sec which guarantees CFL condition (26). This being highly resolved in temporal direction, we use it as the representative solution to compare with all other methods. The implicit schemes BTCS and Crank–Nicolson (CN), being unconditionally stable, have no restriction on the time-step size. The step size is limited only by the accuracy of approximation. By construction, the stabilized Runge–Kutta–Legendre super time-stepping (RKL) scheme is also unconditionally stable. The step size in RKL is determined by the choice of number of substeps in one super time-stepping cycle and, of course, this can be chosen by compromising the accuracy of approximation. As such, to analyze the performance of RKL compared to the implicit schemes BTCS and CN, we run BTCS and CN with Δt=∑k=1sΔτ for different values of s. To deal with the nonlinear dependence of the functional coefficient cϕ in the implicit schemes, we use fix point iteration with maximum allowable iteration MAXIT=10 and relative error tolerance TOL=10−6.


## 4.2. Results and Discussion

Richards equation is a highly nonlinear degenerate partial differential equation. The nonlinear behavior appears from the use of constitutive relationship betweenθ and ψ and K and ψ.The relative curves to the experimental data with appropriate parameter values in the moisture retention function and the hydraulic conductivity model as found in the study by Haverkamp et al. [7] are shown in Figure 1.Figure 1
The shape of functionθψ (a) and Kψ (b).
(a)(b)The numerical solutions of the Richards equation are computed using four different schemes, namely, FTCS, BTCS, CN, and RKL. First, we observed that FTCS scheme is conditionally stable, whereas BTCS, CN, and RKL schemes are unconditionally stable. Thus, to make a long simulation with the explicit scheme, we need to fix the CFL-satisfying time-step size. Since we do not have exact solution to estimate the errors, we use the numerical solution obtained from the fully explicit numerical scheme with a fine mesh (Δz=1 cm and Δt=0.002 sec) as the reference solution.For the performance comparison, we make simulations with a fairly fine mesh ofΔz=2 cm and find Δtexpl using CFL criteria and define Δτ=Δtexpls2+s/2 as the time-step size of a substep in a super time-step. Note that, in general, comparison of explicit and implicit scheme is not realistic. Explicit scheme requires too many time steps, and implicit scheme needs too many iterations. Since we are using linear system solver in the implicit schemes, it is impossible to compare the efficiency using the same time-step size. To make the comparison more realistic, we have used the time-step size for implicit schemes equal to the duration of one super-step in RKL, Δt=sΔτ, the duration of one super time-step size. The simulation results (speedup and the accuracy) are displayed in Tables 1–3 and Figures 2–4. The CPU time (in sec) of all the numerical schemes implemented in this work is shown in Table 1. The CPU timings of the RKL scheme for different values of s (the number of substeps in one super time-step) show the acceleration in the speed of the explicit scheme FTCS (s=1). In Figure 2, the comparison of the CPU timings of the RKL for various values of s with the implicit schemes BTCS and CN is plotted. Tables 2 and 3 show comparison of errors in the final values of moisture content θ at t = 1 hr obtained from RKL and the other three schemes in L∞-norm and L2-norm. Table 2 also shows the stability of the explicit scheme RKL as the overall error in L∞-norm with respect to the reference solution which is not increasing much with the increase in the values of s. Similar behavior can also be observed from Table 3 with the errors measured in L2-norm.Table 1
Timing of four methods.
sEXPLRKLIMPLCN533.1415.40101.23101.251032.668.2628.1828.342032.694.798.418.292533.044.315.926.053032.803.694.444.553532.763.273.603.624032.803.173.083.124532.752.992.672.855032.602.842.412.535532.682.642.222.346032.832.442.082.006532.872.481.892.058032.812.211.681.5310032.772.451.571.71Table 2
Test of RKL scheme with EXPL, IMPL, and CN inL∞.
sCPU time (sec)RKL−EE∞RKL−IE∞RKL−CN∞546.22.25192E−022.24732E−022.24649E−021024.781.37265E−022.24732E−021.34624E−022014.377.33902E−036.61120E−036.48773E−032512.935.80247E−034.63486E−034.43717E−033011.074.20458E−031.93189E−031.63057E−03359.813.27050E−035.05055E−041.19139E−04409.512.53303E−031.19531E−031.78467E−03458.971.61003E−034.21368E−034.69720E−03508.523.80157E−056.90071E−039.08979E−03557.925.36361E−047.11451E−037.43036E−03607.322.17285E−031.13744E−021.53756E−02657.443.23959E−031.15905E−021.69259E−02806.339.77096E−032.24249E−025.21857E−021007.357.96584E−032.15554E−024.52289E−02Table 3
Test of RKL scheme with EXPL, IMPL, and CN inL2.
sCPU time (sec)RKL−EE2RKL−IE2RKL−CN2546.23.23587E−023.22906E−023.22784E−021024.782.02935E−021.99448E−021.98998E−022014.371.10585E−029.97846E−039.79713E−032512.938.77606E−037.04349E−036.75439E−033011.076.38021E−032.99619E−032.55873E−03359.814.97194E−038.79174E−042.55407E−04409.513.84175E−031.69154E−032.53512E−03458.972.43071E−036.18330E−036.84142E−03508.524.63038E−051.01977E−021.32881E−02557.928.80764E−041.05394E−021.08956E−02607.323.32750E−031.67543E−022.22226E−02657.445.05941E−031.72343E−022.45424E−02806.331.35737E−023.17820E−027.01856E−021007.351.27567E−023.22055E−026.30587E−02Figure 2
Convergence of the methods (EXPL, RKL, IMPL, and CN):s = 20 in (a), s = 100 in (b), s = 50 in (c), and s = 35 in (d).
(a)(b)(c)(d)Figure 3
Test comparison of the four methods’ error in infinity norm (a), in 2 norms (b), and speedup with RKL (c).
(a)(b)(c)Figure 4
Variational trend of moisture content in depth (a) and pressure head (b).
(a)(b)Comparison of the solutionθsRKLz,t at t=0.3hr with different values of s to fully explicit θFTCSz,.3, implicit θBTCSz,.3, and Crank–Nicolson θCNz,.3 methods is shown in Figure 3. These values of s are selected randomly to show that as s becomes large, there will be less accurate approximation and as s⟶1, the method becomes the most accurate. From Figure 3, we can see that s=35 gives the best performance. It should be noted that RKL can be run with higher s by compromising with the low accuracy. Finally, we present some numerical simulation results in a very fine spatial mesh Δz=0.5 cm obtained from the RKL method s=100 and CN Δt=sΔτ . In Figures 4 and 5, we describe the profile of moisture content θz,t and pressure head ψz,t along the soil depth at various times t=0.1,0.2,…,1 hr.Figure 5
Variational trend of moisture content in depth and time (a) and pressure head (b).
(a)(b)


## 5. Conclusion

In this work, we considered one-dimensional Richards equation (a highly nonlinear degenerate parabolic partial differential equation) and solved it numerically. We implemented various finite difference schemes FTCS, BTCS, CN, and a stabilized Runge–Kutta–Legendre super time-stepping scheme RKL. The numerical simulations show that the RKL scheme boasts large efficiency gains compared to the standard FTCS scheme and is comparable to CN. Explicit schemes are simple and accurate to implement and convenient for parallelization but suffer severely from stability restriction on the time-step size. In another side, implicit schemes involve iterative approach to solve the nonlinear systems, hence are moderately efficient, can be computationally very costly for highly nonlinear problems, and comparatively difficult to implement. From computational point of view, super time-stepping is more efficient than the standard implicit schemes, in that it runs at least as fast with better accuracy and it is much easier to program; also, it can be easily extended to problems in higher dimensions.



---

*Source: 5573913-2021-04-17.xml*