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.

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:

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 [920] 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 [1020] 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 [2932].

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 [3337] is also called a distance basis function, which is a function of the distance from the node to the node . 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 :where is the Euclidean distance from node to node . Typical RBFs includemultiquadrics function: ,inverse multiquadrics function: ,Gaussians function: ,TPS (thin-plate splines): ,

where is a constant () 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 in the local support region near the point is approximated by the RBFs and its derivative at the Neumann boundary which can be expressed as [22]wherewhere 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. is the interior points, 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. and denoted the direction cosine in and directions on the Neumann boundary, respectively.

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.

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: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 aswhere is the radical basis function centered at source point 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.

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:where is the soil water content, 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), is the moisture content of the surface due to wet condition remain unchanged, 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 can be described aswhere and represents the distance between and . 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 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 and .

Step 4 (discrete the initial and boundary condition). Initial ConditionBoundary Condition
(1) Dirichlet Boundarywhere represents the distance between and .
(2) Neumann Boundarywhere represents the distance between and . 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: wherewherewhere and .

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:

Assuming , are nonlinear term and moment, let , ; use , , as the forecast value of , , ; get K(), D(θin), K(); make it as the forecast value of K(), D(), K(); the solution is solved by the first iteration; then get K(), D(), K() separately, and make it as the correction value of K(), D(), K() and as the next forecast value of . 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 iswhere is the number of iterations and is the allowable relative error, which is determined by the exact requirement of the calculation; for example, . If it satisfied (23), we can replace with .

4. Results and Discussion

Example 1 (one-dimensional unsaturated soil water movement equation). where , , exact solution , the root water uptake item , and the initial value and boundary value , are determined by exact solution. Gaussian function was selected as the radial function. In this paper, we discuss the error of -norm: The comparative result of calculation accuracy about different algorithm is shown in Figures 2, 3 and Tables 1, 2.

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 × 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).

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

Brooks-Corey Model

Van Genuchten Modelwhere in B-C model, in V-G model, is an index parameter representing the soil’s pore size distribution, in B-C model, in V-G model, and 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.

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 -norm error, we can get the -norm error about different models in different wetting front. The -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 -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 , , exact solution , root water absorption item , initial value , and boundary value , are determined by exact solution. The radial basis function chooses the Gaussian function , and the error was estimated using the -norm. Figure 4 gives the comparison result of numerical solution and exact solution in time about different algorithms. Table 3 gives the error accuracy of this paper in and 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.

Example 4 (two-dimensional unsaturated soil water movement equation in regular regional—point source infiltration problem).

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: ,  ,  ,  , and  .

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

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 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). where , , exact solution , root water absorption item , initial value , and boundary value , are determined by exact solution. The radial basis function chooses the Gaussian function , , and the error is estimated using the -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.

Similarly, as shown in Table 4, 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 Table 6 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 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.

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.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

Acknowledgments

This work was cosupported by National Natural Science Foundation of China (NSFC) (51409212, 51409213, 41371239).