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

Abstract

In this work, a hybrid localized meshless method is developed for solving transient groundwater flow in two dimensions by combining the Crank–Nicolson scheme and the generalized finite difference method (GFDM). As the first step, the temporal discretization of the transient groundwater flow 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.

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

As a new approach in recent years, the meshless method is now widely applied in various fields [11,12,13,14,15,16,17,18,19], particularly in computational fluid dynamics (CFD). The developed meshless approaches can be classified into collocation-based and Galerkin-based methods. Compared with the latter, the meshless collocation methods have the advantages of no numerical quadrature and mesh generation, and some of these are the localized method of fundamental solutions (LMFS) [20,21,22], the generalized finite difference method (GFDM) [23,24,25,26,27,28,29,30,31,32,33], the localized Chebyshev collocation method [34], the singular boundary method (SBM) [35,36,37,38,39,40,41,42,43], and the localized knot method (LKM) [44].

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.

2. Problem Definition

The movement of transient groundwater flow of a constant density in a homogeneous and anisotropic two-dimensional (2D) medium with the domain $\mathsf{\Omega }$ can be described by using the following equation:

$$T_x\frac{{\partial }^{2}H(x,y,t)}{\partial x^2}+T_y\frac{{\partial }^{2}H(x,y,t)}{\partial y^2}+W(x,y,t)=u_s\frac{\partial H(x,y,t)}{\partial t},(x,y)\in \mathsf{\Omega },t>0,$$

(1)

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 x- and y-axis, and t is the time.

To obtain the solution of Equation (1), the boundary and initial conditions are imposed as the following:

$$H(x,y,t)=g_1(x,y,t),\hspace{1em}(x,y)\in \partial {\mathsf{\Omega }}_D,t\ge 0,$$

(2)

$$\frac{\partial H(x,y,t)}{\partial \mathit{n}(x,y)}=g_2(x,y,t),\hspace{1em}(x,y)\in \partial {\mathsf{\Omega }}_N,t\ge 0,$$

(3)

$$H(x,y,0)=g_3(x,y),\hspace{1em}(x,y)\in \mathsf{\Omega },$$

(4)

where $g_i(i=1,2,3)$ are given functions, $\partial {\mathsf{\Omega }}_D\cup \partial {\mathsf{\Omega }}_N=\partial \mathsf{\Omega }(\partial {\mathsf{\Omega }}_D\cap \partial {\mathsf{\Omega }}_N=\varnothing )$, and $\mathit{n}$ is unit outward normal vector to $\partial {\mathsf{\Omega }}_N$.

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

3.1. Temporal Discretization by the Crank–Nicolson Scheme

