**October 14, 2016**

Those that are new to groundwater modeling are often overwhelmed with terms such as finite differences, finite elements, grid discretization and triangulation. Once you gain a basic understanding of the terms and the processes involved in groundwater modeling, you are often faced with the ultimate question “Which numerical method should I use for my project?” The answer is not often so much what you should use but often is based on what method or software package is readily available. Like most people, groundwater modelers want to stick with what they know. There are, however, three main numerical methods that should be considered: Finite Difference, Finite Element and Finite Volume.

**Finite Difference**

MODFLOW uses the finite differences method to solve the groundwater flow equation. The grid is created using structured, rectilinear (rectangular) grids. This method is widely used and accepted by many governing bodies. It is the most popular groundwater simulator available, is open source and is well documented. The finite difference solution is easy to understand and calculate, the solutions are mass-conservative and there are several numerical extensions available such as PEST, transport, particle tracking and zone budget. However, the finite difference method is not without weakness. Although grids are easy to create, they cannot be efficiently refined around areas of interest, such as wells and model boundaries. This can be partially addressed by the MODFLOW-LGR (Local Grid Refinement) extension.

Complex geology is difficult to represent particularly if there are discontinuous or pinched out layers; when there are steep gradients in the stratigraphy, this can result in disconnected cells which causes problems with running the model. MODFLOW requires model layers to be continuous across the entire model domain, making it a challenge to model smaller formations of materials (lenses) within a larger body, without adding a large number of unnecessary cells. MODFLOW provides an approximate representation of variable anisotropy in model parameters (e.g. it only uses the principle components of the hydraulic conductivity tensor; Kxx, Kyy and Kzz), and therefore is limited for modeling faults and fractures.

**Finite Element**

The finite element method can also be used to solve the groundwater flow equation. There are several codes available; FEFLOW being the most popular. FEFLOW uses a finite element (triangular) mesh to represent the model domain. The use of triangles allows for a more efficient refinement around wells and boundaries. The triangular mesh can more easily adapt to variable stratigraphy such as sloping or pinch outs, and allows for versatile discretization of non-rectangular model domains. Finite element methods provide a better representation of anisotropy (fully represents conductivity tensor – each triangular element is given its own coordinate system, whereas MODFLOW requires conductivity to be perpendicular to the faces of the finite difference cells (horizontally rectangular)).

No method is perfect; with finite elements the local mass conservation is not guaranteed. There can be discontinuous velocities at the element boundaries which can make it difficult to determine unique pathlines.

FEFLOW is similar to MODFLOW in that model layers must be continuous across the entire model domain. In order to represent smaller formations within larger bodies, an entire model layer must be inserted with a minimum layer thickness enforced in the discontinuous regions.

**Finite Volume**

The finite volume method is now being introduced into groundwater modeling through the MODFLOW-USG (UnStructured Grids) code. MODFLOW-USG was released by the USGS in May 2013 and follows a Control Volume Finite Difference (CVFD) formulation in which a cell can be connected to an arbitrary number of adjacent cells. This allows infinite possibilities for the cell geometry. The model domain can be discretized using triangular, rectangular or Voronoi polygons and allows for effective vertical gridding. This means that the grid can be refined locally around areas of interest, such as wells and boundaries, without adding extra cells outside the areas of interest. MODFLOW-USG does not require layers to be continuous across the model domain, allowing for efficient representation of discontinuous layers (lenses, perched aquifers, etc). The CVFD solution remains mass conservative and efficient.

MODFLOW-USG follows the traditional MODFLOW notation for most packages, reducing the learning curve for those that are already familiar with MODFLOW. Many of the packages follow similar formats, with just a few small adjustments to consider the differences in cell geometry and numbering. The code is open source and readily available.

For both Finite Difference and Finite Element methods, the horizontal discretization (number of cells/elements) must be the same in all model layers; this is quite inefficient, because it means that you end up with a large number of cells/elements in layers outside the area of interest; with Finite Volume, MODFLOW-USG, each layer can have its own discretization; this means that you can refine upper layers around rivers and streams, and have coarser refinement in lower layers; likewise, you can refine around wells screens in lower layers, and have less refinement in layers above. The result with MODFLOW-USG is fewer cells, faster solutions, while maintaining high accuracy.

Essentially, MODFLOW-USG contains the best of the Finite Difference and Finite Element codes. The shortcomings of this solution method do not seem to be technical but related to its ‘newness’. Groundwater modelers will need to learn the ‘ins’ and ‘outs’ of this new code to see how it holds up with the tried and true methods. In addition, avid modelers will need to wait for this code to incorporate add-ons such as transport, particle tracking, and the UZF package. At the end of the day, it is a matter of acceptance; will MODFLOW- USG be adopted by the groundwater modeling community as the next great thing?