A finite point algorithm for soil water-salt movement equation

*Correspondence: tabdeljawad@psu.edu.sa; muhammad.abbas@uos.edu.pk 3Department of Mathematics and General Sciences, Prince Sultan University, Riyadh 11586, Saudi Arabia 6Department of Mathematics, University of Sargodha, Sargodha 40100, Pakistan Full list of author information is available at the end of the article Abstract In this paper, we propose the meshless finite point method for solving a type of fluid flow problem. The moving least square function is combined with the collocation method to treat nonlinear oneand two-dimensional soil water-salt movement equations. An adaptive windward scheme is used to stabilize the numerical solution in regions with a large gradient change. Numerical examples with the comparison among the proposed method, finite element method and characteristic finite element method show that the meshless finite point method is more accurate and is used to eliminate the numerical oscillation phenomenon.


Introduction
Soil salinization caused by unreasonable irrigation has become a major bottleneck restricting the sustainable development of agriculture [1]. With the increase of world population, world food production and water shortages are under greater pressure. Therefore, the rational development of inferior water resources and the improvement of the utilization of saline-alkali land and soil water have become an issue of focus which all countries are concerned about. Studying the process of soil water-salt balance migration and carrying out quantitative numerical simulations of an appropriate mathematical model provide the necessary foundation for monitoring, evaluating and controlling soil salinization in arid areas. It is particularly important to prevent the accumulation of salinity, the loss of water and soil resources, and the desertification of land. And it is also necessary to save water, make irrigation, and increase production and income. The study of soil water-salt movement [2] is derived from the study of soil water movement, which originates from Darcy's law. Since soil water movement and salt transport occur simultaneously and interact with each other, soil salt migration has arisen. As early as the 1960s, the research on the movement of water and salt has gradually developed from simple matter balance calculation to mathematics and computer simulations, which makes the concept of water and salt movement clearer and more quantitative. The numerical simulation can set the initial and boundary conditions flexibly, meanwhile it can describe the process of water and salt mi-gration, which is more consistent with the actual situation. Scholars such as Shi Yuanchun, Li Yunzhu and Jia Dalin have done a great deal of research on the simulation in soil salinization treatment. Zhang Weizhen [3] proposed a preliminary study on the simulation of soil water and salt transportation, therefore, the study of soil water and salt transport has entered a complete new stage. Most of the studies are aimed at simulating the infiltration characteristics of water infiltration and the transport of water [1], however, the salt migration of saline-alkali soil is less simulated. The purpose of this paper is to study the law of water and salt migration under the conditions of water infiltration, the accuracy and validity of this algorithm are verified by practical problems, which provides a scientific basis for reasonable irrigation.
The convection diffusion model is classical and popular among researchers of soil watersalt transport models. However, only a handful of partial differential equations can be obtained with exact solutions due to the heterogeneity of soil and the complexity of the conditions. Therefore, the numerical solution of partial differential equations continues to evolve. At present, the numerical methods for solving the equations of soil water and salt movement are the finite difference method (FDM) [4,5], the finite element method (FEM) [6][7][8] and the characteristic finite element method (CFEM) [9,10]. The central idea of FEM and CFEM is that a problem domain can be split into small, non-overlapping elements so that simple interpolation functions can be used to approximate any field function within each element. In addition, the Newton-Kantorovich method can also better solve the problems associated with partial differential equations [11][12][13][14]. However, shape functions for heavily distorted elements will not produce acceptable numerical solutions. One way of dealing with such distortions is to re-mesh the local domain and apply adaptive techniques, which has a high computational cost and adaptive analysis is difficult.
A meshless method based on nodes can completely or partially eliminate the grid. Therefore it overcomes the dependence on the grid and is superior to FEM and CFEM, because it not only reduces the difficulty of the simulation calculation, but also improves the accuracy of the calculation. There are many different meshless methods [15,16]. The main difference lies in the formation of approximate function and discrete control equation. The shape function is constructed using the moving least squares method (MLS), with its control equation discretized using the collocation method [17]. Without requiring a background grid and numerical integrals, the finite point method has many advantages, including a simple discrete format, short calculation time (CPU) and easy programming. It has been successfully applied to compressible flow [18], sound propagation and other practical problems [19][20][21][22].
In this paper, a numerical solution for the convection dispersion model is calculated using the finite point method. A stabilizing term is added to the control equation before discretizing governing equations. Along with an increase in the Reynolds number, Re, the adaptive upwind support domain is used. This method is called the upwind finite point method (UFPM). This numerical method can eliminate the uncertainty in the boundary area and the case of large gradient changes.
The structure of this paper as follows: the moving least squares approximation, finite point method and its numerical algorithm are given in Sect. 2, the accuracy of the solution is examined with some examples and the actual problem in Sect. 3, and finally, conclusions are summarized in Sect. 4.

