# Hermite Radial Basis Collocation Method for Unsaturated Soil Water Movement Equation

Authors: Jiao Wang; Lijun Su; Xinqiang Qin
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Due to the nonlinear diffusion term, it is hard to use the collocation method to solve the unsaturated soil water movement equation directly. In this paper, a nonmesh Hermite collocation method with radial basis functions was proposed to solve the nonlinear unsaturated soil water movement equation with the Neumann boundary condition. By preprocessing the nonlinear diffusion term and using the Hermite radial basis function to deal with the Neumann boundary, the phenomenon that the collocation method cannot be used directly is avoided. The numerical results of unsaturated soil moisture movement with Neumann boundary conditions on the regular and nonregular regions show that the new method improved the accuracy significantly, which can be used to solve the low precision problem for the traditional collocation method when simulating the Neumann boundary condition problem. Moreover, the effectiveness and reliability of the algorithm are proved by the one-dimensional and two-dimensional engineering problem of soil water infiltration in arid area. It can be applied to engineering problems.


---

## Body

## 1. Introduction

Unsaturated soil moisture movement [1] referring to the water in the soil is not filled with all the pores of the water movement, and it is an important fluid movement form of porous media. It has important practical significance in extensive engineering sciences such as runoff calculation, groundwater resource evaluation, soil science, hydrology, farmland irrigation and drainage, and soil salinization prevention and management [2]. At present, the soil water movement equation mostly refers to the Richards’ equation, which is established based on the Buckingham-Darcy law and continuous equation. The specific expression with water content as dependent variable is shown as follows:(1)∂θ∂t=∂∂xDθ∂θ∂x+∂∂zDθ∂θ∂z-∂Kθ∂z.In recent years, differential method and finite element method have been used to solve this problem; for example, Lei and Yang [3] used the Ritz finite element method to solve the one-dimensional unsaturated soil water flow problem, used finite difference method to establish a numerical model for the one-dimensional flow of unsaturated soil water, and applied it to practical problems [4] successfully; Ma et al. [5] used ADI alternating implicit difference dispersion and Gauss-Seidel iteration to solve the soil water movement equation; Li et al. [6] proposed a radial basis collocation method, combined with the radial basis function and collocation method, for the one-dimensional soil water movement equation; Su et al. [7] solved the two-dimensional unsaturated soil water movement equation by the radial basis collocation with difference method; Hu et al. [8] solved the boundary value problem by weighted radial basis point method. Many scholars [9–20] have solved a series of the numerical solution problems based on the radial basis point method.Compared with the finite element method, the meshless method [21, 22] has the following main advantages: no grid, preprocessing, and simple adaptive analysis. It is easy to expand the low-dimensional problem to the high-dimensional problem and solve the problem that is difficult to solve by the traditional numerical methods. Therefore, it was a hot research at home and abroad in last twenty years and considered as a very promising numerical analysis method. Radial basis function meshless method [23] has two main forms: collocation method (strong form), Galerkin method (weak form). The radial basis function collocation method (RBCM) is a meshless method for solving the numerical solution of partial differential equations. The shape function formed by the interpolation of the radial basis functions (RBFs) satisfies the Kronecker Delta function property, so that the essential boundary condition can be applied very well [9]. For the high order boundary value problem, the Neumann boundary condition affects the numerical calculation result directly. However, the research results in [10–20] showed that the accuracy obtained by using direct collocation scheme is a bit poor especially on boundary. Therefore, in order to satisfy the Neumann boundary condition accurately, the essential boundary condition contains the description of the field function and the description of the derivative of the field function.Approximation of RBFs possesses optimal error estimation and convergence. However, the accuracy of the derivative of the interpolating function is usually very poor on the boundary when collocation method is used. In order to enhance the accuracy of the derivative of the approximation functions, Chu et al. [24] used Hermite radial basis point interpolation method for the nonuniform material gradient plane plate vibration; a meshless method was developed based on the collocation method and ridge basis function interpolation by Wang et al. [25]; Xing [9] constructed the governing equations of laminated plates with large deflection bending problem, approximated these field variables with RBF and HRBF, and discreted the control equations with least square collocation method; Wang et al. [26] proposed an improved meshless method which neither needs the computation of integrals nor requires a partition of the region and its boundary, and this method is applied to elliptic equations for examining its appropriateness; Parzlivand and Shahrezaee [27] mentioned a numerical technique which is a combination of collocation method, radial basis functions, the operational matrix of derivative for radial basis functions, and the new computational technique. Then they got the solution of a parabolic partial differential equation with a time-dependent coefficient subject to an extra measurement; Krowiak [28] proposed a Hermite radial basis function differential quadrature method for higher order equations and so on [29–32].The traditional radial basis point method has a large error at the Neumann boundary condition, so it cannot be used to deal with the Neumann boundary condition well. However, due to the existence of partial guidance in the Hermite radial basis function, it can reduce the error better. In this paper, the Hermite radial basis function combined with the collocation method was used to solve the unsaturated soil water movement equation with Neumann boundary condition. The Hermite radial basis function interpolation not only is satisfied with the Kronecker Delta function, but also can improve the precision of the interpolation function significantly. In addition, due to the nonlinear diffusion term of the second-order partial derivative term in the soil water equation, the chain rule was used to preprocess the nonlinear term of (1) to avoid the phenomenon that the collocation method cannot be used directly. The remainder of this paper is organized as follows: In Section 2, some preliminaries were provided regarding the Hermite radial basis collocation and nonlinear term processing methods. In Section 3, we introduced the Hermite radial basis collocation scheme for the unsaturated soil water movement equation with Neumann boundary condition. In Section 4, the numerical examples were used to examine the appropriateness and efficiency of the method proposed in this paper by comparing it with some traditional methods and applying it to the real background engineering problem. Conclusions are given at the end of the paper.


## 2. Hermite Radial Basis Collocation Method (HRBCM)

### 2.1. Radial Basis Function

RBF [33–37] is also called a distance basis function, which is a function of the distance r=x-xi from the node x to the node xi. The RBFs, which have the advantages of having a simple form, being independent of space dimension and isotropy, and being suitable for applying to the numerical approximation. In general, the RBFs can be divided into two categories. One is compacted supported radial basis function, and the other one is the global radial basis function. In contrast to the compacted supported function, the global radial basis function has better accuracy in the interpolation calculation, but the computational complexity is larger. It is a function that depends on the radial coordinate r:(2)φjr=φx-xj,where r=x-xj is the Euclidean distance from node x to node xj. Typical RBFs includemultiquadrics function:ϕIx=c2+r21/2,inverse multiquadrics function:ϕIx=c2+r2-1/2,Gaussians function:ϕIx=e-cr2,TPS (thin-plate splines):ϕIx=r2βlog⁡r,wherec is a constant (c>0) and β is an integer.


### 2.2. Hermite Radial Basis Function Approximation

Using the radial basis collocation method (RBCM), the derivative of the approximation function will produce a large error at the boundary of the region. Thus, RBCM cannot be applied to solve the differential equation with Neumann boundary condition, but HRBCM can be used to solve this problem well [22]. Hermite interpolation of the field quantity function u(x) in the local support region near the point x(x,z) is approximated by the RBFs and its derivative at the Neumann boundary which can be expressed as [22](3)ux≅u⌢x=∑i=1Nsaiϕi+∑j=1Nbbj∂ϕj∂n,where(4)ϕi=ϕx-xi,ϕj=ϕx-xj,∂ϕj∂n=ljx∂ϕj∂x+ljz∂ϕj∂z,where Ns is the number of nodes in the local support domain and Nb is the number of nodes on the Neumann boundary; let Ns+Nb=N; N represents the total number of nodes. xi is the interior points, xj is the boundary point, ϕi is the RBFs, ai is the undetermined coefficient of the corresponding RBFs, and bj is the number of undetermined coefficients corresponding to the derivatives of the RBFs on the Neumann boundary. ljx=cos(n,x) and ljz=cos(n,z) denoted the direction cosine in x and z directions on the Neumann boundary, respectively.We can obtain algebraic equations by interpolating at thekth point in the local support domain:(5)u^k≅u^xk=∑i=1Nsaiϕki+∑j=1Nbbj∂ϕkj∂n,k=1,2,…,Ns,where ϕki=ϕxk-xi, ϕkj=ϕxk-xj.We also can be obtain algebraic equation by interpolating at themth point on the Neumann boundaries:(6)∂u^m∂n=∂u^xm∂n=∑i=1Nsai∂ϕmi∂n+∑j=1Nbbj∂∂n∂ϕmj∂n,m=1,2,…,Nb.As shown the Figure1, “0” is a collocation point, “1”, “2,” and “3” are the Neumann boundary points, “4” and “5” are the Dirichlet boundary points, and the “hollow” points in the circle are the other collocation points used to construct the approximating function at the “0” point.Figure 1
Hermite collocation figure.


### 2.3. Hermite-Type Collocation

Many practical problems can be solved by the solution of partial differential equations with the initial and boundary conditions—a fixed solution problem. Consider a boundary value problem as follows:(7)Lux=fx,x∈Ω,Bux=gx,x∈Γ,where L is the differential operator defining the governing differential equations to be satisfied in the open domain Ω, B is the differential operator defining the boundary condition on boundary Γ,u is the problem unknown, f is the source term, and g is the term associated with boundary conditions. For the unknown u, the approximation represented by HRBFs can be expressed as(8)ux≈uhx=∑i=1Nsαiφix+∑j=1Nbαj∂∂nφjx=ΦTx,where φi(x) is the radical basis function centered at source point xi and αi is the expansion coefficient. Ns is the number of source points in the domain, and Nb is the number of source points on the boundary.Based on the main principle of the collocation method, the expressions can be obtained by putting (8) into (7) as follows:(9)∑j=1NLφjxkuj=fxk,x∈Ω,k=1,2,…,Ns,∑j=1NB∂∂nφjxkuj=gxk,x∈Γ,k=1,2,…,Nb,when Ns+Nb>N and (9) is the overdetermined system to be solved using the least squares method. When Ns+Nb=N, (9) can be summarized as (10)AU=F,where A is the coefficients matrix formed by the operators L and B under the action, U is the vector formed by the unknown quantity, and F is the source term and the term associated with boundary conditions formed by the vector.The basic principle of the collocation method [9, 22, 25] is to make the nodes satisfied with the differential equation in the domain and the nodes on the boundary satisfied with the boundary condition. Then the nonlinear equations about the initial-boundary value problem of differential equation are obtained.


## 2.1. Radial Basis Function

RBF [33–37] is also called a distance basis function, which is a function of the distance r=x-xi from the node x to the node xi. The RBFs, which have the advantages of having a simple form, being independent of space dimension and isotropy, and being suitable for applying to the numerical approximation. In general, the RBFs can be divided into two categories. One is compacted supported radial basis function, and the other one is the global radial basis function. In contrast to the compacted supported function, the global radial basis function has better accuracy in the interpolation calculation, but the computational complexity is larger. It is a function that depends on the radial coordinate r:(2)φjr=φx-xj,where r=x-xj is the Euclidean distance from node x to node xj. Typical RBFs includemultiquadrics function:ϕIx=c2+r21/2,inverse multiquadrics function:ϕIx=c2+r2-1/2,Gaussians function:ϕIx=e-cr2,TPS (thin-plate splines):ϕIx=r2βlog⁡r,wherec is a constant (c>0) and β is an integer.


## 2.2. Hermite Radial Basis Function Approximation

Using the radial basis collocation method (RBCM), the derivative of the approximation function will produce a large error at the boundary of the region. Thus, RBCM cannot be applied to solve the differential equation with Neumann boundary condition, but HRBCM can be used to solve this problem well [22]. Hermite interpolation of the field quantity function u(x) in the local support region near the point x(x,z) is approximated by the RBFs and its derivative at the Neumann boundary which can be expressed as [22](3)ux≅u⌢x=∑i=1Nsaiϕi+∑j=1Nbbj∂ϕj∂n,where(4)ϕi=ϕx-xi,ϕj=ϕx-xj,∂ϕj∂n=ljx∂ϕj∂x+ljz∂ϕj∂z,where Ns is the number of nodes in the local support domain and Nb is the number of nodes on the Neumann boundary; let Ns+Nb=N; N represents the total number of nodes. xi is the interior points, xj is the boundary point, ϕi is the RBFs, ai is the undetermined coefficient of the corresponding RBFs, and bj is the number of undetermined coefficients corresponding to the derivatives of the RBFs on the Neumann boundary. ljx=cos(n,x) and ljz=cos(n,z) denoted the direction cosine in x and z directions on the Neumann boundary, respectively.We can obtain algebraic equations by interpolating at thekth point in the local support domain:(5)u^k≅u^xk=∑i=1Nsaiϕki+∑j=1Nbbj∂ϕkj∂n,k=1,2,…,Ns,where ϕki=ϕxk-xi, ϕkj=ϕxk-xj.We also can be obtain algebraic equation by interpolating at themth point on the Neumann boundaries:(6)∂u^m∂n=∂u^xm∂n=∑i=1Nsai∂ϕmi∂n+∑j=1Nbbj∂∂n∂ϕmj∂n,m=1,2,…,Nb.As shown the Figure1, “0” is a collocation point, “1”, “2,” and “3” are the Neumann boundary points, “4” and “5” are the Dirichlet boundary points, and the “hollow” points in the circle are the other collocation points used to construct the approximating function at the “0” point.Figure 1
Hermite collocation figure.


## 2.3. Hermite-Type Collocation

Many practical problems can be solved by the solution of partial differential equations with the initial and boundary conditions—a fixed solution problem. Consider a boundary value problem as follows:(7)Lux=fx,x∈Ω,Bux=gx,x∈Γ,where L is the differential operator defining the governing differential equations to be satisfied in the open domain Ω, B is the differential operator defining the boundary condition on boundary Γ,u is the problem unknown, f is the source term, and g is the term associated with boundary conditions. For the unknown u, the approximation represented by HRBFs can be expressed as(8)ux≈uhx=∑i=1Nsαiφix+∑j=1Nbαj∂∂nφjx=ΦTx,where φi(x) is the radical basis function centered at source point xi and αi is the expansion coefficient. Ns is the number of source points in the domain, and Nb is the number of source points on the boundary.Based on the main principle of the collocation method, the expressions can be obtained by putting (8) into (7) as follows:(9)∑j=1NLφjxkuj=fxk,x∈Ω,k=1,2,…,Ns,∑j=1NB∂∂nφjxkuj=gxk,x∈Γ,k=1,2,…,Nb,when Ns+Nb>N and (9) is the overdetermined system to be solved using the least squares method. When Ns+Nb=N, (9) can be summarized as (10)AU=F,where A is the coefficients matrix formed by the operators L and B under the action, U is the vector formed by the unknown quantity, and F is the source term and the term associated with boundary conditions formed by the vector.The basic principle of the collocation method [9, 22, 25] is to make the nodes satisfied with the differential equation in the domain and the nodes on the boundary satisfied with the boundary condition. Then the nonlinear equations about the initial-boundary value problem of differential equation are obtained.


## 3. Hermite Radial Basis Collocation Scheme

### 3.1. Discretization Scheme

Considering the two-dimensional unsaturated vertical infiltration, the Richards’ equation and initial-boundary conditions are shown as follows:(11)∂θ∂t=∂∂xDθ∂θ∂x+∂∂zDθ∂θ∂z-∂Kθ∂z+fx,z,t,x,z∈Ω,t>0;θx,z,0=θ0,x,z∈Ω;θx,z,t=θs,x,z∈Γ1,t>0;∂θx,z,t∂n=θ1,x,z∈Γ,t>0,where θ is the soil water content, θ0 is the initial water content, D(θ) and K(θ) are the unsaturated soil water diffusivity and conductivity, f(x,z,t) is the root water absorption term, t is the time, x is the horizontal distance, z is the underground depth (positive upward), θ1 is the moisture content of the surface due to wet condition remain unchanged, Γ1 is the Dirichlet boundary, and Γ is the Neumann boundary.Step 1 (preprocessing the nonlinear term).
Because of the nonlinear term, the collocation method cannot be used directly. It needs to preprocess the nonlinear term, and the nonlinear diffusion term is considered as a composite function and resolved by the chain rule as follows:(12)∂θ∂t-Dθ·∂2θ∂x2+∂2θ∂z2-Aθ·∂θ∂x-Bθ·∂θ∂z=fx,z,t,where Aθ=∂D(θ)/∂θ·∂θ/∂x and Bθ=∂D(θ)/∂θ·∂θ/∂z-∂K(θ)/∂θ.Step 2 (constructing the approximation function).
By taking the Gaussian function as an example, the approximation function at the pointx(x,z) can be described as(13)φjri=e-cri2c>0,where ri=x-xi2+z-zi2 and ri represents the distance between x and xi. Thus, the approximate solution of the Richards’ equation can be expressed by Hermite radial basis function interpolation as follows:(14)θxj,tm≅∑i=1Nbai∂∂nφjri+∑i=Nb+1Nb+Nsbiφjri,where Nb is the number of boundary points, Ns is the number of inner points, and xj=(xj,zj) at the nodes numbered j.Step 3 (discrete the Richards’ equation).
The forward difference and collocation methods are used for the time discretization and space discretization, respectively. The discretization scheme for the Richards’ equation can be obtained as follows:(15)θjm+1-θjmΔt-Dθjm+1·∑i=1Nbai·∂3φjm+1ri∂x2∂n+∂3φjm+1ri∂z2∂n+∑i=1+NbNb+Nsbi·∂2φjm+1ri∂x2+∂2φjm+1ri∂z2x=xj,z=zj-Aθjm+1·∑i=1Nbai·∂2φjm+1ri∂x∂n+∑i=1+NbNb+Nsbi·∂φjm+1ri∂xx=xj,z=zj-Bθjm+1∑i=1Nbai∂2φjm+1ri∂z∂n+∑i=1+NbNb+Nsbi∂φjm+1ri∂zx=xj,z=zj=fjm+1.Collating (15), (16)
∑
i
=
1
N
b
a
i
m
+
1
·
∂
∂
n
φ
j
m
+
1
r
i
+
D
θ
j
m
+
1
·
Δ
t
·
∂
3
φ
j
m
+
1
r
i
∂
x
2
∂
n
+
∂
3
φ
j
m
+
1
r
i
∂
z
2
∂
n
-
Δ
t
·
A
θ
j
m
+
1
·
∂
2
φ
j
m
+
1
r
i
∂
x
∂
n
-
B
θ
j
m
+
1
·
∂
2
φ
j
m
+
1
r
i
∂
z
∂
n
+
∑
i
=
1
+
N
b
N
b
+
N
s
b
i
m
+
1
φ
j
m
+
1
r
i
+
D
θ
j
m
+
1
·
Δ
t
·
∂
2
φ
j
m
+
1
r
i
∂
x
2
+
∂
2
φ
j
m
+
1
r
i
∂
z
2
-
Δ
t
·
A
θ
j
m
+
1
·
∂
φ
j
m
+
1
r
i
∂
x
-
B
θ
j
m
+
1
·
∂
φ
j
m
+
1
r
i
∂
z
=
θ
j
m
+
Δ
t
·
f
j
m
+
1
,
j
=
N
b
+
1
,
…
,
N
s
+
N
b
,where Aθjm+1=dD(θjm+1)/dθ·∂θjm+1/∂x and B(θjm+1)=dD(θjm+1)/dθ·∂θjm+1/∂z-∂K(θjm+1)/∂θ.Step 4 (discrete the initial and boundary condition).
Initial Condition (17) θ x j , z j , 0 = θ 0 ,  j = 1,2 , … , N b + N s .
Boundary Condition
(  1) Dirichlet Boundary (18) θ x , z , t Γ = θ 1 ⟹ ∑ i = 1 N b α i m + 1 · φ j r i + ∑ i = N b N b i m + 1 · ∂ ∂ n φ j r i = ∑ i = 1 N b α i m + 1 · φ j r i - ∑ i = N b N b i m + 1 · ∂ ∂ x φ j r i = θ 1 m + 1 ,where ri represents the distance between x and xi.
(  2) Neumann Boundary (19) ∂ θ ∂ n x = x j , z = z j = θ 2 ⟹ l j x ∂ θ j b ∂ x + l j z ∂ θ j b ∂ z x = x j , z = z j = θ 2 ⟹ l j x · ∂ ∑ i = 1 N b α i m + 1 · φ j r i + ∑ i = N b + 1 N b i m + 1 · ∂ / ∂ n φ j r i ∂ x + l j z · ∂ ∑ i = 1 N b α i m + 1 · φ j r i + ∑ i = N b + 1 N b i m + 1 · ∂ / ∂ n φ j r i ∂ z = θ 2 m + 1 ,where ri represents the distance between xj and xi. ljx and ljz are the cosine of the Neumann boundary in the x and z directions.Step 5 (constructing nonlinear equations).
Combining (16), (17), (18), and (19), the matrix equation can be written as follows: (20)Hα=F,where(21)H=∂2φ1r1∂y2⋯∂2φ1rn∂y2∂φ1rn+1∂y⋯∂φ1rNb∂y∂φ1rNb+1∂y⋯∂φ1rN∂y⋮⋮⋮⋮⋮⋮⋮⋮⋮∂2φnr1∂y2⋯∂2φnrn∂y2∂φnrn+1∂y⋯∂φnrNb∂y∂φnrNb+1∂y⋯∂φnrN∂y∂φn+1r1∂y⋯∂φn+1rn∂yφn+1rn+1⋯φn+1rNbφn+1rNb+1⋯φn+1rN⋮⋮⋮⋮⋮⋮⋮⋮⋮∂φNbr1∂y⋯∂φNbrn∂yφNbrn+1⋯φNbrNbφNbrNb+1⋯φNbrNΦNb+1r1⋯ΦNb+1rnΨNb+1n+1rn+1⋯ΨNb+1NrNbΨNb+1N+1rNb+1⋯ΨNb+1NrN⋮⋮⋮⋮⋮⋮⋮⋮⋮ΦNr1⋯ΦNrnΨNr1+n⋯ΨNrNbΨNrNb+1⋯ΨNrNα=α1m+1,…,αNbm+1,bNb+1m+1,…,bNm+1Τ,F=θ1m,…,θNbm,θNb+1m+Δt·fNb+1m+1,…,θNm+Δt·fNm+1Τ,where(22)Φjri=∑i=1Nb∂∂nφjm+1ri+Dθjm+1·Δt·∂3φjm+1ri∂x2∂n+∂3φjm+1ri∂z2∂n-Δt·Aθjm+1·∂2φjm+1ri∂x∂n-Bθjm+1·∂2φjm+1ri∂z∂n,Ψjri=φjm+1ri+Dθjm+1·Δt·∂2φjm+1ri∂x2+∂2φjm+1ri∂z2-Δt·Aθjm+1·∂φjm+1ri∂x-Bθjm+1·∂φjm+1ri∂z,where Aθjm+1=dD(θjm+1)/dθ·∂θjm+1/∂x and B(θjm+1)=dD(θjm+1)/dθ·∂θjm+1/∂z-∂K(θjm+1)/∂θ.


### 3.2. Linearization Scheme

The iterative method is used to solve systems of the nonlinear equations. The common method of solving nonlinear equations are direct iteration method, Newton-Raphson method, modified Newton-Raphson method, incremental method, arc length method, and so on. A predictor correction method is adopted in this paper. The specific steps are as follows:AssumingK(θ), D(θ) are nonlinear term and tn+1 moment, let K(θ(xi,zj,tn+1))=K(θin+1), D(θ(xi,zj,tn+1))=D(θin+1); use θi-1n, θin, θi+1nas the forecast value of θi-1n+1, θin+1, θi+1n+1; getK(θi-1n),D(θin),K(θi+1n); make it as the forecast value ofK(θi-1n+1),D(θin+1),K(θi+1n+1); the solution θi-1n+1l(i=1,2,3,…,N) is solved by the first iteration; then getK(θi-1n+11),D(θin+1(1)),K(θi+1n+11) separately, and make it as the correction value ofK(θi-1n+1),D(θin+1),K(θi+1n+1) and as the next forecast value of θin+12(i=1,2,3,…,N). Repeat the above steps until the difference between the two iterations before and after the end of the solution is less than the specified allowable error, which is(23)maxiθin+1q-θin+1q-1θin+1q-1≤e,i=1,2,…,N,where q is the number of iterations and e is the allowable relative error, which is determined by the exact requirement of the calculation; for example, e=10-4. If it satisfied (23), we can replace θin+1q-1 with θin+1qj=1,2,…,N.


## 3.1. Discretization Scheme

Considering the two-dimensional unsaturated vertical infiltration, the Richards’ equation and initial-boundary conditions are shown as follows:(11)∂θ∂t=∂∂xDθ∂θ∂x+∂∂zDθ∂θ∂z-∂Kθ∂z+fx,z,t,x,z∈Ω,t>0;θx,z,0=θ0,x,z∈Ω;θx,z,t=θs,x,z∈Γ1,t>0;∂θx,z,t∂n=θ1,x,z∈Γ,t>0,where θ is the soil water content, θ0 is the initial water content, D(θ) and K(θ) are the unsaturated soil water diffusivity and conductivity, f(x,z,t) is the root water absorption term, t is the time, x is the horizontal distance, z is the underground depth (positive upward), θ1 is the moisture content of the surface due to wet condition remain unchanged, Γ1 is the Dirichlet boundary, and Γ is the Neumann boundary.Step 1 (preprocessing the nonlinear term).
Because of the nonlinear term, the collocation method cannot be used directly. It needs to preprocess the nonlinear term, and the nonlinear diffusion term is considered as a composite function and resolved by the chain rule as follows:(12)∂θ∂t-Dθ·∂2θ∂x2+∂2θ∂z2-Aθ·∂θ∂x-Bθ·∂θ∂z=fx,z,t,where Aθ=∂D(θ)/∂θ·∂θ/∂x and Bθ=∂D(θ)/∂θ·∂θ/∂z-∂K(θ)/∂θ.Step 2 (constructing the approximation function).
By taking the Gaussian function as an example, the approximation function at the pointx(x,z) can be described as(13)φjri=e-cri2c>0,where ri=x-xi2+z-zi2 and ri represents the distance between x and xi. Thus, the approximate solution of the Richards’ equation can be expressed by Hermite radial basis function interpolation as follows:(14)θxj,tm≅∑i=1Nbai∂∂nφjri+∑i=Nb+1Nb+Nsbiφjri,where Nb is the number of boundary points, Ns is the number of inner points, and xj=(xj,zj) at the nodes numbered j.Step 3 (discrete the Richards’ equation).
The forward difference and collocation methods are used for the time discretization and space discretization, respectively. The discretization scheme for the Richards’ equation can be obtained as follows:(15)θjm+1-θjmΔt-Dθjm+1·∑i=1Nbai·∂3φjm+1ri∂x2∂n+∂3φjm+1ri∂z2∂n+∑i=1+NbNb+Nsbi·∂2φjm+1ri∂x2+∂2φjm+1ri∂z2x=xj,z=zj-Aθjm+1·∑i=1Nbai·∂2φjm+1ri∂x∂n+∑i=1+NbNb+Nsbi·∂φjm+1ri∂xx=xj,z=zj-Bθjm+1∑i=1Nbai∂2φjm+1ri∂z∂n+∑i=1+NbNb+Nsbi∂φjm+1ri∂zx=xj,z=zj=fjm+1.Collating (15), (16)
∑
i
=
1
N
b
a
i
m
+
1
·
∂
∂
n
φ
j
m
+
1
r
i
+
D
θ
j
m
+
1
·
Δ
t
·
∂
3
φ
j
m
+
1
r
i
∂
x
2
∂
n
+
∂
3
φ
j
m
+
1
r
i
∂
z
2
∂
n
-
Δ
t
·
A
θ
j
m
+
1
·
∂
2
φ
j
m
+
1
r
i
∂
x
∂
n
-
B
θ
j
m
+
1
·
∂
2
φ
j
m
+
1
r
i
∂
z
∂
n
+
∑
i
=
1
+
N
b
N
b
+
N
s
b
i
m
+
1
φ
j
m
+
1
r
i
+
D
θ
j
m
+
1
·
Δ
t
·
∂
2
φ
j
m
+
1
r
i
∂
x
2
+
∂
2
φ
j
m
+
1
r
i
∂
z
2
-
Δ
t
·
A
θ
j
m
+
1
·
∂
φ
j
m
+
1
r
i
∂
x
-
B
θ
j
m
+
1
·
∂
φ
j
m
+
1
r
i
∂
z
=
θ
j
m
+
Δ
t
·
f
j
m
+
1
,
j
=
N
b
+
1
,
…
,
N
s
+
N
b
,where Aθjm+1=dD(θjm+1)/dθ·∂θjm+1/∂x and B(θjm+1)=dD(θjm+1)/dθ·∂θjm+1/∂z-∂K(θjm+1)/∂θ.Step 4 (discrete the initial and boundary condition).
Initial Condition (17) θ x j , z j , 0 = θ 0 ,  j = 1,2 , … , N b + N s .
Boundary Condition
(  1) Dirichlet Boundary (18) θ x , z , t Γ = θ 1 ⟹ ∑ i = 1 N b α i m + 1 · φ j r i + ∑ i = N b N b i m + 1 · ∂ ∂ n φ j r i = ∑ i = 1 N b α i m + 1 · φ j r i - ∑ i = N b N b i m + 1 · ∂ ∂ x φ j r i = θ 1 m + 1 ,where ri represents the distance between x and xi.
(  2) Neumann Boundary (19) ∂ θ ∂ n x = x j , z = z j = θ 2 ⟹ l j x ∂ θ j b ∂ x + l j z ∂ θ j b ∂ z x = x j , z = z j = θ 2 ⟹ l j x · ∂ ∑ i = 1 N b α i m + 1 · φ j r i + ∑ i = N b + 1 N b i m + 1 · ∂ / ∂ n φ j r i ∂ x + l j z · ∂ ∑ i = 1 N b α i m + 1 · φ j r i + ∑ i = N b + 1 N b i m + 1 · ∂ / ∂ n φ j r i ∂ z = θ 2 m + 1 ,where ri represents the distance between xj and xi. ljx and ljz are the cosine of the Neumann boundary in the x and z directions.Step 5 (constructing nonlinear equations).
Combining (16), (17), (18), and (19), the matrix equation can be written as follows: (20)Hα=F,where(21)H=∂2φ1r1∂y2⋯∂2φ1rn∂y2∂φ1rn+1∂y⋯∂φ1rNb∂y∂φ1rNb+1∂y⋯∂φ1rN∂y⋮⋮⋮⋮⋮⋮⋮⋮⋮∂2φnr1∂y2⋯∂2φnrn∂y2∂φnrn+1∂y⋯∂φnrNb∂y∂φnrNb+1∂y⋯∂φnrN∂y∂φn+1r1∂y⋯∂φn+1rn∂yφn+1rn+1⋯φn+1rNbφn+1rNb+1⋯φn+1rN⋮⋮⋮⋮⋮⋮⋮⋮⋮∂φNbr1∂y⋯∂φNbrn∂yφNbrn+1⋯φNbrNbφNbrNb+1⋯φNbrNΦNb+1r1⋯ΦNb+1rnΨNb+1n+1rn+1⋯ΨNb+1NrNbΨNb+1N+1rNb+1⋯ΨNb+1NrN⋮⋮⋮⋮⋮⋮⋮⋮⋮ΦNr1⋯ΦNrnΨNr1+n⋯ΨNrNbΨNrNb+1⋯ΨNrNα=α1m+1,…,αNbm+1,bNb+1m+1,…,bNm+1Τ,F=θ1m,…,θNbm,θNb+1m+Δt·fNb+1m+1,…,θNm+Δt·fNm+1Τ,where(22)Φjri=∑i=1Nb∂∂nφjm+1ri+Dθjm+1·Δt·∂3φjm+1ri∂x2∂n+∂3φjm+1ri∂z2∂n-Δt·Aθjm+1·∂2φjm+1ri∂x∂n-Bθjm+1·∂2φjm+1ri∂z∂n,Ψjri=φjm+1ri+Dθjm+1·Δt·∂2φjm+1ri∂x2+∂2φjm+1ri∂z2-Δt·Aθjm+1·∂φjm+1ri∂x-Bθjm+1·∂φjm+1ri∂z,where Aθjm+1=dD(θjm+1)/dθ·∂θjm+1/∂x and B(θjm+1)=dD(θjm+1)/dθ·∂θjm+1/∂z-∂K(θjm+1)/∂θ.


## 3.2. Linearization Scheme

The iterative method is used to solve systems of the nonlinear equations. The common method of solving nonlinear equations are direct iteration method, Newton-Raphson method, modified Newton-Raphson method, incremental method, arc length method, and so on. A predictor correction method is adopted in this paper. The specific steps are as follows:AssumingK(θ), D(θ) are nonlinear term and tn+1 moment, let K(θ(xi,zj,tn+1))=K(θin+1), D(θ(xi,zj,tn+1))=D(θin+1); use θi-1n, θin, θi+1nas the forecast value of θi-1n+1, θin+1, θi+1n+1; getK(θi-1n),D(θin),K(θi+1n); make it as the forecast value ofK(θi-1n+1),D(θin+1),K(θi+1n+1); the solution θi-1n+1l(i=1,2,3,…,N) is solved by the first iteration; then getK(θi-1n+11),D(θin+1(1)),K(θi+1n+11) separately, and make it as the correction value ofK(θi-1n+1),D(θin+1),K(θi+1n+1) and as the next forecast value of θin+12(i=1,2,3,…,N). Repeat the above steps until the difference between the two iterations before and after the end of the solution is less than the specified allowable error, which is(23)maxiθin+1q-θin+1q-1θin+1q-1≤e,i=1,2,…,N,where q is the number of iterations and e is the allowable relative error, which is determined by the exact requirement of the calculation; for example, e=10-4. If it satisfied (23), we can replace θin+1q-1 with θin+1qj=1,2,…,N.


## 4. Results and Discussion

