Regional model#

This example shows a simplified script for building a groundwater model in the northeast of the Netherlands. A primary feature of this area is an ice-pushed ridge called the Hondsrug. This examples demonstrates modifying external data for use in a MODFLOW6 model.

In overview, the model features:

  • Thirteen layers: seven aquifers and six aquitards;

  • A dense ditch network in the east;

  • Pipe drainage for agriculture;

  • Precipitation and evapotranspiration.


We’ll start with the usual imports, and an import from scipy.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.ndimage.morphology
import xarray as xr

import imod

Before starting to create the input data, we will create the groundwater model variable (gwf_model) using imod.mf6.GroundwaterFlowModel. The data from all the model packages will be added to this variable.

gwf_model = imod.mf6.GroundwaterFlowModel()

This package allows specifying a regular MODFLOW grid. This grid is assumed to be rectangular horizontally, but can be distorted vertically.

Load data#

We’ll load the data from the examples that come with this package.

layermodel = imod.data.hondsrug_layermodel()
idomain = layermodel["idomain"].astype(int)
top = layermodel["top"]
bot = layermodel["bottom"]

Adding information to the DIS package#

The following step is to add the previously created discretization data to the gwf_model variable. This is done using the function imod.mf6.StructuredDiscretization. The data to include is the top of the model domain, the bottom of the layers, and the idomain. All this information comes from the previously imported tifs (now converted to xarray.DataArray.

gwf_model["dis"] = imod.mf6.StructuredDiscretization(
    top=top, bottom=bot, idomain=idomain
)

Node property flow package - NPF#

This package contains the information related to the aquifer properties used to calculate hydraulic conductance. This package replaces the Layer Property Flow (LPF), Block-Centered Flow (BCF), and Upstream Weighting (UPW) packages from previous MODFLOW versions.

Hydraulic conductivity#

k = layermodel["k"]

icelltype#

The cell type to be used in the model (confined or convertible) can be defined under ICELLTYPE, which is an input to the NPF package. ICELLTYPE == 0: Confined cell - Constant transmissivity ICELLTYPE != 0: Convertible cell - Transmissivity varies depending on the calculated head in the cell (based on the saturated cell thickness)

In this example, all layers have a ICELLTYPE equal to 0, indicating confined cells. This is defined in the following section.

Adding information to the NPF package#

The information for the NPF package is added to the gwf_model variable using imod.mf6.NodePropertyFlow. The information included is the icelltype value (equal to zero), the array for the hydraulic conductivity (considered to be the same for the horizontal and vertical direction) and, optionally, the variable_vertical_conductance, dewatered, perched and save_flows options have been activated. For more details about the meaning of these variables and other variables available to be used within this package, please refer to the documentation.

gwf_model["npf"] = imod.mf6.NodePropertyFlow(
    icelltype=0,
    k=k,
    k33=k,
    variable_vertical_conductance=True,
    dewatered=True,
    perched=True,
    save_flows=True,
)

Initial conditions package - IC#

This package reads the starting heads for a simulation.

Starting heads interpolation#

The starting heads to be used in this model are based on the interpolation of x-y head measurements, which were interpolated on a larger area. This example was created in this example –insert-link-here–

The heads were interpolated on a larger area, therefore these have to be clipped first

initial = imod.data.hondsrug_initial()
interpolated_head_larger = initial["head"]

xmin = 237_500.0
xmax = 250_000.0
ymin = 559_000.0
ymax = 564_000.0

interpolated_head = interpolated_head_larger.sel(
    x=slice(xmin, xmax), y=slice(ymax, ymin)
)

# Plotting the clipped interpolation
fig, ax = plt.subplots()
interpolated_head.plot.imshow(ax=ax)
dx = 25.0, dy = -25.0
<matplotlib.image.AxesImage object at 0x7ff23ffb5a20>

The final step is to assign the 2D heads interpolation to all the model layers (as a reference value) by using the xarray tool xarray.full_like. The 3d idomain array is used as reference for the geometry and then its original values are replaced by NaNs. This array is combined with the interpolated_head array using the xarray DataArray.combine_first option. The final result is an starting_heads xarray where all layers have the 2d interpolated_head information.

# Assign interpolated head values to the all the model layers
like_3d = xr.full_like(idomain, np.nan, dtype=float)
starting_head = like_3d.combine_first(interpolated_head)
starting_head
<xarray.DataArray 'idomain' (layer: 13, y: 200, x: 500)>
array([[[5.6706604 , 5.65319374, 5.63573351, ..., 1.27008811,
         1.27823924, 1.28676046],
        [5.67764711, 5.66021213, 5.64275829, ..., 1.26336139,
         1.27183901, 1.28070347],
        [5.68462064, 5.66720984, 5.64976361, ..., 1.25620192,
         1.26504986, 1.27431479],
        ...,
        [7.28608515, 7.26550211, 7.24273429, ..., 4.67051889,
         4.66997651, 4.66944365],
        [7.2907364 , 7.26985601, 7.24672249, ..., 4.70303061,
         4.70234778, 4.70168152],
        [7.29728577, 7.27649525, 7.25355321, ..., 4.73551559,
         4.73469247, 4.73389453]],

       [[5.6706604 , 5.65319374, 5.63573351, ..., 1.27008811,
         1.27823924, 1.28676046],
        [5.67764711, 5.66021213, 5.64275829, ..., 1.26336139,
         1.27183901, 1.28070347],
        [5.68462064, 5.66720984, 5.64976361, ..., 1.25620192,
         1.26504986, 1.27431479],
...
        [7.28608515, 7.26550211, 7.24273429, ..., 4.67051889,
         4.66997651, 4.66944365],
        [7.2907364 , 7.26985601, 7.24672249, ..., 4.70303061,
         4.70234778, 4.70168152],
        [7.29728577, 7.27649525, 7.25355321, ..., 4.73551559,
         4.73469247, 4.73389453]],

       [[5.6706604 , 5.65319374, 5.63573351, ..., 1.27008811,
         1.27823924, 1.28676046],
        [5.67764711, 5.66021213, 5.64275829, ..., 1.26336139,
         1.27183901, 1.28070347],
        [5.68462064, 5.66720984, 5.64976361, ..., 1.25620192,
         1.26504986, 1.27431479],
        ...,
        [7.28608515, 7.26550211, 7.24273429, ..., 4.67051889,
         4.66997651, 4.66944365],
        [7.2907364 , 7.26985601, 7.24672249, ..., 4.70303061,
         4.70234778, 4.70168152],
        [7.29728577, 7.27649525, 7.25355321, ..., 4.73551559,
         4.73469247, 4.73389453]]])
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int32 1 2 3 4 5 6 7 8 9 10 11 12 13


Adding information to the IC package#

The function for indicating the initial conditions is imod.mf6.InitialConditions. It is necessary to indicate the value(s) to be considered as the initial (starting) head of the simulation. In this case, this value is equal to the previously created starting_head array.

gwf_model["ic"] = imod.mf6.InitialConditions(starting_head)

Constant head package - CHD#

This package allows to indicate if the head varies with time, if it is constant or if it is inactive.

Constant head edge#

The previously interpolated starting_head array will be used to define the constant head value which will be used along the model boundaries. A function is defined to indicate the location of the outer edge (returning a boolean array).

def outer_edge(da):
    data = da.copy()
    from_edge = scipy.ndimage.morphology.binary_erosion(data)
    is_edge = (data == 1) & (from_edge == 0)
    return is_edge.astype(bool)

For the next calculations, it is necessary to create a template array which can be used for assigning the corresponding geometry to other arrays. In this case, a 2d template is created based on the idomain layer information and filled with ones.

like_2d = xr.full_like(idomain.isel(layer=0), 1)
like_2d
<xarray.DataArray 'idomain' (y: 200, x: 500)>
array([[1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       ...,
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1]])
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
    layer    int32 1


Using the previously created function and the 2d template, the outer edge is defined for this example.

edge = outer_edge(xr.full_like(like_2d.drop("layer"), 1))
/builds/deltares/imod/imod-python/examples/mf6/hondsrug.py:200: DeprecationWarning: Please use `binary_erosion` from the `scipy.ndimage` namespace, the `scipy.ndimage.morphology` namespace is deprecated.
  from_edge = scipy.ndimage.morphology.binary_erosion(data)

Adding information to the CHD package#

To add the information to the CHD package within the gwf_model variable, the imod.mf6.ConstantHead function is used. The required information is the head array for this boundary condition. In this example, the starting_head array is selected where the idomain is > 0 (active) and it is located in the edge array.

It is also possible (and optional) to indicate if the CHD information will be written to the listing file after it is read (print_input), if the constant head flow rates will be printed to the listing file for every stress period time step in which “BUDGET PRINT” is specified in Output Control (print_flows) and if the constant head flow terms will be written to the file specified with “BUDGET FILEOUT” in Output Control (save_flows). By default, these three options are set to False.

gwf_model["chd"] = imod.mf6.ConstantHead(
    starting_head.where((idomain > 0) & edge),
    print_input=False,
    print_flows=True,
    save_flows=True,
)

Recharge#

This package is used to represent areally distributed recharge to the groundwater system. To calculate the recharge, the precipitation and evapotranspiration information from the KNMI website has been downloaded for the study area. This information is saved in netCDF files, which have been imported using the xarray function xr.open_dataset, slicing the area to the model’s miminum and maximum dimensions.

Note that the meteorological data has mm/d as unit and this has to be converted to m/d for Modflow 6.

xmin = 230_000.0
xmax = 257_000.0
ymin = 550_000.0
ymax = 567_000.0

meteorology = imod.data.hondsrug_meteorology()
pp = meteorology["precipitation"]
evt = meteorology["evapotranspiration"]

pp = pp.sel(x=slice(xmin, xmax), y=slice(ymax, ymin)) / 1000.0  # from mm/d to m/d
evt = evt.sel(x=slice(xmin, xmax), y=slice(ymax, ymin)) / 1000.0  # from mm/d to m/d

Recharge - Steady state#

For the steady state conditions of the model, the data from the period 2000 to 2009 was considered as reference. The initial information was sliced to this time period and averaged to obtain the a mean value grid. This process was done for both precipitation and evapotranspiration datasets.

Precipitation

pp_ss = pp.sel(time=slice("2000-01-01", "2009-12-31"))
pp_ss_mean = pp_ss.mean(dim="time")

fig, ax = plt.subplots()
pp_ss_mean.plot(ax=ax)
dx = 1e+03, dy = -1e+03
<matplotlib.collections.QuadMesh object at 0x7ff240a05360>

Evapotranspiration

evt_ss = evt.sel(time=slice("2000-01-01", "2009-12-31"))
evt_ss_mean = evt_ss.mean(dim="time")

fig, ax = plt.subplots()
evt_ss_mean.plot(ax=ax)
dx = 1e+03, dy = -1e+03
<matplotlib.collections.QuadMesh object at 0x7ff240bf6ec0>

For the recharge calculation, a first estimate is the difference between the precipitation and evapotranspiration values.

rch_ss = pp_ss_mean - evt_ss_mean

fig, ax = plt.subplots()
rch_ss.plot.imshow(ax=ax)
dx = 1e+03, dy = -1e+03
<matplotlib.image.AxesImage object at 0x7ff240ba9600>

Recharge - Transient#

The transient model will encompass the period from 2010 to 2015. The initial pp and evt datasets have been sliced to this time frame.

pp_trans = pp.sel(time=slice("2010-01-01", "2015-12-31"))
evt_trans = evt.sel(time=slice("2010-01-01", "2015-12-31"))

As previously done, it is assumed that the recharge is equal to the difference between precipitation and evapotranspiration as a first estimate. Furthermore, the negative values found after doing this calculation have been replaced by zeros, as the recharge should not have a negative value.

rch_trans = pp_trans - evt_trans
rch_trans = rch_trans.where(rch_trans > 0, 0)  # check negative values

The original information is on a daily step, so it is going to be resampled to a yearly step by using the xarray function Dataset.resample.

rch_trans_yr = rch_trans.resample(time="A", label="left").mean()
rch_trans_yr
<xarray.DataArray (time: 6, y: 17, x: 27)>
array([[[0.00191725, 0.00191589, 0.00191248, ..., 0.00158879,
         0.00160275, 0.00161772],
        [0.00193096, 0.00193016, 0.0019277 , ..., 0.00159777,
         0.00160856, 0.00162071],
        [0.00194592, 0.00194606, 0.00194384, ..., 0.00160863,
         0.00161658, 0.0016263 ],
        ...,
        [0.00205102, 0.0020491 , 0.00204451, ..., 0.00175834,
         0.00175477, 0.00175247],
        [0.0020481 , 0.00204493, 0.00203995, ..., 0.00177083,
         0.0017666 , 0.00176362],
        [0.00204444, 0.002041  , 0.00203658, ..., 0.00178251,
         0.00177807, 0.00177411]],

       [[0.00175232, 0.0017549 , 0.00175632, ..., 0.00160428,
         0.0016166 , 0.00163142],
        [0.00175832, 0.00176321, 0.00176716, ..., 0.00161909,
         0.00162744, 0.00163787],
        [0.00176593, 0.00177297, 0.00177865, ..., 0.00163613,
         0.00164029, 0.00164618],
...
        [0.00190413, 0.0019    , 0.00189647, ..., 0.00168478,
         0.00167463, 0.00166721],
        [0.00191761, 0.00191336, 0.00190866, ..., 0.00169624,
         0.00168655, 0.0016792 ],
        [0.0019308 , 0.00192603, 0.0019212 , ..., 0.00170712,
         0.00169771, 0.00169032]],

       [[0.00226806, 0.00225167, 0.00223808, ..., 0.00221006,
         0.00219559, 0.00218172],
        [0.0022924 , 0.00227892, 0.00226807, ..., 0.00220936,
         0.00219364, 0.00217812],
        [0.00231699, 0.00230603, 0.00229702, ..., 0.00221065,
         0.0021932 , 0.00217596],
        ...,
        [0.00249076, 0.00248442, 0.00247622, ..., 0.00224485,
         0.00222806, 0.00221361],
        [0.00249192, 0.00248449, 0.00247561, ..., 0.00224585,
         0.00223096, 0.00221743],
        [0.00249155, 0.00248376, 0.00247603, ..., 0.00224736,
         0.00223372, 0.00222133]]], dtype=float32)
Coordinates:
  * x        (x) float64 2.305e+05 2.315e+05 2.325e+05 ... 2.555e+05 2.565e+05
  * y        (y) float64 5.665e+05 5.655e+05 5.645e+05 ... 5.515e+05 5.505e+05
    dx       float64 1e+03
    dy       float64 -1e+03
  * time     (time) datetime64[ns] 2009-12-31 2010-12-31 ... 2014-12-31


To create the final recharge for the transient simulation, the steady state information needs to be concatenated to the transient recharge data. The steady state simulation will be run for one second. This is achieved by using the numpy Timedelta function, first creating a time delta of 1 second, which is assigned to the steady state recharge information. This dataset is then concatenated using the xarray function xarray.concat to the transient information and indicating that the dimension to join is “time”.

starttime = "2009-12-31"

# Add first steady-state
timedelta = np.timedelta64(1, "s")  # 1 second duration for initial steady-state
starttime_steady = np.datetime64(starttime) - timedelta
rch_ss = rch_ss.assign_coords(time=starttime_steady)

rch_ss_trans = xr.concat([rch_ss, rch_trans_yr], dim="time")
rch_ss_trans
<xarray.DataArray (y: 17, x: 27, time: 7)>
array([[[0.00091255, 0.00191725, 0.00175232, ..., 0.00167949,
         0.00161827, 0.00226806],
        [0.0009093 , 0.00191589, 0.0017549 , ..., 0.001689  ,
         0.00161518, 0.00225167],
        [0.00090591, 0.00191248, 0.00175632, ..., 0.00169961,
         0.00161301, 0.00223808],
        ...,
        [0.00078183, 0.00158879, 0.00160428, ..., 0.00174276,
         0.00162122, 0.00221006],
        [0.00078091, 0.00160275, 0.0016166 , ..., 0.00174685,
         0.0016207 , 0.00219559],
        [0.00078071, 0.00161772, 0.00163142, ..., 0.00175276,
         0.00162407, 0.00218172]],

       [[0.000913  , 0.00193096, 0.00175832, ..., 0.00168658,
         0.00164151, 0.0022924 ],
        [0.0009098 , 0.00193016, 0.00176321, ..., 0.00169678,
         0.00163919, 0.00227892],
        [0.0009065 , 0.0019277 , 0.00176716, ..., 0.00170842,
         0.00163841, 0.00226807],
...
        [0.00085378, 0.00177083, 0.00186831, ..., 0.00174713,
         0.00169624, 0.00224585],
        [0.00084055, 0.0017666 , 0.00185661, ..., 0.00173329,
         0.00168655, 0.00223096],
        [0.00082746, 0.00176362, 0.00184691, ..., 0.00172136,
         0.0016792 , 0.00221743]],

       [[0.00098196, 0.00204444, 0.00188612, ..., 0.00182747,
         0.0019308 , 0.00249155],
        [0.00098031, 0.002041  , 0.00189853, ..., 0.00183236,
         0.00192603, 0.00248376],
        [0.00097932, 0.00203658, 0.00191231, ..., 0.00183771,
         0.0019212 , 0.00247603],
        ...,
        [0.00086102, 0.00178251, 0.00188156, ..., 0.00175411,
         0.00170712, 0.00224736],
        [0.00084771, 0.00177807, 0.00187044, ..., 0.0017413 ,
         0.00169771, 0.00223372],
        [0.0008352 , 0.00177411, 0.00186046, ..., 0.00173005,
         0.00169032, 0.00222133]]], dtype=float32)
Coordinates:
  * x        (x) float64 2.305e+05 2.315e+05 2.325e+05 ... 2.555e+05 2.565e+05
  * y        (y) float64 5.665e+05 5.655e+05 5.645e+05 ... 5.515e+05 5.505e+05
    dx       float64 1e+03
    dy       float64 -1e+03
  * time     (time) datetime64[ns] 2009-12-30T23:59:59 2009-12-31 ... 2014-12-31


The data obtained from KNMI has different grid dimensions than the one considered in this example. To fix this, imod-python includes the option imod.prepare.Regridder, which modifies the original grid dimensions to a different one. It is also possible to define the regridding method such as nearest, multilinear, mean, among others. In this case, mean was selected and the 2d template (like_2d) was used as reference, as this is the geometry to be considered in the model.

rch_ss_trans = imod.prepare.Regridder(method="mean").regrid(rch_ss_trans, like_2d)
rch_ss_trans
<xarray.DataArray (y: 200, x: 500, time: 7)>
array([[[0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
         0.00169166, 0.00232221],
        [0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
         0.00169166, 0.00232221],
        [0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
         0.00169166, 0.00232221],
        ...,
        [0.00083373, 0.00163785, 0.00171206, ..., 0.00176909,
         0.00164895, 0.00231468],
        [0.00083373, 0.00163785, 0.00171206, ..., 0.00176909,
         0.00164895, 0.00231468],
        [0.00083373, 0.00163785, 0.00171206, ..., 0.00176909,
         0.00164895, 0.00231468]],

       [[0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
         0.00169166, 0.00232221],
        [0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
         0.00169166, 0.00232221],
        [0.00089943, 0.00190719, 0.00182535, ..., 0.00179266,
         0.00169166, 0.00232221],
...
        [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
         0.00171221, 0.0023817 ],
        [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
         0.00171221, 0.0023817 ],
        [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
         0.00171221, 0.0023817 ]],

       [[0.00091826, 0.00198224, 0.00191962, ..., 0.00183427,
         0.00179237, 0.00242177],
        [0.00091826, 0.00198224, 0.00191962, ..., 0.00183427,
         0.00179237, 0.00242177],
        [0.00091826, 0.00198224, 0.00191962, ..., 0.00183427,
         0.00179237, 0.00242177],
        ...,
        [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
         0.00171221, 0.0023817 ],
        [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
         0.00171221, 0.0023817 ],
        [0.00088318, 0.00169605, 0.00183059, ..., 0.00180298,
         0.00171221, 0.0023817 ]]])
Coordinates:
  * time     (time) datetime64[ns] 2009-12-30T23:59:59 2009-12-31 ... 2014-12-31
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
    dx       float64 25.0
    dy       float64 -25.0
    layer    int32 1


The previously created recharge array is a 2D array that needs to be assigned to a 3D array. This is done using the xarray DataArray.where option, where the recharge values are applied to the cells where the idomain value is larger than zero (that is, the active cells) and for the uppermost active cell (indicated by the minimum layer number).

rch_total = rch_ss_trans.where(
    idomain["layer"] == idomain["layer"].where(idomain > 0).min("layer")
)
rch_total
<xarray.DataArray (y: 200, x: 500, time: 7, layer: 13)>
array([[[[       nan,        nan, 0.00089943, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00190719, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00182535, ...,        nan,
                 nan,        nan],
         ...,
         [       nan,        nan, 0.00179266, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00169166, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00232221, ...,        nan,
                 nan,        nan]],

        [[       nan,        nan, 0.00089943, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00190719, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00182535, ...,        nan,
                 nan,        nan],
...
         [       nan,        nan, 0.00180298, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00171221, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.0023817 , ...,        nan,
                 nan,        nan]],

        [[       nan,        nan, 0.00088318, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00169605, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00183059, ...,        nan,
                 nan,        nan],
         ...,
         [       nan,        nan, 0.00180298, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.00171221, ...,        nan,
                 nan,        nan],
         [       nan,        nan, 0.0023817 , ...,        nan,
                 nan,        nan]]]])
Coordinates:
  * time     (time) datetime64[ns] 2009-12-30T23:59:59 2009-12-31 ... 2014-12-31
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int32 1 2 3 4 5 6 7 8 9 10 11 12 13


Finally, transposing the array dimensions using DataArray.transpose so they are in the correct order.

rch_total = rch_total.transpose("time", "layer", "y", "x")
rch_total

fig, ax = plt.subplots()
rch_total.isel(layer=2, time=6).plot.imshow(ax=ax)
time = 2014-12-31, dx = 25.0, dy = -25.0, layer...
<matplotlib.image.AxesImage object at 0x7ff23feaef50>

Adding information to the RCH package#

The information for the RCH package is added with the function imod.mf6.Recharge. It is required to insert the recharge flux rate, and it is optional to include the print_input, print_flows and save_flows information.

gwf_model["rch"] = imod.mf6.Recharge(rch_total)

Drainage package - DRN#

The drain package is used to simulate features that remove water from the aquifer, such as agricultural drains or springs. This occurs at a rate proportional to the head difference between the head in the aquifer and the drain elevation (the head in the aquifer has to be above that elevation). The conductance is the proportionality constant.

Import drainage information#

drainage = imod.data.hondsrug_drainage()

pipe_cond = drainage["conductance"]
pipe_elev = drainage["elevation"]

pipe_cond
<xarray.DataArray 'conductance' (layer: 4, y: 200, x: 500)>
[400000 values with dtype=float32]
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int32 1 3 5 6


Adding information to the DRN package#

To add the information to the DRN package within the gwf_model variable, the imod.mf6.Drainage function is used. It is required to add the previously created arrays for the drain elevation and the drain conductance. It is optional to insert the information for print_input, print_flows and save_flows which are set to False by default.

gwf_model["drn-pipe"] = imod.mf6.Drainage(conductance=pipe_cond, elevation=pipe_elev)

River package - RIV#

This package simulates the effects of flow between surface-water features and groundwater systems.

Import river information#

river = imod.data.hondsrug_river()
riv_cond = river["conductance"]
riv_stage = river["stage"]
riv_bot = river["bottom"]

Adding information to the RIV package#

The data is assigned to the gwf_model variable by using imod.mf6.River, based on the previously imported conductance, stage and bottom arrays.

gwf_model["riv"] = imod.mf6.River(
    conductance=riv_cond, stage=riv_stage, bottom_elevation=riv_bot
)

Storage package - STO#

When the STO Package is included in a model, storage changes will be calculated, and thus, the model will be transient.

ss = 0.0003
layer = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])
sy = xr.DataArray(
    [0.16, 0.16, 0.16, 0.16, 0.15, 0.15, 0.15, 0.15, 0.14, 0.14, 0.14, 0.14, 0.14],
    {"layer": layer},
    ("layer",),
)
conv = xr.full_like(idomain, 0, dtype=int)

