diff --git a/src/Model/Results/River/River.py b/src/Model/Results/River/River.py
index 51eaf1d19af9a04ef903df2e9e437e5c8245ad3c..7754796dc4c7b9619884f101926991b392c9075b 100644
--- a/src/Model/Results/River/River.py
+++ b/src/Model/Results/River/River.py
@@ -59,6 +59,8 @@ class Profile(object):
     def get_ts_key(self, timestamp, key):
         return self._data[timestamp][key]
 
+    def has_sediment(self):
+        return any(map(lambda ts: "sl" in self._data[ts], self._data))
 
 class Reach(object):
     def __init__(self, reach, study):
@@ -92,6 +94,8 @@ class Reach(object):
     def set(self, profile_id, timestamp, key, data):
         self._profiles[profile_id].set(timestamp, key, data)
 
+    def has_sediment(self):
+        return any(map(lambda profile: profile.has_sediment(), self._profiles))
 
 class River(object):
     def __init__(self, study):
diff --git a/src/View/Results/PlotKPC.py b/src/View/Results/PlotKPC.py
index b54d833ed2af5777a180a996c5c1e1a396ae30bd..c7e288d4b990855651fac4f97931a7518731fd26 100644
--- a/src/View/Results/PlotKPC.py
+++ b/src/View/Results/PlotKPC.py
@@ -18,6 +18,8 @@
 
 import logging
 
+from functools import reduce
+
 from tools import timer
 from View.Tools.PamhyrPlot import PamhyrPlot
 
@@ -86,6 +88,136 @@ class PlotKPC(PamhyrPlot):
         self._init = True
 
     def draw_bottom(self, reach):
+        if reach.has_sediment():
+            self.draw_bottom_with_bedload(reach)
+        else:
+            self.draw_bottom_geometry(reach)
+
+    def draw_bottom_with_bedload(self, reach):
+        kp = reach.geometry.get_kp()
+        z_min = reach.geometry.get_z_min()
+        z_max = reach.geometry.get_z_max()
+
+        initial_sl = self.sl_compute_initial(reach)
+        current_sl = self.sl_compute_current(reach)
+
+        max_sl_num = reduce(
+            lambda acc, sl: max(acc, len(sl)),
+            current_sl, 0
+        )
+
+        sl_init, sl = self.sl_completed_layers(
+            initial_sl, current_sl, max_sl_num
+        )
+
+        z_sl = self.sl_apply_z_min_to_initial(sl_init, z_min)
+        d_sl = self.sl_compute_diff(sl_init, sl)
+
+        final_z_sl = self.sl_apply_diff(reach, z_sl, d_sl)
+        final_z_sl = list(reversed(final_z_sl))
+
+        self.line_kp_sl = []
+        for i, z in enumerate(final_z_sl):
+            self.line_kp_sl.append(None)
+            self.line_kp_sl[i], = self.canvas.axes.plot(
+                kp, z,
+                linestyle="solid" if i == len(final_z_sl) - 1 else "--",
+                lw=1.,
+                color=(
+                    self.color_plot_river_bottom if i == len(final_z_sl) - 1
+                    else self.colors[i]
+                )
+            )
+
+        self._initial_sl = initial_sl
+        self._river_bottom = final_z_sl[-1]
+
+    def sl_compute_initial(self, reach):
+        """
+        Get SL list for profile p at initial time (initial data)
+        """
+        return list(
+            map(
+                lambda p: p.get_ts_key(min(self._timestamps), "sl")[0],
+                reach.profiles
+            )
+        )
+
+    def sl_compute_current(self, reach):
+        """
+        Get SL list for profile p at current time
+        """
+        return list(
+            map(
+                lambda p: p.get_ts_key(self._current_timestamp, "sl")[0],
+                reach.profiles
+            )
+        )
+
+    def sl_completed_layers(self, initial_sl, current_sl, max_sl_num):
+        sl = []
+        sl_init = []
+        for i in range(max_sl_num):
+            cur = []
+            cur_init = []
+            for current, initial in zip(current_sl, initial_sl):
+                if i < len(initial_sl):
+                    cur.append(current[i][0])
+                    cur_init.append(initial[i][0])
+                else:
+                    cur.append(0)
+                    cur_init.append(0)
+            sl.append(cur)
+            sl_init.append(cur_init)
+
+        return sl_init, sl
+
+    def sl_apply_z_min_to_initial(self, sl_init, z_min):
+        """
+        Compute sediment layer from initial data in function to profile
+        z_min
+        """
+        return reduce(
+            lambda acc, v: acc + [
+                list(
+                    map(
+                        lambda x, y: y - x,
+                        v, acc[-1]
+                    )
+                )
+            ],
+            sl_init, [z_min]
+        )
+
+    def sl_compute_diff(self, sl_init, sl):
+        return list(
+            map(
+                lambda ln0, lni: list(
+                    map(
+                        lambda z0, zi: z0 - zi,
+                        ln0, lni
+                    )
+                ),
+                sl_init, sl
+            )
+        )
+
+    def sl_apply_diff(self, reach, z_sl, d_sl):
+        f = list(map(lambda p: 0, reach.profiles))
+
+        return list(
+            map(
+                lambda z, d: list(
+                    map(
+                        lambda zn, dn: zn - dn,
+                        z, d
+                    )
+                ),
+                z_sl, d_sl + [f] # HACK: Add dummy data for last layer
+            )
+        )
+
+    def draw_bottom_geometry(self, reach):
         kp = reach.geometry.get_kp()
         z_min = reach.geometry.get_z_min()
         z_max = reach.geometry.get_z_max()
