Numerical Simulation of Soil Water Movement under Subsurface Irrigation

Abstract

By constructing a radial basis function collocation method combined with a difference method, a two-dimensional mathematical model with boundary conditions of soil water movement under irrigation is proposed. The nonlinear term is dealt with a difference method and the equation is solved using an implicit scheme. In addition, the existence and uniqueness of the solution to the soil water movement equation are proven. Numerical results show that the proposed method has high precision and is easier to use than traditional methods. Moreover, the selection of parameter c plays an important role in guaranteeing calculation precision. It lays the foundation for the numerical solutions to high-dimensional soil water movement equations.

1. Introduction

Between 1960 and 2000, numerical techniques such as boundary element and finite element methods have made spectacular advances in the computation of physical phenomena in engineering and sciences. Particularly, the last two decades have witnessed substantial efforts in developing a new class of numerical methods. These methods, using traditional methods that require certain underlying meshes such as the triangulation of a region for computation, were primarily developed for simulating low-dimensional application problems. However, designing relevant meshes is usually quite a difficult task for two-dimensional regions and might become impossible for higher-dimensional problems.

However, the triangulation process is too time consuming, even if a sophisticated mesh-generator is employed. In fact, for a given distribution of points, it is possible to achieve a mesh quickly, but this always requires a considerable number of iterations including manual interaction before reaching a satisfactory mesh.

Meshless methods were first introduced by Lucyt, and Gingold and Monaghan in 1977 [1, 2]. Further works, such as Nayroles et al. [3], Belytschko et al. [4], Schaback [5], Sukumar et al. [6], Wu [7], and Wendland [8, 9] (see books [10, 11]), attempted to reduce or even eliminate the need for discretization of a domain or surface in the context of numerical solutions for boundary and initial value problems. These rapidly developed methods can not only reduce the large costs of labor, but also spare computational time compared to the finite element method, boundary element method, and other mesh-dependent methods.

The initial idea of meshless methods could date back to the smooth particle hydrodynamics (SPH) method for modeling astrophysical phenomena [1, 2]. In 1992, Nayroles et al. proposed the diffuse element method (DEM) [3]. In 1994, Belytschko et al. improved the DEM and introduced the element-free Galerkin method (EFG) [4]. Since then, there has been a great deal of research into meshless methods. Meshless methods, such as the reproducing kernel particle method (RKPM) by Liu et al. [12], Hp-cloud method by Duarte and Oden [13], natural element method (NEM) by Sukumar et al. [6], partition of unity finite element method (PUFEM) by Babuška and Melenk [14], and meshless Galerkin method using radial basis functions (RBF) by Wendland [9], have also been described in literatures. The major differences between these meshless methods, all of which can be classified as Galerkin methods, come only from the techniques used for interpolating the trial functions. Even though no mesh is required in these methods for the interpolation of the trial or test functions and the solution variables, the use of shadow elements is unavoidable during the integration of the symmetric weak form. Therefore, these methods are not truly meshless.

The finite point method (FPM) introduced by Oñate et al. [15], which is based on the moving least-squares and the collocation method, is a truly meshless method. Other truly meshless methods include the PCM [16], the Hp-meshless clouds [17], the LBIE [18], the LSC [19], and the WLSM [19] among others. In particular, in 1990, Kansa introduced the radial basis function (RBF) collocation method for solving elliptic, hyperbolic, and parabolic partial differential equations (PDEs) [20]. Later, the RBF was further developed by Schaback [5] and Fasshauer [21]. Hon et al. applied the RBF to the numerical computation of the variable normal equations, nonlinear Burgers equation, and shallow water wave equation [22]. In 1999, Rippa carried out work on selecting the correct shape parameter 𝑐 for the RBFs [23]. Another class of RBF was further developed by Wendland [8, 9], Wu [7], and Buhmann [24]. Zhang et al. have made more research on element-free kp-Ritz method and applied it to solve different problems [25–29]. Therefore, the collocation-based meshless method has been an important meshless method in the current literature [30].

