diff --git a/Dockerfile b/Dockerfile
index e9c7a0f..3b50174 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -86,20 +86,15 @@ RUN git clone --recursive https://github.com/TUM-I5/ASAGI.git \
     && CC=mpicc CXX=mpicxx cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DSHARED_LIB=off -DSTATIC_LIB=on -DNONUMA=on \
     && make -j$(nproc) && make install
 
-RUN git clone https://github.com/uphoffc/ImpalaJIT.git \
-    && cd ImpalaJIT \
-    && mkdir build && cd build \
-    && cmake .. && make -j $(nproc) install 
-
 RUN git clone https://github.com/SeisSol/easi \
     && cd easi \
     && mkdir build && cd build \
-    && CC=mpicc CXX=mpicxx cmake .. -DEASICUBE=OFF -DLUA=ON -DCMAKE_PREFIX_PATH=/home/tools -DCMAKE_INSTALL_PREFIX=/home/tools -DASAGI=ON -DIMPALAJIT=ON .. \
+    && CC=mpicc CXX=mpicxx cmake .. -DEASICUBE=OFF -DLUA=ON -DCMAKE_PREFIX_PATH=/home/tools -DCMAKE_INSTALL_PREFIX=/home/tools -DASAGI=ON -DIMPALAJIT=OFF .. \
     && make -j$(nproc) && make install
 
 RUN pip install numpy && docker-clean
 
-RUN git clone --recursive --depth 1 --single-branch --branch v1.1.3 https://github.com/SeisSol/SeisSol.git \
+RUN git clone --recursive --depth 1 --single-branch --branch v1.3.0 https://github.com/SeisSol/SeisSol.git \
     && cd SeisSol \
     && mkdir build_hsw && cd build_hsw \
     && export PATH=$PATH:/home/tools/bin \
@@ -116,13 +111,7 @@ RUN cd SeisSol/preprocessing/science/rconv \
     && CC=mpicc CXX=mpicxx cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_PREFIX_PATH=/home/tools \
     && make -j$(nproc) && cp rconv /home/tools/bin/
 
-RUN git clone --recursive --depth 1 --single-branch --branch v2.2.7 https://github.com/SCOREC/core.git \
-    && cd core \
-    && mkdir build && cd build \
-    && cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release -DSCOREC_CXX_FLAGS="-Wno-error=array-bounds" \
-    && make -j$(nproc) && make install
-
-RUN git clone --recursive --branch v1.0.1 https://github.com/SeisSol/PUMGen.git \
+RUN git clone --recursive --branch v1.1.0 https://github.com/SeisSol/PUMGen.git \
     && cd PUMGen \
     && mkdir build && cd build \
     && cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release \
@@ -163,6 +152,8 @@ RUN conda install \
     matplotlib \
     gmsh \
     python-gmsh \
+    sympy \
+    pandas \
     && docker-clean
 ENV PATH=/home/tools/bin:$PATH
 ENV OMP_PLACES="cores"
diff --git a/Dockerfile_jupyterlab b/Dockerfile_jupyterlab
index 2d71c4f..fa05cfb 100644
--- a/Dockerfile_jupyterlab
+++ b/Dockerfile_jupyterlab
@@ -93,7 +93,7 @@ RUN git clone https://github.com/SeisSol/easi \
     && CC=mpicc CXX=mpicxx cmake .. -DEASICUBE=OFF -DLUA=ON -DCMAKE_PREFIX_PATH=/home/tools -DCMAKE_INSTALL_PREFIX=/home/tools -DASAGI=ON -DIMPALAJIT=OFF .. \
     && make -j$(nproc) && make install
 
-RUN pip install numpy && pip install git+https://github.com/SeisSol/PSpaMM.git@davschneller/compile-fixes && docker-clean
+RUN pip install numpy && pip install git+https://github.com/SeisSol/PSpaMM.git && docker-clean
 
 #RUN git clone --depth 1 --single-branch --branch main https://github.com/libxsmm/libxsmm.git \
 #    && cd libxsmm \