times_sto = [
    np.datetime64("2009-12-30T23:59:59.000000000"),
    np.datetime64("2009-12-31T00:00:00.000000000"),
    np.datetime64("2010-12-31T00:00:00.000000000"),
    np.datetime64("2011-12-31T00:00:00.000000000"),
    np.datetime64("2012-12-31T00:00:00.000000000"),
    np.datetime64("2013-12-31T00:00:00.000000000"),
    np.datetime64("2014-12-31T00:00:00.000000000"),
]

transient = xr.DataArray(
    [False, True, True, True, True, True, True], {"time": times_sto}, ("time",)
)

Adding information to the STO package#

The data is assigned to the gwf_model variable by using imod.mf6.Storage. It is necessary to indicate the values of specific storage, specific yield and if the layers are convertible.

gwf_model["sto"] = imod.mf6.SpecificStorage(
    specific_storage=ss, specific_yield=sy, transient=transient, convertible=conv
)

Output Control package - OC#

This package determines how and when heads are printed to the listing file and/or written to a separate binary output file

Adding information to the OC package#

The function imod.mf6.OutputControl is used to store the information for this package. It is possible to indicate if the heads and budget information is saved at the end of each stress period (last), for all timesteps a stress period (all), or at the start of a stress period (first)

gwf_model["oc"] = imod.mf6.OutputControl(save_head="last", save_budget="last")

