A Hybrid Localized Meshless Method for the Solution of Transient Groundwater Flow in Two Dimensions

: In this work, a hybrid localized meshless method is developed for solving transient groundwater ﬂow in two dimensions by combining the Crank–Nicolson scheme and the generalized ﬁnite difference method (GFDM). As the ﬁrst step, the temporal discretization of the transient groundwater ﬂow equation is based on the Crank–Nicolson scheme. A boundary value problem in space with the Dirichlet or mixed boundary condition is then formed at each time node, which is simulated by introducing the GFDM. The proposed algorithm is truly meshless and easy to program. Four linear or nonlinear numerical examples, including ones with complicated geometry domains, are provided to verify the performance of the developed approach, and the results illustrate the good accuracy and convergency of the method.


Introduction
As an important component of the hydrological cycle system, groundwater is a key source of domestic and industrial water supply. Therefore, the analysis of the groundwater flow has great significance for water supply security. Due to the complexity of the problem, an analytical solution is rarely available for most models of groundwater flow. With the development of computing techniques, more and more numerical approaches have been developed and applied to numerical simulations of science and engineering problems, such as the finite element method (FEM) [1][2][3], the boundary element method (BEM) [4,5], and the meshless method [6][7][8][9][10].
The GFDM, as a popular localized meshless collocation method, employs the Taylor series expansions and moving least squares (MLS) approximations [45,46] to form the system of algebraic equations with a spare matrix [47,48]. Thanks to this spare system, this method is highly efficient and suitable for the numerical simulations of large-scale problems. Many physical applications have been addressed by the GFDM, such as the thin elastic plate bending analysis [49], the electroelastic analysis of 3D piezoelectric structures [50], the acoustic wave propagation [51], the inverse Cauchy problem in 2D elasticity [52], the heat conduction problems [53], and the stationary flow in a dam [54].
A hybrid localized meshless method is proposed in this paper for the solution of transient groundwater flow in a two-dimensional space by combining the Crank-Nicolson scheme and the GFDM. As the first step, the Crank-Nicolson scheme is applied to the temporal discretization of the transient groundwater flow equation. At each time node, a boundary value problem in space is then formed and subsequently solved with the GFDM. Through the above process, a hybrid localized meshless approach is finally established, which is truly meshless and easy to program. The rest of the work is organized as follows. The governing equation of transient groundwater flow with boundary and initial conditions is described in Section 2. The formulations of the hybrid localized meshless method are derived in Section 3. Several linear and nonlinear numerical examples are provided in Section 4 to verify the performance of the developed method. Some conclusions are presented in Section 5.

Problem Definition
The movement of transient groundwater flow of a constant density in a homogeneous and anisotropic two-dimensional (2D) medium with the domain Ω can be described by using the following equation: where H denotes hydraulic head, W is the volumetric flow rate of a source or sink per unit volume, u s means the specific aquifer storativity, T x and T y are hydraulic conductivities along the xand y-axis, and t is the time.

Hybrid Localized Meshless Method
To solve the transient groundwater flow of Equations (1)-(4), the temporal discretization of this system is first made by using the Crank-Nicolson scheme. The spatial equation is then formed at each time node. The GFDM is finally used for the solution of the spatial equation with corresponding boundary conditions.

Temporal Discretization by the Crank-Nicolson Scheme
We insert n nodes t 1 = 0, t 2 , . . . , t n = T f in the time domain [0, T f ] where T f is the final time. By using the Crank-Nicolson scheme, the governing Equation (1) at each time node t i+1 is then recast as the following: where the time step size ∆t i = t i+1 − t i . Then we can reformulate Equation (5) as the following: As a result, a system of spatial equations with the boundary conditions (2) and (3) at time node t i+1 is formed and will be solved by using the GFDM.

