diff --git a/src/analysis/forecast_twl.py b/src/analysis/forecast_twl.py
index cd02071..70edcab 100644
--- a/src/analysis/forecast_twl.py
+++ b/src/analysis/forecast_twl.py
@@ -49,42 +49,7 @@ def forecast_twl(
         df_twl["beta"] = pd.concat(results)
 
     elif slope == "mean":
-        logger.info("Calculating mean (dune toe to MHW) slopes")
-        btm_z = 0.5  # m AHD
-
-        # When calculating mean slope, we go from the dune toe to mhw. However, in some profiles, the dune toe is not
-        # defined. In these cases, we should go to the dune crest. Let's make a temporary dataframe which has this
-        # already calculated.
-        df_top_ele = df_profile_features.xs(profile_type, level="profile_type").copy()
-        df_top_ele.loc[:, "top_ele"] = df_top_ele.dune_toe_z
-        df_top_ele.loc[
-            df_top_ele.top_ele.isnull().values, "top_ele"
-        ] = df_top_ele.dune_crest_z
-
-        n_no_top_ele = len(df_top_ele[df_top_ele.top_ele.isnull()].index)
-        if n_no_top_ele != 0:
-            logger.warning(
-                "{} sites do not have dune toes/crests to calculate mean slope".format(
-                    n_no_top_ele
-                )
-            )
-
-        df_slopes = (
-            df_profiles.xs(profile_type, level="profile_type")
-            .dropna(subset=["z"])
-            .groupby("site_id")
-            .apply(
-                lambda x: slope_from_profile(
-                    profile_x=x.index.get_level_values("x").tolist(),
-                    profile_z=x.z.tolist(),
-                    top_elevation=df_top_ele.loc[x.index[0][0], :].top_ele,
-                    btm_elevation=btm_z,
-                    method="least_squares",
-                )
-            )
-            .rename("beta")
-            .to_frame()
-        )
+        df_slopes = get_mean_slope(df_profile_features, df_profiles, profile_type)
 
         # Merge calculated slopes onto each twl timestep
         df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True)
@@ -92,26 +57,7 @@ def forecast_twl(
     elif slope == "intertidal":
 
         logger.info("Calculating intertidal slopes")
-        top_z = 1.15  # m AHD = HAT from MHL annual ocean tides summary report
-        btm_z = -0.9  # m AHD = HAT from MHL annual ocean tides summary report
-
-        # Calculate slopes for each profile
-        df_slopes = (
-            df_profiles.xs(profile_type, level="profile_type")
-            .dropna(subset=["z"])
-            .groupby("site_id")
-            .apply(
-                lambda x: slope_from_profile(
-                    profile_x=x.index.get_level_values("x").tolist(),
-                    profile_z=x.z.tolist(),
-                    top_elevation=top_z,
-                    btm_elevation=max(min(x.z), btm_z),
-                    method="least_squares",
-                )
-            )
-            .rename("beta")
-            .to_frame()
-        )
+        df_slopes = get_intertidal_slope(df_profiles, profile_type)
 
         # Merge calculated slopes onto each twl timestep
         df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True)
@@ -134,12 +80,83 @@ def forecast_twl(
         df_twl["tide"] + 1.1 * df_twl["setup"] - 1.1 / 2 * df_twl["S_total"]
     )
 
-    # Drop unneeded columns
-    # df_twl.drop(columns=["E", "Exs", "P", "Pxs", "dir"], inplace=True, errors="ignore")
-
     return df_twl
 
 
+def get_intertidal_slope(df_profiles, profile_type, top_z=1.15, btm_z=-0.9):
+    """
+    Gets intertidal slopes
+    :param df_profiles:
+    :param profile_type:
+    :return:
+    """
+
+    # Calculate slopes for each profile
+    df_slopes = (
+        df_profiles.xs(profile_type, level="profile_type")
+        .dropna(subset=["z"])
+        .groupby("site_id")
+        .apply(
+            lambda x: slope_from_profile(
+                profile_x=x.index.get_level_values("x").tolist(),
+                profile_z=x.z.tolist(),
+                top_elevation=top_z,
+                btm_elevation=max(min(x.z), btm_z),
+                method="least_squares",
+            )
+        )
+        .rename("beta")
+        .to_frame()
+    )
+    return df_slopes
+
+
+def get_mean_slope(df_profile_features, df_profiles, profile_type, btm_z=0.5):
+    """
+    Calculates the mean slopes for all profiles
+    :param df_profile_features:
+    :param df_profiles:
+    :param profile_type:
+    :param btm_z: Typically mean high water
+    :return:
+    """
+
+    logger.info("Calculating mean (dune toe to MHW) slopes")
+
+    # When calculating mean slope, we go from the dune toe to mhw. However, in some profiles, the dune toe is not
+    # defined. In these cases, we should go to the dune crest. Let's make a temporary dataframe which has this
+    # already calculated.
+    df_top_ele = df_profile_features.xs(profile_type, level="profile_type").copy()
+    df_top_ele.loc[:, "top_ele"] = df_top_ele.dune_toe_z
+    df_top_ele.loc[
+        df_top_ele.top_ele.isnull().values, "top_ele"
+    ] = df_top_ele.dune_crest_z
+    n_no_top_ele = len(df_top_ele[df_top_ele.top_ele.isnull()].index)
+    if n_no_top_ele != 0:
+        logger.warning(
+            "{} sites do not have dune toes/crests to calculate mean slope".format(
+                n_no_top_ele
+            )
+        )
+    df_slopes = (
+        df_profiles.xs(profile_type, level="profile_type")
+        .dropna(subset=["z"])
+        .groupby("site_id")
+        .apply(
+            lambda x: slope_from_profile(
+                profile_x=x.index.get_level_values("x").tolist(),
+                profile_z=x.z.tolist(),
+                top_elevation=df_top_ele.loc[x.index[0][0], :].top_ele,
+                btm_elevation=btm_z,
+                method="least_squares",
+            )
+        )
+        .rename("beta")
+        .to_frame()
+    )
+    return df_slopes
+
+
 def mean_slope_for_site_id(
     site_id,
     df_twl,