This example illustrates how to setup a very simple unstructured groundwater model using the imod package and associated packages.

In overview, we’ll set the following steps:

  • Create a triangular mesh for a disk geometry.

  • Create the xugrid UgridDataArrays containg the MODFLOW6 parameters.

  • Feed these arrays into the imod mf6 classes.

  • Write to modflow6 files.

  • Run the model.

  • Open the results back into UgridDataArrays.

  • Visualize the results.

We’ll start with the following imports:

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import xugrid as xu
from pandas import isnull

import imod

Create a mesh#

The first steps consists of generating a mesh. In this example, we’ll use data included with iMOD Python for a circular mesh. Note that this is a Ugrid2D object. For more information on working with unstructured grids see the Xugrid documentation

grid =

Dimensions:            (mesh2d_nFaces: 216, mesh2d_nMax_face_nodes: 3,
                        mesh2d_nEdges: 342, two: 2, mesh2d_nNodes: 127)
    mesh2d_node_x      (mesh2d_nNodes) float64 0.0 6.123e-14 ... 331.1 157.5
    mesh2d_node_y      (mesh2d_nNodes) float64 0.0 1e+03 500.0 ... 764.7 818.3
Dimensions without coordinates: mesh2d_nFaces, mesh2d_nMax_face_nodes,
                                mesh2d_nEdges, two, mesh2d_nNodes
Data variables:
    mesh2d             int64 0
    mesh2d_face_nodes  (mesh2d_nFaces, mesh2d_nMax_face_nodes) int32 0 7 ... 1
    mesh2d_edge_nodes  (mesh2d_nEdges, two) int64 0 7 0 17 0 ... 124 126 125 126
    Conventions:  CF-1.9 UGRID-1.0

We can plot this object as follows:

fig, ax = plt.subplots()
xu.plot.line(grid, ax=ax)

Create UgridDataArray#

Now that we have defined the grid, we can start defining the model parameter data.

Our goal here is to define a steady-state model with:

  • Uniform conductivities of 1.0 m/d;

  • Two layers of 5.0 m thick;

  • Uniform recharge of 0.001 m/d on the top layer;

  • Constant heads of 1.0 m along the exterior edges of the mesh.

From these boundary conditions, we would expect circular mounding of the groundwater; with small flows in the center and larger flows as the recharge accumulates while the groundwater flows towards the exterior boundary.

nface = grid.n_face
nlayer = 2

