diff --git a/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py b/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
index 2815053412b3631d20225a57061e945562eec684..1a5e7be16f24a0607e85d2ea2479f22522c8776b 100755
--- a/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
+++ b/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
@@ -402,18 +402,27 @@ if __name__ == '__main__':
         for mesh_resolution, solver_tol in resolutions.items():
             simulation_parameter.update({"solver_tol": solver_tol})
             hlp.info(simulation_parameter["use_case"])
+            processQueue = mp.Queue()
             LDDsim = mp.Process(
                         target=hlp.run_simulation,
                         args=(
                             simulation_parameter,
+                            processQueue,
                             starttime,
                             mesh_resolution
                             )
                         )
             LDDsim.start()
-            LDDsim.join()
+            # LDDsim.join()
             # hlp.run_simulation(
             #     mesh_resolution=mesh_resolution,
             #     starttime=starttime,
             #     parameter=simulation_parameter
             #     )
+
+        LDDsim.join()
+        if mesh_study:
+            simulation_output_dir = processQueue.get()
+            hlp.merge_spacetime_errornorms(isRichards=isRichards,
+                                           resolutions=resolutions,
+                                           output_dir=simulation_output_dir)
diff --git a/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py b/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py
index e10649a0c6c98396dcc23191145741b67524ec24..3861275a01f4d6091ce369980bd3c7808253127f 100755
--- a/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py
+++ b/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py
@@ -417,18 +417,27 @@ if __name__ == '__main__':
         for mesh_resolution, solver_tol in resolutions.items():
             simulation_parameter.update({"solver_tol": solver_tol})
             hlp.info(simulation_parameter["use_case"])
+            processQueue = mp.Queue()
             LDDsim = mp.Process(
                         target=hlp.run_simulation,
                         args=(
                             simulation_parameter,
+                            processQueue,
                             starttime,
                             mesh_resolution
                             )
                         )
             LDDsim.start()
-            LDDsim.join()
+            # LDDsim.join()
             # hlp.run_simulation(
             #     mesh_resolution=mesh_resolution,
             #     starttime=starttime,
             #     parameter=simulation_parameter
             #     )
+
+        LDDsim.join()
+        if mesh_study:
+            simulation_output_dir = processQueue.get()
+            hlp.merge_spacetime_errornorms(isRichards=isRichards,
+                                           resolutions=resolutions,
+                                           output_dir=simulation_output_dir)
diff --git a/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py b/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py
index a47f33fbf12352ee5e90ff67bb740f7c7b65339a..4d1e0bb5ec69b15e31fb327c90e0cc4aa42dd2ef 100755
--- a/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py
+++ b/Usecases/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py
@@ -361,7 +361,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -411,18 +410,27 @@ if __name__ == '__main__':
         for mesh_resolution, solver_tol in resolutions.items():
             simulation_parameter.update({"solver_tol": solver_tol})
             hlp.info(simulation_parameter["use_case"])
+            processQueue = mp.Queue()
             LDDsim = mp.Process(
                         target=hlp.run_simulation,
                         args=(
                             simulation_parameter,
+                            processQueue,
                             starttime,
                             mesh_resolution
                             )
                         )
             LDDsim.start()
-            LDDsim.join()
+            # LDDsim.join()
             # hlp.run_simulation(
             #     mesh_resolution=mesh_resolution,
             #     starttime=starttime,
             #     parameter=simulation_parameter
             #     )
+
+        LDDsim.join()
+        if mesh_study:
+            simulation_output_dir = processQueue.get()
+            hlp.merge_spacetime_errornorms(isRichards=isRichards,
+                                           resolutions=resolutions,
+                                           output_dir=simulation_output_dir)