@@ -96,6 +228,8 @@ class PlotKPC(PamhyrPlot):
             lw=1.
         )
 
+        self._river_bottom = z_min
+
     def draw_water_elevation(self, reach):
         if len(reach.geometry.profiles) != 0:
             kp = reach.geometry.get_kp()
@@ -116,7 +250,7 @@ class PlotKPC(PamhyrPlot):
             )
 
             self.water_fill = self.canvas.axes.fill_between(
-                kp, z_min, water_z,
+                kp, self._river_bottom, water_z,
                 color=self.color_plot_river_water_zone,
                 alpha=0.7,
                 interpolate=True
@@ -214,6 +348,7 @@ class PlotKPC(PamhyrPlot):
         if not self._init:
             self.draw()
 
+        self.update_bottom_with_bedload()
         self.update_water_elevation()
 
         self.update_idle()
@@ -238,7 +373,7 @@ class PlotKPC(PamhyrPlot):
 
         self.water_fill.remove()
         self.water_fill = self.canvas.axes.fill_between(
-            kp, z_min, water_z,
+            kp, self._river_bottom, water_z,
             color=self.color_plot_river_water_zone,
             alpha=0.7, interpolate=True
         )
@@ -248,14 +383,38 @@ class PlotKPC(PamhyrPlot):
         kp = reach.geometry.get_kp()
         z_min = reach.geometry.get_z_min()
         z_max = reach.geometry.get_z_max()
+        cid = self._current_profile_id
+
         self.profile.set_data(
-            [
-                kp[self._current_profile_id],
-                kp[self._current_profile_id]
-            ],
-            [
-                z_max[self._current_profile_id],
-                z_min[self._current_profile_id]
-            ]
+            [kp[cid], kp[cid]],
+            [z_max[cid], z_min[cid]]
         )
         self.canvas.figure.canvas.draw_idle()
+
+    def update_bottom_with_bedload(self):
+        reach = self.results.river.reach(self._current_reach_id)
+        kp = reach.geometry.get_kp()
+        z_min = reach.geometry.get_z_min()
+
+        initial_sl = self._initial_sl
+        current_sl = self.sl_compute_current(reach)
+
+        max_sl_num = reduce(
+            lambda acc, sl: max(acc, len(sl)),
+            current_sl, 0
+        )
+
+        sl_init, sl = self.sl_completed_layers(
+            initial_sl, current_sl, max_sl_num
+        )
+
+        z_sl = self.sl_apply_z_min_to_initial(sl_init, z_min)
+        d_sl = self.sl_compute_diff(sl_init, sl)
+
+        final_z_sl = self.sl_apply_diff(reach, z_sl, d_sl)
+        final_z_sl = list(reversed(final_z_sl))
+
+        for i, z in enumerate(final_z_sl):
+            self.line_kp_sl[i].set_data(kp, z)
+
+        self._river_bottom = final_z_sl[-1]
diff --git a/src/View/Results/Window.py b/src/View/Results/Window.py
index 7a7c49dd5f1a36a3232ab5a8fd3f4544d53a2268..a1b4d086e56674243102bbb1fe0d4c47da8488fd 100644
--- a/src/View/Results/Window.py
+++ b/src/View/Results/Window.py
@@ -230,45 +230,6 @@ class ResultsWindow(PamhyrWindow):
         )
         self.plot_h.draw()
 