@@ -104,7 +104,7 @@ RUN pip install numpy && pip install git+https://github.com/SeisSol/PSpaMM.git@d
 #       -j$(nproc) generator \
 #    && cp bin/libxsmm_gemm_generator /home/tools/bin
 
-RUN git clone --recursive --depth 1 --single-branch --branch v1.1.3 https://github.com/SeisSol/SeisSol.git \
+RUN git clone --recursive --depth 1 --single-branch --branch v1.3.0 https://github.com/SeisSol/SeisSol.git \
     && cd SeisSol \
     && mkdir build_hsw && cd build_hsw \
     && export PATH=$PATH:/home/tools/bin \
@@ -121,13 +121,7 @@ RUN cd SeisSol/preprocessing/science/rconv \
     && CC=mpicc CXX=mpicxx cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_PREFIX_PATH=/home/tools \
     && make -j$(nproc) && cp rconv /home/tools/bin/
 
-RUN git clone --recursive --depth 1 --single-branch --branch v2.2.7 https://github.com/SCOREC/core.git \
-    && cd core \
-    && mkdir build && cd build \
-    && cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release -DSCOREC_CXX_FLAGS="-Wno-error=array-bounds" \
-    && make -j$(nproc) && make install
-
-RUN git clone --recursive https://github.com/SeisSol/PUMGen.git \
+RUN git clone --recursive --branch v1.1.0 https://github.com/SeisSol/PUMGen.git \
     && cd PUMGen \
     && mkdir build && cd build \
     && cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release \
@@ -169,6 +163,8 @@ RUN conda install \
     matplotlib \
     gmsh \
     python-gmsh \
+    sympy \
+    pandas \
     && docker-clean
 ENV PATH=/home/tools/bin:$PATH
 ENV OMP_PLACES="cores"
diff --git a/kaikoura/parametersLSW.par b/kaikoura/parametersLSW.par
index 3949d68..4edf97b 100644
--- a/kaikoura/parametersLSW.par
+++ b/kaikoura/parametersLSW.par
@@ -7,6 +7,8 @@ Tv=0.05
 !Attenuation parameters (ignored if not compiled with attenuation)
 FreqCentral=0.5
 FreqRatio=100
+numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov'
+numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov'
 /
 
 &IniCondition
diff --git a/kaikoura/parametersRS.par b/kaikoura/parametersRS.par
index 92fdbbb..9dd594e 100644
--- a/kaikoura/parametersRS.par
+++ b/kaikoura/parametersRS.par
@@ -7,6 +7,9 @@ Tv=0.05
 !Attenuation parameters (ignored if not compiled with attenuation)
 FreqCentral=0.5
 FreqRatio=100
+
+numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov'
+numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov'
 /
 
 &IniCondition
diff --git a/sulawesi/parametersLSW.par b/sulawesi/parametersLSW.par
index 8f0ca36..f2a5292 100644
--- a/sulawesi/parametersLSW.par
+++ b/sulawesi/parametersLSW.par
@@ -7,6 +7,8 @@ Tv=0.05
 !Attenuation parameters (ignored if not compiled with attenuation)
 FreqCentral=0.3
 FreqRatio=100
+numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov'
+numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov'
 /
 
 &IniCondition
diff --git a/sulawesi/parametersRS.par b/sulawesi/parametersRS.par
index 6bed675..753e173 100644
--- a/sulawesi/parametersRS.par
+++ b/sulawesi/parametersRS.par
@@ -7,6 +7,8 @@ Tv=0.05
 !Attenuation parameters (ignored if not compiled with attenuation)
 FreqCentral=0.3
 FreqRatio=100
+numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov'
+numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov'
 /
 
 &IniCondition
diff --git a/tpv13/parameters.par b/tpv13/parameters.par
index ee3be63..367c296 100644
--- a/tpv13/parameters.par
+++ b/tpv13/parameters.par
@@ -4,6 +4,8 @@ MaterialFileName = 'tpv12_13_material.yaml'
 !enable off-fault plasticity (ignored in Plasticity=0)
 Plasticity = 1
 Tv = 0.03
+numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov'
+numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov'
 /
 
 &IniCondition
@@ -34,6 +36,7 @@ printtimeinterval_sec = 1.0                    ! Time interval at which output w
 OutputMask = 1 1 1 1 1 1 1 1 1 1 1             ! turn on and off fault outputs
 refinement_strategy = 2
 refinement = 1
+!vtkorder = 3
 /
 
 !parameterize ascii fault file outputs
@@ -68,8 +71,9 @@ refinement = 1
 
 ! Free surface output
 SurfaceOutput = 1
-SurfaceOutputRefinement = 1
+SurfaceOutputRefinement = 0
 SurfaceOutputInterval = 2.0
+surfacevtkorder = 2
 
 pickdt = 0.005                       ! Pickpoint Sampling
 RFileName = 'tpv12_13_receivers.dat'      ! Record Points in extra file
diff --git a/tpv13/tpv12_13_fault.yaml b/tpv13/tpv12_13_fault.yaml
index 1eb8f2a..d2141e3 100644
--- a/tpv13/tpv12_13_fault.yaml
+++ b/tpv13/tpv12_13_fault.yaml
@@ -1,26 +1,30 @@
 !Switch
 [s_xx, s_yy, s_zz, s_xy, s_yz, s_xz]: !Include tpv12_13_initial_stress.yaml
-[mu_s, mu_d, d_c, cohesion]: !IdentityMap
-  components:
-    # Inside nucleation patch
-    - !AxisAlignedCuboidalDomainFilter
-      limits:
-        x: [-1500, 1500]
-        y: [-.inf, .inf]
-        z: [-11691.34295108992, -9093.266739736605]
-      components: !ConstantMap
-        map:
-          mu_s:        0.54
-          mu_d:        0.10
-          d_c:         0.50
-          cohesion: -200000
-    # Outside nucleation patch
-    - !ConstantMap
+[mu_d, d_c, cohesion]: !ConstantMap
       map:
-        mu_s:        0.70
         mu_d:        0.10
         d_c:         0.50
         cohesion: -200000
+[mu_s]: !LuaMap
+  returns: [mu_s]
+  function: |
+    function f(x)
+        mu_s = 0.7
+        nucleation_size = 3000.0
+        along_width_nucleation_center = 12000.0
+        sinDip = math.sin(60 * math.pi / 180)
+        x_in_nucleation = (math.abs(x["x"]) <= 0.5 * nucleation_size)
+        z_in_nucleation = (
+            math.abs(x["z"] + along_width_nucleation_center * sinDip) <= 
+            0.5 * nucleation_size * sinDip
+        )     
+        if x_in_nucleation and z_in_nucleation then
+            mu_s = 0.54
+        end
+        
+        return { mu_s = mu_s }
+    end
+        
 [Tnuc_n, Tnuc_s, Tnuc_d, forced_rupture_time]: !ConstantMap
     map:
         Tnuc_n: 0
diff --git a/tpv13/tpv12_13_initial_stress.yaml b/tpv13/tpv12_13_initial_stress.yaml
index 77c517f..3490980 100644
--- a/tpv13/tpv12_13_initial_stress.yaml
+++ b/tpv13/tpv12_13_initial_stress.yaml
@@ -11,9 +11,10 @@
         depth = math.abs(x["z"])
         s_max_minus_Pf = 9.8 * (2700.0 - 1000.0);
         if (depth<=11951.15) then
+          s3_to_s1 = 0.3496
           -- math.floor(a+0.5) is round
-          s_xx =  -0.01 * math.floor(100.0 * (0.5 * (1.0 + 0.3496) * s_max_minus_Pf) + 0.5) * depth
-          s_yy =  -0.01 * math.floor(100.0 *               (0.3496 * s_max_minus_Pf) + 0.5) * depth
+          s_xx =  -0.01 * math.floor(100.0 * (0.5 * (1.0 + s3_to_s1) * s_max_minus_Pf) + 0.5) * depth
+          s_yy =  -0.01 * math.floor(100.0 *               (s3_to_s1 * s_max_minus_Pf) + 0.5) * depth
           s_zz = -s_max_minus_Pf * depth
         else
           s_xx = -s_max_minus_Pf * depth
diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 5fecdcc..3af0c25 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -46,16 +46,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "!gmsh -3 tpv13_training.geo\n",
+    "!gmsh -3 tpv13_training.geo -setnumber dip 60\n",
     "# on Frontera with apptainer, replace with:\n",
     "# !mpirun apptainer run {\"~/my-training.sif\"} gmsh -3 tpv13_training.geo"
    ]
   },
   {
    "cell_type": "markdown",
+   "id": "d67515aa",
    "metadata": {},
    "source": [
-    "You should now have a file called \"tpv13_training.msh\".\n",
+    "You should now have a file called `tpv13_training.msh`.\n",
     "In the next step, we translate this file into the efficient HDF5 format, which is going to be read in by SeisSol. HDF5 stands for the Hierarchical Data Format version 5, and is an open-source file format that supports large, complex, heterogeneous data. HDF5 uses a \"file directory\" like structure that allows you to organize data within the file in many different structured ways, as you might do with files on your computer, and is built for fast I/O processing and storage.\n",
     "\n",
     "We use our tool pumgen for mesh generation for SeisSol, which is [open-source software](  https://github.com/SeisSol/PUMGen), see also [SeisSol's documentation]( https://seissol.readthedocs.io/en/latest/meshing-with-pumgen.html)."
@@ -67,16 +68,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "!mpirun -n 1 pumgen -s msh2 tpv13_training.msh\n",
+    "!mpirun -n 1 pumgen -s msh4 tpv13_training.msh\n",
     "# on Frontera with apptainer, replace with:\n",
-    "# !mpirun -n 1 apptainer run {\"~/my-training.sif\"} pumgen -s msh2 tpv13_training.msh"
+    "# !mpirun -n 1 apptainer run {\"~/my-training.sif\"} pumgen -s msh4 tpv13_training.msh"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The files \"tpv13_training.puml.h5\" and \"tpv13_training.xdmf\" were created.\n",
+    "The files `tpv13_training.puml.h5` and `tpv13_training.xdmf` were created.\n",
     "While the .h5 file is read by SeisSol, the .xdmf file can be used to visualize the mesh, as in the following."
    ]
   },
@@ -108,8 +109,8 @@
    "metadata": {},
    "source": [
     "### Exercises\n",
-    "* Inspect the tpv13_training.geo file and try to understand the meshing script. Consult the [Gmsh documentation](https://gmsh.info/doc/texinfo/gmsh.html) about the BooleanFragments operation (which may not be intuitive but is a powerful operation gmsh provides).\n",
-    "* Create a finer resolved mesh by changing the h_fault parameter. Do not change the .geo file but use the [*-setnumber*](https://gmsh.info/doc/texinfo/gmsh.html#Command_002dline-options) PARAMETER VALUE command line argument of above *gmsh* command. Try, for example, to mesh the nucleation patch with 300 m and the remaining fault with 600 m element edge length.\n",
+    "* Inspect the `tpv13_training.geo` file and try to understand the meshing script. Consult the [Gmsh documentation](https://gmsh.info/doc/texinfo/gmsh.html) about the BooleanFragments operation (which may not be intuitive but is a powerful operation gmsh provides).\n",
+    "* Create a finer resolved mesh by changing the `h_fault` parameter. Do not change the .geo file but use the [`-setnumber PARAMETER VALUE`](https://gmsh.info/doc/texinfo/gmsh.html#Command_002dline-options) command line argument of above *gmsh* command. Try, for example, to mesh the nucleation patch with 300 m and the remaining fault with 600 m element edge length.\n",
     "* Change the dip of the fault by changing the dip parameter."
    ]
   },