By creating a univariate basis function with an Euclidean norm, meshless methods are often naturally radically symmetric, and the high-dimensional problem can be turned into virtually one-dimensional one. Consequently, the study of the numerical solutions of PDEs through radial basis function interpolation has yielded a number of significant results.

1.1. The Problem and the Approach to Solve the Problem

The FEM has difficulty in remeshing and adaptive analysis. In contrast, meshless methods do not require a mesh to discretize the domain, and the approximate solution is constructed entirely with a set of scattered nodes [31]. However, meshless methods may lead to lower computational efficiency than FEM because more computational effort for the meshless interpolation and numerical integrations are required [32]. Hence, the improvement of the computational efficiency of meshless methods targeted at meshless interpolation and numerical integrations becomes an important issue. Other concerns or weaknesses of the existing meshless methods include difficulty in introducing the essential boundary conditions, a greater cost in evaluating the shape function derivatives, problems in handling discontinuities such as those due to heterogeneous material distributions, and the need for complicated node connectivity to ensure accurate results.

A truly meshless method, based on collocation with radial basis functions and radial basis functions are chosen to represent the solutions of PDEs, is the main focus of this paper. Moreover, the collocation-based meshless method is a truly meshless technique without mesh discretization. This method has the advantage of higher accuracy, convenient for computing, and has been successfully applied to numerical solutions of various PDEs. By using collocation with radial basis functions, the partition of the domain is not needed; hence the method can be applied to complex domains and overcomes some drawbacks of the traditional finite element method. It only needs to compute the shape functions and their derivatives, while the finite element method requires a calculation of the relevant integrals that often lowers the efficiency of computation significantly. The boundary conditions are relatively easier to be imposed without special manipulation. Finally, the way of obtaining the solutions to the PDEs is intuitive. From a practical point of view, this approach can achieve higher accuracy with more easily coded computer programs. Therefore, the collocation method with radial basis function interpolation can produce good results.

1.2. Mathematical Model of Soil Water Movement

Assume that permeable pipes in fields are parallel and have the same spacing and depth. Figure 1 shows a soil profile which is perpendicular to the permeable pipeline. Assuming that the water flow in the soil is negligible, then the problem considering irrigation can be solved as a plane problem. Obviously, in the case of Figure 1, it is sufficient to focus on soil water movement in the shaded part abcd, namely, the flow region.

Figure 1 

Schematic diagram of flow region.

In the region of soil water flow, both perpendiculars ad and bc are symmetrical lines, where ad traverses through the center of a permeable pipe and bc is equidistant from two adjacent permeable pipes. The symmetry means that the water flux along the normal direction of two lines equals zero. Line ab, the earth's surface, is the upper boundary of the region of soil water flow. The water flow which crosses line ab is affected by the soil’s conveyance capacity, which is in turn related to the surface soil transpiration and strength of rainfall. By ignoring the transpiration and rainfall effects, the water flux across ab is zero. Line dc is the lower boundary of the region of soil water flow but is deep enough so that water movement caused by irrigation cannot reach it.

Assume that, in the region of soil water flow, the soil water flow is continuous in time and space, obeying the law of mass conservation. Then we can develop mathematical models of the soil water movement based on the physical process through the area.

The mathematical model of 2D water movement in unsaturated soil can be established as 𝜕𝜃=𝜕𝜕𝑡𝐷(𝜃)𝜕𝑥𝜕𝜃+𝜕𝜕𝑥𝐷(𝜃)𝜕𝑧𝜕𝜃+(𝜃)𝜕𝑧𝜕𝐾,𝜕𝑧 (1) where 𝜃 is the volumetric soil water content; 𝐷(𝜃) and 𝐾(𝜃) both are continuous functions which denote the diffusivity and hydraulic conductivity of the unsaturated soil flow respectively; 𝑡 is time and 𝑧 is the distance with the upward direction defined as positive.

The remainder of this paper is organized as follows. In Section 2, we introduce the construction of the meshless method based on the collocation method with radial basis functions, where the explicit schemes for the given equation are proposed and the existence and uniqueness of the solution to the method are proved. In Section 3, we give some applications and apply the method developed in this paper to examine the appropriateness and efficiency of the numerical solutions and to analyze the factors influencing the accuracy. We also compare our method to FEM. Section 4 comprises the conclusions and remarks.

