diff --git a/src/Model/Results/River/River.py b/src/Model/Results/River/River.py
index aaa29a75..6e32b7a3 100644
--- a/src/Model/Results/River/River.py
+++ b/src/Model/Results/River/River.py
@@ -104,7 +104,7 @@ class Profile(SQLSubModel):
         return any(map(lambda ts: "sl" in self._data[ts], self._data))
 
     def has_bedload(self):
-        return any(map(lambda ts: "zfd" in self._data[ts], self._data))
+        return any(map(lambda ts: "sl" in self._data[ts], self._data))
 
     @classmethod
     def _db_create(cls, execute, ext=""):
diff --git a/src/Solver/Mage.py b/src/Solver/Mage.py
index ddb3fcdf..6ecedbe6 100644
--- a/src/Solver/Mage.py
+++ b/src/Solver/Mage.py
@@ -1393,12 +1393,14 @@ class Mage8(Mage):
 
             logger.info(f"compute river bed elevation...")
 
+            zfd_lst = []
+
             for r in reachs:
                 z_min = r.geometry.get_z_min()
-                sls = list(map(
+                sls = map(
                     lambda p: p.get_ts_key(ts_list[0], "sl")[0],
                     r.profiles
-                    ))
+                )
                 z_br = list(map(
                     lambda z, sl: reduce(
                         lambda z, h: z - h[0],
@@ -1407,11 +1409,11 @@ class Mage8(Mage):
                     z_min,          # Original geometry
                     sls             # Original sediment layers
                 ))
-                for t in ts_list:
-                    sls = list(map(
+                for its, t in enumerate(ts_list):
+                    sls = map(
                         lambda p: p.get_ts_key(t, "sl")[0],
                         r.profiles
-                    ))
+                    )
                     zfd = list(map(
                         lambda z, sl: reduce(
                             lambda z, h: z + h[0],
@@ -1420,8 +1422,13 @@ class Mage8(Mage):
                         z_br,       # bedrock
                         sls         # Original sediment layers
                     ))
-                    for i, p in enumerate(r.profiles):
-                        r.set(i, t, "zfd", zfd[i])
+
+                    zfd_lst.append(np.array(zfd))
+
+            table = results.get("table")
+            table["zfd"] = results.new_table_data(
+                "zfd", np.array(zfd_lst)
+            )
 
             results.set("sediment_timestamps", ts)
             logger.info(f"read_gra: ... end with {len(ts)} timestamp read")
diff --git a/src/View/Results/PlotRKC.py b/src/View/Results/PlotRKC.py
index 96c12adb..dc658c65 100644
--- a/src/View/Results/PlotRKC.py
+++ b/src/View/Results/PlotRKC.py
@@ -103,12 +103,15 @@ class PlotRKC(PamhyrPlot):
             self.draw_bottom_geometry(reach)
 
     def draw_bottom_with_bedload(self, reach):
+        results = self.results[self._current_res_id]
         rk = reach.geometry.get_rk()
+
+        table = results.get("table")["zfd"]
+        ts = results.get_timestamp_id(self._current_timestamp)
+
         zfd = list(
             map(
-                lambda p: p.get_ts_key(
-                    self._current_timestamp, "zfd"
-                ),
+                lambda p: table[ts, p.global_index],
                 reach.profiles
             )
         )
@@ -396,14 +399,13 @@ class PlotRKC(PamhyrPlot):
         rk = reach.geometry.get_rk()
         # z = self.sl_compute_current_z(reach)
 
-        zfd = list(
-            map(
-                lambda p: p.get_ts_key(
-                    self._current_timestamp, "zfd"
-                ),
-                reach.profiles
-            )
-        )
+        table = results.get("table")["zfd"]
+        ts = results.get_timestamp_id(self._current_timestamp)
+
+        zfd = list(map(
+            lambda p: table[ts, p.global_index],
+            reach.profiles
+        ))
 
         self.line_bottom.remove()
 
diff --git a/src/View/Results/Window.py b/src/View/Results/Window.py
index c9bf9667..9dcc0bec 100644
--- a/src/View/Results/Window.py
+++ b/src/View/Results/Window.py
@@ -891,17 +891,14 @@ class ResultsWindow(PamhyrWindow):
         z = table["Z"]
         q = table["Q"]
         v = table["V"]
+        zfd = table["zfd"]
 
         if "bed_elevation" in y:
             if reach.has_bedload():
-                zmin = list(
-                        map(
-                            lambda p: p.get_ts_key(
-                                timestamp, "zfd"
-                            ),
-                            reach.profiles
-                        )
-                    )
+                zmin = list(map(
+                    lambda p: zfd[id_ts, p.global_index],
+                    reach.profiles
+                ))
             else:
                 zmin = reach.geometry.get_z_min()
             my_dict[dict_y["bed_elevation"]] = zmin
@@ -1153,6 +1150,7 @@ class ResultsWindow(PamhyrWindow):
         z = table["Z"][:, profile.global_index]
         q = table["Q"][:, profile.global_index]
         v = table["V"][:, profile.global_index]
+        zfd = table["zfd"][:, profile.global_index]
 
         if self._current_results == 2:
             reach1 = self._results[0].river.reach(self._reach)
@@ -1172,7 +1170,7 @@ class ResultsWindow(PamhyrWindow):
 
         if "bed_elevation" in y:
             if reach.has_bedload():
-                z_min = profile.get_key("zfd")
+                z_min = zfd
             else:
                 z_min = [profile.geometry.z_min()] * len(self._timestamps)