@@ -182,33 +183,78 @@
     "    mesh = pv.wrap(reader.GetOutput())\n",
     "    plotter = pv.Plotter(notebook=True)\n",
     "    plotter.set_background('white')\n",
-    "    plotter.add_mesh(mesh, cmap='Blues', scalars=var)\n",
+    "    #fix color range for rupture speed\n",
+    "    clim_arg = {\"clim\":(0, 6000)} if var=='Vr' else {}\n",
+    "    plotter.add_mesh(mesh, cmap='Blues', scalars=var, **clim_arg)\n",
     "    plotter.view_xz()\n",
-    "    plotter.show(jupyter_backend='panel')"
+    "    plotter.show(jupyter_backend='panel')\n",
+    "    plotter.close()"
    ]
   },
   {
    "cell_type": "markdown",
+   "id": "5291f06b-f8c9-473f-9361-eddfa993b0c1",
    "metadata": {},
    "source": [
-    "# Local-time stepping exercises"
+    "# Initial stress and fault strength\n",
+    "\n",
+    "The loading in tpv13 is defined as a stress tensor, which is converted into fault tractions by SeisSol.\n",
+    "In the next cell, we compute analytically the normal and shear tractions, and compute the relative prestress ratio, a key parameter for rupture dynamics, characterizing the initial stress.\n",
+    "By modifying params with the nucleation static friction (mu_s=0.54), verify that this area is critically stressed (R>1)."
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "One of the key features of SeisSol is an efficient local-time stepping scheme, which clusters elements with similar time steps. This method pays off for applications featuring big differences in element sizes, caused by geometric complexity ([Ulrich et al., 2021](https://eartharxiv.org/repository/view/39/)) and/or multi-physics such as fully coupled seismic-acoustic-tsunami simulations ([Krenz et al., 2021](https://arxiv.org/abs/2107.06640)) or poroelastic wave propagation ([Wolf et al., 2021](https://arxiv.org/abs/2108.10565)).  \n",
+    "import sympy as sp\n",
+    "import numpy as np\n",
     "\n",
-    "A time-step cluster $c \\in \\mathbb{N}_0$ has the time-step\n",
-    "$$\\Delta t_c = r^c \\Delta t_{\\text{min}},$$\n",
-    "where the rate $r \\in \\mathbb{N}$. In the parameter file the rate $r$ can be set with the *ClusteredLTS* parameter.\n",
+    "# Define symbolic variables\n",
+    "s3_to_s1, sigma_max, dip, mu_s, mu_d = sp.symbols(\"s3_to_s1 sigma_max dip mu_s mu_d\")\n",
     "\n",
-    "* Run SeisSol once with rate-2 local time-stepping (LTS) and once with global time-stepping (GTS).\n",
-    "* In the log, look for \"Elapsed time (via clock_gettime)\" for both LTS and GTS.\n",
-    "* Compute the speed-up and compare it to the theoretical speed-up due to LTS.\n",
+    "# stress tensor\n",
+    "stress = sp.diag((1 + s3_to_s1) * sigma_max / 2, s3_to_s1 * sigma_max, sigma_max)\n",
+    "# normal and along dip unit vectors\n",
+    "u_n = sp.Matrix([0, -sp.sin(dip), sp.cos(dip)])\n",
+    "udip = sp.Matrix([0, sp.cos(dip), sp.sin(dip)])\n",
     "\n",
-    "*Hint:* You may reduce the *EndTime* parameter for this exercise."
+    "# tractions\n",
+    "traction = stress * u_n\n",
+    "\n",
+    "sigma_n = traction.dot(u_n)\n",
+    "print(\"sigma_n:\", sigma_n)\n",
+    "sigma_d = traction.dot(udip)\n",
+    "print(\"sigma_d:\", sigma_d)\n",
+    "\n",
+    "R = (sigma_d - mu_d * sigma_n) / ((mu_s - mu_d) * sigma_n)\n",
+    "print(\"relative prestress ratio R:\", R)\n",
+    "\n",
+    "params = {\"dip\": 60 * np.pi / 180, \"s3_to_s1\": 0.3495, \"mu_s\": 0.7, \"mu_d\": 0.1}\n",
+    "print(\n",
+    "    f\"relative prestress ration R evaluated with {params}:\",\n",
+    "    R.subs(params),\n",
+    ")\n",
+    "\n",
+    "print(\"sigma_n:\", sigma_n.subs(params))\n",
+    "print(\"sigma_d:\", sigma_d.subs(params))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "de67024b-6854-4fa3-9fe4-13562d828801",
+   "metadata": {},
+   "source": [
+    "# Rupture speed analysis\n",
+    "\n",
+    "We analyze the mechanism causing supershear in this setup and disable it.\n",
+    "\n",
+    "* Examine the earthquake rupture speed and identify the conditions that trigger supershear rupture.\n",
+    "* Change the nucleation size from 3 to 1.2km in both `tpv12_13_fault.yaml` (update parameters nucleation_size) and `tpv13_training.geo` (update parameters `l_n` and `w_n`) to disable the current supershear triggering mechanism.\n",
+    "* You should now observe a delay in the supershear transition.\n",
+    "* Verify that you can suppress the supershear transition mechanism by decreasing the initial shear stress, that is changing `s3_to_s1` from 0.3495 to 0.4 in `tpv12_13_initial_stress.yaml` and reenabling nucleation by decreasing `mu_s` in the nucleation area from 0.54 to 0.43 and increasing nucleation size from 1.2 to 1.8 km (in `tpv12_13_fault.yaml` and `tpv13_training.geo`)."
    ]
   },
   {
@@ -219,9 +265,9 @@
     "Now we compare the subsidence and uplift of dynamic rupture earthquake scenarios varying the dip of the fault\n",
     "\n",
     "To do so, perform the following steps:  \n",
-    "* Generate a new mesh with an alternative dip angle of your choice.\n",
-    "* Run a new simulation using the same command as above but with its output saved in a different folder, and change the prefix that all files will have (OutputFile) in the &Output namelist \n",
-    "* Add the new prefix to the prefixList of the code block below\n",
+    "* Generate a new mesh with an alternative dip angle of your choice. In addition to the mesh, don't forget to update the dip value in `tpv12_13_fault.yaml`, and adapt `mu_s` for the nucleation to remain overstressed (use the python code in the section \"Initial stress\").\n",
+    "* Run a new simulation using the same command as above but with its output saved in a different folder, and change the prefix that all files will have (`OutputFile`) in the `&Output` namelist \n",
+    "* Add the new prefix to the `prefixList` of the code block below\n",
     "* Analyse differences using the visualisation widget above and the cross-sections of surface displacements below\n"
    ]
   },