2. Construction of a Meshless Method Based on Collocation with Radial Basis Functions

2.1. Radial Basis Functions

Radial basis functions (RBFs), also known as distance basis functions, are a type of functions with a basic variable 𝑑_𝑗  (𝑑_𝑗=‖𝑥−𝑥_𝑗‖). They are isotropic and simple in form and can be solved easily with numerical calculation. The method, which combines radial basis functions with collocation, has many advantages when solving the partial differential equation [33, 34]. It has no correlation with the space dimensions and does not need any element or mesh for interpolation. Therefore, it is a truly meshless method. Currently, FEM [35] and FDM [36–38] are mostly used to solve soil moisture movement equations. With the development of meshless methods in recent years, many scholars now solve the soil water seepage problem using the finite volume method [39] and the RBF method [40]. The mathematical models of soil water movement are usually convection-diffusion equations, meaning that numerical oscillation occurs frequently when FDM or FEM is used. However, the radial basis function collocation method can solve many of these problems.

2.2. Construction of Radial Basis Function Collocation Method

The 2D soil water movement equation is 𝜕𝜃=𝜕𝜕𝑡𝐷(𝜃)𝜕𝑥𝜕𝜃+𝜕𝜕𝑥𝐷(𝜃)𝜕𝑧𝜕𝜃+(𝜃)𝜕𝑧𝜕𝐾,()𝜃()𝜕𝑧𝑥,𝑧∈Ω,𝑥,𝑧,0=𝜃₀,()𝜃()𝑥,𝑧∈Ω,𝑥,𝑧,𝑡=𝜃₁,(𝑥,𝑧)∈𝜕Ω, (2) where 𝜃 is the volumetric soil water content, 𝐷(𝜃) and 𝐾(𝜃) are the diffusivity and hydraulic conductivity of the unsaturated soil flow, respectively, 𝑡 is time, and 𝑧 is the distance with the upward direction defined as positive. 𝜃₀ is the initial water content, 𝜃₁ is the constant water content of ground under humid conditions, Ω is the seepage region, and 𝜕Ω is the boundary of Ω.

Because the quadratic term in (2) is nonlinear, we cannot use the collocation method directly. Thus, we apply centered differences to deal with the nonlinear term in (𝑥_𝑖,𝑧_𝑗): 𝜕𝐷(𝜃)𝜕𝑥𝜕𝜃𝜕𝑥_{𝑥=𝑥𝑖}≈1𝐷(𝜃)2Δ𝑥𝜕𝜃𝜕𝑥_{𝑥=𝑥𝑖+1}−𝐷(𝜃)𝜕𝜃𝜕𝑥_{𝑥=𝑥𝑖−1}=1𝐷𝜃2Δ𝑥^𝑖+1,𝑗𝜕𝜃_𝑖+1,𝑗𝜕𝑥_{𝑥=𝑥𝑖+1}−𝐷𝜃_{𝑖−1,𝑗}𝜕𝜃_{𝑖−1,𝑗}𝜕𝑥_{𝑥=𝑥𝑖−1},𝜕𝐷(𝜃)𝜕𝑧𝜕𝜃𝜕𝑧_{𝑧=𝑧𝑗}≈1𝐷(𝜃)2Δ𝑧𝜕𝜃𝜕𝑧_{𝑧=𝑧𝑗+1}−𝐷(𝜃)𝜕𝜃𝜕𝑧_{𝑧=𝑧𝑗−1}=1𝐷𝜃2Δ𝑧^𝑖,𝑗+1𝜕𝜃_𝑖,𝑗+1𝜕𝑧_{𝑧=𝑧𝑗+1}−𝐷𝜃_{𝑖,𝑗−1}𝜕𝜃_{𝑖,𝑗−1}𝜕𝑧_{𝑧=𝑧𝑗−1}, (3) where Δ𝑥 and Δ𝑧 are spatial intervals, respectively. (𝑥_𝑖,𝑧_𝑗) are boundary points when 𝑖,𝑗=1,𝑁, are inner points when 𝑖,𝑗=2,3,…,𝑁−1.