Model simulation#

In MODFLOW 6, the concept of “model” is that part of the program that solves a hydrologic process. MODFLOW 6 documentation supports one type of model: the GWF Model. It is possible within the MODFLOW 6 framewotk to solve multiple, tightly coupled, numerical models in a single system of equation, which may be multiple models of the same type or of different types.

The previously created gwf_model variable now contains the information from all the variables.

gwf_model
{'dis': <imod.mf6.dis.StructuredDiscretization object at 0x7ff249018ee0>, 'npf': <imod.mf6.npf.NodePropertyFlow object at 0x7ff24901b1f0>, 'ic': <imod.mf6.ic.InitialConditions object at 0x7ff2492c2020>, 'chd': <imod.mf6.chd.ConstantHead object at 0x7ff23ffa7b20>, 'rch': <imod.mf6.rch.Recharge object at 0x7ff23fec1fc0>, 'drn-pipe': <imod.mf6.drn.Drainage object at 0x7ff23fd62050>, 'riv': <imod.mf6.riv.River object at 0x7ff23fd611b0>, 'sto': <imod.mf6.sto.SpecificStorage object at 0x7ff23fd77d30>, 'oc': <imod.mf6.oc.OutputControl object at 0x7ff23fd77940>}

Attach the model information to a simulation#

The function imod.mf6.Modflow6Simulation allows to assign models to a simulation (in this case, the gwf_model).

simulation = imod.mf6.Modflow6Simulation("mf6-mipwa2-example")
simulation["GWF_1"] = gwf_model

Solver settings#

The solver settings are indicated using imod.mf6.Solution. If the values are not indicated manually, the defaults values will be considered.

simulation["solver"] = imod.mf6.Solution(
    modelnames=["GWF_1"],
    print_option="summary",
    csv_output=False,
    no_ptc=True,
    outer_dvclose=1.0e-4,
    outer_maximum=500,
    under_relaxation=None,
    inner_dvclose=1.0e-4,
    inner_rclose=0.001,
    inner_maximum=100,
    linear_acceleration="cg",
    scaling_method=None,
    reordering_method=None,
    relaxation_factor=0.97,
)

Assign time discretization#

The time discretization of this model is 6 years.

simulation.create_time_discretization(
    additional_times=["2009-12-30T23:59:59.000000000", "2015-12-31T00:00:00.000000000"]
)

Run the model#

Note

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.

modeldir = imod.util.temporary_directory()
simulation.write(modeldir)
simulation.run()

Results visualization#

The next section indicated how to visualize the model results.

Import heads results#

The heads results are imported using imod.mf6.open_hds and indicating the location of the heads, as well as the DIS file.

