diff --git a/apps/showcases/Antidunes/pyvista.ipynb b/apps/showcases/Antidunes/pyvista.ipynb
index 1f74046eda657ec9e4d2cbb80e94a2ca07b431b3..6cb28649aca82ce6aa152f0b56a3a49552f6c00a 100644
--- a/apps/showcases/Antidunes/pyvista.ipynb
+++ b/apps/showcases/Antidunes/pyvista.ipynb
@@ -119,7 +119,7 @@
    "id": "5d1248ea",
    "metadata": {},
    "source": [
-    "## Plot velocity field "
+    "## Plot velocity field (side view)"
    ]
   },
   {
@@ -222,7 +222,7 @@
    "id": "8c1f1df5",
    "metadata": {},
    "source": [
-    "## Plot 3D velocity field for graphical abstract"
+    "## Plot velocity field (3D view) for graphical abstract"
    ]
   },
   {
@@ -348,7 +348,7 @@
     "## Plot full velocity field at each sampling point (for creating the animations)\n",
     "\n",
     "### Warning: there are severe bugs in pyvista that overcomplicate this script:\n",
-    "1. memory leak in pl.screenshot(); memory leak can be avoided with pl.clear(); \n",
+    "1. memory leak in pl.screenshot(); memory leak is supposed to be circumvented with pl.clear();\n",
     "2. HOWEVER: when pl.clear() was called, pl.view_xz() somehow disrupts pl.camera.zoom() such that the latter command does not work anymore\n",
     "\n",
     "=> new strategy: call python from bash such that the kernel restarts for every file; this avoids the memory leak"
@@ -432,6 +432,7 @@
     "\n",
     "# empirically set after viewing data in ParaView (in lattice units)\n",
     "color_bar_limits_velocity = np.array([-0.3, 1]) # E1 with e_dry=0.97\n",
+    "color_bar_limits_velocity = np.array([-0.3, 1.4]) # E4 with e_dry=0.97\n",
     "\n",
     "fluid_file_directory = os.fsencode(fluid_file_basepath)\n",
     "\n",
@@ -473,37 +474,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "537dc005",
+   "id": "a5b66e66",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "#!/bin/bash -l\n",
-    "\n",
-    "conda activate\n",
-    "\n",
-    "fluid_file_basepath1=\"/simdata/on74yces/2502347/vtk-out/fluid_field/\"\n",
-    "# fluid_file_basepath2=\"/simdata/on74yces/2504475/vtk-out/fluid_field/\"\n",
-    "# fluid_file_basepath3=\"/simdata/on74yces/2509267/vtk-out/fluid_field/\"\n",
-    "# fluid_file_basepath4=\"/simdata/on74yces/2519434/vtk-out/fluid_field/\"\n",
-    "\n",
-    "declare -a PathList=($fluid_file_basepath1 $fluid_file_basepath2 $fluid_file_basepath3 $fluid_file_basepath4)\n",
-    "\n",
-    "for path in ${PathList[@]}; do\n",
-    "   for file in $path*.vtu; do\n",
-    "      python3 create-image.py $file\n",
-    "      # echo $file\n",
-    "   done\n",
-    "done"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "ab05ec15",
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
    "source": [
     "# this Python script (named create_image.py) must be called by the bash script\n",
     "\n",
@@ -515,13 +488,16 @@
     "import pyvista as pv\n",
     "import numpy as np\n",
     "\n",
+    "pv.start_xvfb()\n",
+    "\n",
     "# taken from job output file\n",
-    "dx = 0.00025 \n",
+    "dx = 0.00025\n",
     "dt = 1.38e-5 # E1\n",
     "#dt = 1.08125e-5 # E4\n",
     "\n",
     "# empirically set after viewing data in ParaView (in lattice units)\n",
     "color_bar_limits_velocity = np.array([-0.3, 1]) # E1 with e_dry=0.97\n",
+    "#color_bar_limits_velocity = np.array([-0.3, 1.4]) # E4 with e_dry=0.97\n",
     "\n",
     "# path to fluid file must be given as first command line argument\n",
     "fluid_file_path = str(sys.argv[1])\n",
@@ -551,7 +527,7 @@
     "particle_mesh.point_data[\"diameter_si\"] = diameter_si\n",
     "\n",
     "# create glyphs for particles\n",
-    "sphere_glyphs = particle_mesh.glyph(scale=\"radius\", factor=1, \n",
+    "sphere_glyphs = particle_mesh.glyph(scale=\"radius\", factor=1,\n",
     "                                    geom=pv.Sphere(radius=1, theta_resolution=150, phi_resolution=150))\n",
     "\n",
     "pl = pv.Plotter(lighting='three lights', off_screen=True)\n",
@@ -582,6 +558,249 @@
     "pl.window_size = [10000, 1000]\n",
     "pl.screenshot(output_file_path,transparent_background=False)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "537dc005",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#!/bin/bash -l\n",
+    "\n",
+    "conda activate\n",
+    "\n",
+    "fluid_file_basepath1=\"/simdata/on74yces/2502347/vtk-out/fluid_field/\"\n",
+    "# fluid_file_basepath2=\"/simdata/on74yces/2504475/vtk-out/fluid_field/\"\n",
+    "# fluid_file_basepath3=\"/simdata/on74yces/2509267/vtk-out/fluid_field/\"\n",
+    "# fluid_file_basepath4=\"/simdata/on74yces/2519434/vtk-out/fluid_field/\"\n",
+    "\n",
+    "declare -a PathList=($fluid_file_basepath1 $fluid_file_basepath2 $fluid_file_basepath3 $fluid_file_basepath4)\n",
+    "\n",
+    "for path in ${PathList[@]}; do\n",
+    "   for file in $path*.vtu; do\n",
+    "      python3 create-image.py $file\n",
+    "      # echo $file\n",
+    "   done\n",
+    "done"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ebe75f76",
+   "metadata": {},
+   "source": [
+    "## Plot full bed height elevation at each sampling point (for creating the animations)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9254f001",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "# this Python script (named create_image.py) must be called by the bash script\n",
+    "\n",
+    "#!/usr/bin/env python3\n",
+    "\n",
+    "import os\n",
+    "import sys\n",
+    "from PIL import Image\n",
+    "import pyvista as pv\n",
+    "import numpy as np\n",
+    "\n",
+    "pv.start_xvfb()\n",
+    "\n",
+    "# taken from job output file\n",
+    "dx = 0.00025\n",
+    "\n",
+    "# path to particle file must be given as first command line argument\n",
+    "particle_file_path = \"/simdata/on74yces/2509268/vtk-out/particles/simulation_step_3799000.vtu\" #str(sys.argv[1])\n",
+    "\n",
+    "filenumber = particle_file_path.split('_')[-1]\n",
+    "filenumber = filenumber.split('.')[0]\n",
+    "filenumber = str(filenumber).rjust(10, '0')\n",
+    "output_file_path = \"/simdata/ca36xymo/antidunes/animation-particles-e4/\" + filenumber + \".jpg\"\n",
+    "\n",
+    "particle_reader = pv.get_reader(particle_file_path)\n",
+    "particle_mesh = particle_reader.read()\n",
+    "\n",
+    "# add z-coordinate of particle center in SI units\n",
+    "particle_mesh.point_data[\"center_coordinate_z_si\"] = particle_mesh.points[:,2] * dx\n",
+    "\n",
+    "# create glyphs for particles\n",
+    "sphere_glyphs = particle_mesh.glyph(scale=\"radius\", factor=1, \n",
+    "                                    geom=pv.Sphere(radius=1, theta_resolution=30, phi_resolution=30))\n",
+    "\n",
+    "pl = pv.Plotter(lighting='three lights', off_screen=True)\n",
+    "\n",
+    "# add mesh at periodic sides to enlarge bed region\n",
+    "transform_matrix = np.array([[1, 0, 0, 0],\n",
+    "                                     [0, 1, 0, -60],\n",
+    "                                     [0, 0, 1, 0],\n",
+    "                                     [0, 0, 0, 1]])\n",
+    "particle_mesh_transformed_once = particle_mesh.transform(transform_matrix,inplace=False)\n",
+    "sphere_glyphs_transformed_once = particle_mesh_transformed_once.glyph(scale=\"radius\", factor=1, \n",
+    "                                    geom=pv.Sphere(radius=1, theta_resolution=30, phi_resolution=30))\n",
+    "particle_mesh_transformed_twice = particle_mesh_transformed_once.transform(transform_matrix,inplace=False)\n",
+    "sphere_glyphs_transformed_twice = particle_mesh_transformed_twice.glyph(scale=\"radius\", factor=1, \n",
+    "                                    geom=pv.Sphere(radius=1, theta_resolution=30, phi_resolution=30))\n",
+    "\n",
+    "pl.add_mesh(sphere_glyphs_transformed_once, scalars=\"center_coordinate_z_si\", show_scalar_bar=False)\n",
+    "pl.update_scalar_bar_range(clim=[0.008,0.012])\n",
+    "pl.add_mesh(sphere_glyphs, scalars=\"center_coordinate_z_si\", show_scalar_bar=False)\n",
+    "pl.update_scalar_bar_range(clim=[0.008,0.012])\n",
+    "pl.add_mesh(sphere_glyphs_transformed_twice, scalars=\"center_coordinate_z_si\", show_scalar_bar=False)\n",
+    "pl.add_scalar_bar(title=\"Bed height elevation / m\", vertical=False, width=0.15, height=0.15,\n",
+    "                  position_x=0.425, position_y=0.7, color='black', title_font_size=70, label_font_size=50,\n",
+    "                  fmt=\"%.3f\")\n",
+    "pl.update_scalar_bar_range(clim=[0.008,0.012])\n",
+    "\n",
+    "# pl.add_mesh(pv.Line((0,0,100), (3200,0,100)), color=\"lightgray\", line_width=10)\n",
+    "# pl.add_mesh(pv.Line((0,-60,100), (3200,-60,100)), color=\"lightgray\", line_width=10)\n",
+    "\n",
+    "pl.view_xy()\n",
+    "#pl.enable_parallel_projection()\n",
+    "pl.set_background('white')\n",
+    "#pl.camera.zoom(10)\n",
+    "\n",
+    "# specify camera so that image stays at a fixed position\n",
+    "pl.camera.clipping_range = (5812.79, 6720.05)\n",
+    "pl.camera.distance = 6207.82\n",
+    "pl.camera.focal_point = (1600.44, 29.9889, 38.3558)\n",
+    "pl.camera.parallel_projection = True\n",
+    "pl.camera.parallel_scale = 160.67\n",
+    "pl.camera.thickness = 907.259\n",
+    "pl.camera.view_angle = 30\n",
+    "\n",
+    "pl.window_size = [10000, 1000]\n",
+    "pl.screenshot(output_file_path,transparent_background=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "73d07351",
+   "metadata": {},
+   "source": [
+    "## Plot full free-surface height elevation at each sampling point (for creating the animations)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "95013d8f",
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "# this Python script (named create_image.py) must be called by the bash script\n",
+    "\n",
+    "#!/usr/bin/env python3\n",
+    "\n",
+    "import os\n",
+    "import sys\n",
+    "from PIL import Image\n",
+    "import pyvista as pv\n",
+    "import numpy as np\n",
+    "\n",
+    "pv.start_xvfb()\n",
+    "\n",
+    "# taken from job output file\n",
+    "dx = 0.00025\n",
+    "\n",
+    "# path to fluid file must be given as first command line argument\n",
+    "fluid_file_path = str(sys.argv[1]) #\"/simdata/on74yces/2517354/vtk-out/fluid_field/simulation_step_7551000.vtu\" \n",
+    "\n",
+    "filenumber = fluid_file_path.split('_')[-1]\n",
+    "filenumber = filenumber.split('.')[0]\n",
+    "filenumber = str(filenumber).rjust(10, '0')\n",
+    "output_file_path = \"/simdata/ca36xymo/antidunes/animation-free-surface-e4/\" + filenumber + \".jpg\"\n",
+    "\n",
+    "fluid_reader = pv.get_reader(fluid_file_path)\n",
+    "fluid_mesh = fluid_reader.read()\n",
+    "\n",
+    "# remove gas data from fluid_mesh\n",
+    "fluid_mesh = fluid_mesh.threshold(value=2, scalars=\"mapped_flag\", invert=True)\n",
+    "\n",
+    "# slice fluid mesh to get real 2D data\n",
+    "fluid_mesh = fluid_mesh.slice(normal=[0, 1, 0])\n",
+    "\n",
+    "# extract surface so we apply the extrusion\n",
+    "fluid_mesh = fluid_mesh.extract_surface()\n",
+    "\n",
+    "# extrude fluid_mesh to 3D\n",
+    "fluid_mesh = fluid_mesh.extrude([0, -60, 0], capping=False)\n",
+    "\n",
+    "# add z-coordinate of cell center in SI units\n",
+    "fluid_mesh.point_data[\"coordinate_z_si\"] = fluid_mesh.points[:,2] * dx\n",
+    "\n",
+    "pl = pv.Plotter(lighting='three lights', off_screen=True)\n",
+    "\n",
+    "pl.add_mesh(fluid_mesh, scalars=\"coordinate_z_si\", component=0, show_scalar_bar=False)\n",
+    "pl.add_scalar_bar(title=\"Free-surface elevation / m\", vertical=False, width=0.15, height=0.15,\n",
+    "                  position_x=0.425, position_y=0.7, color='black', title_font_size=70, label_font_size=50,\n",
+    "                  fmt=\"%.3f\")\n",
+    "#pl.update_scalar_bar_range(clim=[0.018,0.022]) # E1\n",
+    "pl.update_scalar_bar_range(clim=[0.020,0.024]) # E4\n",
+    "\n",
+    "pl.view_xy()\n",
+    "#pl.enable_parallel_projection()\n",
+    "pl.set_background('white')\n",
+    "#pl.camera.zoom(10)\n",
+    "\n",
+    "# specify camera so that image stays at a fixed position\n",
+    "pl.camera.clipping_range = (5812.79, 6720.05)\n",
+    "pl.camera.distance = 6207.82\n",
+    "pl.camera.focal_point = (1600.44, 29.9889, 38.3558)\n",
+    "pl.camera.parallel_projection = True\n",
+    "pl.camera.parallel_scale = 160.67\n",
+    "pl.camera.thickness = 907.259\n",
+    "pl.camera.view_angle = 30\n",
+    "\n",
+    "pl.window_size = [10000, 1000]\n",
+    "\n",
+    "pl.screenshot(output_file_path,transparent_background=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b516bfd5",
+   "metadata": {},
+   "source": [
+    "## Plot coordinate systems\n",
+    "This is very troublesome to do in the actual plots => add them afterwards using ffmpeg"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e46cfd29",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pyvista as pv\n",
+    "\n",
+    "output_file_path = \"/simdata/ca36xymo/antidunes/animation-free-surface-e1/coordinate-system.jpg\"\n",
+    "\n",
+    "pl = pv.Plotter(lighting='three lights', off_screen=True)\n",
+    "pl.view_xz()\n",
+    "pl.add_axes(color=\"black\", x_color=\"black\", y_color=\"black\", z_color=\"black\", viewport=(0, 0, 1, 1))\n",
+    "pl.set_background('white')\n",
+    "\n",
+    "pl.screenshot(output_file_path,transparent_background=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "63082336",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {