diff --git a/Tests/Data/PhaseField/Kregime_Propagating_jupyter_notebook/Kregime_Propagating_jupyter.ipynb b/Tests/Data/PhaseField/Kregime_Propagating_jupyter_notebook/Kregime_Propagating_jupyter.ipynb
index 1ff1f2395d12ff2f995fd1d9d83f53d30aedeb2c..d53861a1a65e4c70bec8fc1530a38fcbf07740cc 100644
--- a/Tests/Data/PhaseField/Kregime_Propagating_jupyter_notebook/Kregime_Propagating_jupyter.ipynb
+++ b/Tests/Data/PhaseField/Kregime_Propagating_jupyter_notebook/Kregime_Propagating_jupyter.ipynb
@@ -185,8 +185,11 @@
     "prj_name = \"Kregime_Propagating.prj\"\n",
     "meshname = \"mesh_full_pf\"\n",
     "\n",
-    "out_dir = os.environ.get(\"OGS_TESTRUNNER_OUT_DIR\", \"_out\")\n",
-    "os.makedirs(out_dir, exist_ok=True)"
+    "from pathlib import Path\n",
+    "\n",
+    "out_dir = Path(os.environ.get(\"OGS_TESTRUNNER_OUT_DIR\", \"_out\"))\n",
+    "if not out_dir.exists():\n",
+    "    out_dir.mkdir(parents=True)"
    ]
   },
   {
@@ -582,7 +585,7 @@
     "fluid_volume = []\n",
     "surface_energy = []\n",
     "# Open the file for reading\n",
-    "with open(f\"{out_dir}/log.txt\") as fd:\n",
+    "with (out_dir / \"log.txt\").open() as fd:\n",
     "    # Iterate over the lines\n",
     "    for _i, line in enumerate(fd):\n",
     "        match_surface_energy = re.search(\n",
@@ -610,7 +613,7 @@
     "fluid_volume = []\n",
     "pressure = []\n",
     "# Open the file for reading\n",
-    "with open(f\"{out_dir}/log.txt\") as fd:\n",
+    "with (out_dir / \"log.txt\").open() as fd:\n",
     "    # Iterate over the lines\n",
     "    for _i, line in enumerate(fd):\n",
     "        match_pressure = re.search(\n",
diff --git a/Tests/Data/PhaseField/PForthotropy_jupyter_notebook/sen_shear.ipynb b/Tests/Data/PhaseField/PForthotropy_jupyter_notebook/sen_shear.ipynb
index 3e56bdfb07554062abedd5b347909dc9b42ed2ec..06271057ae121df064f7e09ac070ae4b6b87dabc 100644
--- a/Tests/Data/PhaseField/PForthotropy_jupyter_notebook/sen_shear.ipynb
+++ b/Tests/Data/PhaseField/PForthotropy_jupyter_notebook/sen_shear.ipynb
@@ -147,7 +147,7 @@
     "    print(\n",
     "        f\"> Running single edge notched shear test {phasefield_model} - {energy_split_model} ... <\"\n",
     "    )\n",
-    "    logfile = f\"{out_dir}/log_{phasefield_model}_{energy_split_model}.txt\"\n",
+    "    logfile = f\"{out_dir}/log_{phasefield_model}_{energy_split_model}.txt\"  # noqa: F841\n",
     "    model = ogs.OGS(INPUT_FILE=prj_name, PROJECT_FILE=f\"{out_dir}/{prj_name}\", MKL=True)\n",
     "\n",
     "    # generate prefix from properties\n",
@@ -206,7 +206,7 @@
     "    # run ogs\n",
     "    t0 = time.time()\n",
     "    if MPI:\n",
-    "        print(\"  > OGS started execution with MPI - \" f\"{ncores} cores...\")\n",
+    "        print(f\"  > OGS started execution with MPI - {ncores} cores...\")\n",
     "        ! mpirun -np {ncores} ogs {out_dir}/{prj_name} -o {output_dir} >> {logfile}\n",
     "    else:\n",
     "        print(\"  > OGS started execution - \")\n",
@@ -309,8 +309,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import pyvista as pv\n",
-    "\n",
     "reader = pv.get_reader(f\"{out_dir}/AT2_OrthoMasonry.pvd\")\n",
     "\n",
     "plotter = pv.Plotter()\n",
@@ -432,7 +430,7 @@
     "\n",
     "\n",
     "def force_displ_from_pvd(pvd):\n",
-    "    doc = minidom.parse(pvd)\n",
+    "    doc = minidom.parse(str(pvd))\n",
     "    DataSets = doc.getElementsByTagName(\"DataSet\")\n",
     "    vtu_files = [x.getAttribute(\"file\") for x in DataSets]\n",
     "    forces_sum = [force_midpoint(f\"{out_dir}/{x}\") for x in vtu_files]\n",
@@ -457,8 +455,8 @@
     "plt.rc(\"text\", usetex=True)\n",
     "fig.set_size_inches(18.5, 10.5)\n",
     "for i, pre in enumerate(prefixes):\n",
-    "    pvd = f\"{out_dir}/{pre}.pvd\"\n",
-    "    if os.path.isfile(pvd):\n",
+    "    pvd = out_dir / f\"{pre}.pvd\"\n",
+    "    if pvd.is_file():\n",
     "        curve = force_displ_from_pvd(pvd)\n",
     "        ax.plot(\n",
     "            curve[0],\n",
diff --git a/Tests/Data/PhaseField/beam_jupyter_notebook/beam.ipynb b/Tests/Data/PhaseField/beam_jupyter_notebook/beam.ipynb
index 1a201beecd618d6e33e0ecedefc3d132e60877a2..0f4bb75b9c77782a41d81c0fa491a6bcc569445b 100644
--- a/Tests/Data/PhaseField/beam_jupyter_notebook/beam.ipynb
+++ b/Tests/Data/PhaseField/beam_jupyter_notebook/beam.ipynb
@@ -150,14 +150,14 @@
     "    with_hypre = \"-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_strong_threshold 0.7 -ksp_atol 1e-8 -ksp_rtol 1e-8\"\n",
     "    # file's name\n",
     "    prj_name = \"beam.prj\"\n",
-    "    print(f\"> Running beam model {phasefield_model} – {energy_split_model} ... <\")\n",
-    "    logfile = f\"{out_dir}/log_{phasefield_model}_{energy_split_model}.txt\"\n",
+    "    print(f\"> Running beam model {phasefield_model} - {energy_split_model} ... <\")\n",
+    "    logfile = f\"{out_dir}/log_{phasefield_model}_{energy_split_model}.txt\"  # noqa: F841\n",
     "    # beam dimensions\n",
     "    beam_height = 0.05\n",
-    "    beam_depth = beam_height\n",
-    "    beam_length = 1.0\n",
+    "    beam_depth = beam_height  # noqa: F841\n",
+    "    beam_length = 1.0  # noqa: F841\n",
     "    # mesh properties\n",
-    "    h = mesh_size  # distance between nodes\n",
+    "    h = mesh_size  # noqa: F841, distance between nodes\n",
     "    ls = length_scale\n",
     "    # generate prefix from properties\n",
     "    if energy_split_model == \"VolumetricDeviatoric\":\n",
@@ -365,7 +365,7 @@
     "\n",
     "\n",
     "def force_displ_from_pvd(pvd):\n",
-    "    doc = minidom.parse(pvd)\n",
+    "    doc = minidom.parse(str(pvd))\n",
     "    DataSets = doc.getElementsByTagName(\"DataSet\")\n",
     "    vtu_files = [x.getAttribute(\"file\") for x in DataSets]\n",
     "    forces_right = [force_right(f\"{out_dir}/{x}\") for x in vtu_files]\n",
@@ -424,8 +424,8 @@
     "plt.rc(\"text\", usetex=True)\n",
     "fig.set_size_inches(18.5, 10.5)\n",
     "for i, pre in enumerate(prefixes):\n",
-    "    pvd = f\"{output_dir}/{pre}.pvd\"\n",
-    "    if os.path.isfile(pvd):\n",
+    "    pvd = output_dir / f\"{pre}.pvd\"\n",
+    "    if pvd.is_file():\n",
     "        curve = force_displ_from_pvd(pvd)\n",
     "        ax.plot(\n",
     "            curve[1],\n",
diff --git a/Tests/Data/PhaseField/surfing/Surfing_python.py b/Tests/Data/PhaseField/surfing/Surfing_python.py
index e99c3c562ad0f897f14d8aba969d05721c8f0854..a8e240266813b4e772cc83c7ee8258649deaa398 100644
--- a/Tests/Data/PhaseField/surfing/Surfing_python.py
+++ b/Tests/Data/PhaseField/surfing/Surfing_python.py
@@ -55,7 +55,7 @@ def AsymptoticSolution_Y(x, y, t):
 
 # Dirichlet BCs X-dir
 class BCTop_X(OpenGeoSys.BoundaryCondition):
-    def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+    def getDirichletBCValue(self, t, coords, _node_id, _primary_vars):
         x, y, z = coords
         assert y == b
         value = AsymptoticSolution_X(x, y, t)
@@ -63,7 +63,7 @@ class BCTop_X(OpenGeoSys.BoundaryCondition):
 
 
 class BCLeft_X(OpenGeoSys.BoundaryCondition):
-    def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+    def getDirichletBCValue(self, t, coords, _node_id, _primary_vars):
         x, y, z = coords
         assert x == 0.0
         value = AsymptoticSolution_X(x, y, t)
@@ -71,7 +71,7 @@ class BCLeft_X(OpenGeoSys.BoundaryCondition):
 
 
 class BCBottom_X(OpenGeoSys.BoundaryCondition):
-    def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+    def getDirichletBCValue(self, t, coords, _node_id, _primary_vars):
         x, y, z = coords
         assert y == 0.0
         value = AsymptoticSolution_X(x, y, t)
@@ -79,7 +79,7 @@ class BCBottom_X(OpenGeoSys.BoundaryCondition):
 
 
 class BCRight_X(OpenGeoSys.BoundaryCondition):
-    def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+    def getDirichletBCValue(self, t, coords, _node_id, _primary_vars):
         x, y, z = coords
         assert x == a
         value = AsymptoticSolution_X(x, y, t)
@@ -88,7 +88,7 @@ class BCRight_X(OpenGeoSys.BoundaryCondition):
 
 # Dirichlet BCs Y-dir
 class BCTop_Y(OpenGeoSys.BoundaryCondition):
-    def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+    def getDirichletBCValue(self, t, coords, _node_id, _primary_vars):
         x, y, z = coords
         assert y == b
         value = AsymptoticSolution_Y(x, y, t)
@@ -96,7 +96,7 @@ class BCTop_Y(OpenGeoSys.BoundaryCondition):
 
 
 class BCLeft_Y(OpenGeoSys.BoundaryCondition):
-    def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+    def getDirichletBCValue(self, t, coords, _node_id, _primary_vars):
         x, y, z = coords
         assert x == 0.0
         value = AsymptoticSolution_Y(x, y, t)
@@ -104,7 +104,7 @@ class BCLeft_Y(OpenGeoSys.BoundaryCondition):
 
 
 class BCBottom_Y(OpenGeoSys.BoundaryCondition):
-    def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+    def getDirichletBCValue(self, t, coords, _node_id, _primary_vars):
         x, y, z = coords
         assert y == 0.0
         value = AsymptoticSolution_Y(x, y, t)
@@ -112,7 +112,7 @@ class BCBottom_Y(OpenGeoSys.BoundaryCondition):
 
 
 class BCRight_Y(OpenGeoSys.BoundaryCondition):
-    def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+    def getDirichletBCValue(self, t, coords, _node_id, _primary_vars):
         x, y, z = coords
         assert x == a
         value = AsymptoticSolution_Y(x, y, t)