Discretizing the left side of (2) with forward differences gives 𝜕𝜃𝜕𝑡^{𝑡=𝑡𝑛+1}_{(𝑥𝑖,𝑧𝑗)}≈𝜃^𝑛+1−𝜃^𝑛Δ𝑡_{(𝑥𝑖,𝑧𝑗)}=𝜃^𝑛+1_𝑖,𝑗−𝜃^𝑛_𝑖,𝑗.Δ𝑡 (4)

Discretizing the third term in the right side of (2) with central differences gives 𝜕𝐾(𝜃)𝜕𝑧^{𝑡=𝑡𝑛+1}_{(𝑥𝑖,𝑧𝑗)}≈𝐾𝜃^𝑛+1_𝑖,𝑗+1−𝐾𝜃^𝑛_{𝑖,𝑗−1}.2Δ𝑧 (5)

Let function 𝜃(𝑋,𝑡^𝑛) be an approximation of 𝜃(𝑋,𝑡^𝑛): 𝜃𝑋,𝑡^𝑛=^𝑁𝐼_𝑖=1𝛼^𝑛_𝑖𝜑𝑋−𝐼_𝑖+^𝑁𝑏_𝑖=1𝛽^𝑛_𝑖𝜑𝑋−𝐵_𝑖, (6) where 𝑋=(𝑥,𝑧), 𝐼_𝑖∈Ω, and 𝐵∈𝜕Ω. 𝑁_𝐼 is the number of the nodes in the region and 𝑁_𝑏 is the number of the nodes on the border.

Applying collocation method, (6) should satisfy the differential equation for the region Ω in (7) and boundary conditions on the borders 𝜕Ω in (8). That is, ^𝑁𝐼_𝑖=1𝛼^𝑛_𝑖𝜑𝐼_𝑗−𝐼_𝑖−𝑑𝐷𝜃𝑑𝑥^𝑛𝐼_𝑗𝑑𝜑𝐼_𝑗−𝐼_𝑖−𝑑𝑑𝑥𝐷𝜃𝑑𝑧^𝑛𝐼_𝑗𝑑𝜑𝐼_𝑗−𝐼_𝑖+𝑑𝑧_{𝑁𝑏𝑖=1}𝛽^𝑛_𝑖𝜑𝐼_𝑗−𝐵_𝑖−𝑑𝐷𝜃𝑑𝑥^𝑛𝐼_𝑗𝑑𝜑𝐼_𝑗−𝐵_𝑖−𝑑𝑑𝑥𝐷𝜃𝑑𝑧^𝑛𝐼_𝑗𝑑𝜑𝐼_𝑗−𝐵_𝑖𝑑𝑧=𝜃^𝑛_𝑖,𝑗+𝜕𝐾𝜃^𝑛𝐼_𝑗𝜕𝑧,𝑗=1,2,…,𝑁_𝐼,_{𝑁𝐼𝑖=1}𝛼^𝑛_𝑖𝜑𝐵_𝑗−𝐼_𝑖+^𝑁𝑏_𝑖=1𝛽^𝑛_𝑖𝜑𝐵_𝑗−𝐵_𝑖=𝜃^𝑛₁𝐵_𝑗.𝑗=1,2,…,𝑁_𝑏. (7)(8)

