Next Article in Journal
Application of Response Surface Methodology on Brewery Wastewater Treatment Using Chitosan as a Coagulant
Next Article in Special Issue
Integrated Management and Environmental Impact Assessment of Sustainable Groundwater-Dependent Development in Toshka District, Egypt
Previous Article in Journal
Evaluation of Social Vulnerability to Flood Hazard in Basilicata Region (Southern Italy)
Previous Article in Special Issue
Analytical Method for Groundwater Seepage through and Beneath a Fully Penetrating Cut-off Wall Considering Effects of Wall Permeability and Thickness
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Deforming Mixed-Hybrid Finite Element Model for Robust Groundwater Flow Simulation in 3D Unconfined Aquifers with Unstructured Layered Grids

1
Independent Researcher, 45000 Orléans, France
2
Regional Water Centre of Maghreb, LIMEN, Ecole Mohammadia d’ingénieurs, Université Mohammed V de Rabat, Avenue Ibn Sina B.P. 765, Rabat 10090, Morocco
*
Authors to whom correspondence should be addressed.
Water 2023, 15(6), 1177; https://doi.org/10.3390/w15061177
Submission received: 29 January 2023 / Revised: 7 March 2023 / Accepted: 16 March 2023 / Published: 18 March 2023
(This article belongs to the Special Issue Flow and Transport Processes in Groundwater Systems)

Abstract

:
Determining the water table shape and position in unconfined aquifers is fundamental to many groundwater flow assessment studies. The commonly used industry-standard fixed mesh models, contrary to popular belief, do not provide an accurate description of the phreatic surface. When using such models, the water table position is post-processed from the simulated groundwater heads, leading to an approximation error. This error becomes larger for coarse vertical grids. This paper introduces a novel moving mesh technique to simulate the groundwater table in three-dimensional unconfined aquifers under steady-state or transient conditions. We adopt the face-based mixed-hybrid finite element discretization approach in space, leading to a more accurate approximation of the specific discharge field. The model uses an adaptive unstructured but layered mesh which is iteratively adjusted until its top fits the phreatic surface. The developed algorithm accounts for a linearized form of the kinematic boundary condition prescribed on the moving boundary and also supports usual boundary conditions as well. The model was compared to the existing analytical, fixed mesh, and previously published solutions. The obtained results show that the developed model is superior in terms of its numerical stability, convergence behavior, and accuracy. Furthermore, the simulated phreatic surface is free from a cellwise interpolation error and independent of the vertical grid size as used in fixed mesh methods. We also found that the robustness of the moving mesh method cannot be surpassed by a fixed mesh alternative. The model’s efficiency is supported by an almost quadratic rate of convergence of the outer iteration loop.

1. Introduction

Nearly half of the world population is supplied by groundwater aquifers [1], while the percentage of the U.S. population relying exclusively on groundwater is estimated to be between 40% and 50% [2]. Groundwater flow modeling is a widespread technique to sustainably manage these resources [3]. Indeed, many groundwater modeling tools are readily available [4,5]. The numerical modeling of groundwater flow in unconfined aquifers is much more involved than in confined aquifers. This is because the governing equation (i.e., Richard’s equation) is highly nonlinear and is subject to nonlinear boundary conditions as well. The nonlinearity in Richard’s equation is related to the dependence of the relative permeability and the water retention in the unsaturated zone on the pressure head [6,7]. Moreover, fully saturated models are typically associated with boundary conditions such as constant head, pumping/injection flow rates, and leakage flux through a semi-confining bed. These are essentially linear and do not pose additional challenges to standard numerical solution techniques. This is not the case for unconfined flow models where some boundary conditions are nonlinear and therefore are unknowns a priori; thus, they are an integral part of the numerical solution. A typical example is the free and moving water table boundary. Another complication in modeling unconfined groundwater flow is to locate the position of seepage faces when the water table reaches the land surface. The seepage boundary faces are not known prior to the numerical solution and are, therefore, not fixed as typically done with other boundary conditions. Additionally, wells pumping groundwater from unconfined aquifers might become dry as the water table drops below the well screens, and natural or artificial drains may stop draining groundwater as the levels drop below predefined elevations. These are typical examples that necessitate additional bookkeeping procedures during the nonlinear iterative solution process when solving for an unconfined groundwater flow model.
Natural groundwater flow occurs mostly in highly heterogeneous and anisotropic materials, thus increasing the nonlinearity of the resulting discrete algebraic systems of the equations. Hence, robust and advanced numerical methods are needed in this context.
For unconfined aquifers, a decision has to be, very often, made for the choice of the groundwater modeling approach. The most rigorous approach selects a numerical discretization of Richard’s equation which explicitly accounts for the processes in the unsaturated zone (Figure 1B). These are, for instance, the occurrence of a capillary barrier above the saturated zone where the soil moisture is variably distributed, and gravity routing of infiltrated surface water to confined storage. The second choice is to use an approach that simplifies the underlying physics in the porous medium. Hence, the water table is viewed as a sharp interface whose position in space and time is sought, and the flow in the unsaturated zone is neglected. The last approach, although not free from limitations, is by far the most used in practice when dealing with regional-scale modeling studies. This is because the more detailed and accurate approach requires much more field data which are frequently not available on a large scale. Hence, groundwater modeling of the unsaturated zone is, very often, performed on a local scale of a site where experimental data are available. For some projects, the two approaches might advantageously be used concurrently to understand the local behavior in the unsaturated zone, thus enabling a better parametrization of an unconfined flow model built with the simplified approach. This paper deals exclusively with the second approach; modeling the unsaturated zone will not be discussed further.
Modeling unconfined groundwater flow subject to free and moving boundaries is performed using two broad classes of methods. The first category assumes a fixed mesh as illustrated in Figure 1C. Within this class, an algorithm should determine the spatial and temporal distribution of three types of cells: wet, dry, and partially wet cells (i.e., crossed by the water table). The second category considers a moving mesh that adapts to the dynamics of the water table boundary (Figure 1D). Hence, unlike for a fixed mesh method, the dry zone (i.e., approximate unsaturated zone) is not part of the computational domain. The main advantage of an adaptive mesh technique is to solve a saturated groundwater flow equation, hence a linear system of algebraic equations, during each iteration of the mesh adaption process. Meanwhile, more elaborate algorithms are needed for mesh adaption and the interpolation of hydraulic properties between moving cells.
The fixed mesh approach is probably the most widespread owing to the availability of many computer implementations such as MODFLOW developed by the USGS [4,8,9]. Legacy versions of MODFLOW [8] were criticized for an unphysical representation of unconfined groundwater flow dynamics, leading to a numerical instability for highly nonlinear problems [10,11]. This was due, in part, to the heuristics upon which the decision to rewet a previously dry cell are based. Second, Picard iteration for some highly nonlinear problems may converge at a very slow rate or not at all. The latest releases include, however, a new Newton–Raphson-based formulation that substantially improves the computational stability [4,9]. Another known difficulty within the fixed mesh approach derives from the way to artificially hide the flow patterns in the dry cells. Because this is not possible as the mesh is fixed, different approaches have been developed to tackle this problem. Desai [12] developed a two-dimensional finite element model that linearizes the relative permeability function and recast it as a function of hydraulic head. The same approach was used by [13,14] in the framework of three-dimensional finite element models. Other authors [15] used a sharp representation of this function; that is, the hydraulic conductivity in all dry cells is multiplied by a small residual factor (i.e., ~ 10−2–10−3). New developments presented the numerical solution of 3D unconfined seepage problems with the smoothed finite element method borrowed from solid mechanics [16]. A limitation of the fixed mesh approach is that the exact position and shape of the water table could only be determined by post-processing computed groundwater heads. Hence, it is expected that a moving mesh technique will always lead to a more accurate position and a smooth shape of the phreatic surface when compared with a fixed mesh approach using the same resolution.
The second class of methods solve a series of linear equations on a deforming mesh. The mesh is adapted iteratively to fit the water table position. Here, again, many sub-approaches have been developed that could be categorized as rigorous or simplified. As described in standard textbooks [6], the phreatic surface does not only fulfill the zero-pressure condition, but it is a kinematic boundary condition too. In other words, it is an interface through which the effects of unconfined storage, net recharge, and nonlinear effects related to the changing shape of the surface geometry are balanced. To the best of our knowledge, the only published numerical model taking into account both conditions was presented by [10]. As we will demonstrate later, even under steady-state conditions, the kinematic boundary condition cannot be neglected.
Previous works in the literature using a moving mesh approach were developed using standard finite difference (FDM), finite volume (FVM), and finite element (FEM) methods. This is because, within a FDM framework, the grid should remain orthogonal and fixed, so the groundwater flow equation is rewritten in a rectangular computational domain using a coordinates transformation technique. An example of this is the so-called boundary-fitted coordinates (BFC) method [17,18]. The same method was also used by [10] for an application to the Delaware Basin of south-eastern New Mexico. This approach leads to a transformed groundwater flow equation which is of the advection-diffusion type, leading to more computational storage as the resulting system matrix is unsymmetric. Additionally, the numerical solution of this equation type is not only more involved than the standard groundwater flow equation, but it sounds unnatural to many groundwater modelers. Another method’s disadvantage is using ghost cells close to the water table boundary because it is not possible to discretize the kinematic condition exactly on the phreatic surface.
Owing to the limitations of the FDM method, the standard FEM is a more attractive alternative. An FEM-based moving mesh approach has been introduced by [19,20] for groundwater applications for the first time. A two-dimensional cross-section model using triangular finite elements was developed for heterogeneous alluvial aquifers and applied to simulating irrigation and drainage in California’s western San Joaquin Valley [21]. In this model, the stratigraphic domain was cut from the top along the new position of the water table curve, and the updated computational domain was re-meshed prior to finding a new groundwater head distribution. In [22], a nonlinear programming minimization algorithm was coupled with an adaptive finite element procedure to determine the location of the phreatic line in earth dams. This model solves a steady-state two-dimensional seepage problem using bilinear triangular finite elements. A parallel three-dimensional adaptive mesh finite element model was presented for a steady-state unconfined flow in subsurface aquifers [23]. An initially coarse unstructured mesh is iteratively adapted along the horizontal direction to refine the mesh in areas of a high hydraulic gradient, but also along the vertical direction to fit the upper slice of the mesh to the phreatic surface.
In [24], a moving curvilinear finite volume grid was used for solving steady-state seepage flow in two-dimensional homogeneous and isotropic earth dams subject to nonlinear boundary conditions. The developed algorithm was favorably validated against the experimental results.
As described earlier, solving unconfined groundwater flow with a moving mesh technique using a finite difference method has some limitations. These problems do not exist when using a conforming (i.e., nodal-based) finite element technique such as the Galerkin weighted residual approach [25,26]. However, while the groundwater head and the water table shape obtained by the conforming FEM using a moving mesh method might be accurate, it is well known that this is not the case for the specific discharge field. Moreover, the conforming FEM is not exactly mass conservative, which is problematic for nodes sharing elements with a high contrast of hydraulic conductivities. Post-processing techniques were introduced to derive continuous nodal-specific discharge fields, but it was reported that this procedure has a high computational cost [27,28,29,30]. Mixed finite element methods were introduced to simultaneously approximate a scalar variable and its first-order derivative fields [31]. Promising results highlighting a higher accuracy of this approximation technique when compared with alternative techniques were obtained for flow problems [32]. The mixed finite element method leads, however, to an indefinite system of algebraic equations which is difficult to solve with direct or iterative methods. A hybridization technique avoids this problem by reformulating the mixed problem with primary unknowns on the faces of the mesh and using local algebraic relationships to recover the cell-centered heads and the normal specific discharge components [33]. This gives rise to a method formally known as the mixed-hybrid finite element method (MHFEM). Previous works established the outstanding behavior of this method to tackle groundwater flow problems in highly heterogeneous and anisotropic formations [28,34,35]. While the MHFEM approximation was used for many classes of groundwater flow problems such as unsaturated flow [36,37,38,39] and two-phase flow [40,41], to the best of the authors’ knowledge, this is the first time that it has been used in combination with a deforming mesh approach to solve an unconfined flow in phreatic aquifers.
The objective of this paper is to introduce a new state-of-the-art three-dimensional MHFEM numerical model for confined–unconfined groundwater flow in complex aquifer systems characterized by irregular geometries, heterogeneous and anisotropic materials, and nonlinear boundary conditions of a practical relevance. These include, in particular, the nonlinear kinematic and seepage face boundary conditions generally ignored or weakly dealt with in previous studies. The model formulation and numerical discretization on layered unstructured grids are presented in detail in the next sections. Several theoretical and realistic examples are provided to demonstrate the model’s accuracy, efficiency, and capability to successfully simulate unconfined groundwater flow with a moving mesh technique.
The manuscript is organized as follows. We start by exposing the theory of unconfined groundwater flow in a deforming mesh, namely, how the free and moving boundary conditions are formulated in this particular case. The next section shows the discrete form of the three-dimensional mixed-hybrid finite element model when a moving mesh approach is used. Notably, the numerical implementation of the associated boundary conditions, which is tricky in some cases, is presented in sufficient detail. Section 4 presents the numerical algorithm that iteratively solves a steady-state or transient unconfined flow problem with a deforming mesh. The next section validates the developed model on a series of available analytical and numerical solutions of an increasing complexity. They demonstrate the model’s accuracy, efficiency, robustness, and applicability to a wide range of hydrogeological problems typically encountered in practice. Finally, a summary and concluding remarks section closes the paper.

2. Theory

When modeling unconfined groundwater flow in three-dimensions using a moving mesh approach, there are two domains of interest. The first is the stratigraphic domain (known also as the physical or geological model) where the hydrogeological units and their hydraulic properties are defined. The second is the moving computational domain where the top surface is adjusted to fit the unknown water table position. The stratigraphic domain is inherently static and is used mainly in the background to map the hydraulic properties into the deforming domain when needed. One advantage of the moving mesh approach is to solve a series of linear problems on the deforming domain using a saturated flow equation. The non-linearity is, therefore, moved from the governing equation to the process of mesh adaption.

2.1. Governing Equations

The governing equation for confined groundwater flow in the porous medium filling the deforming domain is presented as in [6,7]:
S p h t + . q = q s
where S p [L−1] is the specific storage coefficient, h [L] is the hydraulic head, t [T] is the time, q s [T−1] denotes the source/sink strength term positive for sources and negative for sinks, [L−1] is the differential operator, and q [LT−1] is the specific groundwater discharge, following Darcy’s law:
q = K . h
where K [LT−1] is the hydraulic conductivity tensor. In the presented model, we use a full tensor for each soil type in the stratigraphic domain:
K = [ k x x k x y k x z k y x k y y k y z k z x k z y k z z ]
Supporting a full hydraulic conductivity tensor allows for modeling of the groundwater flow in cases where the stratigraphic domain is highly heterogeneous and anisotropic, but also in the frequent situation where there is a mismatch between the scales of the geological model and the computational domain. Using static or dynamic upscaling approaches may produce non-zero off-diagonal equivalent hydraulic conductivities even if the original (i.e., fine scale) hydraulic conductivity tensor is aligned with the principal directions [42].
By inserting Equation (2) in Equation (1), we obtain the governing equation for three-dimensional groundwater flow in the heterogeneous and anisotropic deforming domain:
S p h t = . ( K . h ) + q s
This equation is subject to the initial and boundary conditions of Dirichlet, Neumann, and Cauchy types, as will be further discussed in the next subsection.

2.2. Boundary Conditions

In the framework of a moving mesh approach, some boundary conditions require a special treatment. The conceptual approaches behind such boundary conditions when using the fixed mesh approach are discussed in many references [3,4,8]. Special attention is paid to these boundary conditions because they are part of the solution and, second, they affect the final solution inside the computational domain. The numerical implementation of these boundary conditions is discussed in the following section once the discrete numerical model has been obtained.

2.2.1. Phreatic Surface

The sharp phreatic surface, Γ W T , delimiting the water and air phases, is of special interest here. This is because the shape and location of this interface is unknown a priori. Hence, additional constraints on the free surface are needed. The theory of unconfined groundwater flow [6] (pp. 252–259) stipulates the existence of two distinct boundary conditions. The first states that the pressure at the water table is atmospheric, hence:
h ( x , y , z W T , t ) = z W T ( x , y , t )
where z W T [L] is the water table elevation. Second, the kinematic boundary condition is a statement of the continuity of the mass flux across the phreatic surface. According to [6], it is provided as:
S y h t = ( K . h + N ) . ( h z )
where S y [L−1] is the specific yield and N ( x , y , t ) = N z [LT−1] is the flux associated with net recharge (i.e., inflows such as recharge minus outflows such as evapotranspiration) across the surface. The kinematic condition is nonlinear, owing, first, to the occurrence of the quadratic terms in h and, second, because the conductivity coefficients depend upon the location of the water table inside the stratigraphic domain. It indicates that the flux vector across the water table is not entirely vertical. It is a combination of linear and nonlinear flux terms. The nonlinearity derives from the changing shape of the deforming computational domain volume. A direct numerical discretization of Equation (4) subject to boundary conditions (5) and (6), in particular, increases the model’s nonlinearity and impacts its efficiency. Alternatively, we use a linearized version of Equation (6) borrowed from Dagan [43] that was used also by Neuman [44]:
S y h t = ( K z + N ) h z + N
where K z is the vertical component of a diagonal hydraulic conductivity tensor. The linearized flux in Equation (7) does not only neglect the quadratic terms, but it also neglects the horizontal flux components because a diagonal hydraulic conductivity tensor was adopted to derive this expression. The simplified kinematic boundary condition in Equation (7) stipulates that the head rate of change at the water table equals N / S y , excluding a term which is proportional to the vertical hydraulic gradient h / z .
When a full hydraulic conductivity tensor is used as adopted herein, the development of Equation (6), while still neglecting the quadratic terms, yields:
q z = N + N h z + S y h t
q z = ( K z x h x + K z y h y + K z z h z )
where q z is the vertical component of the specific discharge vector. Equation (8) relates the vertical specific discharge to net recharge, the vertical hydraulic head gradient, and the unconfined storage owing to the water table displacement. This is in sharp contrast to an unconfined flow model relying on a fixed mesh where q z would be, naturally, equal to N . Prescribing the net recharge flux for a moving mesh approach in this way would result in an erroneous estimation of the specific discharge at the interface. This error equals to | N h z + S y h t | , which accounts for contributions from the vertical hydraulic head gradient and the unconfined storage. It vanishes only in the particular case when considering a steady-state regime with a zero net recharge.
Under steady-state conditions, Equation (8) becomes:
q z = N + N h z
enforcing the vertical specific discharge component to scale with the net recharge across the water table interface. Here, again, prescribing the recharge directly on the moving water table would result in an erroneous estimation of the vertical flux component of the specific discharge.
The form of the kinematic condition provided in Equation (8) is suitable for a face-based finite element technique such as the mixed-hybrid finite element method. The numerical implementation of the kinematic condition is described in detail in the next section.

2.2.2. Seepage Face

Another complication arises when modeling unconfined flow with a moving water table approach. Rising groundwater levels may seep through the land surface as often happens in earth dams and trench drains. This interaction between the free-surface and the basin topography is approached using the seepage boundary condition [6]:
h ( x , y , z W T , t ) = z T O P ( x , y )
where z T O P [L] is the land surface elevation. This is essentially a switching Dirichlet boundary condition that applies to a set of unknown locations, adding further nonlinearity to the model. Indeed, a numerical implementation iteratively checks for the occurrence of this condition. When the computed free-surface elevation is above z T O P , Equation (11) replaces Equations (5) and (8) as they apply only when the water table is below the land surface. The kinematic condition provided by Equation (8) is not prescribed at seepage faces. The complete numerical implementation of the seepage face boundary conditions is provided in the next section. The numerical treatment of seepage faces is identical under steady-state and transient conditions.

2.2.3. Recharge

The downward recharge flux, R E ( x , y , t ) [LT−1], is prescribed directly on the moving water table boundary, as presented in Equations (8) and (10), for transient and steady-state regimes, respectively. Under transient conditions, the recharge may take a non-negligible time to reach the saturated zone from the land surface. This is known as the unsaturated zone time lag, t u ( x , y ) , which is site-specific and depends on many parameters, such the water table depth, layered heterogeneity in the unsaturated zone, rainfall intensity, crops water demand, etc. Many approaches have been developed in the literature [45,46,47] to provide first-order estimates of this parameter which can be provided as a model input to delay the recharge flux prescribed at the horizontal coordinates ( x , y ) and time t as R E ( x , y , t t u ( x , y ) ) . Obviously, a delayed recharge flux is not needed under steady-state conditions. Based on observed rainfall and water levels time series in six chalk aquifers, ref. [48] concluded that for shallower water tables, a linear relationship exists between the lag time and the depth to the water table. In [49], simple analytical methods based on a steady-state assumption and rigorous one-dimensional transient simulations were compared to assess the travel time in unsaturated zone profiles. The authors found that none of the results obtained with the simplified methods matched the transient simulations.

2.2.4. Evapotranspiration

