diff --git a/docs/user-guide/data-structures.qmd b/docs/user-guide/data-structures.qmd
index f5bcd5b52..5cc6e312f 100644
--- a/docs/user-guide/data-structures.qmd
+++ b/docs/user-guide/data-structures.qmd
@@ -20,31 +20,37 @@ The primary data of ModelSkill are the data that needs to be compared: observati
 
 * `point`: 0D time series data
 * `track`: 0D time series data at moving locations (trajectories)
+* `vertical`: 1D depth profile time series at a fixed location
 * `grid`: gridded 2D data
 * `dfsu`: flexible mesh 2D data
 
-Point and track data are both `TimeSeries` objects, while grid and dfsu data are both `SpatialField` objects. `TimeSeries` objects are ready to be compared whereas data from `SpatialField` object needs to be *extracted* first (the extracted object will be of the `TimeSeries` type).
+Point, track and vertical data are all `TimeSeries` objects, while grid and dfsu data are both `SpatialField` objects. `TimeSeries` objects are ready to be compared whereas data from `SpatialField` object needs to be *extracted* first (the extracted object will be of the `TimeSeries` type).
 
 `TimeSeries` objects contains its data in an [](`xarray.Dataset`) with the actual data in the first DataArray and optional auxilliary data in the following DataArrays. The DataArrays have a `kind` attribute with either `observation` or `model`.
 
+Vertical data is stored in *long* format: a single `time` dimension in which each timestamp is repeated once per depth, with the depth stored in a non-index `z` coordinate.
+
 
 ## Comparer objects
 
 Comparer objects are results of a matching procedure (between observations and model results) or constructed directly from already matched data. A comparison of a *single* observation and one or more model results are stored in a [](`~modelskill.Comparer`) object. A comparison of *multiple* observations and one or more model results are stored in a [](`~modelskill.ComparerCollection`) object which is a collection of [](`~modelskill.Comparer`) objects.
 
-The matched data in a [](`~modelskill.Comparer`) is stored in an [](`xarray.Dataset`) which can be accessed via the `data` attribute. The Dataset has an attribute `gtype` which is a string describing the type of data (e.g. `point`, `track`). The first DataArray in the Dataset is the observation data, the next DataArrays are model result data and optionally additional DataArrays are auxilliarye data. Each of the DataArrays have a `kind` attribute with either `observation`, `model` or `aux`.
+The matched data in a [](`~modelskill.Comparer`) is stored in an [](`xarray.Dataset`) which can be accessed via the `data` attribute. The Dataset has an attribute `gtype` which is a string describing the type of data (e.g. `point`, `track`, `vertical`). The first DataArray in the Dataset is the observation data, the next DataArrays are model result data and optionally additional DataArrays are auxilliarye data. Each of the DataArrays have a `kind` attribute with either `observation`, `model` or `aux`.
 
 Both [](`~modelskill.Comparer`) and [](`~modelskill.ComparerCollection`) have a `plot` accessor for plotting the data (e.g. `cmp.plot.timeseries()` or `cmp.plot.scatter()`).
 
+A [](`~modelskill.Comparer`) with `gtype="vertical"` additionally has a `vertical` accessor (e.g. `cmp.vertical.skill()` or `cmp.vertical.plot.hovmoller()`) for depth-binned skill and profile plots.
+
 
 
 ## Skill objects
 
-Calling a skill method on a comparer object will return a skill object with skill scores (statistics) from comparing observation and model result data using different metrics (e.g. root mean square error). Two skill objects are currently implemented: [](`~modelskill.SkillTable`) and [](`~modelskill.skill_grid.SkillGrid`). The first is relevant for all ModelSkill users while the latter is relevant for users of the track data (e.g. MetOcean studies using satellite altimetry data).
+Calling a skill method on a comparer object will return a skill object with skill scores (statistics) from comparing observation and model result data using different metrics (e.g. root mean square error). Three skill objects are currently implemented: [](`~modelskill.SkillTable`), [](`~modelskill.skill_grid.SkillGrid`) and [](`~modelskill.skill_profile.SkillProfile`). The first is relevant for all ModelSkill users, while the second is relevant for users of track data (e.g. MetOcean studies using satellite altimetry data) and the third for users of vertical profile data.
 
 If `c` is a comparer object, then the following skill methods are available:
 
 * `c.skill()` -> `SkillTable`
 * `c.mean_skill()` -> `SkillTable`
 * `c.gridded_skill()` -> `SkillGrid`
