diff --git a/Makefile b/Makefile
index 8205197..3dc2c8c 100644
--- a/Makefile
+++ b/Makefile
@@ -100,6 +100,7 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
 	--tides-csv "./data/interim/tides.csv" \
 	--profiles-csv "./data/interim/profiles.csv" \
 	--profile-features-csv "./data/interim/profile_features_crest_toes.csv" \
+	--grain-size-csv "./data/interim/sites_grain_size.csv" \
 	--runup-function "sto06" \
 	--slope "foreshore" \
 	--profile-type "prestorm" \
@@ -111,6 +112,7 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
 	--tides-csv "./data/interim/tides.csv" \
 	--profiles-csv "./data/interim/profiles.csv" \
 	--profile-features-csv "./data/interim/profile_features_crest_toes.csv" \
+	--grain-size-csv "./data/interim/sites_grain_size.csv" \
 	--runup-function "sto06" \
 	--slope "mean" \
 	--profile-type "prestorm" \
@@ -122,6 +124,7 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
 	--tides-csv "./data/interim/tides.csv" \
 	--profiles-csv "./data/interim/profiles.csv" \
 	--profile-features-csv "./data/interim/profile_features_crest_toes.csv" \
+	--grain-size-csv "./data/interim/sites_grain_size.csv" \
 	--runup-function "hol86" \
 	--slope "mean" \
 	--profile-type "prestorm" \
@@ -133,6 +136,7 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
 	--tides-csv "./data/interim/tides.csv" \
 	--profiles-csv "./data/interim/profiles.csv" \
 	--profile-features-csv "./data/interim/profile_features_crest_toes.csv" \
+	--grain-size-csv "./data/interim/sites_grain_size.csv" \
 	--runup-function "nie91" \
 	--slope "mean" \
 	--profile-type "prestorm" \
diff --git a/src/analysis/forecast_twl.py b/src/analysis/forecast_twl.py
index 8815d4c..60e58e1 100644
--- a/src/analysis/forecast_twl.py
+++ b/src/analysis/forecast_twl.py
@@ -20,6 +20,7 @@ def forecast_twl(
     df_profiles,
     df_waves,
     df_profile_features,
+    df_grain_size,
     runup_function,
     n_processes=MULTIPROCESS_THREADS,
     slope="foreshore",
@@ -82,6 +83,7 @@ def forecast_twl(
         Hs0=df_twl["Hs0"].tolist(),
         Tp=df_twl["Tp"].tolist(),
         beta=df_twl["beta"].tolist(),
+        r=df_twl.merge(df_grain_size, on="site_id").r.tolist(),
     )
 
     df_twl["R2"] = R2
@@ -340,6 +342,7 @@ def slope_from_profile(
     type=click.Choice(["prestorm", "poststorm"]),
 )
 @click.option("--output-file", required=True, help="")
+@click.option("--grain-size-csv", required=False, help="")
 def create_twl_forecast(
     waves_csv,
     tides_csv,
@@ -349,6 +352,7 @@ def create_twl_forecast(
     slope,
     profile_type,
     output_file,
+    grain_size_csv,
 ):
     logger.info("Creating forecast of total water levels")
     logger.info("Importing data")
@@ -356,6 +360,7 @@ def create_twl_forecast(
     df_tides = pd.read_csv(tides_csv, index_col=[0, 1])
     df_profiles = pd.read_csv(profiles_csv, index_col=[0, 1, 2])
     df_profile_features = pd.read_csv(profile_features_csv, index_col=[0, 1])
+    df_grain_size = pd.read_csv(grain_size_csv, index_col=[0])
 
     logger.info("Forecasting TWL")
     df_twl = forecast_twl(
@@ -363,6 +368,7 @@ def create_twl_forecast(
         df_profiles,
         df_waves,
         df_profile_features,
+        df_grain_size,
         runup_function=getattr(runup_models, runup_function),
         slope=slope,
         profile_type=profile_type,
diff --git a/src/analysis/runup_models.py b/src/analysis/runup_models.py
index dc0a06d..b7ae476 100644
--- a/src/analysis/runup_models.py
+++ b/src/analysis/runup_models.py
@@ -2,7 +2,8 @@ import numpy as np
 import pandas as pd
 
 
-def sto06(Hs0, Tp, beta):
+
+def sto06(Hs0, Tp, beta, **kwargs):
     """
     :param Hs0: List or float of offshore significant wave height values
     :param Tp: List or float of peak wave period
@@ -43,8 +44,7 @@ def sto06(Hs0, Tp, beta):
     )
 
 
-def hol86(Hs0, Tp, beta):
-
+def hol86(Hs0, Tp, beta, **kwargs):
     df = pd.DataFrame(
         {"Hs0": Hs0, "Tp": Tp, "beta": beta}, index=[x for x in range(0, np.size(Hs0))]
     )
@@ -67,8 +67,7 @@ def hol86(Hs0, Tp, beta):
     )
 
 
-def nie91(Hs0, Tp, beta):
-
+def nie91(Hs0, Tp, beta, **kwargs):
     df = pd.DataFrame(
         {"Hs0": Hs0, "Tp": Tp, "beta": beta}, index=[x for x in range(0, np.size(Hs0))]
     )