Code
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.
This data is distributed with the xugrid package. You need an internet connection to download the data.
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
[5248 values with dtype=float32]
[5248 values with dtype=float64]
[5248 values with dtype=float64]
array([ 0, 1, 2, ..., 5245, 5246, 5247])
PandasIndex(RangeIndex(start=0, stop=5248, step=1, name='mesh2d_nFaces'))
Let’s plot the data
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.
<xarray.DataArray (layer: 1)> Size: 8B array([1]) Coordinates: * layer (layer) int64 8B 1
array([1])
array([1])
PandasIndex(Index([1], dtype='int64', name='layer'))
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.
<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
array([1, 1, 1, ..., 1, 1, 1])
[5248 values with dtype=float64]
[5248 values with dtype=float64]
array([ 0, 1, 2, ..., 5245, 5246, 5247])
PandasIndex(RangeIndex(start=0, stop=5248, step=1, name='mesh2d_nFaces'))
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.
We can now use the template grid to transform the 2D elevation grid to a 3D elevation grid, required for the drainage package.
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.
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.
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.
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!
Initialization of the main MODFLOW 6 object. We give it a name.
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.
The DIS package with vertices discretization needs input on tops, bottom and idomain.
Assign storage package:
Assign initial conditions, we’ll set them equal to drainage level.
The Node Property Flow is most noteworthy where we specify the hydraulic conductivities. For a simple homogeneous, isotropic medium of 10 m/d:
Assign recharge package:
Assign drainage package:
Now the ground water model object ‘gwf_model’ is ready.
Let’s add this model to the MODFLOW 6 simulation object.
The final things we have to define are specifying a solver settings and time discretization to the Modflow 6 simulation.
Now run our simulation:
Let’s look at the results:
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.
Let’s split our simulation with the label array.
Now write the split simulation and run it
And plot the results. Note that iMOD Python takes care of merging the 4 partitioned models into one single head grid.