Spatial Discretization by the GFDM
For the GFDM, some collocation nodes are first distributed in the computational domain Ω and its boundary Γ. A supporting domain called a star for each node x 0 = (x 0 , y 0 ) is defined by collecting m nearest nodes x j = (x j , y j )(j = 1, 2, . . . , m) as shown in Figure 1. In the star, x 0 and x j (j = 1, 2, . . . , m) are respectively named as the central node and the supporting nodes. This distance criterion is the simplest way of star selection. However, it should be noted that distorted (ill-conditioned) stars may be formed based on this distance criterion, particularly for cases with very irregular distributions of collocation nodes. To overcome the above drawback, the "cross" and the "Voronoi neighbours" criterions discussed in Refs. [55,56] can be used to form more reasonable stars.
As a result, a system of spatial equations with the boundary conditions (2) and (3) at time node 1 i t  is formed and will be solved by using the GFDM. are respectively named as the central node and the supporting nodes. This distance criterion is the simplest way of star selection. However, it should be noted that distorted (ill-conditioned) stars may be formed based on this distance criterion, particularly for cases with very irregular distributions of collocation nodes. To overcome the above drawback, the "cross" and the "Voronoi neighbours" criterions discussed in Refs. [55,56] can be used to form more reasonable stars.  To conveniently derive the system of algebraic equations, the H j (j = 0, 1, . . . , m) are employed to denote the values of hydraulic head at nodes x j (j = 0, 1, . . . , m) in the star of central node x 0 . Based on Taylor series expansions, H j can be given as

Spatial Discretization by the GFDM
with a j = x j − x 0 , and b j = y j − y 0 .
By truncating the expansion (7) after second-order derivatives of hydraulic head H, we can define a residual function R(H) as the following: with the following weighting function κ j [57,58]: in which λ j = x j − x 0 , and λ m = max λ j , j = 1, 2, . . . , m . It should be noted that the weighting function in the GFDM should be a monotonic decreasing function of λ j . Because the Taylor series approximation becomes more accurate when the distance λ j is smaller, which should have a higher weight κ j in residual function R(H) of Equation (9). Some other weighting functions can be found in [53,59]. A vector P H is defined by the following: Minimizing the residual function R(H) with respect to each element in the vector P H , i.e.,

∂R(H)
we can have a system of linear equations as the following: with the following: where is a 5 × 5 square matrix that all elements are one, diag E (j) 1 and diag E (j) 2 are both diagonal matrices with their diagonal elements as the following: and With the help of Equation (13) and Equation (15), the vector P H can be formulated as the following: As a result, all the first-and second-order derivatives of hydraulic head H at central node x 0 are expressed as the linear combinations of values H j (j = 0, 1, . . . , m). By substituting the corresponding second-order derivatives of Equation (17) into Equation (6), we can easily recast Equation (6) as the following: where α j (j = 0, 1, . . . , m) and Q are obviously known and can be determined by Equations (6) and (17). For each collocation node excepting that of satisfying boundary condition (2), one equation can be obtained by using the above similar derivation. It should be noted that H 0 = g(x 0 , y 0 , t i+1 ) is directly used as one equation for collocation nodes of satisfying boundary condition (2). Finally, by using the GFDM, the spatial equation with the corresponding boundary conditions has been transformed into a system of linearly algebraic equations with a sparse matrix. The values of hydraulic head H at all collocation nodes can be calculated once this system is solved. In addition, we consider a case with a nonlinear volumetric flow rate in the following numerical example 4 to further verify the proposed approach. Through a similar derivation process, a nonlinear system of algebraic equations is established and solved by using "fsolve" function of MATLAB.

Numerical Examples
In this section, four numerical examples with square and complicated domains are given to test the accuracy and stability of the proposed method. To preferably estimate the precision of numerical results, two different error definitions are provided as the following [60,61]: (20) in which N denotes the number of collocation nodes, H i and H i represents the exact and numerical results of hydraulic head at the i-th node. Unless otherwise specified, the number of supporting nodes in a star is set to m = 12 in all numerical examples.

