Journal of Hydrology 580 (2020) 124240

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

A control volume scheme using compact integrated radial basis function
stencils for solving the Richards equationT

Duc Ngo-Conga, , Nam Mai-Duya,b, Diogenes L. Antillec,d, Martinus Th. van Genuchtene,f⁎aUniversity of Southern Queensland, Institute for Advanced Engineering and Space Sciences, Toowoomba, QLD 4350, Australia
University of Southern Queensland, School of Mechanical and Electrical Engineering, Faculty of Health, Engineering and Sciences, Toowoomba, QLD 4350, Australiac
CSIRO Agriculture and Food, Canberra, ACT 2601, AustraliadUniversity of Southern Queensland, Centre for Agricultural Engineering, Toowoomba, QLD 4350, AustraliaeFederal University of Rio de Janeiro, Department of Nuclear Engineering, Rio de Janeiro, BrazilfSão Paulo State University, Center for Environmental Studies, CEA, Rio Claro, Brazilb

ARTICLE INFOABSTRACTThis manuscript was handled by Corrado
Corradini, Editor-in-Chief, with the assistanceof Zhiming Lu, Associate Editor.

A new control volume approach is developed based on compact integrated radial basis function (CIRBF) stencils
for solution of the highly nonlinear Richards equation describing transient water flow in variably saturated soils.
Unlike the conventional control volume method, which is regarded as second-order accurate, the proposed
approach has high-order accuracy owing to the use of a compact integrated radial basis function approximationthat enables improved flux predictions. The method is used to solve the Richards equation for transient flow in
1D homogeneous and heterogeneous soil profiles. Numerical results for different boundary conditions, initialconditions and soil types are shown to be in good agreement with Warrick’s semi-analytical solution and simulations using the HYDRUS-1D software package. Results obtained with the proposed method were far lessdependent upon the grid spacing than the HYDRUS-1D finite element solutions.

Keywords:
Richards equation
Finite volume method
Integrated radial basis function
Compact stencil
Unsaturated flow

1. Introduction
The complex problem of transient water flow in unsaturated orvariably-saturated porous media has been widely modelled using thehighly nonlinear Richards equation. Neglecting sink/source terms (e.g,root water uptake) and soil and water compressibilities, the Richards
equation in its mixed form can be written as (Celia et al., 1990):(h )t·[K (h) h]

K (h)
=0
z

(1)

where θ and K are, respectively, the volumetric water content andhydraulic conductivity, both being a function of the pressure head h, isthe vertical coordinate (assumed positive upward here), and t is time. It
is inherently difficult to obtain numerical solutions of Eq. (1) due to thehighly nonlinear nature of the soil hydraulic functions θ(h) and K(h).
Several approaches have been used to deal with nonlinearity of the
Richards equation, such as applying Picard iteration (e.g., Celia et al.,
1990; Huang et al., 1996) and Newton Raphson (e.g., Lehmann andAckerer, 1998) methods, the method of lines (e.g., Fahs et al., 2009), or
certain transformations (Kirkland et al., 1992; Pan and Wierenga,
1995). The nonlinearities pose especially difficulties under some⁎extreme conditions (e.g., sharp wetting fronts, very dry initial conditions, extreme soil textures, heterogeneous and/or hysteretic soils)
using traditional analytical or numerical techniques (Farthing andOgden, 2017). For example, finite element solutions often suffer fromnumerical oscillations near sharp wetting fronts unless very dense grids
are used (Celia et al., 1990; Pan and Wierenga, 1995; van Genuchten,1982).
In the past few years, many numerical schemes have been developed
for solving the Richards equation (e.g., Fahs et al., 2009; Farthing et al.,
2003; Kavetski et al., 2001; Lehmann and Ackerer, 1998; McBride et al.,2006; Zha et al., 2019). For example, Lehmann and Ackerer (1998)
used the modified Picard and Newton iteration schemes to solve boththe h-based and mixed forms of the Richards equation governing infiltration in homogeneous and heterogeneous dry porous media. Theirresults showed that, for the same numerical accuracy, a combination of
the modified Picard and Newton scheme required fewer iterations thanthe modified Picard or Newton schemes as such. Kavetski et al. (2001)
introduced an adaptive time stepping with embedded error control tosolve the mixed form of the Richards equation. Their second-orderadaptive time stepping algorithm improved the performance of thenonlinear solver by providing accurate initial solution estimates for the

Corresponding author.
E-mail address: duc.ngo@usq.edu.au (D. Ngo-Cong).

https://doi.org/10.1016/j.jhydrol.2019.124240
Received 10 July 2019; Received in revised form 7 October 2019; Accepted 14 October 2019
Available online 17 October 20190022-1694/ © 2019 Elsevier B.V. All rights reserved.

Journal of Hydrology 580 (2020) 124240

D. Ngo-Cong, et al.

