diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
index b991ae3e8a8eb923645e84cd292bcbd108d443d4..405bc4251d86ee640a956a0b052909fcf8742387 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
@@ -59,7 +59,7 @@ resolutions = {
 #  for each element t_0 in starttimes.
 starttimes = [0.0]
 timestep_size = 0.001
-number_of_timesteps = 1
+number_of_timesteps = 5
 
 # LDD scheme parameters  ######################################################
 
@@ -953,9 +953,9 @@ for subdomain in isRichards.keys():
 
 # RUN #########################################################################
 def run_simulation(
-        **kwargs
-        # starttime: float = 0.0,
-        # mesh_resolution: float = 32
+        parameter,
+        starttime: float = 0.0,
+        mesh_resolution: float = 32,
         ) -> None:
     """Setup and run LDD simulation
 
@@ -967,53 +967,51 @@ def run_simulation(
                                     of the simulation.
     mesh_resolution #type float     mesh resolution parameter for Fenics.
     """
-    for key, value in kwargs.items():
-        print ("%s == %s" %(key, value))
-
+    # parameter = kwargs
     # initialise LDD simulation class
     simulation = ldd.LDDsimulation(
-        tol=kwargs["tol"],
-        LDDsolver_tol=kwargs["solver_tol"],
-        debug=kwargs["debugflag"],
-        max_iter_num=kwargs["max_iter_num"],
-        FEM_Lagrange_degree=kwargs["FEM_Lagrange_degree"],
-        mesh_study=kwargs["mesh_study"]
+        tol=parameter["tol"],
+        LDDsolver_tol=parameter["solver_tol"],
+        debug=parameter["debugflag"],
+        max_iter_num=parameter["max_iter_num"],
+        FEM_Lagrange_degree=parameter["FEM_Lagrange_degree"],
+        mesh_study=parameter["mesh_study"]
         )
 
     simulation.set_parameters(
-        use_case=kwargs["use_case"],
-        output_dir=kwargs["output_string"],
-        subdomain_def_points=kwargs["subdomain_def_points"],
-        isRichards=kwargs["isRichards"],
-        interface_def_points=kwargs["interface_def_points"],
-        outer_boundary_def_points=kwargs["outer_boundary_def_points"],
-        adjacent_subdomains=kwargs["adjacent_subdomains"],
-        mesh_resolution=kwargs["mesh_resolution"],
-        viscosity=kwargs["viscosity"],
-        porosity=kwargs["porosity"],
-        L=kwargs["L"],
-        lambda_param=kwargs["lambda_param"],
-        relative_permeability=kwargs["relative_permeability"],
-        saturation=kwargs["sat_pressure_relationship"],
-        starttime=kwargs["starttime"],
-        number_of_timesteps=kwargs["number_of_timesteps"],
-        number_of_timesteps_to_analyse=kwargs["number_of_timesteps_to_analyse"],
-        plot_timestep_every=kwargs["plot_timestep_every"],
-        timestep_size=kwargs["timestep_size"],
-        sources=kwargs["source_expression"],
-        initial_conditions=kwargs["initial_condition"],
-        dirichletBC_expression_strings=kwargs["dirichletBC"],
-        exact_solution=kwargs["exact_solution"],
-        densities=kwargs["densities"],
-        include_gravity=kwargs["include_gravity"],
-        gravity_acceleration=kwargs["gravity_acceleration"],
-        write2file=kwargs["write_to_file"],
+        use_case=parameter["use_case"],
+        output_dir=parameter["output_string"],
+        subdomain_def_points=parameter["subdomain_def_points"],
+        isRichards=parameter["isRichards"],
+        interface_def_points=parameter["interface_def_points"],
+        outer_boundary_def_points=parameter["outer_boundary_def_points"],
+        adjacent_subdomains=parameter["adjacent_subdomains"],
+        mesh_resolution=mesh_resolution,  # parameter["mesh_resolution"],
+        viscosity=parameter["viscosity"],
+        porosity=parameter["porosity"],
+        L=parameter["L"],
+        lambda_param=parameter["lambda_param"],
+        relative_permeability=parameter["relative_permeability"],
+        saturation=parameter["sat_pressure_relationship"],
+        starttime=starttime,  # parameter["starttime"],
+        number_of_timesteps=parameter["number_of_timesteps"],
+        number_of_timesteps_to_analyse=parameter["number_of_timesteps_to_analyse"],
+        plot_timestep_every=parameter["plot_timestep_every"],
+        timestep_size=parameter["timestep_size"],
+        sources=parameter["source_expression"],
+        initial_conditions=parameter["initial_condition"],
+        dirichletBC_expression_strings=parameter["dirichletBC"],
+        exact_solution=parameter["exact_solution"],
+        densities=parameter["densities"],
+        include_gravity=parameter["include_gravity"],
+        gravity_acceleration=parameter["gravity_acceleration"],
+        write2file=parameter["write_to_file"],
         )
 
     simulation.initialise()
     output_dir = simulation.output_dir
     # simulation.write_exact_solution_to_xdmf()
-    output = simulation.run(analyse_condition=kwargs["analyse_condition"])
+    output = simulation.run(analyse_condition=parameter["analyse_condition"])
 
     for subdomain_index, subdomain_output in output.items():
         mesh_h = subdomain_output['mesh_size']
@@ -1068,7 +1066,6 @@ def info(title):
 
 # MAIN ########################################################################
 if __name__ == '__main__':
-    # mp.set_start_method('spawn')
     # dictionary of simualation parameters to pass to the run function.
     # mesh_resolution and starttime are excluded, as they get passed explicitly
     # to achieve parallelisation in these parameters in these parameters for
@@ -1110,24 +1107,22 @@ if __name__ == '__main__':
     }
     for starttime in starttimes:
         for mesh_resolution, solver_tol in resolutions.items():
-            simulation_parameter.update(
-                {"solver_tol": solver_tol,
-                 "mesh_resolution": mesh_resolution,
-                 "starttime": starttime}
-            )
-            # info(simulation_parameter["use_case"])
-            # LDDsim = mp.Process(
-            #             target=run_simulation,
-            #             args=(
-            #                 starttime,
-            #                 mesh_resolution,
-            #                 simulation_parameter
-            #                 )
-            #             )
-            # LDDsim.start()
-            # LDDsim.join()
-            run_simulation(
-                # simulation_parameter,
-                # starttime=starttime,
-                kwargs=mesh_resolution
-                )
+            simulation_parameter.update({"solver_tol": solver_tol})
+            simulation_parameter.update({"mesh_resolution": mesh_resolution})
+            simulation_parameter.update({"starttime": starttime})
+            info(simulation_parameter["use_case"])
+            LDDsim = mp.Process(
+                        target=run_simulation,
+                        args=(
+                            simulation_parameter,
+                            starttime,
+                            mesh_resolution
+                            )
+                        )
+            LDDsim.start()
+            LDDsim.join()
+            # run_simulation(
+            #     mesh_resolution=mesh_resolution,
+            #     starttime=starttime,
+            #     parameter=simulation_parameter
+            #     )