Example 1: Hydraulic Head Distribution in a Square Domain
As the first numerical example, we consider the distribution of hydraulic heads in a unit square domain with its central point (0.5, 0.5). The specific aquifer storativity is u s = 1, and the hydraulic conductivities are T x = 1 and T y = 3. The volumetric flow rate W is given as the following: The Dirichlet boundary condition H(x, y, t) = 0 is imposed on the boundaries x = 0, and y = 0. The Neumann boundary condition ∂H(x,y,t) ∂n = 0 is applied to the boundaries x = 1, and y = 1. The initial condition is H(x, y, 0) = 0. Based on these, the exact solution can be determined as the following: In the numerical simulation, the final time is set to T f = 2, and the time step size is ∆t = 0.05. For the spatial discretization, 396 collocation nodes are distributed in the domain and its boundary, which have the following two different patterns: regular distribution (case 1) and irregular distribution (case 2), as shown in Figure 2a. From t = 0 to t = 2, the variations of global and max errors of hydraulic head H calculated by the GFDM with the Crank-Nicolson scheme (CN-GFDM) are plotted in Figure 2b. As we can see from this figure, the developed method has good performance for two different collocation node distributions. We also find that errors obtained by employing the irregular distribution (case 2) are larger than those obtained by using the regular distribution (case 1).  (20) in which N denotes the number of collocation nodes, i The Dirichlet boundary condition In the numerical simulation, the final time is set to 2 f T  , and the time step size is 0.05 t   . For the spatial discretization, 396 collocation nodes are distributed in the domain and its boundary, which have the following two different patterns: regular distribution (case 1) and irregular distribution (case 2), as shown in   Figure 2b. As we can see from this figure, the developed method has good performance for two different collocation node distributions. We also find that errors obtained by employing the irregular distribution (case 2) are larger than those obtained by using the regular distribution (case 1). To investigate the convergence behavior of the CN-GFDM, the program is rerun by only changing the number of collocation nodes compared with the previous setting. Here, the distribution of nodes is regular. By using the developed method with a different number of collocation nodes (or a different mean distance of collocation nodes), Table 1 provides the max and global errors of the numerical results of the hydraulic head at time t = 2. It can be obviously observed that the errors decay with an increasing collocation node number, which indicates the CN-GFDM has a good convergence property for this case.

Example 2: Hydraulic Head Distribution in a Heart-Shaped Domain
The second numerical example is a hydraulic head distribution in a heart-shaped domain, and the dimension of this domain is shown in Figure 3. The specific aquifer storativity and hydraulic conductivity are assumed to be as the following: The volumetric flow rate is W = 0. In this case, the Dirichlet boundary and initial conditions are imposed as the following: H(x, y, t) = e −t cos x cos y + c, (x, y) ∈ Γ, H(x, y, 0) = cos x cos y + c, (x, y) ∈ Ω, where c = 0.2. The exact solution for this example is determined as H(x, y, t) = e −t cos x cos y + c. To investigate the convergence behavior of the CN-GFDM, the program is rerun by only changing the number of collocation nodes compared with the previous setting. Here, the distribution of nodes is regular. By using the developed method with a different number of collocation nodes (or a different mean distance of collocation nodes), Table 1 provides the max and global errors of the numerical results of the hydraulic head at time 2 t  . It can be obviously observed that the errors decay with an increasing collocation node number, which indicates the CN-GFDM has a good convergence property for this case.