iterative process. Farthing et al. (2003) formulated a mixed hybrid finite element (MHFE) method and an enhanced cell-centred difference
scheme in combination with the method of lines for both the h-based
and mixed forms of the Richards equation. Fahs et al. (2009) combined
a lumped MHFE method with the method of lines for a more accurate
resolution of the Richards equation. They used the lumped MHFE approach (Younes et al., 2006) for approximating the flux to avoid unphysical oscillations for triangular meshes without additional numericalerrors. Other recent developments in numerical methods and software
for solving the Richards equation can be found in a review paper by Zha
et al. (2019).
Warrick et al. (1985) presented a semi-analytical solution to the
Richards equation subject to the van Genuchten (1980) soil hydraulicfunctions. Their solution has been used in many studies as a benchmark
for verifying numerical solutions (Caviedes-Voullième et al., 2013;
Phoon et al., 2007). Caviedes-Voullième et al. (2013) presented a finite
volume method in combination with four temporal integration schemes
to solve the 1D Richards equation. The four schemes included explicit
mixed (EM), implicit mixed (IM), explicit pressure head (EP), and implicit pressure head (IP) schemes. They found that the EM and IM
schemes produced accurate and conservative solutions, but not the EP
and IP schemes. The EM scheme, furthermore, was not able to solve theRichards equation under saturated conditions, while the IM scheme wasable to do so. Hayek (2016) further developed an analytical solution of
the 1D Richards equation for semi-infinite soils using the special case ofan exponential water retention curve and a power law for the hydraulicconductivity. Analytical and semi-analytical solutions of this type areimportant to test the accuracy of proposed numerical schemes.
In this study we focus on the use of schemes using radial basis
functions (RBFs). Such a scheme was first proposed by Kansa (1990) for
solving partial differential equations (PDEs) with a high level of accuracy for both gridded and scattered data. Unlike Kansa’s approach, MaiDuy and Tran-Cong (2001) suggested the use of integration to construct
the global RBF approximations, which improved the accuracy of numerical solutions. In related work, Wright and Fornberg (2006) presented a generalised finite difference formula based on compact radialbasis function approximations for solving partial differential equations
on irregular domains with scattered nodes. Šarler and Vertnik (2006)
developed a meshless local RBF collocation method for the solution ofdiffusion problems. The method overcame ill-conditioning problems ofthe global RBF collocation method (Kansa, 1990) and mitigated thesensitivity problem of choosing the RBF shape parameter. The scheme
has been used to solve several industrial problems such as thermodriven fluid flow problems with a free surface (Hon et al., 2015) andturbulent molten steel flow and solidification (Vertnik et al., 2019).
Divo and Kassab (2007) developed an efficient localised meshless collocation method using the Hardy Multiquadrics RBF for solution of
conjugate heat transfer problems.
Schemes using compact integrated radial basis function (CIRBF)
stencils were originally proposed by Mai-Duy and Tran-Cong (2011,2013) for solving partial differential equations. Their CIRBF schemes
incorporate not only nodal function values but also nodal values ofderivatives into the local integrated radial basis function (IRBF) approximations. Ngo-Cong et al. (2017) presented a generalised finite
difference scheme based on the CIRBF approximation to simulate waterflow in 1D and 2D heterogeneous soils. Hoang-Trieu et al. (2013)
subsequently incorporated the CIRBF approximation into the control
volume formulation (called CV-CIRBF) for solving second-order differential problems in one and two dimensions, with application to naturalconvection flow. Their numerical results show that the CV-CIRBF
method yields more accurate results than the standard finite volume
method (FVM) and a point collocation CIRBF-based method.
In this article, we further extend the CV-CIRBF method to solve thehighly nonlinear Richards equation governing water flow in unsaturated or variably saturated soils. Instead of using the nodal valuesof the governing equation as in Hoang-Trieu et al. (2013), the nodal

second derivatives of a field variable are added into the local IRBF
stencils when converting the radial basis function network-weight space
into physical space. In the CIRBF scheme of Ngo-Cong et al. (2017), theextra information of second derivatives at some nodal point i were
implicitly calculated from the function values at all nodal points over a
grid line passing through point i. As a result, this scheme (referred to asa global CIRBF scheme) yielded a fully populated system matrix associated with a grid line. In contrast, the present CV-CIRBF scheme produced a tridiagonal system matrix, which is computationally far more
efficient than the global CIRBF scheme. In Section 3.1, we will showthat calculations of internodal fluid fluxes using the global CIRBF
scheme (instead of the local CIRBF) can improve the solution accuracy
and efficiency in solving the Richards equation. We note that this internodal flux treatment does not affect the tridiagonal system matrix
produced by the present CV-CIRBF scheme. A rational function transformation technique (Pan and Wierenga, 1995) in combination with
Picard iterations is subsequently employed to treat the nonlinearity of
the Richards equation. Our proposed CV-CIRBF method is first illustrated and verified in Section 3.1 through comparisons with the semianalytical solution of Warrick et al. (1985). The numerical approach isthen used in Sections 3.2 and 3.3 to solve several unsaturated flow
problems to demonstrate its performance.
2. Control volume – Compact integrated radial basis function
discretisation of the Richards equation
In this section, a control volume approach based on compact integrated radial basis function (CV-CIRBF) stencils is used to discretise
the Richards equation governing moisture motion in soils. Consider Nznodal points on a z grid line. Each node z i is surrounded by a control
volume i . For interior nodes z i with 2 i Nz 1, i is defined as[z i 1/2, z i + 1/2]. For boundary nodes z1 andz Nz , the control volumes are[z1, z 3/2] and [z Nz 1/2, z Nz ], respectively. In the compact integrated radialbasis function (CIRBF) approximation as used by Mai-Duy and TranCong (2011, 2013), the first and second spatial derivatives of the field
variable u at a nodal point z(i) are computed as follows:
2u(i)

z2

=

