diff --git a/Makefile b/Makefile
index 3bbede6..4fea72d 100644
--- a/Makefile
+++ b/Makefile
@@ -121,6 +121,17 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
 	--profile-type "prestorm" \
 	--output-file "./data/interim/twl_mean_slope_hol86.csv"
 
+./data/interim/twl_mean_slope_nie91.csv: ./data/interim/waves.csv ./data/interim/tides.csv ./data/interim/profiles.csv ./data/interim/sites.csv ./data/interim/profile_features_crest_toes.csv
+	$(PYTHON_CLI) create-twl-forecast \
+	--waves-csv "./data/interim/waves.csv" \
+	--tides-csv "./data/interim/tides.csv" \
+	--profiles-csv "./data/interim/profiles.csv" \
+	--profile-features-csv "./data/interim/profile_features_crest_toes.csv" \
+	--runup-function "nie91" \
+	--slope "mean" \
+	--profile-type "prestorm" \
+	--output-file "./data/interim/twl_mean_slope_nie91.csv"
+
 ### IMPACTS
 
 ./data/interim/impacts_observed.csv: ./data/interim/profiles.csv ./data/interim/profile_features_crest_toes.csv ./data/raw/profile_features_chris_leaman/profile_features_chris_leaman.csv
@@ -148,9 +159,15 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
 	--forecasted-twl-csv "./data/interim/twl_mean_slope_hol86.csv" \
 	--output-file "./data/interim/impacts_forecasted_mean_slope_hol86.csv"
 
+./data/interim/impacts_forecasted_mean_slope_nie91.csv: ./data/interim/profile_features_crest_toes.csv ./data/interim/twl_mean_slope_nie91.csv
+	$(PYTHON_CLI) create-forecasted-impacts \
+	--profile-features-csv "./data/interim/profile_features_crest_toes.csv" \
+	--forecasted-twl-csv "./data/interim/twl_mean_slope_nie91.csv" \
+	--output-file "./data/interim/impacts_forecasted_mean_slope_nie91.csv"
+
 ### GEOJSONs
 
-geojsons: ./data/interim/impacts_forecasted_mean_slope_hol86.geojson ./data/interim/impacts_forecasted_mean_slope_hol86_R_high.geojson ./data/interim/impacts_forecasted_mean_slope_sto06.geojson ./data/interim/impacts_forecasted_mean_slope_sto06_R_high.geojson ./data/interim/profile_features_crest_toes.geojson ./data/interim/sites.geojson
+geojsons: ./data/interim/impacts_forecasted_mean_slope_nie91.geojson ./data/interim/impacts_forecasted_mean_slope_nie91_R_high.geojson ./data/interim/impacts_forecasted_mean_slope_hol86.geojson ./data/interim/impacts_forecasted_mean_slope_hol86_R_high.geojson ./data/interim/impacts_forecasted_mean_slope_sto06.geojson ./data/interim/impacts_forecasted_mean_slope_sto06_R_high.geojson ./data/interim/profile_features_crest_toes.geojson ./data/interim/sites.geojson
 
 ./data/interim/impacts_forecasted_mean_slope_sto06.geojson: ./data/interim/impacts_forecasted_mean_slope_sto06.csv ./data/interim/impacts_observed.csv
 	$(PYTHON_CLI) impacts-to-geojson \
@@ -182,6 +199,21 @@ geojsons: ./data/interim/impacts_forecasted_mean_slope_hol86.geojson ./data/inte
 	--impacts-csv "./data/interim/impacts_forecasted_mean_slope_hol86.csv" \
 	--output-geojson "./data/interim/impacts_forecasted_mean_slope_hol86_R_high.geojson"
 
