Numerical model of seepage flows by reformulating finite element method based on new spherical Hankel shape functions

The water penetration in soil is investigated numerically using the finite element method (FEM) in a novel way. In the suggested method, new spherical Hankel shape functions are used and the finite element method is reformulated based on them. These new functions are obtained from the first and second kind of Bessel functions. The properties of Hankel shape functions lead to having more accuracy and robustness for the proposed method with low number of elements. To validate the suggested approach, at first, a boundary value problem is solved and the results are compared with the available analytical solution. Then, in order to prove the efficiency and applicability of the present model in the seepage problems, five examples including saturated and unsaturated flow in porous media are studied and the hydraulic head is calculated. Afterward, the results obtained from the classical and new method are compared together. The comparisons indicate that the suggested method with the low number of elements is more precise than the classic FEM with the same mesh.


Introduction
The seepage phenomenon is an important issue in soil mechanics and water engineering and generally is defined as the slow moving of a fluid in a porous medium. This problem is related to various fields such as water seepage through body of a dam and its foundation and also groundwater flows. To study these problems, the governing equation of the boundary value problems should be solved. This equation is obtained by combining Darcy's and continuity laws. By solving the mentioned equation, the potential variable or hydraulic head is calculated.
The numerical study of the seepage problems has been considered by many researchers. The analytical methods are applied when one problem has simple geometry and boundary conditions. However, when the computational domain boundary conditions are complicated, the numerical methods are more efficient and usually have suitable precision for the modeling. In the literature, several works have been done in the field of leakage flow problems by the numerical methods. Juanes (2005) presented a variational multiscale finite element method for multiphase flow in a porous media. Fu and Sheng (2009) presented a numerical model for simulation of the unsteady seepage flow through dam. These authors verified the reliability and accuracy of this model by laboratory data. A threedimensional numerical manifold method for unconfined seepage analysis was proposed by Jiang et al. (2010). In this method, these researchers used tetrahedral finite element meshes. Bereslavskii (2011) simulated the plan steady state seepage in a homogeneous isotropic ground from channel as a hydrodynamic formulation. Alyavuz (2012) used a new computational analysis method known as discrete singular convolution (DSC) method for solving seepage flow problems. A numerical modeling was carried out for groundwater flow in Daxing, China, by Yang et al. (2012). Kazemzadeh-parsi and Daneshmand (2013) solved the 3D unconfined seepage problems using smoothed fixed grid finite element method. The transient seepage in an unsaturated porous media was solved by Pedroso (2015) using finite element solution. He also predicted the free surface location in the time-dependent problems. Athani  et al. (2015) presented the results of seepage and stability analysis of an earth dam using finite element method. A groundwater model was developed by Karay and Hajnal (2015) for flow in fractured rocks. These authors considered the non-laminar flow in discontinuities. Fukuchi (2016) studied the steady state seepage problems using interpolation finite difference method. He was able to solve two-and three-dimensional elliptic partial differential equations over the complicated domain. He also achieved results with high accuracy. The 2D subsurface seepage flow in an anisotropic porous medium was simulated by Lande et al. (2016). In this simulation, a new analytical solution was developed to estimate the transient behavior of the phreatic surface. Yuan and Zhong (2016) proposed a weak form quadrature element formulation in three-dimension for analysis of unconfined seepage in an earth dam. These authors validated the numerical solution by comparing with the available analytical and numerical solutions. Jiang et al. (2017) simulated the groundwater bypass seepage at the dam site of Dongzhiang hydro-junction. They considered Hydrological conditions include underground water divides, a hanging river and hydraulic drops along faults.
In this paper, it is developed a new finite element method for modeling the seepage problems. In this present approach, spherical Hankel shape functions are used instead of Lagrange shape functions in finite element method. Generally, in the classical FEM, when the number of elements are increased or in other words, when the mesh used is refined, the numerical solution become close to the exact solution, but by doing so, the CPU time and storage space are enhanced. Therefore, by using the new suggested shape functions, the degrees of freedom can be decreased while high accuracy can be resulted.
In order to estimate multivariable functions by linear combinations of terms depending on a univariate function, radial basis functions (RBFs) are used (Buhmann 2003). RBFs are categorized into two classes: oscillatory and nonoscillatory. The Gaussian functions, the thin plate splines, compact supported functions, multiquadrics and inverse multiquadrics are in the non-oscillatory class, while real and complex Fourier and J-Bessel RBFs are placed in oscillatory category (Hamzeh Javaran et al. 2011a, 2011bKhaji 2012, 2014;Khaji and Hamzehei Javaran 2013;Hamzehei Javaran and Shojaee 2017;Bahrampour et al. 2018;Farmani et al. 2018).
Complex Fourier and J-Bessel RBFs were proposed by Hamzeh Javaran et al. (2011a). The J-Bessel RBFs only have the properties of the first kind of Bessel function, while spherical Hankel shape functions contain the features of both first and second kind of Bessel functions, simultaneously. Therefore, this is the reason of the robustness of Hankel functions for increasing the solution accuracy. In the area of the application of Hankel shape functions, Mousavi et al. (2021) solved fracture mechanics problem using boundary element method based on these functions recently.
In order to apply Hankel functions in the fluid mechanics, first, Farmani et al. (2018) made an improvement in the solution of Navier-Stokes equations for the benchmark fluid mechanics problems. Then, Hankel shape functions are used for boundary value problems (Farmani et al. 2019a) and for the sloshing behavior in 2-D tanks (Farmani et al. 2019b). Recently, Hankel shape functions are used for modeling the performance of Tuned Liquid Dampers (TLD) and also the interaction force between TLD-structure (Farmani et al. 2021). According to the past researches, it is useful to investigate the performance of the Hankel shape functions in the seepage and groundwater flow modeling in this paper.
The present paper consists of the following parts: at first, the governing equation and FEM formulation are stated. Then, spherical Hankel shape function and enrichment procedure are investigated. After that, to validate the proposed method for the leakage flows problems, one boundary value problem with different boundary conditions is solved and the results obtained are compared with the analytical solution. Also, in order to show applicability and reliability, five numerical problems are studied including saturated and unsaturated flow in porous media. Finally, the results are analyzed and the accuracy and precision of the present model are demonstrated.

Methods
The governing equations and FEM formulation are discussed in this part. Then, the suggested functions and the enrichment procedure are stated.

Governing equation and FEM formulation
In order to study the behavior of water in the porous media, two-dimensional seepage flow is considered. Therefore, the governing equation can be written as Eq. (1) which is obtained from Darcy's and continuity laws (Lam et al. 1987).
The above relation is a general equation that takes into account the transient condition. k x and k y are hydraulic permeability coefficients in x and y directions, respectively. h(x ⋅ y ⋅ t) denote hydraulic head. Also, w and g indicate the density of water and gravitational acceleration, respectively. m w 2 is the slop of water retention curve. There are two types of boundary conditions for seepage problems: boundary with specific head (Dirichlet boundary condition) and impermeable boundary (Neumann boundary conditions) which are expressed as h = h 0 and h∕ n , respectively. Of course, for the groundwater problems, if there is a river in the aquifer domain, Neumann boundary conditions in the form of h∕ n = q 0 , can be considered. The purpose of numerical modeling is finding the unknown variable h(x, y) . Therefore, in the finite element method, the nodal variables are expressed by the following relation: where j denotes the interpolation (shape) function. Using the FEM in solving Eq. (1), the below relation can be resulted: According to employment of the four-node elements, it can be written as: in which m = w g m w 2 and ̇ n ′ is time derivative of nodal hydraulic head.
It should be noted that can be classic Lagrange or Hankel shape functions. In the following section, the development procedure of Hankel functions is stated.

Suggested new shape functions
The J-Bessel functions that have been used in the other researches contain just the properties of the first kind of Bessel functions. Equation (5) expresses the Bessel equation of order u: The solution of Eq. (5) is in the form of the following form: in which A and B are constant and J u and Y u show the first and second kind Bessel of order u , respectively. It should be noted that the first kind of Bessel functions is not capable of describing all properties of a physical phenomenon, while, by combining the first and second kind of it, Hankel RBFs are obtained and by using them, the numerical modeling is improved. In the following section, the enrichment process is stated completely.

Enrichment procedure of spherical Hankel shape functions
In this part, the enrichment procedure of Hankel RBFs is stated for an arbitrary n-node element in direction, which is shown in Fig. 1.
At first, the polynomial terms and functional expansion of redial basis functions are combined: in which n and m show the number of nodes and basis polynomial terms, respectively. Also, the following relation indicates the parameters and functions used in relation (7).
If the nodal points are replaced in Eq. (7), relation (9) can be resulted: It should be noted that the number of unknowns in Eq. (9) is equal to n + m, while there are only n equations, and thus in order to balance the number of knowns and unknowns, additional conditions are essential [for more details see Farmani et al. (2018) and Bahrampour et al. (2018)].
Finally, the set of equations can be resulted as the following form: By using some algebraic manipulations, relation (12) is resulted: in which, If c and d are replaced into Eq. (7), the following relation can be obtained: The final matrix of the shape functions is suggested as below: Consequently, the following relation expresses the intended RBF: Since there is a singularity at imaginary part of h (1) n ( r) = j̃n( r) + iỹn( r) , the term of ( r)̃n +1 is applied for resolving this problem. Thus the limiting values of the RBF can be obtained as: The intended RBF for n-node element shown in Fig. 1 is calculated as below: In the enrichment process of Hankel RBF, the matrices can be expressed as: , which is restricted by deleting the singularity in Eq. (17).
Finally, the shape functions for an n-node element are written as relation (22): where F, G and E are expressed as: that j, o ∈ A and the set A is expressed as A = {1, 2, ..., n} (Farmani et al. 2018).
It should be noted that the spherical Hankel shape functions have some properties such as Kronecker delta property, Partition of unity, Linear independence property and Infinite piecewise continuity which are stated completely in Farmani et al. (2018) and Bahrampour et al. (2018).
The above process can be easily extended to 2-D elements. In this paper, four-node quadrilateral elements are used for solving the seepage problems by the classic finite element and new methods.