hds = imod.mf6.open_hds(
    modeldir / "GWF_1" / "GWF_1.hds", modeldir / "GWF_1" / "dis.dis.grb"
)
hds
<xarray.DataArray 'head' (time: 7, layer: 13, y: 200, x: 500)>
dask.array<stack, shape=(7, 13, 200, 500), dtype=float64, chunksize=(1, 13, 200, 500), chunktype=numpy.ndarray>
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int64 1 2 3 4 5 6 7 8 9 10 11 12 13
  * time     (time) float64 1.157e-05 365.0 730.0 ... 1.826e+03 2.191e+03


We can plot the data of an individual layer as follows

fig, ax = plt.subplots()
hds.sel(layer=3).isel(time=3).plot(ax=ax)
dx = 25.0, dy = -25.0, layer = 3, time = 1.096e+03
<matplotlib.collections.QuadMesh object at 0x7ff23fdd0c10>

As you can see layer 3 has some missing cells in the west Whereas layer 4 only contains active cells in the eastern peatland area

fig, ax = plt.subplots()
hds.sel(layer=4).isel(time=3).plot(ax=ax)
dx = 25.0, dy = -25.0, layer = 4, time = 1.096e+03
<matplotlib.collections.QuadMesh object at 0x7ff23fc32890>

Layer 5 contains more data towards the west, but has no active cells in the centre.

