diff --git a/docs/examples/00_model_from_scratch.ipynb b/docs/examples/00_model_from_scratch.ipynb
index da381157..b043358b 100644
--- a/docs/examples/00_model_from_scratch.ipynb
+++ b/docs/examples/00_model_from_scratch.ipynb
@@ -122,7 +122,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds = nlmod.time.set_ds_time(ds, time=pd.Timestamp.today())"
+    "ds = nlmod.time.set_ds_time(ds, time=pd.Timestamp.today(), start=1)"
    ]
   },
   {
diff --git a/docs/examples/01_basic_model.ipynb b/docs/examples/01_basic_model.ipynb
index 7d1699a7..66b7e294 100644
--- a/docs/examples/01_basic_model.ipynb
+++ b/docs/examples/01_basic_model.ipynb
@@ -99,7 +99,7 @@
     "\n",
     "# add time discretisation\n",
     "ds = nlmod.time.set_ds_time(\n",
-    "    ds, start_time=start_time, steady_state=steady_state, perlen=365 * 5\n",
+    "    ds, start=start_time, steady=steady_state, perlen=365 * 5\n",
     ")\n",
     "\n",
     "if add_northsea:\n",
diff --git a/docs/examples/02_surface_water.ipynb b/docs/examples/02_surface_water.ipynb
index cdb38806..8e705bd7 100644
--- a/docs/examples/02_surface_water.ipynb
+++ b/docs/examples/02_surface_water.ipynb
@@ -351,7 +351,7 @@
     "ds = nlmod.to_model_ds(layer_model, model_name, model_ws, delr=delr, delc=delc)\n",
     "\n",
     "# create model time dataset\n",
-    "ds = nlmod.time.set_ds_time(ds, start_time=start_time, steady_state=True)\n",
+    "ds = nlmod.time.set_ds_time(ds, start=start_time, steady=True, perlen=1)\n",
     "\n",
     "ds"
    ]
diff --git a/docs/examples/03_local_grid_refinement.ipynb b/docs/examples/03_local_grid_refinement.ipynb
index 53fb8379..7e567600 100644
--- a/docs/examples/03_local_grid_refinement.ipynb
+++ b/docs/examples/03_local_grid_refinement.ipynb
@@ -98,11 +98,10 @@
     "# add time discretisation\n",
     "ds = nlmod.time.set_ds_time(\n",
     "    ds,\n",
-    "    start_time=start_time,\n",
-    "    steady_state=steady_state,\n",
+    "    start=start_time,\n",
+    "    steady=steady_state,\n",
     "    steady_start=steady_start,\n",
-    "    transient_timesteps=transient_timesteps,\n",
-    "    perlen=perlen,\n",
+    "    perlen=[perlen] * (transient_timesteps+1),\n",
     ")"
    ]
   },
diff --git a/docs/examples/08_gis.ipynb b/docs/examples/08_gis.ipynb
index 358dee67..cde523a3 100644
--- a/docs/examples/08_gis.ipynb
+++ b/docs/examples/08_gis.ipynb
@@ -71,7 +71,7 @@
    "outputs": [],
    "source": [
     "ds_struc = xr.load_dataset(\n",
-    "    os.path.join(model_ws, \"cache\", f\"{model_name}.nc\"), mask_and_scale=False\n",
+    "    os.path.join(model_ws, f\"{model_name}.nc\"), mask_and_scale=False\n",
     ")\n",
     "ds_struc"
    ]
@@ -111,7 +111,7 @@
    "outputs": [],
    "source": [
     "ds_vert = xr.load_dataset(\n",
-    "    os.path.join(model_ws, \"cache\", f\"{model_name}.nc\"), mask_and_scale=False\n",
+    "    os.path.join(model_ws, f\"{model_name}.nc\"), mask_and_scale=False\n",
     ")\n",
     "\n",
     "# get modelgrid\n",
diff --git a/docs/examples/09_schoonhoven.ipynb b/docs/examples/09_schoonhoven.ipynb
index b9d87e0f..530f0852 100644
--- a/docs/examples/09_schoonhoven.ipynb
+++ b/docs/examples/09_schoonhoven.ipynb
@@ -226,7 +226,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds = nlmod.time.set_ds_time(ds, time=time)\n"
+    "ds = nlmod.time.set_ds_time(ds, time=time, start=3652)\n"
    ]
   },
   {
@@ -338,6 +338,7 @@
     "    \"identificatie\"\n",
     "].isin(ids_oude_haven)\n",
     "lakes = bgt_grid[mask].copy()\n",
+	"lakes[\"name\"] = \"\"\n",
     "lakes.loc[lakes[\"identificatie\"].isin(ids_grote_gracht), \"name\"] = \"grote gracht\"\n",
     "lakes.loc[lakes[\"identificatie\"].isin(ids_oude_haven), \"name\"] = \"oude haven\"\n",
     "bgt_grid = bgt_grid[~mask]\n"
@@ -502,7 +503,7 @@
     "    colorbar_label=\"head [m NAP]\",\n",
     "    title=\"mean head\",\n",
     ")\n",
-    "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\");"
+    "bgt.plot(ax=pc.axes, edgecolor=\"k\", facecolor=\"none\");"
    ]
   },
   {
@@ -552,7 +553,7 @@
     "x = 118228\n",
     "y = 439870\n",
     "head_point = nlmod.gwf.get_head_at_point(head, x=x, y=y, ds=ds)\n",
-    "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n",
+    "fig, ax = plt.subplots(1, 1, figsize=(10, 8))\n",
     "handles = head_point.plot.line(ax=ax, hue=\"layer\")\n",
     "ax.set_ylabel(\"head [m NAP]\");"
    ]
