Building a model with local grid refinement

This notebook shows how nlmod can be used to create a model with local grid refinement.

import os
import geopandas as gpd

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

Modflow 6 makes it possible to use locally refined grids. In nlmod you can use a shapefile and a number of levels to specify where and how much you want to use local grid refinement. Below we use a shapefile of the Planetenweg in IJmuiden and set the refinement levels at 2. This well create a grid with cells of 100x100m except at the Planetenweg where the cells will be refined to 25x25m. See figures below.

# model settings vertex
model_ws = "03_local_grid_refinement"
model_name = "IJm_planeten"
figdir, cachedir = nlmod.util.get_model_dirs(model_ws)
refine_shp_fname = os.path.abspath(os.path.join("data", "planetenweg_ijmuiden"))
levels = 2
extent = [95000.0, 105000.0, 494000.0, 500000.0]
delr = 100.0
delc = 100.0
steady_state = False
steady_start = True
transient_timesteps = 5
perlen = 1.0
start_time = "2015-1-1"
add_northsea = True
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")

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
[                                        ] | 0% Completed | 136.65 us
[                                        ] | 0% Completed | 101.45 ms
[                                        ] | 0% Completed | 202.22 ms
[                                        ] | 0% Completed | 302.88 ms
[                                        ] | 0% Completed | 403.53 ms
[                                        ] | 0% Completed | 504.16 ms
[                                        ] | 0% Completed | 604.94 ms
[                                        ] | 0% Completed | 705.57 ms
[                                        ] | 0% Completed | 806.21 ms
[                                        ] | 0% Completed | 906.83 ms
[                                        ] | 0% Completed | 1.01 s
[###                                     ] | 8% Completed | 1.11 s
[###                                     ] | 8% Completed | 1.21 s
[######                                  ] | 17% Completed | 1.31 s
[######                                  ] | 17% Completed | 1.41 s
[######                                  ] | 17% Completed | 1.51 s
[##########                              ] | 26% Completed | 1.61 s
[##########                              ] | 26% Completed | 1.71 s
[###################                     ] | 47% Completed | 1.81 s
[###################                     ] | 47% Completed | 1.91 s
[###################                     ] | 47% Completed | 2.01 s
[###################                     ] | 47% Completed | 2.12 s
[###################                     ] | 47% Completed | 2.22 s
[######################                  ] | 56% Completed | 2.32 s
[######################                  ] | 56% Completed | 2.42 s
[##########################              ] | 65% Completed | 2.52 s
[##########################              ] | 65% Completed | 2.62 s
[##########################              ] | 65% Completed | 2.72 s
[#################################       ] | 82% Completed | 2.82 s
[#################################       ] | 82% Completed | 2.92 s
[########################################] | 100% Completed | 3.02 s

# 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.54s/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

Modflow 6 makes it possible to use locally refined grids. In nlmod you can use a shapefile and a number of levels to specify where and how much you want to use local grid refinement. Below we use a shapefile of the Planetenweg in IJmuiden and set the refinement levels at 2. This well create a grid with cells of 100x100m except at the Planetenweg where the cells will be refined to 25x25m. See figures below.

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

# create a model ds by changing grid of 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,
    steady_start=steady_start,
    perlen=[perlen] * (transient_timesteps + 1),
)

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

ds

<xarray.Dataset> Size: 6MB
Dimensions:      (layer: 40, y: 60, x: 100, time: 6)
Coordinates:
  * layer        (layer) <U8 1kB 'AAOP' 'NASC' 'ONAWA' ... 'PZc' 'PZWAz4' 'MSz1'
  * y            (y) float64 480B 5e+05 4.998e+05 ... 4.942e+05 4.94e+05
  * x            (x) float64 800B 9.505e+04 9.515e+04 ... 1.048e+05 1.05e+05
  * time         (time) datetime64[ns] 48B 2015-01-02 2015-01-03 ... 2015-01-07
    spatial_ref  int64 8B 0
Data variables:
    top          (y, x) float64 48kB -18.75 -18.75 -18.75 ... 1.25 0.25 0.25
    botm         (layer, y, x) float64 2MB -18.75 -18.75 ... -218.2 -218.5
    kh           (layer, y, x) float64 2MB 1.0 1.0 1.0 1.0 ... 6.1 6.1 6.11 6.11
    kv           (layer, y, x) float64 2MB 0.1 0.1 0.1 0.1 ... 0.61 0.611 0.611
    area         (y, x) float64 48kB 1e+04 1e+04 1e+04 ... 1e+04 1e+04 1e+04
    steady       (time) int64 48B 1 0 0 0 0 0
    nstp         (time) int64 48B 1 1 1 1 1 1
    tsmult       (time) float64 48B 1.0 1.0 1.0 1.0 1.0 1.0
Attributes:
    extent:      [95000.0, 105000.0, 494000.0, 500000.0]
    gridtype:    structured
    model_name:  IJm_planeten
    mfversion:   mf6
    created_on:  20260513_15:11:59
    exe_name:    /home/docs/checkouts/readthedocs.org/user_builds/nlmod/envs/...
    model_ws:    03_local_grid_refinement
    figdir:      03_local_grid_refinement/figure
    cachedir:    03_local_grid_refinement/cache
    transport:   0

Local grid refinement

The code below applies a local grid refinement to the layer model. The local grid refinement is based on the shapefile ‘planetenweg_ijmuiden.shp’, which contains a line shape of the Planetenweg, and the levels, which is 2. This means that the model cells at the Planetenweg will get a size of 25 x 25m because we halving the cell size twice (100 / (2^2) = 25).

# use gridgen to create vertex grid
ds = nlmod.grid.refine(ds, refinement_features=[(refine_shp_fname, "line", levels)])

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

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

INFO:nlmod.dims.grid.refine:create vertex grid using gridgen
INFO:nlmod.dims.grid.ds_to_gridprops:resample model Dataset to vertex modelgrid
INFO:nlmod.cache._same_function_arguments:cache was created using different function argument 'ds' not in cached arguments, do not use cached data
INFO:nlmod.cache.decorator:caching data -> rws_surface_water.pklz
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
disv = nlmod.gwf.disv(ds, gwf)

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

# 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._disv:creating mf6 DISV
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
da_name = "rws_oppwater"
rws_ds = nlmod.read.rws.discretize_surface_water(
    ds, gdf=gdf_surface_water, da_basename=da_name
)

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

# build ghb package
ghb = nlmod.gwf.ghb(ds, gwf, bhead=f"{da_name}_stage", cond=f"{da_name}_cond")

# 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.ghb:creating mf6 GHB
INFO:nlmod.gwf.gwf.chd:creating mf6 CHD
INFO:nlmod.dims.grid.da_to_reclist:ignore 7304 out of 12640 cells because idomain is inactive

# 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

ds

<xarray.Dataset> Size: 11MB
Dimensions:             (layer: 40, icell2d: 6231, time: 6, iv: 6459, icv: 7)
Coordinates:
  * layer               (layer) <U8 1kB 'AAOP' 'NASC' ... 'PZWAz4' 'MSz1'
  * icell2d             (icell2d) int64 50kB 0 1 2 3 4 ... 6227 6228 6229 6230
    x                   (icell2d) float64 50kB 9.505e+04 9.515e+04 ... 1.05e+05
    y                   (icell2d) float64 50kB 5e+05 5e+05 ... 4.94e+05 4.94e+05
  * time                (time) datetime64[ns] 48B 2015-01-02 ... 2015-01-07
    spatial_ref         int64 8B 0
Dimensions without coordinates: iv, icv
Data variables: (12/20)
    top                 (icell2d) float64 50kB -18.75 -18.75 ... 0.25 0.25
    botm                (layer, icell2d) float64 2MB -18.75 -18.75 ... -218.5
    kh                  (layer, icell2d) float64 2MB 1.0 1.0 1.0 ... 6.11 6.11
    kv                  (layer, icell2d) float64 2MB 0.1 0.1 0.1 ... 0.611 0.611
    steady              (time) int64 48B 1 0 0 0 0 0
    nstp                (time) int64 48B 1 1 1 1 1 1
    ...                  ...
    rws_oppwater_area   (icell2d) float64 50kB 1e+04 1e+04 1e+04 ... 0.0 0.0 0.0
    rws_oppwater_cond   (icell2d) float64 50kB 1e+03 1e+03 1e+03 ... 0.0 0.0 0.0
    rws_oppwater_stage  (icell2d) float64 50kB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    ahn                 (icell2d) float64 50kB nan nan nan ... 0.8357 0.1605
    edge_mask           (layer, icell2d) int64 2MB 1 1 1 1 1 1 1 ... 1 1 1 1 1 1
    recharge            (time, icell2d) float64 299kB 0.0019 0.0019 ... 0.0014
Attributes:
    extent:                  [95000.0, 105000.0, 494000.0, 500000.0]
    gridtype:                vertex
    model_name:              IJm_planeten
    mfversion:               mf6
    created_on:              20260513_15:11:59
    exe_name:                /home/docs/checkouts/readthedocs.org/user_builds...
    model_ws:                03_local_grid_refinement
    figdir:                  03_local_grid_refinement/figure
    cachedir:                03_local_grid_refinement/cache
    transport:               0
    surface_drn_resistance:  10.0

Write and Run

nlmod.sim.write_and_run(sim, ds, write_ds=True)

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 IJm_planeten...
    writing model name file...
    writing package disv...
    writing package npf...
    writing package ic...
    writing package oc...
    writing package ghb...
    writing package drn...
    writing package chd...
    writing package rch...
    writing package ts_0...
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:12:30
 
 Writing simulation list file: mfsim.lst
 Using Simulation name file: mfsim.nam
    Solving:  Stress period:     1    Time step:     1
    Solving:  Stress period:     2    Time step:     1
    Solving:  Stress period:     3    Time step:     1
    Solving:  Stress period:     4    Time step:     1
    Solving:  Stress period:     5    Time step:     1
    Solving:  Stress period:     6    Time step:     1
 
 Run end date and time (yyyy/mm/dd hh:mm:ss): 2026/05/13 15:12:33
 Elapsed run time:  2.813 Seconds
 
 Normal termination of simulation.

Visualise

Using the ds and gwf variables it is quite easy to visualise model data. Below the modelgrid together with the surface water is shown.

plan_weg_gdf = gpd.read_file(refine_shp_fname + ".shp")

# plot modelgrid
ax = nlmod.plot.modelgrid(ds)
nlmod.plot.surface_water(ds, ax=ax)
plan_weg_gdf.plot(ax=ax, color="r", label="Planetenweg")
ax.legend()

# plot zoomed modelgrid
ax = nlmod.plot.modelgrid(ds)
nlmod.plot.surface_water(ds, ax=ax)
ax.set_title("Planetenweg")
plan_weg_gdf.plot(ax=ax, color="r", label="Planetenweg")
ax.set_xlim(100000, 103000)
ax.set_ylim(495000, 497500)
ax.legend();
../_images/af5edeeefb4af1c2f3a9617898711a0909a8a4143b857634af3c9725eaf5a1c5.png ../_images/bade0871aa14ddcf81a0285cee129d26331227d1b59e793eaa822d4d8e9ff148.png

The model dataset of a vertex model differs from a structured model dataset. The data is stored relative to the cell-id instead of the row and column number. Therefore the model dataset has the dimension icell2d instead of the dimensions x and y.

ds

<xarray.Dataset> Size: 11MB
Dimensions:             (layer: 40, icell2d: 6231, time: 6, iv: 6459, icv: 7)
Coordinates:
  * layer               (layer) <U8 1kB 'AAOP' 'NASC' ... 'PZWAz4' 'MSz1'
  * icell2d             (icell2d) int64 50kB 0 1 2 3 4 ... 6227 6228 6229 6230
    x                   (icell2d) float64 50kB 9.505e+04 9.515e+04 ... 1.05e+05
    y                   (icell2d) float64 50kB 5e+05 5e+05 ... 4.94e+05 4.94e+05
  * time                (time) datetime64[ns] 48B 2015-01-02 ... 2015-01-07
    spatial_ref         int64 8B 0
Dimensions without coordinates: iv, icv
Data variables: (12/20)
    top                 (icell2d) float64 50kB -18.75 -18.75 ... 0.25 0.25
    botm                (layer, icell2d) float64 2MB -18.75 -18.75 ... -218.5
    kh                  (layer, icell2d) float64 2MB 1.0 1.0 1.0 ... 6.11 6.11
    kv                  (layer, icell2d) float64 2MB 0.1 0.1 0.1 ... 0.611 0.611
    steady              (time) int64 48B 1 0 0 0 0 0
    nstp                (time) int64 48B 1 1 1 1 1 1
    ...                  ...
    rws_oppwater_area   (icell2d) float64 50kB 1e+04 1e+04 1e+04 ... 0.0 0.0 0.0
    rws_oppwater_cond   (icell2d) float64 50kB 1e+03 1e+03 1e+03 ... 0.0 0.0 0.0
    rws_oppwater_stage  (icell2d) float64 50kB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    ahn                 (icell2d) float64 50kB nan nan nan ... 0.8357 0.1605
    edge_mask           (layer, icell2d) int64 2MB 1 1 1 1 1 1 1 ... 1 1 1 1 1 1
    recharge            (time, icell2d) float64 299kB 0.0019 0.0019 ... 0.0014
Attributes: (12/14)
    extent:                            [95000.0, 105000.0, 494000.0, 500000.0]
    gridtype:                          vertex
    model_name:                        IJm_planeten
    mfversion:                         mf6
    created_on:                        20260513_15:11:59
    exe_name:                          /home/docs/checkouts/readthedocs.org/u...
    ...                                ...
    cachedir:                          03_local_grid_refinement/cache
    transport:                         0
    surface_drn_resistance:            10.0
    model_dataset_written_to_disk_on:  20260513_15:12:27
    model_data_written_to_disk_on:     20260513_15:12:30
    model_ran_on:                      20260513_15:12:33

To plot the same rasters as for the previous model we can use the nlmod.plot.data_array() function.

fig, axes = nlmod.plot.get_map(extent, nrows=2, ncols=2, figsize=14)

nlmod.plot.data_array(ds["ahn"], ds, ax=axes[0][0])
nlmod.plot.data_array(ds["botm"][0], ds, ax=axes[0][1])
nlmod.plot.data_array(nlmod.layers.get_idomain(ds)[0], ds, ax=axes[1][0])
nlmod.plot.data_array(ds["edge_mask"][0], ds, ax=axes[1][1])

fig, axes = nlmod.plot.get_map(extent, nrows=2, ncols=2, figsize=(14, 11))
nlmod.plot.data_array(ds["bathymetry"], ds, ax=axes[0][0])
nlmod.plot.data_array(ds["northsea"], ds, ax=axes[0][1])
nlmod.plot.data_array(ds["kh"][1], ds, ax=axes[1][0])
nlmod.plot.data_array(ds["recharge"][0], ds, ax=axes[1][1]);
../_images/36433ee03c871c9c3367c92daaf0e8f64fbd9669519948e95c1f2548640df8a2.png ../_images/80fca7ee46c228c195bc112803ea431c24e1b521ed492896ee29668f7698b15c.png

We can save the entire model as a UGRID NetCDF-file. This can be opened in qgis, as a ‘Mesh Layer’. For more information see https://docs.qgis.org/3.16/en/docs/user_manual/working_with_mesh/mesh_properties.html

fname = os.path.join(ds.figdir, "results.nc")
nlmod.gis.ds_to_ugrid_nc_file(ds, fname)

<xarray.Dataset> Size: 11MB
Dimensions:                 (icell2d: 6231, time: 6, iv: 6459, icv: 6)
Coordinates:
  * icell2d                 (icell2d) int64 50kB 0 1 2 3 ... 6227 6228 6229 6230
    x                       (icell2d) float64 50kB 9.505e+04 ... 1.05e+05
    y                       (icell2d) float64 50kB 5e+05 5e+05 ... 4.94e+05
  * time                    (time) datetime64[ns] 48B 2015-01-02 ... 2015-01-07
    spatial_ref             int64 8B 0
Dimensions without coordinates: iv, icv
Data variables: (12/216)
    top                     (icell2d) float64 50kB -18.75 -18.75 ... 0.25 0.25
    botm_AAOP               (icell2d) float64 50kB -18.75 -18.75 ... -0.25 -0.25
    botm_NASC               (icell2d) float64 50kB -18.75 -18.75 ... -0.25 -0.25
    botm_ONAWA              (icell2d) float64 50kB -18.75 -18.75 ... -1.75 -1.75
    botm_NAZA               (icell2d) float64 50kB -18.75 -18.75 ... -4.25 -4.75
    botm_NAWA               (icell2d) float64 50kB -18.75 -18.75 ... -4.25 -4.75
    ...                      ...
    edge_mask_PZWAz3        (icell2d) int64 50kB 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1
    edge_mask_PZc           (icell2d) int64 50kB 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1
    edge_mask_PZWAz4        (icell2d) int64 50kB 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1
    edge_mask_MSz1          (icell2d) int64 50kB 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1
    recharge                (time, icell2d) float64 299kB 0.0019 ... 0.0014
    mesh_topology           int64 8B 0
Attributes: (12/14)
    extent:                            [95000.0, 105000.0, 494000.0, 500000.0]
    gridtype:                          vertex
    model_name:                        IJm_planeten
    mfversion:                         mf6
    created_on:                        20260513_15:11:59
    exe_name:                          /home/docs/checkouts/readthedocs.org/u...
    ...                                ...
    cachedir:                          03_local_grid_refinement/cache
    transport:                         0
    surface_drn_resistance:            10.0
    model_dataset_written_to_disk_on:  20260513_15:12:27
    model_data_written_to_disk_on:     20260513_15:12:30
    model_ran_on:                      20260513_15:12:33