Let 𝜓−𝑑(‖𝑋−𝐼‖)=𝜑(‖𝑋−𝐼‖)𝐷(𝜃(𝑋𝑑𝑥))𝑑𝜑(‖𝑋−𝐼‖)−𝑑𝑑𝑥𝐷(𝜃(𝑋𝑑𝑧))𝑑𝜑(‖𝑋−𝐼‖)𝑓𝑑𝑧_𝑛𝐼_𝑗=𝜃^𝑛−1𝐼_𝑗+𝜕𝐾𝜃^𝑛𝐼_𝑗.𝜕𝑧 (9) Then (7) and (8) can be expressed with the following matrix equation: 𝐇𝐔=𝐅, (10) where𝜓(0)𝐼𝐇=⋯𝜓₁−𝐼_𝑁𝑖𝜓𝐼₁−𝐵₁⋯𝜓𝐼₁−𝐵_𝑁𝑏𝜓𝐼⋮⋮⋮⋮⋮⋮_𝑁𝑖−𝐼₁⋯𝜓(0)𝜓𝐼_𝑁𝑖−𝐵₁⋯𝜓𝐼_𝑁𝑖−𝐵_𝑁𝑏𝜑𝐵₁−𝐼₁⋯𝜑𝐵₁−𝐼_𝑁𝑖𝜑(0)⋯𝜑𝐵₁−𝐵_𝑁𝑏𝜑𝐵⋮⋮⋮⋮⋮⋮_𝑁𝑏−𝐼₁⋯𝜑𝐵_𝑁𝑏−𝐼_𝑁𝑖𝜑𝐵_𝑁𝑏−𝐵₁(0),𝛼⋯𝜑𝐔=^𝑛₁,𝛼^𝑛₂,…,𝛼^𝑛_𝑁𝑖,𝛽^𝑛₁,𝛽^𝑛₂,…,𝛽^𝑛_𝑁𝑏^𝑇,𝐅=𝑓^𝑛₁(𝐼₁),𝑓^𝑛₁(𝐼₂),…,𝑓^𝑛₁(𝐼_𝑁𝑖),𝜃^𝑛₁(𝐵₁),𝜃^𝑛₁𝐵₂,…,𝜃^𝑛₁𝐵_𝑁𝑏^𝑇. (11)

2.3. The Existence and Uniqueness of the Solution

Let 𝐀₁=𝜓(0)⋯𝜓𝐼₁−𝐼_𝑁𝑖𝜓𝐼⋮⋱⋮_𝑁𝑖−𝐼₁(0),𝐀⋯𝜓₂=𝜓𝐼₁−𝐵₁⋯𝜓𝐼₁−𝐵_𝑁𝑏𝜓𝐼⋮⋱⋮_𝑁𝑖−𝐵₁⋯𝜓𝐼_𝑁𝑖−𝐵_𝑁𝑏,𝐀₃=𝜑𝐵₁−𝐼₁⋯𝜑𝐵₁−𝐼_𝑁𝑖𝜑𝐵⋮⋱⋮_𝑁𝑏−𝐼₁⋯𝜑𝐵_𝑁𝑏−𝐼_𝑁𝑖,𝐀₄=𝜑(0)⋯𝜑𝐵₁−𝐵_𝑁𝑏𝜑𝐵⋮⋱⋮_𝑁𝑏−𝐵₁(0).⋯𝜑 (12) Then 𝐻=^𝐴1𝐴2_𝐴3𝐴4. For (10), the following theorem is obtained.

Theorem 1. If the Fourier transform 𝐹[𝜙] of 𝜙(𝜔) is almost everywhere larger than 0 and (𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄)⁻¹ exists, then 𝐻 is invertible. Consequently, the matrix equation (10) has a unique solution.

