Building a model anywhere in the Netherlands

This example notebook shows a basic example of a model created using online data with the nlmod package. nlmod contains functions to create modflow models anywhere in the Netherlands.

import nlmod

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

Model settings

We create a modflow model with the name ‘IJmuiden’. This model has the following properties:

an extent that covers part of the Northsea, Noordzeekanaal and the small port city IJmuiden.
a structured grid with cells of 100 x 100 m2
the model is a steady state model with a single time step.
starting heads of 1 m NAP in every cell.

# model settings
model_ws = "01_basic_model"
model_name = "IJmuiden"
figdir, cachedir = nlmod.util.get_model_dirs(model_ws)
extent = [95000.0, 105000.0, 494000.0, 500000.0]
delr = 100.0
delc = 100.0
steady_state = True
start_time = "2015-1-1"
starting_head = 1.0

Download data

First, download the data needed for the model from online sources. It is recommended to store a copy somewhere in case the online dataset changes. For this model we download data from these sources:

regis (layer model)
geotop (layer model)
rijkswaterstaat (surface water)
jarkus (bathymetry)
ahn (digital elevation model)
knmi (precipitation and evaporation)

# layer models
regis_ds = nlmod.read.regis.download_regis(extent, botm_layer="MSz1",
                                            cachedir=cachedir,
                                            cachename="regis.nc")

geotop_ds = nlmod.read.geotop.download_geotop(extent,
                                              cachedir=cachedir,
                                              cachename="geotop.nc",
                                              chunks=None)

WARNING:nlmod.dims.layers.remove_layer_dim_from_top:Botm of layer is not equal to top of deeper layer in 3 cells
INFO:nlmod.cache.wrapper:caching data -> regis.nc
INFO:nlmod.cache.wrapper:caching data -> geotop.nc

# read surface water data
gdf_surface_water = nlmod.read.rws.get_gdf_surface_water(
    extent=extent, cachedir=cachedir, cachename="rws_surface_water.pklz"
)

# bathymetry
bathymetry_da = nlmod.read.jarkus.download_bathymetry(extent=extent,
                                                      cachedir=cachedir,
                                                      cachename="bathymetry.nc")

INFO:nlmod.cache.decorator:caching data -> rws_surface_water.pklz
INFO:nlmod.cache.wrapper:caching data -> bathymetry.nc

# Digital elevation model
ahn_da = nlmod.read.ahn.download_ahn(extent, cachedir=cachedir, cachename="ahn")

INFO:nlmod.cache.wrapper:caching data -> ahn.nc

Downloading tiles of AHN4:   0%|          | 0/2 [00:00<?, ?it/s]
Downloading tiles of AHN4: 100%|██████████| 2/2 [00:05<00:00,  2.78s/it]

# knmi data
oc_knmi = nlmod.read.knmi.download_knmi(extent=extent, delr=delr, delc=delc, start='2000-1-1', end='2020-1-1')

