Modifying model layers

This notebook shows methods for modyfying model layers for MODFLOW models, for example by combining layers and splitting layers. Multiple layers can be combined into one layer or one layer can be split into sub-layers based on a fraction of the original thickness.

import matplotlib.pyplot as plt
import numpy as np
from shapely.geometry import LineString

import nlmod
from nlmod.plot import DatasetCrossSection

nlmod.util.get_color_logger("INFO")
nlmod.show_versions()

Python version     : 3.11.14
NumPy version      : 2.4.4
Xarray version     : 2026.4.0
Matplotlib version : 3.10.9
Flopy version      : 3.10.0

nlmod version      : 0.11.3dev

def compare_layer_models(
    ds1,
    line,
    colors,
    ds2=None,
    zmin=-200,
    zmax=10,
    min_label_area=1000,
    title1="REGIS original",
    title2="Modified layers",
    xlabel="Distance along x-sec (m)",
    ylabel="m NAP",
):
    if ds2 is None:
        fig, ax1 = plt.subplots(1, 1, figsize=(14, 6))
    else:
        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 12), sharex=True)
    dcs1 = DatasetCrossSection(ds1, line=line, ax=ax1, zmin=zmin, zmax=zmax)
    polys2 = dcs1.plot_layers(colors=colors, min_label_area=min_label_area)
    dcs1.plot_grid(linewidth=0.5, vertical=False)
    ax1.set_ylabel(ylabel)

    if ds2 is not None:
        ax1.set_title(title1)
        dcs2 = DatasetCrossSection(ds2, line=line, ax=ax2, zmin=zmin, zmax=zmax)
        polys1 = dcs2.plot_layers(colors=colors, min_label_area=min_label_area)
        dcs2.plot_grid(linewidth=0.5, vertical=False)
        ax2.set_ylabel(ylabel)
        ax2.set_xlabel(xlabel)
        ax2.set_title(title2)
    else:
        ax1.set_xlabel(xlabel)

Set new model top

The nlmod.layers.set_model_top changes the top of the model. When the new top is lower than the old top, the new top is burned in the layer model, lowering the top of the top layer(s). Top layers can become incactive, when the thickness is reduced to 0. When the new top is higher than the old top, the thickness of the most upper active layer (not necessarily the first) is increased. This method can be used to change the model top to a digital terrain model with a higher accuracy.

First transform the regis-date to a model Dataset, as the next methods need a model Dataset.

ds = nlmod.to_model_ds(regis)

INFO:nlmod.dims.base.to_model_ds:resample layer model data to structured modelgrid
INFO:nlmod.dims.layers.get_kh_kv:kv and kh both undefined in layer HLc
INFO:nlmod.dims.layers.get_kh_kv:kv and kh both undefined in layer DTc
INFO:nlmod.dims.layers._fill_var:Filling 11886 values in active cells of kh by multipying kv with an anisotropy of 10
INFO:nlmod.dims.layers._fill_var:Filling 51924 values in active cells of kv by dividing kh by an anisotropy of 10
INFO:nlmod.dims.layers._fill_var:Filling 2485 values in active cells of kh with a value of 1.0 m/day
INFO:nlmod.dims.layers._fill_var:Filling 2485 values in active cells of kv with a value of 0.1 m/day

ds_new = nlmod.layers.set_model_top(ds.copy(deep=True), 5.0)
compare_layer_models(ds, line, colors, ds2=ds_new, title2="New top")

../_images/dfd87c8ef403cc9c0e858b5ca9c3483c1daaf597141792a35cf29f6992a4f1ad.png

Set layer top

nlmod.layers.set_layer_top sets the layer top to a specified value or array.

This method only changes the shape of the layer, and does not check if all hydrological properties are defined for cells that had a thickness of 0 before.

ds_new = nlmod.layers.set_layer_top(ds.copy(deep=True), "WAk1", -40.0)
compare_layer_models(ds, line, colors, ds2=ds_new, title2="Modified")

../_images/23b5a6773f6ef0bcf52c58ed086889ecc89277e7b5972246594bc046cec17c2f.png

Set layer bottom

nlmod.layers.set_layer_botm sets the layer botm to a specified value or array.

This method only changes the shape of the layer, and does not check if all hydrological properties are defined for cells that had a thickness of 0 before.

# set the botm of 'WAk1' to -70 m NAP
ds_new = nlmod.layers.set_layer_botm(ds.copy(deep=True), "WAk1", -70.0)
compare_layer_models(ds, line, colors, ds2=ds_new, title2="Modified")

../_images/2e9503752b9356a4dad710b88a4098c1727df7f0f0c8936c3102cbd7b90c749f.png

Set layer thickness

nlmod.layers.set_layer_thickness sets the thickness of a layer to a specified value or array. With a parameter called ‘change’ you can specify in which direction the layer is changed. The only supported option for now is ‘botm’, which changes the layer botm.

This method only changes the shape of the layer, and does not check if all hydrological properties are defined for cells that had a thickness of 0 before.

# set the thickness of 'WAk1' to 20 m NAP
ds_new = nlmod.layers.set_layer_thickness(ds.copy(deep=True), "WAk1", 20)
compare_layer_models(ds, line, colors, ds2=ds_new, title2="Modified")

../_images/9b7eddb5ec9fca93b7c874ccbbcd773638dcd0410db6c38b4627111b403004d0.png

Set minimum layer thickness

nlmod.layers.set_minimum layer_thickness increases the thickness of a layer if the thickness is less than a specified value. With a parameter called ‘change’ you can specify in which direction the layer is changed. The only supported option for now is ‘botm’, which changes the layer botm.

This method only changes the shape of the layer, and does not check if all hydrological properties are defined for cells that had a thickness of 0 before.

# set the minimum thickness of 'PZWAz2' to 20 m
ds_new = nlmod.layers.set_minimum_layer_thickness(ds.copy(deep=True), "PZWAz2", 20.0)
compare_layer_models(ds, line, colors, ds2=ds_new, title2="Modified")

../_images/b9574e40784bef1ca314a82424e4cbb7474188077f3b97e3968d1f9010d25400.png

Modifying model layers

Get data

Split layers

Combine layers

Set new model top

Set layer top

Set layer bottom

Set layer thickness

Set minimum layer thickness