Proof. First, we can show that 𝐀₁⁻¹ and 𝐀₄⁻¹ exist. According to the characteristics of radial basis function, matrix 𝐀₁ is symmetric. Then we show that it is positive definiteness. For ∀𝜉(𝜉≠0)∈𝑅^𝑁𝑏, 𝐀₁𝜉,𝜉=^𝑁𝑏𝑁_𝑗=1^𝑏_𝑘=1𝜉_𝑗𝜉_𝑘𝜑_𝑗𝑥^′_𝑗=^𝑁𝑏𝑁_𝑗=1^𝑏_𝑘=1𝜉_𝑗𝜉_𝑘𝜑𝑥_𝑗−𝑥_𝑘=12𝜋^𝑁𝑏+∞_−∞𝐹𝜑𝜉(𝜔)²𝑑𝜔>0, (13) where 𝜉(𝜔)=∑^𝑁𝑏_𝑗=1𝜉_𝑗𝑒^{𝑖<𝜔,𝑥𝑗>}. So 𝐀₁ is a symmetric positive definite matrix; that is, 𝐀₁⁻¹ exists.
Because 𝜑(𝑟) is an even function, then −𝜑(𝑟) is an even function too. Moreover, 𝐹[−𝜑]=−(𝑖𝜔)²𝐹[𝜑]=𝜔²𝐹[𝜑]>0. Similarly, 𝐀₄ is also a symmetric positive definite matrix. Hence, 𝐀₄⁻¹ exits. Now, we have 𝐄−𝐀₂𝐀₄⁻¹𝐀𝟎𝐄¹𝐀₂𝐀₃𝐀₄=𝐀₁−𝐀₂𝐀₄⁻¹𝐀₃𝟎𝐀₃𝐀₄. (14) Computing the determinants of both sides of (14), we have |𝐇|=𝐀₁𝐀₂𝐀₃𝐀₄=𝐀₁−𝐀₂𝐀₄⁻¹𝐀₃⋅𝐀₄. (15) As (𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄)⁻¹ exists, 𝐀₁−𝐀₂𝐀₄⁻¹𝐀₃×𝐀₁⁻¹−𝐀₁⁻¹𝐀₂(𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄)⁻¹𝐀₃𝐀₁⁻¹=𝐄−𝐀₂(𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄)⁻¹𝐀₃𝐀₁⁻¹−𝐀₂𝐀₄⁻¹𝐀₃𝐀₁⁻¹+𝐀₂𝐀₄⁻¹𝐀₃𝐀₁⁻¹𝐀₂×𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄⁻¹𝐀₃𝐀₁⁻¹. (16) In addition, 𝐀₂𝐀₄⁻¹𝐀₃𝐀₁⁻¹𝐀₂(𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄)⁻¹𝐀₃𝐀₁⁻¹−𝐀₂(𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄)⁻¹𝐀₃𝐀₁⁻¹=𝐀₂𝐀₄⁻¹𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄×(𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄)⁻¹𝐀₃𝐀₁⁻¹=𝐀₂𝐀₄⁻¹𝐀₃𝐀₁⁻¹. (17) By (16) and (17), we obtain 𝐀₁−𝐀₂𝐀₄⁻¹𝐀₃×𝐀₁⁻¹−𝐀₁⁻¹𝐀₂𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄⁻¹𝐀₃𝐀₁⁻¹=𝐄. (18) So (𝐀₁−𝐀₂𝐀₄⁻¹𝐀₃)⁻¹=𝐀₁⁻¹−𝐀₁⁻¹𝐀₂𝐀₃𝐀₁⁻¹𝐀₂−𝐀₄⁻¹𝐀₃𝐀₁⁻¹. (19)
According to (15) and (19), we know that |𝐇|≠𝟎. Thus, the matrix equation (10) has a unique solution. Now the numerical solution to problem (1), based on the radial basis function collocation method, can be calculated.

3. Numerical Examples

Example 1. Solve the following linear model using the radial basis function collocation method: 𝜕𝜃=𝜕𝜕𝑡𝐷(𝜃)𝜕𝑥𝜕𝜃+𝜕𝜕𝑥𝐷(𝜃)𝜕𝑧𝜕𝜃()[]×[],𝜃()𝜕𝑧+𝑓,𝑥,𝑧∈Ω=0,10,1𝑥,𝑧,0=𝜃₀,()𝜃()𝑥,𝑧∈Ω,𝑥,𝑧,𝑡=𝜃₁,(𝑥,𝑧)∈𝜕Ω,𝑡>0. (20)

Let 𝐷(𝜃)=1; then the analytic solution of (20) is 𝜃(𝑧,𝑡)=(𝑥²−𝑥)(𝑧²−𝑧)𝑡 and 𝑓=(𝑥²−𝑥)(𝑧²−𝑧)−2𝑡(𝑥²+𝑧²−𝑥−𝑧). The values of 𝜃₀ and 𝜃₁ are determined by analytic solution using a spatial step ℎ=0.1 and a time step Δ𝑡=0.01 from 𝑡=0 to 𝑡=10. The Gaussian function exp(−𝑐𝑟²) is selected as the radial basis function and the error estimate is based on 𝐿₂-norm. The numerical and exact solutions of Example 1 are shown in Figures 2 and 3, respectively.

