MODFLOW–USG (UnStructured Grids) is capable of simulating groundwater flow using a wide range of grid geometries. All of the grid information is characterized by the connectivity, connection lengths, connection flow areas, and cell areas and volumes. This information is provided to MODFLOW–USG in the unstructured discretization file (.DISU). For accurate flux calculations, the Control Volume Finite Difference (CVFD) method requires certain geometrical properties pertaining to the cell connections. For most grid types, this means that the line connecting the centers of two cells should bisect the shared edge at right angles. This condition is satisfied by equilateral triangles, rectangles, regular polygons, as shown below. (Panday et. al., 2013)



Voronoi Tesselations: An Ideal Grid Geometry for MODFLOW-USG

Another ideal geometry that satisfies the MODFLOW-USG CVFD geometry constraints is a Voronoi Tesselation. When you generate a Delaunay Triangulation from a set of points, a Voronoi Tesselation is created as the “dual” to this triangulation. (A more thorough explanation of this concept is available in the References section below). The process works as follows: You start with a set of shapes that represent your conceptual boundary conditions; these include points for pumping/observation well locations, and polylines and polygons representing rivers/streams, surface water bodies, or other boundary conditions (or geological features). The vertices from polylines and polygons are extracted and used as points for the grid design. A simple example of the conceptual representations is shown below: constant heads are defined on the east and west edges, one pumping well is located in the northeast corner, another is located in the southwest corner, and a meandering river passes through the middle of the model domain.


These points are inputs (control points) for a Delaunay triangulation; in finite element modeling, specifically FEFLOW, these objects would define the super-element mesh. A resulting triangulation is shown below, with refinement and smoothing around the wells and boundaries.



Delaunay triangulation with refinement around around the constant heads (on west and east boundaries),
the two pumping wells, and the river.

During the triangulation, the “by-product” Voronoi Tesselation is created. The points/vertices from the conceptual model objects form the node generators for the Voronoi Tesselation. The Voronoi Tesselation consists of a set of polygons that can be used as an unstructured grid for MODFLOW-USG.



MODFLOW-USG grid with refinement around around the constant heads (on west and east boundaries),
the two pumping wells, and the river.


A main feature of the Voronoi Tesselation is that the lines connecting the node generators are always perpendicular to the cell edges, making it ideal for MODFLOW-USG and CVFD requirements.


Benefits of Using Voronoi Tesselation Grids

Voronoi Tessellations provide similar capabilities as Finite Element (FE) meshes: efficient refinement around wells and other boundary condition objects, ideal representation of variable stratigraphy and pinchout layers, and suitable representation of irregularly shaped model domains (i.e. non-rectangular shaped model boundary). You can have higher accuracy “where it counts” and coarser cell size outside the areas of interest.

MODFLOW USG-Benchmark Tests in Visual MODFLOW Flex

We benchmarked several Voronoi Tesselation Unstructured Grids in Visual MODFLOW Flex. The full results will be presented at the MODFLOW 2013 Conference in June. The simple benchmark model is a 2D unconfined aquifer with the following characteristics:

  • approximately 20 km * 19 km
  • flat top and bottom elevation with an aquifer thickness of 50 m
  • hydraulic conductivity of 1 m/d (for Kx, Ky, Kx)
  • constant heads on the east and west edges of the model
  • two pumping wells, each with an extraction rate of 500 m3/day
  • a meandering river through the middle (with river stage=45 m, riverbed bottom=42 m, riverbed thickness=1 m, width=10 m, river bed conductivity=1 m/day)

For this conceptual model, we generated a Delaunay Triangulation as a mesh for FEFLOW and the corresponding Voronoi Tesselation as an unstructured grid for MODFLOW-USG. We compared the MODFLOW-USG grid to the equivalent FEFLOW model, and also compared both models to a finely-spaced MODFLOW-2005 model (with a uniform grid spacing of 25 m, resulting in a total of ~600,000 cells). FEFLOW calculates the heads/fluxes at the nodes, whereas MODFLOW-USG calculates heads/fluxes at the Voronoi cell generators. Since the generators of the MODFLOW-USG cells exactly coincide with the nodes of the finite element mesh, this makes it possible to do a direct comparison between results from MODFLOW-USG with results from FEFLOW. The Ghost Node Correction (GNC) package of MODFLOW-USG was not used for these examples. (for more details on the GNC package, please refer to p.24 of the MODFLOW-USG User Manual)

Three different MODFLOW-USG grids were created, with varying degrees of refinement:
• Coarse Grid: approx. 200 m USG cell spacing around the wells and boundaries (total of 1495 nodes)
• Medium Grid: approx. 100 m USG cell spacing around the wells and boundaries (total of 4571 nodes)
• Fine Grid: approx. 25 m USG cell spacing around the wells and boundaries (total of 19048 nodes)

The xMD solver was used for MODFLOW-USG with the default parameters; the PCG solver was used for FEFLOW with the default settings (with FEFLOW v.5.4)


The “Fine” MODFLOW-USG grid yielded the best results when looking at the head observed at Pumping Well 2 (pw2). A head value of 24.5 m was observed which matches closely with the FEFLOW results (24.67 m) and the MODFLOW-2005 results (24.55 m). For the coarse and medium sized MODFLOW-USG grids, the observed head at PW2 was higher (approx. 32 m), which is to be expected given the larger cell sizes.



Results from MODFLOW-USG run with the coarse grid. Cells are colored by the calculated head value.



Results from the FEFLOW run with coarse mesh, with continous color map of heads.


The mass balance for the “Fine USG” grid closely matches the results observed with the “Fine Mesh” FEFLOW run and the MODFLOW-2005 run (0.00% mass balance discrepency in all cases).



The run times are very quick for the fine grid for both MODFLOW-USG and FEFLOW (less than 1.0 second). The MODFLOW-2005 runtime was approximately 60 seconds (which is to be expected given that the MODFLOW-2005 model required much more cells in order to obtain the same level of cell resolution around the pumping wells).



The example below shows the results from the “Fine Grid” MODFLOW-USG run, with detailed refinement around the pumping wells.



Results from the MODFLOW-USG “Fine Grid” with zoomed in windows showing the
detailed refinement around the pumping wells.


Future Work

Future work will include testing with more complex 3D problems, including transient boundary conditions, and benchmarking accuracy and performance to MODFLOW-2005 and FEFLOW. We will also implement Ghost Node Correction (GNC) for the above grids to improve accuracy.



We would like to thank Sorab Panday for his technical guidance and troubleshooting support during this exercise.



MODFLOW–USG Version 1: An Unstructured Grid Version of MODFLOW for Simulating Groundwater Flow and Tightly Coupled Processes Using a Control Volume Finite-Difference Formulation. USGS Techniques and Methods 6–A45.  More details:

Voronoi Tesselation:

Delaunay Triangulation: