import imod
import xarray as xr
import xugrid as xu
In this example, we’ll work with an unstructured grid. We’ll create a very simple unstructured model of the Netherlands from scratch. We’ll use a digital elevation model (DEM) to set as drain level to simulate overland flow. To this we add constant recharge.
Let’s start with importing the required packages. These are iMOD Python, xarray, and xugrid.
Load the data
This data is distributed with the xugrid package. You need an internet connection to download the data.
= xu.data.elevation_nl()
elevation
elevation
Downloading file 'elevation_nl.nc' from 'https://github.com/deltares/xugrid/raw/main/data/elevation_nl.nc' to '/home/runner/.cache/xugrid'.
<xarray.DataArray 'elevation' (mesh2d_nFaces: 5248)> Size: 21kB [5248 values with dtype=float32] Coordinates: mesh2d_face_x (mesh2d_nFaces) float64 42kB ... mesh2d_face_y (mesh2d_nFaces) float64 42kB ... * mesh2d_nFaces (mesh2d_nFaces) int64 42kB 0 1 2 3 4 ... 5244 5245 5246 5247 Attributes: unit: m NAP
Let’s plot the data
=-20, vmax=90, cmap="terrain") elevation.ugrid.plot(vmin
Template grid
It is convenient to create a template grid, which has the right shape, dimensions, and coordinates. This can be used to construct DataArrays for different variables.
We’ll start off with creating a layer
DataArray, which is a 1D DataArray, containing all layer information required.
= xr.DataArray(data=[1], coords={"layer": [1]}, dims=("layer",))
layer
layer
<xarray.DataArray (layer: 1)> Size: 8B array([1]) Coordinates: * layer (layer) int64 8B 1
Next, we’ll create a 2D planar template grid. The xarray function .ones_like
returns a new object of ones with the same shape and type as a given dataarray or dataset.
= imod.typing.grid.ones_like(elevation, dtype=int)
template_2d
template_2d
<xarray.DataArray 'elevation' (mesh2d_nFaces: 5248)> Size: 42kB array([1, 1, 1, ..., 1, 1, 1]) Coordinates: mesh2d_face_x (mesh2d_nFaces) float64 42kB ... mesh2d_face_y (mesh2d_nFaces) float64 42kB ... * mesh2d_nFaces (mesh2d_nFaces) int64 42kB 0 1 2 3 4 ... 5244 5245 5246 5247 Attributes: unit: m NAP
Finally, we’ll create our unstructured template 3D grid.
We’ll multiply our template 2d grid with the layer
template to construct the full 3D grid. Pay attention to the dimension order, iMOD Python is very specific about this. We’ll put the dimensions in the right order immediately, to save ourselves headaches later.
# iMOD Python requires the layer dimension first, faces dimension second
= (template_2d * layer).transpose("layer", "mesh2d_nFaces")
template # plot
=1).ugrid.plot(vmin=0, vmax=2, cmap="viridis") template.sel(layer
Drainage
We can now use the template grid to transform the 2D elevation grid to a 3D elevation grid, required for the drainage package.
# Drain info
= template * elevation # requires layer dimension
drainage_lvl # plot
=1).ugrid.plot(vmin=-5, vmax=40, cmap="viridis") drainage_lvl.sel(layer
Compute conductance. The function .ugrid.grid.area
yields a numpy array with for each unstructured element its area (m2). Multiplied with the unstructured grid with value 1 yiels grid with area.
= template * template.ugrid.grid.area # [m2]
cell_area = 1 # [d]
resistance = cell_area / resistance # [m2/d]
conductance # Plot
=1).ugrid.plot() conductance.sel(layer
Discretization
Prepare our model discretization. idomain
specifies whether a call is active [1], inactive [0], or vertical pass through [-1]. In this case we’ll go for a grid that is active everywhere.
= template.copy() idomain
We’ll also define the bottom of the model. To keep things simple, we’ll subtract 300 from the drainage level, which will also be our top. In this way, we’ll construct a model with thickness 300 everywhere.
= drainage_lvl - 300.0
bottom # Plot
=1).ugrid.plot() bottom.sel(layer
Recharge
Let’s prepare our recharge. Note that Modflow 6 requires consistent length units, so we have to provide recharge in m/d instead of mm/d!
# Define recharge
= template * 0.002
recharge_rate # Plot
=1).ugrid.plot() recharge_rate.sel(layer
Construct simulation
Initialization of the main MODFLOW 6 object. We give it a name.
= imod.mf6.Modflow6Simulation("nl_elevation") simulation
Within a Modflow6Simulation
object one or more models can be defined. In this example there is just one ground water flow model.
In the next cells we will assemble our groundwater model in the object “gwf_model”. Further down, we add it to the object “simulation”.
In the next cell we define that option newton is chosen.
= imod.mf6.GroundwaterFlowModel(newton=True) gwf_model
The DIS package with vertices discretization needs input on tops, bottom and idomain.
"dis"] = imod.mf6.VerticesDiscretization(
gwf_model[=elevation, idomain=idomain, bottom=bottom
top )
Assign storage package:
"sto"] = imod.mf6.SpecificStorage(
gwf_model[=1e-5, specific_yield=0.01, transient=False, convertible=1
specific_storage )
Assign initial conditions, we’ll set them equal to drainage level.
"ic"] = imod.mf6.InitialConditions(start=drainage_lvl) gwf_model[
The Node Property Flow is most noteworthy where we specify the hydraulic conductivities. For a simple homogeneous, isotropic medium of 10 m/d:
"npf"] = imod.mf6.NodePropertyFlow(icelltype=1, k=10.0, k33=10.0) gwf_model[
Assign recharge package:
"rch"] = imod.mf6.Recharge(rate=recharge_rate) gwf_model[
Assign drainage package:
"drn"] = imod.mf6.Drainage(elevation=drainage_lvl, conductance=conductance) gwf_model[
"oc"] = imod.mf6.OutputControl(save_head="last", save_budget="last") gwf_model[
Now the ground water model object ‘gwf_model’ is ready.
Let’s add this model to the MODFLOW 6 simulation object.
"gwf"] = gwf_model simulation[
The final things we have to define are specifying a solver settings and time discretization to the Modflow 6 simulation.
"solver"] = imod.mf6.SolutionPresetModerate(modelnames=["gwf"]) simulation[
=["2000-01-01", "2001-01-01"]) simulation.create_time_discretization(additional_times
Write model
= imod.util.temporary_directory()
tmpdir
= tmpdir / "reference"
modeldir simulation.write(modeldir)
Now run our simulation:
# mf6_path = "path/to/mf6.exe"
# If you installed Modflow6 in your PATH environment variable, you can use the
# following argument:
= "mf6"
mf6_path
simulation.run(mf6_path)
Let’s look at the results:
= simulation.open_head().load().isel(time=-1, layer=0)
head # calculate ground water below surface level.
= elevation - head
groundwater_depth # Plot
groundwater_depth.ugrid.plot()
iMOD Python also supports partitioning a simulation, for parallel computation. It offers convenience functions to split existing simulations. For this we require a label array first.
from imod.mf6.multimodel.partition_generator import get_label_array
= get_label_array(simulation, npartitions=4)
label_array
label_array.ugrid.plot()
Let’s split our simulation with the label array.
= simulation.split(label_array) simulation_split
Now write the split simulation and run it
= tmpdir / "partitioned"
modeldir
simulation_split.write(modeldir) simulation_split.run(mf6_path)
And plot the results. Note that iMOD Python takes care of merging the 4 partitioned models into one single head grid.
= simulation_split.open_head()["head"].isel(time=-1, layer=0)
head_merged # plot
=-20, vmax=90) head_merged.ugrid.plot(vmin