diff --git a/docs/examples/howto_meshplotlib/plot_animation.py b/docs/examples/howto_meshplotlib/plot_animation.py
index 929c96620f440e27d8619d38eed3be22a3d29483..02ecb6ba0bb4ef76f19cdbc0d97855200a669718 100644
--- a/docs/examples/howto_meshplotlib/plot_animation.py
+++ b/docs/examples/howto_meshplotlib/plot_animation.py
@@ -4,22 +4,43 @@ How to create Animations
 
 .. sectionauthor:: Florian Zill (Helmholtz Centre for Environmental Research GmbH - UFZ)
 
-For this example we load a 2D meshseries dataset from within the ``meshplotlib`` examples.
-...
+To demonstrate the creation of an animated plot we use a component transport
+example from the ogs benchmark gallery
+(https://www.opengeosys.org/docs/benchmarks/hydro-component/elder/).
 """
 # %%
 import numpy as np
 
 from ogstools.meshplotlib import examples, setup
 from ogstools.meshplotlib.animation import animate
-from ogstools.propertylib.defaults import temperature
+from ogstools.propertylib import ScalarProperty
 
-setup.p_max = 48.5
-setup.num_levels = 25
-# TODO: needs better example
-mesh_series = examples.meshseries_THM_2D
+mesh_series = examples.meshseries_CT_2D
+# %%
+# Let's use fixed scale limits to prevent rescaling during the animation.
+setup.p_min = 0
+setup.p_max = 1
+
+# %%
+# You can choose which timesteps to render by passing either an int array
+# corresponding to the indices, or a float array corresponding to the timevalues
+# to render. If a requested timevalue is not part of the timeseries it will be
+# interpolated.
+timevalues = np.linspace(
+    mesh_series.timevalues[0], mesh_series.timevalues[-1], num=25
+)
 
-timevalues = np.linspace(0, mesh_series.timevalues[-1], num=5)
-titles = [f"{tv/(365.25*86400):.1f} yrs" for tv in timevalues]
-anim = animate(mesh_series, temperature, timevalues, titles)
 # %%
+# Now, let's animate the saturation solution. A timescale at the top
+# indicates existing timesteps and the position of the current timevalue.
+titles = [f"{tv/(365.25*86400):.1f} yrs" for tv in timevalues]
+si = ScalarProperty("Si", "", "%", "Saturation")
+anim = animate(mesh_series, si, timevalues, titles)
+
+# sphinx_gallery_start_ignore
+# note for developers:
+# unfortunately when creating the documentation the animation is run twice:
+# once for saving a thumbnail and another time for generating the html_repr
+# see .../site-packages/sphinx_gallery/scrapers.py", line 234, 235
+# I don't see a way to circumvent this.
+# sphinx_gallery_end_ignore
diff --git a/ogstools/meshplotlib/animation.py b/ogstools/meshplotlib/animation.py
index f4515e56ae7df4a6ed0e0d56833ed74995d9bff7..0012275618592186b48b1207bdf14886c68ccdb3 100644
--- a/ogstools/meshplotlib/animation.py
+++ b/ogstools/meshplotlib/animation.py
@@ -17,15 +17,15 @@ from .mesh_series import MeshSeries
 
 
 def timeline(ax: plt.Axes, x: float, xticks: list[float]) -> None:
-    y = 1.07
+    y = 1.1
     ap = {"arrowstyle": "-", "color": "k"}
     axfr = "axes fraction"
     ax.annotate("", (1, y), (0, y), axfr, arrowprops=ap)
     align = dict(ha="center", va="center")  # noqa: C408: noqa
     for xt in xticks:
-        ax.annotate("|", (xt, y), (xt, y), axfr, **align, size=6)
+        ax.annotate("|", (xt, y), (xt, y), axfr, **align, size=28)
     align = dict(ha="center", va="center")  # noqa: C408: noqa
-    style = dict(color="g", size=18, weight="bold")  # noqa: C408: noqa
+    style = dict(color="g", size=36, weight="bold")  # noqa: C408: noqa
     ax.annotate("|", (x, y), (x, y), axfr, **align, **style)
 
 
@@ -49,7 +49,11 @@ def animate(
     start_time = time.time()
     ts = mesh_series.timesteps if timesteps is None else timesteps
     tv = mesh_series.timevalues
-    xticks = [t / tv[-1] for t in tv]
+
+    def t_frac(t):
+        return (t - tv[0]) / (tv[-1] - tv[0])
+
+    xticks = [t_frac(t) for t in tv]
     fig = plot(mesh_series.read(0, False), property)
 
     def init() -> None:
@@ -72,11 +76,11 @@ def animate(
         if isinstance(i, int):
             mesh = mesh_series.read(i)
         else:
-            mesh = mesh_series.read_interp(i, False)
+            mesh = mesh_series.read_interp(i, True)
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             fig = _plot(fig, mesh, property)
-            x = i / len(ts) if isinstance(i, int) else i / tv[-1]
+            x = i / len(ts) if isinstance(i, int) else t_frac(i)
             timeline(fig.axes[0], x, xticks)
 
     _func = partial(animate_func, fig=fig)