diff --git a/pyproject.toml b/pyproject.toml
index 890fc7574..bf4bc089f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,7 +12,7 @@ version = "1.4.0a3"
 dependencies = [
     "numpy  > 1.24.4",
     "pandas >= 1.4, < 3.0", # TODO remove upper limit
-    "mikeio >= 1.2",
+    "mikeio >= 3.2",
     "matplotlib",
     "xarray",
     "netCDF4",
diff --git a/requirements_min.txt b/requirements_min.txt
deleted file mode 100644
index 0491e51c0..000000000
--- a/requirements_min.txt
+++ /dev/null
@@ -1,3 +0,0 @@
-numpy==1.23.0
-pandas==2.0.0
-xarray==2023.1.0
\ No newline at end of file
diff --git a/src/modelskill/model/dfsu.py b/src/modelskill/model/dfsu.py
index 7fec7bb4e..86d757cf6 100644
--- a/src/modelskill/model/dfsu.py
+++ b/src/modelskill/model/dfsu.py
@@ -1,5 +1,4 @@
 from __future__ import annotations
-import inspect
 from pathlib import Path
 from typing import Literal, get_args, cast
 
@@ -195,9 +194,7 @@ def _extract_vertical(
             if isinstance(self.data, mikeio.Dataset):
                 ds_column = self.data.sel(x=x, y=y)
             elif isinstance(self.data, mikeio.dfsu.Dfsu3D):
-                # FIXME: open with specific element instead...make sure mikeio fix is in place before
-                # ds_column = self.data.read(elements=elemids) # when bug is fixed in mikeio
-                ds_column = self.data.read(items=self.sel_items.all).sel(x=x, y=y)
+                ds_column = self.data.read(elements=elemids, items=self.sel_items.all)
             else:
                 raise ValueError(
                     "Unsupported data type for vertical profile extraction."
@@ -207,8 +204,8 @@ def _extract_vertical(
                 "Only spatial_method='contained' is currently implemented for vertical profile extraction from DfsuModelResult. "
             )
 
-        # get layer depth info
-        layer_boundaries = ds_column.geometry.calc_ze(ds_column._zn)
+        # z-coordinates at element centers, one row per timestep
+        element_depths = ds_column.z.elements
 
         item_name = self.sel_items.values
 
@@ -218,7 +215,7 @@ def _extract_vertical(
 
         # Create flattened arrays # Repeat each timestamp n_layers times
         time_flat = np.repeat(ds_column.time, n_layers)
-        z_flat = layer_boundaries.flatten()  # Flatten z-coordinates
+        z_flat = element_depths.flatten()  # Flatten z-coordinates
         item_values_1d = ds_column[item_name].to_numpy().flatten()  # Flatten item?
         # aux_items_flat = {aux_item: ds_column[aux_item].to_numpy().flatten() for aux_item in self.sel_items.aux}
 
@@ -261,8 +258,8 @@ def _extract_point(
             )
 
         if method == "contained":
-            signature = inspect.signature(self.data.geometry.find_index)
-            if "z" in signature.parameters and z is not None:
+            # only layered geometries can be indexed by z
+            if z is not None and hasattr(self.data.geometry, "n_layers"):
                 elemids = self.data.geometry.find_index(x=x, y=y, z=z)
             else:
                 elemids = self.data.geometry.find_index(x=x, y=y)
diff --git a/tests/model/test_vertical.py b/tests/model/test_vertical.py
index 7c86fe504..7957559fa 100644
--- a/tests/model/test_vertical.py
+++ b/tests/model/test_vertical.py
@@ -229,7 +229,7 @@ def test_extract_from_dfsu_correct_layers(self, dfsu_ds):
         # ===
 
         # expected element depths from dfsu geometry
-        element_depths_expected = dfsu_col.geometry.calc_ze(dfsu_col._zn)
+        element_depths_expected = dfsu_col.z.elements
 
         # element depths at first timestep from VerticalModelResult
         element_depths_t0 = vmr.data.sel(time=vmr.data.time.values[0]).z.values
@@ -239,6 +239,40 @@ def test_extract_from_dfsu_correct_layers(self, dfsu_ds):
         assert np.allclose(element_depths_t0, element_depths_expected[0, :])
         assert np.allclose(element_depths_tend, element_depths_expected[-1, :])
 
+    def test_extract_from_dfsu_file_matches_dataset(self, dfsu_fpath, dfsu_ds):
+        # A dfsu file takes the Dfsu3D branch, which reads only the observation
+        # column; a Dataset takes the sel(x, y) branch. Both must give the same
+        # z and values.
+        dfsu_mr = ms.DfsuModelResult(dfsu_fpath, item=0, name="test")
+        assert isinstance(dfsu_mr.data, mikeio.dfsu.Dfsu3D)
+
+        dummy_obs = pd.DataFrame(
+            {"z": [-5.0, -4.0, -3.0], "salt": [30.0, 31.0, 32.0]},
+            index=pd.to_datetime(["2022-06-14 00:00:00"] * 3),
+        )
+        XPOS = 6.575e5
+        YPOS = 6.55e6
+        vo = ms.VerticalObservation(dummy_obs, x=XPOS, y=YPOS, item="salt", z_item="z")
+
+        vmr = dfsu_mr.extract(vo, spatial_method="contained")
+        vmr_from_ds = ms.DfsuModelResult(dfsu_ds, item=0, name="test").extract(
+            vo, spatial_method="contained"
+        )
+
+        assert np.allclose(vmr.data.z.values, vmr_from_ds.data.z.values)
+        assert np.allclose(vmr.data["test"].values, vmr_from_ds.data["test"].values)
+
+        # ...and against the dfsu column itself, so both branches breaking the
+        # same way is still a failure
+        dfsu_col = dfsu_ds.sel(x=XPOS, y=YPOS)
+        item_name = dfsu_ds.items[0].name
+        assert np.allclose(vmr.data.z.values, dfsu_col.z.elements.flatten())
+        assert np.allclose(
+            vmr.data["test"].values, dfsu_col[item_name].to_numpy().flatten()
+        )
+        assert vmr.x == pytest.approx(dfsu_col.geometry.element_coordinates[0, 0])
+        assert vmr.y == pytest.approx(dfsu_col.geometry.element_coordinates[0, 1])
+
     @pytest.mark.parametrize("spatial_method", ["nearest", "inverse_distance"])
     def test_extract_from_dfsu_unsupported_spatial_methods_raise(
         self, dfsu_ds, spatial_method