diff --git a/docs/examples/01_basic_model.ipynb b/docs/examples/01_basic_model.ipynb
index a1ee2ef6..fa46f0d7 100644
--- a/docs/examples/01_basic_model.ipynb
+++ b/docs/examples/01_basic_model.ipynb
@@ -24,7 +24,7 @@
     "import flopy\n",
     "import geopandas as gpd\n",
     "import matplotlib.pyplot as plt\n",
-    "import nlmod"
+    "import nlmod\n"
    ]
   },
   {
@@ -35,7 +35,7 @@
    "source": [
     "print(f\"nlmod version: {nlmod.__version__}\")\n",
     "\n",
-    "nlmod.util.get_color_logger(\"INFO\");"
+    "nlmod.util.get_color_logger(\"INFO\")"
    ]
   },
   {
@@ -76,7 +76,7 @@
     "steady_state = True\n",
     "start_time = \"2015-1-1\"\n",
     "add_northsea = True\n",
-    "starting_head = 1.0"
+    "starting_head = 1.0\n"
    ]
   },
   {
@@ -103,7 +103,7 @@
     ")\n",
     "\n",
     "if add_northsea:\n",
-    "    ds = nlmod.read.rws.add_northsea(ds, cachedir=cachedir)"
+    "    ds = nlmod.read.rws.add_northsea(ds, cachedir=cachedir)\n"
    ]
   },
   {
@@ -134,7 +134,7 @@
     "ic = nlmod.gwf.ic(ds, gwf, starting_head=starting_head)\n",
     "\n",
     "# Create the output control package\n",
-    "oc = nlmod.gwf.oc(ds, gwf)"
+    "oc = nlmod.gwf.oc(ds, gwf)\n"
    ]
   },
   {
@@ -152,7 +152,7 @@
     "ds.update(rws_ds)\n",
     "\n",
     "# build ghb package\n",
-    "ghb = nlmod.gwf.ghb(ds, gwf, bhead=f\"{da_name}_stage\", cond=f\"{da_name}_cond\")"
+    "ghb = nlmod.gwf.ghb(ds, gwf, bhead=f\"{da_name}_stage\", cond=f\"{da_name}_cond\")\n"
    ]
   },
   {
@@ -167,7 +167,7 @@
     "ds.update(ahn_ds)\n",
     "\n",
     "# build surface level drain package\n",
-    "drn = nlmod.gwf.surface_drain_from_ds(ds, gwf, resistance=10.0)"
+    "drn = nlmod.gwf.surface_drain_from_ds(ds, gwf, resistance=10.0)\n"
    ]
   },
   {
@@ -178,7 +178,7 @@
    "source": [
     "# add constant head cells at model boundaries\n",
     "ds.update(nlmod.grid.mask_model_edge(ds, ds[\"idomain\"]))\n",
-    "chd = nlmod.gwf.chd(ds, gwf, mask=\"edge_mask\", head=\"starting_head\")"
+    "chd = nlmod.gwf.chd(ds, gwf, mask=\"edge_mask\", head=\"starting_head\")\n"
    ]
   },
   {
@@ -193,7 +193,7 @@
     "ds.update(knmi_ds)\n",
     "\n",
     "# create recharge package\n",
-    "rch = nlmod.gwf.rch(ds, gwf)"
+    "rch = nlmod.gwf.rch(ds, gwf)\n"
    ]
   },
   {
@@ -209,7 +209,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds"
+    "ds\n"
    ]
   },
   {
@@ -228,7 +228,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "nlmod.sim.write_and_run(sim, ds, write_ds=True, nb_path=\"01_basic_model.ipynb\")"
+    "nlmod.sim.write_and_run(sim, ds, write_ds=True, script_path=\"01_basic_model.ipynb\")\n"
    ]
   },
   {
@@ -246,7 +246,7 @@
    "outputs": [],
    "source": [
     "ax = nlmod.plot.modelgrid(ds)\n",
-    "nlmod.plot.surface_water(ds, ax=ax)"
+    "nlmod.plot.surface_water(ds, ax=ax)\n"
    ]
   },
   {
@@ -279,7 +279,7 @@
     "\n",
     "for axes1 in axes:\n",
     "    for ax in axes1:\n",
-    "        ax.axis(\"scaled\")"
+    "        ax.axis(\"scaled\")\n"
    ]
   }
  ],
