diff --git a/docs/examples/18_observations.ipynb b/docs/examples/18_observations.ipynb
index c7026de4..627491da 100644
--- a/docs/examples/18_observations.ipynb
+++ b/docs/examples/18_observations.ipynb
@@ -351,14 +351,6 @@
     "ax.plot(obswell.x, obswell.y, \"ro\", markersize=10, label=obswell.name.item())\n",
     "ax.legend(loc=(0, 1), frameon=False, ncol=2, fontsize=\"small\");"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "9198f61f",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/docs/examples/18_particle_tracking_prt.ipynb b/docs/examples/19_particle_tracking_prt.ipynb
similarity index 100%
rename from docs/examples/18_particle_tracking_prt.ipynb
rename to docs/examples/19_particle_tracking_prt.ipynb
diff --git a/docs/examples/19_voronoi.ipynb b/docs/examples/20_voronoi.ipynb
similarity index 100%
rename from docs/examples/19_voronoi.ipynb
rename to docs/examples/20_voronoi.ipynb
diff --git a/nlmod/dims/layers.py b/nlmod/dims/layers.py
index d505a13c..cb83fddf 100644
--- a/nlmod/dims/layers.py
+++ b/nlmod/dims/layers.py
@@ -2278,7 +2278,7 @@ def get_modellayers_indexer(
     screen_top="screen_top",
     screen_bottom="screen_bottom",
     full_output=False,
-    drop_nan_layers=False,
+    drop_nan_layers=True,
 ):
     """Get a model layer indexer (dataset) for a dataframe with observation wells.
 
@@ -2308,7 +2308,7 @@ def get_modellayers_indexer(
         If False, return only the variables needed to index a data array.
     drop_nan_layers : bool, optional
         If True, drop the observation wells for which the model layer cannot be
-        determined. These probably lie above or below the model. The default is False.
+        determined. These probably lie above or below the model. The default is True.
         If False, the model layer indexer will contain NaN values for these wells.
 
     Returns
@@ -2404,9 +2404,12 @@ def get_modellayers_indexer(
         ].values
     elif drop_nan_layers:
         obs_ds = obs_ds.dropna(dim, subset=["modellayer"])
+        obs_ds["modellayer"].values = ds["layer"][
+            obs_ds["modellayer"].astype(int)
+        ].values
     else:
         logger.warning(
-            "There are observation wells that lie above/below the model. "
+            "There are observation wells above/below the model top/bottom. "
             "The model layer indexer will contain NaN values for these wells."
         )
 
diff --git a/tests/test_009_layers.py b/tests/test_009_layers.py
index 5d520504..209236d2 100644
--- a/tests/test_009_layers.py
+++ b/tests/test_009_layers.py
@@ -287,12 +287,12 @@ def test_insert_layer():
 
 
 def test_remove_thin_layers():
-    # %% download regis and define min_thickness
+    # download regis and define min_thickness
     regis = get_regis_horstermeer()
 
     min_thickness = 1.0
 
-    # %% test update_thickness_every_layer = False
+    # test update_thickness_every_layer = False
     ds_new = nlmod.layers.remove_thin_layers(regis, min_thickness)
 
     th = nlmod.layers.calculate_thickness(regis)
@@ -307,7 +307,7 @@ def test_remove_thin_layers():
     # test if all active cells in the new dataset were active in the original dataset
     assert (th.data[th_new > 0] > 0).all()
 
-    # %% test update_thickness_every_layer = True
+    # test update_thickness_every_layer = True
     ds_new2 = nlmod.layers.remove_thin_layers(
         regis, min_thickness, update_thickness_every_layer=True
     )
@@ -358,27 +358,31 @@ def test_get_modellayers_indexer():
     # structured grid
     idx = nlmod.layers.get_modellayers_indexer(ds, df)
     # check result
-    assert np.isnan(idx["layer"].values[0])
-    assert idx["layer"].values[1] == 1.0
-    assert idx["layer"].values[2] == ds.sizes["layer"] - 1
+    assert idx["layer"].values[0] == ds['layer'].values[1]
+    assert idx["layer"].values[1] == ds['layer'].values[ds.sizes["layer"] - 1]
+    # test getting bottom elevations using indexer
+    _ = ds["botm"].sel(**idx)
+
+    # structured grid keep nan
+    idx2 = nlmod.layers.get_modellayers_indexer(ds, df, drop_nan_layers=False)
+    # check result
+    assert np.isnan(idx2["layer"].values[0])
 
     # drop nans
-    idx = idx.dropna("index")
-    # get layer names (annoying step that maybe should be performed in the function)
-    idx["layer"].values = ds.layer[idx["layer"].astype(int)].values
+    idx2 = idx2.dropna("index")
+    # get layer names (step is unnecesary if you keep drop_nan_layers=True)
+    idx2["layer"].values = ds.layer[idx2["layer"].astype(int)].values
     # test getting bottom elevations using indexer
-    _ = ds["botm"].sel(**idx)
+    _ = ds["botm"].sel(**idx2)
 
     # vertex grid
     ds_ref = nlmod.grid.refine(ds, refinement_features=[])
-    idx2 = nlmod.layers.get_modellayers_indexer(ds_ref, df)
-    idx2 = idx2.dropna("index")
-    idx2["layer"].values = ds_ref.layer[idx2["layer"].astype(int)].values
-    _ = ds_ref["botm"].sel(**idx2)
+    idx3 = nlmod.layers.get_modellayers_indexer(ds_ref, df)
+    _ = ds_ref["botm"].sel(**idx3)
 
-    assert (idx2["layer"] == idx["layer"]).all()
+    assert (idx3["layer"] == idx["layer"]).all()
 
     # full output
     idxfull = nlmod.layers.get_modellayers_indexer(ds, df, full_output=True)
-    assert (idxfull["modellayer_top"] == np.array([np.inf, 1, 4, 0])).all()
-    assert (idxfull["modellayer_bot"] == np.array([np.inf, 2, 4, 4])).all()
+    assert (idxfull["modellayer_top"] == np.array([1, 4, 0])).all()
+    assert (idxfull["modellayer_bot"] == np.array([2, 4, 4])).all()