INFO:nlmod.dims.base.to_model_ds:resample layer model data to structured modelgrid
INFO:hydropandas.io.knmi.get_knmi_obs:get data from station 201 and variable RD from 2000-01-01 to 2020-01-01
INFO:hydropandas.io.knmi.fill_missing_measurements:station 201 has no measurements between 2000-01-01 00:00:00 and 2020-01-01 00:00:00 trying to get measurements from nearest station -> 226
INFO:hydropandas.io.knmi.get_knmi_obs:get data from station 226 and variable RD from 2000-01-01 to 2020-01-01
INFO:hydropandas.io.knmi.get_knmi_obs:get data from station 231 and variable RD from 2000-01-01 to 2020-01-01
INFO:hydropandas.io.knmi.fill_missing_measurements:station 231 has no measurements between 2000-01-01 00:00:00 and 2020-01-01 00:00:00 trying to get measurements from nearest station -> 210
INFO:hydropandas.io.knmi._add_missing_indices:station 210 has no measurements before 2000-01-14 09:00:00
INFO:hydropandas.io.knmi._add_missing_indices:station 210 has no measurements after 2006-02-27 09:00:00
INFO:hydropandas.io.knmi.fill_missing_measurements:Filled 5048 observations from station 263 ASSENDELFT -> 210 BEVERWIJK
INFO:hydropandas.io.knmi.fill_missing_measurements:Filled 171 observations from station 226 WIJK-AAN-ZEE -> 210 BEVERWIJK
INFO:hydropandas.io.knmi.get_knmi_obs:get data from station 210 and variable RD from 2000-01-01 to 2020-01-01
INFO:hydropandas.io.knmi._add_missing_indices:station 210 has no measurements before 2000-01-14 09:00:00
INFO:hydropandas.io.knmi._add_missing_indices:station 210 has no measurements after 2006-02-27 09:00:00
INFO:hydropandas.io.knmi.fill_missing_measurements:Filled 5048 observations from station 263 ASSENDELFT -> 210 BEVERWIJK
INFO:hydropandas.io.knmi.fill_missing_measurements:Filled 171 observations from station 226 WIJK-AAN-ZEE -> 210 BEVERWIJK
INFO:hydropandas.io.knmi.get_knmi_obs:get data from station 254 and variable RD from 2000-01-01 to 2020-01-01
INFO:hydropandas.io.knmi.fill_missing_measurements:station 254 has no measurements between 2000-01-01 00:00:00 and 2020-01-01 00:00:00 trying to get measurements from nearest station -> 210
INFO:hydropandas.io.knmi._add_missing_indices:station 210 has no measurements before 2000-01-14 09:00:00
INFO:hydropandas.io.knmi._add_missing_indices:station 210 has no measurements after 2006-02-27 09:00:00
INFO:hydropandas.io.knmi.fill_missing_measurements:Filled 5048 observations from station 263 ASSENDELFT -> 210 BEVERWIJK
INFO:hydropandas.io.knmi.fill_missing_measurements:Filled 171 observations from station 226 WIJK-AAN-ZEE -> 210 BEVERWIJK
INFO:hydropandas.io.knmi.get_knmi_obs:get data from station 257 and variable EV24 from 2000-01-01 to 2020-01-01
INFO:hydropandas.io.knmi._add_missing_indices:station 257 has no measurements before 2001-05-03 01:00:00
INFO:hydropandas.io.knmi.fill_missing_measurements:Filled 488 observations from station 240 SCHIPHOL -> 257 WIJK-AAN-ZEE

Create model

Now create the modflow model:

Build a layer model using Regis and Geotop. Use MSz1 as the bottom layer; all deeper Regis layers are excluded. Use nlmod.read.regis.get_layer_names() to list available layers.
Replace the Holocene layer in Regis with Geotop due to missing hydraulic properties.
Extend Regis/Geotop layers into the North Sea using shoreline extrapolation and Jarkus bathymetry data.
Add large water bodies (e.g. North Sea, IJsselmeer) as general head boundaries.
Add surface drainage using the AHN with a conductance of 1000 m²/d.
Calculate recharge from precipitation and evaporation time series
Add constant head boundaries on all model edges using the starting head.

# create layer model
layer_model = nlmod.read.regis.add_geotop_to_regis_layers(regis_ds,
                                                          geotop_ds,
                                                          layers="HLc")

# create a dataset from they layer model
ds = nlmod.to_model_ds(layer_model, model_name, model_ws, delr=delr, delc=delc)

# add time discretisation
ds = nlmod.time.set_ds_time(ds, start=start_time, steady=steady_state, perlen=365 * 5)

# add northsea and bathymetry to modelgrid
ds.update(nlmod.read.rws.discretize_northsea(ds, gdf=gdf_surface_water))
ds.update(nlmod.read.jarkus.discretize_bathymetry(ds, bathymetry_da.drop_vars("time"), cachedir=cachedir, cachename="bathymetry.nc"))
ds = nlmod.dims.add_bathymetry_to_layer_model(ds)

