A RetroSearch Logo

Home - News ( United States | United Kingdom | Italy | Germany ) - Football scores

Search Query:

Showing content from https://deltares.github.io/imod-python/user-guide/09-topsystem.html below:

from 2D map to 3D model — iMOD Python 1.0.0rc4 documentation

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 MODFLOW 6 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.

Example data#

Let’s load the data first. We have a layer model containing a basic hydrogeological schemitization of our model, so the tops and bottoms of model layers, the hydraulic conductivity (k), and which cells are active (idomain=1) or vertical passthrough (idomain=-1).

import imod

layer_model = imod.data.hondsrug_layermodel_topsystem()

layer_model

<xarray.Dataset> Size: 64MB
Dimensions:  (layer: 20, y: 200, x: 500)
Coordinates:
  * layer    (layer) int64 160B 1 2 3 4 5 6 7 8 9 ... 12 13 14 15 16 17 18 19 20
  * x        (x) float64 4kB 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 2kB 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 8B 25.0
    dy       float64 8B -25.0
Data variables:
    k        (layer, y, x) float64 16MB nan nan nan nan ... 15.39 15.38 15.44
    idomain  (layer, y, x) float64 16MB -1.0 -1.0 -1.0 -1.0 ... 1.0 1.0 1.0 1.0
    top      (layer, y, x) float64 16MB 6.329 6.329 6.329 ... -16.92 -16.92
    bottom   (layer, y, x) float64 16MB 6.329 6.329 6.329 ... -19.68 -19.68

Furthermore we have planar grid of river, containing a river stage, bed elevation and conductance.

planar_river = imod.data.hondsrug_planar_river()

planar_river

<xarray.Dataset> Size: 1MB
Dimensions:      (x: 500, y: 200)
Coordinates:
  * x            (x) float64 4kB 2.375e+05 2.375e+05 ... 2.5e+05 2.5e+05
  * y            (y) float64 2kB 5.64e+05 5.64e+05 ... 5.59e+05 5.59e+05
    dx           float64 8B ...
    dy           float64 8B ...
    layer        int32 4B 1
Data variables:
    stage        (y, x) float32 400kB nan nan nan 5.014 ... nan nan nan nan
    bottom       (y, x) float32 400kB nan nan nan 3.7 3.7 ... nan nan nan nan
    conductance  (y, x) float32 400kB nan nan nan 8.21 6.02 ... nan nan nan nan

Let’s plot the top elevation of the model on a map. You can see we have a ridge roughly the centre of the model, sided by two low-lying areas.

import numpy as np

imod.visualize.plot_map(
    layer_model["top"].sel(layer=1), "viridis", np.linspace(1, 20, 11)
)

(<Figure size 640x480 with 2 Axes>, <Axes: >)

Let’s plot the river stages on a map. You can see most rivers are located in the low-lying areas.

imod.visualize.plot_map(planar_river["stage"], "viridis", np.linspace(-1, 19, 9))

(<Figure size 640x480 with 2 Axes>, <Axes: >)

Allocate river cells#

Let’s allocate river cells across model layers to cells.

from imod.prepare import ALLOCATION_OPTION, allocate_riv_cells

riv_allocated, _ = allocate_riv_cells(
    allocation_option=ALLOCATION_OPTION.at_elevation,
    active=layer_model["idomain"] == 1,
    top=layer_model["top"],
    bottom=layer_model["bottom"],
    stage=planar_river["stage"],
    bottom_elevation=planar_river["bottom"],
)

Let’s take a look at what we just produced. Since we are dealing with information with depth, it is simplest to make a crosssection plot. For that, we first have to select a crosssection.

import geopandas as gpd
from shapely.geometry import LineString

geometry = LineString([[238725, 561800], [241050, 560350]])

# Define overlay
overlays = [
    {"gdf": gpd.GeoDataFrame(geometry=[geometry]), "edgecolor": "black", "linewidth": 3}
]
# Plot
imod.visualize.plot_map(
    planar_river["stage"], "viridis", np.linspace(-1, 19, 9), overlays
)

(<Figure size 640x480 with 2 Axes>, <Axes: >)

Select a cross section. The plot also requires top and bottom information which we will add first as coordinates before selecting.

riv_allocated.coords["top"] = layer_model["top"]
riv_allocated.coords["bottom"] = layer_model["bottom"]