Validation
In this section, to validate and verify the numerical modeling of leakage problems, Laplace equation ( 2 x 2 + 2 y 2 = 0 ) is solved for the domain shown in Fig. 2. In this problem, both boundary conditions, namely Dirichlet and Neumann conditions, exist. Thus, it can be a suitable representative for the seepage flow problems. For the modeling, the numerical results are obtained with high number of elements and then 8 elements for both methods (Fig. 3). Then, these results are compared with the analytical solution [Eq. (26)] and presented in Table 1.
According to the results in Table 1, it can be justified that firstly, Hankel shape functions with 8 elements present

Applicability
In order to demonstrate the applicability and robustness of the suggested formulation in the seepage flows, various numerical models are carried out. Then the present model results with coarse meshes are compared with the traditional FEM with the same and refined meshes. The results are presented in the framework of tables and graphs. It should be noted that, for the results in the tables, the relative error L 2 relative error norm and the computational time are presented.

First problem: groundwater flow under a coffer dam
The first problem is the steady state groundwater flow beneath a coffer dam. The computational domain with boundary conditions is shown in Fig. 4 (Reddy 2004). The Laplace equation is solved by the present model and also the classic method. For this purpose, a coarse mesh (10 elements) and then a mesh with high number of elements are chosen (Fig. 5). The results obtained from both methods are given in Table 2. According to the Table 2, it is inferred that with fewer number of elements in the present formulation, better accuracy than the classic FEM is achieved and the calculated errors are also less. Also, the contours of hydraulic head are indicated in Fig. 6 for this problem with 10-element mesh.

Second problem: flow through the foundation of a dam
In this case, the steady state seepage flow of soil foundation of a dam is studied. It is assumed that the soil under the dam is isotropic and heterogeneous (Reddy 2004). Figure 7 displays the geometry and boundary conditions for this problem. The governing equation is solved by both methods with refined meshes and then the meshes with low number of elements. The grids are shown in Fig. 8. To compare the proposed method results with the classic finite element method, the results are indicated in Table 3. The errors are also computed based on the numerical solution with high degrees of freedom.
As it can be observed in Table 3, the proposed formulation has better agreement with the results by high elements compared with the classic method with the same number of elements. For this case, Fig. 9 shows the contours of hydraulic head for the mesh with 12 elements.

Third problem: groundwater flow in an aquifer with river boundary, pumping and recharges wells
In this part, a homogeneous and anisotropic aquifer is investigated. Equation (1) with steady state condition is solved using the proposed and classical methods. For this problem, the permeability coefficients are considered k x = 1.75 m/ day and k y = 1 m/day in the x and y-direction, respectively. Furthermore, it is assumed that a river passes through the aquifer with infiltrating rate of q 0 = 0.5 m 3 /day/m. Two pumping wells and one recharges well with rates of Q 1 = 900 m 3 /day, Q 2 = − 220 m 3 /day and Q 3 = − 450 m 3 /day are also considered (Reddy 2004). The geometry, boundary conditions, the locations of river and wells are illustrated in Fig. 10. The grids with 100 × 50 and 4 × 4 meshes are used which are depicted in Fig. 11. Also, Figs. 12 and 13 show the results obtained from the both methods at y = 0 and y = According to Figs. 11 and 12, it can be found that using spherical Hankel shape functions with 16 elements leads to having better results than Lagrange shape functions with the same mesh.

Fourth problem: groundwater flow in anisotropic heterogeneous region with river boundary
A domain with different permeability in two directions is considered in this section. Then, the steady state   groundwater flow is studied for the intended aquifer. In the top boundary, the hydraulic head is varied as linear and in the function of x. There is a river at the left boundary with rate of q 0 = 0.24 m 3 /day/m. Other boundaries are impermeable (Reddy 2004). The computational domain and boundary conditions are also shown in Fig. 15. Also, the meshes used are illustrated in Fig. 16. The numerical results are also exhibited in Fig. 17 and Table 4. The contours of hydraulic head are presented in Fig. 18.
It can be concluded from Fig. 17 and Table 4 that the suggested formulation improves the accuracy in the numerical modeling of the mentioned problem. The computational time of Hankel functions is marginally more than Lagrange functions with the same mesh for all problems. However, this time is less than Lagrange functions with high number of elements. Also, the improved accuracy using Hankel functions is significant.