INFO:nlmod.read.geotop.add_kh_and_kv:Determining kh and kv of geotop-data based on lithoclass
INFO:nlmod.read.geotop.add_kh_and_kv:Setting kv equal to kh / 1.0
INFO:nlmod.dims.base.to_model_ds:resample layer model data to structured modelgrid
INFO:nlmod.dims.layers._fill_var:Filling 8032 values in active cells of kh by multipying kv with an anisotropy of 10
INFO:nlmod.dims.layers._fill_var:Filling 69674 values in active cells of kv by dividing kh by an anisotropy of 10
INFO:nlmod.cache._same_function_arguments:cache was created using different function argument 'arg1' not in cached arguments, do not use cached data
INFO:nlmod.dims.resample.structured_da_to_ds:No crs in da. Setting crs equal to ds: EPSG:28992
INFO:nlmod.cache.wrapper:caching data -> bathymetry.nc
INFO:nlmod.dims.layers.add_bathymetry_to_layer_model:Filling NaN values in top/botm and kh/kv in North Sea using bathymetry data from jarkus

# create simulation
sim = nlmod.sim.sim(ds)

# create time discretisation
tdis = nlmod.sim.tdis(ds, sim)

# create ims
ims = nlmod.sim.ims(sim)

# create groundwater flow model
gwf = nlmod.gwf.gwf(ds, sim)

# Create discretization
dis = nlmod.gwf.dis(ds, gwf)

# create node property flow
npf = nlmod.gwf.npf(ds, gwf)

# Create the initial conditions package
ic = nlmod.gwf.ic(ds, gwf, starting_head=starting_head)

# Create the output control package
oc = nlmod.gwf.oc(ds, gwf)

INFO:nlmod.sim.sim.sim:creating mf6 SIM
INFO:nlmod.sim.sim.tdis:creating mf6 TDIS
INFO:nlmod.sim.sim.ims:creating mf6 IMS
INFO:nlmod.gwf.gwf.gwf:creating mf6 GWF
INFO:nlmod.gwf.gwf._dis:creating mf6 DIS
INFO:nlmod.gwf.gwf.npf:creating mf6 NPF
INFO:nlmod.gwf.gwf.ic:creating mf6 IC
INFO:nlmod.gwf.gwf.ic:adding 'starting_head' data array to ds
INFO:nlmod.gwf.gwf.oc:creating mf6 OC

# discretize surface water bodies from geodataframe
rws_ds = nlmod.read.rws.discretize_surface_water(
    ds, gdf=gdf_surface_water, da_basename="rws_oppwater"
)

# add data to model dataset
ds.update(rws_ds)

# build ghb package
ghb = nlmod.gwf.ghb(ds, gwf, bhead="rws_oppwater_stage", cond="rws_oppwater_cond")

INFO:nlmod.gwf.gwf.ghb:creating mf6 GHB

# discretize ahn
ahn_ds = nlmod.read.ahn.discretize_ahn(ds, ahn_da)

# add data to model dataset
ds.update(ahn_ds)

# build surface level drain package
drn = nlmod.gwf.surface_drain_from_ds(ds, gwf, resistance=10.0)

# add constant head cells at model boundaries
ds.update(nlmod.grid.mask_model_edge(ds))
chd = nlmod.gwf.chd(ds, gwf, mask="edge_mask", head="starting_head")

INFO:nlmod.gwf.gwf.chd:creating mf6 CHD

# discretize knmi recharge data
knmi_ds = nlmod.read.knmi.discretize_knmi(ds, oc_knmi, cachedir=cachedir, cachename="recharge")

# update model dataset
ds.update(knmi_ds)

# create recharge package
rch = nlmod.gwf.rch(ds, gwf)

WARNING:nlmod.read.knmi.discretize_knmi:The default of hourly_precision=False will be changed to True in a future version of nlmod. Pass hourly_precision=False to retain current behavior or hourly_precision=True to adopt the future default and silence this warning.
INFO:nlmod.cache.wrapper:caching data -> recharge.nc
INFO:nlmod.gwf.gwf.rch:creating mf6 RCH