xsection_allocated = imod.select.cross_section_linestring(riv_allocated, geometry)

xsection_allocated

<xarray.DataArray (layer: 20, s: 150)> Size: 24kB
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]], shape=(20, 150))
Coordinates:
    x        (s) float64 1kB 2.387e+05 2.388e+05 2.388e+05 ... 2.41e+05 2.41e+05
    y        (s) float64 1kB 5.618e+05 5.618e+05 ... 5.604e+05 5.604e+05
    dx       (s) float64 1kB 25.0 15.09 9.914 25.0 ... 25.0 9.914 15.09 25.0
    dy       (s) float64 1kB -15.59 -9.409 -6.183 ... -6.183 -9.409 -15.59
  * layer    (layer) int64 160B 1 2 3 4 5 6 7 8 9 ... 12 13 14 15 16 17 18 19 20
    top      (layer, s) float64 24kB 5.76 5.76 5.76 5.76 ... 2.875 2.875 2.875
    bottom   (layer, s) float64 24kB 5.76 5.76 5.76 5.76 ... 1.679 1.679 1.679
  * s        (s) float64 1kB 14.73 38.35 53.08 ... 2.687e+03 2.702e+03 2.725e+03
    ds       (s) float64 1kB 29.46 17.78 11.68 29.46 ... 29.46 11.68 17.78 29.46

Now that we have selected our data we can plot it.

imod.visualize.cross_section(xsection_allocated, "viridis", [0, 1])

(<Figure size 640x480 with 2 Axes>, <Axes: >)

Let’s plot the locations of our river stages and bottom elevations.

import matplotlib.pyplot as plt

stage_line = imod.select.cross_section_linestring(planar_river["stage"], geometry)
stage_bottom = imod.select.cross_section_linestring(planar_river["bottom"], geometry)

fig, ax = plt.subplots()

imod.visualize.cross_section(xsection_allocated, "viridis", [0, 1], fig=fig, ax=ax)
x_line = stage_line["s"] + stage_line["ds"] / 2
ax.scatter(x_line, stage_line.values, marker=7, c="k")
ax.scatter(x_line, stage_bottom.values, marker=6, c="k")

<matplotlib.collections.PathCollection object at 0x0000016F37D04980>

The above plot might look a bit off. Let’s plot the layer numbers, so that we can identify where model layers are located.

import xarray as xr

layer_grid = layer_model.layer * xr.ones_like(layer_model["top"])
layer_grid.coords["top"] = layer_model["top"]
layer_grid.coords["bottom"] = layer_model["bottom"]
xsection_layer_nr = imod.select.cross_section_linestring(layer_grid, geometry)

imod.visualize.cross_section(xsection_layer_nr, "tab20", np.arange(21))
ax.scatter(x_line, stage_line.values, marker=7, c="k")
ax.scatter(x_line, stage_bottom.values, marker=6, c="k")

<matplotlib.collections.PathCollection object at 0x0000016F40112450>

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()

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.

Distribute conductance#

Next, we’ll take a look at distributing conductances, as there are multiple ways to distribute conductances over layers. For example, it is possible to distribute conductances equally across layers, weighted by layer thickness, or by transmissivity.

from imod.prepare import DISTRIBUTING_OPTION, distribute_riv_conductance

Here’s a map of how the conductances are distributed in our dataset.

imod.visualize.plot_map(
    planar_river["conductance"], "magma", np.logspace(-2, 3, 11), overlays
)

(<Figure size 640x480 with 2 Axes>, <Axes: >)

First compute the allocated river cells for stage to river bottom elevation again. This time we’ll use the stage_to_riv_bot option.

riv_allocated, _ = allocate_riv_cells(
    allocation_option=ALLOCATION_OPTION.stage_to_riv_bot,
    active=layer_model["idomain"] == 1,
    top=layer_model["top"],
    bottom=layer_model["bottom"],
    stage=planar_river["stage"],
    bottom_elevation=planar_river["bottom"],
)

riv_allocated

