diff --git a/docs/examples/17_unsaturated_zone_flow.ipynb b/docs/examples/17_unsaturated_zone_flow.ipynb
index 9dbd49eb..68c675ad 100644
--- a/docs/examples/17_unsaturated_zone_flow.ipynb
+++ b/docs/examples/17_unsaturated_zone_flow.ipynb
@@ -38,6 +38,18 @@
     "nlmod.show_versions()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "34f28281-fa68-4618-b063-8b9604306974",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Ignore warning about 1d layers, in numpy 1.26.1 and flopy 3.4.3\n",
+    "import warnings\n",
+    "warnings.filterwarnings(\"ignore\", message=\"Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future.\")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "9813a7da",
@@ -241,7 +253,7 @@
    "outputs": [],
    "source": [
     "# run model\n",
-    "nlmod.sim.write_and_run(sim, ds)\n",
+    "nlmod.sim.write_and_run(sim, ds, silent=True)\n",
     "\n",
     "# get heads\n",
     "head_rch = nlmod.gwf.get_heads_da(ds)\n",
@@ -311,7 +323,7 @@
    "outputs": [],
    "source": [
     "# run model\n",
-    "nlmod.sim.write_and_run(sim, ds)\n",
+    "nlmod.sim.write_and_run(sim, ds, silent=True)\n",
     "\n",
     "# get heads\n",
     "head_rch = nlmod.gwf.get_heads_da(ds)\n",
diff --git a/nlmod/dims/layers.py b/nlmod/dims/layers.py
index 4e58d775..c3a6e4b5 100644
--- a/nlmod/dims/layers.py
+++ b/nlmod/dims/layers.py
@@ -689,6 +689,8 @@ def set_model_top(ds, top, min_thickness=0.0):
                 "Either add attribute or use nlmod functions to build the dataset."
             )
         )
+    if "layer" in ds["top"].dims:
+        raise (ValueError("set_model_top does not support top with a layer dimension"))
     if isinstance(top, (float, int)):
         top = xr.full_like(ds["top"], top)
     if not top.shape == ds["top"].shape:
@@ -712,6 +714,8 @@ def set_model_top(ds, top, min_thickness=0.0):
 def set_layer_top(ds, layer, top):
     """Set the top of a layer."""
     assert layer in ds.layer
+    if "layer" in ds["top"].dims:
+        raise (ValueError("set_layer_top does not support top with a layer dimension"))
     lay = np.where(ds.layer == layer)[0][0]
     if lay == 0:
         # change the top of the model
@@ -735,6 +739,8 @@ def set_layer_top(ds, layer, top):
 def set_layer_botm(ds, layer, botm):
     """Set the bottom of a layer."""
     assert layer in ds.layer
+    if "layer" in ds["top"].dims:
+        raise (ValueError("set_layer_botm does not support top with a layer dimension"))
     lay = np.where(ds.layer == layer)[0][0]
     # if lay > 0 and np.any(botm > ds["botm"][lay - 1]):
     #    raise (Exception("set_layer_botm cannot change botm of higher layers yet"))
@@ -782,8 +788,11 @@ def set_minimum_layer_thickness(ds, layer, min_thickness, change="botm"):
 def remove_thin_layers(ds, min_thickness=0.1):
     """Remove layers from cells with a thickness less than min_thickness
 
-    THe thickness of the removed cells is added to the first active layer below
+    The thickness of the removed cells is added to the first active layer below
     """
+    if "layer" in ds["top"].dims:
+        msg = "remove_thin_layers does not support top with a layer dimension"
+        raise (ValueError(msg))
     thickness = calculate_thickness(ds)
     for lay_org in range(len(ds.layer)):
         # determine where the layer is too thin
diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py
index 1768b9fe..50a17cfe 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -446,7 +446,7 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
     Parameters
     ----------
     da : xarray.DataArray
-        THe data-array to be resampled, with dimensions x and y.
+        The data-array to be resampled, with dimensions x and y.
     ds : xarray.Dataset
         The model dataset.
     method : string or rasterio.enums.Resampling, optional
@@ -457,7 +457,7 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
         method is 'linear' or 'nearest' da.interp() is used. Otherwise
         da.rio.reproject_match() is used. The default is "average".
     nodata : float, optional
-        THe nodata value in input and output. THe default is np.NaN.
+        The nodata value in input and output. The default is np.NaN.
 
     Returns
     -------
diff --git a/nlmod/dims/time.py b/nlmod/dims/time.py
index 4bbf6015..b7906d9d 100644
--- a/nlmod/dims/time.py
+++ b/nlmod/dims/time.py
@@ -503,7 +503,9 @@ def dataframe_to_flopy_timeseries(
     filename=None,
     time_series_namerecord=None,
     interpolation_methodrecord="stepwise",
+    append=False,
 ):
+    assert not df.isna().any(axis=None)
     if ds is not None:
         # set index to days after the start of the simulation
         df = df.copy()
@@ -519,7 +521,10 @@ def dataframe_to_flopy_timeseries(
 
     if isinstance(interpolation_methodrecord, str):
         interpolation_methodrecord = [interpolation_methodrecord] * len(df.columns)
-    package.ts.initialize(
+
+    # initialize or append a new package
+    method = package.ts.append_package if append else package.ts.initialize
+    method(
         filename=filename,
         timeseries=timeseries,
         time_series_namerecord=time_series_namerecord,
diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py
index 9a2ed383..abc3d30e 100644
--- a/nlmod/gwf/gwf.py
+++ b/nlmod/gwf/gwf.py
@@ -215,10 +215,6 @@ def _disv(ds, model, length_units="METERS", pname="disv", **kwargs):
         xorigin = ds.attrs["xorigin"]
         yorigin = ds.attrs["yorigin"]
         angrot = ds.attrs["angrot"]
-    elif "extent" in ds.attrs.keys():
-        xorigin = ds.attrs["extent"][0]
-        yorigin = ds.attrs["extent"][2]
-        angrot = 0.0
     else:
         xorigin = 0.0
         yorigin = 0.0
diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py
index 1cf237ab..aed11a13 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -1015,7 +1015,13 @@ def gdf_to_seasonal_pkg(
         **kwargs,
     )
     # add timeseries for the seasons 'winter' and 'summer'
-    add_season_timeseries(ds, package, summer_months=summer_months, seasons=seasons)
+    add_season_timeseries(
+        ds,
+        package,
+        summer_months=summer_months,
+        winter_name="winter",
+        summer_name="summer",
+    )
     return package
 
 
@@ -1024,8 +1030,26 @@ def add_season_timeseries(
     package,
     summer_months=(4, 5, 6, 7, 8, 9),
     filename="season.ts",
-    seasons=("winter", "summer"),
+    winter_name="winter",
+    summer_name="summer",
 ):
+    """Add time series indicating which season is active (e.g. summer/winter).
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        xarray dataset used for time discretization
+    package : flopy.mf6 package
+        Modflow 6 package to add time series to
+    summer_months : tuple, optional
+        summer months. The default is (4, 5, 6, 7, 8, 9), so from april to september.
+    filename : str, optional
+        name of time series file. The default is "season.ts".
+    winter_name : str, optional
+        The name of the time-series with ones in winter. The default is "winter".
+    summer_name : str, optional
+        The name of the time-series with ones in summer. The default is "summer".
+    """
     tmin = pd.to_datetime(ds.time.start)
     if tmin.month in summer_months:
         ts_data = [(0.0, 0.0, 1.0)]
@@ -1048,20 +1072,14 @@ def add_season_timeseries(
     package.ts.initialize(
         filename=filename,
         timeseries=ts_data,
-        time_series_namerecord=seasons,
+        time_series_namerecord=[winter_name, summer_name],
         interpolation_methodrecord=["stepwise", "stepwise"],
     )
 
 
 def rivdata_from_xylist(gwf, xylist, layer, stage, cond, rbot):
-    # TODO: temporary fix until flopy is patched
-    if gwf.modelgrid.grid_type == "structured":
-        gi = flopy.utils.GridIntersect(gwf.modelgrid, rtree=False)
-        cellids = gi.intersect(xylist, shapetype="linestring")["cellids"]
-    else:
-        gi = flopy.utils.GridIntersect(gwf.modelgrid)
-        cellids = gi.intersects(xylist, shapetype="linestring")["cellids"]
-
+    gi = flopy.utils.GridIntersect(gwf.modelgrid, method="vertex")
+    cellids = gi.intersect(xylist, shapetype="linestring")["cellids"]
     riv_data = []
     for cid in cellids:
         if len(cid) == 2:
diff --git a/nlmod/modpath/modpath.py b/nlmod/modpath/modpath.py
index 4b3a4588..20db899e 100644
--- a/nlmod/modpath/modpath.py
+++ b/nlmod/modpath/modpath.py
@@ -430,7 +430,15 @@ def pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5):
 
 
 def sim(
-    mpf, particlegroups, direction="backward", gwf=None, ref_time=None, stoptime=None
+    mpf,
+    particlegroups,
+    direction="backward",
+    gwf=None,
+    ref_time=None,
+    stoptime=None,
+    simulationtype="combined",
+    weaksinkoption="pass_through",
+    weaksourceoption="pass_through",
 ):
     """Create a modpath backward simulation from a particle group.
 
@@ -472,10 +480,10 @@ def sim(
 
     mpsim = flopy.modpath.Modpath7Sim(
         mpf,
-        simulationtype="combined",
+        simulationtype=simulationtype,
         trackingdirection=direction,
-        weaksinkoption="pass_through",
-        weaksourceoption="pass_through",
+        weaksinkoption=weaksinkoption,
+        weaksourceoption=weaksourceoption,
         referencetime=ref_time,
         stoptimeoption=stoptimeoption,
         stoptime=stoptime,
diff --git a/nlmod/read/__init__.py b/nlmod/read/__init__.py
index 63dbf4f4..efbcfac4 100644
--- a/nlmod/read/__init__.py
+++ b/nlmod/read/__init__.py
@@ -13,6 +13,7 @@
     rws,
     waterboard,
     webservices,
+    nhi,
 )
 from .geotop import get_geotop
 from .regis import get_regis
diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index c0c0e571..9499817f 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -131,7 +131,9 @@ def _add_ts_to_ds(timeseries, loc_sel, variable, ds):
     """Add a timeseries to a variable at location loc_sel in model DataSet."""
     end = pd.Timestamp(ds.time.data[-1])
     if timeseries.index[-1] < end:
-        raise ValueError(f"no recharge available at {timeseries.name} for date {end}")
+        raise ValueError(
+            f"no data available for time series'{timeseries.name}' on date {end}"
+        )
 
     # fill recharge data array
     model_recharge = pd.Series(index=ds.time, dtype=float)
diff --git a/nlmod/read/nhi.py b/nlmod/read/nhi.py
new file mode 100644
index 00000000..858e3a16
--- /dev/null
+++ b/nlmod/read/nhi.py
@@ -0,0 +1,175 @@
+import logging
+import os
+
+import numpy as np
+import requests
+import rioxarray
+
+from ..dims.resample import structured_da_to_ds
+
+logger = logging.getLogger(__name__)
+
+
+def download_file(url, pathname, filename=None, overwrite=False, timeout=120.0):
+    """
+    Download a file from the NHI website.
+
+    Parameters
+    ----------
+    url : str
+        The url of the file to download.
+    pathname : str
+        The pathname to which the file is downloaded.
+    filename : str, optional
+        The name of the file to contain the downloadded data. When filename is None, it
+        is derived from url. The default is None.
+    overwrite : bool, optional
+        Overwrite the file if it allready exists. The default is False.
+    timeout : float, optional
+        How many seconds to wait for the server to send data before giving up. The
+        default is 120.
+
+    Returns
+    -------
+    fname : str
+        The full path of the downloaded file.
+
+    """
+    if filename is None:
+        filename = url.split("/")[-1]
+    fname = os.path.join(pathname, filename)
+    if overwrite or not os.path.isfile(fname):
+        logger.info(f"Downloading {filename}")
+        r = requests.get(url, allow_redirects=True, timeout=timeout)
+        with open(fname, "wb") as file:
+            file.write(r.content)
+    return fname
+
+
+def download_buisdrainage(pathname, overwrite=False):
+    """
+    Download resistance and depth of buisdrainage from the NHI website
+
+    Parameters
+    ----------
+    pathname : str
+        The pathname to which the files are downloaded.
+    overwrite : bool, optional
+        Overwrite the files if they allready exists. The default is False.
+
+    Returns
+    -------
+    fname_c : str
+        The full path of the downloaded file containing the resistance of buisdrainage.
+    fname_d : str
+        The full path of the downloaded file containing the depth of buisdrainage.
+
+    """
+    url_bas = "https://thredds.data.nhi.nu/thredds/fileServer/opendap/models/nhi3_2/25m"
+
+    # download resistance
+    url = f"{url_bas}/buisdrain_c_ras25/buisdrain_c_ras25.nc"
+    fname_c = download_file(url, pathname, overwrite=overwrite)
+
+    # download drain depth
+    url = f"{url_bas}/buisdrain_d_ras25/buisdrain_d_ras25.nc"
+    fname_d = download_file(url, pathname, overwrite=overwrite)
+
+    return fname_c, fname_d
+
+
+def add_buisdrainage(
+    ds,
+    pathname=None,
+    cond_var="buisdrain_cond",
+    depth_var="buisdrain_depth",
+    cond_method="average",
+    depth_method="mode",
+):
+    """
+    Add data about the buisdrainage to the model Dataset.
+
+    This data consists of the conductance of buisdrainage (m2/d) and the depth of
+    buisdrainage (m to surface level). With the default settings for `cond_method` and
+    `depth_method`, the conductance is calculated from the area weighted average of the
+    resistance in a cell, while the depth is set to the value which appears most often
+    in a cell. Cells without data (or 0 in the case of the resistance) are ignored in
+    these calculations.
+
+    The original data of the resistance and the depth of buisdrainage at a 25x25 m scale
+    is downloaded to `pathname` when not found.
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        The model Dataset.
+    pathname : str, optional
+        The pathname containing the downloaded files or the pathname to which the files
+        are downloaded. When pathname is None, it is set the the cachedir. The default
+        is None.
+    cond_var : str, optional
+        The name of the variable in ds to contain the data about the conductance of
+        buisdrainage. The default is "buisdrain_cond".
+    depth_var : str, optional
+        The name of the variable in ds to contain the data about the depth of
+        buisdrainage. The default is "buisdrain_depth".
+    cond_method : str, optional
+        The method to transform the conductance of buisdrainage to the model Dataset.
+        The default is "average".
+    depth_method : str, optional
+        The method to transform the depth of buisdrainage to the model Dataset. The
+        default is "mode".
+
+    Returns
+    -------
+    ds : xr.Dataset
+        The model dataset with added variables with the names `cond_var` and
+        `depth_var`.
+
+    """
+    if pathname is None:
+        pathname = ds.cachedir
+    # download files if needed
+    fname_c, fname_d = download_buisdrainage(pathname)
+
+    # make sure crs is set on ds
+    if ds.rio.crs is None:
+        ds = ds.rio.write_crs(28992)
+
+    # use cond_methd for conductance
+    # (default is "average" to account for locations without pipe drainage, where the
+    # conductance is 0)
+    buisdrain_c = rioxarray.open_rasterio(fname_c, mask_and_scale=True)[0]
+    # calculate a conductance (per m2) from a resistance
+    cond = 1 / buisdrain_c
+    # set conductance to 0 where resistance is infinite or 0
+    cond = cond.where(~(np.isinf(cond) | np.isnan(cond)), 0.0)
+    cond = cond.rio.write_crs(buisdrain_c.rio.crs)
+    # resample to model grid
+    ds[cond_var] = structured_da_to_ds(cond, ds, method=cond_method)
+    # multiply by area to get a conductance
+    ds[cond_var] = ds[cond_var] * ds["area"]
+
+    # use depth_method to retrieve the depth
+    # (default is "mode" for depth that occurs most in each cell)
+    mask_and_scale = False
+    buisdrain_d = rioxarray.open_rasterio(fname_d, mask_and_scale=mask_and_scale)[0]
+    if mask_and_scale:
+        nodata = np.nan
+    else:
+        nodata = buisdrain_d.attrs["_FillValue"]
+    # set buisdrain_d to nodata where it is 0
+    mask = buisdrain_d != 0
+    buisdrain_d = buisdrain_d.where(mask, nodata).rio.write_crs(buisdrain_d.rio.crs)
+    # resample to model grid
+    ds[depth_var] = structured_da_to_ds(
+        buisdrain_d, ds, method=depth_method, nodata=nodata
+    )
+    if not mask_and_scale:
+        # set nodata values to NaN
+        ds[depth_var] = ds[depth_var].where(ds[depth_var] != nodata)
+
+    # from cm to m
+    ds[depth_var] = ds[depth_var] / 100.0
+
+    return ds
diff --git a/pyproject.toml b/pyproject.toml
index e3de3d4a..d63b685c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta"
 [project]
 name = "nlmod"
 dynamic = ["version"]
-description = "nlmod is an open-source Python package for building Modflow 6 groundwater models from online data sources in The Netherlands"
+description = "Python package to build, run and visualize MODFLOW 6 groundwater models in the Netherlands."
 license = { file = "LICENSE" }
 readme = "README.md"
 authors = [
diff --git a/tests/test_005_external_data.py b/tests/test_005_external_data.py
index 8019b418..da72fe7a 100644
--- a/tests/test_005_external_data.py
+++ b/tests/test_005_external_data.py
@@ -1,3 +1,5 @@
+import pandas as pd
+import pytest
 import test_001_model
 
 import nlmod
@@ -11,9 +13,10 @@ def test_get_recharge():
     ds.update(nlmod.read.knmi.get_recharge(ds))
 
 
-def test_get_reacharge_most_common():
+def test_get_recharge_most_common():
     # model with sea
-    ds = test_001_model.get_ds_from_cache("basic_sea_model")
+    ds = nlmod.get_ds([100000, 110000, 420000, 430000])
+    ds = nlmod.time.set_ds_time(ds, start="2021", time=pd.date_range("2022", "2023"))
 
     # add knmi recharge to the model dataset
     ds.update(nlmod.read.knmi.get_recharge(ds, most_common_station=True))
@@ -31,6 +34,14 @@ def test_get_recharge_steady_state():
     ds.update(nlmod.read.knmi.get_recharge(ds))
 
 
+def test_get_recharge_not_available():
+    ds = nlmod.get_ds([100000, 110000, 420000, 430000])
+    time = pd.date_range("2010", pd.Timestamp.now())
+    ds = nlmod.time.set_ds_time(ds, start="2000", time=time)
+    with pytest.raises(ValueError):
+        ds.update(nlmod.read.knmi.get_recharge(ds))
+
+
 def test_ahn_within_extent():
     extent = [95000.0, 105000.0, 494000.0, 500000.0]
     da = nlmod.read.ahn.get_ahn_from_wcs(extent)
diff --git a/tests/test_021_nhi.py b/tests/test_021_nhi.py
new file mode 100644
index 00000000..af768339
--- /dev/null
+++ b/tests/test_021_nhi.py
@@ -0,0 +1,22 @@
+import os
+import numpy as np
+import tempfile
+import nlmod
+import pytest
+
+tmpdir = tempfile.gettempdir()
+
+
+@pytest.mark.slow
+def test_buidrainage():
+    model_ws = os.path.join(tmpdir, "buidrain")
+    ds = nlmod.get_ds([110_000, 130_000, 435_000, 445_000], model_ws=model_ws)
+    ds = nlmod.read.nhi.add_buisdrainage(ds)
+
+    # assert that all locations with a specified depth also have a positive conductance
+    mask = ~ds["buisdrain_depth"].isnull()
+    assert np.all(ds["buisdrain_cond"].data[mask] > 0)
+
+    # assert that all locations with a positive conductance also have a specified depth
+    mask = ds["buisdrain_cond"] > 0
+    assert np.all(~np.isnan(ds["buisdrain_depth"].data[mask]))