The evapotranspiration flux, E T ( x , y , t ) [LT−1], is prescribed directly on the moving water table boundary, such that the net recharge in Equations (8) and (10) is N = R E E T . This condition is characterized by two conceptual parameters. The extinction elevation, z E X T , below which evapotranspiration from the water table ceases. Next, the ET surface elevation, z E T > z E X T , above which the evapotranspiration loss from the water table occurs at a fixed rate, E T m a x . In between, it is assumed that the E T flux varies linearly with the water table elevation. Hence,
E T ( x , y , t ) = { 0 z W T z E X T z W T z E X T z E T z E X T E T m a x z E X T < z W T < z E T E T m a x z W T z E T

2.2.5. Pumping and Injection Wells

Any vertical well is provided by the top, z t o p w , and bottom, z b o t w , vertical positions of the strainer. The prescribed flow rate, Q w [L3T−1], depends on the position of the water table relative to the well screen. When the water table position is higher than z t o p w , the flow rate is distributed among the formations perforated by the well screen; otherwise, when it intersects, the well screen length the flow rate is distributed among the formations lying between the water table position and the bottom of the well screen. A third case corresponds to the situation when the water table is below the bottom of the well screen, leading to shutting-down the well (i.e., Q w = 0 ). Additionally, we have implemented an option for a gradual decrease in the prescribed flow rate of a pumping well when the water table position starts to decrease below a threshold elevation, z t h r w ,   ( z b o t w < z t h r w < z t o p w ) , as was introduced in [9]. To summarize, the well boundary condition is formulated as follows:
Q w = { f Q Q w 0 Q w 0 < 0 Q w 0 Q w 0 > 0  
f Q = { s w 2 ( 2 ζ 3 s w + 3 ζ 2 ) z b o t w z W T z t h r w 1 z W T > z t h r w 0 z W T < z b o t w
where Q w 0 is the initially prescribed flow rate, f Q is a flow rate scaling factor for a pumping well, s w = z W T z b o t w is the saturated length in the well, and ζ = z t h r w z b o t w is the critical well length triggering a declining pumping flow rate.

2.2.6. River

The conceptualization of river–aquifer interactions assumes that a gaining or losing reach of the river is characterized by three constant parameters defining flow into or from the aquifer through a leaky riverbed layer. These are the river stage, h R , the bottom elevation of the riverbed layer, z R B , and the hydraulic conductance of the riverbed material, C R . When the water table elevation is above z R B , flow is directed towards the river using a mixed boundary condition where the negative flow rate depends linearly on the water table elevation. Otherwise, the flow rate Q R is directed towards the aquifer and a constant positive flow rate is prescribed at the water table directly, such that N = R E E T + Q R is just beneath the riverbed. To sum up, this boundary condition is expressed as:
Q R = { C R ( h R z W T ) z W T > z R B C R ( h R z R B ) z W T z R B  
An inherent assumption when the water table is below the riverbed is that the depth of the unsaturated zone is relatively thin, such that surface water enters the aquifer without a significant delay. Otherwise, similar to the recharge condition, we can estimate the unsaturated zone time lag, and, consequently, the time-dependent flow rate, Q R , is delayed.

2.2.7. Drain

The conceptualization of aquifer drainage using natural or artificial structures is somewhat similar in nature to the preceding condition type. A mixed boundary condition is usually adopted to estimate a drain flow rate, Q D R N , coming from the aquifer as the water table elevation remains above the drain elevation z D R N . When the water table elevation drops below the drain, this one exerts no action on the groundwater flow such that Q D R N = 0 . A mathematical statement of this condition reads:
Q D R N = { C D R N ( z D R N z W T ) z W T > z D R N 0 z W T z D R N
where C D R N is the hydraulic conductance of the drain.

3. The Discrete Moving Mesh Model

3.1. Mixed-Hybrid Finite Element Approximation (MHFEM)

We consider a spatial discretization of the groundwater flow, in Equation (4), using a mixed-hybrid finite element approximation on unstructured hexahedral grids. Let nf and ne be the total number of faces and elements in the mesh, respectively. The basic idea of the MHFEM is to simultaneously approximate the dependent variable and its first-order derivative [31,33]. Hence, a simultaneous discrete representation of the groundwater head and the specific discharge will be considered. According to [33], . q is constant inside each element, e, and facet, f, of the mesh. q e is uniquely determined when all normal fluxes through the six faces of a hexahedral element are perfectly known, such that:
q e = j = 1 6 q e , j   b j
where q e , j is the normal component of the groundwater discharge across the jth face of the hexahedral element e, and b j are the lowest-order Raviart–Thomas (RT0) basis functions, provided in Appendix A, and satisfying the fundamental relationship:
e . b j   d e = f i b j . n i   d f i = δ i j       i ,   j = 1 ,   ,   6
where n i is the outward normal unit vector to the face, fi, of element e. This implies that the basis function b i contributes a unit flux to face number i, and a null flux to all others.
The variational formulation of Darcy′s law is obtained by taking the scalar product of Equation (2) with test functions b i , integrating over the element, and using Green’s formula [33], yielding:
e ( K e 1 . q e ) . b i   d e = h e e . b i   d e j = 1 6 T h e , j   f j b i . n j   d f j       i = 1 ,   ,   6
where T h e , j is the trace of the groundwater head at face number j of element e, and he is the cell-centered groundwater head.
Inserting Equations (17) and (18) into Equation (19) yields:
j = 1 6 q e , j e ( K e 1 . b j ) . b i   d e = h e T h e , i       i = 1 ,   ,   6
Next, we introduce the local symmetric matrix [ A e ] 6 x 6 with entries:
A i j = e ( K e 1 . b j ) . b i   d e
Appendix D provides a pseudo-code for assembly of the local mixed-hybrid finite element matrix [ A e ] . Hence, Equation (19) is rewritten as:
q e , i = α e , i h e j = 1 6 A i j 1   T h j
α e , i = j = 1 6 A i j 1
where the shape factor α e , i is the sum at row number i of the inverse of the local resistance matrix Ae. From Equation (21), it is clear that this shape factor depends only on the geometry and the hydraulic conductivity components of the hexahedral element e.
The continuity of the normal component of the specific discharge is expressed for all internal faces, fi, shared by two neighbor elements, e1 and e2, such that:
q e 1 , i + q e 2 , i = 0
Equation (22) substituted twice (for e1 and e2) into Equation (24) yields:
α e 1 , i h e 1 + α e 2 , i h e 2 j = 1 6 [ ( A i j 1 ) e 1 ( A i j 1 ) e 2 ]   T h j = 0
Notably, when a constant flux q N is prescribed on Neumann boundary faces lying in a portion of the domain boundary Γ N (i.e., Q w in Equation (13) or Q D R N in Equation (16)), Equation (25) reduces to:
α e , i h e j = 1 6 A i j 1   T h j = q N
Likewise, we consider a Cauchy-type or general head boundary condition provided as q C = λ ( h h 0 ) , where q C is a head-dependent flux, λ is a hydraulic conductance (i.e., C R in Equation (15)), and h 0 is a prescribed head [L] (i.e., z W T or z R B in Equation (15)). q C is prescribed on all Cauchy boundary faces lying in a portion of the domain boundary Γ C . In this case, Equation (25) is rewritten as:
α e , i h e j = 1 6 A i j 1   T h j λ i T h i = λ i T h i , 0
where λ i and T h i , 0 are the prescribed conductance and head values at face i being part of Γ C .
Dirichlet boundary conditions of faces lying in Γ D (a portion of the domain boundary) are handled by simply moving the Dirichlet terms in the second term of Equation (27) to the right-hand side for rows belonging to non-Dirichlet faces, and directly substituting Equation (27) by the Dirichlet conditions (i.e., T h i = T h i , 0 ) for all Dirichlet rows.
Merging the continuity of the specific discharge provided by Equation (25), and the specific equations for Dirichlet, Neumann, and Cauchy boundary conditions presented by Equations (26) and (27), leads to the following matrix equation:
[ G ] { h } [ H ] { Th } = { D } + { N } + { C }
Entries of the globally assembled conductance matrices [ G ] n f x n e and [ H ] n f x n f are:
G i j = { α e j , i f i e j 0 f i e j
H i j = { f i ,   f j e ( A i j 1 ) e λ i f i Γ C f i ,   f j e ( A i j 1 ) e f i Γ C
Entries of the right-hand side column vectors { D } , { N } , and { C } resulting from prescribed Dirichlet, Neumann, and Cauchy boundary conditions on faces Γ D , Γ N , and Γ C , respectively, are:
D i = { H i j   T h i , 0 f i Γ D 0 f i Γ D  
N i = { q i f i Γ N 0 f i Γ N
C i = { λ i T h i , 0 f i Γ C 0 f i Γ C
Integration of Equation (1), using a first-order fully implicit Euler discretization in time, and making use of the special properties of the RT0 basis functions as presented in Equation (18), yields:
| V e | S p , e h e n + 1 h e n Δ t n + j = 1 6 q e , j n + 1 = Q s n + 1
where V e [L3] is the volume of element e, n is a time index, Δ t n = t n + 1 t n [T] is the time step at time level n , and Q s n + 1 = q s n + 1 V e [L3T1] is a volumetric source/sink term.
The insertion of Equation (22) into Equation (34) yields:
h e n + 1 ( 1 β e ) h e n β e j = 1 6 α e , j α e T h j n + 1 = β e Q s n + 1 α e
where β e = ω e 1 + ω e , ω e = α e Δ t n | V e | S p , e , and α e = i = 1 6 α e , i are the sum of all the entries of the inverse of the local matrix Ae. Equation (35) is rewritten in globally assembled matrix form as:
{ h } n + 1 [ M ] { h } n [ T ] { Th } n + 1 = { R } n + 1 + { Q } n + 1
where entries of the diagonal mass matrix [ M ] n e x n e and the trace matrix [ T ] n e x n f are:
M i j = { 1 β i e i = e j 0 e i e j
T i j = { β i α e i , j α i f j e i 0 f j e i
Entries of the right-hand side column vector { R } n e result from the Dirichlet face boundaries in Γ D whose terms are moved to the right-hand side, such that:
R i = f j ( e i Γ D ) β i α e i , j α i T h j , 0
The source–sink column vector { Q } n e has the following entries:
Q i = Q S i β e i α e i
Combining Equations (28) and (36) eliminates the cell-centered head unknowns and yields a linear system where the head traces at the faces of the mesh are the primary unknowns:
( [ H ] [ G ] [ T ] ) { Th } n + 1 = [ G ] ( [ M ] { h } n + { R } n + 1 + { Q } n + 1 ) { D } n + 1 { N } n + 1 { C } n + 1
The global sparse matrix ( [ H ] [ G ] [ T ] ) n f x n f is symmetric. Appendix E provides a pseudo-code to assembly of the global mixed-hybrid finite element matrix. A standard preconditioned conjugate gradient technique, such as ILU(0) preconditioning, is a suitable solver [50]. In situations where the matrix is not positive-definite owing to highly stretched elements or the occurrence of small dihedral angles, an M-matrix transformation always guaranties the existence of an ILU(0) preconditioner [50,51].
Once the linear system in Equation (41) is solved, the obtained trace heads vector, { T h } n + 1 , is substituted into Equation (36) to obtain the cell-centered heads at the current time level, { h } n + 1 , and, consequently, the specific discharge on all faces of the mesh using Equation (22) on a cell-by-cell basis. This post-processing step does not require, therefore, any linear system solution. Appendix F provides a pseudo-code for the hybridization procedure.
Notably, the assembled linear system for steady-state groundwater flow could be obtained simply by letting Δ t + , thus the coefficients β i are equal to 1 in Equations (37)–(39), and consequently the mass matrix, M, vanishes. Therefore, under steady-state conditions, Equation (41) simplifies to:
( [ H ] [ G ] [ T ] ) { Th } = [ G ] ( { R } + { Q } ) { D } { N } { C }
where the β i coefficients cancel out in the entries of the trace matrix [ T ] and the right-hand side vector { R } .
The global matrix ( [ H ] [ G ] [ T ] ) resulting from the MHFEM approximation is sparser than the global matrix resulting from a conforming finite element approximation. To illustrate this, we first recall that the non-zero entries of row number i in the global MHFEM and FEM assembled matrices result from face/node neighbors to face/node number i, respectively. In the finite element method, this neighborhood relationship is naturally defined as the faces or nodes shared by the same elements containing the target face or node. For a structured three-dimensional problem considering the conforming FEM on a hexahedral mesh, each node has 27 neighbors, including itself. Owing to the symmetry of the global conductance matrix, at most, 14 non-zeros are stored for each row of the sparse matrix. Alternatively, when considering an MHFEM for the same problem, each face has 11 neighbors, including itself. Therefore, because the matrix is symmetric, only (at most) six non-zero entries in each row of the resulting sparse matrix need to be stored. We conclude that even if the number of unknowns for an MHFEM approximation are greater than those for a conforming FEM for the same problem (as 2   n n < n f < 3   n n for a three-dimensional hexahedral mesh), the resulting global matrix associated with an MHFEM approximation is sparser and, therefore, not necessarily slower to process with sparse linear algebra operations such as matrix–vector multiplication, triangular solves, incomplete factorization, etc.
When the spatial grid is orthogonal and the hydraulic conductivity tensor is diagonal, the MHFEM approximation becomes equivalent to the standard FDM approximation (See Appendix B and Appendix C for demonstration).

3.2. Mesh Adaption to Water Table Movement

The computational stencil is a three-dimensional hexahedral mesh that is optionally unstructured in the plane and made from many layers to concept the aquifer geometry. Generally, then the given slices of the mesh conform with the stratigraphy of the aquifers’ bottoms and the land surface as much as possible. Hence, the mesh is composed of vertical columns and the mesh adaption proceeds independently in each column. This simplifies the interpolation (linear or harmonic) of the materials’ properties from the background geological mesh to the moving computational mesh. This also makes this process fast and practically tractable.
The phreatic surface position is locally updated using the boundary condition provided by Equation (5) at the nodal points in the surface. Prior to this, the computed head traces at the face centers of the upper slice of the mesh are interpolated to the nodal points of the same slice. This is because, usually, the coordinates of a finite element mesh are generally provided at the nodal points and not the face centers. If this updated value is above the land surface elevation, then the corresponding face is marked as a seepage face and Equation (11) takes place.
After updating the phreatic surface, it is necessary to do so for other moving nodes in the mesh. Only a set of the upper mesh slices is allowed to change their position while the others remain fixed. As suggested by other authors [10,23], a linear interpolation proves to work satisfactorily to find the position of the moving mesh slices. Optionally, a Laplacian smoothing technique is applied to the free nodes of the mesh. Those are the nodes not subject to any boundary condition. This operation is necessary for extreme situations where the vertical-to-horizontal scales ratio is higher than 1, which is typical for civil engineering structures such as earth dams (see Section 4.3).

3.3. Numerical Implementation of Boundary Conditions

The discrete systems of Equations (41) and (42) for groundwater flow under transient and steady-state conditions, respectively, have provisions for Dirichlet (i.e., constant or fixed head), Neuman, and Cauchy boundary conditions. Herein, we will discuss the numerical implementation for all the boundary conditions encountered in Section 2.2.
Notably, all discrete boundary conditions can use time series functions to allow for time-dependent variables. Spatial variability of any boundary condition parameter is dictated by the mesh resolution. Hence, each boundary face may have a different set of boundary conditions parameters.

3.3.1. Phreatic Surface

Any face on the top slice of the mesh is a phreatic surface boundary governed by Equation (8). At these boundaries, we have h z | Γ W T = K z z 1 | Γ W T q z , where K z z 1 is the component of the inverse of the hydraulic conductivity tensor, K , along the third dimension. Notably, this hydraulic conductivity tensor belongs to the stratigraphic element, which is crossed by the moving surface. Additionally, the z-component of the specific discharge, q z , is related to the normal component through the face, q , by the relationship q z = q cos θ , where θ is the cosine angle between the face normal and the vertical axis. Inserting the two previous relationships into Equation (8) yields:
q = 1 cos θ   ( 1 + N K z z 1 ) ( N + S y h t )  
Using an implicit first-order discretization of the time derivative term in Equation (43) as ( h t ) n + 1 T h n + 1 T h n Δ t n shows that a kinematic boundary face is a Cauchy-type boundary condition where the face hydraulic conductance and specified head are:
λ = S y Δ t n cos θ   ( 1 + N K z z 1 )  
T h 0 = N + S y T h n Δ t n cos θ   ( 1 + N K z z 1 )
Under steady-state conditions, the kinematic boundary condition reduces to a fixed flux on the given face, whose magnitude is presented as:
q = N cos θ   ( 1 + N K z z 1 )  
As stated earlier, a kinematic condition is dependent on other boundary conditions which are the components of the net recharge through recharge, evapotranspiration, and river faces.

3.3.2. Seepage Face

All the free top faces of the mesh are marked as potential candidates as seepage faces. The model implements equally an option where only a restricted set of these faces is included. This helps the algorithm to be faster, as areas where the groundwater depth is known to be significant are safely excluded (i.e., in high valleys, hills, etc.). Next, at each outer iteration, marked faces retain their phreatic surface conditions if the head is lower than the land surface position, z T O P , at the same horizontal coordinates. If the head exceeds this position, a constant head boundary condition, T h i = z T O P , is prescribed, where i is the face number. The numerical treatment of a constant head boundary is described in the next sub-section.

3.3.3. Constant Head

The head in any boundary face of the finite element mesh may be fixed during the simulation. The fixed head terms of a Dirichlet face, whose number is i , are moved to the right-hand side in all matrix rows, j i , according to Equation (39). It turns out that these are the rows for faces, j , who can share the same element with face i and that are not Dirichlet boundaries too. Furthermore, the coefficient matrix of the ith diagonal of the global matrix is set equal to one. In this way, the size of the matrix equations does not need to be modified prior or next to an iterative sparse matrix solution.

3.3.4. Recharge and Evapotranspiration

Recharge and evapotranspiration boundary conditions are directly applied on the water table boundary to calculate the net recharge through all its faces. Those are Neuman-type boundary conditions. The flux is internally multiplied by the face surface area prior to prescribing the boundary condition. Any recharge boundary condition is optionally delayed in time, while an evapotranspiration condition is applicable only when the calculated head in the face is higher than the specified extinction elevation at the horizontal coordinates of this face center.

3.3.5. Pumping/Injection Well

There are two alternatives to setup boundary conditions for a pumping/injection well when using an MHFEM numerical approximation. First, the flow rate boundary condition could be assigned to a group of boundary faces from which the fluid is extracted/injected. This is the case for a tubular geometry that represents a fictive well diameter. Another representation is possible by adding a source/sink term to an element containing the well. The first schematization is more flexible when greater accuracy at the proximity of the well is needed, such as when interpreting pumping tests. The second option allows the well effects to be incorporated at the regional scale without additional constraints on the mesh stencil near the wells, which is practical for field studies with a large number of wells. Both well’s representations are used in this paper. However, only the second option is supported within a moving mesh approach.
When representing a well boundary condition with a source/sink term, the difficulty with the moving mesh approach is to iteratively relocate the cells lying between the bottom and top screens of the well along a vertical column of cells. Provision is made to fit the bottom of the lower element in this column with the bottom of the well screen. Additionally, the flow rate is distributed among the cells in this column lying in the saturated interval of the well screen. In particular, for a pumping well, when the position of the water table in the well (i.e., the position of the upper surface in the well vertical column of cells) is below a given fraction of the well screen length, the flow rate is smoothly decreased according to Equation (14). In the preceding two situations, a source/sink term is added to the right-hand side vector { Q } in the discrete mass balance in Equation (36) according to Equation (40). Because the vector { Q } is multiplied by the matrix [ G ] in Equation (41) to obtain the final system of the equations on the mesh faces, during the assembly of the right-hand side of this system, a source/sink term such as in Equation (41) is distributed among the six j-faces of cell number i with the weighting shape factor α e i , j . When the water table position at the well falls below its screen bottom, this condition is deactivated.
The cellwise injection/pumping rates are adjusted, therefore, iteratively during each outer iteration of the overall algorithm.

3.3.6. River

A river boundary condition is applied to any face in the upper slice of the finite element mesh. When the water table position is above the riverbed elevation, z R B , of a river face, whose number is i , this results in a third type of boundary condition. The prescribed hydraulic conductance, C R , at this face is added to the ith diagonal of the global matrix according to Equation (30), and the product C R h R is added to the ith entry of the right-hand side column vector according to Equation (33). Otherwise, when the water table falls below the riverbed elevation of the face, we switch the boundary condition to the Neuman type with a prescribed flow rate that equals C R ( h R z R B ) > 0 . Hence, a river boundary condition is essentially a switching boundary condition. This status is updated at each outer (i.e., mesh level) iteration prior to solving a new linear system of equations and after updating the water table elevations.

3.3.7. Drain

Conceptually, there is two distinct ways to setup drains within an MHFEM approach. First, boundary faces could drain groundwater when the explicit geometric details of the structures, such as underground tunnels, are a part of the mesh. A second representation of this boundary condition may use a sink term to model drained water, such that all the drain geometric details are lumped at the cell center of an element. The first representation, while offering more flexibility, is more appropriate at detailed local scale models. The second representation is more flexible for regional scale studies. The first approach has some limitations when using the moving mesh method because of the difficulty to adjust or complete the details of the drain geometry when crossed by the water table. Hence, this representation could be used only when the drain is guaranteed to remain in the saturated zone during the simulation. Both methods are implemented in the developed simulator, but only the second one is used in this paper, and is therefore described as follows. When considering the second method, we search for the element center with the closest coordinates to the drain at ( x D R N ,   y D R N ,   z D R N ) . If this element exists (i.e., is inside the actual saturated region), a prescribed negative flow rate equal to C D R N ( z D R N z W T ) < 0 is attributed to this element, where z W T is the water table elevation at coordinates ( x D R N ,   y D R N ) , that is the highest vertical elevation of the mesh column containing the drain element, and C D R N is the hydraulic conductance of the drain. Otherwise, when the element is not found (i.e., is inside the virtual unsaturated region), this condition is inactive (i.e., no prescribed flow rate).

3.4. Computer Code

HydroTec© is the computer code implementing the adaptive MHFEM method presented in this paper. It is entirely programmed using the object-oriented (OO) C++ language. Although there was an increasing support of OO methodologies in other widely used computer languages by the groundwater modeling community, such as the Fortran 2003 standard, C++ was preferred because it is tightly integrated to high-level scripting languages (i.e., such as Python), which are planned to be used in the future to accelerate the addition of new features. Furthermore, the extended support of OO features and template-based metaprogramming were other features needed in this project. The core computational framework includes classes for sparse linear algebra, structured/unstructured mesh management, iterative sparse linear solvers, and boundary condition packages. Similar to MODFLOW, all supported types of boundary conditions could be loaded on-the-fly from their input files to the discrete sparse system of the equations. This flexibility supports the flexible and modular testing of groundwater flow models using different sets of hypothetical boundary conditions when using several conceptual models.
Finite element analysis is supported by a hierarchy of polymorphic face and element classes. For instance, any instance of an element class refers to either a quadrilateral or hexahedral object depending on the dimension of the problem being solved. The FEM assembly routines directly use the element class, enhancing the code genericity and leading to a less duplicate code.
A set of utility classes was designed to support input/output to different file formats, the management of time stress periods, and time series datasets (i.e., piezometric records, time-dependent boundary conditions, and time-dependent hydraulic conductivity or storage coefficients).
The code is under active development to support many groundwater modeling projects. A new MHFEM simulator for solute/heat transport was recently developed using the same framework and will be documented elsewhere.

4. Results

Five examples are presented to illustrate and validate the deforming mesh MHFEM model developed in this paper. The first evaluates the model by comparison with a two-dimensional analytical solution for drawdown in an infinite anisotropic confined aquifer. The second problem compares the model’s accuracy with MODFLOW + MODFPATH simulations for a two-dimensional highly heterogeneous formation. The third problem simulates the water table and seepage face on a rectangular earth dam, where the obtained results are compared with those of other authors. The fourth example uses a two-dimensional cross-section model of an unconfined aquifer to compare the free surface computed by HydroTec© with the same solution obtained by MODFLOW. The last example is an adapted three-dimensional model of an alluvial aquifer which is used as a test bed to compare the two models’ performance when simulating an unconfined groundwater flow.

4.1. Pumping in a Two-Dimensional Anisotropic Aquifer

Papadopulos provided an analytical solution for the hydraulic head drawdown, s, in an infinite anisotropic confined aquifer [52]. This analytical solution features a full anisotropic transmissivity tensor, T , and is used in this section for model verification. The radial drawdown resulting from the pumping well action depends on the pumping flow rate, Q w , the components of the transmissivity tensor, T , and the well function, W ( u ) as:
s ( x , y , t ) = Q w 4 π T x x T y y T x y 2 W ( u )
where u depends on the observation point coordinates (x, y), the aquifer specific storage, S p , the components of the transmissivity tensor, T , and the pumping time, t , such that:
u ( x , y , t ) = S p ( T x x y 2 + T y y x 2 2 T x y x y ) 4 t ( T x x T y y T x y 2 )
The well function, known also as the negative exponential integral, reads:
W ( u ) = u + e ξ ξ   d ξ
Notably, when T x x = T y y , and T x y = 0 , Equation (47) reduces to the well-known analytical solution provided by Theis for non-leaky confined aquifers [53].
The computational domain has an extension of 100 m in each side of the XY plane and has a uniform depth at 1 m. The physical parameters of the problem are provided in Table 1. Two cases are selected. The first corresponds to a diagonal transmissivity tensor, while the second case uses a full symmetric transmissivity tensor. Transmissivity components in the second case are obtained by rotating the main axes of the tensor ellipsoid, taken in the first case, by 30° in the anti-clockwise direction.
An automatically generated quality single-layer hexahedral mesh honoring the internal geometrical features of the model and the target maximal subdomain size is shown in Figure 2a. These mesh constraints are 1 m in the area bounded by a 10 m circle centered on the pumping well position, as shown in Figure 2b. Next, a second refinement zone is considered in the area enclosed by a circle with a radius of 1 m, as depicted in Figure 2c. Notably, a well radius of 0.1 m is adopted, and the well circumference is approximated with a regular pentagon (Figure 2c) and is located at the center of the domain. The hexahedral mesh contains a total of 3778 facets and 1876 elements.
The two transient simulations span a single stress period whose length equals 103 s. A constant time step of 2 s is selected. The hydraulic head drawdown is monitored at two piezometers, P1 and P2, as shown in Figure 2b. Figure 3a,b show the evolution of the numerically computed drawdowns with those analytically evaluated using Equation (47) at piezometers P1 and P2 for the diagonal and fully symmetric transmissivity tensors, respectively. In all cases, the MHFEM and analytical solutions are in good agreement. For the second case, it is clear that the hydraulic head drawdown is higher in P2 and lower in P1 than in the first case as expected, at least qualitatively, from the change in the transmissivity values provided in Table 1.

4.2. Flow in a Heterogeneous Porous Medium

Here, we consider a previously published model from the literature [28]. The authors illustrated the accuracy of the MHFEM approximation when compared to the conforming FEM using a hypothetical two-dimensional problem as shown in Figure 4. Their results suggested using particle tracking shape as a metric to evaluate the accuracy of the underlying numerical methods. In this work, we additionally use the particles forward travel time, τ , as a second metric for these accuracy assessments.
Figure 4 shows the 100 m by 100 m conceptual model considered in this test case. The domain comprises four hydraulic conductivity zones ranging from 102 m.d−1 to 10−3 m.d−1. Groundwater flows from north to south between two fixed head boundaries at 100 m and 99 m, respectively. A total of 19 particles were placed at an equal distance of 5 m at the top boundary tracing their flow paths. MODPATH simulation uses a semi-analytical particle tracking method [54]. Meanwhile, particle tracking simulation using the continuous specific discharge fields obtained by the MHFEM approximation is based on previously published algorithms [34].
The steady-state flow and particle tracking simulations are approximated with the MHFEM and the finite difference method (FDM) approximations. Their results refer to implementations in the HydroTec© and MODFLOW/MODPATH computer codes, respectively. In the first column of Figure 5, the simulated path lines using the MHFEM approximation with increasing grid resolutions from top to bottom are shown for the 10 × 10, 20 × 20, and 30 × 30 grid resolutions. Similarly, the second column shows the path lines simulated with MODPATH with identical grid resolutions. Notably, in all subfigures, the reference path lines are plotted with a red color for a direct comparison of both MHFEM and FDF approximations. These reference path lines are identical to the reference solution presented by the authors in the original publication [28].
Comparing Figure 5a,d indicates that the path lines resulting from the MHFEM approximation are closer to the reference solution than those computed with the FDM approximation. Particularly, it can be noticed that all flow paths which start traveling from the points located above the northmost least permeable block are penetrating it. This should not be the case according to the reference solution.
The same above-mentioned remarks hold for the other grid resolutions considered in this study. This can be concluded by comparing Figure 5b,e, then Figure 5c,f.
The improvement in the path lines shape is much better for the MHFEM approximation when the grid resolution steadily increases. For instance, Figure 5c shows that a 30 × 30 grid resolution using the MHFEM method is deemed sufficient for accurate particle tracking as the resulting flow paths follow very closely the reference solution that requires a much higher resolution (i.e., 320 × 320 grid). Notably, the mixed approximation yields a more accurate solution at a 20 × 20 grid resolution than a finite difference approximation at a 30 × 30 resolution. Hence, this confirms the superiority of the MHFEM method as a more appropriate choice when using coarse grids for particle tracking simulations, as concluded in previous works [28,34].
It Is worth noting that the FDM approximation necessitates a higher resolution than a 30 × 30 grid resolution to simulate accurate particle tracks. This is because one particle path is still penetrating the uppermost most impermeable block (Figure 5f), which is physically unrealistic.
As mentioned earlier, in reference [28], the authors did not consider the forward travel time as a metric for comparative assessments of the MHFEM with an alternative spatial approximation technique. This is performed in this work by calculating the percent relative error in the simulated travel time of a flowing particle at the exit boundary as:
error ( % ) = 100 × | τ r e f τ | / τ r e f
where τ r e f is the travel time of the reference flow path at the exit boundary simulated on the 320 × 320 fine grid resolution with a mixed approximation, while τ is a similar quantity evaluated at the grid resolution of interest.
Figure 6a,b depict the relative errors obtained from the MHFEM and FDM approximations for different grid resolutions for particles released at coordinates (25, 100) m and (95, 100) m, respectively. They show that travel time errors for the MHFEM approximation are small even for the lowest resolution. Indeed, they do not exceed 2.5% and 1.7% for the considered flow lines, respectively. These errors are much smaller than those obtained when using an FDM approximation. Moreover, the relative error in the MHFEM decreases monotonically as the grid resolution increases, confirming the suitability and superiority of this numerical scheme for particle tracking analysis on coarse grids. This is unlikely the case with the FDM approximation as the travel time convergence when increasing the grid resolutions is erratic.

4.3. Two-Dimensional Free Surface in a Homogeneous Earth Dam

In this third example, a homogeneous cross-section model of a regular earth dam is considered. The model domain has a length of 5 m and a height of 10 m. This area is divided into uniform quadrilateral elements at 0.1 m on each side. Therefore, the mesh stencil contains 5000 elements and 10,150 edges. The hydraulic conductivity equals 1 m/d everywhere. The flow occurs between two fixed head boundaries. On the left, the head was fixed to 10 m, while this was set to 2 m on the right side. All nodes in the mesh, except those lying at the bottom, and fixed head faces were allowed to change their position. Notably, in this problem, Laplacian mesh smoothing was enabled because simply moving the nodes along the vertical lines could result in highly distorted elements or even elements with self-intersecting edges owing to the high vertical-to-horizontal aspect ratio of the problem.
Figure 7 shows the newly obtained mesh and the position of the steady-state phreatic surface. The length of the seepage face results directly from the water table height on the right side. There was no need to specify potential seepage faces at this boundary prior to the simulation as would be the case within a fixed mesh approach. The simulated water table height is compared with other numerical solutions gathered from the literature. There was a good agreement with the solutions provided by [55,56]. The first solution developed in [55] uses a bilinear fixed mesh conforming FEM model. The second solution from [56] was computed with an FDM technique implemented in a spreadsheet.

4.4. Two-Dimensional Free Surface in a Pumped Phreatic Aquifer

In this third example, we consider a 6 Km long and 50 m deep vertical section of a phreatic aquifer bounded by two constant head boundaries at 50 m and 35 m from west to east, respectively. The aquifer is assumed to be homogeneous and isotropic with a hydraulic conductivity at 5 m/d. The mean groundwater recharge was assumed to linearly increase from 0 to 0.05 m/d between the ground surface and a depth of 35 m following previous observations [48]. A partially penetrating well located at the center of the model extension withdraws groundwater with a flow rate equal to 62.5 m3/d. The well is screened in the lower half of the aquifer thickness. Owing to the aquifer homogeneity, the flow rate is equally divided among the cells lying along the well perforation.
A structured mesh was adopted to enable a fair comparison with MODFLOW 2005 [8]. The grid has 120 columns and 10 layers. Hence, the finite-difference simulation with MODFLOW solves a nonlinear problem with 1200 unknowns, while the MHFEM simulation with HydroTec© solves a problem with 2530 unknowns using the deforming mesh algorithm developed in this paper.
The two compared codes use the same preconditioned iterative method, and the same convergence control parameters, for the inner linear systems in the outer iteration loop, facilitating a fair intercomparison. Namely, the maximum number of iterations (i.e., 900) and the tolerance of inner convergence of groundwater heads (i.e., 10−12). In the FDM- based simulation, a Picard iteration is adopted with a convergence closure equal to 10−3 m. The same parameter was used as a convergence closure for the heads in the mesh adaption loop. Likewise, a maximum number of 100 outer iterations was similarly taken for both algorithms. All 10 upper horizontal lines of the 2D mesh were allowed to move during the iterative adaption process.
Figure 8 shows the final mesh shape and head distribution after the convergence of the adaptive MHFEM solution. In particular, the groundwater surface is smooth and all its nodes (i.e., the med-edges) explicitly satisfy the prescribed external convergence criteria. This is not the case for MODFLOW’s FDM solution, because the water table is inside the grid blocks and must be found by post-processing (i.e., interpolation) from the simulated groundwater heads in the center of their cells. This explains why the moving mesh solution always provides directly a much more accurate position of the water table.
The convergence behavior of mesh adaptation is shown in Figure 9. This is about 3.2 and is therefore faster than a linear trend. More importantly, it outperforms the convergence rate of Picard’s method and even the well-known quadratic rate of Newton’s iteration method. Only 8 iterations were needed in the adaptive mesh MHFEM model, while 32 Picard iterations were performed by MODFLOW. In this case, even for this simple problem, convergence was more difficult to reach with the chosen external stopping criterion. It was necessary to move to a small under-relaxation factor (i.e., 0.25) for MODFLOW to converge.
On the other hand, the MHFEM model took less than 2 s CPU time for this steady-state simulation, while the FDM solution with a fixed mesh took about 11 s.
Figure 10 shows the simulated water table with the MHFEM and FDM approximations. While there is a good agreement between the drawdowns calculated with the two methods, small discrepancies are observed far from the well. The maximum error does not exceed 0.74 m, which amounts to 15% of the vertical grid size (i.e., 5 m). As concluded earlier, the differences are due to interpolation artifacts inside the fixed mesh. This error would be greater as the mesh becomes coarser along the vertical direction. Furthermore, the MHFEM approximation with an adaptive mesh is always more accurate in the particular case of groundwater abstraction. This is because the mesh density becomes higher than the initial fixed mesh when it contracts due to pumping. This observation would also apply to the groundwater head distribution in the moving saturated region.
We conclude from the obtained results relative to this third test problem that the moving mesh method is an attractive technique to accurately simulate groundwater flow in unconfined aquifers subject to stresses such as pumping, recharge, etc. The next test problem could be regarded as a three-dimensional extension that would serve as a basis to consolidate these findings.

4.5. Three-Dimensional Free Surface in a Pumped Phreatic Aquifer

This problem is adapted from a field case study in a small alluvial aquifer. The domain geometry is shown in the inset of Figure 11. Except for the left constant head boundary, this three-dimensional problem inherits many characteristics of the previous example. A second difference concerns the well properties including its position, well rate, and perforated length. The well is positioned at coordinates x = 3700 m and y = 4000 m. The pumping rate equals 6000 m3/d. The well is screened in the lower 10 m of the aquifer.
As in the previous example, a structured and relatively coarse mesh was adopted to enable a direct comparison with MODFLOW 2005. The grid has 60 columns, 70 rows, and 4 layers. Hence, when the mesh is fixed, the horizontal cell dimension is close to 100 m and the vertical cell size is uniform at 12.5 m. The finite-difference simulation with MODFLOW solves a nonlinear problem with 16,284 unknowns, while the MHFEM simulation with HydroTec© solves a problem with 53,435 unknowns.
The convergence control parameters for the two simulations are identical to those selected in the previous problem, except for the nonlinear/mesh convergence closure criterion, which is equal to 10−2 m. This is because the MODFLOW simulation does not converge when the last parameter was increased up to 10−3 m. Contrarily, a higher outer convergence criterion does not pose any challenge to the adaptive mesh model as only three more mesh iterations were required. All four upper horizontal slices of the 3D mesh were allowed to move during the iterative adaption process.
Figure 11 shows a cutaway view of the final mesh shape and head distribution at steady-state conditions. Even with a coarse mesh, the model successfully computes and visualizes the cone of depression with a rapidly changing curvature of the free surface. As expected, and also depicted in Figure 11, the hydraulic gradient is the highest in this zone.
The convergence behavior of mesh adaptation is shown in Figure 12. This is about 1.9 and is therefore faster than a linear trend and close to the quadratic rate of Newton’s iteration method. Only 13 iterations were needed by the adaptive mesh MHFEM model, while 54 Picard iterations were performed by MODFLOW. The MHFEM model took 32 s CPU time for this steady-state simulation, while the FDM solution with a fixed mesh took about 48 s.
Figure 13 shows the simulated water table profile with the MHFEM and FDM approximations. This profile is parallel to the X-axis and crosses the well position. The drawdowns calculated with the two methods are not in good agreement at the well as the difference amounts to 2.65 m. This is a consequence of the coarse vertical discretization of the fixed mesh. Hence, post-processing the water table position, including the drawdown at the well, can lead to higher errors as the vertical cell size increases. The results of this test problem confirm the ability of the MHFEM approximation with an adaptive mesh to deliver an accurate representation of the water table geometry, which is free from interpolation errors and less prone to the vertical discretization of the mesh.

5. Discussion

The adaptive meshing technique developed in this paper is relatively simple to implement into the framework of an existing code for modeling saturated groundwater flow. This upgrading effort is somewhat equivalent to upgrading an FEM fixed mesh model to support unconfined flow using a Picard loop. For the moving mesh technique, this outer loop corresponds to mesh iteration. Hence, the nonlinearity is removed from the flow process in favor of mesh adaption. The provided examples demonstrate that this approach enhances the model’s robustness, provided that the geological model complexity can survive the layered mesh adaption methodology. When this process succeeds, the moving mesh approach is unbeatable. This is because the underlying preconditioned linear systems using the M-matrix transformation for the intermediate linear systems are always guaranteed to converge. Another robustness aspect of the moving mesh approach is its resilience with respect to the outer tolerance or mesh convergence closure. As can be seen from the provided examples, convergence is no problem for values up to 10−3 m, while this impacts even the simplest models with a fixed mesh approach.
The provided problems in this work support that the moving mesh technique is more efficient than the fixed mesh. Indeed, the convergence rate is higher than linear and is close to quadratic for three-dimensional models. There is even room for improvement as the PCG iteration could be replaced with a more efficient algebraic multigrid preconditioner when targeting large-scale discrete problems [57].
The chief advantage of the moving mesh approach is obtaining an accurate water table interface shape, which is free from the uncertainties related to coarse discretization and spatial interpolation from groundwater heads. The water table interface is derived naturally because it is part of the adapted mesh geometry.
The moving mesh method has few limitations with regard to the current state-of-the-art method. The most challenging obstacle to its universal application is the occurrence of the random heterogeneity of the porous medium. This is the case of hydraulic conductivity realizations mapped from geostatistical techniques. When the hydraulic conductivity is different from cell to cell along all directions, it becomes almost impossible to find any practical guidelines to correctly adapt the mesh. Another minor limitation is the inability to simulate perched water table aquifers. This is because the continuity of the saturated and unsaturated zones is implicitly taken as an initial assumption. Furthermore, the simple mesh adaption scheme used in this work cannot handle the extreme situation when the water table becomes nearly vertical, spanning a large number of vertical cells in the model. This is the case of a heterogeneous earth dam where, for instance, two pile sheet materials are placed along a vertical discontinuous interface. However, this event is impossible to occur in realistic unconfined aquifers characterized by their large horizontal-to-vertical aspect ratio.
Overall, we think that any widely used groundwater modeling software should support not only the fixed mesh approach. This will not only enlarge its scope and usability, but it will advantageously enable users to cross-check the two approaches for the same application.
An interesting extension of the moving mesh approach is handling the general case of multiple free and moving groundwater interfaces. A typical example is the situation of two exploited aquifers separated by an impermeable layer where a dewatered zone results from extensive pumping of the lower aquifer. In this case, not only the mesh geometry accounts but also its topology (i.e., connectivity) should be adapted accordingly to exclude the dewatered volume from the mesh. Advances in this direction are awaited in future works.

6. Conclusions

This paper deals with a novel moving mesh approach for modeling unconfined groundwater flow with a mixed-hybrid finite element approximation. A layered mesh is adjusted to fit with the water table interface modeled as a kinematic boundary condition. The model accounts for many types of boundary conditions needed in practical applications. This approach has been validated using analytical solutions and numerical solutions computed with industry-standard groundwater FDM modeling codes such as MODFLOW. The following conclusions can be drawn from this study:
  • The mixed-hybrid finite element approximation is suitable for highly heterogeneous and anisotropic aquifers. It provides accurate solutions for the discharge fields even for coarser grids than the FDMs. Therefore, it is a suitable approximation technique prior to advective particle tracking or solute/heat transport simulations.
  • The moving mesh technique based on the MHFEM approximations provides an accurate solution for the water table interface, groundwater heads, and specific discharge fields. The water table is accurate because it is free from cellwise interpolation and independent of the vertical grid size, as used in fixed mesh methods.
  • The robustness of the moving mesh method cannot be surpassed by a fixed mesh alternative. The nonlinearity associated with the latter is replaced by a series of saturated flow problems which are always guaranteed to converge for some preconditioners. Moreover, the moving mesh outer iteration loop accepts higher convergence closure parameters.
  • Owing to the simplicity of the mesh adaption scheme, the moving mesh model proves to be equally efficient. For all of the studied problems, it was faster than the Picard loop-based FDM solutions when using the same resolution.
  • The efficiency of the moving mesh MHFEM model is supported by an almost quadratic rate of convergence.
  • The presented samples prove that it is possible to confidently develop practical groundwater applications using HydroTec© code. Nevertheless, additional benchmarks and testing will strengthen users’ confidence.

Author Contributions

Conceptualization, M.A.S.; methodology, M.A.S. and A.L.; software, M.A.S.; validation, M.A.S. and A.L.; formal analysis, M.A.S.; investigation, M.A.S.; writing—original draft preparation, M.A.S.; writing—review and editing, A.L.; visualization, M.A.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to acknowledge the guest editors of the special issue “Flow and Transport Processes in Groundwater Systems” for the invitation to write this paper. We are also thankful to all the reviewers for their constructive comments.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Expressions of the Vector Basis Functions

There exist six vector basis functions for a local interpolation of the specific discharge inside a hexahedral element as in Equation (17). Herein, we use the lowest-order Raviart–Thomas (RT0) basis functions defined in the local coordinates ( ξ , η , ζ ) as listed below:
b 1 ( ξ , η , ζ ) = 1 8   [ ξ 1 ,   0 ,   0 ] T
b 2 ( ξ , η , ζ ) = 1 8   [ ξ + 1 ,   0 ,   0 ] T
b 3 ( ξ , η , ζ ) = 1 8   [ 0 ,   η 1 ,   0 ] T
b 4 ( ξ , η , ζ ) = 1 8   [ 0 ,   η + 1 ,   0 ] T
b 5 ( ξ , η , ζ ) = 1 8   [ 0 ,   0 ,   ζ 1 ] T
b 6 ( ξ , η , ζ ) = 1 8   [ 0 ,   0 ,   ζ + 1 ] T
In local coordinates, the original hexahedral element is mapped into a cube (i.e., the reference element) whose corners are located at ξ = ± 1 , η = ± 1 , ζ = ± 1 . It is quite straightforward to show that all these basis functions satisfy Equation (18).

Appendix B. MHFEM Approximation on a Rectangular Element with a Diagonal Hydraulic Conductivity Tensor

This appendix derives an expression for the mixed-hybrid finite element approximation on a rectangular element with lengths Δ x , Δ y , and Δ z along the X, Y, and Z axes, respectively. We assume, furthermore, that the diagonal hydraulic conductivity tensor is diagonal and presented as diag ( K x , K y , K z ) . Owing to the orthogonality of the vector basis functions, it is easy to show from Equation (21) that the matrix [ A e ] is also diagonal in this case. It is presented as:
A e = V   diag ( 1 K x , 1 K x , 1 K y , 1 K y , 1 K z , 1 K z )
where V = Δ x Δ y Δ z is the element volume and local face numbers are ordered along the faces orthogonal to X, Y, and Z axes, respectively.
When the hydraulic conductivity tensor is diagonal, Equation (22) stipulates that the normal specific discharge component on any face depends only on the heads at the cell center and at the same face. This reduces Equation (22) to a form of finite difference (i.e., two points) approximation. The only difference with a standard FDM approximation is that the latter does not consider the trace heads as additional variables, and the local specific discharge fluxes are written exclusively in terms of cell centers. For instance, the flux along a face orthogonal to the X-axis is deduced from Equations (22) and (23) combined with Equation (A1) as:
q x = K x L x ( h e T h x )
where L x is the distance between the face and the cell center.
Equation (A8) is a simple statement of Darcy’s law between two specific points in the element. Hence, this establishes the equivalence between the MHFEM and the FDM approximations for a rectangular mesh and a diagonal hydraulic conductivity tensor. However, the MHFEM has an interesting plus as the specific discharge does not result from post-processing the cell-centered heads, meaning they are, therefore, more accurate.
Notably, under the particular conditions of mesh orthogonality and a diagonal hydraulic conductivity tensor, there is no need to use a numerical quadrature integration formula to assemble the local element matrix, as discussed in Appendix D, as Equation (A7) would be used directly.

Appendix C. Equivalence with the FDM Approximation for a Diagonal Hydraulic Conductivity Tensor

We consider the two elements shared by a face orthogonal to the X-axis. We denote by h 1 and h 2 the heads at the centers of cells number 1 and 2, respectively. T h x is the trace of the head at the face. As in the previous appendix, we assume a rectangular cell and a diagonal hydraulic conductivity tensor. Inserting Equation (A8) twice in Equation (24) yields:
T h x = K x , 1 L x , 1 h 1 + K x , 2 L x , 2 h 2 K x , 1 L x , 1 + K x , 2 L x , 2
Inserting Equation (A9) into Equation (A8) provides an expression of the specific discharge as a function of the heads in the cell centers:
q x , 1 = K e q ( h 1 h 2 ) L x , 1
K e q = K x , 1 L x , 1 K x , 2 L x , 2 K x , 1 L x , 1 + K x , 2 L x , 2
where K e q is a weighted harmonic mean equivalent hydraulic conductivity. The weight is the inverse length between the cell center and the face.
Hence, the MHFEM and FDM approximations are equivalents for a rectangular mesh with a diagonal hydraulic conductivity tensor. In particular, when L x , 1 = L x , 2 (i.e., uniform X grid size), the equivalent hydraulic conductivity in Equation (A11) is simply the harmonic mean of the hydraulic conductivities of the two cells traditionally used in the FDM.

Appendix D. Assembly of the Local Mixed-Hybrid Finite Element Matrix

This appendix provides and documents a generic pseudo-code of the function or subroutine that assembles the local matrix for element, e, for a mixed-hybrid finite element approximation. Inputs for this function are the element structure, e, holding the local numbering of its six faces, and the 3-by-3 hydraulic conductivity matrix, Ke, for the soil material filling the same element. Outputs are the six-by-six inverse of the local mixed-hybrid matrix, A e 1 , the six-by-one array of its row-sums (i.e., the shape factors α e , i   i = 1 , , 6 ), the sum of all entries in A e 1 : α e , and the hexahedral element volume Ve.
1 Begin Function AssembleLocalMatrix (Input e, Ke) Returns A e 1 , α e , i , α e , Ve
2      Obtain K−1 the inverse of Ke
3      Set Ve to 0
4      Set Ae to 0
5      Begin Cycle for Gauss in [1,6]
7         Obtain the transpose of the Jacobian JT and its determinant |J|
8         Obtain the mixed basis function bGauss at Gauss point
9         Add WGauss x |J| to Ve
10       Begin Cycle for face i in [1,6]
11         Obtain Wi = K−1 x JT x bGauss(i)
12         Begin Cycle for face j in [1,6]
13         If j is lower than or equals i Add ((Wi x bGauss(i)) . (JT x bGauss(j)))/|J| to Ae(i,j)
14         End Cycle
15       End Cycle
16      End Cycle
17      Begin Cycle for face i in [1,6]
18         Begin Cycle for face j in [1,6]
19           If j is greater than i Set Ae(j,i) to Ae(i,j)
20         End Cycle
21      End Cycle
22      Obtain A e 1 the inverse of Ae
23      Set α e , i to j = 1 6 A e 1 ( i , j ) for i in [1,6]
24      Set α e to i = 1 6 α e , i
25  End Function
The direct translation of this pseudo-code to any high-level computer programming language, such as Fortran or C/C++, will result in an efficient implementation.
The first loop between lines 5 and 16 calculates the lower triangular part of the local mixed-hybrid matrix by adding the respective contributions from the six Gauss points of the integration. This loop expects three external function calls to obtain (i) the Gauss point weight, (ii) the three-by-three transpose of the Jacobian matrix transformation, JT, from the real to the reference element and its determinant |J|, and (iii) the vector basis function at the Gauss point. This loop also calculates the element volume by adding each Gauss point contribution, as shown in line nine.
Because the matrix is symmetric, the loop between lines 17 and 21 simply calculate the strictly upper triangular entries of this matrix from its lower triangular part.
Finally, the loop between lines 23 and 24 simply calculate the shape factors α e , i   i = 1 , , 6 and their sum α e , respectively.

Appendix E. Assembly of the Global Mixed-Hybrid Finite Element Matrix

This appendix provides and documents a generic pseudo-code of the function or subroutine that assembles the global matrix resulting from a mixed-hybrid finite element approximation of a groundwater flow problem. Inputs for this function are the current time step size (only for a transient simulation), Δ t , the logical nf-by-1 array FixedHead indicating whether a face is of a fixed head type or not, the ne-by-1 array CellHead of cell-centered heads at the previous time step, and the ne-by-1 array for the cellwise source/sink values Q at the current time step. For steady-state conditions, the CellHead array may take any initial set of head values. The number of elements ne in line 2 may be deduced from the shape of the CellHead array.
1 Begin Function AssembleGlobalMatrix (Input Δ t , FixedHead, CellHead, Q) Returns Matrix, Rhs
2      Begin Cycle for e in [1,ne]
3            Obtain the element hydraulic conductivity tensor Ke
4            AssembleLocalMatrix (e, Ke) Returns A e 1 , α e , i , α e , Ve
5            Obtain S e the specific storage coefficient in element e
6            If the flow is steady-state Set β e to 1
7            If the flow in transient Set β e to ( α e Δ t / S e Ve + α e Δ t )
8            Begin Cycle for face i in [1,6]
9                  Obtain fi the global number of face i
10                 If fi is a constant head face Cycle
11                 Begin Cycle for face j in [1,6]
12                         Obtain fj the global number of face j
13                         If fj is not a constant head face
14                            If fj is lower than or equals fi Add A e 1 (i,j) − β e α e , i α e , j α e to Matrix(fi,fj)
15                         Else
16                            Add -( A e 1 (i,j)- β e α e , i α e , j α e ) FixedHead(fj) to Rhs(fi)
17                         End If
18                 End Cycle
19                 Add Q(e) β e α e , i α e   + α e , i ( 1 β e )   CellHead(e) to Rhs(fi)
20           End Cycle
21     End Cycle
22   End Function
Outputs are the nf-by-nf sparse global mixed-hybrid finite element Matrix corresponding to ( [ H ] [ G ] [ T ] ) in Equations (41) and (42), the nf-by-1 array Rhs corresponding to all right-hand-side vectors in Equation (41), excluding Neuman and Cauchy arrays, which are separately assembled when processing the boundary conditions.
The most outer loop between lines 2 and 21 sequentially processes all elements to fill contributions to the global matrix from local hydraulic interactions between any couple of two faces among the six faces. To do this, two inner loops between lines 8 and 20 for the first, and lines 11 and 18 for the second, record these local interactions. For each local couple of faces i and j, their global numbers fi and fj are gathered in lines 9 and 12, respectively. Then, only for the lower triangular entries of the sparse matrix (i.e., when fj <= fi) is the value in line 14 added to the position at row fi and column fj. This operation is performed only when both faces fi and fj are not of the Dirichlet type, as indicated by lines 10 and 13. For constant head faces, their contributions are moved to the right-hand-side vector, as stated in line 16. Line 19 adds contributions deriving from the weighted contributions of the element source–sink and the cell-centered head at the previous time step. Other contributions to the right-hand-side vector from Neuman and Cauchy boundary conditions are added when processing the boundary conditions. The code works for steady-state and transient conditions have the provision to adjust the local coefficient β e , as shown in lines 6 and 7.
Notably, the global matrix is stored in a suitable sparse matrix format, such as the compressed sparse row (CSR) storage. Only the lower triangular part is stored.

Appendix F. Retrieving Cell-Centered Heads and Normal Face Fluxes from Head Traces

This appendix provides and documents a generic pseudo-code of the function or subroutine that implements hybridization for a mixed-hybrid finite element approximation.
1 Begin Function Hybridization (Input Δ t , FaceHead, CellHead, Q) Returns CellHead, FaceFlux
2    Begin Cycle for e in [1,ne]
3      Obtain the element hydraulic conductivity tensor Ke
4      AssembleLocalMatrix (e, Ke) Returns A e 1 , α e , i , α e , Ve
5      Obtain S e the specific storage coefficient in element e
6      If the flow is steady-state Set β e to 1
7      If the flow in transient Set β e to ( α e Δ t / S e Ve + α e Δ t )
8      Set sum to 0
9      Begin Cycle for face i in [1,6]
10           Obtain fi the global number of face i
11           Add α e , i FaceHead(fi) to sum
12        End Cycle
13        Set CellHead(e) to ( 1 β e )   CellHead(e) + (Q(e)+sum) β e α e
14        Begin Cycle for face i in [1,6]
15           Obtain fi the global number of face i
16           Set FaceFlux(i,e) to α e , i CellHead(e)
17           Begin Cycle for face j in [1,6]
18           Obtain fj the global number of face j
19           Add - A e 1 (i,j)FaceHead(fj) to FaceFlux(i,e)
20        End Cycle
21     End Cycle
22    End Cycle
23   End Function
The hybridization process is known as the local operations where the cell-centered heads and the normal fluxes to faces are retrieved once the head traces on the faces are determined. Inputs for this function are the current time step size (only for a transient simulation), Δ t , the nf-by-1 array FaceHead storing the computed heads at the faces of the mesh, the ne-by-1 array CellHead where cell-centered heads at the previous time step are stored, and he ne-by-1 array for the cellwise source/sink values Q at the current time step.
Outputs are the ne-by-1 array CellHead for cell-centered heads at the current time step, and the six-by-ne array FaceFlux where the normal components of the specific discharge at each face of an element for the current time step are stored.
The most outer loop between lines 2 and 22 sequentially processes the elements to calculate the heads at their centers and to add each elemental contributions to the water flux through connected faces in the mesh. The loop between lines 9 and 11 simply calculates the weighted sum of trace heads in all the faces connected to an element. This term is substituted to calculate the cell-centered head at the current time step, as stated in line 13. Additionally, the first term in this line corresponds to the cell-centered head at the previous time step. The loop between lines 14 and 21 adds elemental contributions to the face fluxes according to the expression provided in Equation (22).

References

  1. Smith, M.; Cross, K.; Paden, M.; Laban, P. (Eds.) Spring—Managing Groundwater Sustainably; IUCN: Gland, Switzerland, 2016. [Google Scholar] [CrossRef] [Green Version]
  2. Johnson, T.D.; Belitz, K.; Kauffman, L.J.; Watson, E.; Wilson, J.T. Populations using public-supply groundwater in the conterminous U.S. 2010; Identifying the wells, hydrogeologic regions, and hydrogeologic mapping units. Sci. Total Environ. 2022, 806, 150618. [Google Scholar] [CrossRef]
  3. Anderson, M.P.; Woessner, W.W.; Hunt, R.J. Applied Groundwater Modeling—Simulation of Flow and Advective Transport, 2nd ed.; Elsevier: London, UK, 2015. [Google Scholar] [CrossRef]
  4. Langevin, C.D.; Hughes, J.D.; Provost, A.M.; Banta, E.R.; Niswonger, R.G.; Panday, S. Documentation for the MODFLOW 6 Groundwater Flow (GWF) Model; U.S. Geological Survey Techniques and Methods, Book 6, Chapter A55; U.S. Geological Survey: Reston, VA, USA, 2017; 197p. [Google Scholar] [CrossRef] [Green Version]
  5. Diersch, H.-J.G. FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media; Springer: Berlin/Heidelberg, Germany, 2013. [Google Scholar] [CrossRef]
  6. Bear, J. Dynamics of Fluids in Porous Media; Elsevier: New York, NY, USA, 1972. [Google Scholar] [CrossRef] [Green Version]
  7. Bear, J.; Cheng, A.H.D. Modeling Groundwater Flow and Contaminant Transport; Springer: Berlin/Heidelberg, Germany, 2010. [Google Scholar] [CrossRef]
  8. Harbaugh, A.W. MODFLOW-2005, the U.S. Geological Survey Modular Ground-Water Model—The Ground-Water Flow Process; U.S. Geological Survey Techniques and Methods, Book 6, Chapter A16, variously paged; U.S. Geological Survey: Reston, VA, USA, 2015. Available online: https://pubs.usgs.gov/tm/2005/tm6A16/ (accessed on 2 December 2022).
  9. Niswonger, R.G.; Panday, S.; Ibaraki, M. MODFLOW-NWT, A Newton Formulation for MODFLOW-2005; U.S. Geological Survey Techniques and Methods, Book 6, Chapter A37; U.S. Geological Survey: Reston, VA, USA, 2011; 44p. Available online: https://pubs.er.usgs.gov/publication/tm6A37 (accessed on 2 December 2022).
  10. Knupp, P. A moving mesh algorithm for 3-D regional groundwater flow with water table and seepage face. Adv. Water Resour. 1996, 19, 83–95. [Google Scholar] [CrossRef]
  11. Keating, E.; Zyvoloski, G. A Stable and efficient numerical algorithm for unconfined aquifer analysis. Groundwater 2009, 47, 569–579. [Google Scholar] [CrossRef] [PubMed]
  12. Desai, C.S.; Li, G.C. A residual flow procedure and application for free surface flow in porous media. Adv. Water Resour. 1983, 6, 27–35. [Google Scholar] [CrossRef]
  13. Desai, C.S.; Baseghi, B. Theory and verification of residual flow procedure for 3-D free surface seepage. Adv. Water Resour. 1988, 11, 192–203. [Google Scholar] [CrossRef]
  14. Sbai, M.A. Modelling Three-Dimensional Groundwater Flow and Transport by Hexahedral Finite Elements. Ph.D. Thesis, Free University of Brussels, Belgium, Brussels, 1999. [Google Scholar]
  15. Larabi, A.; De Smedt, F. Numerical solution of 3-D groundwater flow involving free boundaries by a fixed finite element method. J. Hydrol. 1997, 201, 161–182. [Google Scholar] [CrossRef]
  16. Kazemzadeh-Parsi, M.J.; Daneshmand, F. Three-dimensional smoothed fixed grid finite element method for the solution of unconfined seepage problems. Finite Elem. Anal. Des. 2013, 64, 24–35. [Google Scholar] [CrossRef]
  17. Koo, M.-H.; Leap, D.I. Modeling three-dimensional groundwater flows by the body-fitted coordinate (BFC) method: I. Stationary boundary problems. Transp. Porous Media 1998, 30, 217–239. [Google Scholar] [CrossRef]
  18. Koo, M.-H.; Leap, D.I. Modeling three-dimensional groundwater flows by the body-fitted coordinate (BFC) method: II. Free and moving boundary problems. Transp. Porous Media 1998, 30, 345–362. [Google Scholar] [CrossRef]
  19. Neuman, S.P.; Witherspoon, P.A. Variational principles for confined and unconfined flow of groundwater. Water Resour. Res. 1970, 6, 889–897. [Google Scholar] [CrossRef] [Green Version]
  20. Neuman, S.P.; Witherspoon, P.A. Analysis of nonsteady flow with a free surface using the finite element method. Water Resour. Res. 1971, 7, 611–623. [Google Scholar] [CrossRef]
  21. Purkey, D.R.; Wallender, W.W.; Islam, N.; Fogg, G.E.; Sivakumar, B. Describing near surface, transient flow processes in unconfined aquifers below irrigated lands: A deforming finite element model for heterogeneous aquifers. J. Hydrol. 2006, 330, 435–443. [Google Scholar] [CrossRef]
  22. Toufigh, V. Constrained optimization based F.E. mesh deforming algorithm for unconfined seepage problems. Applied Math. Model. 2016, 40, 6754–6765. [Google Scholar] [CrossRef]
  23. Kourakos, G.; Harter, T. Simulation of unconfined aquifer flow based on parallel adaptive mesh refinement. Water Resour. Res. 2021, 57, e2020WR029354. [Google Scholar] [CrossRef]
  24. Darbandi, M.; Torabi, S.O.; Saadat, M.; Daghighi, Y.; Jarrahbashi, D. A moving-mesh finite-volume method to solve free-surface seepage problem in arbitrary geometries. Int. J. Num. Anal. Meth. Geomech. 2007, 31, 1609–1629. [Google Scholar] [CrossRef]
  25. Huyakorn, P.S.; Pinder, G.F. Computational Methods in Subsurface Flow; Academic Press: Orlando, FL, USA, 1983. [Google Scholar] [CrossRef]
  26. Istok, J. Groundwater Modeling by the Finite Element Method; Water Resources Monograph Series; American Geophysical Union: Washington, DC, USA, 1989; Volume 13. [Google Scholar] [CrossRef]
  27. Cordes, C.; Kinzelbach, W. Continuous groundwater velocity fields and path lines in linear, bilinear, and trilinear finite elements. Water Resour. Res. 1992, 28, 2903–2911. [Google Scholar] [CrossRef]
  28. Mosé, R.; Siegel, P.; Ackerer, P.; Chavent, G. Application of the mixed hybrid finite element approximation in a groundwater flow model: Luxury or necessity? Water Resour. Res. 1994, 30, 3001–3012. [Google Scholar] [CrossRef]
  29. Cordes, C.; Kinzelbach, W. Comment on “Application of the mixed hybrid finite element approximation in a groundwater flow model: Luxury or necessity?” by R. Mosé, P. Siegel, P. Ackerer, and G. Chavent. Water Resour. Res. 1996, 32, 1905–1909. [Google Scholar] [CrossRef]
  30. Ackerer, P.; Mosé, R.; Siegel, P.; Chavent, G. Reply [to “Comment on ‘Application of the mixed hybrid finite element approximation in a groundwater flow model: Luxury or necessity? ’by R. Mosé, P. Siegel, P. Ackerer, and G. Chavent”]. Water Resour. Res. 1996, 32, 1911–1913. [Google Scholar] [CrossRef]
  31. Raviart, P.A.; Thomas, J.M. A mixed finite element method for 2-nd order elliptic problems. In Mathematical Aspects of Finite Element Methods; Galligani, I., Magenes, E., Eds.; Lecture Notes in Mathematics; Springer: Berlin/Heidelberg, Germany, 1977; Volume 606, pp. 292–315. [Google Scholar] [CrossRef]
  32. Durlofsky, L. Accuracy of mixed and control volume finite element approximations to Darcy velocity and related quantities. Water Resour. Res. 1994, 30, 965–973. [Google Scholar] [CrossRef]
  33. Chavent, G.; Roberts, J.E. A unified physical presentation of mixed, mixed-hybrid finite elements and standard finite difference approximations for the determination of velocities in waterflow problems. Adv. Water Resour. 1991, 14, 329–348. [Google Scholar] [CrossRef]
  34. Kaasschieter, E.F. Mixed finite elements for accurate particle tracking in saturated groundwater flow. Adv. Water Resour. 1995, 18, 277–294. [Google Scholar] [CrossRef] [Green Version]
  35. Traverso, L.; Phillips, T.N.; Yang, Y. Mixed finite element methods for groundwater flow in heterogeneous aquifers. Comput. Fluids 2013, 88, 60–80. [Google Scholar] [CrossRef]
  36. Bergamaschi, L.; Putti, M. Mixed finite elements and Newton-type linearizations for the solution of Richards’ equation. Int. J. Num. Meth. Eng. 1999, 45, 1025–1046. [Google Scholar] [CrossRef]
  37. Farthing, M.W.; Kees, C.E.; Miller, C.T. Mixed finite element methods and higher order temporal approximations for variably saturated groundwater flow. Adv. Water Resour. 2003, 26, 373–394. [Google Scholar] [CrossRef]
  38. Fahs, M.; Younes, A.; Lehmann, F. An easy and efficient combination of the mixed finite element method and the method of lines for the resolution of Richards’ equation. Env. Model. Soft. 2009, 24, 1122–1126. [Google Scholar] [CrossRef]
  39. Belfort, B.; Ramasomanana, F.; Younes, A.; Lehmann, F. An efficient lumped mixed hybrid finite element formulation for variably saturated groundwater flow. Vadose Zone J. 2009, 8, 352–362. [Google Scholar] [CrossRef] [Green Version]
  40. Nayagum, D.; Schäfer, G.; Mosé, R. Modelling two-phase incompressible flow in porous media using mixed hybrid and discontinuous finite elements. Comput. Geosci. 2004, 8, 49–73. [Google Scholar] [CrossRef]
  41. Fučík, R.; Klinkovský, J.; Solovský, J.; Oberhuber, T.; Mikyška, J. Multidimensional mixed-hybrid finite element method for compositional two-phase flow in heterogeneous porous media and its parallel implementation on GPU. Comput. Phys. Commun. 2019, 238, 165–180. [Google Scholar] [CrossRef]
  42. Durlofsky, L.J. Upscaling and gridding of fine scale geological models for flow simulation. In Proceedings of the 8th International Forum on Reservoir Simulation Iles Borromees, Stresa, Italy, 20–24 June 2005; Volume 2024. [Google Scholar]
  43. Dagan, G. Second order linearized theory of free surface flow in porous media. La Houille Blanche 1964, 8, 901–910. [Google Scholar] [CrossRef] [Green Version]
  44. Neuman, S.P. Theory of flow in unconfined aquifers considering delayed response of the water table. Water Resour. Res. 1972, 8, 1031–1045. [Google Scholar] [CrossRef]
  45. Sousa, M.R.; Jones, J.P.; Frind, E.O.; Rudolph, D.L. A simple method to assess unsaturated zone time lag in the travel time from ground surface to receptor. J. Contam. Hydrol. 2013, 144, 138–151. [Google Scholar] [CrossRef]
  46. Wossenyeleh, B.K.; Verbeiren, B.; Diels, J.; Huysmans, M. Vadose zone lag time effect on groundwater drought in a temperate climate. Water 2020, 12, 2123. [Google Scholar] [CrossRef]
  47. Teramoto, E.H.; Crioni, P.L.B.; Chang, H.K. Daily time series of groundwater recharge derived from temporal variation of water level. Sustain. Water Resour. Manag. 2021, 7, 67. [Google Scholar] [CrossRef]
  48. Lee, L.J.E.; Lawrence, D.S.L.; Price, M. Analysis of water-level response to rainfall and implications for recharge pathways in the Chalk aquifer, SE England. J. Hydrol. 2006, 330, 604–620. [Google Scholar] [CrossRef]
  49. Szymkiewicz, A.; Gumuła-Kawecka, A.; Potrykus, D.; Jaworska-Szulc, B.; Pruszkowska-Caceres, M.; Gorczewska-Langner, W. Estimation of conservative contaminant travel time through vadose zone based on transient and steady flow approaches. Water 2018, 10, 1417. [Google Scholar] [CrossRef] [Green Version]
  50. Saad, Y. Iterative Methods for Sparse Linear Systems; SIAM: Philadelphia, PA, USA, 2003. [Google Scholar] [CrossRef]
  51. Larabi, A.; De Smedt, F. Solving three-dimensional hexahedral finite element groundwater models by preconditioned conjugate gradient methods. Water Resour. Res. 1994, 30, 509–521. [Google Scholar] [CrossRef]
  52. Papadopulos, I.S. Nonsteady flow to a well in an infinite anisotropic aquifer. In Proceedings of the Dubrovnik Symposium on the Hydrology of Fractured Rocks, Dubrovnik, Croatia, 7–14 October 1965; International Association of Scientific Hydrology: Wallingford, UK, 1965; pp. 21–31. [Google Scholar]
  53. Theis, C.V. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage. AGU Trans. 1935, 16, 519–524. [Google Scholar] [CrossRef]
  54. Pollock, D.W. Semi-analytical computation of path lines for finite difference models. Groundwater 1988, 26, 743–750. [Google Scholar] [CrossRef]
  55. Lacy, S.J.; Prevost, J.H. Flow through porous media: A procedure for locating the free surface. Int. J. Numer. Anal. Met. Geomech. 1987, 11, 508–601. [Google Scholar] [CrossRef]
  56. Bardet, J.-P.; Tobita, T. A practical method for solving free-surface seepage problems. Comput. Geotech. 2002, 29, 451–475. [Google Scholar] [CrossRef]
  57. Sbai, M.A.; Larabi, A. On solving groundwater flow and transport models with algebraic multigrid preconditioning. Groundwater 2021, 59, 100–108. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Illustration of different approaches for modeling unconfined groundwater flow in subsurface aquifers. (A) An example of flow processes occurring in an unconfined aquifer such as recharge, abstraction, and aquifer-river interactions. (B) A discrete variably saturated groundwater flow model seeks to determine the saturation state in all cells; the water table shape and position are post-processed from simulated pressure heads. (C) A discrete unconfined flow model neglecting processes in the unsaturated zone using a fixed mesh. The numerical model identifies wet, partially wet, and dry cells. The water table is not represented explicitly and is post-processed from simulated heads. (D) A discrete unconfined flow model neglecting processes in the unsaturated zone using a moving mesh. The numerical model works exclusively with fully saturated cells. The water table is explicitly represented by the upper mesh slice and is, therefore, not post-processed.
Figure 1. Illustration of different approaches for modeling unconfined groundwater flow in subsurface aquifers. (A) An example of flow processes occurring in an unconfined aquifer such as recharge, abstraction, and aquifer-river interactions. (B) A discrete variably saturated groundwater flow model seeks to determine the saturation state in all cells; the water table shape and position are post-processed from simulated pressure heads. (C) A discrete unconfined flow model neglecting processes in the unsaturated zone using a fixed mesh. The numerical model identifies wet, partially wet, and dry cells. The water table is not represented explicitly and is post-processed from simulated heads. (D) A discrete unconfined flow model neglecting processes in the unsaturated zone using a moving mesh. The numerical model works exclusively with fully saturated cells. The water table is explicitly represented by the upper mesh slice and is, therefore, not post-processed.
Water 15 01177 g001
Figure 2. (a) Plan view of the unstructured mesh used for simulating the hydraulic head drawdown in an anisotropic aquifer. (b) Enlarged view around the first level of mesh refinement within a 10 m radius from the domain center. The two piezometers are shown. (c) Enlarged view around the second level of mesh refinement within a 1 m radius from the domain center. The well is approximated with a regular pentagon and has a radius of 0.1 m.
Figure 2. (a) Plan view of the unstructured mesh used for simulating the hydraulic head drawdown in an anisotropic aquifer. (b) Enlarged view around the first level of mesh refinement within a 10 m radius from the domain center. The two piezometers are shown. (c) Enlarged view around the second level of mesh refinement within a 1 m radius from the domain center. The well is approximated with a regular pentagon and has a radius of 0.1 m.
Water 15 01177 g002
Figure 3. The hydraulic head drawdown in an anisotropic aquifer: comparison of the Papadopulos analytical solution with results of the MHFEM model at piezometers P1 and P2 with (a) a diagonal transmissivity tensor, and (b) a full symmetric transmissivity tensor.
Figure 3. The hydraulic head drawdown in an anisotropic aquifer: comparison of the Papadopulos analytical solution with results of the MHFEM model at piezometers P1 and P2 with (a) a diagonal transmissivity tensor, and (b) a full symmetric transmissivity tensor.
Water 15 01177 g003
Figure 4. Conceptual model of the second test problem. Groundwater is flowing between two fixed head boundaries, through a heterogeneous porous medium. The hydraulic conductivity distribution spans four orders of magnitude.
Figure 4. Conceptual model of the second test problem. Groundwater is flowing between two fixed head boundaries, through a heterogeneous porous medium. The hydraulic conductivity distribution spans four orders of magnitude.
Water 15 01177 g004
Figure 5. Simulated path lines (black lines) with both the MHFEM code HydroTec© and the finite-difference-based MODFLOW/MODPATH codes are presented in the first (ac) and second columns (df), respectively, and for the 10 × 10, 20 × 20, and 30 × 30 grid resolutions starting from the top. The reference path lines (red lines) are computed using a mixed approximation with a 320 × 320 grid resolution. The MHFEM approximation yields a more accurate solution at a 20 × 20 resolution (b) than a finite difference approximation at a higher resolution (f).
Figure 5. Simulated path lines (black lines) with both the MHFEM code HydroTec© and the finite-difference-based MODFLOW/MODPATH codes are presented in the first (ac) and second columns (df), respectively, and for the 10 × 10, 20 × 20, and 30 × 30 grid resolutions starting from the top. The reference path lines (red lines) are computed using a mixed approximation with a 320 × 320 grid resolution. The MHFEM approximation yields a more accurate solution at a 20 × 20 resolution (b) than a finite difference approximation at a higher resolution (f).
Water 15 01177 g005
Figure 6. Calculated relative errors in simulated residence times of two particles released at (25, 100) (a) and (95, 100) (b) using the MHFEM and FDM approximations.
Figure 6. Calculated relative errors in simulated residence times of two particles released at (25, 100) (a) and (95, 100) (b) using the MHFEM and FDM approximations.
Water 15 01177 g006
Figure 7. Simulated steady-state water table position using the adaptive mesh MHFEM technique for the third test problem. All horizontal lines of the two-dimensional mesh were allowed to move during the iterative procedure. The water table is favorably compared with numerical solutions provided by other authors [55,56].
Figure 7. Simulated steady-state water table position using the adaptive mesh MHFEM technique for the third test problem. All horizontal lines of the two-dimensional mesh were allowed to move during the iterative procedure. The water table is favorably compared with numerical solutions provided by other authors [55,56].
Water 15 01177 g007
Figure 8. Simulated steady-state water table position and groundwater heads using the adaptive mesh MHFEM technique for the fourth test problem. All horizontal lines of the two-dimensional mesh were allowed to move during the iterative procedure (vertical aspect ratio = 80).
Figure 8. Simulated steady-state water table position and groundwater heads using the adaptive mesh MHFEM technique for the fourth test problem. All horizontal lines of the two-dimensional mesh were allowed to move during the iterative procedure (vertical aspect ratio = 80).
Water 15 01177 g008
Figure 9. Performance of the adaptive mesh iteration for the fourth test problem. The water table error decreases with a rate greater than a linear trend. The bars show the number of linear iterations at each outer iteration. Eight mesh iteration steps were sufficient to achieve accuracy with a 10−3 m closure criterion.
Figure 9. Performance of the adaptive mesh iteration for the fourth test problem. The water table error decreases with a rate greater than a linear trend. The bars show the number of linear iterations at each outer iteration. Eight mesh iteration steps were sufficient to achieve accuracy with a 10−3 m closure criterion.
Water 15 01177 g009
Figure 10. Comparison between the water table simulated with the adaptive MHFEM approximation and the post-processed FDM approximation. There is a good match for the well drawdown. The observed discrepancies far from the well are mainly attributed to post-processing interpolation errors.
Figure 10. Comparison between the water table simulated with the adaptive MHFEM approximation and the post-processed FDM approximation. There is a good match for the well drawdown. The observed discrepancies far from the well are mainly attributed to post-processing interpolation errors.
Water 15 01177 g010
Figure 11. Simulated steady-state water table position and groundwater heads using the adaptive mesh MHFEM technique for the three-dimensional problem. The cutaway view shows the pumping well’s depression cone and the head contours at the aquifer base (vertical aspect ratio = 50).
Figure 11. Simulated steady-state water table position and groundwater heads using the adaptive mesh MHFEM technique for the three-dimensional problem. The cutaway view shows the pumping well’s depression cone and the head contours at the aquifer base (vertical aspect ratio = 50).
Water 15 01177 g011
Figure 12. Performance of the adaptive mesh iteration for the three-dimensional problem. The water table error decreases with a rate greater than a linear trend. The bars show the number of linear iterations at each outer iteration. Thirteen mesh iterations were sufficient to achieve accuracy with a 10−2 m closure criterion.
Figure 12. Performance of the adaptive mesh iteration for the three-dimensional problem. The water table error decreases with a rate greater than a linear trend. The bars show the number of linear iterations at each outer iteration. Thirteen mesh iterations were sufficient to achieve accuracy with a 10−2 m closure criterion.
Water 15 01177 g012
Figure 13. Comparison between the water table profile (parallel to X-axis and crossing the well position) simulated with the adaptive MHFEM approximation and the post-processed FDM approximation for the three-dimensional problem. The maximum difference occurs at the well due to the coarse vertical discretization of the FDM fixed mesh.
Figure 13. Comparison between the water table profile (parallel to X-axis and crossing the well position) simulated with the adaptive MHFEM approximation and the post-processed FDM approximation for the three-dimensional problem. The maximum difference occurs at the well due to the coarse vertical discretization of the FDM fixed mesh.
Water 15 01177 g013
Table 1. Physical parameters used in the pumping test problem. Case 1 and 2 correspond to diagonal and fully symmetric transmissivity tensors, respectively.
Table 1. Physical parameters used in the pumping test problem. Case 1 and 2 correspond to diagonal and fully symmetric transmissivity tensors, respectively.
ParameterValue—Case 1Value—Case 2
Transmissivity   along   X ,   T x x m2.s−12.93 × 10−52.57 × 10−5
Transmissivity   along   Y ,   T y y m2.s−11.47 × 10−51.84 × 10−5
Transmissivity   along   XY ,   T x y = T y x m2.s−106.36 × 10−6
Specific   storage ,   S p m−11.957 × 10−31.957 × 10−3
Pumping   flow   rate ,   Q w m3.s−15 × 10−35 × 10−3
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Sbai, M.A.; Larabi, A. A Deforming Mixed-Hybrid Finite Element Model for Robust Groundwater Flow Simulation in 3D Unconfined Aquifers with Unstructured Layered Grids. Water 2023, 15, 1177. https://doi.org/10.3390/w15061177

AMA Style

Sbai MA, Larabi A. A Deforming Mixed-Hybrid Finite Element Model for Robust Groundwater Flow Simulation in 3D Unconfined Aquifers with Unstructured Layered Grids. Water. 2023; 15(6):1177. https://doi.org/10.3390/w15061177

Chicago/Turabian Style

Sbai, Mohammed Adil, and Abdelkader Larabi. 2023. "A Deforming Mixed-Hybrid Finite Element Model for Robust Groundwater Flow Simulation in 3D Unconfined Aquifers with Unstructured Layered Grids" Water 15, no. 6: 1177. https://doi.org/10.3390/w15061177

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop