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.
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]:where is the unsaturated hydraulic conductivity and is the vertical dimension (downward positive).
The constitutive relationship between and allows to express the Richards equation (1) with a single dependent variable or .
The hydraulic conductivity 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 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 and represent the saturated and residual moisture content, respectively, corresponds to the saturated hydraulic conductivity, and , , , and are dimensionless soil parameters.
Both the hydraulic conductivity and the moisture content functions 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 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.
2. 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). 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:
2.1. Kirchhof Integral Transform
To apply Kirchhof integral transformation to equation (3), we let and define
Since from (2), the function is strictly increasing with . 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 formwith . The corresponding initial and boundary conditions for the transformed equation (8) take the following form:
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
Differentiating (2) and (9) with respect to , we get
Using (10) and (11), the transformed Richards equation (8) takes the formwhere the functional coefficient depends on through as
3.1. Finite Difference Discretization
Let and . We construct a grid , with and . Let denote . The partial differential equation (12) can be approximated using forward difference in time and central difference in space as
Let . Using a weighted average of the derivative at two time levels, and , equation (12) can be discretized aswhere .
Equation (15) is used to update the values of for the internal nodes. Using the second order central difference with ghost node approach at the upper boundary, we get
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 , , and , respectively [14]. The error associated with this approximation is for all . In the case of Crank–Nicolson, it is , second order accurate in both space and time.
We can express the above numerical procedure in a tridiagonal matrix system as
The numerical scheme (17) can be used to update the transformed variable to its value in the next time level . But we cannot advance the algorithm to the next time level without evaluating the function which requires computing the intermediate variable . For this, we employ equation (11) which can be approximated as
3.2. Von Neumann Stability Analysis
Substituting Fourier component in equation (15) and rearranging the terms, we get
Utilizing the relations and , equation (19) becomes
Defining the amplification factor bywe write equation (20) as
For a stable solution, the absolute value of must be less or equal to 1 for all values of . Since at , and , the stability requirement can be expressed as
Since is always , stability condition (24) reduces to
It is easy to see that if , inequality (24) is always satisfied and if , inequality (24) is satisfied only if
Hence, we see that if , our scheme (15) is unconditionally stable and for scheme is stable only if which puts a severe restriction on the time-step size known as Courant–Friedrichs–Lewy (CFL) condition, as
Note that when in (12) becomes purely a diffusion equation. Thus, to avoid degeneracy, we need to adjust a small parameter in the calculation of .
3.3. Super Time-Stepping Scheme
When , equation (17) takes the formwhich can be expressed in the following form:where the amplification factor with the tridiagonal matrix 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 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 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 substeps, rather than at the end of each time-step , thus leading to a Runge–Kutta-like method with 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 is defined asand the intermediate substeps 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 and ensures stability condition 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 is defined in terms of the Legendre polynomials [13, 16].
The parameters , , , and are chosen to satisfy the consistency conditions at first order and . Thus, the s-stage RK scheme takes the form
Using the three-point recursion property of the Legendre polynomials,
RKL scheme (30) can be written aswhere
The maximum stable super time-step for the above -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.
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 depth cm in a time period of hr.
Following Haverkamp et al. [7], we use the soil parameters and the characteristics relationship between the soil moisture content and the hydraulic conductivity as follows:
The simulation starts with a uniform saturation cm/cm and a constant water head cm is maintained at the bottom boundary . At the upper boundary (the soil surface), a constant flux cm/hr for hr and zero normal flux condition for hr.
To compute the approximate solution using the explicit scheme FTCS, we have used a uniform spatial step size cm and the time-step size 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 for different values of . To deal with the nonlinear dependence of the functional coefficient in the implicit schemes, we use fix point iteration with maximum allowable iteration and relative error tolerance .
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 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.
(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 ( cm and sec) as the reference solution.
For the performance comparison, we make simulations with a fairly fine mesh of cm and find using CFL criteria and define 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, , 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 (the number of substeps in one super time-step) show the acceleration in the speed of the explicit scheme FTCS (). In Figure 2, the comparison of the CPU timings of the RKL for various values of 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 = 1 hr obtained from RKL and the other three schemes in -norm and -norm. Table 2 also shows the stability of the explicit scheme RKL as the overall error in -norm with respect to the reference solution which is not increasing much with the increase in the values of . Similar behavior can also be observed from Table 3 with the errors measured in -norm.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(a)
(b)
Comparison of the solution at with different values of to fully explicit , implicit , and Crank–Nicolson methods is shown in Figure 3. These values of are selected randomly to show that as becomes large, there will be less accurate approximation and as , the method becomes the most accurate. From Figure 3, we can see that gives the best performance. It should be noted that RKL can be run with higher by compromising with the low accuracy. Finally, we present some numerical simulation results in a very fine spatial mesh cm obtained from the RKL method and CN . In Figures 4 and 5, we describe the profile of moisture content and pressure head along the soil depth at various times hr.
(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.
Data Availability
The data used for supporting the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest.