Hermite Radial Basis Collocation Method for Unsaturated Soil Water Movement Equation

Due to the nonlinear diffusion term, it is hard to use the collocationmethod 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.


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: 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][10][11][12][13][14][15][16][17][18][19][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 highdimensional 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][11][12][13][14][15][16][17][18][19][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][30][31][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.

Hermite Radial Basis Collocation
Method (HRBCM) We can obtain algebraic equations by interpolating at the kth point in the local support domain: where We also can be obtain algebraic equation by interpolating at the mth point on the Neumann boundaries: As shown the Figure 1, "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.

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: where  is the differential operator defining the governing differential equations to be satisfied in the open domain Ω,  is the differential operator defining the boundary condition on boundary Γ,  is the problem unknown,  is the source term, and  is the term associated with boundary conditions.For the unknown , the approximation represented by HRBFs can be expressed as where   (x) is the radical basis function centered at source point x  and   is the expansion coefficient.  is the number of source points in the domain, and   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: when   +   >  and ( 9) is the overdetermined system to be solved using the least squares method.When   +  = , (9) can be summarized as where  is the coefficients matrix formed by the operators  and  under the action,  is the vector formed by the unknown quantity, and  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.

Hermite Radial Basis Collocation Scheme
where  is the soil water content,  0 is the initial water content, () and () are the unsaturated soil water diffusivity and conductivity, (, , ) is the root water absorption term,  is the time,  is the horizontal distance,  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: where () = ()/⋅/ and () = ()/⋅/− ()/.
Step 2 (constructing the approximation function).By taking the Gaussian function as an example, the approximation function at the point x(, ) can be described as where and   represents the distance between x and x i .Thus, the approximate solution of the Richards' equation can be expressed by Hermite radial basis function interpolation as follows: where   is the number of boundary points,   is the number of inner points, and x  = (  ,   ) at the nodes numbered .
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: Collating (15), where ( +1 Step 4 (discrete the initial and boundary condition).

Initial Condition
Boundary Condition (1) Dirichlet Boundary where   represents the distance between x and x i . ( where   represents the distance between x  and x  .   and    are the cosine of the Neumann boundary in the  and  directions. Step 5 (constructing nonlinear equations).Combining ( 16), ( 17), (18), and ( 19), the matrix equation can be written as follows: where where where ( +1

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: Assuming where  is the number of iterations and  is the allowable relative error, which is determined by the exact requirement of the calculation; for example,  = 10

Results and Discussion
Example 1 (one-dimensional unsaturated soil water movement equation).
It can be seen from Figures 2, 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  will increase with increasing the space nodes when we choose the different time steps, space distributions, and parameters .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 Two most common infiltration models used to calculate the unsaturated conductivity () and diffusivity () 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).

Van Genuchten Model
where 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 onedimensional vertical (Figures 4(b) and 4(d)) and horizontal (Figures 4(a) and 4(c)) infiltration for the two soils at different times.
Figures 4(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  2norm error, we can get the  2 -norm error about different models in different wetting front.The  2 -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  2 -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 Figures 2, 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).
where () = , () = 0.01 + 0.01, exact solution  =   (1 + )(1 + ), root water absorption item (, , ), initial value  0 (), and boundary value  1 ,  2 are determined by exact solution.The radial basis function chooses the Gaussian function () =  − 2 ( > 0), and the error was estimated using the  2 -norm.Figure 4 gives the comparison result of numerical solution and exact solution in  = 1 time about different algorithms.Table 3 gives the error accuracy of this paper in  = 0.5 and  = 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 Figure 5 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  increases.
Table 4 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 .
According to the data in Table 4, we can fit the relationship between parameter  and the number of nodes  and give the empirical formula as follows: It can be seen from the above equation that, with the increase of the spatial nodes , the parameter  increases.
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:   = 0.041,   = 0.412,  = 0.322,   = 0.04317, and ℎ  = 14.66.
In 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 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 × 10 3 and 1.1646 × 10 3 when the vertical wetting front was 48 cm.The error is very small, which is 1.828%.
From Figures 5 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).time for simulating the soil moisture equation with Neumann boundary condition when we choose optimal collocations points, the appropriate time step, and parameter C.
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 Table 6.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  of the basis function remains unchanged and the spatial node increased.The error accuracy is the same at the same moment.

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 onedimensional 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.

Figure 2 :
Figure 2: The comparison between numerical solution and exact solution about different algorithm in  = 0.5 time.

Table 1 : 9 Table 2 :FDMFigure 3 :
Figure 3: The error results between numerical solution and exact solution about different algorithm in  = 0.5 time.

Figure 4 :
Figure 4: Time step is 1, comparison of numerical solution and simulation solution of HYDRUS software about the sandy loam and silt at different times.

Figure 5 :
Figure 5: Time step is 0.1, comparative results of numerical solution and exact solution in  = 1 time.
is the number of nodes in the local support domain and   is the number of nodes on the Neumann boundary; let   +  = ;  represents the total number of nodes.x  is the interior points, x  is the boundary point,   is the RBFs,   is the undetermined coefficient of the corresponding RBFs, and   is the number of undetermined coefficients corresponding to the derivatives of the RBFs on the Neumann boundary.   = cos(, ) and    = cos(, ) denoted the direction cosine in  and  directions on the Neumann boundary, respectively.

Table 4 :
Time step is 0.01, error accuracy of this paper in  = 0.5 and  = 1 time with the same collocation points number and the calculation points number.

Table 5 :
Comparison of the results of this paper and other algorithms using the 16 collocation points in Table4to construct the approximate function to calculate more node values.

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.