Model data is stored in the variable ds which is an xarray.Dataset.

Write and Run

Write the model files and run the model using the function nlmod.sim.write_and_run) as shown below. This function has two additional options:

Saving the model dataset (write_ds=True) for faster future loading.
Saving a copy of the current notebook (nb_path=”<notebook_name>.ipynb”) alongside the model files.

nlmod.sim.write_and_run(sim, ds, write_ds=True, script_path="01_basic_model.ipynb")

INFO:nlmod.sim.sim.write_and_run:write script 20260513_01_basic_model.ipynb to model workspace
INFO:nlmod.sim.sim.write_and_run:write model dataset to cache
INFO:nlmod.sim.sim.write_and_run:write modflow files to model workspace
writing simulation...
  writing simulation name file...
  writing simulation tdis package...
  writing solution package ims...
  writing model IJmuiden...
    writing model name file...
    writing package dis...
    writing package npf...
    writing package ic...
    writing package oc...
    writing package ghb...
    writing package drn...
    writing package chd...
    writing package rch...
INFO:nlmod.sim.sim.write_and_run:run model
FloPy is using the following executable to run the model: ../../../../../envs/latest/lib/python3.11/site-packages/nlmod/bin/mf6
                                   MODFLOW 6
                U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
                            VERSION 6.6.3 09/29/2025

   MODFLOW 6 compiled Oct 07 2025 22:51:46 with Intel(R) Fortran Intel(R) 64
   Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0
                             Build 20220726_000000

This software has been approved for release by the U.S. Geological 
Survey (USGS). Although the software has been subjected to rigorous 
review, the USGS reserves the right to update the software as needed 
pursuant to further analysis and review. No warranty, expressed or 
implied, is made by the USGS or the U.S. Government as to the 
functionality of the software and related material nor shall the 
fact of release constitute any such warranty. Furthermore, the 
software is released on condition that neither the USGS nor the U.S. 
Government shall be held liable for any damages resulting from its 
authorized or unauthorized use. Also refer to the USGS Water 
Resources Software User Rights Notice for complete use, copyright, 
and distribution information.

 
 MODFLOW runs in SEQUENTIAL mode
 
 Run start date and time (yyyy/mm/dd hh:mm:ss): 2026/05/13 15:11:24
 
 Writing simulation list file: mfsim.lst
 Using Simulation name file: mfsim.nam
    Solving:  Stress period:     1    Time step:     1
 
 Run end date and time (yyyy/mm/dd hh:mm:ss): 2026/05/13 15:11:24
 Elapsed run time:  0.716 Seconds
 
 Normal termination of simulation.

Visualise

Plot the modelgrid and surface water

ax = nlmod.plot.modelgrid(ds)
nlmod.plot.surface_water(ds, ax=ax)

<Axes: >

../_images/ce6c7e868ee4a54521e337007c328924f76d043643b609875243af10498f61bc.png

Plot model input data.

fig, axes = nlmod.plot.get_map(ds.extent, nrows=2, ncols=2, figsize=14)
ds["ahn"].plot(ax=axes[0][0])
ds["botm"][0].plot(ax=axes[0][1])
nlmod.layers.get_idomain(ds)[0].plot(ax=axes[1][0])
ds["edge_mask"][0].plot(ax=axes[1][1])

fig, axes = nlmod.plot.get_map(ds.extent, nrows=2, ncols=2, figsize=14)
ds["bathymetry"].plot(ax=axes[0][0])
ds["northsea"].plot(ax=axes[0][1])
ds["kh"][1].plot(ax=axes[1][0])
ds["recharge"].plot(ax=axes[1][1]);

../_images/307e15604f2f1f60766a3b8d1fd8ac7994e58d4d8e30520a0a90f226b6c75340.png

../_images/ce7279f8b92d315d8d099e7736c1d6cafdf4af089aedf49272fa1b9a14afcf89.png