idomain = xu.UgridDataArray(
        np.ones((nlayer, nface), dtype=np.int32),
        coords={"layer": [1, 2]},
        dims=["layer", grid.face_dimension],
icelltype = xu.full_like(idomain, 0)
k = xu.full_like(idomain, 1.0, dtype=float)
k33 = k.copy()
rch_rate = xu.full_like(idomain.sel(layer=1), 0.001, dtype=float)
bottom = idomain * xr.DataArray([5.0, 0.0], dims=["layer"])

All the data above have been constants over the grid. For the constant head boundary, we’d like to only set values on the external border. We can py:method:xugrid.UgridDataset.binary_dilation to easily find these cells:

chd_location = xu.zeros_like(idomain.sel(layer=2), dtype=bool).ugrid.binary_dilation(
constant_head = xu.full_like(idomain.sel(layer=2), 1.0, dtype=float).where(chd_location)

fig, ax = plt.subplots()
xu.plot.line(grid, ax=ax, color="black")
layer = 2

Write the model#

The first step is to define an empty model, the parameters and boundary conditions are added in the form of the familiar MODFLOW packages.

gwf_model = imod.mf6.GroundwaterFlowModel()
gwf_model["disv"] = imod.mf6.VerticesDiscretization(
    top=10.0, bottom=bottom, idomain=idomain
gwf_model["chd"] = imod.mf6.ConstantHead(
    constant_head, print_input=True, print_flows=True, save_flows=True
gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0)
gwf_model["npf"] = imod.mf6.NodePropertyFlow(
gwf_model["sto"] = imod.mf6.SpecificStorage(
gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all")
gwf_model["rch"] = imod.mf6.Recharge(rch_rate)

simulation = imod.mf6.Modflow6Simulation("circle")
simulation["GWF_1"] = gwf_model
simulation["solver"] = imod.mf6.Solution(
simulation.create_time_discretization(additional_times=["2000-01-01", "2000-01-02"])

We’ll create a new directory in which we will write and run the model.

modeldir = imod.util.temporary_directory()

Run the model#


The following lines assume the mf6 executable is available on your PATH. The Modflow 6 examples introduction shortly describes how to add it to yours.

Open the results#

First, we’ll open the heads (.hds) file.

head = simulation.open_head()

<xarray.DataArray 'head' (time: 1, layer: 2, mesh2d_nFaces: 216)>
dask.array<stack, shape=(1, 2, 216), dtype=float64, chunksize=(1, 2, 216), chunktype=numpy.ndarray>
  * layer          (layer) int64 1 2
  * time           (time) float64 1.0
  * mesh2d_nFaces  (mesh2d_nFaces) int64 0 1 2 3 4 5 ... 210 211 212 213 214 215

For a DISV MODFLOW6 model, the heads are returned as a UgridDataArray. While all layers are timesteps are available, they are only loaded into memory as needed.

We may also open the cell-by-cell flows (.cbc) file.

cbc = simulation.open_flow_budget()

dict_keys(['flow-horizontal-face', 'flow-horizontal-face-x', 'flow-horizontal-face-y', 'flow-lower-face', 'chd'])

The flows are returned as a dictionary of UgridDataArrays. This dictionary contains all entries that are stored in the CBC file, but like for the heads file the data are only loaded into memory when needed.

The horizontal flows are stored on the edges of the UgridDataArray topology. The other flows are generally stored on the faces; this includes the flow-lower-face.

We’ll create a dataset for the horizontal flows for further analysis.

ds = xu.UgridDataset(grids=grid)
ds["u"] = cbc["flow-horizontal-face-x"]
ds["v"] = cbc["flow-horizontal-face-y"]

Visualize the results#

We can quickly and easily visualize the output with the plotting functions provided by xarray and xugrid. We’ll add some some edge coordinates to the dataset so that they can be used to place the arrows in the quiver plot.

ds = ds.ugrid.assign_edge_coords()
fig, ax = plt.subplots()
head.isel(time=0, layer=0).compute().ugrid.plot(ax=ax)
ds.isel(time=0, layer=0).plot.quiver(
    x="mesh2d_edge_x", y="mesh2d_edge_y", u="u", v="v", color="white"
layer = 1, time = 1.0

As would be expected from our model input, we observe circular groundwater mounding and increasing flows as we move from the center to the exterior.

Slice the model domain#

We may also quickly setup a smaller model. We’ll select half of the original domain. To set up the boundary conditions on the clipped edges you can provide a states_for_boundary dictionary. In this case we add the head values of the computed full domain simulation as the clipped boundary values

states_for_boundary = {
    "GWF_1": head.compute(),

half_simulation = simulation.clip_box(
    x_max=0.0, states_for_boundary=states_for_boundary

Let’s run the model, read the results, and visualize.

modeldir = imod.util.temporary_directory()
head = half_simulation.open_head()

Let’s add constant head boundaries together and plot them

half_simulation_constant_head = half_simulation["GWF_1"]["chd"]["head"]

clipped_half_simulation_constant_head = (

all_boundaries_constant_head = half_simulation_constant_head.where(
    ~isnull(half_simulation_constant_head), clipped_half_simulation_constant_head

# plot boundary conditions
fig, ax = plt.subplots()
layer = 2, time = 1.0

plot computed heads

fig, ax = plt.subplots()
head.isel(time=0, layer=0).compute().ugrid.plot(ax=ax)
layer = 1, time = 1.0

Total running time of the script: (0 minutes 3.258 seconds)

Gallery generated by Sphinx-Gallery