@@ -231,6 +277,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import pandas as pd\n",
+    "\n",
     "prefixList = ['tpv13']\n",
     "nSim = len(prefixList)\n",
     "# coordinates of the line end points\n",
@@ -238,21 +287,20 @@
     "b = [0, 10e3, 0]\n",
     "# create line polydata for vizualizing\n",
     "line = pv.Line(a, b)\n",
-    "p = pv.Plotter(shape=(1, nSim), notebook=True)\n",
+    "p = pv.Plotter(shape=(1, nSim), notebook=True, window_size=[1200, 500])\n",
     "\n",
     "def load_mesh(prefix):\n",
-    "    # Read xdmf data\n",
-    "    reader = vtk.vtkXdmfReader()\n",
-    "    reader.SetFileName(f'output/{prefix}-surface.xdmf')\n",
+    "    # Read free-surface output data\n",
+    "    reader = vtk.vtkHDFReader()\n",
+    "    reader.SetFileName(f'output/{prefix}-free-surface-4.vtkhdf')\n",
     "    reader.Update()\n",
-    "    reader.UpdateTimeStep(5.)\n",
     "    return pv.wrap(reader.GetOutput())\n",
     "\n",
     "# plot the line of the elevation profile over the surface displacement \n",
     "for i, pref in enumerate(prefixList):\n",
     "    mesh = load_mesh(pref)\n",