1 ū +2u(i )= µ1 u¯ + µ 2zwhere u¯ = (u(i 1)((2u(i 1)2u(i + 1)z2

z2

2u(i 1)2u(i + 1)

z2

z2

)

(2)

)

(3)

T

T

and µ are known vectors of dimensionof 1 × 5; 1 = (1: 3); 2 = (4: 5); µ1 = µ (1: 3) and µ 2 = µ (4: 5) .
Readers are referred to the Appendix for further detail on the CIRBF
approximation. Application of Eqs. (2) and (3) at nodal pointsz = ( z (1) z (2) . . . z (Nz ) )T on the z-grid line yields
2u (m + 1)

z2

u(i)

T
u(i + 1) ) ;m)m)= D2z u (m + 1) + k 2(zk, k 2(zk= D2zku (m + 1)m)m)= D1z u (m + 1) + k1(zk, k1(zk= D1zkz

2u (m )z2
2u ( m )

z2

,

(4)

,

(5)

where u = (u(1) u(2) . . . u(Nz ) )T ; D1z , D1zk , D2z and D2zk are known
matrices of dimension of Nz × Nz ; and (m) and (m + 1) are the previous
and current time levels for the Picard iteration, respectively.
Similarly,
we
apply
Eqs.
(4)
and
(5)
at
T
T
z + = ( z (3/2) z (5/2) . . . z (Nz 1/2) z (Nz ) ) and z =( z (1) z (3/2) . . . z (Nz 3/2) z (Nz 1/2) )to calculate the second and first derivatives of the function u at thesurfaces of each control volume as follows:
2u (m + 1)

z2

+

m)
(m )
= D2z +u (m + 1) + k 2(zk+, k 2zk + = D 2zk +

2u (m)z2

2u (m )u (m + 1)m)(m )= D1z +u (m + 1) + k1(zk,+, k1zk + = D1zk +
z +
z2

2

,

(6)
(7)

Journal of Hydrology 580 (2020) 124240

D. Ngo-Cong, et al.

and
2u (m + 1)

2u (m)m)m)= D2z u (m + 1) + k2(zk, k 2(zk= D2zkz2u (m + 1)zm)= D1z u (m + 1) + k1(zk
,

m)
k1(zk

z2
2 u (m )= D1zkz2,

Table 1
Van Genuchten-Mualem (VG) soil hydraulic parameters of the three soils considered in this study.

(8).

Soil texture(9)
CoarseMediumFine

Here, we adopt the rational function transformation method proposed by Pan and Wierenga (1995) to transfer the pressure head h in
Eq. (1) into a transformed function f as follows:
h
1+ h

f=
where

if h

1

h=

f

if h

(11)

0.

(1 +
h
=
1
f

if h < 0
if h 0.

·(K(13)·(K (h) f )

K
= 0.
z

f)

C

(14)

(n )

t

z

K (n+ 1, m)z

K (n + 1, m)=0zf (n+ 1, m + 1)

= K+

where

(n + 1, m) + C (n + 1, m )C (n + 1, m)
t
=

z
z

f (n + 1, m)

K (n + 1, m)
f (n+ 1, m)z

K (n+ 1, m)

z

K

z

C (n + 1, m)
t

K
+
z

=

z

+

(

+ K (n + 1, m)
=

z

(K

) (z

(n + 1, m)(n + 1, m + 1)

z

(n )

f

( n)
i

dz

+

+

2f (m + 1)

z2

fi + 1/2)(n + 1, m )t(n + 1, m)

D1z + + K

(n + 1, m )K(m + 1)(m + 1)f (n + 1, m) (zi 1/2 )z(n + 1, m)

D1z.

K (z i + 1/2)

(21)

K (z i 1/2)

(22)

(n + 1, m + 1)

k1z
(n)

(n + 1, m)zt(m + 1)m)(m )+ k 2(zk+, k 2zk + = D2zk +m)(m )+ k1(zk+, k1zk + = D1zk +(m + 1)= D2z +f

(m + 1)z(19)= D1z +f

= D1z f(m + 1)(m + 1)= D2z f(m + 1)zi 1/2)qz(n + 1, m) (z i + 1/2) + qz(n+ 1, m) (z i 1/2 ) +

K+

= D2z

= D1z

z2

(n + 1, m)

t

) (z

= D1z +

2f (m + 1)

(18)

t

qz(n + 1, m)

f (n + 1, m) (z i + 1/2 )
z

(23)

m)
m)
+ k 2(zk, k 2(zk= D2zkm)m)+ k1(zk, k1(zk= D1zk2 (m )(24)z2
2 (m )(25)z2
2 (m )(26)z2
2 (m )

(27)

z2

and

(n + 1, m)(n )K (n+ 1, m) (z i 1/2 )k1z +(m + 1)zz(n + 1, m + 1)+z(n + 1, m)(n + 1, m + 1) dz+

z2. Substitution of Eq. (17) into (15) yields+

K (n+ 1, m) (z i + 1/2 )= D2z +

2 (m + 1)

The Richards equation (19) is now integrated over each control
volume i as follows:i0.24310.02700.0043(m + 1)(17)+ O ( 2)(n + 1, m)z2which can be rewritten in the form:

C (n+ 1, m)
t

2.28
2.00
1.31

z. I

2 (m + 1)(16)= f (n+ 1, m + 1)

0.12400.01890.0190where I is the identity matrix. The first and second spatial derivatives ofand f at the control surfaces z i + 1/2 and z i 1/2 are discretised by using
the CIRBF scheme (Eqs. (6)-(9)) as follows:(15)

Application of a Taylor series expansion about the point f (n + 1, m)
results in(n + 1, m + 1) =

0.410
0.420
0.410

qz(n + 1, m) (z i + 1/2) + qz(n + 1, m) (z i 1/2 ) +where the superscript n refers to a physical time level, and m to thePicard iteration. The specific moisture capacity function C (h) = d /dh
is transformed as follows:

h
d
C =C
=
.
f
df

(n + 1, m)t

The vertical flux is defined as qz = K f / z K . The nonlinear
term in Eq. (14) is now linearised using the Picard iteration method,which for a one-dimensional system leads to(n + 1, m + 1)

0.057
0.061
0.095

respectively, in which the internodal conductivity given by the arithmetic
mean
of
the
nodal
values,i.e.,
K*(zi+1/
2) = [K*(zi) + K*(zi+1)]/2, is assumed constant over each control volume. The integrals in Eq. (20) are approximated using the midpoint
rule. Application of Eq. (20) to all control volumes i on the z-axis
results in(12)

where K = K h/ f . Substituting Eq. (13) into Eq. (1) givest

Ks (cm/min)qz(n + 1, m) (z i 1/2) =Making use of the chain rule of partial derivatives, we obtain·[K (h) h] =η (−)

qz(n + 1, m) (z i + 1/2) =Taking the derivative of both sides of Eq. (11) and using Eq. (10) toeliminate f results inh) 2

α (cm−1)

where the vertical fluxes at the control surfaces z i + 1/2 and z i 1/2 aregiven byif h < 0fθs (cm3/cm3)

(10)

0

is a transformation constant. Rearranging Eq. (10) givesfθr (cm3/cm3)

θr: residual water content; θs: saturated water content; α: VG multiplier of thepressure head; η: VG exponent, Ks saturated hydraulic conductivity.

if h < 0hvan Genuchten hydraulic parameters(m + 1)(m + 1)m)(m )+ k 2(zkf+, k 2zkf + = D2zk +m)(m )+ k1(zkf+, k1zkf + = D1zk +m)m)+ k2(zkf, k 2(zkf= D2zkm)m)+ k1(zkf, k1(zkf= D1zk

2f (m )

(28)z2
2 f (m )(29)

z2
2f (m)z2
2f ( m )

z2

.

(30)
(31)

Substitution of these equations into Eq. (23) leads to a system of
linear equations, which is solved to obtain at each physical time step.
Note that the resulting system matrix (for the 1D case) has only 3 non-(20)
3

Journal of Hydrology 580 (2020) 124240

D. Ngo-Cong, et al.

(a)

(b)

Fig. 1. Calculated pressure head (h) distributions versus depth (z) obtained with the CV-CIRBF method using different time steps t (a) and nodal spacings (b),compared with the semi-analytical solution of Warrick et al. (1985) (a, b) and the finite difference solution (FDM) of Celia et al. (1990) (b). Calculations are for the
medium-textured soil subject to ponded infiltration.
Table 2
Example 1: Comparisons of the relative error norm (Ne), total number of iterations, and total CPU time required to obtain converged solutions with TOL = 10−5
between the FDM and CV-CIRBF using Approaches 1 and 2. The time step was set at 0.5 min for all cases. FDM, App.1 and App. 2 stand for the finite difference
method of Celia et al. (1990), Approach 1 and Approach 2, respectively.
Grid

21
31
41
51
61

Ne(z)

Number of iterations

CPU time (s)FDMApp. 1

App. 2

FDM

App. 1

App. 2

FDM

App. 1

App. 2

2.25E−01
1.91E−01
1.80E−01
1.58E−01
1.44E−−01

4.63E−02
5.91E−02
7.77E−02
9.53E−02
9.82E−02

4.20E−02
5.91E−02
7.77E−02
9.53E−02
9.81E−02

4735
4788
4847
4937
5104

7947
7451
7537
7591
7638

4737
4789
4846
4898
4929

3.86
4.03
4.28
4.73
5.10

6.27
6.66
7.21
7.95
8.34

4.16
4.37
4.69
5.22
5.42

3. Numerical results and discussion
In this section, the CV-CIRBF method in combination with Picard
iteration and the rational function transformation is used to solve theRichards equation in one dimension (1D) for different boundary conditions (constant pressure head and constant flux), initial pressure
heads (i.e., h 0 = 800, 10000 and 50000 cm ) and soil types (relatively coarse-, medium- and fine-textured media). We used the standard
van Genuchten-Mualem (VG) model to describe the soil water retention
curve, θ(h), and the relative hydraulic conductivity function, K(h), as
follows (van Genuchten, 1980):s(h ) =r[1 + ( | h |) ]µ

+

< TOL

0
(33)

1

1

Seµ

µ 2

if h

if h > 00(34)where θs and θr are, respectively, the saturated and residual watercontents, α and η are fitting parameters, Ks is the saturated hydraulic1/ ,conductivity, Se = (r )/( sr ) is effective saturation, µ = 1
and L = 0.5. Table 1 presents the van Genuchten-Mualem (VG) hydraulic parameters of the three soils we considered. For the relatively
coarse-textured soil we used the hydraulic parameters of a loamy sandas provided by Carsel and Parrish (1988). For the medium soil we took
the hydraulic parameters from a study by Vogel and Hopmans (1992)
on furrow irrigation. The soil in their study was categorized as a sandbut in reality had hydraulic parameters more typical of a mediumtextured soil as indicated by the class-pedotransfer functions of Carsel
and Parrish (1988) and Schaap et al. (2001). And for the fine-textured

zero entries on each row. At each physical time step, the Picard iteration is terminated when the following criterion is fulfilledNz(f (m + 1)
f i(m) )2i=1 i
Nz(m + 1) 2
(f
)
i=1 i

Ks SeL 1
Ks

Fig. 2. Calculated pressure head (h) distributions versus depth at selected times
as obtained with the CV-CIRBF scheme, the HYDRUS-1D numerical solution(Šimůnek et al., 2008) and the semi-analytical solution of Warrick et al. (1985).
Results are for the medium-textured soil subject to ponded infiltration.

h

if h > 0s

K (h) =

r if

(32)

where TOL is a prescribed tolerance, presently set at 10 5 .
4Journal of Hydrology 580 (2020) 124240

D. Ngo-Cong, et al.

(a)

(b)

Fig. 3. Calculated pressure head (h) distributions versus depth at selected times for initial conditions of −10,000 cm (a) and −50,000 cm (b) as obtained with the
CV-CIRBF and HYDRUS-1D numerical solutions assuming 21 and 51 nodes (grid spacings of 5 and 2 cm, respectively). Calculations are for the medium-textured soil
subject to ponded infiltration.

• Approach 1: Internodal fluxes are calculated using the present local
CIRBF scheme, Eq. (5).
• Approach 2: Internodal fluxes are calculated using the global CIRBF

soil we used the hydraulic parameters of a clay loam soil as categorized
by Carsel and Parrish (1988).
For all simulations we considered vertical downward infiltration
into a 100-cm deep soil profile with the lower boundary at z = 0 cm andthe upper boundary at z = 100 cm. In this study we limited ourselves to
fixed nodal spacings (Δz).

scheme (Ngo-Cong et al., 2017), Eq. (A.11).

The CV-CIRBF formulation was verified by comparing results withthe semi-analytical solution (Celia et al., 1990). Performance was assessed based on relative error norms between the numerical solutions
and the semi-analytical solution of Warrick et al. (1985) at three marks(Fig. 1) corresponding to Se = 0.25, 0.5 and 0.75 where
Se = (Se Se0)/(1 Se0) and Se0 is the initial effective saturation.
Table 2 shows comparisons of the relative error norm, Ne(z), the total
number of iterations, and the total CPU time required to obtain theconverged solution with TOL = 10−5. The relative error norm, Ne(z), is
calculated as

3.1. Example 1: Constant pressure head at upper and lower boundaries
Initial and boundary conditions for this first example areh (0, t ) = hbottom = 800 cm ,h (z, t = 0) = h 0 = 800 cm ,andh (100, t ) = htop = 0 cm . In this example, we used fixed time steps (Δt).
Internodal fluxes are known to have errors when changes in the watercontent are significant, especially for heterogeneous soils due to discontinuities in the water content at the interface of different soil layers.
In this example, we used two different approaches to calculate the internodal fluxes as given by Eqs. (21) and (22):

(a)

(b)

Fig. 4. Calculated pressure head (h) distributions versus depth at selected times for initial conditions of −800 cm (a) and −50,000 cm (b) as obtained with the CVCIRBF and HYDRUS-1D numerical solutions assuming 21 and 51 nodes (grid spacings of 5 and 2 cm, respectively). Calculations are for the coarse-textured soil subject
to ponded infiltration.

5

Journal of Hydrology 580 (2020) 124240

D. Ngo-Cong, et al.

(a)

(b)

Fig. 5. Calculated pressure head (h) distributions versus depth at selected times for initial conditions of −800 cm (a) and −50,000 cm (d) as obtained with the CVCIRBF and HYDRUS-1D numerical solutions assuming 21 and 51 nodes (grid spacings of 5 and 2 cm, respectively). Calculations are for the fine-textured soil subject toponded infiltration.

influence of the time step on the CV-CIRBF results for the mediumtextured soil using t = {0.5, 1, 10, 30} min and a grid of 51 nodes.
The Warrick semi-analytical solution (Phoon et al., 2007; Warrick et al.,1985) and the finite difference method (FDM) results of Celia et al.
(1990) are also depicted for comparison purposes. Distributions of thepressure head for different grid sizes Nz = {21, 51} and t = 0.5 min are
described in Fig. 1b. The results indicate that the present method yields
more accurate solutions when refining the grid and time step sizes andthat the CV-CIRBF solutions are more accurate than those based on the
FDM scheme.
Fig. 2 compares similar pressure head distributions versus depth,
now at four different times, with results obtained using both the semianalytical solution of Warrick et al. (1985) and the HYDRUS-1D software of Šimůnek et al. (2008). Results show a slightly better performance of the CV-CIRBF scheme as compared to HYDRUS-1D. Calculations in Fig. 2 were obtained for an initial condition of −800 cm.
Results for more extreme initial conditions of −10,000 cm and−50,000 cm are shown in Fig. 3 assuming nodal spacings of 5 and 2 cm
(i.e., Nz values of 21 and 51, respectively). The distributions indicate
only a marginal effect of nodal spacing on the CV-CIRBF results, butwith the HYDRUS-1D distributions depending far more on the numberof finite element nodes.
The CV-CIRBF numerical approach was also used to study infiltration in more coarse- and fine-textured soils whose hydraulic parameters
are presented in Table 1. Figs. 4 and 5 show the CV-CIRBF and HYDRUS
results at different times for the coarser and finer soils, respectively,
using initial pressure heads (h0) of −800 and −50,000 cm. Similarly as
the results shown in Figs. 2 and 3, the CV-CIRBF calculations agreed
well with the HYDRUS-1D results, with HYDRUS-1D again showing
more of an effect of nodal spacing on the computed pressure head
distributions.

Fig. 6. Calculated pressure head (h) distributions versus depth at timet = 780 min for initial conditions of −800 cm as obtained with the CV-CIRBF
and HYDRUS-1D numerical solutions assuming 51 nodes (grid spacings of
2 cm). Calculations are for the medium-textured soil subject to constant infiltration.
3

Ne (z ) =

3

(z k
k=1

z¯k )2 /
k=1

z¯k2

(35)

where z k and z̄ k are the numerical and semi-analytical values at the k-th
mark, respectively. The time step was set at 0.5 min for all cases. For the
purpose of CPU time comparisons, all related computations were carried out on a single 3.10 GHz processor machine with 8 GB RAM.
The CV-CIRBF methods yielded more accurate solutions than the
FDM of Celia et al. (1990), with Approach 2 solutions being slightly
more accurate than those of Approach 1. Approach 2 was also more
efficient than Approach 1 in terms of the total number of iterations andthe total CPU time. Since Approach 2 outperformed Approach 1 interms of accuracy and efficiency, we will use here further only the CVCIRBF method with Approach 2. For a given grid size, Approach 2 wasslower than the FDM. However, Approach 2 achieved a given level ofaccuracy with a much coarser grid and hence was more efficient. For
instance, Approach 2 with a grid of 21 nodes produced better results in
4.16 s than the FDM with a grid of 61 nodes in 5.1 s. Fig. 1a shows the

3.2. Example 2: Constant flux at upper boundary and constant pressure
head at lower boundary
In this example, we apply a constant flux at the upper boundary (qA)
and again a constant pressure head at the lower boundary (hC). The
constant flux at the upper boundary was set at half of the saturated
hydraulic conductivity, i.e., qA = 0.5Ks. Fig. 6 shows the influence ofthe time step (Δt) on the CV-CIRBF results for infiltration in themedium-textured soil, in comparison also with the HYDRUS-1D solutions, assuming an initial pressure heads of −800 cm. The CV-CIRBF
result became closer to the HYDRUS-1D result with reducing the time6

Journal of Hydrology 580 (2020) 124240

D. Ngo-Cong, et al.

(a)

(b)

Fig. 7. Calculated pressure head (h) distributions versus depth at time t = 780 min for initial conditions of −800 cm (a) and −50,000 cm (d) as obtained with the
CV-CIRBF and HYDRUS-1D numerical solutions assuming 21 and 51 nodes (grid spacings of 5 and 2 cm, respectively). Calculations are for the medium-textured soil
subject to constant infiltration.

step size. We note that an adaptive time-stepping scheme was used toobtain the HYDRUS-1D result.
Fig. 7 presents the CV-CIRBF and HYDRUS-1D results for infiltration
into the medium-textured soil assuming initial pressure heads (h0) of−800 and −50,000 cm, again for 21 and 51 nodes (grid spacings of 5
and 2 cm, respectively). For the calculations here we now used anadaptive time-stepping scheme to improve the computational efficiency
and robustness. Following Kirkland et al. (1992), the time step (Δt) was
increased by 10% when the number of Picard iterations at the previous
time step was less than 4, and decreased by 10% when the number of
iterations was greater than 8. The plots in Fig. 7 again show excellent
agreement between the CV-CIRBF and HYDRUS-1D results, especiallyfor the finer nodal spacing, with the CV-CIRBF distributions being lessaffected by the number of nodes as compared to HYDRUS-1D (especially for infiltration in the drier soil). Our CV-CIRBF solver produced amass-conservative solution as demonstrated in Fig. 8 showing the massbalance measure (MB) versus time t. The mass balance measure is defined as MB = m1/m2 with m1 and m2 being the total additional mass inthe domain and the total net flux into the domain, respectively. Except
for some minor deviations at early times, MB values were always at orvery close to 1.00. Results for the coarse- and fine-textured soils were
found to be very similar and are not further shown here.

Fig. 8. Calculated mass balance measure (MB) versus time t obtained with the
CV-CIRBF scheme.
Table 3
Flow in a layered soil: initial and boundary conditions, and elapsed times. Notethat h0 is the initial pressure head and qA is the vertical flux at the top boundary.
Caseh0 (cm)qA (cm/h)

Elapsed Time (h)

1
2
3
4

−200
−50000−200−50000

0.3
0.3
1.25
1.25

4
12
3.8
6

3.3. Example 3: Flow in a layered soil with a constant flux at upperboundary and zero-flux at lower boundary
This example considers infiltration into a layered soil assuming a
constant flux at the top boundary, a zero-flux at the bottom boundary,and different initial pressure heads as shown in Table 3. For the threelayered soil profile we used Berino loamy fine sand from 0 to 50 cm and
90 to 100 cm, and Glendale clay loam from 50 to 90 cm. The VG hydraulic parameters of the two soils are given in Table 4. Fig. 9 presents
the FDM and CV-CIRBF results for infiltration assuming an initialpressure head of −200 cm and a flux of 0.3 cm/h at the top boundary.
The CV-CIRBF yielded more accurate solution than the FDM using agrid of 51 nodes. The CV-CIRBF solution using 101 nodes was almost
the same as the FVM solution of McBride et al. (2006) who used aconventional finite-volume-based CFD code (Fig. 9b). Fig. 10 similarlyshows excellent agreement between the CV-CIRBF results and those of
McBride et al. (2006) for qA = 0.3 and 1.25 cm/h, assuming again an
initial condition (h0) of −200 cm, but also for the far more extreme casewhere h0 = −50,000 cm.

Table 4
Flow in a layered soil: van Genuchten soil hydraulic parameters (McBride et al.,2006).
Soil

Berino loamy finesand
Glendale clay loam

van Genuchten hydraulic parametersθr (cm3/
cm3)

θs (cm3/cm3)

α (cm−1)η (–)

Ks (cm/min)

0.02860.36580.02802.23900.37570.10600.46860.0104

1.39540.0091

θr: residual water content; θs: saturated water content; α: VG multiplier of thepressure head; η: VG exponent, Ks saturated hydraulic conductivity.
7

Journal of Hydrology 580 (2020) 124240

D. Ngo-Cong, et al.

(a)

(b)

Fig. 9. Calculated pressure head (a) and volumetric water content (b) distributions versus depth for Case 1 as obtained with the CV-CIRBF and FDM schemes, and theFVM solution of McBride et al. (2006) assuming 51 and 101 nodes (grid spacings of 2 and 1 cm, respectively).

(a)

(b)

Fig. 10. Calculated volumetric water content distributions versus depth for all cases as obtained with the CV-CIRBF scheme, and FVM solution of McBride et al.
(2006) using 101 nodes (grid spacing of 1 cm) for the flux of 0.3 cm/h (a) and 1.25 cm/h (b) at the top boundary.

4. Conclusions

upon the nodal spacing. The CV-CIRBF scheme was also used successfully to simulate flow in a layered soil profile, with its solutions being in
excellent agreement with the standard FVM solutions of McBride et al.
(2006). Finally, we acknowledge that in this study we limited ourselvesto constant nodal spacing (Δz). Extending the solutions to irregular grid
systems is relatively straightforward and has not been further addressedin this paper.

We successfully developed a control volume technique based on acompact integrated radial basis function approximation (CV-CIRBF) for
solving the nonlinear Richards equation governing water flow in unsaturated and variably saturated soils. Performance of the proposed
method appears satisfactory as demonstrated with ponded and constantflux infiltration in a one-dimensional vertical soil profile. The local
CIRBF stencil was constructed based on information associated withthree nodal points, which produced a tridiagonal system matrix similarly as with the conventional control volume method. Results showedthat our calculation of internodal fluid fluxes using the global CIRBF
scheme was more accurate than the approach used earlier for the local
CIRBF scheme. The CV-CIRBF outperformed the FDM scheme used by
Celia et al. (1990), especially when a coarser grid was used (Example1). The proposed numerical solutions agreed well with the semi-analytical solution of Warrick et al. (1985), as well as with HYDRUS-1D
results, for different boundary conditions (constant pressure heads orfluxes), initial pressure heads (h 0 = 800, 10000 and 50000 cm ) anda range of soil textural classes exhibiting contrasting hydraulic behaviour. The HYDRUS-1D results closely duplicated the CV-CIRBF results
only when relatively small nodal spacing were used in the HYDRUS-1D
finite element solutions, with the CV-CIRBF being far less dependent

Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
This research was funded in part by the University of Southern
Queensland, Australia, through the University's Strategic Research
Fund initiative, and in part by the Basin Studies Laboratory (LEBAC) of
the Department of Geology, associated with the Center for
Environmental Studies (CEA) of UNESP in Brazil. The authors would
like to thank the editor and reviewers for their helpful comments.

8

Journal of Hydrology 580 (2020) 124240

D. Ngo-Cong, et al.

AppendixCompact integrated radial basis function approximation{i

The field variable u can be approximated over a compact integrated radial basis function (CIRBF) stencil corresponding to three nodal points
1, i , i + 1} as follows (Mai-Duy and Tran-Cong, 2013):
2u (z )z2

i+1

w (k ) G (k) (z ),

=

(A1)

k=i 1
i+1

u (z )(k )=w (k) H[1](z ) + c1,
z
k=i 1
i+1

u (z ) =

(A2)

(k )
w (k) H[0](z ) + c1 z + c2,

(A3)

k=i 1

where
are radial basis function (RBF) weights that need to be determined;
are known RBFs; H[1] (z ) =G (z ) dz; H[0] (z ) = H[1] (z ) dz ; and c1 and c2 are integration constants, which are also unknown. The RBFs G (k) (z ) are defined as{w (k ) }ik+=1i 1

G (k ) ( z )G (k) (z ) = (z z (k) ) 2 + a(k )2 where a(k ) is the RBF width determined presently as a(k ) = g , in which is a positive factor, and g the grid size. The
influence of on the CIRBF solution accuracy was investigated in previously published papers (Hoang-Trieu et al., 2013; Mai-Duy and Tran-Cong,
2011; Mai-Duy and Tran-Cong, 2013; Ngo-Cong et al., 2017). Based on our numerical experiments, the value of
was chosen to be 20 for the
calculations in this paper. We used the nodal function values and the nodal second derivatives as extra information to construct the stencil approximations for the first and second derivatives. Application of Eq. (A.3) at three nodal points {i 1, i , i + 1} and Eq. (A.1) at two nodal points{i 1, i + 1} leads to the following system:
(i 1)
(i )
(i + 1)H[0](z i 1 ) H[0](z i 1 ) H[0](z i 1 ) z i 1 1u¯2u(i 1)

z2

(i 1)
H[0]
(z i )=

2u(i + 1)

z2

(i )
H[0]
(z i )(i + 1)H[0](z i )

zi

1

(i 1)
(i )
(i + 1)H[0](z i + 1 ) H[0](z i + 1 ) H[0](z i + 1 ) z i + 1 1
(i 1)
G[0]
(z i 1 )(i )G[0](z i 1 )(i + 1)G[0](z i 1 )

0

0

(i 1)
G[0]
(z i + 1 )

(i )
G[0]
(z i + 1 )(i + 1)G[0](z i + 1 )00Tu(i + 1) ) ; w¯ = ( w (i 1)w (i )( wc¯¯ )
(A4)

C

or

( wc¯¯ ) = C (u¯12u(i 1)2u(i + 1)z2

z2

in which u¯ = (u(i 1)

u(i)

)

T

(A5)
T
w (i + 1) ) and c = ( c1c2 )T . Inserting Eq. (A.5) into Eqs. (A.1) and (A.2), and applying them at

nodal point i , expressions for first and second derivatives of the field variable u in physical space can be obtained:
2u(i)

z2

(

= G (z (i) ) C 1 u¯

2u(i + 1)z2

z2

2u(i 1)2u(i + 1)

z2

z2

(

u(i )
= H[1] (z (i) ) C 1 u¯zµ)2u(i 1)

T

(A6)

)

T

(A7)

or

(
u
= µ u¯ + µ (
z
2u(i)

z2

=

1 ū +

2

(i )

1

2

2u(i 1)2u(i + 1)z2

z2

2u(i 1)2u(i + 1)

z2

z2

)
)

T

(A8)

T

(A9)

where and µ are known vectors of dimension of 1 × 5; 1 = (1: 3); 2 = (4: 5); µ1 = µ (1: 3) and µ 2 = µ (4: 5) . Eqs. (A.8) and (A.9) were solved
with the Picard iteration method, with the extra information of the second derivatives obtained from the previous iteration.
Instead of using Picard iteration, the extra information of the second derivatives in (A.8) and (A.9) can be implicitly calculated using the function
values at all nodes on a grid line as shown in Ngo-Cong et al. (2017). The values of the second and first derivatives of u against z at the nodal points
on the z-grid line are then given by
2u (m + 1)

z2
u (m + 1)z= D2z u (m + 1) ,(A10)

= D1z u (m + 1) ,(A11)

where D1z and D2z are fully populated Nz × Nz matrices.

9

Journal of Hydrology 580 (2020) 124240

D. Ngo-Cong, et al.

References

second-order elliptic differential problems. J. Comput. Phys. 230 (12), 4772–4794.
https://doi.org/10.1016/j.jcp.2011.03.002.
Mai-Duy, N., Tran-Cong, T., 2013. A compact five-point stencil based on integrated RBFs
for 2D second-order differential problems. J. Comput. Phys. 235, 302–321. https://doi.org/10.1016/j.jcp.2012.10.048.
McBride, D., Cross, M., Croft, N., Bennett, C., Gebhardt, J., 2006. Computational modelling of variably saturated flow in porous media with complex three-dimensional
geometries. Int. J. Numer. Methods Fluids 50 (9), 1085–1117. https://doi.org/10.
1002/fld.1087.
Ngo-Cong, D., et al., 2017. A generalised finite difference scheme based on compact integrated radial basis function for flow in heterogeneous soils. Int. J. Numer. Methods
Fluids 85 (7), 404–429. https://doi.org/10.1002/fld.4386.
Pan, L., Wierenga, P.J., 1995. A transformed pressure head-based approach to solve
Richards' equation for variably saturated soils. Water Resour. Res. 31 (4), 925–931.
https://doi.org/10.1029/94WR03291.
Phoon, K.-K., Tan, T.-S., Chong, P.-C., 2007. Numerical simulation of Richards equation inpartially saturated porous media: under-relaxation and mass balance. Geotech. Geol.
Eng. 25 (5), 525–541. https://doi.org/10.1007/s10706-007-9126-7.
Šarler, B., Vertnik, R., 2006. Meshfree explicit local radial basis function collocation
method for diffusion problems. Comput. Math. Appl. 51 (8), 1269–1282. https://doi.
org/10.1016/j.camwa.2006.04.013.
Schaap, M.G., Leij, F.J., Van Genuchten, M.T., 2001. Rosetta: A computer program forestimating soil hydraulic parameters with hierarchical pedotransfer functions. J.
Hydrol. 251 (3–4), 163–176. https://doi.org/10.1016/S0022-1694(01)00466-8.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., Van Genuchten, M.T., 2008. The HYDRUS-1D
software package for simulating the movement of water, heat, and multiple solutes invariably saturated media, version 4.0: HYDRUS Software Series 3. Department ofEnvironmental Sciences, University of California Riverside, Riverside, California,
USA, pp. 315.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. https://doi.org/
10.2136/sssaj1980.03615995004400050002x.
van Genuchten, M.T., 1982. A comparison of numerical solutions of the one-dimensional
unsaturated—saturated flow and mass transport equations. Adv. Water Resour. 5 (1),47–55. https://doi.org/10.1016/0309-1708(82)90028-8.
Vertnik, R., Mramor, K., Šarler, B., 2019. Solution of three-dimensional temperature andturbulent velocity field in continuously cast steel billets with electromagnetic stirring
by a meshless method. Eng. Anal. Bound. Elem. 104, 347–363. https://doi.org/10.
1016/j.enganabound.2019.03.026.
Vogel, T., Hopmans, J.W., 1992. Two-dimensional analysis of furrow infiltration. J. Irrig.
Drain Eng. 118 (5), 791–806. https://doi.org/10.1061/(ASCE)0733-9437(1992)
118:5(791).
Warrick, A., Lomen, D., Yates, S., 1985. A generalized solution to infiltration. Soil Sci.
Soc. Am. J. 49 (1), 34–38. https://doi.org/10.2136/sssaj1985.
03615995004900010006x.
Wright, G.B., Fornberg, B., 2006. Scattered node compact finite difference-type formulas
generated from radial basis functions. J. Comput. Phys. 212 (1), 99–123. https://doi.
org/10.1016/j.jcp.2005.05.030.
Younes, A., Ackerer, P., Lehmann, F., 2006. A new mass lumping scheme for the mixed
hybrid finite element method. Int. J. Numer. Methods Eng. 67 (1), 89–107. https://doi.org/10.1002/nme.1628.
Zha, Y., et al., 2019. Review of numerical solution of Richardson-Richards equation forvariably saturated flow in soils. Wiley Interdiscip. Rev. Water e1364. https://doi.org/
10.1002/wat2.1364.

Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24 (5), 755–769. https://doi.org/10.
1029/WR024i005p00755.
Caviedes-Voullième, D., García-Navarro, P., Murillo, J., 2013. Verification, conservation,
stability and efficiency of a finite volume method for the 1D Richards equation. J.
Hydrol. 480, 69–84. https://doi.org/10.1016/j.jhydrol.2012.12.008.
Celia, M.A., Bouloutas, E.T., Zarba, R.L., 1990. A general mass-conservative numericalsolution for the unsaturated flow equation. Water Resour. Res. 26 (7), 1483–1496.
https://doi.org/10.1029/WR026i007p01483.
Divo, E., Kassab, A.J., 2007. An efficient localized radial basis function meshless method
for fluid flow and conjugate heat transfer. J. Heat Transf. 129 (2), 124–136. https://doi.org/10.1115/1.2402181.
Fahs, M., Younes, A., Lehmann, F., 2009. An easy and efficient combination of the mixed
finite element method and the method of lines for the resolution of Richards’ equation. Environ. Model. Softw. 24 (9), 1122–1126. https://doi.org/10.1016/j.envsoft.
2009.02.010.
Farthing, M.W., Kees, C.E., Miller, C.T., 2003. Mixed finite element methods and higher
order temporal approximations for variably saturated groundwater flow. Adv. Water
Resour. 26 (4), 373–394. https://doi.org/10.1016/S0309-1708(02)00187-2.
Farthing, M.W., Ogden, F.L., 2017. Numerical solution of Richards’ equation: a review ofadvances and challenges. Soil Sci. Soc. Am. J. 81 (6), 1257–1269. https://doi.org/10.
2136/sssaj2017.02.0058.
Hayek, M., 2016. An exact explicit solution for one-dimensional, transient, nonlinear
Richards’ equation for modeling infiltration with special hydraulic functions. J.
Hydrol. 535, 662–670. https://doi.org/10.1016/j.jhydrol.2016.02.021.
Hoang-Trieu, T.-T., Mai-Duy, N., Tran, C.-D., Tran-Cong, T., 2013. A finite-volume
method based on compact local integrated radial basis function approximations forsecond-order differential problems. CMES: Comput. Model. Eng. Sci. 91 (6), 485–516.
https://doi.org/10.3970/cmes.2013.091.485.
Hon, Y.-C., Šarler, B., Yun, D.-F., 2015. Local radial basis function collocation method forsolving thermo-driven fluid-flow problems with free surface. Eng. Anal. Bound. Elem.
57, 2–8. https://doi.org/10.1016/j.enganabound.2014.11.006.
Huang, K., Mohanty, B., Van Genuchten, M.T., 1996. A new convergence criterion for themodified Picard iteration method to solve the variably saturated flow equation. J.
Hydrol. 178 (1–4), 69–91. https://doi.org/10.1016/0022-1694(95)02799-8.
Kansa, E.J., 1990. Multiquadrics - A scattered data approximation scheme with applications to computational fulid-dynamics - II: Solutions to parabolic, hyperbolic andelliptic partial differential equations. Comput. Math. Appl. 19 (8–9), 147–161.
https://doi.org/10.1016/0898-1221(90)90271-K.
Kavetski, D., Binning, P., Sloan, S., 2001. Adaptive time stepping and error control in a
mass conservative numerical solution of the mixed form of Richards equation. Adv.
Water Resour. 24 (6), 595–605. https://doi.org/10.1016/S0309-1708(00)00076-2.
Kirkland, M.R., Hills, R., Wierenga, P., 1992. Algorithms for solving Richards' equationfor variably saturated soils. Water Resour. Res. 28 (8), 2049–2058. https://doi.org/
10.1029/92WR00802.
Lehmann, F., Ackerer, P., 1998. Comparison of iterative methods for improved solutionsof the fluid flow equation in partially saturated porous media. Transp. Porous Media31 (3), 275–292. https://doi.org/10.1023/A:1006555107450.
Mai-Duy, N., Tran-Cong, T., 2001. Numerical solution of Navier-Stokes equations using
multiquadric radial basis function networks. Int. J. Numer. Methods Fluids 37 (1),65–86. https://doi.org/10.1002/fld.165.
Mai-Duy, N., Tran-Cong, T., 2011. Compact local integrated-RBF approximations for10