Construction of the finite point algorithm 2.1 Moving least squares approximation
Consider an approximation function u(x) in closed domain ∈ R d (d = 1, 2 or 3). Divide into subdomain i (i = 1, 2, . . . , n), and use i to represent the coverage of , so as to obtain the local approximation of u(x). A point x i and a set of points x j (j = 1, 2, . . . , m) surrounding x i form a local point, thus completing i . Assuming that u(x) is sufficiently smooth in i , it has the following form: In Eq. (2.1), the weighted least squares are used to approximate the local approximation to obtain the coefficients a j (x). Define the following form: The coefficients a(x) can be found through the extreme value of J, and the value is as follows:

4)
Thus, u h (x) is As the weight function is compact, the global approximation function of the unknown function can be obtained: where the shape function is given by (2.9) Let γ T = A -1 p, then Eq. (2.9) can be expressed as = γ T B. The derivatives of the shape functions are, respectively,

Collocation method
The governing equation with boundary conditions is assumed to be of the following form: where L is a differential operator defining the governing equations in the domain , and B is the differential operator defining the boundary condition at boundary . The collocation method uses N points selected in the domain and N points selected in the boundary to satisfy Eqs. (2.12) and (2.13), respectively. The unknown function u(x) may be approximated at points x i (i = 1, 2, . . . , N ) with a compact support function ϕ j (x) as follows: (2.14) This is the main idea of the collocation method, put Eq. (2.14) into Eqs. (2.12) and (2.13).
When N + N = N , Eq. (2.14) is the determined system needed to approximate the value of each node x i (i = 1, 2, . . . , N ). The method studied in this paper is the problem N + N = N . If L, B are nonlinear operators, we can obtain the nonlinear equation group. An iterative method is needed to solve the node approximation.
There are two forms of the support field given in Figs. 1-2. One is completely skewed to the upwind side of the support domain; the other is the adaptive upwind support domain. d u is the distance between the center of the old and new support areas, i.e., the adaptive upwind support region allows the collocation point to shift in the opposite direction of  A schematic diagram of the two support domains is shown below.

Finite point algorithm for soil water-salt transport equation
Water is the carrier of soluble salt. In addition, it also reflects the macro regulation in a saline-alkali area. Corresponding models have been established for the water movement and salt transport equations. The dynamic change in soil water content has been simulated using the finite point method. In this paper, the algorithm developed is based on the stable finite point method, in which the large gradient change making its support domain gradually tends to be windward. That is to say that we are to use the adaptive windward format to solve the stability of numerical solution. The algorithm is given now: Step 1. The nodes are divided in the solution domain according to certain rules.
Step 2. A stability term is applied to the control equation.
Step 3. The moving least squares method is used to construct the shape function and the approximate expression of the unknown function is obtained. Step 5. Generate and solve a system of algebraic equations.