Fifth problem: unsaturated-saturated soil with transient condition
A time-dependent seepage problem is investigated in this section. The studied soil consists of two unsaturated and saturated regions. In order to model the unsaturated part, the permeability coefficient (or hydraulic conductivity coefficient) is assumed to be a function of pore-water pressure. The permeability function for different soils is shown in Fig. 19.
The governing equation for seepage in unsaturated-saturated soil with transient condition is expressed as the following relation:  It should be noted that the soil is assumed sandy for this problem. There is the anisotropic condition in the unsaturated part where the angle between the major permeability coefficient and x-axis is considered = 30 • . In Eq. (27),   (Lam et al, 1987) k xx = k 1 cos 2 + k 2 sin 2 , k yy = k 1 sin 2 + k 2 cos 2 a n d k xy = k yx = (k 1 − k 2 ) sin cos . That, k 1 and k 2 are the major and minor permeability coefficient, respectively.
It can be mentioned that matrix C is defined as = k xx k xy k yx k yy for this problem.
In order to approximate the hydraulic head at the different time step, the central finite difference approximation is used as the following relation: where Δt is the time step that is considered equal to 0.5 min.
It should be noted that the transient seepage equation in the unsaturated region is nonlinear. Because the permeability coefficient is a function of pore-water pressure, an iterative scheme is needed to obtain the correct nodal hydraulic heads. First, an initial permeability is assumed and nodal heads are calculated. Then, the water-pore pressure is obtained by (u w ) t = h n � t − z n � w g that z n ′ indicates the elevation at the nodes of the elements. Therefore, using the pressure and permeability function, a modified permeability value is considered to calculate the new nodal heads. This process is repeated until the difference of the head and permeability for two iterations are smaller than 0.0001.
To model this numerical example, the geometry with boundary conditions is considered according to illustrated conditions in Fig. 21. The problem discussed here is a transient groundwater flow in two unsaturated-saturated regions. The initial depth of water in pond is assumed 1 m. Also, a lining with permeability coefficient equal to 1 × 10 −6 m/s is considered on the bottom of pond. At the left and right boundaries, there are hydrostatic conditions.
In order to investigate the applicability of this problem, the present model results are compared with the results of Lam et al. (1987). A grid with 20 elements (Fig. 22) is studied for the numerical modeling. Tables 5, 6, 7 and 8 present the nodal hydraulic head at different times using finite element method based on both Lagrange and Hankel functions. The negative sign for some values in the tables represents the suction pressure. Also, the contours of hydraulic head at the points when the time is equal to 1316 min are shown in Fig. 23 using 20 elements based on the proposed functions.
It is observed that Hankel functions present the results with more precision than Lagrange functions in comparison with numerical results of Lam et al. (1987). Although the computational cost of Hankel functions is slightly more, the increased accuracy is significant.

Conclusion
Leakage flow in soil was studied numerically by suggesting new shape functions based on the finite element method. In the proposed formulation, the new functions   L 2 relative error norm = 6.01% L 2 relative error norm = 2.62% Table 6 Comparison of the results obtained from Hankel shape function, Lagrange shape function and numerical solution of Lam et al. (1987) at t = 66.5 min for the fifth problem L 2 relative error norm = 9.08% L 2 relative error norm = − 5.88% were derived by combining polynomials and spherical Hankel radial basis functions, and therefore, they contain the advantages of both functions. These advantages also provide robustness and strength for the offered approach. In addition, these properties improve the solution accuracy with fewer degrees of freedom. In order to validate the present method, one boundary value problem with different boundary conditions was studied and then the results were compared with analytical solution. Also, five seepage numerical problems for finding hydraulic head were investigated, and through them, the applicability and reliability of the proposed method were examined. The results indicated that the present model with the low number of elements is more efficient than the classic FEM with the same mesh.
Funding The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript. The authors received no specific funding for this work.

Availability of data and material Not applicable.
Code availability Not applicable.

Declarations
Conflict of interest This article has been submitted to the journal with the full authors consent and agreement. All authors have participated in (a) analysis and interpretation of the data, (b) drafting the article or revising it critically for important intellectual content and (c) approval of the final version. This manuscript has not been submitted to nor is under review at, another journal or other publishing venue. On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.  Lam et al. (1987) at t = 1316 min for the fifth problem Node Lagrange shape function (20 elements) (relative error%) Hankel shape function (20 elements) (relative error%)