-    "\n",
     "    p.subplot(0, i)\n",
+    "    p.add_text(pref, font_size=10, position=\"upper_left\")\n",
     "    p.set_background('lightgrey')\n",
     "    p.add_mesh(mesh, scalars='u3')\n",
     "    p.add_mesh(line, color=\"white\", line_width=10)\n",
@@ -261,17 +309,56 @@
     "    )\n",
     "    p.camera.zoom(3)\n",
     "p.show(jupyter_backend='panel')\n",
-    "    \n",
-    "# plot subsidence/uplift profile\n",
+    "p.close()\n",
+    "\n",
+    "# plot fault output\n",
+    "p = pv.Plotter(shape=(1, nSim), notebook=True, window_size=[1200, 350])\n",
     "for i, pref in enumerate(prefixList):\n",
-    "    load_mesh(pref).plot_over_line(\n",
-    "        a,\n",
-    "        b,\n",
-    "        resolution=10000,\n",
-    "        title=f\"subsidence/uplift profile for {pref}\",\n",
-    "        ylabel=\"vectical displacement (m)\",\n",
-    "        scalars='u3'\n",
-    "    )"
+    "    reader = vtk.vtkXdmfReader()\n",
+    "    reader.SetFileName(f'output/{pref}-fault.xdmf')\n",
+    "    reader.Update()\n",
+    "    reader.UpdateTimeStep(9)\n",
+    "    mesh = pv.wrap(reader.GetOutput())\n",
+    "    p.subplot(0, i)\n",
+    "    p.add_text(pref, font_size=10, position=\"upper_left\")\n",
+    "    p.set_background('white')\n",
+    "    p.add_mesh(mesh, cmap='plasma', scalars=\"ASl\")\n",
+    "    p.view_xz()\n",
+    "    p.camera.zoom(1.6)  \n",
+    "p.show(jupyter_backend='panel')\n",
+    "p.close()\n",
+    "\n",
+    "# plot subsidence/uplift profiles\n",
+    "def plot_combined_profiles(prefixList, a, b, resolution=10000):\n",
+    "    plt.figure(figsize=(6, 4))\n",
+    "    for prefix in prefixList:\n",
+    "        mesh = load_mesh(prefix)\n",
+    "        sampled_data = mesh.sample_over_line(a, b, resolution=resolution)\n",
+    "        distance = sampled_data['Distance']\n",
+    "        displacement = sampled_data['u3']\n",
+    "        plt.plot(distance/1000, displacement, label=prefix, linewidth=2)\n",
+    "    plt.axhline(y=0,color='k', linestyle='--', linewidth=0.5)\n",
+    "    plt.title(\"Comparison of subsidence/uplift profiles\")\n",
+    "    plt.xlabel(\"Distance (km)\")\n",
+    "    plt.ylabel(\"Vertical Displacement (m)\")\n",
+    "    plt.legend()\n",
+    "    plt.tight_layout()\n",
+    "    plt.show()\n",
+    "plot_combined_profiles(prefixList, a, b)\n",
+    "\n",
+    "# plot moment rate release\n",
+    "fig, ax = plt.subplots(figsize=(6, 4))\n",
+    "for i, pref in enumerate(prefixList):\n",
+    "    df = pd.read_csv(f'output/{pref}-energy.csv')\n",
+    "    df = df.pivot_table(index=\"time\", columns=\"variable\", values=\"measurement\")\n",
+    "    df[\"seismic_moment_rate\"] = np.gradient(df[\"seismic_moment\"], df.index[1])\n",
+    "    ax.plot(df.index, df[\"seismic_moment_rate\"], label=pref)\n",
+    "\n",
+    "ax.set_xlabel('Time (s)')\n",
+    "ax.set_ylabel('Seismic moment rate (Nm/s)')\n",
+    "ax.set_title('Seismic moment rate comparison')\n",
+    "ax.legend()\n",
+    "plt.show()"
    ]
   },
   {
@@ -322,6 +409,30 @@
     "\n",
     "![](paraview.jpg)\n"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Local-time stepping exercises"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "One of the key features of SeisSol is an efficient local-time stepping scheme, which clusters elements with similar time steps. This method pays off for applications featuring big differences in element sizes, caused by geometric complexity ([Ulrich et al., 2021](https://eartharxiv.org/repository/view/39/)) and/or multi-physics such as fully coupled seismic-acoustic-tsunami simulations ([Krenz et al., 2021](https://arxiv.org/abs/2107.06640)) or poroelastic wave propagation ([Wolf et al., 2021](https://arxiv.org/abs/2108.10565)).  \n",
+    "\n",
+    "A time-step cluster $c \\in \\mathbb{N}_0$ has the time-step\n",
+    "$$\\Delta t_c = r^c \\Delta t_{\\text{min}},$$\n",
+    "where the rate $r \\in \\mathbb{N}$. In the parameter file the rate $r$ can be set with the *ClusteredLTS* parameter.\n",
+    "\n",
+    "* Run SeisSol once with rate-2 local time-stepping (LTS) and once with global time-stepping (GTS).\n",
+    "* In the log, look for \"Elapsed time (via clock_gettime)\" for both LTS and GTS.\n",
+    "* Compute the speed-up and compare it to the theoretical speed-up due to LTS.\n",
+    "\n",
+    "*Hint:* You may reduce the *EndTime* parameter for this exercise."
+   ]
   }
  ],
  "metadata": {
diff --git a/tpv13/tpv13_training.geo b/tpv13/tpv13_training.geo
index b08e8cf..6865c03 100644
--- a/tpv13/tpv13_training.geo
+++ b/tpv13/tpv13_training.geo
@@ -10,7 +10,7 @@ SetFactory("OpenCASCADE");
 
 // Length of the fault
 l_f = 30e3;
-// Depth of the fault
+// Width of the fault (along-dip)
 w_f = 15e3;
 dip_rad = dip*Pi/180.;
 
@@ -18,7 +18,7 @@ dip_rad = dip*Pi/180.;
 // Hypocenter, z_n is here the depth along-dip
 x_n = 0e3;
 z_n = -12e3;
-// size of the nucleaton patch
+// Dimensions of the nucleation patch
 l_n = 3e3;
 w_n = 3e3;
 
@@ -32,7 +32,7 @@ Z0 = -42e3;
 // Create the domain as a box
 domain = newv; Box(domain) = {X0, Y0, Z0, X1-X0, Y1-Y0, -Z0};
 
-// Create the fault as a vertically dipping rectangle, centered in x at the hypocenter
+// Create the fault as a rectangle in the x-y plane, centered in x at the hypocenter
 fault = news; Rectangle(fault) = {-l_f/2, -w_f, 0, l_f, w_f};
 
 // Create the nucleation patch as a smaller rectangle
@@ -44,7 +44,7 @@ Rotate{ {1, 0, 0}, {0, 0, 0}, dip_rad } { Surface{fault, nucl}; }
 // Intersect the domain box with the fault rectangle at the free surface
 v() = BooleanFragments{ Volume{domain}; Delete; }{ Surface{fault, nucl}; Delete; };
 
-// Update all coordinates that define important surfaces within the mesh
+// Extract surfaces within the mesh for defining boundary conditions
 eps = 1e-3;
 fault_final[] = Surface In BoundingBox{-l_f/2-eps, -w_f-eps, -w_f-eps, l_f/2+eps, w_f+eps, w_f+eps};
 top[] = Surface In BoundingBox{X0-eps, Y0-eps, -eps, X1+eps, Y1+eps, eps};
@@ -66,4 +66,3 @@ Physical Surface(3) = {fault_final[]};
 Physical Surface(5) = {other[]};
 
 Physical Volume(1) = {domain};
-Mesh.MshFileVersion = 2.2;