<xarray.DataArray (layer: 20, y: 200, x: 500)> Size: 2MB
array([[[False, False, False, ...,  True,  True,  True],
        [False, False, False, ..., False,  True,  True],
        [False, False, False, ..., False,  True,  True],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],

       [[False, False, False, ...,  True,  True,  True],
        [False, False, False, ..., False,  True,  True],
        [False, False, False, ..., False,  True,  True],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],

       [[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False,  True],
        ...,
...
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],

       [[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],

       [[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]]],
      shape=(20, 200, 500))
Coordinates:
  * x        (x) float64 4kB 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 2kB 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 8B 25.0
    dy       float64 8B -25.0
  * layer    (layer) int64 160B 1 2 3 4 5 6 7 8 9 ... 12 13 14 15 16 17 18 19 20

Distribute river conductance over model layers. There are multiple options available, which are fully described in imod.prepare.DISTRIBUTING_OPTION(). We can print all possible options as follows:

for option in DISTRIBUTING_OPTION:
    print(option.name)

by_corrected_transmissivity
by_corrected_thickness
equally
by_crosscut_thickness
by_layer_thickness
by_crosscut_transmissivity
by_conductivity
by_layer_transmissivity

To reduce duplicate code, we are going to store all input data in this dictionary which we can provide further as keyword arguments.

distributing_data = {
    "allocated": riv_allocated,
    "conductance": planar_river["conductance"],
    "top": layer_model["top"],
    "bottom": layer_model["bottom"],
    "k": layer_model["k"],
    "stage": planar_river["stage"],
    "bottom_elevation": planar_river["bottom"],
}

Let’s keep things simple first and distribute conductances across layers equally.

riv_conductance = distribute_riv_conductance(
    distributing_option=DISTRIBUTING_OPTION.by_layer_thickness, **distributing_data
)
riv_conductance.coords["top"] = layer_model["top"]
riv_conductance.coords["bottom"] = layer_model["bottom"]

Lets repeat the earlier process to produce a nice cross-section plot.

# Select the conductance over the cross section again.
xsection_distributed = imod.select.cross_section_linestring(riv_conductance, geometry)

fig, ax = plt.subplots()

# Plot grey background of active cells
is_active = ~np.isnan(xsection_distributed.coords["top"])
imod.visualize.cross_section(
    is_active,
    "Greys",
    [0, 1, 2, 3],
    kwargs_colorbar={"plot_colorbar": False},
    fig=fig,
    ax=ax,
)
# Plot conductances
imod.visualize.cross_section(
    xsection_distributed, "magma", np.logspace(-2, 3, 11), fig=fig, ax=ax
)

(<Figure size 640x480 with 2 Axes>, <Axes: >)

Let’s compare the results of all possible options visually. On the top left we’ll plot the hydraulic conductivity, as we haven’t looked at that yet. The other plots show the effects of different settings. Again, distributing options are described in more detail in imod.prepare.DISTRIBUTING_OPTION()

fig, axes = plt.subplots(3, 3, figsize=[13, 11], sharex=True, sharey=True)
axes = np.ravel(axes)

k = distributing_data["k"].copy()
k.coords["top"] = layer_model["top"]
k.coords["bottom"] = layer_model["bottom"]

xsection_k = imod.select.cross_section_linestring(k, geometry)

ax = axes[0]
imod.visualize.cross_section(
    xsection_k,
    "magma",
    np.logspace(-2, 3, 11),
    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("hydraulic conductivity")

for i, option in enumerate(DISTRIBUTING_OPTION, start=1):
    ax = axes[i]
    riv_conductance = distribute_riv_conductance(
        distributing_option=option, **distributing_data
    )
    riv_conductance.coords["top"] = layer_model["top"]
    riv_conductance.coords["bottom"] = layer_model["bottom"]

    xsection_distributed = imod.select.cross_section_linestring(
        riv_conductance, geometry
    )

    if (i % 3) == 2:
        kwargs_colorbar = {"plot_colorbar": True}
    else:
        kwargs_colorbar = {"plot_colorbar": False}

    # Plot grey background of active cells
    is_active = ~np.isnan(xsection_distributed.coords["top"])
    imod.visualize.cross_section(
        is_active,
        "Greys",
        [0, 1, 2, 3],
        kwargs_colorbar={"plot_colorbar": False},
        fig=fig,
        ax=ax,
    )
    # Plot conductances
    imod.visualize.cross_section(
        xsection_distributed,
        "magma",
        np.logspace(-2, 3, 11),
        kwargs_colorbar=kwargs_colorbar,
        fig=fig,
        ax=ax,
    )

    ax.set_title(f"option: {option.name}")

plt.tight_layout()

You can see quite some wildly varying conductances with depth. First, the most simple algorithm equally keeps conductance constant with depth. Second, correcting by_conductivity (hydraulic) decreases the conductance in the deepest layer where rivers occur in the centre of the plot, as conductivity is lower over there. Correcting by_layer_thickness, however, increases conductance in this deep layer, as this is a thicker layer. These differences even out when we correct by_layer_transmissivity (k * thickness). The by_crosscut_thickness algorithm accounts for how far the river bottom penetrates a layer. You can see this reduces conductance in the deep layer compared to distributing by_layer_thickness. by_crosscut_transmissivity uses the crosscut thickness instead of the layer thickness and therefore shows a lower conductance in the deeper layer compared to by_layer_transmissivity. Finally by_corrected_transmissivity and by_corrected_thickness also correct for the displacement of the midpoint over the length where crosscut transmissivity is computed over ([layer top - river bottom]/2) compared to the model cell centre. This further reduces the conductance in the deeper layer.

MODFLOW 6 package#

The data created can now be used to create a MODFLOW 6 package. To construct 3D grids from planar grids for the stages, we can utilize xarrays broadcasting:

from imod.typing.grid import enforce_dim_order

riv_stage = planar_river["stage"].where(riv_allocated)
# Use this function to enforce the right dimension order for iMOD Python.
riv_stage = enforce_dim_order(riv_stage)

riv_stage

<xarray.DataArray 'stage' (layer: 20, y: 200, x: 500)> Size: 8MB
array([[[      nan,       nan,       nan, ..., 1.9508   , 1.9358001,
         1.9258001],
        [      nan,       nan,       nan, ...,       nan, 1.9258001,
         1.9258001],
        [      nan,       nan,       nan, ...,       nan, 1.7708   ,
         1.6558001],
        ...,
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan]],

       [[      nan,       nan,       nan, ..., 1.9508   , 1.9358001,
         1.9258001],
        [      nan,       nan,       nan, ...,       nan, 1.9258001,
         1.9258001],
        [      nan,       nan,       nan, ...,       nan, 1.7708   ,
         1.6558001],
...
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan]],

       [[      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        ...,
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan]]], shape=(20, 200, 500), dtype=float32)
Coordinates:
  * x        (x) float64 4kB 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 2kB 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 8B 25.0
    dy       float64 8B -25.0
  * layer    (layer) int64 160B 1 2 3 4 5 6 7 8 9 ... 12 13 14 15 16 17 18 19 20