diff --git a/docs/examples/02_surface_water.ipynb b/docs/examples/02_surface_water.ipynb
index f96a20bd..4c8e3c00 100644
--- a/docs/examples/02_surface_water.ipynb
+++ b/docs/examples/02_surface_water.ipynb
@@ -42,7 +42,7 @@
    "source": [
     "print(f\"nlmod version: {nlmod.__version__}\")\n",
     "\n",
-    "nlmod.util.get_color_logger(\"INFO\");"
+    "nlmod.util.get_color_logger(\"INFO\")\n"
    ]
   },
   {
@@ -152,7 +152,7 @@
    "outputs": [],
    "source": [
     "f, ax = nlmod.plot.get_map(extent)\n",
-    "bgt.plot(\"bronhouder\", legend=True, ax=ax);"
+    "bgt.plot(\"bronhouder\", legend=True, ax=ax)\n"
    ]
   },
   {
@@ -453,7 +453,7 @@
     "xmin, xmax, ymin, ymax = extent\n",
     "offset = 100\n",
     "ax.set_xlim(xmin - offset, xmax + offset)\n",
-    "ax.set_ylim(ymin - offset, ymax + offset);"
+    "ax.set_ylim(ymin - offset, ymax + offset)\n"
    ]
   },
   {
@@ -501,7 +501,7 @@
     "gwf.modelgrid.plot(ax=ax)\n",
     "ax.set_xlim(xlim[0], xlim[0] + ds.delr * 1.1)\n",
     "ax.set_ylim(ylim)\n",
-    "ax.set_title(f\"Surface water shapes in cell: {cid}\");"
+    "ax.set_title(f\"Surface water shapes in cell: {cid}\")\n"
    ]
   },
   {
@@ -713,7 +713,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "nlmod.sim.write_and_run(sim, ds, write_ds=True, nb_path=\"02_surface_water.ipynb\")"
+    "nlmod.sim.write_and_run(sim, ds, write_ds=True, script_path=\"02_surface_water.ipynb\")"
    ]
   },
   {
@@ -757,7 +757,7 @@
     "qm = mv.plot_array(head[-1], cmap=\"RdBu\")  # last timestep\n",
     "mv.plot_ibound()  # plot inactive cells in red\n",
     "fig.colorbar(qm, shrink=1.0)\n",
-    "ax.set_title(f\"Heads top-view, layer {ilay}\");"
+    "ax.set_title(f\"Heads top-view, layer {ilay}\")\n"
    ]
   },
   {
@@ -783,7 +783,7 @@
     "xs.plot_ibound()  # plot inactive cells in red\n",
     "fig.colorbar(qm, shrink=1.0)\n",
     "col = gwf.modelgrid.ncol // 2\n",
-    "ax.set_title(f\"Cross-section along column {col}\");"
+    "ax.set_title(f\"Cross-section along column {col}\")\n"
    ]
   }
  ],
diff --git a/docs/examples/03_local_grid_refinement.ipynb b/docs/examples/03_local_grid_refinement.ipynb
index 79d7f9cd..d6461114 100644
--- a/docs/examples/03_local_grid_refinement.ipynb
+++ b/docs/examples/03_local_grid_refinement.ipynb
@@ -40,7 +40,7 @@
    "source": [
     "print(f\"nlmod version: {nlmod.__version__}\")\n",
     "\n",
-    "nlmod.util.get_color_logger(\"INFO\");"
+    "nlmod.util.get_color_logger(\"INFO\")\n"
    ]
   },
   {
@@ -233,7 +233,7 @@
    "outputs": [],
    "source": [
     "nlmod.sim.write_and_run(\n",
-    "    sim, ds, write_ds=True, nb_path=\"03_local_grid_refinement.ipynb\"\n",
+    "    sim, ds, write_ds=True, script_path=\"03_local_grid_refinement.ipynb\"\n",
     ")"
    ]
   },