-        # self.canvas_5 = MplCanvas(width=5, height=4, dpi=100)
-        # self.canvas_5.setObjectName("canvas_5")
-        # self.toolbar_5 = PamhyrPlotToolbar(
-        #     self.canvas_5, self
-        # )
-        # self.plot_layout_5 = self.find(QVBoxLayout, "verticalLayout_sed_reach")
-        # self.plot_layout_5.addWidget(self.toolbar_5)
-        # self.plot_layout_5.addWidget(self.canvas_5)
-
-        # if self._study.river.has_sediment():
-        #     self.plot_sed_reach = PlotSedReach(
-        #         canvas=self.canvas_5,
-        #         results=self._results,
-        #         reach_id=0,
-        #         profile_id=0,
-        #         toolbar=self.toolbar_5
-        #     )
-        #     self.plot_sed_reach.draw()
-
-        # self.canvas_6 = MplCanvas(width=5, height=4, dpi=100)
-        # self.canvas_6.setObjectName("canvas_6")
-        # self.toolbar_6 = PamhyrPlotToolbar(
-        #     self.canvas_6, self
-        # )
-        # self.plot_layout_6 = self.find(
-        #     QVBoxLayout, "verticalLayout_sed_profile")
-        # self.plot_layout_6.addWidget(self.toolbar_6)
-        # self.plot_layout_6.addWidget(self.canvas_6)
-
-        # if self._study.river.has_sediment():
-        #     self.plot_sed_profile = PlotSedProfile(
-        #         canvas=self.canvas_6,
-        #         results=self._results,
-        #         reach_id=0,
-        #         profile_id=0,
-        #         toolbar=self.toolbar_6
-        #     )
-        #     self.plot_sed_profile.draw()
-
     def closeEvent(self, event):
         try:
             self._timer.stop()
@@ -395,10 +356,6 @@ class ResultsWindow(PamhyrWindow):
             self.plot_kpc.set_reach(reach_id)
             self.plot_h.set_reach(reach_id)
 
-            if self._study.river.has_sediment():
-                self.plot_sed_reach.set_reach(reach_id)
-                self.plot_sed_profile.set_reach(reach_id)
-
             for plot in self._additional_plot:
                 self._additional_plot[plot].set_reach(reach_id)
 
@@ -411,10 +368,6 @@ class ResultsWindow(PamhyrWindow):
             self.plot_kpc.set_profile(profile_id)
             self.plot_h.set_profile(profile_id)
 
-            if self._study.river.has_sediment():
-                self.plot_sed_reach.set_profile(profile_id)
-                self.plot_sed_profile.set_profile(profile_id)
-
             for plot in self._additional_plot:
                 self._additional_plot[plot].set_profile(profile_id)
 
@@ -426,10 +379,6 @@ class ResultsWindow(PamhyrWindow):
             self.plot_kpc.set_timestamp(timestamp)
             self.plot_h.set_timestamp(timestamp)
 
-            if self._study.river.has_sediment():
-                self.plot_sed_reach.set_timestamp(timestamp)
-                self.plot_sed_profile.set_timestamp(timestamp)
-
             for plot in self._additional_plot:
                 self._additional_plot[plot].set_timestamp(timestamp)
 
@@ -500,19 +449,11 @@ class ResultsWindow(PamhyrWindow):
         self.plot_kpc.results = self._results
         self.plot_h.results = self._results
 
-        if self._study.river.has_sediment():
-            self.plot_sed_reach.results = self._results
-            self.plot_sed_profile.results = self._results
-
         self.plot_xy.draw()
         self.plot_ac.draw()
         self.plot_kpc.draw()
         self.plot_h.draw()
 
-        if self._study.river.has_sediment():
-            self.plot_sed_reach.draw()
-            self.plot_sed_profile.draw()
-
     def _reload_slider(self):
         self._slider_time = self.find(QSlider, f"horizontalSlider_time")
         self._slider_time.setMaximum(len(self._timestamps) - 1)