We can do the same for the river bottom and construct a river package. Note that we use the previously distributed conductance which we assigned to the variable riv_conducance.

riv_bottom = planar_river["bottom"].where(riv_allocated)
riv_bottom = enforce_dim_order(riv_bottom)

# Remove coordinates that were added for cross-section plots previously
riv_conductance = riv_conductance.drop_vars(["bottom", "top"])

riv = imod.mf6.River(
    stage=riv_stage, conductance=riv_conductance, bottom_elevation=riv_bottom
)

riv

River

<xarray.Dataset> Size: 32MB
Dimensions:           (x: 500, y: 200, layer: 20)
Coordinates:
  * x                 (x) float64 4kB 2.375e+05 2.375e+05 ... 2.5e+05 2.5e+05
  * y                 (y) float64 2kB 5.64e+05 5.64e+05 ... 5.59e+05 5.59e+05
    dx                float64 8B 25.0
    dy                float64 8B -25.0
  * layer             (layer) int64 160B 1 2 3 4 5 6 7 ... 14 15 16 17 18 19 20
Data variables:
    stage             (layer, y, x) float32 8MB nan nan nan nan ... nan nan nan
    conductance       (layer, y, x) float64 16MB nan nan nan nan ... nan nan nan
    bottom_elevation  (layer, y, x) float32 8MB nan nan nan nan ... nan nan nan
    print_input       bool 1B False
    print_flows       bool 1B False
    save_flows        bool 1B False
    observations      object 8B None
    repeat_stress     object 8B None

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

Gallery generated by Sphinx-Gallery

RetroSearch is an open source project built by @garambo | Open a GitHub Issue

Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo

HTML: 3.2 | Encoding: UTF-8 | Version: 0.7.4