diff --git a/docs/examples/09_schoonhoven.ipynb b/docs/examples/09_schoonhoven.ipynb
index e07f069c..8e305a80 100644
--- a/docs/examples/09_schoonhoven.ipynb
+++ b/docs/examples/09_schoonhoven.ipynb
@@ -29,7 +29,7 @@
     "import geopandas as gpd\n",
     "from nlmod.plot import DatasetCrossSection\n",
     "from shapely.geometry import LineString, Point\n",
-    "import warnings"
+    "import warnings\n"
    ]
   },
   {
@@ -41,7 +41,7 @@
    "source": [
     "print(f\"nlmod version: {nlmod.__version__}\")\n",
     "\n",
-    "nlmod.util.get_color_logger(\"INFO\");"
+    "nlmod.util.get_color_logger(\"INFO\")"
    ]
   },
   {
@@ -64,7 +64,7 @@
     "model_ws = \"schoonhoven\"\n",
     "figdir, cachedir = nlmod.util.get_model_dirs(model_ws)\n",
     "extent = [116_500, 120_000, 439_000, 442_000]\n",
-    "time = pd.date_range(\"2020\", \"2023\", freq=\"MS\")  # monthly timestep"
+    "time = pd.date_range(\"2020\", \"2023\", freq=\"MS\")  # monthly timestep\n"
    ]
   },
   {
@@ -98,7 +98,7 @@
     "            f\"{fname_bgt} not found. Please run notebook 02_surface_water.ipynb first\"\n",
     "        )\n",
     "    )\n",
-    "bgt = gpd.read_file(fname_bgt)"
+    "bgt = gpd.read_file(fname_bgt)\n"
    ]
   },
   {
@@ -121,7 +121,7 @@
     "norm = matplotlib.colors.Normalize(vmin=-3, vmax=1)\n",
     "cmap = \"viridis\"\n",
     "bgt.plot(\"summer_stage\", ax=ax, norm=norm, cmap=cmap)\n",
-    "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)"
+    "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)\n"
    ]
   },
   {
@@ -141,7 +141,7 @@
    "source": [
     "f, ax = nlmod.plot.get_map(extent)\n",
     "bgt.plot(\"ahn_min\", ax=ax, norm=norm, cmap=cmap)\n",
-    "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)"
+    "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)\n"
    ]
   },
   {
@@ -167,7 +167,7 @@
     "    cachedir=cachedir,\n",
     "    cachename=\"layer_model.nc\",\n",
     ")\n",
-    "layer_model"
+    "layer_model\n"
    ]
   },
   {
@@ -186,7 +186,7 @@
    "outputs": [],
    "source": [
     "ds = nlmod.to_model_ds(layer_model, model_name, model_ws, delr=100.0, delc=100.0)\n",
-    "ds"
+    "ds\n"
    ]
   },
   {
@@ -208,7 +208,7 @@
    "outputs": [],
    "source": [
     "refinement_features = [(bgt[bgt[\"bronhouder\"] == \"L0002\"], 2)]\n",
-    "ds = nlmod.grid.refine(ds, refinement_features=refinement_features)"
+    "ds = nlmod.grid.refine(ds, refinement_features=refinement_features)\n"
    ]
   },
   {
@@ -226,7 +226,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds = nlmod.time.set_ds_time(ds, time=time)"
+    "ds = nlmod.time.set_ds_time(ds, time=time)\n"
    ]
   },
   {
@@ -245,7 +245,7 @@
    "outputs": [],
    "source": [
     "knmi_ds = nlmod.read.knmi.get_recharge(ds, cachedir=cachedir, cachename=\"recharge.nc\")\n",
-    "ds.update(knmi_ds)"
+    "ds.update(knmi_ds)\n"
    ]
   },
   {
@@ -289,7 +289,7 @@
     "oc = nlmod.gwf.oc(ds, gwf)\n",
     "\n",
     "# create storagee package\n",
-    "sto = nlmod.gwf.sto(ds, gwf)"
+    "sto = nlmod.gwf.sto(ds, gwf)\n"
    ]
   },
   {
@@ -340,7 +340,7 @@
     "lakes = bgt_grid[mask].copy()\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]"
+    "bgt_grid = bgt_grid[~mask]\n"
    ]
   },
   {
@@ -362,7 +362,7 @@
     "lek[\"stage\"] = 0.0\n",
     "lek[\"rbot\"] = -3.0\n",
     "spd = nlmod.gwf.surface_water.build_spd(lek, \"RIV\", ds)\n",
-    "riv = flopy.mf6.ModflowGwfriv(gwf, stress_period_data={0: spd})"
+    "riv = flopy.mf6.ModflowGwfriv(gwf, stress_period_data={0: spd})\n"
    ]
   },
   {
@@ -381,7 +381,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "drn = nlmod.gwf.surface_water.gdf_to_seasonal_pkg(bgt_grid, gwf, ds);"
+    "drn = nlmod.gwf.surface_water.gdf_to_seasonal_pkg(bgt_grid, gwf, ds)"
    ]
   },
   {
@@ -425,7 +425,7 @@
     "] = 1.0  # overstort hoogte\n",
     "\n",
     "# add lake to groundwaterflow model\n",
-    "nlmod.gwf.lake_from_gdf(gwf, lakes, ds, boundname_column=\"name\");"
+    "nlmod.gwf.lake_from_gdf(gwf, lakes, ds, boundname_column=\"name\")"
    ]
   },
   {
@@ -436,7 +436,7 @@
    "outputs": [],
    "source": [
     "# create recharge package\n",
-    "rch = nlmod.gwf.rch(ds, gwf)"
+    "rch = nlmod.gwf.rch(ds, gwf)\n"
    ]
   },
   {
@@ -454,7 +454,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "nlmod.sim.write_and_run(sim, ds)"
+    "nlmod.sim.write_and_run(sim, ds)\n"
    ]
   },
   {
@@ -473,7 +473,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "head = nlmod.gwf.get_heads_da(ds)"
+    "head = nlmod.gwf.get_heads_da(ds)\n"
    ]
   },
   {
@@ -497,7 +497,7 @@
     "    head.sel(layer=\"HLc\").mean(\"time\"), ds=ds, edgecolor=\"k\", norm=norm\n",
     ")\n",
     "cbar = nlmod.plot.colorbar_inside(pc)\n",
-    "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")"
+    "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")\n"
    ]
   },
   {
@@ -525,7 +525,7 @@
     "cbar = nlmod.plot.colorbar_inside(pc, bounds=[0.05, 0.05, 0.02, 0.9])\n",
     "dcs.plot_grid()\n",
     "dcs.plot_layers(alpha=0.0, min_label_area=1000)\n",
-    "f.tight_layout(pad=0.0)"
+    "f.tight_layout(pad=0.0)\n"
    ]
   },
   {
@@ -546,7 +546,7 @@
     "x = 118228\n",
     "y = 439870\n",
     "head_point = nlmod.gwf.get_head_at_point(head, x=x, y=y, ds=ds)\n",
-    "head_point.plot.line(hue=\"layer\", size=10);"
+    "head_point.plot.line(hue=\"layer\", size=10)"
    ]
   },
   {
@@ -566,7 +566,7 @@
    "source": [
     "df = pd.read_csv(os.path.join(model_ws, \"lak_STAGE.csv\"), index_col=0)\n",
     "df.index = ds.time.values\n",
-    "ax = df.plot(figsize=(10, 3));"
+    "ax = df.plot(figsize=(10, 3))"
    ]
   },
   {
@@ -679,7 +679,7 @@
     "        ha=\"center\",\n",
     "        va=\"top\",\n",
     "        transform=ax.transAxes,\n",
-    "    )"
+    "    )\n"
    ]
   },
   {
@@ -712,7 +712,7 @@
     "pg = nlmod.modpath.pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5)\n",
     "\n",
     "# create the modpath simulation file\n",
-    "mpsim = nlmod.modpath.sim(mpf, pg, \"forward\", gwf=gwf)"
+    "mpsim = nlmod.modpath.sim(mpf, pg, \"forward\", gwf=gwf)\n"
    ]
   },
   {
@@ -723,7 +723,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\")\n"
    ]
   },
   {
@@ -733,7 +733,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pdata = nlmod.modpath.load_pathline_data(mpf)"
+    "pdata = nlmod.modpath.load_pathline_data(mpf)\n"
    ]
   },
   {
@@ -776,7 +776,7 @@
     "line = ax.add_collection(lc)\n",
     "nlmod.plot.colorbar_inside(line, label=\"Travel time (years)\")\n",
     "\n",
-    "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\");"
+    "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")"
    ]
   },
   {
@@ -812,13 +812,19 @@
     "dcs.plot_grid()\n",
     "# add labels with layer names\n",
     "dcs.plot_layers(alpha=0.0, min_label_area=1000)\n",
-    "f.tight_layout(pad=0.0)"
+    "f.tight_layout(pad=0.0)\n"
    ]
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "nlmod",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
-   "name": "python"
+   "name": "python",
+   "version": "3.9.4"
   }
  },
  "nbformat": 4,
diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py
index 3cdbc39d..3e544dd4 100644
--- a/nlmod/gwf/recharge.py
+++ b/nlmod/gwf/recharge.py
@@ -39,19 +39,25 @@ 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
-    rch_name_arr, rch_unique_dic = _get_unique_series(ds, "recharge", pname)
-    ds["rch_name"] = ds["top"].dims, rch_name_arr
-    if mask is not None:
-        mask = (ds["rch_name"] != "") & mask
+    if ds.time.steady_state:
+        recharge = "recharge"
+        if "time" in ds["recharge"].dims:
+            mask = ds["recharge"].isel(time=0) != 0
+        else:
+            mask = ds["recharge"] != 0
     else:
-        mask = ds["rch_name"] != ""
-
-    recharge = "rch_name"
+        recharge = "rch_name"
+        rch_name_arr, rch_unique_dic = _get_unique_series(ds, "recharge", pname)
+        ds[recharge] = ds["top"].dims, rch_name_arr
+        if mask is not None:
+            mask = (ds["rch_name"] != "") & mask
+        else:
+            mask = ds["rch_name"] != ""
 
     spd = da_to_reclist(
         ds,
         mask,
-        col1=recharge,
+        col1=ds[recharge].squeeze(),
         first_active_layer=True,
         only_active_cells=False,
     )
@@ -67,6 +73,9 @@ def ds_to_rch(gwf, ds, mask=None, pname="rch", **kwargs):
         **kwargs,
     )
 
+    if ds.time.steady_state:
+        return rch
+
     # create timeseries packages
     _add_time_series(rch, rch_unique_dic, ds)
 
@@ -118,6 +127,7 @@ def ds_to_evt(gwf, ds, pname="evt", nseg=1, surface=None, depth=None, **kwargs):
         logger.info("Setting extinction depth to 1 meter below surface")
         depth = 1.0
 
+    # check for nan values
     if ds["evaporation"].isnull().any():
         raise ValueError("please remove nan values in evaporation data array")