Finite point algorithm for water movement equation
The Richards equation with water content as the variable is Step 1: Construct the approximation function Step 2: Discretize the control equation.
The governing equations in (2.15) have the following form: where u = (-∂K(θ) ∂θ , -∂K(θ) ∂θ ) is the coefficient of convection term and h is the characteristic length. Assuming K(θ ) = θ , D(θ ) is a linear function of θ . Putting Eqs. (2.17) and (2.18) into (2.19), To deal with the problem of two-dimensional partial differential equations, a bottomup principle is used for rearranging the calculation points (x, z) in two dimensions. The space variables in the equation are semi-discretized. The time variable is fully discrete. The value of the shape function in each node is solved in the support domain and put into the approximation expression for the unknown function. The following algebraic equations can be obtained: There are two ways to solve the equations. One method is the forward iteration method when the right side of the equation is all about the value of the m layer of θ . Another method is the predictor corrector method when the equation gets this form: θ m+1 = A m+1 θ m + F m+1 .

Finite point algorithm for the soil salt transport equation
The convection dispersion equation for non-volatile, non-reactive, non-adsorptive soil water and salt transport can be written as where the hydrodynamic dispersion coefficient D 1h is formed by the superposition of the molecular diffusion coefficient and mechanical dispersion coefficient, and J w is the soil moisture flux. In the case of saturated flow, the water content θ is constant, and Eq. (2.22) gets the following form: For a steady flow condition, Eq. (2.23) gets the following form: where D v = D L /θ is the convection dispersion coefficient, and v = J w /θ is the mean pore water velocity. The finite point algorithm is shown for the migration equation of nonvolatile, non-reactive, non-adsorptive soil water and salt: where c is the salt concentration, c 0 (x) is the initial value, ψ 1 (t), ψ 2 (t) are boundary values.
Step 1: Construct the approximation function where ϕ ip (x) is the shape function of the point x i , obtained for x p in the support domain, ϕ ip (x p ) is the value of ϕ ip (x) at point x p , and n is the number of nodes in the local support domain for the shape function at point x p .
Step 2: Discretize the control equation A stability term is applied to governing Eq. (2.25) as follows: The space variables in the equation are semi-discrete and the interval I = [0, L] has the number of N nodes with the x = x i (i = 1, 2, . . . , N ) distribution. Discretization at the node x p is as follows: We can remove the term ∂c ∂t | x=x p at the time layer m t in Eq. (2.29) by using the forward difference method, (2.32) Equation (2.31) can be written as