+* `c.vertical.skill()` -> `SkillProfile` (only for vertical data)
 
diff --git a/docs/user-guide/matching.qmd b/docs/user-guide/matching.qmd
index 8b0b583fd..fec696a77 100644
--- a/docs/user-guide/matching.qmd
+++ b/docs/user-guide/matching.qmd
@@ -91,7 +91,7 @@ cc = ms.match([o1, o2], mr_dfsu, spatial_method='contained')
 
 ::: callout-note
 * Extraction of *track* data does not currently support the "contained" method.
-* Extraction of point data from 3D dfsu files is not yet fully supported. It is recommended to extract the data "offline" prior to using ModelSkill.
+* Extraction of *point* data from 3D dfsu files is not yet fully supported. It is recommended to extract the data "offline" prior to using ModelSkill. Extraction of whole *profiles* from 3D dfsu files is supported, see [matching of vertical profiles](#matching-of-vertical-profiles).
 :::
 
 
@@ -104,6 +104,43 @@ cc = ms.match([o1, o2], mr_nc, spatial_method='nearest')
 ```
 
 
+## Matching of vertical profiles
+
+A [](`~modelskill.VerticalObservation`) (e.g. a CTD cast) is matched with either a
+[](`~modelskill.VerticalModelResult`) - a column extracted beforehand - or a *3D*
+[](`~modelskill.DfsuModelResult`), in which case the column at the observation
+position is extracted as part of the matching.
+
+::: {.callout-note}
+Vertical profile support is still in development and the API may change.
+:::
+
+```{python}
+vo = ms.VerticalObservation("../data/vertical/VerticalProfile_obs1.dfs0",
+                            item="Salinity", z_item="z", x=657500, y=6553600)
+mr_3d = ms.model_result("../data/vertical/sigma_z_coast.dfsu", item="Salinity")
+cmp = ms.match(vo, mr_3d)
+cmp
+```
+
+Matching is done in two steps:
+
+1. The model profile is taken from the model time *closest* to each observation time. The tolerance defaults to half the median model time step, so observation times without a nearby model profile are dropped.
+2. The model values are interpolated *linearly in depth* onto the observation depths. Observation depths outside the model depth range are discarded - no extrapolation is performed.
+
+Note that the temporal step is a nearest-neighbour selection rather than the
+linear interpolation in time used for point and track data. This preserves each
+modelled profile as a whole instead of blending two profiles.
+
+The result is a [](`~modelskill.Comparer`) with `gtype="vertical"` in *long* format:
+one row per (time, depth) pair.
+
+::: callout-note
+* Extraction from a 3D dfsu file only supports `spatial_method="contained"` (the default); the other methods are not implemented for vertical profiles.
+* The `max_model_gap` argument is not yet supported for vertical matching.
+:::
+
+
 ## Event-based matching and handling of gaps
 
 If the model result data contains gaps either because only events are stored or because of missing data, the `max_model_gap` argument of the [](`~modelskill.match`) function can be used to specify the maximum allowed gap (in seconds) in the model result data. This will avoid interpolating model data over long gaps in the model result data!
diff --git a/docs/user-guide/overview.qmd b/docs/user-guide/overview.qmd
index e531d03dc..c9f11a2a3 100644
--- a/docs/user-guide/overview.qmd
+++ b/docs/user-guide/overview.qmd
@@ -18,7 +18,7 @@ they will need to be defined and then matched in space and time with the `match(
 
 ### Define observations
 
-The observations can be defined as either a [`PointObservation`](`modelskill.PointObservation`) or a [`TrackObservation`](`modelskill.TrackObservation`) (a moving point).
+The observations can be defined as a [`PointObservation`](`modelskill.PointObservation`), a [`TrackObservation`](`modelskill.TrackObservation`) (a moving point) or a [`VerticalObservation`](`modelskill.VerticalObservation`) (a depth profile at a fixed position).
 
 ```{python}
 import modelskill as ms
@@ -34,8 +34,9 @@ A model result will either be a simple point/track like the observations, or spa
 
 * [`PointModelResult`](`modelskill.PointModelResult`) - a point result from a dfs0/nc file or a DataFrame
 * [`TrackModelResult`](`modelskill.TrackModelResult`) - a track result from a dfs0/nc file or a DataFrame
+* [`VerticalModelResult`](`modelskill.VerticalModelResult`) - a depth profile from a dfs0 file or a DataFrame
 * [`GridModelResult`](`modelskill.GridModelResult`) - a spatial field from a dfs2/nc file or a Xarray Dataset
-* [`DfsuModelResult`](`modelskill.DfsuModelResult`) - a spatial field from a dfsu file
+* [`DfsuModelResult`](`modelskill.DfsuModelResult`) - a spatial field from a dfsu file (2D or 3D)
 
 ```{python}
 mr1 = ms.PointModelResult("../data/model.dfs0", item="WL")
@@ -51,6 +52,24 @@ cc
 ```
 
 
+## Vertical profiles
+
+Model results can also be validated against vertical profiles, e.g. CTD casts of
+salinity or temperature. The observation is a
+[`VerticalObservation`](`modelskill.VerticalObservation`) and the model result is either a
+[`VerticalModelResult`](`modelskill.VerticalModelResult`) (a pre-extracted column) or a
+3D [`DfsuModelResult`](`modelskill.DfsuModelResult`), from which the column at the
+observation position is *extracted* automatically (see [matching](matching.qmd)).
+
+The result is an ordinary [`Comparer`](`modelskill.Comparer`) with `gtype="vertical"`,
+so `cmp.skill()` and all the usual plots work as before. In addition, a `vertical`
+accessor gives access to depth-aware skill and plots (see [plotting](plotting.qmd)).
+The comparer can be sliced by depth with `cmp.sel(z=-2)` or `cmp.sel(z=slice(-10, 0))`
+(see [selecting data](selecting-data.qmd)), and the profile can be collapsed to a
+single value per timestep with `cmp.vertical.mean()`, `.min()` or `.max()` (which
+return point comparers).
+
+
 ## Analysis
 
 Once the observations and model results are matched, the `Comparer` object can be used for analysis and plotting. 
diff --git a/docs/user-guide/plotting.qmd b/docs/user-guide/plotting.qmd
index 5fdce145d..4ded728d7 100644
--- a/docs/user-guide/plotting.qmd
+++ b/docs/user-guide/plotting.qmd
@@ -72,6 +72,55 @@ cmp.plot.scatter();
 ```
 
 
+## Vertical profiles
+
+A [](`~modelskill.Comparer`) holding vertical profile data (`gtype="vertical"`)
+has an extra `vertical` accessor with depth-aware plots, available through
+[`cmp.vertical.plot`](`modelskill.comparison.VerticalPlotter`).
+
+The `profile` plot shows model and observed values against depth. It is meant
+for a single point in time, so select a time first:
+
+```{python}
+vo = ms.VerticalObservation("../data/vertical/VerticalProfile_obs1.dfs0",
+                            item="Salinity", z_item="z", x=657500, y=6553600)
+mr3d = ms.DfsuModelResult("../data/vertical/sigma_z_coast.dfsu", item="Salinity")
+vcmp = ms.match(vo, mr3d)
+
+vcmp.sel(time="2022-06-15").vertical.plot.profile();
+```
+
+The model is drawn as a line (the full modelled column) and the observations as
+markers. Pass `show_matched_model=True` to also mark the model values
+interpolated to the observation depths. `vertical.plot()` is a shorthand for
+`vertical.plot.profile()`.
+
+The `hovmoller` plot shows the whole comparison as a function of depth and time,
+with the model as a filled contour and the observations overlaid as markers:
+
+```{python}
+vcmp.vertical.plot.hovmoller(figsize=(12, 4));
+```
+
+With more than one model result, the model to contour must be named explicitly,
+e.g. `vcmp.vertical.plot.hovmoller(model="m1")`.
+
+Skill computed in depth bins with
+[`vertical.skill()`](`modelskill.comparison.VerticalAccessor.skill`) returns a
+[](`~modelskill.skill_profile.SkillProfile`), whose metrics can be plotted as a
+horizontal bar plot per depth bin:
+
+```{python}
+sk = vcmp.vertical.skill(bins=4)
+sk.rmse.plot();
+```
+
+The bins can also be given explicitly, e.g. `bins=[-20, -15, -10, -5, 0]`.
+
+Note that all the ordinary Comparer plots shown above work on vertical data too;
+they simply ignore the depth dimension.
+
+
 ## Taylor diagrams
 
 A Taylor diagram shows how well a model result matches an observation in terms of correlation, standard deviation and root mean squared error. The `taylor` plot can be accessed through the Comparer [`plot`](`modelskill.comparison.ComparerPlotter`) accessor or the ComparerCollection [`plot`](`modelskill.comparison.ComparerCollectionPlotter`) accessor:
diff --git a/docs/user-guide/selecting-data.qmd b/docs/user-guide/selecting-data.qmd
index 609a0b370..6d95383a9 100644
--- a/docs/user-guide/selecting-data.qmd
+++ b/docs/user-guide/selecting-data.qmd
@@ -53,13 +53,14 @@ cmp = ms.match(o, [m1, m2])
 
 The [](`~modelskill.Comparer.sel`) method allows you to select data based on specific criteria such as time, model name, or spatial area. It returns a new `Comparer` object with the selected data. This method is highly versatile and supports multiple selection parameters, which can be combined.
 
-**Syntax:** `Comparer.sel(model=None, time=None, area=None)`
+**Syntax:** `Comparer.sel(model=None, time=None, area=None, z=None)`
 
 | Parameter    | Type                        | Description                                             | Default |
 |--------------|-----------------------------|---------------------------------------------------------|---------|
 | `model`      | str, int, or list           | Model name or index. Selects specific models.           | None    |
 | `time`       | str, datetime, or slice     | Specific time or range for selection.                  | None    |
 | `area`       | list of float or Polygon    | Bounding box [x0, y0, x1, y1] or a polygon area filter. | None    |
+| `z`          | float or slice              | Depth or depth range (vertical data only).              | None    |
 
 **Example 1: Selecting data by time**
 ```{python}
@@ -84,6 +85,43 @@ cmp_area = cmp.sel(area=[4.0, 52.5, 5.0, 53.0])
 
 This filters the data within the bounding box defined by `[x0, y0, x1, y1]`.
 
+**Example 4: Selecting by depth**
+
+For vertical profile data (`gtype="vertical"`), `sel()` additionally accepts a
+`z` argument, given either as a single depth or as a slice:
+
+```{python}
+#| code-fold: true
+#| code-summary: "Vertical observation and model data"
+vo = ms.VerticalObservation("../data/vertical/VerticalProfile_obs1.dfs0",
+                            item="Salinity", z_item="z", x=657500, y=6553600)
+mr3d = ms.DfsuModelResult("../data/vertical/sigma_z_coast.dfsu", item="Salinity")
+vcmp = ms.match(vo, mr3d)
+```
+
+```{python}
+vcmp.sel(z=slice(-10, 0))
+```
+
+A slice keeps all observation depths within the range. A single value instead
+selects the nearest depth:
+
+```{python}
+vcmp.sel(z=-2)
+```
+
+Selections can be combined and given a descriptive name, which is convenient
+when comparing skill between depth ranges:
+
+```{python}
+cmp_surface = vcmp.sel(z=slice(-10, 0)).rename({vcmp.name: "surface"})
+cmp_deep = vcmp.sel(z=slice(None, -10)).rename({vcmp.name: "deep"})
+ms.ComparerCollection([cmp_surface, cmp_deep]).skill()
+```
+
+`z` also works on a [](`~modelskill.ComparerCollection`), and depth can equally
+be filtered with `query()`, e.g. `vcmp.query("z > -10")`.
+
 ### `where()` method
 
 The [](`~modelskill.Comparer.where`) method is used to filter data conditionally. It works similarly to `xarray`'s `where` method and returns a new `Comparer` object with values satisfying a given condition. Other values will be masked (set to `NaN`).
@@ -94,14 +132,14 @@ The [](`~modelskill.Comparer.where`) method is used to filter data conditionally
 |-----------|------------------------------|-----------------------------------------------|
 | `cond`    | bool, np.ndarray, or xr.DataArray | Condition to filter values (True or False). |
 
-**Example 4: Filtering data conditionally**
+**Example 5: Filtering data conditionally**
 ```{python}
 cmp.where(cmp.data.Observation > 3)
 ```
 
 This filters out any rows where the observation values are not greater than 3.
 
-**Example 5: Multiple conditions**
+**Example 6: Multiple conditions**
 ```{python}
 cmp.where((cmp.data.m1 < 2.9) & (cmp.data.Observation > 3))
 ```
@@ -120,7 +158,7 @@ The [](`~modelskill.Comparer.query`) method uses a [](`pandas.DataFrame.query`)-
 |-----------|--------|---------------------------------------------|
 | `query`   | str    | Query string for filtering data.           |
 
-**Example 6: Querying data**
+**Example 7: Querying data**
 ```{python}
 cmp.query("Observation > 3.0 and m1 < 2.9")
 ```
@@ -143,14 +181,14 @@ sk = cmp.skill(metrics=["rmse", "mae", "si"])
 sk
 ```
 
-**Example 7: Select model**
+**Example 8: Select model**
 ```{python}
 sk.sel(model='m1')
 ```
 
 Here, `sk` contains skill scores for all models, and `sk_m1` filters the results to include only model "m1". Observations can be selected in the same way.
 
-**Example 8: Querying skill scores**
+**Example 9: Querying skill scores**
 ```{python}
 sk_high_rmse = sk.query("rmse > 0.3")
 sk_high_rmse
@@ -158,7 +196,7 @@ sk_high_rmse
 
 This filters the `SkillTable` to include only rows where the root mean square error (RMSE) exceeds 0.3.
 
-**Example 9: Accessing and visualizing specific metrics**
+**Example 10: Accessing and visualizing specific metrics**
 ```{python}
 sk_rmse = sk.rmse
 sk_rmse
diff --git a/src/modelskill/comparison/_vertical_comparison.py b/src/modelskill/comparison/_vertical_comparison.py
index a272d2dfe..552e4d31c 100644
--- a/src/modelskill/comparison/_vertical_comparison.py
+++ b/src/modelskill/comparison/_vertical_comparison.py
@@ -23,7 +23,7 @@ def __init__(self, comparer):
         self.comparer = comparer
 
     def __call__(self, *args, **kwargs) -> matplotlib.axes.Axes:
-        """Plot scatter plot of modelled vs observed data"""
+        """Plot vertical profiles of model and observations (alias for `profile`)"""
         return self.profile(*args, **kwargs)
 
     def profile(
@@ -43,6 +43,9 @@ def profile(
             Matplotlib Axes to plot on (if None, a new figure and axes will be created).
         figsize : tuple, optional
             Size of the figure (only used if ax is None).
+        show_matched_model : bool, optional
+            Also show the model values interpolated to the observation depths
+            as markers, by default False (only the raw model profile is drawn).
 
         Returns
         -------
@@ -284,15 +287,27 @@ def _agg(self, agg_func: str = "mean"):
         return Comparer(ds, raw_mod_data=raw_mod_data)
 
     def mean(self):
-        """Aggregate the comparison vertically using a specified aggregation function."""
+        """Aggregate the profile to a single mean value per timestep.
+
+        Returns a point `Comparer`. The raw model profile is restricted to the
+        depth range covered by the observation before aggregating.
+        """
         return self._agg(agg_func="mean")
 
     def min(self):
-        """Aggregate the comparison vertically using a specified aggregation function."""
+        """Aggregate the profile to a single min value per timestep.
+
+        Returns a point `Comparer`. The raw model profile is restricted to the
+        depth range covered by the observation before aggregating.
+        """
         return self._agg(agg_func="min")
 
     def max(self):
-        """Aggregate the comparison vertically using a specified aggregation function."""
+        """Aggregate the profile to a single max value per timestep.
+
+        Returns a point `Comparer`. The raw model profile is restricted to the
+        depth range covered by the observation before aggregating.
+        """
         return self._agg(agg_func="max")
 
     def skill(
diff --git a/src/modelskill/model/vertical.py b/src/modelskill/model/vertical.py
index 5b571469a..214086e2c 100644
--- a/src/modelskill/model/vertical.py
+++ b/src/modelskill/model/vertical.py
@@ -12,31 +12,77 @@
 
 
 class VerticalModelResult(TimeSeries):
-    """Model result for a vertical column.
+    """Model result for a vertical profile at a fixed (x, y) location.
 
-    Construct a VerticalColumnModelResult from a dfs0 file,
-    mikeio.Dataset, pandas.DataFrame or a xarray.Datasets
+    Construct a VerticalModelResult from a dfs0 file, mikeio.Dataset,
+    pandas.DataFrame or xarray.Dataset in long format: one row per
+    (time, z) pair, with one item/column holding the vertical coordinate
+    and another holding the modelled value. At least two items are required
+    (z + value); if more are present, `item` must be given.
 
     Parameters
     ----------
-    data : str, Path, pd.DataFrame, mikeio.Dfs0, mikeio.Dfs0, xr.Dataset
-        The input data or file path
-    name : str | None, optional
-        The name of the model result,
-        by default None (will be set to file name or item name)
-    item : str | int | None, optional
-        If multiple items/arrays are present in the input an item
-        must be given (as either an index or a string), by default None
-    z_item : str | int | None, optional
-        Item of the first coordinate of positions, by default None
-    x : float, optional
-        lateral coordinate of point position, inferred from data if not given, else None
-    y : float, optional
-        zonal coordinate of point position, inferred from data if not given, else None
+    data : str, Path, pd.DataFrame, mikeio.Dfs0, mikeio.Dataset, xr.Dataset
+        Input data or path to a dfs0 file.
+    name : str, optional
+        Name of the model result, by default the file or item name.
+    item : str or int, optional
+        Index or name of the value item. Required if the input has more
+        than two items.
     quantity : Quantity, optional
         Model quantity, for MIKE files this is inferred from the EUM information
-    aux_items : list[int | str] | None, optional
-        Auxiliary items, by default None
+    z_item : str or int, optional
+        Index or name of the item holding the vertical coordinate, by default 0.
+    x : float, optional
+        x-coordinate of the profile location, inferred from data when possible.
+    y : float, optional
+        y-coordinate of the profile location, inferred from data when possible.
+    aux_items : list[int | str], optional
+        Auxiliary items to keep alongside the value item, by default None.
+
+    Notes
+    -----
+    A dfs0 with N depth levels has its profile timestamps repeated N times on
+    a non-equidistant time axis. Duplicate (time, z) pairs are not allowed and
+    will raise a ValueError.
+
+    Examples
+    --------
+    From a `pandas.DataFrame` in long format:
+
+    ```{python}
+    import modelskill as ms
+    import pandas as pd
+
+    times = pd.to_datetime(["2010-01-01 01:00"] * 3 + ["2010-01-01 02:00"] * 3)
+    df = pd.DataFrame(
+        {
+            "z": [0.0, -5.0, -10.0, 0.0, -5.0, -10.0],
+            "Salinity": [30.1, 30.3, 30.4, 30.5, 30.3, 30.3],
+        },
+        index=times,
+    )
+    ms.VerticalModelResult(
+        df,
+        item="Salinity",
+        z_item="z",
+        x=12.0,
+        y=55.0,
+        quantity=ms.Quantity("Salinity", "PSU"),
+    )
+    ```
+
+    From a dfs0 file (with z, Salinity and Temperature items):
+
+    ```{python}
+    ms.VerticalModelResult(
+        "../data/vertical/VerticalModel_at_obs.dfs0",
+        item="Salinity",
+        z_item="z",
+        x=657500,
+        y=6553600,
+    )
+    ```
     """
 
     def __init__(
diff --git a/src/modelskill/obs.py b/src/modelskill/obs.py
index e06b79858..fa2cc01ca 100644
--- a/src/modelskill/obs.py
+++ b/src/modelskill/obs.py
@@ -367,70 +367,81 @@ def __init__(
 
 
 class VerticalObservation(Observation):
-    """Class for observations of vertical profiles.
+    """Observation of a vertical profile at a fixed (x, y) location.
 
-    Create a VerticalObservation from a dfs0/nc file or tabular data
-    containing time, vertical coordinate, and observed values.
+    Create a VerticalObservation from a dfs0 file, mikeio.Dataset,
+    pandas.DataFrame or xarray.Dataset in long format: one row per
+    (time, z) pair, with one item/column holding the vertical coordinate
+    and another holding the observed value. At least two items are required
+    (z + value); if more are present, `item` must be given.
 
     Parameters
     ----------
-    data : (str, Path, pd.DataFrame, mikeio.Dfs0, mikeio.Dataset, xr.Dataset)
-        Input data with vertical profile observations.
+    data : str, Path, pd.DataFrame, mikeio.Dfs0, mikeio.Dataset, xr.Dataset
+        Input data or path to a dfs0 file.
     item : int or str, optional
-        Index or name of the primary observation item.
-        If the input contains more than one candidate value item,
-        this argument must be provided.
+        Index or name of the value item. Required if the input has more
+        than two items.
     x : float, optional
-        x-coordinate of the observation location. If not provided,
-        it is inferred from data when possible.
+        x-coordinate of the profile location, inferred from data when possible.
     y : float, optional
-        y-coordinate of the observation location. If not provided,
-        it is inferred from data when possible.
+        y-coordinate of the profile location, inferred from data when possible.
     z_item : int or str, optional
-        Index or name of the vertical coordinate item, by default 0.
+        Index or name of the item holding the vertical coordinate, by default 0.
     name : str, optional
-        User-defined name for identification in plots and summaries.
+        Name of the observation, by default the file or item name.
     weight : float, optional
-        Weighting factor for skill scores, by default 1.0.
+        Weighting factor for skill scores in a ComparerCollection, by default 1.0.
     quantity : Quantity, optional
-        Physical quantity metadata used for validation against model results.
+        Observed quantity, for MIKE files this is inferred from the EUM information
     aux_items : list[int | str], optional
-        List of auxiliary item names or indices to keep in the dataset.
+        Auxiliary items to keep alongside the value item, by default None.
     attrs : dict, optional
         Additional attributes to be added to the underlying dataset.
 
+    Notes
+    -----
+    A dfs0 with N depth levels has its profile timestamps repeated N times on
+    a non-equidistant time axis. Duplicate (time, z) pairs are not allowed and
+    will raise a ValueError.
+
     Examples
     --------
-    >>> import modelskill as ms
-    >>> import pandas as pd
-    >>> df = pd.DataFrame(
-    ...     {
-    ...         "z": [0.0, -5.0, -10.0, 0.0, -5.0, -10.0],
-    ...         "value": [0.1, 0.3, 0.4, 0.5, 0.3, 0.3],
-    ...     },
-    ...     index=pd.to_datetime(
-    ...         [
-    ...             "2010-01-01 01:00:00",
-    ...             "2010-01-01 01:00:00",
-    ...             "2010-01-01 01:00:00",
-    ...             "2010-01-01 02:00:00",
-    ...             "2010-01-01 02:00:00",
-    ...             "2010-01-01 02:00:00",
-    ...         ]
-    ...     ),
-    ... )
-    >>> df.index.name = "t"
-    >>> print(df.to_string())
-                           z  value
-    t
-    2010-01-01 01:00:00   0.0    0.1
-    2010-01-01 01:00:00  -5.0    0.3
-    2010-01-01 01:00:00 -10.0    0.4
-    2010-01-01 02:00:00   0.0    0.5
-    2010-01-01 02:00:00  -5.0    0.3
-    2010-01-01 02:00:00 -10.0    0.3
-
-    >>> o = ms.VerticalObservation(df, item="value", z_item="z", x=12.0, y=55.0)
+    From a `pandas.DataFrame` in long format:
+
+    ```{python}
+    import modelskill as ms
+    import pandas as pd
+
+    times = pd.to_datetime(["2010-01-01 01:00"] * 3 + ["2010-01-01 02:00"] * 3)
+    df = pd.DataFrame(
+        {
+            "z": [0.0, -5.0, -10.0, 0.0, -5.0, -10.0],
+            "Salinity": [30.0, 30.2, 30.3, 30.4, 30.2, 30.2],
+        },
+        index=times,
+    )
+    ms.VerticalObservation(
+        df,
+        item="Salinity",
+        z_item="z",
+        x=12.0,
+        y=55.0,
+        quantity=ms.Quantity("Salinity", "PSU"),
+    )
+    ```
+
+    From a dfs0 file (with z and Salinity items):
+
+    ```{python}
+    ms.VerticalObservation(
+        "../data/vertical/VerticalProfile_obs1.dfs0",
+        item="Salinity",
+        z_item="z",
+        x=657500,
+        y=6553600,
+    )
+    ```
     """
 
     def __init__(