Example 1 (one-dimensional unsaturated soil water movement equation).
(24) ∂ θ ∂ t = ∂ ∂ z D θ ∂ θ ∂ z + ∂ K θ ∂ z + f z , t ,  0 ≤ z ≤ 1 ,   t > 0 ; θ z , 0 = θ 0 z ,  t = 0 ,   0 ≤ z ≤ 1 ; θ 0 , t = g 1 ,  z = 0 ,   t > 0 ; ∂ ∂ n θ 1 , t = g N ,  z = 1 ,   t > 0 ,where K(θ)=θ, D(θ)=0.001θ+0.001, exact solution θ=e-t(z-z2), the root water uptake item f(z,t), and the initial value θ0(z) and boundary value g1, g2 are determined by exact solution. Gaussian function φ(r)=e-cr2(c>0) was selected as the radial function. In this paper, we discuss the error of L2-norm: (25)θ-θhL2=∫Ωθ-θh2dΩ1/2=∑i=1Nθi-θih2∑i=1Nθi2.The comparative result of calculation accuracy about different algorithm is shown in Figures 2, 3 and Tables 1, 2.Table 1
Number of nodes is 11 and 51; the comparison of relative error about different algorithm.
Nodes number Time step ParameterC  (HRBCM/RBCM) Time HRBCM error FDM error RBCM error 11 0.01 2.78/1.6 0.1 2.7507 e - 05 1.0850 e - 02 9.3827 e - 04 0.5 2.5677 e - 04 1.1018 - 02 3.5145 e - 03 1 6.3712 e - 03 1.0850 e - 02 4.6151 e - 02 0.05 1.8/1.6 0.1 5.3214 e - 04 1.0489 e - 02 1.9501 e - 03 0.5 1.6053 e - 03 1.0381 e - 02 2.2310 e - 02 1 4.2426 e - 02 1.0344 e - 02 6.6038 e - 02 0.1 2.1/1.6 0.1 3.9982 e - 03 9.7463 e - 03 3.2716 e - 03 0.5 2.6891 e - 02 9.7037 e - 03 7.5105 e - 02 1 3.5269 e - 02 9.7029 e - 03 1.3290 e - 01 51 0.01 225/205 0.1 7.1412 e - 05 1.1087 e - 02 2.2450 e - 03 0.5 3.1904 e - 04 1.0933 e - 02 1.2102 e - 02 1 8.4911 e - 03 1.0876 e - 02 3.3314 e - 02 0.05 226/205 0.1 1.2680  e - 04 1.0537 e - 02 2.2680 e - 03 0.5 2.6000  e - 03 1.0524 e - 02 7.3873 e - 03 1 1.4026  e - 02 1.0451 e - 02 2.8443 e - 02 0.1 215/205 0.1 5.1545 e - 04 9.9190 e - 03 5.5687 e - 03 0.5 3.1763 e - 03 1.0384 e - 02 2.7842 e - 02 1 5.9878 e - 03 1.1884 e - 02 7.5204 e - 02Table 2
The comparison of relative error about different algorithm under different number of nodes.
Time step The number of nodes ParameterC (HRBCM/RBCM) Time HRBCM error FDM error RBCM error 0.01 11 2.75/1.6 0.1 2.7507 e - 05 1.0850 e - 02 9.3827 e - 04 0.5 2.5677 e - 04 1.1018 - 02 3.5145 e - 03 1 6.3712 e - 03 1.0850 e - 02 4.6151 e - 02 21 25/30 0.1 3.7880 e - 05 1.1078 e - 02 2.7659 e - 03 0.5 2.9703 e - 04 1.0924 e - 02 1.1632 e - 02 1 8.1927 e - 03 1.0867 e - 02 5.8906 e - 02 51 225/205 0.1 7.1412 e - 05 1.1087 e - 02 2.2450 e - 03 0.5 3.1904 e - 04 1.0933 e - 02 1.2102 e - 02 1 8.4911 e - 03 1.0876 e - 02 3.3314 e - 02 FDM: finite difference method, RBCM: radial basis collocation method.Figure 2
The comparison between numerical solution and exact solution about different algorithm inT=0.5 time.
(a)
Number of nodes is 21; time step is 0.01 (b)
Number of nodes is 51; time step is 0.01 (c)
Number of nodes is 51; time step is 0.01 (d)
Number of nodes is 51; time step is 0.01Figure 3
The error results between numerical solution and exact solution about different algorithm inT=0.5 time.
(a)
Number of nodes is 21; time step is 0.01 (b)
Number of nodes is 51; time step is 0.01 (c)
Number of nodes is 51; time step is 0.01 (d)
Number of nodes is 51; time step is 0.01It can be seen from Figures2, 3 and Tables 1, 2 that the algorithm of HRBCM has the highest calculation accuracy among the other methods for the Neumann boundary problem. It is obvious that the parameter C will increase with increasing the space nodes when we choose the different time steps, space distributions, and parameters C. From Table 1, we can see that the HRBCM error is better. For the same time steps 0.01 and the number of spatial nodes 11, 21, and 51, the HRBCM error, RBCM error, and FDM error range from 2.5677 × 10−4 to 3.1904 × 10−4, from 3.5145 × 10−3 to 2.1202 × 10−2, and from 1.0018 × 10−2 to 1.1933 × 10−2 when the computation time is 0.5, respectively. The RBCM and the FDM have larger error at the Neumann boundary condition, and the HRBCM reduces the error at the Neumann boundary condition (see Figures 3(b), 3(c), and 3(d)). Therefore, this method is relatively superior, and the error accuracy is higher compared with other algorithms.Example 2 (simulation of the soil water infiltration process).
(26) ∂ θ ∂ t = ∂ ∂ z D θ ∂ θ ∂ z + ∂ K θ ∂ z , θ z , 0 = θ 0 z , θ 0 , t = θ s , ∂ ∂ n θ L , t = 0 .Two most common infiltration models used to calculate the unsaturated conductivityK(θ) and diffusivity D(θ) in practice are Brooks-Corey (B-C) model and van Genuchten (V-G) model [1, 38, 39]. According to the B-C model and V-G model, we can obtain the expressions of nonlinear terms (27) and (28).Brooks-Corey Model (27) K θ = k s θ - θ r θ s - θ r m / n D θ = K θ · d h d θ = k s · h d · - 1 n · θ s - θ r 1 - m / n · θ - θ r m - 1 - n / n .Van Genuchten Model (28) K θ = k s θ - θ r θ s - θ r 1 / 2 1 - 1 - θ - θ r θ s - θ r 1 / m m 2 D θ = K θ · d h d θ = k s 1 - m θ - θ r / θ s - θ r 1 / 2 - 1 / m α m θ s - θ r · 1 - θ - θ r θ s - θ r 1 / m - m + 1 - θ - θ r θ s - θ r 1 / m m - 2 ,where θ-θr/θs-θr=hd/hn in B-C model, θ-θr/θs-θr=1/1+(αh)nm in V-G model, n is an index parameter representing the soil’s pore size distribution, m=2+3n in B-C model, m=1-1/n in V-G model, and hd=1/α is bubbling pressure (cm).The HYDRUS-1D 4.xx [38, 39] software was used to solve (26), and the relevant data pertaining to the infiltration of one-dimensional soil under a fixed water head was simulated. Table 3 shows the values of the physical parameters for two soils. Figure 4 shows the simulated results of the one-dimensional vertical (Figures 4(b) and 4(d)) and horizontal (Figures 4(a) and 4(c)) infiltration for the two soils at different times.Table 3
Hydraulic property parameters of two types of soil.
Soil θ r (cm3⋅min−3) θ s (cm3⋅min−3) θ 0 (cm3⋅min−3) K s (cm⋅min−1) α (cm) n Sandy loam (V-G) 0.065 0.41 0.18 0.0736806 0.075 1.89 Silt (B-C) 0.015 0.501 0.15 0.0113335 0.0482 0.211Figure 4
Time step is 1, comparison of numerical solution and simulation solution of HYDRUS software about the sandy loam and silt at different times.
(a)
The horizontal infiltration (b)
The vertical infiltration (c)
The horizontal infiltration (d)
The vertical infiltrationFigures4(a), 4(b), 4(c), and 4(d) are the infiltration processes of the B-C model and V-G model in the horizontal and vertical directions, respectively. According to the L2-norm error, we can get the L2-norm error about different models in different wetting front. The L2-norm error for the vertical infiltration was 3.8396% and 5.0855% when the wetting front depth is 23.5 cm and 48 cm, and it was 2.7664% and 1.1824% for the horizontal infiltration when wetting front depth is 20 cm and 40 cm for silt soil, respectively. Similarly, the L2-norm error for the vertical infiltration was 4.5637% and 1.6325% when the wetting front depth is 20 cm and 40 cm, and it was 2.9743% and 2.1441% for the horizontal infiltration when the wetting front depth is 25 cm and 35 cm for sandy loam soil, respectively.From Figures2, 3, and 4, we can see that the new method proposed in this paper adapts to the numerical examples, and it is feasible and effective in simulating the actual background problem within the allowable range of error. So the new method proposed in this paper can be applied to solve practical one-dimensional infiltration problems in engineering.Example 3 (two-dimensional unsaturated soil water movement equation in regular region).
(29) ∂ θ ∂ t = ∂ ∂ x D θ ∂ θ ∂ x + ∂ ∂ z D θ ∂ θ ∂ z - ∂ K θ ∂ z + f x , z , t ,  x , z ∈ Ω ,   Ω = 0,1 × 0,1 ; θ x , z , 0 = θ 0 ,  x , z ∈ Ω ; θ x , z , t = g 1 ,  x = 1 ; θ x , z , t = g 2 ,  x = 0 ; θ x , z , t = g 3 ,  z = 1 ; ∂ ∂ z θ x , z , t = g 4 ,  z = 0 ,where K(θ)=θ, D(θ)=0.01θ+0.01, exact solution θ=et(1+x)(1+z), root water absorption item f(x,z,t), initial value θ0(z), and boundary value g1, g2 are determined by exact solution. The radial basis function chooses the Gaussian function φ(r)=e-cr2(c>0), and the error was estimated using the L2-norm. Figure 4 gives the comparison result of numerical solution and exact solution in T=1 time about different algorithms. Table 3 gives the error accuracy of this paper in T=0.5 and T=1 time when the time step is 0.01.The main advantage of the nonmesh method is that less points were used as the collocation points to calculate the numerical solutions at all nodes in the whole region. It can be seen from Figure5 that the numerical solution and the exact solution fit are better and the error precision is well high. With the increase of the spatial nodes, the parameter C increases.Figure 5
Time step is 0.1, comparative results of numerical solution and exact solution inT=1 time.
(a)
Number of nodes is 441 in this method (b)
Number of nodes is 2551 in this method (c)
Comparison of FDM and RBCMTable4 shows the optimal numerical error. It can be seen that the HRBCM error is 9.0137 × 10−6 when the collocation points number is 16 and the calculated time is 1. Thus, the numerical solutions at the more nodes were simulated by the 16th collocation points, and the numerical errors were shown in Table 5. For the same time steps 0.01 and the number of spatial nodes 16, 36, 121, and 441, the HRBCM error, RBCM error, and FDM error range from 8.1384 × 10−6 to 9.4099 × 10−6, from 6.1375 × 10−5 to 7.3059 × 10−4, and from 3.0860 × 10−3 to 3.6695 × 10−2 when the computation time is 0.5, respectively. Therefore, the method of this paper is relatively superior, and the error accuracy is higher compared with other algorithm. The optimal numerical error was observed by choosing the optimal collocation points, the appropriate time step, and parameter C.Table 4
Time step is 0.01, error accuracy of this paper inT=0.5 and T=1 time with the same collocation points number and the calculation points number.
The collocation points number ParameterC Time HRBCM error 16 0.007 0.5 4.1384 e - 06 1 9.0137 e - 06 25 0.022 0.5 4.4073 e - 05 1 9.0726 e - 05 36 1.85 0.5 2.2788 e - 04 1 4.2127 e - 03 81 10 0.5 9.4896 e - 03 1 2.0206 e - 02 121 35.6 0.5 8.6331 e - 03 1 2.8443 e - 02 256 106.4 0.5 9.1025 e - 03 1 3.1263 e - 02 441 399.8 0.5 7.2096 e - 03 1 3.4099 e - 02Table 5
Comparison of the results of this paper and other algorithms using the 16 collocation points in Table4 to construct the approximate function to calculate more node values.
The calculation points number Time step ParameterC Time HRBCM error FDM error RBCM error 16 0.01 0.0046 0.1 2.0432 e - 06 7.4438 e - 04 7.4283 e - 06 0.5 8.1384 e - 06 3.0860 e - 03 6.1375 e - 05 1 9.0137 e - 06 4.6655 e - 03 2.5356 e - 04 0.1 0.007 0.1 7.9557 e - 06 7.2009 e - 03 2.4414 e - 05 0.5 6.6330 e - 06 3.0898 e - 02 4.5256 e - 05 1 5.1927 e - 06 4.7364 e - 02 7.0511 e - 05 36 0.01 0.0046 0.1 2.2891 e - 06 8.5735 e - 04 7.7374 e - 06 0.5 8.7090 e - 06 3.4395 e - 03 2.5466 e - 04 1 3.3153 e - 05 4.9781 e - 03 7.1789 e - 04 0.1 0.007 0.1 5.4913 e - 06 8.4151 e - 03 2.5976 e - 05 0.5 9.1136 e - 06 3.4894 e - 02 4.7182 e - 05 1 7.7744 e - 06 4.9995 e - 02 6.4569 e - 05 121 0.01 0.0046 0.1 2.1382 e - 06 8.9728 e - 04 3.3954 e - 06 0.5 9.1683 e - 06 3.5253 e - 03 7.2631 e - 05 1 3.3885 e - 05 4.9750 e - 03 2.5684 e - 04 0.1 0.007 0.1 5.5878 e - 06 8.8523 e - 03 2.9172 e - 05 0.5 1.0040 e - 05 3.6006 e - 02 4.6436 e - 05 1 8.6827 e - 05 5.0086 e - 02 5.1277 e - 05 441 0.01 0.0046 0.1 2.2096 e - 06 9.7449 e - 04 4.6080 e - 06 0.5 9.4099 e - 06 3.5397 e - 03 7.3059 e - 04 1 3.4266 e - 05 4.6280 e - 03 2.5792 e - 03 0.1 0.007 0.1 5.6725 e - 06 9.7607 e - 03 3.0843 e - 05 0.5 6.7735 e - 06 3.6695 e - 02 5.3469 e - 05 1 9.1582 e - 06 5.1487 e - 02 7.1210 e - 05According to the data in Table4, we can fit the relationship between parameter C and the number of nodes N and give the empirical formula as follows:(30)C=2×10-6×N6.665/2.It can be seen from the above equation that, with the increase of the spatial nodesN, the parameter C increases.Example 4 (two-dimensional unsaturated soil water movement equation in regular regional—point source infiltration problem).
(31) ∂ θ ∂ t = ∂ ∂ x D θ ∂ θ ∂ x + ∂ ∂ z D θ ∂ θ ∂ z - ∂ K θ ∂ z , θ x , z , 0 = θ 0 ,  x , z ∈ 0,100 × - 100,0 ; θ x , z , t = θ s ,  x = 0 ,   z = 0 ; ∂ ∂ z θ x , z , t = 0 ,  x , z ∈ 0,100 × - 100,0 .Similarly, according to the Brooks-Corey model about (31), we can also get the actual process of the two-dimensional soil water content infiltration about different type soils. Figure 6 gives the single point infiltration of two-dimensional unsaturated soil water movement equation when the wetting front is 48 cm. The physical parameters values [38, 39] of the sandy loam soils are shown as follows: θr=0.041,  θs=0.412,  n=0.322,  ks=0.04317, and  hd=14.66.Figure 6
Single point infiltration of two-dimensional unsaturated soil water movement equation.
(a)
HRBCM (b)
HYDRUS softwareIn order to verify the effectiveness of the method in this paper, we obtained the cumulative infiltration according to the given model parameters through the HYDRUS software to simulate the single point infiltration process. The cumulative infiltration can be obtained by the double integral of the water content calculated in this paper.(32)I=∫-1000∫0100θdxdz.The double integral of the above formula can be calculated by complex trapezoidal formula. The cumulative infiltration volumes of (32) and the HYDRUS software were 1.1859 × 103 and 1.1646 × 103 when the vertical wetting front was 48 cm. The error is very small, which is 1.828%.From Figures5 and 6, we can see that the algorithm proposed in this paper is also suitable for two-dimensional actual background problem within the allowable range of error. Therefore, the algorithm proposed in this paper can be applied to solve practical engineering problems.Example 5 (two-dimensional unsaturated soil water movement equation in irregular regional).
(33) ∂ θ ∂ t - ∂ ∂ x D θ ∂ θ ∂ x - ∂ ∂ z D θ ∂ θ ∂ z + ∂ K θ ∂ z = f x , z , t ,  x , z ∈ Ω ,   Ω = 0,1 × 0,1 ; θ x , z , 0 = g 0 ,  x , z ∈ Ω ; ∂ ∂ n θ x , z , t = g 1 ,  x = 0 ; ∂ ∂ n θ x , z , t = g 2 ,  z = 0 ; θ x , z , t = g 3 ,  z = - tan ⁡ α · x - 1 ,   0 ≤ x ≤ 1 ,where K(θ)=θ, D(θ)=0.01θ+0.01, exact solution θ=(1+t)(x-x2)(z-z2), root water absorption item f(x,z,t), initial value θ0(z), and boundary value g1, g2 are determined by exact solution. The radial basis function chooses the Gaussian function φ(r)=e-cr2(c>0), α=π/4, and the error is estimated using the L2-norm. Figure 7 shows the typical distribution of a triangle area, which is the sparse distribution on the left side and the dense distribution based on the actual situation of soil water infiltration in the top of the slope. Comparison results of this paper and other algorithms are shown in Figures 8, 9 and Table 6.Table 6
Comparison of the results of this paper and other algorithms using the 16 collocation points to construct the approximate function to calculate more node values.
The calculation points number Time step ParameterC Time RBCM error FDM error RBCM error 16 0.01 4 0.1 7.8912 e - 04 2.1524 e - 02 4.0892 e - 03 0.5 3.2868 e - 03 2.5384 e - 02 5.8206 e - 02 1 7.7251 e - 03 2.7810 e - 02 1.6605 e - 01 0.1 4.5 0.1 4.1504 e - 03 1.0864 e - 01 1.7071 e - 02 0.5 8.2285 e - 03 2.8542 e - 01 7.1026 e - 02 1 1.0176 e - 02 2.7239 e - 01 1.0936 e - 01 55 0.01 4 0.1 6.0602 e - 04 1.8701 e - 02 4.8589 e - 03 0.5 2.6782 e - 03 2.3108 e - 02 8.2325 e - 02 1 9.4395 e - 03 2.5500 e - 02 1.1613 e - 01 0.1 4.5 0.1 3.2473 e - 03 3.4521 e - 01 1.1552 e - 02 0.5 9.7611 e - 03 5.6369 e - 01 4.7614 e - 01 1 2.4074 e - 02 6.9379 e - 01 8.8699 e - 01 155 0.01 4 0.1 6.6656 e - 04 1.7437 e - 02 1.6260 e - 03 0.5 5.1966 e - 03 2.1955 e - 02 9.6772 e - 02 1 1.5110 e - 02 2.4220 e - 02 1.8763 e - 01 0.1 4.5 0.1 2.8588 e - 03 1.7646 e - 01 4.3647 e - 03 0.5 8.6347 e - 03 2.3540 e - 01 7.9979 e - 02 1 2.1302 e - 02 4.1907 e - 01 3.9264 e - 01 1275 0.01 4 0.1 5.4615 e - 04 2.9934 e - 02 1.2475 e - 03 0.5 5.5338 e - 03 2.1229 e - 02 5.4090 e - 03 1 9 .  2195 e - 03 2.3499 e - 02 9.7695 e - 02 0.1 4.5 0.1 2.6199 e - 03 1.6632 e - 01 4.0319 e - 03 0.5 7.9396 e - 03 2.9909 e - 01 7.7675 e - 02 1 1.9591 e - 02 2.2210 e - 01 2.6655 e - 01Figure 7
The stationing ways of triangle area.Similarly, as shown in Table4, the optimal numerical error was observed when the collocation points number was 16. Thus, the numerical solutions at the more nodes were simulated by the 16th collocation points and the numerical errors were shown in Table 6. For the same time steps 0.01 and the number of spatial nodes 16, 55, 155, and 1275, the HRBCM error, RBCM error, and FDM error range from 3.2868 × 10−3 to 7.9396 × 10−3, from 5.8206 × 10−2 to 9.7695 × 10−2, and from 2.5384 × 10−2 to 2.1299 × 10−2 when the time is 0.5, respectively.The comparison results between the RBCM and FDM in Table6 and Figure 9 have shown that the algorithm in this paper has better effect in dealing with Neumann boundary condition. Therefore, it is convenient to save time for simulating the soil moisture equation with Neumann boundary condition when we choose optimal collocations points, the appropriate time step, and parameterC.We can learn the calculation accuracy at more nodes, such as 55, 155, 1275, which are the same as 16 points when the 16 points are used as collocation points from Table6. As can be seen from Figures 8 and 9, the numerical solution is in good agreement with the exact solution, and the error precision is well high. The time step has little influence on the numerical results when the parameter C of the basis function remains unchanged and the spatial node increased. The error accuracy is the same at the same moment.Figure 8
Δ t = 0.01, comparison of numerical solution and exact solution in T=0.2.
(a)
Number of nodes is 155 (b)
Number of nodes is 1275Figure 9
Δ t = 0.01; number of nodes is 1275, comparison of RBCM and FDM in T=0.2.


## 5. Conclusions

Since soil water movement equation is a nonlinear partial differential equation, the existence of nonlinear diffusion terms often leads to the collocation method is hard to be used in discreting partial differential equations directly. This paper focused on the nonlinear problem of the Neumann boundary and used the chaining principle of the composite function derivation to preprocess the nonlinear term in the equation firstly. Then the HRBCM was used to construct algorithm for the Richards’ equation with the Neumann boundary condition. The new method is a completely meshless method without any differential processing. It can be seen that this method can be used to solve the low precision problem for the traditional collocation method when simulating the Neumann boundary condition problem from the one-dimensional and two-dimensional nonlinear examples. In addition, the algorithm has better error accuracy in dealing with Neumann boundary condition compared with the other algorithms. The numerical solution of more nodes in the whole region was simulated by using the optimal collocation points accurately, which can greatly reduce the amount of computation and save time. The validity and accuracy of the proposed method were verified by the one-dimensional and two-dimensional real background engineering problems.



---

*Source: 8298915-2018-04-08.xml*