A Solution of Richards’ Equation by Generalized Finite Differences for Stationary Flow in a Dam

: The accurate description of the ﬂow of water in porous media is of the greatest importance due to its numerous applications in several areas (groundwater, soil mechanics, etc.). The nonlinear Richards equation is often used as the governing equation that describes this phenomenon and a large number of research studies aimed to solve it numerically. However, due to the nonlinearity of the constitutive expressions for permeability, it remains a challenging modeling problem. In this paper, the stationary form of Richards’ equation used in saturated soils is solved by two numerical methods: generalized ﬁnite differences, an emerging method that has been successfully applied to the transient case, and a ﬁnite element method, for benchmarking. The nonlinearity of the solution in both cases is handled using a Newtonian iteration. The comparative results show that a generalized ﬁnite difference iteration yields satisfactory results in a standard test problem with a singularity at the boundary.


Introduction
Water flow is one of the most relevant phenomena in geotechnical engineering because it modifies the mechanical properties of the soil. In the field of civil engineering, there are many structures where the soil is used as a construction material; among the most important ones are dams and embankments. They are subject to environmental effects since the water enters into structures that are subject to the environment and eventually modify its mechanical behavior. Among the catastrophic events that can be caused by infiltration processes are sinkhole collapses, Ref. [1], and landslides [2,3].
To describe the dynamic mentioned above, Richards' equation is adopted [4]. It is a degenerate elliptic-parabolic equation and is considered an adequate formulation of the problems since it takes into account the unsaturated condition of the soil in a realistic way. Although, an accurate numerical formulation requires an adequate discretization of the strong nonlinear relation between the pressure head, the water content, and the hydraulic conductivity. This nonlinearity must be included both in the temporal and spatial discretizations.
The numerical problem must be addressed from two points of view: the discretization of the differential equation and the solution of the discrete system. This paper aims to address the discretization of Richards' equation. Due to its importance and its wide range of applications, several authors have solved it using different techniques, such as the Finite Element Method (FEM), the Finite Volume Method (FV) and the Finite Differences Method (FD).Extensive reviews on the subject have been conducted, among other authors, by Zha et al. [5], Farthing and Ogden [6], and List and Radu [7]. Nevertheless, despite the intensive research conducted in recent decades, there is still a strong need for robust numerical schemes for multiphase flow in porous media.
FEM has been widely used and is a standard reference method for benchmarking due to many studies where successful results were obtained in one-, two-, and threedimensional problems [8][9][10][11][12][13][14][15]. However, in some transient problems, the numerical results for infiltration were found to exhibit strong oscillations, a phenomenon that has been associated with the mass matrix in the formulation [8].
The finite volume has been used to solve Richards' equation in one [15][16][17][18][19] and two spatial dimensions [17,20]. One of the advantageous features of the finite volume method is that water balance is taken into account in the scheme at every time step since it is calculated at the borders of a control volume as a part of the algorithm, implementing the mass-conservative iterations straightforwardly [17,18]. Another advantage is that the method has applicability in unstructured meshes, which endows it with flexibility [17,21].
In the case of finite differences, most of the papers analyze the one-dimensional problem [8,[22][23][24][25] to assess the accuracy, stability, and convergence of the method [21]. One of its main limitations is the lack of schemes in irregular domains because the stencils developed for structured meshes could not be used in 2D complex domains with unstructured grids. Nevertheless, some interesting adaptations were proposed, such as that of Hsieh [26], which was designed for block-rectangular domains.
Differences in block domains and symmetric nonrectangular domains prepared the way for more general schemes. Nowadays, there are new formulations of finite differences that can be applied to complex domains with structured and unstructured meshes. In recent times, the Generalized Finite Differences Method (GFDM) [27,28], which is an extended case of the classic finite differences, has shown to be especially useful for irregular boundaries without losing the simplicity of the traditional method. The method has shown the capabilities to numerically solve the transient differential equations which describe water flow in complex boundaries with structured grids [29][30][31][32]. The simplicity behind the formulation of the GFDM is one of its main advantages. In this case, the idea is to obtain weights for each node, as a function of the relative distance of the surrounding nodes to a central point using the Taylor expansion. This yields a method that can even be applied as a meshless one, because a specific grid structure is not required to assign the weights. For instance, in Milewski's work [33], an unstructured, arbitrarily distributed cloud of nodes is used, as well as an ad hoc node selection, to produce accurate results. Certainly, generalized finite differences can also be applied to structured meshes and irregular domains [27], which has an advantage over the meshless methods that the grid structure for the calculation of the coefficients uses.
Due to this reason, the GFDM has also been used for solving parabolic and elliptic partial differential equations successfully. In the present work, the GFDM is tested against FEM to solve Richards' equation in stationary formulation for a dam case with boundary singularities [34]. This form of the Richards' equation is an elliptical PDE with highly nonlinear constitutive relations. The hydraulic conductivity k(θ) and suction s(θ) functions have nearly zero values that can produce a degenerate form of the equation [6], which makes the problem especially complex due to the dependence of the hydraulic permeability on the solution, so a Newton-Raphson (NR) iteration needs to be applied to solve the nonlinear discrete problem. This paper is organized as follows: Section 1 presents the introduction to the problem. The materials and methods are discussed in Section 2 and the results and discussion are addressed in Section 3. Finally, the conclusions are presented in Section 4.

Richards' Equation
Richards' equation [4] for the case of the stationary condition in two dimensions can be written as [35] ∂ ∂x k x ∂h w ∂x where: h w = y + u w γ w is the total hydraulic head; y is the gravitational head; u w is the pore water pressure; γ w is the unit weight of water; k x (s) and k y (s) are the water permeability coefficients along the xand y-axis; s = u a − u w is the matrix suction relative to the reference air pressure u a , which was set as equal to zero in this paper.
In this formulation, the cornerstone to produce realistic results is an adequate definition of the permeability coefficients. Under isotropy considerations, the same coefficient k can be used in both coordinate directions. The van Genuchten-Mualem equation [36], which can be written as is often used for that purpose. Here: θ n = θ w −θ r θ s −θ r is the normalized volumetric water content; θ w is the actual volumetric water content; θ r is the residual volumetric water content; θ s is the saturated volumetric water content; n, m are fitting parameters, related as m = 1 − 1 n . On the other hand, the normalized volumetric water content θ n can be represented by the following equation proposed by van Genuchten [37] θ n (s) = 1 ((αs) n + 1) m where m and n are the fitting parameters defined in Equation (2).

Finite Difference Formulation
Under smoothness assumptions, Richards' Equation (1) can be written in a differential form as follows For this work, bearing in mind its successful application to solve an unsteady Richards equation [28], the formulation of generalized finite differences is used to approximate the solution of the Equation (4).
The procedure to approximate differential operators using generalized finite differences for nonrectangular regions considers a linear combination of the values of the function h w at a finite set of nodes in such a way that the local consistency condition is fulfilled. Thus, the truncation error at τ(p 0 ) is used to define a second-order approximation to the linear operator at the point p 0 . Here, the values Γ i , i = 1, ..., q, are the weights in the generalized finite difference scheme, and the functions A(x, y), B(x, y), C(x, y), D(x, y) and E(x, y) are coordinate dependent functions which define the second-order operator L in Equation (6). The approximation sought is given by the second term on the right-hand side of the Equation (5), which is a linear combination of h w , evaluated at p 0 , p 1 , ..., p k , while the first term is the evaluation of the operator L in the central node p 0 . Consistency implies that the coefficients Γ l can be obtained from the Taylor expansion of Equation (5) at p 0 by imposing the following conditions where ∆x l and ∆y l denote the x and y components of p l − p 0 (see [28]). System (7) defines the generalized finite differences. The classical stencils and schemes are characterized by weights Γ i , which are a subset of its solutions; however, by choosing a large enough number of nodes so that system has a full row rank, it is possible to define complete manifolds where the weights can be chosen. Naturally, the full rank condition can be weakened in the presence of symmetries of the nodes, as happens in the rectangular stencils.
From the numerical point of view, since Expression (7) is graded, its solutions must be thoroughly calculated. Given that it is a block system, the first equation can be separated to solve the system [28].
independently. It is important to note that the latter equation must have a full row rank, but at least 5 independent columns to avoid overdetermination and in consequence lack of consistency. Each column corresponds to each one of the points p 1 , ..., p q , so these nodes must be chosen in a balanced way around p 0 . Several authors use radial weighted functions on clouds of points as meshless methods for that purpose. Very remarkable results have been obtained with this approach, but the accuracy of the results is strongly dependent on an adequate ad hoc weight function selection [31,32,34,38]. On the other hand, the node selection proposed in this paper is one of its distinctive features, since it makes use of the known edges and nodes of a grid as a simple criterion to select the neighbors.
Once the system (9) has been solved, the remaining weight Γ 0 is obtained from Equation (8). In a rather general grid, where the positions of the neighbor nodes p i , i = 1, · · · , q relative to p 0 change from one inner node to another, the coefficients Γ l also change, since they depend on its position. However, two facts must be emphasized: • If the grid is regular, for instance, an actual rectangular grid, the matrix coefficient in Equation (9) will be equal for all the inner grid nodes. Additionally, if the right-hand side of that equation is constant, which is the case in several important equations (v.gr. Poisson equation), the coefficients will be the same for all the inner grid nodes.
This implies that in specific modeling applications, the calculation of the coefficients can be very efficient. • Furthermore, if the positions of the neighbors are symmetrical to the central node and the differential operator is self-adjoint, the number of different coefficients required decreases (v.gr. four nodes is the discretization of the Poisson equation). This is precisely the case of the standard differences in rectangular grids. However, the analogous case occurs in other regular structures, such as those used in [39].
Equation (7) is used for approximating the solution of Equation (4) considering that, at the nodes of a grid, the dependence of k x , k y , h w and their partial derivatives on the nodes and their corresponding neighbors, which yields a nonlinear discrete system of the form where H(ĥ w ) is the assembled matrix of coefficients Γ 0 , Γ 1 , ..., Γ k obtained at each inner grid node from the solution of Equation (7),ĥ w is the vector of unknown hydraulic head values and f is the vector of boundary values assembled with the known Dirichlet boundary conditions. The matrix H depends on the hydraulic conductivity, which is a function of the unknown variable h w (through Equations (2) and (3)); thus, the problem becomes nonlinear and an iterative scheme solution must be applied.
Regarding the discrete assembled systems obtained using the proposed discretization, since H shares matrix structures that are common to other numerical methods, namely, finite elements, finite volumes, the finite pointset method, and the radial basis function method, and the solution of such nonlinear systems and their linearizations has been thouroughly studied [7], the focus on the papers remains on the new discretization. Thus, the solution of Equation (10) can be approximated from the following equation: the residual of the equation, is used to correct iteratively the initial assumed value ofĥ n w aŝ h n+1 w =ĥ n w + α k n ∆ĥ n w (13) where α k n is the adaptive Armijo's step length and J is the Jacobian is defined as J(ĥ n w ) = H(ĥ n w ) + H (ĥ n w ). The last term in the Jacobian corresponds to nonsymmetric terms; if it is neglected, the conventional iterative method is obtained. The solution will be reached until the value of ψ(h n w ) is less than the given tolerance [40]. The adaptive step length α k n was initially set to 0.7. Whether a sufficient decrease in the residual is not obtained, the step ∆ĥ n w is multiplied by α k n . This Modified Newton-Raphson Method (MNR) was used in the numerical tests of the next section.

Problem Domain and Boundary Conditions
The domain for the case study is a typical cross-section of a dam of 12 m height; it is shown in Figure 1 along with its boundary conditions. The dam is built of only one kind of soil; the water level is located at 11 m on the upstream side (h w = 11). The dam bed is considered impervious except for a drain of 8 m of length placed on the downstream side, where the pore water pressure is set as equal to zero. The downstream boundary condition is also set to zero, with the intensity direction of the normal flux defined by the direction cosines x , y . The same case is presented in the upper part of the upstream boundary above the water level and at the top of the dam. It must be noted that in the case where the water level exceeds the downstream boundary height, the water pressure equals the atmospheric pressure, and then the boundary conditions must be corrected accordingly.

Grid Independence
In order to assess the grid independence of the generalized finite difference discretization as defined by Equation (7), a sequence of four grids for the problem domain was generated. The first one is a coarse mesh with 992 triangles; it was generated using Delaunay triangulation for a set of 553 inner grids nodes obtained from a given initial node distribution. No smoothing process was applied. The next three grids were obtained by a process of refining, on which a splitting procedure was selected; in other words, each triangle was subdivided into four triangles by splitting each triangle side by his midpoint [41].
Following this procedure, grids with 3968, 15,872, and 63,488 triangles were generated (see Figure 2). Using this grids, Equation (10) was solved using the parameters described in Section 2.3.3. The phreatic line for the finest mesh (d) is shown in Figure 3. A more detailed view of the phreatic lines is displayed in Figure 4. The fast approximation of the consecutive phreatic lines to the one corresponding to the finest mesh can be seen, even though the meshes were not smoothed or preprocessed.
In Table 1, a summary of the root mean square errors (RMSE) between the different phreatic levels obtained is shown. The values in the column labeled as A correspond to RMSE errors between the phreatic curve for a mesh and the phreatic level for its refinement. In the third column, labeled as B, the phreatic level obtained for Mesh 3 (see Figure 2) was taken as the reference to calculate the RMSE errors. The values in the table are consistent with the fact that the phreatic levels in the coarser grids tend to the last one, thus assessing grid independence for the proposed discretization.

Stationary Flow in a Dam
Regarding the problem of interest, the modeling stationary flow in a dam was solved with two numerical methods: the generalized finite differences method described in the previous section and the finite element method using linear elements as a reference. Both of them were codified in Julia 1.5. Having the shape of the phreatic levels obtained in the previous subsection, two new different meshes were used in the numerical experiments: a rather coarse one (Figure 5a) and another with refinement in the zones where abrupt changes in ∂h w /∂x (Figure 5b) are expected. The first mesh has 778 nodes and 4221 triangular elements with an approximate side length of 0.75 m. The refined mesh has 2559 nodes and 14,661 elements and was obtained from the coarse one by refining around the boundary sides where Dirichlet and Neumann boundary conditions meet. The soil parameters used in the tests are shown in Table 2. The finite element method solution in the refined mesh is shown in Figure 6, where the line at u w = 0 is the phreatic line. For the next analysis, it will be taken as a reference for the comparison of the different solutions presented.  As mentioned, the Newton-Rapson method was used to solve the discretized nonlinear Richards equations. To evaluate the performance of the GFDM compared to that of the FEM method, the number of iterations required to achieve convergence in the solution in both methods with the two meshes is presented. The algorithm uses two-stop criteria: the relative tolerance error (set to 1.0 × 10 −8 ) as the convergence criterion and a maximum number of iterations (100) to avoid endless or slow loops; the first occurrence breaks the cycle.
The root mean square error (RMSE) along the reference phreatic line was calculated with the Equation (15), where: h f d is the finite differences approximation to the hydraulic head; h f e is the finite element approximation to the hydraulic head; N is the number of measurements along the reference line.

Results and Discussion
A zoomed in view of the solution near the drain zone can be seen in Figure 7; the high slope can be easily noticed. In a two-dimensional view, some additional remarkable details can be noted. In Figure 8a, the phreatic line obtained with both methods in the refined mesh is presented. The value of the area enclosed by the solutions is equal to 2.2396 m 2 in the coarse mesh and the value of the RMSE is equal to 0.1152 m. However, if one looks at Figure 8a in detail, near the base there is a separation of the two calculated lines. The one in blue (FEM method) has a slight shift to the right, an effect that occurs due to the mesh, whereas the red line (GFDM method) is completely vertical. This suggests that the GFDM solution is smoother than the FEM solution near the base. To highlight the effect of the mesh mentioned above, in Figure 8b, the solution with the coarse mesh is shown. The separation of the two lines increases and this causes the RMSE error and enclosed areas to increment to 0.1969 and 2.2414 m 2 , respectively, due to abrupt movements of the blue line (FEM results) caused by the singularity at the boundary. As a reference, the maximum distance between the two lines is printed in the figures (max error), and it can be noted that, for coarse mesh, its value is 0.7268 m, and for the refined one, it is equal to 0.2772 m; these numbers confirm the observation mentioned above.
The analysis of the performance of the Newtonian iteration (13) provides some extra insight into the solution process. The semilogarithmic Figures 9 and 10 present the convergence history for solving the problem on both grids; Figure 9 shows results on the coarse grid. Clearly, on the coarse grid, even when in both cases convergence is achieved, the FEM iteration shows a better performance. On the other hand, on the refined mesh, adapting the step size length under the Armijo rule was necessary to achieve convergence, the GFDM being notably faster in the latter case.   Additionally, Figures 11 and 12 show the evolution of the Jacobian condition number on the iterations. In both cases, the condition numbers are relatively large, although the Jacobian on the FEM iteration is in a worse condition than in the case of the GFDM iteration even when the same Newtonian procedure is applied. Nevertheless, when the Armijo rule is applied, convergence was obtained in all the tests.  In Table 3, a summary of the convergence history for each test is shown. According to the graphs, although in most of the tests there is convergence, the error sequences behave quite differently. In the coarse grid, numbered as 1, the FEM iteration is faster, but the estimated convergence order in the table in both cases is essentially the same. It must be pointed out that to produce a fair comparison, the reported value of the estimated empirical order corresponds to the average of the last 40 iterations. Thus, the oscillations shown in the FEM iterations, which do not avoid convergence, are smoothed for the analysis. However, despite the oscillations, there is a clear graphical trend that shows a better performance. The value of the Least Squares Slope (LS Slope), taken in the same last 40 iterations, corroborates this fact.
On the other hand, when the grid is refined, convergence cannot be achieved without the Armijo adapting step size length strategy. It can be seen that, in both methods, convergence is accelerated, with the value α = 0.75 being the one the produces the best results. The convergence order is notably larger, but now the LS Slopes show that the GFDM iterations produce a faster iteration.

Conclusions
As the main conclusion, it is possible to state that the GFDM is a useful numerical technique to solve several differential equations of interest-in the present case, the steady Richards' equation was used, which is a nonlinear expression. As a numerical method, GFDM keeps most of the simplicity of the standard formulation on rectangular grids, but can also be applied to nonstructured meshes. In the proposed tests, carried out under the same conditions and taking advantage of a suitable Newtonian iteration and an adaptive step length strategy, the GFDM was capable of achieving convergence faster than FEM in the refined mesh. Additionally, the calculated GFD solution shows a much smoother convergence history than the FEM solution, even in the presence of a singularity at the boundary, a characteristic that may be useful to improve the iterative solution of the discrete assembled systems on infiltration modeling and similar problems.
Considering the diverse applications addressed in the introduction, as well as the results shown in this paper, GFDM stands as a reliable alternative numerical method to solve some differential equations used in applied engineering. Nevertheless, it is convenient to test this numerical technique with other problems to gain a deeper understanding of its capabilities.