diff --git a/docs/examples/10_modpath.ipynb b/docs/examples/10_modpath.ipynb
index d29b43ca..0d0b2b33 100644
--- a/docs/examples/10_modpath.ipynb
+++ b/docs/examples/10_modpath.ipynb
@@ -63,7 +63,7 @@
     "model_ws = \"ijmuiden\"\n",
     "model_name = \"IJm_planeten\"\n",
     "\n",
-    "ds = xr.open_dataset(os.path.join(model_ws, \"cache\", f\"{model_name}.nc\"))"
+    "ds = xr.open_dataset(os.path.join(model_ws, f\"{model_name}.nc\"))"
    ]
   },
   {
@@ -128,7 +128,7 @@
    "outputs": [],
    "source": [
     "# run modpath model\n",
-    "nlmod.modpath.write_and_run(mpf, nb_path=\"10_modpath.ipynb\")"
+    "nlmod.modpath.write_and_run(mpf, script_path=\"10_modpath.ipynb\")"
    ]
   },
   {
@@ -228,7 +228,7 @@
    "outputs": [],
    "source": [
     "# run modpath model\n",
-    "nlmod.modpath.write_and_run(mpf, nb_path=\"10_modpath.ipynb\")"
+    "nlmod.modpath.write_and_run(mpf, script_path=\"10_modpath.ipynb\")"
    ]
   },
   {
@@ -332,7 +332,7 @@
    "outputs": [],
    "source": [
     "# run modpath model\n",
-    "nlmod.modpath.write_and_run(mpf, nb_path=\"10_modpath.ipynb\")"
+    "nlmod.modpath.write_and_run(mpf, script_path=\"10_modpath.ipynb\")"
    ]
   },
   {
@@ -440,7 +440,7 @@
    "outputs": [],
    "source": [
     "# run modpath model\n",
-    "nlmod.modpath.write_and_run(mpf, nb_path=\"10_modpath.ipynb\")"
+    "nlmod.modpath.write_and_run(mpf, script_path=\"10_modpath.ipynb\")"
    ]
   },
   {
diff --git a/docs/examples/11_grid_rotation.ipynb b/docs/examples/11_grid_rotation.ipynb
index c58d4c28..852f2efb 100644
--- a/docs/examples/11_grid_rotation.ipynb
+++ b/docs/examples/11_grid_rotation.ipynb
@@ -83,7 +83,7 @@
     "# this period will determine the recharge-value and the ratio between summer and\n",
     "# winter stage of the surface water.\n",
     "\n",
-    "ds = nlmod.time.set_ds_time(ds, time=\"2023-1-1\")"
+    "ds = nlmod.time.set_ds_time(ds, time=\"2023-1-1\", start=\"2013\")"
    ]
   },
   {
@@ -358,8 +358,14 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "artesia",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
-   "name": "python"
+   "name": "python",
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,
diff --git a/docs/examples/13_plot_methods.ipynb b/docs/examples/13_plot_methods.ipynb
index 26d3faeb..ebc79dcc 100644
--- a/docs/examples/13_plot_methods.ipynb
+++ b/docs/examples/13_plot_methods.ipynb
@@ -73,8 +73,7 @@
    "source": [
     "model_name = \"Schoonhoven\"\n",
     "model_ws = \"schoonhoven\"\n",
-    "figdir, cachedir = nlmod.util.get_model_dirs(model_ws)\n",
-    "ds = xr.open_dataset(os.path.join(cachedir, f\"{model_name}.nc\"))\n",
+    "ds = xr.open_dataset(os.path.join(model_ws, f\"{model_name}.nc\"))\n",
     "\n",
     "# add calculated heads\n",
     "ds[\"head\"] = nlmod.gwf.get_heads_da(ds)\n",
diff --git a/docs/examples/14_stromingen_example.ipynb b/docs/examples/14_stromingen_example.ipynb
index b7098431..41f4ab96 100644
--- a/docs/examples/14_stromingen_example.ipynb
+++ b/docs/examples/14_stromingen_example.ipynb
@@ -164,7 +164,7 @@
    "source": [
     "# set model time settings\n",
     "t = date_range(tmin, tmax, freq=freq)\n",
-    "ds = nlmod.time.set_ds_time(ds, time=t, steady_start=True)\n",
+    "ds = nlmod.time.set_ds_time(ds, start=3652, time=t, steady_start=True)\n",
     "\n",
     "# build the modflow6 gwf model\n",
     "gwf = nlmod.gwf.ds_to_gwf(ds)"
diff --git a/docs/examples/16_groundwater_transport.ipynb b/docs/examples/16_groundwater_transport.ipynb
index e3828855..1278e55a 100644
--- a/docs/examples/16_groundwater_transport.ipynb
+++ b/docs/examples/16_groundwater_transport.ipynb
@@ -131,12 +131,10 @@
     "# add time discretisation\n",
     "ds = nlmod.time.set_ds_time(\n",
     "    ds,\n",
-    "    start_time=start_time,\n",
-    "    steady_state=False,\n",
+    "    start=start_time,\n",
+    "    steady=False,\n",
     "    steady_start=True,\n",
-    "    steady_start_perlen=1,\n",
-    "    transient_timesteps=10,\n",
-    "    perlen=365.0,\n",
+    "    perlen=[365.0]*10,\n",
     ")\n",
     "\n",
     "if ds.transport == 1:\n",
diff --git a/nlmod/dims/time.py b/nlmod/dims/time.py
index 583496a8..cb869b51 100644
--- a/nlmod/dims/time.py
+++ b/nlmod/dims/time.py
@@ -1,14 +1,16 @@
 import datetime as dt
 import logging
+import warnings
 
 import numpy as np
 import pandas as pd
+import xarray as xr
 from xarray import IndexVariable
 
 logger = logging.getLogger(__name__)
 
 
-def set_ds_time(
+def set_ds_time_deprecated(
     ds,
     time=None,
     steady_state=False,
@@ -71,6 +73,13 @@ def set_ds_time(
     ds : xarray.Dataset
         dataset with time variant model data
     """
+
+    warnings.warn(
+        "this function is deprecated and will eventually be removed, "
+        "please use nlmod.time.set_ds_time() in the future.",
+        DeprecationWarning,
+    )
+
     # checks
     if time_units.lower() != "days":
         raise NotImplementedError()
@@ -124,9 +133,192 @@ def set_ds_time(
     ds.time.attrs["steady_start"] = int(steady_start)
     ds.time.attrs["steady_state"] = int(steady_state)
 
+    # add to ds (for new version nlmod)
+    # add steady, nstp and tsmult to dataset
+    steady = int(steady_state) * np.ones(len(time_dt), dtype=int)
+    if steady_start:
+        steady[0] = 1
+    ds["steady"] = ("time",), steady
+
+    if isinstance(nstp, (int, np.integer)):
+        nstp = nstp * np.ones(len(time), dtype=int)
+    ds["nstp"] = ("time",), nstp
+
+    if isinstance(tsmult, float):
+        tsmult = tsmult * np.ones(len(time))
+    ds["tsmult"] = ("time",), tsmult
+
     return ds
 
 
+def set_ds_time(
+    ds,
+    start,
+    time=None,
+    steady=False,
+    steady_start=True,
+    time_units="DAYS",
+    perlen=None,
+    nstp=1,
+    tsmult=1.0,
+):
+    """Set time discretisation for model dataset.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        model dataset
+    start : int, float, str or pandas.Timestamp, optional
+        model start. When start is an integer or float it is interpreted as the number
+        of days of the first stress-period. When start is a string or pandas Timestamp
+        it is the start datetime of the simulation.
+    time : float, int or array-like, optional
+        float(s) (indicating elapsed time) or timestamp(s) corresponding to the end of
+        each stress period in the model. When time is a single value, the model will
+        have only one stress period. When time is None, the stress period lengths have
+        to be supplied via perlen. The default is None.
+    steady : arraylike or bool, optional
+        arraylike indicating which stress periods are steady-state, by default False,
+        which sets all stress periods to transient with the first period determined by
+        value of `steady_start`.
+    steady_start : bool, optional
+        whether to set the first period to steady-state, default is True, only used
+        when steady is passed as single boolean.
+    time_units : str, optional
+        time units, by default "DAYS"
+    perlen : float, int or array-like, optional
+        length of each stress-period. Only used when time is None. When perlen is a
+        single value, the model will have only one stress period. The default is None.
+    nstp : int or array-like, optional
+        number of steps per stress period, stored in ds.attrs, default is 1
+    tsmult : float, optional
+        timestep multiplier within stress periods, stored in ds.attrs, default is 1.0
+
+    Returns
+    -------
+    ds : xarray.Dataset
+        model dataset with added time coordinate
+
+    """
+    logger.info(
+        "Function set_ds_time() has changed since nlmod version 0.7."
+        " For the old behavior, use `nlmod.time.set_ds_time_deprecated()`."
+    )
+
+    if time is None and perlen is None:
+        raise (ValueError("Please specify either time or perlen in set_ds_time"))
+    elif perlen is not None:
+        if time is not None:
+            msg = f"Cannot use both time and perlen. Ignoring perlen: {perlen}"
+            logger.warning(msg)
+        else:
+            if isinstance(perlen, (int, np.integer, float)):
+                perlen = [perlen]
+            time = np.cumsum(perlen)
+
+    if isinstance(time, str) or not hasattr(time, "__iter__"):
+        time = [time]
+
+    # parse start
+    if isinstance(start, (int, np.integer, float)):
+        if isinstance(time[0], (int, np.integer, float)):
+            raise (ValueError("Make sure start or time contains a valid TimeStamp"))
+        start = time[0] - pd.to_timedelta(start, "D")
+    elif isinstance(start, str):
+        start = pd.Timestamp(start)
+    elif isinstance(start, (pd.Timestamp, np.datetime64)):
+        pass
+    else:
+        raise TypeError("Cannot parse start datetime.")
+
+    # convert time to Timestamps
+    if isinstance(time[0], (int, np.integer, float)):
+        time = pd.Timestamp(start) + pd.to_timedelta(time, time_units)
+    elif isinstance(time[0], str):
+        time = pd.to_datetime(time)
+    elif isinstance(time[0], (pd.Timestamp, np.datetime64, xr.core.variable.Variable)):
+        pass
+    else:
+        raise TypeError("Cannot process 'time' argument. Datatype not understood.")
+
+    if time[0] <= start:
+        msg = (
+            "The timestamp of the first stress period cannot be before or "
+            "equal to the model start time! Please modify `time` or `start`!"
+        )
+        logger.error(msg)
+        raise ValueError(msg)
+
+    ds = ds.assign_coords(coords={"time": time})
+
+    # add steady, nstp and tsmult to dataset
+    if isinstance(steady, bool):
+        steady = int(steady) * np.ones(len(time), dtype=int)
+        if steady_start:
+            steady[0] = 1
+    ds["steady"] = ("time",), steady
+
+    if isinstance(nstp, (int, np.integer)):
+        nstp = nstp * np.ones(len(time), dtype=int)
+    ds["nstp"] = ("time",), nstp
+
+    if isinstance(tsmult, float):
+        tsmult = tsmult * np.ones(len(time))
+    ds["tsmult"] = ("time",), tsmult
+
+    if time_units == "D":
+        time_units = "DAYS"
+    ds.time.attrs["time_units"] = time_units
+    ds.time.attrs["start"] = str(start)
+
+    return ds
+
+
+def ds_time_idx_from_tdis_settings(start, perlen, nstp=1, tsmult=1.0, time_units="D"):
+    """Get time index from TDIS perioddata: perlen, nstp, tsmult.
+
+
+    Parameters
+    ----------
+    start : str, pd.Timestamp
+        start datetime
+    perlen : array-like
+        array of period lengths
+    nstp : int, or array-like optional
+        number of steps per period, by default 1
+    tsmult : float or array-like, optional
+        timestep multiplier per period, by default 1.0
+    time_units : str, optional
+        time units, by default "D"
+
+    Returns
+    -------
+    IndexVariable
+        time coordinate for xarray data-array or dataset
+    """
+    deltlist = []
+    for kper, delt in enumerate(perlen):
+        if not isinstance(nstp, int):
+            kstpkper = nstp[kper]
+        else:
+            kstpkper = nstp
+
+        if not isinstance(tsmult, float):
+            tsm = tsmult[kper]
+        else:
+            tsm = tsmult
+
+        if tsm > 1.0:
+            delt0 = delt * (tsm - 1) / (tsm**kstpkper - 1)
+            delt = delt0 * tsm ** np.arange(kstpkper)
+        else:
+            delt = np.ones(kstpkper) * delt / kstpkper
+        deltlist.append(delt)
+
+    dt_arr = np.cumsum(np.concatenate(deltlist))
+    return ds_time_idx(dt_arr, start_datetime=start, time_units=time_units)
+
+
 def estimate_nstp(
     forcing, perlen=1, tsmult=1.1, nstp_min=1, nstp_max=25, return_dt_arr=False
 ):
@@ -214,6 +406,15 @@ def estimate_nstp(
 
 
 def ds_time_from_model(gwf):
+    warnings.warn(
+        "this function was renamed to `ds_time_idx_from_model`. "
+        "Please use the new function name.",
+        DeprecationWarning,
+    )
+    return ds_time_idx_from_model(gwf)
+
+
+def ds_time_idx_from_model(gwf):
     """Get time index variable from model (gwf or gwt).
 
     Parameters
@@ -227,10 +428,19 @@ def ds_time_from_model(gwf):
         time coordinate for xarray data-array or dataset
     """
 
-    return ds_time_from_modeltime(gwf.modeltime)
+    return ds_time_idx_from_modeltime(gwf.modeltime)
 
 
 def ds_time_from_modeltime(modeltime):
+    warnings.warn(
+        "this function was renamed to `ds_time_idx_from_model`. "
+        "Please use the new function name.",
+        DeprecationWarning,
+    )
+    return ds_time_idx_from_modeltime(modeltime)
+
+
+def ds_time_idx_from_modeltime(modeltime):
     """Get time index variable from modeltime object.
 
     Parameters
diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py
index 4e5962c2..2017abe5 100644
--- a/nlmod/gwf/gwf.py
+++ b/nlmod/gwf/gwf.py
@@ -578,18 +578,15 @@ def sto(
     """
     logger.info("creating mf6 STO")
 
-    if ds.time.steady_state:
+    if "time" not in ds or ds["steady"].all():
+        logger.warning("Model is steady-state, no STO package created.")
         return None
     else:
-        if ds.time.steady_start:
-            sts_spd = {0: True}
-            trn_spd = {1: True}
-        else:
-            sts_spd = None
-            trn_spd = {0: True}
+        sts_spd = {iper: bool(b) for iper, b in enumerate(ds["steady"])}
+        trn_spd = {iper: not bool(b) for iper, b in enumerate(ds["steady"])}
 
         sy = _get_value_from_ds_datavar(ds, "sy", sy, default=0.2)
-        ss = _get_value_from_ds_datavar(ds, "ss", ss, default=0.000001)
+        ss = _get_value_from_ds_datavar(ds, "ss", ss, default=1e-5)
 
         sto = flopy.mf6.ModflowGwfsto(
             gwf,
diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py
index 51242dc1..751751dc 100644
--- a/nlmod/gwf/recharge.py
+++ b/nlmod/gwf/recharge.py
@@ -35,7 +35,7 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", **kwargs):
         raise ValueError("please remove nan values in recharge data array")
 
     # get stress period data
-    if ds.time.steady_state:
+    if ds["steady"].all():
         recharge = "recharge"
         if "time" in ds["recharge"].dims:
             mask = ds["recharge"].isel(time=0) != 0
@@ -69,7 +69,7 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", **kwargs):
         **kwargs,
     )
 
-    if ds.time.steady_state:
+    if ds["steady"].all():
         return rch
 
     # create timeseries packages
@@ -128,7 +128,7 @@ def ds_to_evt(gwf, ds, pname="evt", nseg=1, surface=None, depth=None, **kwargs):
         raise ValueError("please remove nan values in evaporation data array")
 
     # get stress period data
-    if ds.time.steady_state:
+    if ds["steady"].all():
         if "time" in ds["evaporation"].dims:
             mask = ds["evaporation"].isel(time=0) != 0
         else:
@@ -163,7 +163,7 @@ def ds_to_evt(gwf, ds, pname="evt", nseg=1, surface=None, depth=None, **kwargs):
         **kwargs,
     )
 
-    if ds.time.steady_state:
+    if ds["steady"].all():
         return evt
 
     # create timeseries packages
diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py
index 4180714a..4198fc6a 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -509,7 +509,7 @@ def add_info_to_gdf(
         inds = s.query(geom_to)
         if len(inds) == 0:
             continue
-        overlap = gdf_from.geometry[inds].intersection(geom_to)
+        overlap = gdf_from.geometry.iloc[inds].intersection(geom_to)
         if geom_type is None:
             geom_type = overlap.geom_type.iloc[0]
         if geom_type in ["Polygon", "MultiPolygon"]:
diff --git a/nlmod/modpath/modpath.py b/nlmod/modpath/modpath.py
index f605071d..5850dba9 100644
--- a/nlmod/modpath/modpath.py
+++ b/nlmod/modpath/modpath.py
@@ -47,7 +47,7 @@ def write_and_run(mpf, remove_prev_output=True, script_path=None, silent=False):
     logger.info("write modpath files to model workspace")
 
     # write modpath datasets
-    mpf.write_input(silent=silent)
+    mpf.write_input()
 
     # run modpath
     logger.info("run modpath model")
@@ -196,12 +196,6 @@ def mpf(gwf, exe_name=None, modelname=None, model_ws=None):
             "the save_flows option of the npf package should be True not None"
         )
 
-    # check if the tdis has a start_time
-    if gwf.simulation.tdis.start_date_time.array is not None:
-        logger.warning(
-            "older versions of modpath cannot handle this, see https://github.com/MODFLOW-USGS/modpath-v7/issues/31"
-        )
-
     # get executable
     if exe_name is None:
         exe_name = util.get_exe_path("mp7_2_002_provisional")
diff --git a/nlmod/plot/plotutil.py b/nlmod/plot/plotutil.py
index d2fadae0..f41941d3 100644
--- a/nlmod/plot/plotutil.py
+++ b/nlmod/plot/plotutil.py
@@ -6,6 +6,7 @@
 from ..dims.resample import get_affine_mod_to_world
 from ..epsg28992 import EPSG_28992
 
+
 def get_patches(ds, rotated=False):
     """Get the matplotlib patches for a vertex grid."""
     assert "icell2d" in ds.dims
diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index 14a1d817..f104e819 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -85,7 +85,9 @@ def get_recharge(ds, method="linear", most_common_station=False):
         unique_combinations = locations.drop_duplicates(["stn_rd", "stn_ev24"])[
             ["stn_rd", "stn_ev24"]
         ].values
-
+        if unique_combinations.shape[1] > 2:
+            # bug fix for pandas 2.1 where three columns are returned
+            unique_combinations = unique_combinations[:, :2]
         for stn_rd, stn_ev24 in unique_combinations:
             # get locations with the same prec and evap station
             mask = (locations["stn_rd"] == stn_rd) & (locations["stn_ev24"] == stn_ev24)
diff --git a/nlmod/sim/sim.py b/nlmod/sim/sim.py
index 6e98a55e..911d949d 100644
--- a/nlmod/sim/sim.py
+++ b/nlmod/sim/sim.py
@@ -64,7 +64,7 @@ def write_and_run(sim, ds, write_ds=True, script_path=None, silent=False):
     ds.attrs["model_ran_on"] = dt.datetime.now().strftime("%Y%m%d_%H:%M:%S")
 
 
-def get_tdis_perioddata(ds):
+def get_tdis_perioddata(ds, nstp="nstp", tsmult="tsmult"):
     """Get tdis_perioddata from ds.
 
     Parameters
@@ -92,15 +92,15 @@ def get_tdis_perioddata(ds):
     if len(ds["time"]) > 1:
         perlen.extend(np.diff(ds["time"]) / deltat)
 
-    if "nstp" in ds:
-        nstp = ds["nstp"].values
-    else:
-        nstp = [ds.time.nstp] * len(perlen)
+    nstp = util._get_value_from_ds_datavar(ds, "nstp", nstp, return_da=False)
 
-    if "tsmult" in ds:
-        tsmult = ds["tsmult"].values
-    else:
-        tsmult = [ds.time.tsmult] * len(perlen)
+    if isinstance(nstp, (int, np.integer)):
+        nstp = [nstp] * len(perlen)
+
+    tsmult = util._get_value_from_ds_datavar(ds, "tsmult", tsmult, return_da=False)
+
+    if isinstance(tsmult, float):
+        tsmult = [tsmult] * len(perlen)
 
     tdis_perioddata = list(zip(perlen, nstp, tsmult))
 
@@ -144,7 +144,7 @@ def sim(ds, exe_name=None):
     return sim
 
 
-def tdis(ds, sim, pname="tdis"):
+def tdis(ds, sim, pname="tdis", nstp="nstp", tsmult="tsmult", **kwargs):
     """create tdis package from the model dataset.
 
     Parameters
@@ -156,6 +156,8 @@ def tdis(ds, sim, pname="tdis"):
         simulation object.
     pname : str, optional
         package name
+    **kwargs
+        passed on to flopy.mft.ModflowTdis
 
     Returns
     -------
@@ -166,7 +168,7 @@ def tdis(ds, sim, pname="tdis"):
     # start creating model
     logger.info("creating mf6 TDIS")
 
-    tdis_perioddata = get_tdis_perioddata(ds)
+    tdis_perioddata = get_tdis_perioddata(ds, nstp=nstp, tsmult=tsmult)
 
     # Create the Flopy temporal discretization object
     tdis = flopy.mf6.ModflowTdis(
@@ -176,6 +178,7 @@ def tdis(ds, sim, pname="tdis"):
         nper=len(ds.time),
         start_date_time=pd.Timestamp(ds.time.start).isoformat(),
         perioddata=tdis_perioddata,
+        **kwargs,
     )
 
     return tdis
diff --git a/pyproject.toml b/pyproject.toml
index 891b6b45..8464577d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -55,13 +55,20 @@ repository = "https://github.com/ArtesiaWater/nlmod"
 documentation = "https://nlmod.readthedocs.io/en/latest/"
 
 [project.optional-dependencies]
-full = ["nlmod[knmi]", "gdown", "geocube", "bottleneck", "contextily", "scikit-image"]
+full = [
+    "nlmod[knmi]",
+    "gdown",
+    "geocube",
+    "bottleneck",
+    "contextily",
+    "scikit-image",
+]
 knmi = ["h5netcdf", "nlmod[grib]"]
 grib = ["cfgrib", "ecmwflibs"]
 test = ["pytest>=7", "pytest-cov", "pytest-dependency"]
 nbtest = ["nbformat", "nbconvert>6.4.5"]
 lint = ["flake8", "isort", "black[jupyter]"]
-ci = ["nlmod[full,lint,test,nbtest]", "netCDF4>=1.6.3"]
+ci = ["nlmod[full,lint,test,nbtest]", "netCDF4>=1.6.3", "pandas<2.1.0"]
 rtd = [
     "nlmod[full]",
     "ipython",
diff --git a/tests/test_001_model.py b/tests/test_001_model.py
index 9c9c1362..f8fcf25b 100644
--- a/tests/test_001_model.py
+++ b/tests/test_001_model.py
@@ -1,6 +1,8 @@
 import os
 import tempfile
 
+import numpy as np
+import pandas as pd
 import pytest
 import xarray as xr
 
@@ -18,20 +20,17 @@ def test_model_directories(tmpdir):
 def get_ds_time_steady(tmpdir, modelname="test"):
     model_ws = os.path.join(tmpdir, "test_model")
     ds = nlmod.base.set_ds_attrs(xr.Dataset(), modelname, model_ws)
-    ds = nlmod.time.set_ds_time(ds, start_time="2015-1-1", steady_state=True)
+    ds = nlmod.time.set_ds_time(ds, time=["2015-1-2"], start="2015-1-1", steady=True)
     return ds
 
 
 def get_ds_time_transient(tmpdir, modelname="test"):
     model_ws = os.path.join(tmpdir, "test_model")
     ds = nlmod.base.set_ds_attrs(xr.Dataset(), modelname, model_ws)
-    ds = nlmod.time.set_ds_time(
-        ds,
-        start_time="2015-1-1",
-        steady_state=False,
-        steady_start=True,
-        transient_timesteps=10,
-    )
+    nper = 11
+    time = pd.date_range(start="2015-1-2", periods=nper, freq="D")
+    steady = np.zeros(nper)
+    ds = nlmod.time.set_ds_time(ds, time=time, start="2015-1-1", steady=steady)
     return ds
 
 
@@ -79,12 +78,14 @@ def test_create_small_model_grid_only(tmpdir, model_name="test"):
     )
     assert ds.dims["layer"] == 5
 
+    nper = 11
+    steady = np.zeros(nper, dtype=int)
+    steady[0] = 1
     ds = nlmod.time.set_ds_time(
         ds,
-        start_time="2015-1-1",
-        steady_state=False,
-        steady_start=True,
-        transient_timesteps=10,
+        time=pd.date_range("2015-1-2", periods=nper, freq="D"),
+        start="2015-1-1",
+        steady=steady,
     )
 
     # create simulation
@@ -118,12 +119,14 @@ def test_create_sea_model_grid_only(tmpdir, model_name="test"):
         regis_geotop_ds, model_name, model_ws, delr=100.0, delc=100.0
     )
 
+    nper = 11
+    steady = np.zeros(nper, dtype=int)
+    steady[0] = 1
     ds = nlmod.time.set_ds_time(
         ds,
-        start_time="2015-1-1",
-        steady_state=False,
-        steady_start=True,
-        transient_timesteps=10,
+        time=pd.date_range("2015-1-2", periods=nper, freq="D"),
+        start="2005-1-1",
+        steady=steady,
     )
 
     # save ds
@@ -206,9 +209,10 @@ def test_create_sea_model_perlen_list(tmpdir):
     ds = nlmod.base.set_ds_attrs(ds, ds.model_name, model_ws)
 
     # create transient with perlen list
+    start = ds.time.start
     perlen = [3650, 14, 10, 11]  # length of the time steps
-    transient_timesteps = 3
-    start_time = ds.time.start
+    steady = np.zeros(len(perlen), dtype=int)
+    steady[0] = 1
 
     # drop time dimension before setting time
     ds = ds.drop_dims("time")
@@ -216,11 +220,9 @@ def test_create_sea_model_perlen_list(tmpdir):
     # update current ds with new time dicretisation
     ds = nlmod.time.set_ds_time(
         ds,
-        start_time=start_time,
-        steady_state=False,
-        steady_start=True,
-        perlen=perlen,
-        transient_timesteps=transient_timesteps,
+        time=np.cumsum(perlen),
+        start=start,
+        steady=steady,
     )
 
     # create simulation
@@ -277,21 +279,24 @@ def test_create_sea_model_perlen_14(tmpdir):
     ds = nlmod.base.set_ds_attrs(ds, ds.model_name, model_ws)
 
     # create transient with perlen list
-    perlen = 14  # length of the time steps
-    transient_timesteps = 3
-    start_time = ds.time.start
+    perlen = 14  # length of the transient time steps
+    nper = 4
+    start = ds.time.start
+    perlen = perlen * np.ones(nper)
+    perlen[0] = 3652.0  # length of the steady state step
+    steady = np.zeros(nper, dtype=int)
+    steady[0] = 1
+    time = nlmod.time.ds_time_idx_from_tdis_settings(start, perlen=perlen)
 
     # drop time dimension before setting time
     ds = ds.drop_dims("time")
 
-    # update current ds with new time dicretisation
+    # update current ds with new time discretization
     ds = nlmod.time.set_ds_time(
         ds,
-        start_time=start_time,
-        steady_state=False,
-        steady_start=True,
-        perlen=perlen,
-        transient_timesteps=transient_timesteps,
+        time=time,
+        start=start,
+        steady=steady,
     )
 
     # create simulation
diff --git a/tests/test_005_external_data.py b/tests/test_005_external_data.py
index 42d32323..8019b418 100644
--- a/tests/test_005_external_data.py
+++ b/tests/test_005_external_data.py
@@ -25,7 +25,7 @@ def test_get_recharge_steady_state():
 
     # modify mtime
     ds = ds.drop_dims("time")
-    ds = nlmod.time.set_ds_time(ds, start_time="2000-1-1", perlen=3650)
+    ds = nlmod.time.set_ds_time(ds, time=[3650], start="2000-1-1")
 
     # add knmi recharge to the model dataset
     ds.update(nlmod.read.knmi.get_recharge(ds))
diff --git a/tests/test_008_waterschappen.py b/tests/test_008_waterschappen.py
index 1c5d18d1..bd9a6a0b 100644
--- a/tests/test_008_waterschappen.py
+++ b/tests/test_008_waterschappen.py
@@ -1,5 +1,6 @@
-import pytest
 import matplotlib
+import pytest
+
 import nlmod
 
 
diff --git a/tests/test_010_wells.py b/tests/test_010_wells.py
index 70c382d2..3e704db4 100644
--- a/tests/test_010_wells.py
+++ b/tests/test_010_wells.py
@@ -18,7 +18,7 @@ def get_model_ds():
         model_name="from_scratch",
     )
 
-    ds = nlmod.time.set_ds_time(ds, time=pd.Timestamp.today())
+    ds = nlmod.time.set_ds_time(ds, time=[1], start=pd.Timestamp.today())
 
     return ds
 
diff --git a/tests/test_013_surface_water.py b/tests/test_013_surface_water.py
index 8a10b8ce..d05ed8d0 100644
--- a/tests/test_013_surface_water.py
+++ b/tests/test_013_surface_water.py
@@ -12,7 +12,7 @@ def test_gdf_to_seasonal_pkg():
     ds = nlmod.get_ds(
         [170000, 171000, 550000, 551000], model_ws=model_ws, model_name=model_name
     )
-    ds = nlmod.time.set_ds_time(ds, time=pd.Timestamp.today())
+    ds = nlmod.time.set_ds_time(ds, time=[365.0], start=pd.Timestamp.today())
     gdf = nlmod.gwf.surface_water.get_gdf(ds)
 
     sim = nlmod.sim.sim(ds)
@@ -33,7 +33,7 @@ def test_gdf_lake():
     ds = nlmod.get_ds(
         [170000, 171000, 550000, 551000], model_ws=model_ws, model_name=model_name
     )
-    ds = nlmod.time.set_ds_time(ds, time=pd.Timestamp.today())
+    ds = nlmod.time.set_ds_time(ds, time=[1], start=pd.Timestamp.today())
     ds = nlmod.dims.refine(ds)
 
     sim = nlmod.sim.sim(ds)
diff --git a/tests/test_015_gwf_output.py b/tests/test_015_gwf_output.py
index 7f4a0f73..2173948a 100644
--- a/tests/test_015_gwf_output.py
+++ b/tests/test_015_gwf_output.py
@@ -31,13 +31,7 @@ def test_create_small_model_grid_only(tmpdir, model_name="test"):
     )
     assert ds.dims["layer"] == 5
 
-    ds = nlmod.time.set_ds_time(
-        ds,
-        start_time="2015-1-1",
-        steady_state=False,
-        steady_start=True,
-        transient_timesteps=2,
-    )
+    ds = nlmod.time.set_ds_time(ds, time=[1, 2, 3], start="2015-1-1", steady=[1, 0, 0])
 
     # create simulation
     sim = nlmod.sim.sim(ds)
diff --git a/tests/test_016_time.py b/tests/test_016_time.py
index a606e727..1e3339c5 100644
--- a/tests/test_016_time.py
+++ b/tests/test_016_time.py
@@ -1,3 +1,5 @@
+import numpy as np
+
 import nlmod
 
 
@@ -17,3 +19,12 @@ def test_estimate_nstp():
     assert nstp[-1] == nstp_min
     assert max(nstp) == nstp_max
     assert min(nstp) == nstp_min
+
+
+def test_ds_time_from_tdis_settings():
+    tidx = nlmod.time.ds_time_idx_from_tdis_settings(
+        "2000", [100, 100, 100], nstp=[1, 2, 2], tsmult=[1.0, 1.0, 2.0]
+    )
+
+    elapsed = (tidx.to_numpy() - np.datetime64("2000")) / np.timedelta64(1, "D")
+    assert np.allclose(elapsed, [100, 150, 200, 233.33333333, 300.0])