In order to understand and reproduce the movement of groundwater, it is necessary to solve the mathematical equations that represent groundwater flow.

Mathematical equations governing flow are partial differential equations whose solution requires properties of the medium (e.g. hydraulic conductivity and storage coefficients) as well as boundary conditions (e.g. recharge, well flow, impermeable boundaries and bodies of water).

These equations have varying degrees of complexity and can represent homogeneous or heterogeneous media, saturated and unsaturated zones, porous media or fractured media, uniform or variable fluid density and viscosity, among others.

Likewise, the solution of the equations can vary in complexity. For homogeneous media with simple boundary conditions, analytical solutions can be used. Classic examples of analytical solutions are the Theis, Hantush or Neuman solutions, widely used in the interpretation of pumping tests in confined, semi-confined and free aquifers, respectively.

On the other hand, for systems with more realistic geology, aquifer thicknesses and boundary conditions, it is necessary to use numerical solutions, where the domain is compartmentalized into small pieces (called cells, elements or control volumes depending on the method used.

Within each volume, the property is uniform and the method makes a mass balance of all the inputs and outputs of each neighboring volume. The most commonly used numerical methods for groundwater are the Finite Difference Method, used by the USGS program MODFLOW, the Finite Volume Method, also an option in the more recent MODFLOW USG or MODFLOW 6 versions, and the Finite Element Method, used by the DHI program FEFLOW.

The aim of these methods is to obtain an approximate solution to the governing equations through temporal and spatial discretization, where the primary variable is obtained at points defined in space and time. In the finite element method, each element has nodes where the hydraulic load is calculated and Gauss points where a secondary variable can be calculated. In the context of flow in porous media, this secondary variable corresponds to the Darcy flow or Darcy velocity.

In FEFLOW software, we can create 2D and 3D numerical meshes.

In 3D, you can use prismatic triangular elements. In this case, the mesh is called structured, with several layers subdivided into triangles connected by vertical lines between the different layers, or tetrahedral elements. When using tetrahedral elements, there is no longer the concept of a layer and the mesh is called unstructured.

The use of tetrahedral elements allows for a more faithful representation of more complex geologies or structures, such as a non-vertical intrusive fault or dyke, or the filters and elements of an earth dam. However, there is an increase in the complexity of constructing and editing the mesh and its properties, given the current restrictions of FEFLOW’s graphical interface, which is essentially designed to deal with layered systems (slices and layers in FEFLOW terminology).

Unstructured meshes offer the flexibility needed to capture geological and hydrogeological heterogeneity, as well as the complexity of the interfaces between different layers, geological formations and important components of geotechnical structures (dykes, dams, piles, etc.). The use of tetrahedral elements in such meshes is particularly advantageous, as these elements are capable of representing arbitrary volumes, making them suitable for modeling irregular geometry with greater precision. In addition, the adaptability of unstructured meshes allows for greater resolution in specific areas of interest, such as wells, sources of contamination or areas of high hydraulic conductivity (filters, drainage mats, etc.). In this way, it is possible to concentrate computational resources where they are really needed, optimizing the performance of hydrogeological models and reducing computational costs.

It is important to note that the implementation of unstructured meshes composed of tetrahedral elements requires care in the generation and optimization of the mesh, as well as the appropriate choice of numerical methods used in hydrogeological models. It is essential to guarantee the quality of the mesh, avoiding poor quality elements that could lead to inaccurate results or numerical instabilities.

The parameters that control the quality of tetrahedral elements include the Maximum dihedral angles (maximum 150 degrees) and the Aspect ratio of tetrahedra (less than 2.5).

One application implemented within the Numerical Modeling Group at Water Services and Technologies is the use of unstructured meshes for hydrogeological models applied to the study of geotechnical structures. This was adopted because the position of the water level in the reservoirs and in the massifs of the geotechnical structures, through hydrogeological modeling, is of paramount importance to the client since these variables directly influence the safety levels of the geotechnical structure.

Therefore, accurate representation of the complex geometry of the massifs and the components that control flow (mainly internal drainage) is fundamental to understanding the behavior of groundwater flow and assessing the hydrogeological risks involved.

The choice of tetrahedral elements to compose unstructured meshes is particularly advantageous in hydrogeological modeling of dams, as these elements are capable of representing arbitrary volumes in three dimensions. This allows the meshes to adequately adjust to the irregular geometry of the subsoil, including the challenges of representing dam walls and their associated components with high precision, especially those elements that are more conductive (internal drainage).

The biggest challenge, however, is how to generate a numerical tetrahedral mesh that complies with FEFLOW’s constraints in the least time-consuming way possible. Some programs can be used outside of FEFLOW, such as Autodesk CIvil3D, ANSYS and others. At Water Services and Technologies, we use Bentley’s PLAXIS Designer. With software designed for geotechnical simulations, it offers a range of facilities for constructing the mesh for this purpose. Finally, we found that the most efficient way to develop and subsequently manipulate the mesh and properties is to combine mesh types in the same model. We started by building a structured mesh and, in areas where a more detailed definition of the structures is required, we added an unstructured mesh.

Below is an example of a representation of an earth dam with its components perfectly represented in the unstructured mesh portion.