diff --git a/src/View/Results/CustomPlot/Plot.py b/src/View/Results/CustomPlot/Plot.py
index cea91e0d8b2edf12978d9ac1ace818ecf8c14e99..9f3a005ed967549a15598544f94a7f94d8d15e64 100644
--- a/src/View/Results/CustomPlot/Plot.py
+++ b/src/View/Results/CustomPlot/Plot.py
@@ -215,13 +215,34 @@ class CustomPlot(PamhyrPlot):
         if "bed_elevation" in self._y:
 
             ax = self._axes[unit["bed_elevation"]]
-            line = ax.plot(
-                rk, z_min,
-                color='grey', lw=1.,
-            )
+            if self._current_res_id != 2:
+                line = ax.plot(
+                    rk, z_min,
+                    color='grey', lw=1.,
+                )
+            else:
+                if reach.has_sediment():
+                    z_min1 = self.draw_bottom_with_bedload(reach1)
+                    z_min2 = self.draw_bottom_with_bedload(reach2)
+                else:
+                    z_min1 = reach1.geometry.get_z_min()
+                    z_min2 = reach2.geometry.get_z_min()
+                dz = list(
+                    map(
+                        lambda x, y: x - y,
+                        z_min2, z_min2
+                        )
+                    )
+                line = ax.plot(
+                    rk, dz,
+                    color='grey', lw=1.,
+                )
+
             self.lines["bed_elevation"] = line
 
-            if self._envelop and reach.has_sediment():
+            if (self._envelop and
+                reach.has_sediment() and
+                self._current_res_id != 2):
 
                 ax = self._axes[unit["bed_elevation_envelop"]]
 
@@ -619,7 +640,22 @@ class CustomPlot(PamhyrPlot):
             )
         )
         if "bed_elevation" in self._y:
-            self.lines["bed_elevation"][0].set_ydata(z_min)
+            if self._current_res_id != 2:
+                self.lines["bed_elevation"][0].set_ydata(z_min)
+            else:
+                if reach.has_sediment():
+                    z_min1 = self.draw_bottom_with_bedload(reach1)
+                    z_min2 = self.draw_bottom_with_bedload(reach2)
+                else:
+                    z_min1 = reach1.geometry.get_z_min()
+                    z_min2 = reach2.geometry.get_z_min()
+                dz = list(
+                    map(
+                        lambda x, y: x - y,
+                        z_min2, z_min2
+                        )
+                    )
+            self.lines["bed_elevation"][0].set_ydata(dz)
 
         if "water_elevation" in self._y:
             self.lines["water_elevation"][0].set_ydata(z)
@@ -827,6 +863,8 @@ class CustomPlot(PamhyrPlot):
             )
 
         self.canvas.axes.set_xticks(ticks=fx, labels=xt, rotation=45)
+        self.canvas.axes.relim(visible_only=True)
+        self.canvas.axes.autoscale_view()
 
     def _draw_time(self):
         results = self.data[self._current_res_id]
@@ -863,15 +901,16 @@ class CustomPlot(PamhyrPlot):
         z = profile.get_key("Z")
         v = profile.get_key("V")
         z_min = profile.geometry.z_min()
-        if reach.has_sediment():
-            ts_z_min = self.get_ts_zmin(self._profile)
-        else:
-            ts_z_min = list(
-                map(
-                    lambda ts:  z_min,
-                    ts
+        if self._current_res_id != 2:
+            if reach.has_sediment():
+                ts_z_min = self.get_ts_zmin(self._profile)
+            else:
+                ts_z_min = list(
+                    map(
+                        lambda ts:  z_min,
+                        ts
+                    )
                 )
-            )
 
         self.lines = {}
         if "bed_elevation" in self._y:
@@ -879,6 +918,26 @@ class CustomPlot(PamhyrPlot):
 
             ax = self._axes[unit["bed_elevation"]]
 
+            if self._current_res_id == 2:
+                if reach.has_sediment():
+                    ts_z_min1 = self.get_ts_zmin(self._profile1)
+                    ts_z_min2 = self.get_ts_zmin(self._profile2)
+                    ts_z_min = list(
+                        map(
+                            lambda x, y: x - y,
+                            ts_z_min1, ts_z_min2
+                            )
+                        )
+                else:
+                    z_min1 = profile1.geometry.z_min()
+                    z_min2 = profile2.geometry.z_min()
+                    ts_z_min = list(
+                        map(
+                            lambda ts:  z_min1 - z_min2,
+                            ts
+                        )
+                    )
+
             line = ax.plot(
                 ts, ts_z_min,
                 color='grey', lw=1.
@@ -991,16 +1050,16 @@ class CustomPlot(PamhyrPlot):
                     map(lambda z, v:
                         v /
                         sqrt(9.81 * (
-                            profile.geometry.wet_area(z) /
-                            profile.geometry.wet_width(z))
+                            profile1.geometry.wet_area(z) /
+                            profile1.geometry.wet_width(z))
                         ), z1, v1)
                     )
                 d2 = list(
                     map(lambda z, v:
                         v /
                         sqrt(9.81 * (
-                            profile.geometry.wet_area(z) /
-                            profile.geometry.wet_width(z))
+                            profile2.geometry.wet_area(z) /
+                            profile2.geometry.wet_width(z))
                         ), z2, v2)
                     )
                 d = list(
@@ -1024,10 +1083,10 @@ class CustomPlot(PamhyrPlot):
                 )
             else:
                 d1 = list(
-                    map(lambda z: profile.geometry.wet_area(z), z1)
+                    map(lambda z: profile1.geometry.wet_area(z), z1)
                 )
                 d2 = list(
-                    map(lambda z: profile.geometry.wet_area(z), z2)
+                    map(lambda z: profile2.geometry.wet_area(z), z2)
                 )
                 d = list(
                     map(
@@ -1077,17 +1136,37 @@ class CustomPlot(PamhyrPlot):
         q = profile.get_key("Q")
         z = profile.get_key("Z")
         v = profile.get_key("V")
-        if reach.has_sediment():
-            ts_z_min = self.get_ts_zmin(self._profile)
-        else:
-            z_min = profile.geometry.z_min()
-            ts_z_min = list(
-                map(
-                    lambda ts:  z_min,
-                    ts
+        if self._current_res_id != 2:
+            if reach.has_sediment():
+                ts_z_min = self.get_ts_zmin(self._profile)
+            else:
+                z_min = profile.geometry.z_min()
+                ts_z_min = list(
+                    map(
+                        lambda ts:  z_min,
+                        ts
+                    )
                 )
-            )
         if "bed_elevation" in self._y:
+            if self._current_res_id == 2:
+                if reach.has_sediment():
+                    ts_z_min1 = self.get_ts_zmin(self._profile1)
+                    ts_z_min2 = self.get_ts_zmin(self._profile2)
+                    ts_z_min = list(
+                        map(
+                            lambda x, y: x - y,
+                            ts_z_min1, ts_z_min2
+                            )
+                        )
+                else:
+                    z_min1 = profile1.geometry.z_min()
+                    z_min2 = profile2.geometry.z_min()
+                    ts_z_min = list(
+                        map(
+                            lambda ts:  z_min1 - z_min2,
+                            ts
+                        )
+                    )
             self.lines["bed_elevation"][0].set_ydata(ts_z_min)
 
         if "water_elevation" in self._y:
@@ -1164,16 +1243,16 @@ class CustomPlot(PamhyrPlot):
                     map(lambda z, v:
                         v /
                         sqrt(9.81 * (
-                            profile.geometry.wet_area(z) /
-                            profile.geometry.wet_width(z))
+                            profile1.geometry.wet_area(z) /
+                            profile1.geometry.wet_width(z))
                         ), z1, v1)
                     )
                 d2 = list(
                     map(lambda z, v:
                         v /
                         sqrt(9.81 * (
-                            profile.geometry.wet_area(z) /
-                            profile.geometry.wet_width(z))
+                            profile2.geometry.wet_area(z) /
+                            profile2.geometry.wet_width(z))
                         ), z2, v2)
                     )
                 d = list(
@@ -1192,10 +1271,10 @@ class CustomPlot(PamhyrPlot):
                 )
             else:
                 d1 = list(
-                    map(lambda z: profile.geometry.wet_area(z), z1)
+                    map(lambda z: profile1.geometry.wet_area(z), z1)
                 )
                 d2 = list(
-                    map(lambda z: profile.geometry.wet_area(z), z2)
+                    map(lambda z: profile2.geometry.wet_area(z), z2)
                 )
                 d = list(
                     map(
@@ -1271,6 +1350,7 @@ class CustomPlot(PamhyrPlot):
                 self._axes[axes].set_ylabel(
                     "Δ " + self._trad[axes],
                     color='black', fontsize=10
+                )
 
         if self._x == "rk":
             self._draw_rk()
diff --git a/src/View/Results/Window.py b/src/View/Results/Window.py
index f01313a9affc2662bc988997283762253dc0e29d..0fa5f8283eac586967a45315e2def9b386f9c909 100644
--- a/src/View/Results/Window.py
+++ b/src/View/Results/Window.py
@@ -334,6 +334,9 @@ class ResultsWindow(PamhyrWindow):
             # "action_export": self.export_current,
         }
 
+        if len(self._results) > 1:
+            self.find(QAction, "action_reload").setEnabled(False)
+
         for action in actions:
             self.find(QAction, action).triggered.connect(
                 actions[action]
@@ -533,12 +536,13 @@ class ResultsWindow(PamhyrWindow):
 
     def _reload(self):
         logger.debug("Reload results...")
-        self._results = self._results.reload()
+        if len(self._results) == 1:
+            self._results[0] = self._results[0].reload()
 
-        self._timestamps = sorted(list(self._results.get("timestamps")))
+            self.get_timestamps()
 
-        self._reload_plots()
-        self._reload_slider()
+            self._reload_plots()
+            self._reload_slider()
 
     def _add_custom_plot(self):
         dlg = CustomPlotValuesSelectionDialog(parent=self)
@@ -680,9 +684,9 @@ class ResultsWindow(PamhyrWindow):
         )
 
     def export_to(self, filename, x, y, envelop):
-        timestamps = sorted(self._results.get("timestamps"))
-        reach = self._results.river.reachs[self._get_current_reach()]
-        first_line = [f"Study: {self._results.study.name}",
+        results = self._results[self._current_results]
+        reach = results.river.reachs[self._get_current_reach()]
+        first_line = [f"Study: {results.study.name}",
                       f"Reach: {reach.name}"]
         if x == "rk":
             timestamp = self._get_current_timestamp()
@@ -740,7 +744,8 @@ class ResultsWindow(PamhyrWindow):
         )
 
     def export_current_to(self, directory):
-        reach = self._results.river.reachs[self._get_current_reach()]
+        results = self._results[self._current_results]
+        reach = results.river.reachs[self._get_current_reach()]
         self.export_all(reach, directory, [self._get_current_timestamp()])
 
     def delete_tab(self, index):
@@ -749,15 +754,31 @@ class ResultsWindow(PamhyrWindow):
         tab_widget.removeTab(index)
 
     def _export_rk(self, timestamp, y, envelop):
-        reach = self._results.river.reachs[self._get_current_reach()]
+        results = self._results[self._current_results]
+        reach = results.river.reachs[self._get_current_reach()]
         dict_x = self._trad.get_dict("values_x")
         dict_y = self._trad.get_dict("values_y")
+        if self._current_results == 2:
+            envelop = False
+            reach1 = results[0].river.reachs[self._get_current_reach()]
+            reach2 = results[1].river.reachs[self._get_current_reach()]
         if envelop:
             dict_y.update(self._trad.get_dict("values_y_envelop"))
         my_dict = {}
         my_dict[dict_x["rk"]] = reach.geometry.get_rk()
         if "bed_elevation" in y:
-            my_dict[dict_y["bed_elevation"]] = reach.geometry.get_z_min()
+            if self._current_results != 2:
+                zmin = reach.geometry.get_z_min()
+            else:
+                z_min1 = reach1.geometry.get_z_min()
+                z_min2 = reach2.geometry.get_z_min()
+                zmin = list(
+                    map(
+                        lambda x, y: x - y,
+                        z_min2, z_min2
+                        )
+                    )
+            my_dict[dict_y["bed_elevation"]] = zmin
             # if envelop and reach.has_sediment():
         if "discharge" in y:
             my_dict[dict_y["discharge"]] = list(
@@ -817,13 +838,36 @@ class ResultsWindow(PamhyrWindow):
                 my_dict[dict_y["max_velocity"]] = [max(v) for v in velocities]
 
         if "depth" in y:
-            my_dict[dict_y["depth"]] = list(
-                map(
-                    lambda p: p.geometry.max_water_depth(
-                        p.get_ts_key(timestamp, "Z")),
-                    reach.profiles
+            if self._current_results != 2:
+                d = my_dict[dict_y["depth"]] = list(
+                    map(
+                        lambda p: p.geometry.max_water_depth(
+                            p.get_ts_key(timestamp, "Z")),
+                        reach.profiles
+                    )
                 )
-            )
+            else:
+                d1 = list(
+                    map(
+                        lambda p: p.geometry.max_water_depth(
+                            p.get_ts_key(timestamp, "Z")),
+                        reach1.profiles
+                    )
+                )
+                d2 = list(
+                    map(
+                        lambda p: p.geometry.max_water_depth(
+                            p.get_ts_key(timestamp, "Z")),
+                        reach2.profiles
+                    )
+                )
+                d = list(
+                    map(
+                        lambda x, y: x - y,
+                        d1, d2
+                        )
+                    )
+            my_dict[dict_y["depth"]] = d
             if envelop:
                 my_dict[dict_y["min_depth"]] = list(map(
                     lambda p1, p2: p1 - p2, map(
@@ -841,53 +885,157 @@ class ResultsWindow(PamhyrWindow):
                 )
 
         if "mean_depth" in y:
-            my_dict[dict_y["mean_depth"]] = list(
-                map(
-                    lambda p: p.geometry.mean_water_depth(
-                        p.get_ts_key(timestamp, "Z")),
-                    reach.profiles
+            if self._current_results != 2:
+                d = list(
+                    map(
+                        lambda p: p.geometry.mean_water_depth(
+                            p.get_ts_key(timestamp, "Z")),
+                        reach.profiles
+                    )
                 )
-            )
+            else:
+                d1 = list(
+                    map(
+                        lambda p: p.geometry.mean_water_depth(
+                            p.get_ts_key(timestamp, "Z")),
+                        reach1.profiles
+                    )
+                )
+                d2 = list(
+                    map(
+                        lambda p: p.geometry.mean_water_depth(
+                            p.get_ts_key(timestamp, "Z")),
+                        reach2.profiles
+                    )
+                )
+                d = list(
+                    map(
+                        lambda x, y: x - y,
+                        d1, d2
+                        )
+                    )
+            my_dict[dict_y["mean_depth"]] = d
         if "froude" in y:
-            my_dict[dict_y["froude"]] = list(
-                map(
-                    lambda p:
-                        p.get_ts_key(timestamp, "V") /
-                        sqrt(9.81 * (
-                            p.geometry.wet_area(
-                                p.get_ts_key(timestamp, "Z")) /
-                            p.geometry.wet_width(
-                                p.get_ts_key(timestamp, "Z"))
-                            )),
-                    reach.profiles
+            if self._current_results != 2:
+                fr = my_dict[dict_y["froude"]] = list(
+                    map(
+                        lambda p:
+                            p.get_ts_key(timestamp, "V") /
+                            sqrt(9.81 * (
+                                p.geometry.wet_area(
+                                    p.get_ts_key(timestamp, "Z")) /
+                                p.geometry.wet_width(
+                                    p.get_ts_key(timestamp, "Z"))
+                                )),
+                        reach.profiles
+                    )
                 )
-            )
+            else:
+                fr1 = list(
+                    map(
+                        lambda p:
+                            p.get_ts_key(timestamp, "V") /
+                            sqrt(9.81 * (
+                                p.geometry.wet_area(
+                                    p.get_ts_key(timestamp, "Z")) /
+                                p.geometry.wet_width(
+                                    p.get_ts_key(timestamp, "Z"))
+                                )),
+                        reach1.profiles
+                    )
+                )
+                fr2 = list(
+                    map(
+                        lambda p:
+                            p.get_ts_key(timestamp, "V") /
+                            sqrt(9.81 * (
+                                p.geometry.wet_area(
+                                    p.get_ts_key(timestamp, "Z")) /
+                                p.geometry.wet_width(
+                                    p.get_ts_key(timestamp, "Z"))
+                                )),
+                        reach2.profiles
+                    )
+                )
+                fr =  list(
+                    map(
+                        lambda x, y: x - y,
+                        fr1, fr2
+                        )
+                    )
+            my_dict[dict_y["froude"]] = fr
         if "wet_area" in y:
-            my_dict[dict_y["wet_area"]] = list(
-                map(
-                    lambda p: p.geometry.wet_area(
-                        p.get_ts_key(timestamp, "Z")),
-                    reach.profiles
+            if self._current_results != 2:
+                wa = list(
+                    map(
+                        lambda p: p.geometry.wet_area(
+                            p.get_ts_key(timestamp, "Z")),
+                        reach.profiles
+                    )
                 )
-            )
+            else:
+                wa1 = list(
+                    map(
+                        lambda p: p.geometry.wet_area(
+                            p.get_ts_key(timestamp, "Z")),
+                        reach1.profiles
+                    )
+                )
+                wa2 = list(
+                    map(
+                        lambda p: p.geometry.wet_area(
+                            p.get_ts_key(timestamp, "Z")),
+                        reach2.profiles
+                    )
+                )
+                wa = list(
+                    map(
+                        lambda x, y: x - y,
+                        wa1, wa2
+                        )
+                    )
+            my_dict[dict_y["wet_area"]] = wa
 
         return my_dict
 
     def _export_time(self, profile, y):
-        reach = self._results.river.reachs[self._get_current_reach()]
+        results = self._results[self._current_results]
+        reach = results.river.reachs[self._get_current_reach()]
         profile = reach.profile(profile)
-        ts = list(self._results.get("timestamps"))
-        ts.sort()
         dict_x = self._trad.get_dict("values_x")
         dict_y = self._trad.get_dict("values_y")
         my_dict = {}
-        my_dict[dict_x["time"]] = ts
+        my_dict[dict_x["time"]] = self._timestamps
         z = profile.get_key("Z")
         q = profile.get_key("Q")
         v = profile.get_key("V")
+        if self._current_results == 2:
+            reach1 = self._results[0].river.reach(self._reach)
+            reach2 = self._results[1].river.reach(self._reach)
+            profile1 = reach1.profile(self._profile)
+            profile2 = reach2.profile(self._profile)
+
+            q1 = profile1.get_key("Q")
+            z1 = profile1.get_key("Z")
+            v1 = profile1.get_key("V")
+
+            q2 = profile2.get_key("Q")
+            z2 = profile2.get_key("Z")
+            v2 = profile2.get_key("V")
+
         if "bed_elevation" in y:
-            my_dict[dict_y["bed_elevation"]] = [
-                profile.geometry.z_min()] * len(ts)
+            if self._current_results != 2:
+                z_min = [profile.geometry.z_min()] * len(self._timestamps)
+            else:
+                    z_min1 = profile1.geometry.z_min()
+                    z_min2 = profile2.geometry.z_min()
+                    z_min = list(
+                        map(
+                            lambda ts:  z_min1 - z_min2,
+                            self._timestamps
+                        )
+                    )
+            my_dict[dict_y["bed_elevation"]] = z_min
         if "discharge" in y:
             my_dict[dict_y["discharge"]] = q
         if "water_elevation" in y:
@@ -895,25 +1043,95 @@ class ResultsWindow(PamhyrWindow):
         if "velocity" in y:
             my_dict[dict_y["velocity"]] = v
         if "depth" in y:
-            my_dict[dict_y["depth"]] = list(
-                map(lambda z: profile.geometry.max_water_depth(z), z)
-            )
+            if self._current_results != 2:
+                d = list(
+                    map(lambda z: profile.geometry.max_water_depth(z), z)
+                )
+            else:
+                d1 = list(
+                    map(lambda z: profile1.geometry.max_water_depth(z), z1)
+                )
+                d2 = list(
+                    map(lambda z: profile2.geometry.max_water_depth(z), z2)
+                )
+                d = list(
+                    map(
+                        lambda x, y: x - y,
+                        d1, d2
+                        )
+                    )
+            my_dict[dict_y["depth"]] = d
         if "mean_depth" in y:
-            my_dict[dict_y["mean_depth"]] = list(
-                map(lambda z: profile.geometry.mean_water_depth(z), z)
-            )
-        if "froude" in y:
-            my_dict[dict_y["froude"]] = list(
-                map(lambda z, v:
-                    v / sqrt(9.81 * (
-                        profile.geometry.wet_area(z) /
-                        profile.geometry.wet_width(z))
-                    ), z, v)
+            if self._current_results != 2:
+                d = list(
+                    map(lambda z: profile.geometry.mean_water_depth(z), z)
                 )
+            else:
+                d1 = list(
+                    map(lambda z: profile1.geometry.mean_water_depth(z), z1)
+                )
+                d2 = list(
+                    map(lambda z: profile2.geometry.mean_water_depth(z), z2)
+                )
+                d = list(
+                    map(
+                        lambda x, y: x - y,
+                        d1, d2
+                        )
+                    )
+            my_dict[dict_y["mean_depth"]] = d
+        if "froude" in y:
+            if self._current_results != 2:
+                fr = my_dict[dict_y["froude"]] = list(
+                    map(lambda z, v:
+                        v / sqrt(9.81 * (
+                            profile.geometry.wet_area(z) /
+                            profile.geometry.wet_width(z))
+                        ), z, v)
+                    )
+            else:
+                fr1 = list(
+                    map(lambda z, v:
+                        v /
+                        sqrt(9.81 * (
+                            profile1.geometry.wet_area(z) /
+                            profile1.geometry.wet_width(z))
+                        ), z1, v1)
+                    )
+                fr2 = list(
+                    map(lambda z, v:
+                        v /
+                        sqrt(9.81 * (
+                            profile2.geometry.wet_area(z) /
+                            profile2.geometry.wet_width(z))
+                        ), z2, v2)
+                    )
+                fr = list(
+                    map(
+                        lambda x, y: x - y,
+                        fr1, fr2
+                        )
+                    )
+            my_dict[dict_y["froude"]] = fr
         if "wet_area" in y:
-            my_dict[dict_y["wet_area"]] = list(
-                map(lambda z: profile.geometry.wet_area(z), z)
-            )
+            if self._current_results != 2:
+                wa = list(
+                    map(lambda z: profile.geometry.wet_area(z), z)
+                )
+            else:
+                wa1 = list(
+                    map(lambda z: profile1.geometry.wet_area(z), z1)
+                )
+                wa2 = list(
+                    map(lambda z: profile2.geometry.wet_area(z), z2)
+                )
+                wa = list(
+                    map(
+                        lambda x, y: x - y,
+                        wa1, wa2
+                        )
+                    )
+            my_dict[dict_y["wet_area"]] = wa
 
         return my_dict