Topsystem: from 2D map to 3D model#

iMOD Python has multiple features to help you define the topsystem of the groundwater system. With “topsystem” we mean all forcings which act on the top of the groundwater system. These are usually either meteorological (precipitation & evapotranspiration) or hydrological (rivers, ditches, lakes, sea) in nature. In MODFLOW6 these are usually simulated with the DRN, RIV, GHB, and RCH package. This data is usually provided as planar grids (x, y) without any vertical dimension. This user guide will show you how to allocate these forcings across model layers to grid cells and how to distribute conductances for Robin-like boundary conditions (RIV, DRN, GHB) over model layers. We will demonstrate this with the RIV package, as this package supports all options for allocating cells and distributing conductances. A subset of options are available for the DRN, GHB, and RCH packages.

Overview allocation options#

There are multiple options available to allocate rivers. For a full description of all options, see the documentation of imod.prepare.ALLOCATION_OPTION(). We can print all possible options as follows:

for option in ALLOCATION_OPTION:
    print(option.name)

stage_to_riv_bot
first_active_to_elevation
stage_to_riv_bot_drn_above
at_elevation
at_first_active

Let’s make plots for each option to visualize the effect of each choice.

# Create grid for plots
fig, axes = plt.subplots(3, 2, sharex=True, sharey=True, figsize=(11, 11))
axes = np.ravel(axes)

# The top left plot shows the elevation of the river bottom (upward triangle)
# and river stage (downward triangle).
ax = axes[0]
imod.visualize.cross_section(
    xsection_layer_nr,
    "tab20",
    np.arange(21),
    kwargs_colorbar={"plot_colorbar": False},
    fig=fig,
    ax=ax,
)
ax.scatter(x_line, stage_line.values, marker=7, c="k")
ax.scatter(x_line, stage_bottom.values, marker=6, c="k")
ax.set_title("stage and bottom elevation")

# Loop over allocation options, and plot the allocated cells as a polygon,
# using the "aquitard" feature of the cross_section plot.
for i, option in enumerate(ALLOCATION_OPTION, start=1):
    riv_allocated, _ = allocate_riv_cells(
        allocation_option=option,
        active=layer_model["idomain"] == 1,
        top=layer_model["top"],
        bottom=layer_model["bottom"],
        stage=planar_river["stage"],
        bottom_elevation=planar_river["bottom"],
    )
    riv_allocated.coords["top"] = layer_model["top"]
    riv_allocated.coords["bottom"] = layer_model["bottom"]

    xsection_allocated = imod.select.cross_section_linestring(riv_allocated, geometry)
    ax = axes[i]
    if (i % 2) == 0:
        kwargs_colorbar = {"plot_colorbar": False}
    else:
        kwargs_colorbar = {"plot_colorbar": True}
    kwargs_aquitards = {"edgecolor": "k", "facecolor": "grey"}
    imod.visualize.cross_section(
        xsection_layer_nr,
        "tab20",
        np.arange(21),
        aquitards=xsection_allocated,
        kwargs_aquitards=kwargs_aquitards,
        kwargs_colorbar=kwargs_colorbar,
        fig=fig,
        ax=ax,
    )
    ax.set_title(f"option: {option.name}")

# Enforce tight layout to remove whitespace inbetween plots.
plt.tight_layout()

stage and bottom elevation, option: stage_to_riv_bot, option: first_active_to_elevation, option: stage_to_riv_bot_drn_above, option: at_elevation, option: at_first_active

You can see the chosen option matters quite a lot. at_elevation allocates cells in the model layer containing the river bottom elevation. at_first_active allocates only at the first active model layer. first_active_to_elevation allocates cells from first active model layer to the model layer containing the river bottom elevation. stage_to_riv_bot allocates cells from the model layer containing river stage up until the model layer containing bottom elevation. Finally stage_to_riv_bot_drn_above allocates river cells from the model layer containing river stage to the model layer containing the river bottom elevation and allocates drain cells from the first active model layer to the model layer containing the river stage elevation. The allocated drain cells are not shown in the plot.