+./data/interim/impacts_forecasted_mean_slope_nie91.geojson: ./data/interim/impacts_forecasted_mean_slope_nie91.csv ./data/interim/impacts_observed.csv
+	$(PYTHON_CLI) impacts-to-geojson \
+	--sites-csv "./data/interim/sites.csv" \
+	--observed-impacts-csv "./data/interim/impacts_observed.csv" \
+	--forecast-impacts-csv "./data/interim/impacts_forecasted_mean_slope_nie91.csv" \
+	--output-geojson "./data/interim/impacts_forecasted_mean_slope_nie91.geojson"
+
+./data/interim/impacts_forecasted_mean_slope_nie91_R_high.geojson: ./data/interim/impacts_forecasted_mean_slope_nie91.csv
+	$(PYTHON_CLI) r-high-to-geojson \
+	--sites-csv "./data/interim/sites.csv" \
+	--profiles-csv "./data/interim/profiles.csv" \
+	--crest-toes-csv "./data/interim/profile_features_crest_toes.csv" \
+	--impacts-csv "./data/interim/impacts_forecasted_mean_slope_nie91.csv" \
+	--output-geojson "./data/interim/impacts_forecasted_mean_slope_nie91_R_high.geojson"
+
 ./data/interim/profile_features_crest_toes.geojson: ./data/interim/profile_features_crest_toes.csv
 	$(PYTHON_CLI) profile-features-crest-toes-to-geojson \
 	--sites-csv "./data/interim/sites.csv" \
diff --git a/src/analysis/forecast_twl.py b/src/analysis/forecast_twl.py
index 9badbe2..8815d4c 100644
--- a/src/analysis/forecast_twl.py
+++ b/src/analysis/forecast_twl.py
@@ -325,7 +325,10 @@ def slope_from_profile(
 @click.option("--profiles-csv", required=True, help="")
 @click.option("--profile-features-csv", required=True, help="")
 @click.option(
-    "--runup-function", required=True, help="", type=click.Choice(["sto06", "hol86"])
+    "--runup-function",
+    required=True,
+    help="",
+    type=click.Choice(["sto06", "hol86", "nie91"]),
 )
 @click.option(
     "--slope", required=True, help="", type=click.Choice(["foreshore", "mean"])
diff --git a/src/analysis/runup_models.py b/src/analysis/runup_models.py
index 374858f..dc0a06d 100644
--- a/src/analysis/runup_models.py
+++ b/src/analysis/runup_models.py
@@ -67,6 +67,52 @@ def hol86(Hs0, Tp, beta):
     )
 
 
+def nie91(Hs0, Tp, beta):
+
+    df = pd.DataFrame(
+        {"Hs0": Hs0, "Tp": Tp, "beta": beta}, index=[x for x in range(0, np.size(Hs0))]
+    )
+
+    df["Lp"] = 9.8 * df["Tp"] ** 2 / 2 / np.pi
+    df["Ls"] = df["Lp"]  # Need to make this approximation, refer to Atkinson 2017
+    df["Hrms"] = df["Hs0"] / np.sqrt(
+        2
+    )  # Acceptable approximation, refer to Atkinson 2017
+
+    df.loc[df.beta >= 0.1, "LR"] = 0.6 * df["beta"] * np.sqrt(df["Hrms"] * df["Ls"])
+    df.loc[df.beta < 0.1, "LR"] = 0.06 * np.sqrt(df["Hrms"] * df["Ls"])
+
+    df["R2"] = 1.98 * df["LR"]
+
+    # Note that this should be the level above Z100%, which in this case is taken as the time varying tide level,
+    # even though minimum run-down can still occur below tide SWL.
+
+    df["setup"] = np.nan
+    df["S_ig"] = np.nan
+    df["S_inc"] = np.nan
+    df["S_total"] = np.nan
+
+    return (
+        float_or_list(df["R2"].tolist()),
+        float_or_list(df["setup"].tolist()),
+        float_or_list(df["S_total"].tolist()),
+        float_or_list(df["S_inc"].tolist()),
+        float_or_list(df["S_ig"].tolist()),
+    )
+
+
+def atk18(Hs0, Tp, beta):
+    pass
+
+
+def pow18(Hs0, Tp, beta):
+    pass
+
+
+def beu(Hs0, Tp, beta):
+    pass
+
+
 def float_or_list(a):
     """
     If only one value in the array, return the float, else return a list