fig, ax = plt.subplots()
hds.sel(layer=5).isel(time=3).plot(ax=ax)
dx = 25.0, dy = -25.0, layer = 5, time = 1.096e+03
<matplotlib.collections.QuadMesh object at 0x7ff23fb79450>

As you can see the data is individual layers have lots of inactive in different places.

It is difficult for this model to get a good idea what is happening across the area based on 1 layer alone. Luckily xarray allows us to compute the mean across a selection of layers and plot this.

By first selecting 3 layers with sel`, and then computing the mean across the layer dimension with mean(dim="layer").

fig, ax = plt.subplots()
hds.sel(layer=slice(3, 5)).mean(dim="layer").isel(time=3).plot(ax=ax)

# Assign dates to head
# --------------------
#
# MODFLOW6 has no concept of a calendar, so the output is not labelled only
# in terms of "time since start" in floating point numbers. For this model
# the time unit is days and we can assign a date coordinate as follows:

starttime = pd.to_datetime("2000-01-01")
timedelta = pd.to_timedelta(hds["time"], "D")
hds = hds.assign_coords(time=starttime + timedelta)
dx = 25.0, dy = -25.0, time = 1.096e+03

Extract head at points#

A typical operation is to extract simulated heads at point locations to compare them with measurements. In this example, we select the heads at two points:

x = [240_000.0, 244_000.0]
y = [560_000.0, 562_000.0]
selection = imod.select.points_values(hds, x=x, y=y)

The result can be converted into a pandas dataframe for timeseries analysis, or written to a variety of tabular file formats.

dataframe = selection.to_dataframe().reset_index()
dataframe = dataframe.rename(columns={"index": "id"})
dataframe
time layer id x y dx dy head
0 2000-01-01 00:00:00.999993600 1 0 240012.5 560012.5 25.0 -25.0 NaN
1 2000-01-01 00:00:00.999993600 1 1 244012.5 562012.5 25.0 -25.0 10.534013
2 2000-01-01 00:00:00.999993600 2 0 240012.5 560012.5 25.0 -25.0 NaN
3 2000-01-01 00:00:00.999993600 2 1 244012.5 562012.5 25.0 -25.0 NaN
4 2000-01-01 00:00:00.999993600 3 0 240012.5 560012.5 25.0 -25.0 7.682111
... ... ... ... ... ... ... ... ...
177 2005-12-31 00:00:00.999993600 11 1 244012.5 562012.5 25.0 -25.0 4.787419
178 2005-12-31 00:00:00.999993600 12 0 240012.5 560012.5 25.0 -25.0 NaN
179 2005-12-31 00:00:00.999993600 12 1 244012.5 562012.5 25.0 -25.0 NaN
180 2005-12-31 00:00:00.999993600 13 0 240012.5 560012.5 25.0 -25.0 7.498579
181 2005-12-31 00:00:00.999993600 13 1 244012.5 562012.5 25.0 -25.0 4.787407

182 rows × 8 columns



Total running time of the script: ( 1 minutes 56.040 seconds)

Gallery generated by Sphinx-Gallery