We insert $n$ nodes $\{t_1=0,t_2,\dots ,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:

$$\begin{array}{l}\frac{1}{2}\left[T_x\frac{{\partial }^{2}H(x,y,t_{i+1})}{\partial x^2}+T_y\frac{{\partial }^{2}H(x,y,t_{i+1})}{\partial y^2}+W(x,y,t_{i+1})+T_x\frac{{\partial }^{2}H(x,y,t_i)}{\partial x^2}+T_y\frac{{\partial }^{2}H(x,y,t_i)}{\partial y^2}+W(x,y,t_i)\right]\\ =u_s\frac{H(x,y,t_{i+1})-H(x,y,t_i)}{\Delta t_i},i=1,2,\dots ,n-1,\end{array}$$

(5)

where the time step size $\Delta t_i=t_{i+1}-t_i$. Then we can reformulate Equation (5) as the following:

$$\begin{array}{l}\frac{\Delta t_i}{2}\left[T_x\frac{{\partial }^{2}H(x,y,t_{i+1})}{\partial x^2}+T_y\frac{{\partial }^{2}H(x,y,t_{i+1})}{\partial y^2}\right]-u_{s}H(x,y,t_{i+1})=\\ =-u_{s}H(x,y,t_i)-\frac{\Delta t_i}{2}\left[T_x\frac{{\partial }^{2}H(x,y,t_i)}{\partial x^2}+T_y\frac{{\partial }^{2}H(x,y,t_i)}{\partial y^2}+W(x,y,t_i)+W(x,y,t_{i+1})\right],i=1,2,\dots ,n-1,\end{array}$$

(6)

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.

3.2. Spatial Discretization by the GFDM

For the GFDM, some collocation nodes are first distributed in the computational domain $\mathsf{\Omega }$ and its boundary Γ. A supporting domain called a star for each node ${\mathit{x}}_0=(x_0,y_0)$ is defined by collecting $\mathit{m}$ nearest nodes ${\mathit{x}}_j=(x_j,y_j)(j=1,2,\dots ,m)$ as shown in Figure 1. In the star, ${\mathit{x}}_0$ and ${\mathit{x}}_j(j=1,2,\dots ,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.

To conveniently derive the system of algebraic equations, the $H_j(j=0,1,\dots ,m)$ are employed to denote the values of hydraulic head at nodes ${\mathit{x}}_j(j=0,1,\dots ,m)$ in the star of central node ${\mathit{x}}_0$. Based on Taylor series expansions, $H_j$ can be given as

$$H_j=H_0+a_j\frac{\partial H_0}{\partial x}+b_j\frac{\partial H_0}{\partial y}+\frac{1}{2}\left(a_j^2\frac{{\partial }^2{H}_0}{\partial x^2}+b_j^2\frac{{\partial }^2{H}_0}{\partial y^2}\right)+a_j{b}_j\frac{{\partial }^2{H}_0}{\partial x\partial y}+\dots \dots ,\hspace{1em}j=1,2,\dots ,m,$$

(7)

with

$$a_j=x_j-x_0,\mathrm{and}{b}_j=y_j-y_0.$$

(8)

By truncating the expansion (7) after second-order derivatives of hydraulic head $H$, we can define a residual function $R(H)$ as the following:

$$R(H)={{\displaystyle \sum _{j=1}^m\left[\left(H_0-H_j+a_j\frac{\partial H_0}{\partial x}+b_j\frac{\partial H_0}{\partial y}+\frac{a_j^2}{2}\frac{{\partial }^2{H}_0}{\partial x^2}+\frac{b_j^2}{2}\frac{{\partial }^2{H}_0}{\partial y^2}+a_j{b}_j\frac{{\partial }^2{H}_0}{\partial x\partial y}\right){\kappa }_j\right]}}^2,$$

(9)

with the following weighting function ${\kappa }_j$ [57,58]:

$${\kappa }_j=\frac{\exp \left(-{\lambda }_j^2\right)-\exp \left(-{\lambda }_m^2\right)}{1-\exp \left(-{\lambda }_m^2\right)},\hspace{1em}j=1,2,\dots ,m,$$

(10)

in which ${\lambda }_j=\left|{\mathit{x}}_j-{\mathit{x}}_0\right|$, and ${\lambda }_m=\max \left\{{\lambda }_j,j=1,2,\dots ,m\right\}$. It should be noted that the weighting function in the GFDM should be a monotonic decreasing function of ${\lambda }_j$. Because the Taylor series approximation becomes more accurate when the distance ${\lambda }_j$ is smaller, which should have a higher weight ${\kappa }_j$ in residual function $R(H)$ of Equation (9). Some other weighting functions can be found in [53,59].

A vector ${\mathit{P}}_H$ is defined by the following:

$${\mathit{P}}_H={\left[\frac{\partial H_0}{\partial x},\frac{\partial H_0}{\partial y},\frac{{\partial }^2{H}_0}{\partial x^2},\frac{{\partial }^2{H}_0}{\partial y^2},\frac{{\partial }^2{H}_0}{\partial x\partial y}\right]}^T$$

(11)

Minimizing the residual function $R(H)$ with respect to each element in the vector ${\mathit{P}}_H$, i.e.,

$$\frac{\partial R(H)}{\partial \left\{\frac{\partial H_0}{\partial x}\right\}}=0,\frac{\partial R(H)}{\partial \left\{\frac{\partial H_0}{\partial y}\right\}}=0,\frac{\partial R(H)}{\partial \left\{\frac{{\partial }^2{H}_0}{\partial x^2}\right\}}=0,\frac{\partial R(H)}{\partial \left\{\frac{{\partial }^2{H}_0}{\partial y^2}\right\}}=0,\frac{\partial R(H)}{\partial \left\{\frac{{\partial }^2{H}_0}{\partial x\partial y}\right\}}=0,$$

(12)

we can have a system of linear equations as the following:

$$\mathit{A}{\mathit{P}}_H=\mathit{b},$$

(13)

with the following:

$$\mathit{A}={\displaystyle \sum _{j=1}^{m}diag\left({\mathit{E}}_1^{(j)}\right)\coprod diag\left({\mathit{E}}_2^{(j)}\right)},$$

(14)

$$\mathit{b}=\mathit{B}\mathit{H}={\left(-{\displaystyle \sum _{j=1}^m{\kappa }_j{\mathit{E}}_2^{(j)}},{\kappa }_1{\mathit{E}}_2^{(1)},{\kappa }_2{\mathit{E}}_2^{(2)},\cdots ,{\kappa }_m{\mathit{E}}_2^{(m)}\right)}_{5\times (m+1)}{\left(\begin{array}{l}{H}_0\\ H_1\\ H_2\\ ⋮\\ H_m\end{array}\right)}_{(m+1)\times 1},$$

(15)

where $\coprod$ is a $5\times 5$ square matrix that all elements are one, $diag\left({\mathit{E}}_1^{(j)}\right)$ and $diag\left({\mathit{E}}_2^{(j)}\right)$ are both diagonal matrices with their diagonal elements as the following:

$${\mathit{E}}_1^{(j)}={\kappa }_j{\left[\begin{array}{ccccc}{a}_j& b_j& a_j^2& b_j^2& a_j{b}_j\end{array}\right]}^T,\mathrm{and}{\mathit{E}}_2^{(j)}=\frac{{\kappa }_j}{2}{\left[\begin{array}{ccccc}2a_j& 2b_j& a_j^2& b_j^2& 2a_j{b}_j\end{array}\right]}^T,$$

(16)

and $\mathit{H}={\left[H_0,H_1,H_2,\dots ,H_m\right]}^T$.

With the help of Equation (13) and Equation (15), the vector ${\mathit{P}}_H$ can be formulated as the following:

$${\mathit{P}}_H={\mathit{A}}^{-1}\mathit{B}{\left[H_0,H_1,H_2,\dots ,H_m\right]}^T.$$

(17)

As a result, all the first- and second-order derivatives of hydraulic head $H$ at central node ${\mathit{x}}_0$ are expressed as the linear combinations of values $H_j(j=0,1,\dots ,m)$. By substituting the corresponding second-order derivatives of Equation (17) into Equation (6), we can easily recast Equation (6) as the following:

$${\displaystyle \sum _{j=0}^m{\alpha }_j{H}_j}=Q,$$

(18)

where ${\alpha }_j(j=0,1,\dots ,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.

4. 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]:

$$\mathrm{Global}\mathrm{error}=\sqrt{{\displaystyle \sum _{i=1}^N{({\overline{H}}_i-H_i)}^2}}/\sqrt{{\displaystyle \sum _{i=1}^N{H}_i{}^2}},$$

(19)

$$\mathrm{Max}\mathrm{error}=\underset{1\le i\le N}{\max }\left\{\frac{\left|{\overline{H}}_i-H_i\right|}{\left|H_i\right|}\right\}$$

(20)

in which $N$ denotes the number of collocation nodes, $H_i$ and ${\overline{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.

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

$$W(x,y,t)=\sin \left(0.5\pi x\right)\sin \left(0.5\pi y\right)\left({\pi }^2\sin t+\cos t\right).$$

(21)

The Dirichlet boundary condition $H(x,y,t)=0$ is imposed on the boundaries $x=0,\mathrm{and}y=0$. The Neumann boundary condition $\frac{\partial H(x,y,t)}{\partial \mathit{n}}=0$ is applied to the boundaries $x=1,\mathrm{and}y=1$. The initial condition is $H(x,y,0)=0$. Based on these, the exact solution can be determined as the following:

$$H(x,y,t)=\sin (0.5\pi x)\sin (0.5\pi y)\sin t.$$

(22)

In the numerical simulation, the final time is set to $T_f=2$, and the time step size is $\Delta 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).

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. Max and global errors of hydraulic head H with different number of collocation nodes.

Number of collocation nodes	21	96	192	285	396
Mean distance of collocation nodes	$2.50\times {10}^{-1}$	$1.11\times {10}^{-1}$	$7.69\times {10}^{-2}$	$6.25\times {10}^{-2}$	$5.26\times {10}^{-2}$
Max error	$1.01\times {10}^{-2}$	$3.75\times {10}^{-3}$	$2.05\times {10}^{-3}$	$1.42\times {10}^{-3}$	$1.04\times {10}^{-3}$
Global error	$5.94\times {10}^{-3}$	$2.74\times {10}^{-3}$	$1.57\times {10}^{-3}$	$1.11\times {10}^{-3}$	$8.18\times {10}^{-4}$

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.

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

$$u_s=3,T_x=2T_y=2.$$

(23)

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,\hspace{1em}(x,y)\in \Gamma ,$$

(24)

$$H(x,y,0)=\cos x\cos y+c,\hspace{1em}(x,y)\in \mathsf{\Omega },$$

(25)

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 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 $\Delta 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\times {10}^{-4}$ even only using 218 collocation nodes. Moreover, all numerical errors overall decreased with an increasing number of collocation nodes.

Next, we keep the above-mentioned setting and distribute 1700 collocation nodes. The mean distance of these nodes is 0.101.

Table 2. Max and global errors of hydraulic head H with different time step size.

Time step size	0.625	0.50	0.25	0.125	0.05
Max error	$2.51\times {10}^{-3}$	$1.60\times {10}^{-3}$	$4.05\times {10}^{-4}$	$8.42\times {10}^{-5}$	$1.17\times {10}^{-5}$
Global error	$1.20\times {10}^{-3}$	$7.80\times {10}^{-4}$	$1.91\times {10}^{-4}$	$3.89\times {10}^{-5}$	$4.97\times {10}^{-6}$

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 $\Delta t=0.25$.

Table 3. Max and global errors of hydraulic head H and CPU time with different supporting node numbers.

m	12	16	20	24	28
Max error	$4.05\times {10}^{-4}$	$3.53\times {10}^{-4}$	$3.45\times {10}^{-4}$	$3.38\times {10}^{-4}$	$3.23\times {10}^{-4}$
Global error	$1.91\times {10}^{-4}$	$1.71\times {10}^{-4}$	$1.66\times {10}^{-4}$	$1.61\times {10}^{-4}$	$1.55\times {10}^{-4}$
CPU time (s)	0.27	0.34	0.38	0.44	0.49

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.

4.3. 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\times 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:

$$T_x=e^{\left|x\right|},\mathrm{and}{T}_y=e^{\left|y\right|}.$$

(26)

The volumetric flow rate is given as the following:

$$W(x,y,t)=e^{x+y}\left[u_s\cos t-(T_x+T_y)\sin t\right]$$

(27)

The Dirichlet boundary condition is imposed as the following:

$$H(x,y,t)=e^{x+y}\sin t,\hspace{1em}(x,y)\in \Gamma ,$$

(28)

and initial condition is $H=0.$ The exact solution for this example is $H(x,y,t)=e^{x+y}\sin t.$

We first consider the simulation from $t=0$ to $t=5$, and set the time step size as $\Delta t=0.1$. By using the present approach with 3398 (see Figure 7) and 8918 collocation nodes, Figure 8 and Figure 9 respectively plot the contours of relative errors of hydraulic head H and its flux $\frac{\partial H}{\partial 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 $\Delta 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.

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

$$T_x=x^2,\mathrm{and}{T}_y=y^2.$$

(29)

The volumetric flow rate is a nonlinear term given as the following:

$$W(x,y,t)=H^2(x,y,t)-{\left[\left(x^4+y^4\right)e^{\frac{t}{100}}+0.1\right]}^2$$

(30)

The Dirichlet boundary condition is imposed as the following:

$$H(x,y,t)=\left(x^4+y^4\right)e^{\frac{t}{100}}+0.1,\hspace{1em}(x,y)\in \Gamma ,$$

(31)

and initial condition is the following:

$$H(x,y,t)=x^4+y^4+0.1,\hspace{1em}(x,y)\in \mathsf{\Omega }.$$

(32)

The exact solution for this case is determined as the following:

$$H(x,y,t)=\left(x^4+y^4\right)e^{\frac{t}{100}}+0.1$$

(33)

The numerical simulation for this case is run from $t=0$ to $t=10$. The time step size is set to $\Delta 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\times {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.

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