Example 2: Hydraulic Head Distribution in a Heart-Shaped Domain
The second numerical example is a hydraulic head distribution in a heart-shaped domain, and the dimension of this domain is shown in Figure 3. The specific aquifer storativity and hydraulic conductivity are assumed to be as the following: The volumetric flow rate is 0 W  . In this case, the Dirichlet boundary and initial conditions are imposed as the following:  To simulate the solution from t = 0 to t = 5, the developed method employs 218, 464 (see Figure 4), and 1700 collocation nodes. The time step size is set to ∆t = 0.05. Figure 5 provides the contours of relative errors (RE) of the hydraulic head at final time t = 5 by using these three distributions of collocation nodes. It can be found that the max relative error of numerical results in the computational domain is less than 3 × 10 −4 even only using   Next, we keep the above-mentioned setting and distribute 1700 collocation nodes. The mean distance of these nodes is 0.101. Table 2 lists the max and global errors of hydraulic head H at final time 5 t  , which are calculated by the CN-GFDM with different time step sizes.
As we can observe, the errors decay rapidly when decreasing the time step size. . Table 3 provides the variations of two kinds of errors and CPU time with different numbers of supporting nodes. As we can observe, the accuracy of numerical results obtained by the present method is relatively insensitive to the number of supporting   Next, we keep the above-mentioned setting and distribute 1700 collocation nodes. The mean distance of these nodes is 0.101. Table 2 lists the max and global errors of hydraulic head H at final time 5 t  , which are calculated by the CN-GFDM with different time step sizes.
As we can observe, the errors decay rapidly when decreasing the time step size. . Table 3 provides the variations of two kinds of errors and CPU time with different numbers of supporting nodes. As we can observe, the accuracy of numerical results obtained by the present method is relatively insensitive to the number of supporting Next, we keep the above-mentioned setting and distribute 1700 collocation nodes. The mean distance of these nodes is 0.101. Table 2 lists the max and global errors of hydraulic head H at final time t = 5, which are calculated by the CN-GFDM with different time step sizes. As we can observe, the errors decay rapidly when decreasing the time step size. Finally, we investigate the influence of supporting node number m on the precision and computational efficiency of the developed approach. The time step size is reset as ∆t = 0.25. Table 3 provides the variations of two kinds of errors and CPU time with different numbers of supporting nodes. As we can observe, the accuracy of numerical results obtained by the present method is relatively insensitive to the number of supporting nodes. To have a higher computing efficiency of the CN-GFDM, we should choose a relatively small number of supporting nodes when the numerical results satisfy the accuracy requirement.

Example 3: Hydraulic Head Distribution in a Complicated Domain
As the third numerical example, the distribution of hydraulic head in a complicated domain is considered. The dimension of this domain is 2.7 × 1.4, as shown in Figure 6. The specific aquifer storativity is set to u s = 0.3, and hydraulic conductivities are assumed to be functions as the following: The volumetric flow rate is given as the following: W(x, y, t) = e x+y u s cos t − (T x + T y ) sin t The Dirichlet boundary condition is imposed as the following: H(x, y, t) = e x+y sin t, (x, y) ∈ Γ, and initial condition is H = 0. The exact solution for this example is H(x, y, t) = e x+y sin t.

Example 3: Hydraulic Head Distribution in a Complicated Domain
As the third numerical example, the distribution of hydraulic head in a comp domain is considered. The dimension of this domain is 2.7 1.4  , as shown in F The specific aquifer storativity is set to 0.3 s u  , and hydraulic conductivities sumed to be functions as the following: , and .
x y x y T e T e   The volumetric flow rate is given as the following: The Dirichlet boundary condition is imposed as the following: ( , , ) sin , ( , ) , sin .
x y H x y t e t   Figure 6. The dimension of a complicated domain.
We first consider the simulation from 0 t  to 5 t  , and set the time step 0.1 t   . By using the present approach with 3398 (see Figure 7) and 8918 coll nodes, Figures 8 and 9 respectively plot the contours of relative errors of hydraul  We first consider the simulation from t = 0 to t = 5, and set the time step size as ∆t = 0.1. By using the present approach with 3398 (see Figure 7) and 8918 collocation nodes, Figures 8 and 9 respectively plot the contours of relative errors of hydraulic head H and its flux ∂H ∂x at final time t = 5. The numerical results in these figures illustrate the availability and convergency of the developed CN-GFDM.
Finally, a long-time simulation of hydraulic head H from t = 0 to t = 100 is considered. The number of collocation nodes is 8918, and the time step size is ∆t = 0.1. Figure 10 shows the max and global errors of hydraulic head H, which are changed as functions of time. As we can see from this figure, the two kinds of errors both remain stable in this simulation. . Figure 10 shows the max and global errors of hydraulic head H, which are changed as functions of time. As we can see from this figure, the two kinds of errors both remain stable in this simulation.  Figure 10 shows the max and global errors of hydraulic head H, which are changed as functions of time. As we can see from this figure, the two kinds of errors both remain stable in this simulation.

Example 4: Nonlinear Hydraulic head Distribution in a Gear Domain
As the final numerical example, we consider the distribution of nonlinea head in a gear domain. Figure 11 shows the dimension of the gear domain. T aquifer storativity is assumed to be 1200 s u  , and hydraulic conductivities lowing: 2 2 , and .
x y T x T y   The volumetric flow rate is a nonlinear term given as the following:  

Example 4: Nonlinear Hydraulic Head Distribution in a Gear Domain
As the final numerical example, we consider the distribution of nonlinear hydraulic head in a gear domain. Figure 11 shows the dimension of the gear domain. The specific aquifer storativity is assumed to be u s = 1200, and hydraulic conductivities are the following: , and .
x y T x T y    Figure 11. The dimension of a gear domain.
The numerical simulation for this case is run from size is set to 0.05 t   , and 1186 collocation nodes are distributed in the The volumetric flow rate is a nonlinear term given as the following: The Dirichlet boundary condition is imposed as the following: and initial condition is the following: H(x, y, t) = x 4 + y 4 + 0.1, (x, y) ∈ Ω.
The exact solution for this case is determined as the following: H(x, y, t) = x 4 + y 4 e t 100 + 0.1 The numerical simulation for this case is run from t = 0 to t = 10. The time step size is set to ∆t = 0.05, and 1186 collocation nodes are distributed in the gear domain and its boundary as shown in Figure 12. The contours of relative errors of the hydraulic head H at t = 5 and t = 10 are plotted in Figure 13. As we can see from this figure, the present approach obtains the satisfied numerical results for this nonlinear case, and max relative error is less than 4 × 10 −3 . Figure 14 provides the max and global errors of hydraulic head H at each time node, which illustrates that the developed CN-GFDM yields accurate numerical results as an increasing time.
Mathematics 2022, 10, 515 12 of 15 and its boundary as shown in Figure 12. The contours of relative errors of the hydraulic head H at 5 t  and 10 t  are plotted in Figure 13. As we can see from this figure, the present approach obtains the satisfied numerical results for this nonlinear case, and max relative error is less than 3 4 10   . Figure 14 provides the max and global errors of hydraulic head H at each time node, which illustrates that the developed CN-GFDM yields accurate numerical results as an increasing time.    and its boundary as shown in Figure 12. The contours of relative errors of the hydraulic head H at 5 t  and 10 t  are plotted in Figure 13. As we can see from this figure, the present approach obtains the satisfied numerical results for this nonlinear case, and max relative error is less than 3 4 10   . Figure 14 provides the max and global errors of hydraulic head H at each time node, which illustrates that the developed CN-GFDM yields accurate numerical results as an increasing time.

Concluding Remarks
To simulate the transient groundwater flow in homogeneous and anisotropic twodimensional mediums, a hybrid localized meshless method is constructed based on the Crank-Nicolson scheme for temporal discretization and the GFDM for spatial discretization. The present approach is truly meshless and easy to program.
To fully investigate the performance of the developed method, the max and global errors of hydraulic head are both provided for numerical examples with different boundary conditions, complicated geometry domains, and several kinds of hydraulic conductivities. Numerical results indicate that the hybrid localized meshless method developed in this work obtains the satisfied accuracy and convergency in time and space.

Funding:
The work described in this paper was supported by the Government of Shandong's Education System for studying abroad, and the key R & D plan of Zibo city: Intelligent Water IoT cloud service platform of "Zishui Online" (NO. 2019ZBXC246).
Institutional Review Board Statement: Not applicable. Figure 14. Two types of error curves of hydraulic head H from t = 0 to t = 10.

Concluding Remarks
To simulate the transient groundwater flow in homogeneous and anisotropic two-dimensional mediums, a hybrid localized meshless method is constructed based on the Crank-Nicolson scheme for temporal discretization and the GFDM for spatial discretization. The present approach is truly meshless and easy to program.
To fully investigate the performance of the developed method, the max and global errors of hydraulic head are both provided for numerical examples with different boundary conditions, complicated geometry domains, and several kinds of hydraulic conductivities. Numerical results indicate that the hybrid localized meshless method developed in this work obtains the satisfied accuracy and convergency in time and space.