diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 78921ffd83504dd8dc5d4745ca33edcb53d5a0b6..1c678b0722b614ee8a5b44c807c4f505fb7f1cad 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -491,30 +491,6 @@ class LDDsimulation(object):
             # the future.
             for sd_index, subdomain in self.subdomain.items():
                 # set subdomain iteration to current iteration
-                # # mp.set_start_method('fork')
-                # processQueue = mp.Queue()
-                # Lsolver_step = mp.Process(
-                #             target=self.run_Lsolver_step,
-                #             args=(
-                #                 iteration,
-                #                 sd_index,
-                #                 processQueue,
-                #                 solution_over_iteration_within_timestep,
-                #                 debug,
-                #                 analyse_timestep,
-                #                 analyse_condition,
-                #             )
-                #         )
-                # Lsolver_step.start()
-                # max_subsequent_error[sd_index-1] = processQueue.get()
-                # print(Lsolver_step.exitcode)
-                # Lsolver_step.join()
-                # print(Lsolver_step.exitcode)
-                # # update interface iteration numbers
-                # # for index in subdomain.has_interface:
-                # #     interface = self.interface[index]
-                # #     interface.current_iteration[sd_index] = iteration
-                # # Lsolver_step.terminate()
                 max_subsequent_error[sd_index - 1] = self.run_Lsolver_step(
                     iteration=iteration,
                     subdomain_index=sd_index,
diff --git a/LDDsimulation/helpers.py b/LDDsimulation/helpers.py
index 7f3696b47096540efdb966293c79de9da8d3a193..d4c3a70e1529d291e276994958f317afec4de30f 100644
--- a/LDDsimulation/helpers.py
+++ b/LDDsimulation/helpers.py
@@ -13,11 +13,13 @@ import sympy as sym
 import os
 import pandas as pd
 import LDDsimulation as ldd
+import csv
 
 
 # RUN #########################################################################
 def run_simulation(
         parameter,
+        processQueue,
         starttime: float = 0.0,
         mesh_resolution: float = 32,
         ) -> None:
@@ -78,6 +80,7 @@ def run_simulation(
 
     simulation.initialise()
     output_dir = simulation.output_dir
+    processQueue.put(output_dir)
     # simulation.write_exact_solution_to_xdmf()
     output = simulation.run(analyse_condition=parameter["analyse_condition"])
 
@@ -309,3 +312,47 @@ def generate_exact_symbolic_pc(
                 )
 
     return symbolic_pc
+
+
+def merge_spacetime_errornorms(isRichards: tp.Dict[int, bool],
+                               resolutions: tp.Dict[int, float],
+                               output_dir: str):
+    """Merge spacetime error norms for parallel mesh study computation.
+
+    Since spacetime errornorms for each mesh resolution get written into a
+    separate file (to make mesh study computation parallelisable), they need to
+    be merged for ease of plotting.
+    """
+    for subdomain, isR in isRichards.items():
+        subdomain_has_phases = ["wetting"]
+        if not isR:
+            subdomain_has_phases = ["wetting", "nonwetting"]
+
+        for phase in subdomain_has_phases:
+            merge_filename = output_dir + f"subdomain{subdomain}" + \
+                             f"-space-time-errornorm-{phase}-phase.csv"
+
+            with open(merge_filename, 'w', newline='') as mergefile:
+                mergefile_writer = csv.writer(mergefile, delimiter=' ',)
+                # write header into merge five
+                mergefile_writer.writerow(
+                    ['mesh_parameter', 'mesh_h', 'Linf', 'L1', 'L2']
+                    )
+                for meshrez in resolutions.keys():
+                    input_filename = output_dir + f"subdomain{subdomain}" + \
+                                     f"-space-time-errornorm-{phase}-phase" + \
+                                     f"_meshrez{meshrez}.csv"
+                    try:
+                        # print(f"found {input_filename}")
+                        with open(input_filename, 'r', newline='') as inputfile:
+                            inputfile_reader = csv.reader(inputfile,
+                                                          delimiter=' '
+                                                          )
+                            # skip header
+                            next(inputfile_reader)
+                            for row in inputfile_reader:
+                                mergefile_writer.writerow(row)
+                    except IOError:
+                        print("File not accessible or nonexistant")
+
+    print("Merged all spacetime errornorm files")