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

'e 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. 'ere 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. 'e 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]:


Introduction
e 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.
ere 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.
e 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]: where K(ψ) is the unsaturated hydraulic conductivity and z is the vertical dimension (downward positive). e constitutive relationship between θ � θ(z, t) and ψ � ψ(z, t) allows to express the Richards equation (1) with a single dependent variable ψ or θ. e hydraulic conductivity K(ψ) 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]. ere 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 ψ.
where θ s and θ r represent the saturated and residual moisture content, respectively, K s corresponds to the saturated hydraulic conductivity, and A, α, β, and c are dimensionless soil parameters.
Both the hydraulic conductivity and the moisture content functions K 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][5][6][7][8][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]. e 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. e 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). e 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. e 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.
is paper is organized as follows: in Section 2, 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.

Simplified One-Dimensional Model
We consider Richards equation (1) in one space dimension with no source term.
Richards equation (3) is typically used to simulate infiltration experiments (in both lab and field scale). ese 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. at is, we use the following initial and boundary conditions:

Kirchhof Integral Transform.
To apply Kirchhof integral transformation to equation (3), we let h � ψ − z and define Since K(h) > 0 from (2), the function ϕ(h) is strictly increasing with K(h) � K(ψ). Taking derivative of both sides of the transformation, we obtain Again, taking derivative of equation (6), Using equation (7), the Richards equation (3) takes the form with θ(ϕ) � θ(h). e corresponding initial and boundary conditions for the transformed equation (8) e 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.

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 Differentiating (2) and (9) with respect to h, we get Using (10) and (11), the transformed Richards equation (8) takes the form where the functional coefficient c depends on ϕ through h as 3.1. Finite Difference Discretization. Let Δz � L/M and Δt � T/N. We construct a grid (z j , t n ), with z j � jΔz, j � 0, 1, 2, . . . , M and t n � nΔt, n � 0, 1, 2, . . . , N. Let ϕ n j denote ϕ(z j , t n ). e partial differential equation (12) can be approximated using forward difference in time and central difference in space as Let 0 ≤ λ ≤ 1. Using a weighted average of the derivative (z 2 ϕ/zz 2 ) at two time levels, t n and t n+1 , equation (12) can be discretized as where σ n j � Δt/(c n j Δz 2 ). Equation (15) is used to update the values of ϕ n+1 j for the internal nodes. Using the second order central difference with ghost node approach at the upper boundary, we get e 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]. e error associated with this approximation is O(Δz 2 + Δt) for all λ ≠ (1/2). In the case of Crank-Nicolson, it is O(Δz 2 + Δt 2 ), second order accurate in both space and time.
We can express the above numerical procedure in a tridiagonal matrix system as e numerical scheme (17) can be used to update the transformed variable ϕ n j to its value in the next time level ϕ n+1 j . But we cannot advance the algorithm to the next time level ϕ n+2 j without evaluating the function c(ϕ n+1 j ) which requires computing the intermediate variable h n+1 j . For this, we employ equation (11) which can be approximated as

Von Neumann Stability Analysis.
Substituting Fourier component ϕ n j � ϕ n e ji2πhξ in equation (15) and rearranging the terms, we get Utilizing the relations 2hπξ � u and cos u � (e iu + e − iu /2), equation (19) becomes Defining the amplification factor G by we write equation (20) as For a stable solution, the absolute value of G must be less or equal to 1 for all values of u. Since G ′ (u) � 0 at u � 0, ± π, and G(0) � 1, the stability requirement can be expressed as Since (1 − 4σ n j (1 − λ)/1 + 4λσ n j ) is always ≤1, stability condition (24) reduces to It is easy to see that if λ ≥ (1/2), inequality (24) is always satisfied and if λ < (1/2), inequality (24) is satisfied only if Hence, we see that if λ ≥ (1/2), our scheme (15) is unconditionally stable and for λ < (1/2) scheme is stable only if max σ n j ≤ (1/2) which puts a severe restriction on the timestep size known as Courant-Friedrichs-Lewy (CFL) condition, as Note that when c(ϕ) � 0 in (12) becomes purely a diffusion equation.
us, to avoid degeneracy, we need to adjust a small parameter in the calculation of c(ϕ n j ).

Super Time-Stepping Scheme.
When λ � 0, equation (17) takes the form which can be expressed in the following form: where the amplification factor R � I − ΔtA with the tridiagonal matrix A appropriately defined from the above linear system. e explicit scheme (26) is easy to implement as the right hand side of this equation is explicitly known at the time level t n 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-Kuttalike method with s stages. e 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 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 Δτ � s k�1 τ k and ensures stability condition (ρ(R s ) < 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 R s (ΔτA) is defined in terms of the Legendre polynomials [13,16].
us, the s-stage RK scheme takes the form Φ n+1 � P s I − w 1 ΔτA Φ n + ΔτF n . (31) Using the three-point recursion property of the Legendre polynomials, RKL scheme (30) can be written as e maximum stable super time-step for the above s-stage RKL scheme is given by 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.

Numerical Setup.
e 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 Mathematical Problems in Engineering 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 depth L � 70 cm in a time period of t max � 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: 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). is being highly resolved in temporal direction, we use it as the representative solution to compare with all other methods. e implicit schemes BTCS and Crank-Nicolson (CN), being unconditionally stable, have no restriction on the time-step size. e step size is limited only by the accuracy of approximation. By construction, the stabilized Runge-Kutta-Legendre super timestepping (RKL) scheme is also unconditionally stable. e 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 � s k�1 Δτ 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 .

Results and Discussion.
Richards equation is a highly nonlinear degenerate partial differential equation.
e nonlinear behavior appears from the use of constitutive relationship between θ and ψ and K and ψ.
e 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. e 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. us, 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 Δt expl using CFL criteria and define Δτ � Δt expl (s 2 + 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 timestep size for implicit schemes equal to the duration of one super-step in RKL, Δt � sΔτ, the duration of one super timestep size. e simulation results (speedup and the accuracy) are displayed in Tables 1-3 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 Figure 3. ese 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.

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 timestepping scheme RKL. e 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.
Data Availability e data used for supporting the findings of this study are included within the article.