diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
index 7b96ff47ebfc5f148caa738f9bd8a771df071986..9134d00955f127d744028450f7a8a4f86d487845 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
@@ -55,7 +55,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.00001
 number_of_timesteps = 5
 
@@ -278,7 +278,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -321,22 +320,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
index 2b9d2cd8c4f2a30e39d2872a2aba5b2f7e6ccc25..689b8f01cecce7fa7c5bb78596d17332cba8c473 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
@@ -55,7 +55,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.005
 number_of_timesteps = 5
 
@@ -412,7 +412,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -455,22 +454,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py
index fca9fe24f8fe9047647cd5abd6d052ef68ad0476..352e07f6eb631b78b11deae1e8435ec3590f2c37 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py
@@ -56,7 +56,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.001
 number_of_timesteps = 800
 
@@ -267,7 +267,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -310,22 +309,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py
index a2929e89bc6d38d9ea8d257936d208a3d7a804f5..9152aeed5743310ca20320809e81d0cf8a95aead 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py
@@ -56,7 +56,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.001
 number_of_timesteps = 800
 
@@ -267,7 +267,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -310,22 +309,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py
index 1b40691474c2305a209859a9a8288b6f227bf446..dcc14cc30e4508e226aca3bce0568c03979f2ae4 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py
@@ -56,7 +56,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.01
 number_of_timesteps = 8
 
@@ -267,7 +267,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -310,22 +309,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py
index 6b16a1c022ef2f9256a99c46912a7a36f6097a4e..74776902a9f6d49ed6340267e32ccc117d3543d7 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py
@@ -56,7 +56,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.001
 number_of_timesteps = 5
 
@@ -267,7 +267,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -310,22 +309,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index 09106bdcd97cc86a865d580128ff3ae2bdc0ffd1..84ff68c86fe762b731ce2e37a3c580554eeddfe0 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -57,7 +57,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.001
 number_of_timesteps = 20
 
@@ -416,7 +416,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -442,6 +441,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
@@ -458,22 +458,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py
index d6bc8fe1f6d45036e486fe85f2ecba960f562bee..4d05dc91e07a095fda3d19616032d6e2bb01ed24 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py
@@ -56,7 +56,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.001
 number_of_timesteps = 800
 
@@ -267,7 +267,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -310,22 +309,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py
index c613e0bfc7add32c4b989f2f1d18b51365f88255..2cff89a78f76ad195c6be5314897a09416a0febe 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py
@@ -56,7 +56,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.001
 number_of_timesteps = 800
 
@@ -267,7 +267,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -310,22 +309,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
index 73de79ca9cabfbfa75f9e54a65a0db7d400e179c..fa7034891e8da625c7a4b45bc6205389d3a3764f 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
@@ -315,18 +315,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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py
index 877dae6b28bb540247178a59fa36adcd563fe8a7..ca31d632887efa37c5b086da8e37518c6abc8af8 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py
@@ -56,7 +56,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.001
 number_of_timesteps = 800
 
@@ -267,7 +267,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -310,22 +309,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py
index 249ab3f371d28ad5ce167e4d17b9e4e478541583..f124ebf490a6377e5a0020ca6436d948cf719d71 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py
@@ -56,7 +56,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.001
 number_of_timesteps = 800
 
@@ -267,7 +267,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -310,22 +309,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py
index 54910805fab46bd492de45807caa0bb5c6c47156..6943cfe34cbb32e8ce1f56c00255fdf8da8903c5 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py
@@ -56,7 +56,7 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = {0: 0.0}
 timestep_size = 0.001
 number_of_timesteps = 800
 
@@ -267,7 +267,6 @@ f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
-
 # MAIN ########################################################################
 if __name__ == '__main__':
     # dictionary of simualation parameters to pass to the run function.
@@ -310,22 +309,34 @@ if __name__ == '__main__':
         "write_to_file": write_to_file,
         "analyse_condition": analyse_condition
     }
-    for starttime in starttimes:
+    for number_shift, starttime in starttimes.items():
+        simulation_parameter.update(
+            {"starttime_timestep_number_shift": number_shift}
+        )
         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)