Numerical simulation 3.1 Linear soil salt transport equation
Example 1 The one-dimensional steady-state salt transport is  Table 1 shows that the finite point computation error with the stability term only, the finite point error of the adaptive upwind scheme and the error of finite element method are given, respectively, in the case of convection dominance.
As can be seen from Table 1, the proposed method is very good at eliminating error when convection prevails. Numerical simulation results of the FPM are better than that   where e 1 , e 2 are the maximum errors between the exact solution and the numerical solution at the nodes, and p indicates the convergence order. As can be seen from Table 2, the higher number of nodes, the higher the calculation precision. As the convergence order is positive, the proposed algorithm in this paper is convergent.
According to Fig. 3 and Fig. 4, the gradient change is relatively large in the vicinity of x = 1. This is because the diffusion term is ignored. A thin boundary layer is formed at the downstream boundary where x = 1, which is difficult to replicate using a standard numerical algorithm. So it can lead to the shocked solution in the boundary layer. To solve this kind of problem, several methods can be adopted including encrypting the node, expanding the local support domain, and using an adaptive upwind scheme. The method adopted in this paper is the adaptive upwind scheme. As the Reynolds number (Re) increases, the support field gradually inclines to the upwind side, so that the solution can capture more upstream information. It can get the numerical solution close to the exact solution, which shows the effectiveness of this method.   Example 2 The one-dimensional non-steady-state salt transport is The exact solution is c(x, t) = (e -x + x -x e -1)e (D v +v)t -(x -x e -1)e (D v +v)t . Let v = 0.1, D v = 0.001, and use a Gauss function as the weight function. According to Table 3, Table 4 and Table 5, the comparison of the proposed method with FEM shows that UFPM is highly accurate. As can be seen from Table 6, the higher the number of nodes, the higher the calculation precision when time step and time are determined. As the convergence order is positive, UFPM in this paper is convergent.
According to Fig. 5 and Fig. 6, by changing the time step, the number of nodes and using the Gauss function, numerical solutions of the proposed algorithm fit well with the exact solution.
According to Table 7, the error between the numerical and exact solution is very small when the time is calculated to be 6.0. As can be seen from Table 8, the comparison of the proposed algorithm with CFEM shows that the proposed algorithm is highly accurate.
According to Fig. 7 and Fig. 8, the results of CFEM are not good when the number of nodes is 51. CFEM needs to increase the number of nodes to improve accuracy. However, the proposed algorithm can obtain good simulation results whether the number of nodes is 51 or 71.
where D(θ ) = 0.001θ + 0.001, K(θ ) = θ , the exact solution is θ = (t + 1)(xx 2 )(zz 2 ). Initial condition: θ 0 = (xx 2 )(yy 2 ), the resource item f , the initial value θ 0 (x, z) and the boundary value ψ 1 (z, t), ψ 2 (z, t), ψ 3 (x, t), ψ 4 (x, t) are determined by exact solution. Table 9 gives the relative error and calculation time for the proposed method when the time step is 0.  According to Table 9, the size of scale has an impact on the error. That is because the size of support domain directly affects the number of nodes in the support domain. On the whole, the error between the numerical and exact solutions is very small. From the perspective of computing time, the computational efficiency of the proposed algorithm is quite good. As can be seen from Table 10, the proposed algorithm is very stable. But FEM becomes unstable with the increase of time.
Example 5 Simulation and analysis of soil water movement equation under the background of engineering water conservancy.
The one-dimensional soil water movement equation has the following form: Here K(θ ) is the unsaturated hydraulic conductivity, h is the soil-water suction, and x is the horizontal position. θ is the soil water content, θ s is the saturated soil water content, and θ 0 is the initial soil water content. D(θ ) is the diffusion rate. Equation (3.5) repre- Figure 11 The comparison of this paper solution and simulation solution of Hydrus software at t = 500 min and t = 1 sents the vertical infiltration when α = 0, and the horizontal infiltration when α = π/2. The paper takes the Brooks-Corey [19] model as an example in engineering to discuss the feasibility of this method. So the expressions of D(θ ) and K(θ ) are as follows: Here n is an index parameter, m = 2 + 3n, h d is the bubbling pressure (cm). The physical parameters of a given soil are as follows: θ r = 0.015, θ s = 0.486, θ 0 = 0.2, K s = 0.0113333, n = 0.211, h d = 20.747. Figure 11 shows the simulated results of the one-dimensional vertical (right) and horizontal (left) infiltration for the given soil. The solid red line is the solution of this method, the dotted line represents the solution of the Hydrus software [13,14].
To be compared with other algorithms, the errors of vertical and horizontal infiltration are 1.141% and 0.5% when the infiltration time is 500 min, respectively. As can be seen from Table 9 and Fig. 11, to a certain degree, this method is not only suitable for numerical examples, but also feasible for the actual background problem.

Conclusions
Through fully discussing the effects of support domain size, node number and time step size on the simulation results, it can be concluded that the number of nodes has a direct impact on the calculation accuracy. Under normal circumstances, the higher the number of nodes, the higher the accuracy of the numerical results. In addition, the smaller the time step, the smaller the calculation error. It can be seen from the convergence order that the proposed method is feasible. The proposed method applies stability terms to the governing equations before the discretization of the governing equations and makes its support domain gradually tend to windward in the area of large gradient change so that it can effectively avoid numerical shock and obtain good simulation results when dealing with nonlinear water movement equation. Numerical examples with comparison between the FEM and CFEM indicate that the proposed method has high precision. The numerical examples and the problems of hydraulic engineering in the background indicate that this method can be used to solve the problem of alkaline areas. The accuracy and validity of this method have been verified, which provides a reasonable scientific basis for reasonable irrigation saline-alkali land. In addition, it has a certain positive effect on the development of agriculture. We will focus on studying the finite point procedure for solving three-dimensional problems in future work.