Figure 2 

The exact solutions of Example 1.

Figure 3 

Numerical solutions using the radial basis function collocation method.

Table 1 shows a comparison of this method with FEM, where ℎ=0.1 and 𝑡=10.

Table 1 

The results of the new method compared with FEM.

Next, we consider the time step Δ𝑡=0.01, a different spatial step, and different parameters 𝑐 in the Gaussian function exp(−𝑐𝑟²). The results are shown in Table 2.

Table 2 

Numerical solutions of Example 1 related to the parameter and spatial step.

For the linear model, Figures 1–4 and Table 1 show that the method presented in this paper is feasible. And comparing it with the traditional method, it showed good accuracy and rapid convergence rate. And Table 2 shows that the computational accuracy and computing time are both related to the selections of parameter 𝑐 in the radial basis functions and the spatial step.

Figure 4 

The error graph of the new method.

Example 2. We consider the 2D horizontal absorption soil water movement equation: 𝜕𝜃=𝜕𝜕𝑡𝐷(𝜃)𝜕𝑥𝜕𝜃+𝜕𝜕𝑥𝐷(𝜃)𝜕𝑧𝜕𝜃()[]×[],𝜃()𝜕𝑧+𝑓,𝑥,𝑧∈Ω=0,10,1𝑥,𝑧,0=𝜃₀,()𝜃()𝑥,𝑧∈Ω,𝑥,𝑧,𝑡=𝜃₁,(𝑥,𝑧)∈𝜕Ω,𝑡>0. (21) Here 𝐷(𝜃)=0.01𝑒^𝜃, 𝑓=−0.02𝑡²(𝑥−𝑥²)(𝑧−𝑧²)(𝑥²−𝑥+𝑧²−𝑧)−0.01𝑡²(𝑧²(2𝑥𝑧−𝑧−2𝑥+1)²) +𝑥²(2𝑥𝑧−2𝑧−𝑥+1)²+(𝑥−𝑥²)(𝑧²−𝑧) and the analytic solution is 𝜃(𝑧,𝑡)=(𝑥²−𝑥)(𝑧²−𝑧)𝑡. The values of 𝜃₀ and 𝜃₁ are determined by analytic solution. It is a nonlinear equation with the spatial step ℎ=0.1 and time step Δ𝑡=0.01 from 𝑡=0 to 𝑡=10. The radial basis function is also a Gaussian function exp(−𝑐𝑟²)  (𝑐>0). The error estimate is based on 𝐿₂-norm. The numerical and exact solutions of Example 2, along with the error, are shown in Figures 5, 6, and 7, respectively.

Figure 5 

The exact solutions of Example 2.

Figure 6 

Numerical solutions using the radial basis function collocation method.

Figure 7 

The error graph.

Table 3 shows the comparison of this method with FEM, where ℎ=0.1 and 𝑡=10

Table 3 

The results of the new method compared with FEM for Example 2.

For the time step Δ𝑡=0.01, we give the results of the new method based on different spatial steps and different parameters 𝑐. These are shown in Table 4.

Table 4 

Numerical solution based on different spatial steps and different parameters in Example 2.

For the nonlinear soil water movement equation, we can get the same conclusions as in linear equation.

4. Conclusions

In this paper, a mathematical model with the boundary conditions for soil water movement under irrigation has been developed by constructing a radial basis function collocation method. The existence and uniqueness of the solution were proven. Several numerical examples show that the proposed method yields higher precision and is easier to solve 2D soil water movement equations than traditional methods. Moreover, the selection of the time, spatial steps, and parameter 𝑐 has a direct influence on calculation accuracy. Therefore, it is necessary to study the combination of steps, radial basis function, and parameter 𝑐 to obtain the numerical solutions. In addition, it lays the foundation for the numerical solutions to high-dimensional soil water movement equations, which is very important.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This paper is supported by National Natural Science Foundation of China (NSFC 11101330, 51305344) and Scientific Foundation of Shaanxi Province (2014K05-22).