diff --git a/R/plot_biomass.R b/R/plot_biomass.R
index 45dbf3f4..725d3a36 100644
--- a/R/plot_biomass.R
+++ b/R/plot_biomass.R
@@ -141,37 +141,18 @@ plot_biomass <- function(
     ...
   )
   # Add reference line
-  # getting data set - an ifelse statement in the fxn wasn't working
   if (relative) {
     # don't add any reference line here and just add theme for final plot
     final <- plt + theme_noaa()
   } else {
-    if ("unfished" %in% c(names(ref_line), ref_line)) {
-      # find the minimum x axis value from the plot
-      min_year <- min <- ggplot2::ggplot_build(plt)[["data"]][[2]] |>
-        as.data.frame() |>
-        dplyr::pull(x) |>
-        min() |>
-        round(digits = 2)
-      # find the reference point value for unfished
-      ref_point <- calculate_reference_point(
-        dat = stockplotr::example_data,
-        reference_name = "biomass_unfished"
-      ) / scale_amount
-      # add point to plot and add theme
-      final <- plt +
-        ggplot2::geom_point(ggplot2::aes(x = min_year - 1, y = ref_point)) + # should I keep -1 or set as first year?
-        theme_noaa()
-    } else {
-      final <- reference_line(
-        plot = plt,
-        dat = rp_dat,
-        label_name = "biomass",
-        reference = ref_line,
+    final <- plt +
+      add_reference_line(
+        dat = dat, 
+        ref_line = ref_line, 
+        label = "biomass",
         scale_amount = scale_amount
       ) +
-        theme_noaa()
-    }
+      theme_noaa()
   }
 
   ### Make RDA ----
diff --git a/R/plot_fishing_mortality.R b/R/plot_fishing_mortality.R
index 56b01401..93b73e1f 100644
--- a/R/plot_fishing_mortality.R
+++ b/R/plot_fishing_mortality.R
@@ -76,14 +76,19 @@ plot_fishing_mortality <- function(
     facet = facet,
     ...
   )
-  # Add reference line and theme
-  final <- reference_line(
-    plot = plt,
-    dat = dat,
-    label_name = "fishing_mortality",
-    reference = ref_line,
-    scale_amount = 1
-  ) + theme_noaa()
+  # Add reference line
+  # Conditions for ref line
+  # 1. all are comparing same value = msy, target, unfished...
+  # 2. custom input vector = c("sable23_msy"=30, "sable25_msy"=40) -- user must indicate which model in label otherwise it will assign in order of dat
+  # 3. input vector of labels = c("msy", "target")
+  # getting data set - an ifelse statement in the fxn wasn't working
+  final <- plt +
+      add_reference_line(
+        dat = dat, 
+        ref_line = ref_line, 
+        label = "fishing_mortality"
+      ) +
+      theme_noaa()
 
   ### Make RDA ----
   if (make_rda) {
diff --git a/R/plot_spawning_biomass.R b/R/plot_spawning_biomass.R
index 33a67fca..cc12c70c 100644
--- a/R/plot_spawning_biomass.R
+++ b/R/plot_spawning_biomass.R
@@ -145,12 +145,6 @@ plot_spawning_biomass <- function(
       )
     }
   )
-  # Pull first df if in a list to find reference point
-  if (!is.data.frame(dat)) {
-    rp_dat <- dat[[1]]
-  } else {
-    rp_dat <- dat
-  }
 
   if (relative & scale_amount > 1) {
     cli::cli_alert_warning("Scale amount is not applicable when relative = TRUE. Resetting scale_amount to 1.")
@@ -202,41 +196,23 @@ plot_spawning_biomass <- function(
     ylab = spawning_biomass_label,
     group = group,
     # add check in case facet is returned as character(0)
-    facet = if (length(facet) > 0) facet else NULL,
-    ...
+    facet = if (length(facet) > 0) facet else NULL#,
+    # ...
   )
   # Add reference line
-  # getting data set - an ifelse statement in the fxn wasn't working
   if (relative) {
     # don't add any reference line here and just add theme for final plot
     final <- plt + theme_noaa()
   } else {
-    if ("unfished" %in% c(names(ref_line), ref_line)) {
-      # find the minimum x axis value from the plot
-      min_year <- min <- ggplot2::ggplot_build(plt)[["data"]][[2]] |>
-        as.data.frame() |>
-        dplyr::pull(x) |>
-        min() |>
-        round(digits = 2)
-      # find the reference point value for unfished
-      ref_point <- calculate_reference_point(
-        dat = rp_dat,
-        reference_name = "spawing_biomass_unfished"
-      ) / scale_amount
-      # add point to plot and add theme
-      final <- plt +
-        ggplot2::geom_point(ggplot2::aes(x = min_year - 1, y = ref_point)) + # should I keep -1 or set as first year?
-        theme_noaa()
-    } else {
-      final <- reference_line(
-        plot = plt,
-        dat = rp_dat,
+    final <- plt + 
+      add_reference_line(
+        dat = dat, 
+        ref_line = ref_line, 
+        label = "spawning_biomass", 
         lbs = lbs,
-        label_name = "spawning_biomass",
-        reference = ref_line,
         scale_amount = scale_amount
-      ) + theme_noaa()
-    }
+      ) + 
+      theme_noaa()
   }
 
   ### Make RDA ----
diff --git a/R/utils_plot.R b/R/utils_plot.R
index 44817291..eeaffbd1 100644
--- a/R/utils_plot.R
+++ b/R/utils_plot.R
@@ -71,30 +71,19 @@ plot_timeseries <- function(
 ) {
   # Start plot
   plot <- ggplot2::ggplot()
-  # make into new geom?
-  # more defaults and fxnality for ggplot
 
   # Add geom
   plot <- switch(geom,
     "point" = {
-      # point_size = ifelse(
-      #   is.null(list(...)$size),
-      #   2.0,
-      #   list(...)$size
-      # )
       plot +
         ggplot2::geom_point(
           data = dat,
           ggplot2::aes(
             .data[[x]],
             .data[[y]],
-            # TODO: add more groupings
-            # shape = ifelse(any(grepl("shape", names(group))), .data[[group[[grep("shape", names(group))]]]], 1),
-            # color = ifelse(any(grepl("color", names(group))), .data[[group[[grep("color", names(group))]]]], "black")
             color = model,
             shape = group_var
           ),
-          # size = point_size,
           ...
         )
     },
@@ -119,7 +108,8 @@ plot_timeseries <- function(
               }
             }
           ),
-          alpha = 0.3
+          alpha = 0.3 #,
+          # show.legend = ifelse(all(is.na(dat$estimate_lower)), FALSE, TRUE)
         ) +
         # }
         ggplot2::geom_line(
@@ -127,8 +117,6 @@ plot_timeseries <- function(
           ggplot2::aes(
             x = .data[[x]],
             y = .data[[y]],
-            # linetype = group_var,
-            # linetype = ifelse(!is.null(group), group_var, "solid"),
             color = {
               if (length(unique(.data[["model"]])) > 1) {
                 interaction(model, group_var)
@@ -137,7 +125,7 @@ plot_timeseries <- function(
               }
             }
           ),
-          # linewidth = 1.0,
+          # show.legend = FALSE, #ifelse(all(is.na(dat$estimate_lower)), TRUE, FALSE),
           ...
         )
     },
@@ -148,7 +136,7 @@ plot_timeseries <- function(
           ggplot2::aes(
             x = .data[[x]],
             y = .data[[y]],
-            fill = model
+            fill = interaction(model, group_var) # model
           ),
           ...
         )
@@ -160,10 +148,6 @@ plot_timeseries <- function(
     labs <- plot + ggplot2::labs(
       x = xlab,
       y = ylab
-      # color = "Model",
-      # linetype = cap_first_letter(group),
-      # fill = cap_first_letter(group),
-      # shape = cap_first_letter(group)
     ) +
       ggplot2::theme(legend.title = ggplot2::element_blank())
   } else {
@@ -195,15 +179,7 @@ plot_timeseries <- function(
       # return plot if option beyond line and point for now
       labs
     )
-  } # else if (length(unique(dat$model)) == 1) {
-  #   labs <- switch(geom,
-  #     "line" = labs + ggplot2::guides(color = "none"),
-  #     "point" = labs + ggplot2::guides(color = "none"),
-  #     "area" = labs + ggplot2::guides(fill = "none"),
-  #     # return plot if option beyond line and point for now
-  #     labs
-  #   )
-  # }
+  }
 
   # Calc axis breaks
   x_n_breaks <- axis_breaks(dat[[x]])
@@ -529,11 +505,9 @@ cohort_line <- function(
 #' Preformatted reference line
 #'
 #' @inheritParams plot_spawning_biomass
-#' @param plot Plot object. A ggplot2 object where the reference line will be added
-#' @param dat Data frame. Standard data frame where reference point should be extracted
 #' @param label_name String. Name of the quantity that users want to
 #' extract the reference point from
-#' @param reference String of the reference point
+#' @param ref_line String of the reference point
 #'
 #' Options: Including, but not limited to: "msy", "unfished", "target"
 #'
@@ -546,49 +520,39 @@ cohort_line <- function(
 #' reference_line(dat, "biomass", "msy")
 #' }
 reference_line <- function(
-  plot,
-  dat,
   label_name,
-  reference,
+  ref_line,
   scale_amount = 1,
-  lbs = FALSE
+  model_name = "1"
 ) {
-  if (!is.null(names(reference))) {
-    ref_line_val <- reference[[1]]
-    reference <- names(reference)
-  } else {
-    # calculate reference point value
-    ref_line_val <- calculate_reference_point(
-      dat = dat,
-      reference_name = glue::glue("{label_name}_{reference}"),
-      lbs = lbs
-    )
-  }
+  # Extract reference line values and labels
+  ref_line_val <- ref_line[[1]]
+  reference <- names(ref_line)
 
-  # Add geom for ref line
-  if (is.null(ref_line_val)) {
-    plot
-  } else {
-    plot +
-      ggplot2::geom_hline(
-        # ggplot2::aes(
+  list(
+    # Add geom for ref line
+    ggplot2::geom_hline(
+      data = data.frame(yintercept = ref_line_val, model = model_name, group_var = "1"),
+      ggplot2::aes(
         yintercept = ref_line_val / scale_amount,
-        # ),
-        color = "black",
-        linetype = "dashed"
-      ) +
-      ggplot2::annotate(
-        geom = "text",
-        # TODO: need to change this for general process
-        x = as.numeric(max(ggplot2::ggplot_build(plot)[["data"]][[2]][["x"]], na.rm = TRUE)), # - as.numeric(max(dat$year[dat$era == "time"], na.rm = TRUE))/200,
-        y = ref_line_val / scale_amount,
-        label = glue::glue("{stringr::str_replace_all(label_name, '_', '~')}[{reference}]"), # list(bquote(label_name[.(reference)])),
-        parse = TRUE,
-        hjust = 1,
-        vjust = 0,
-        size = 5 # this is not foolproof
-      )
-  }
+        color = interaction(model, group_var) # model
+      ),
+      linetype = "dashed",
+      show.legend = FALSE
+    ),
+    # add line annotation
+    ggplot2::annotate(
+      geom = "text",
+      # TODO: need to change this for general process
+      x = -Inf,
+      y = ref_line_val / scale_amount,
+      label = glue::glue("{reference}"), # glue::glue("{stringr::str_replace_all(label_name, '_', '~')}[{reference}]"), # list(bquote(label_name[.(reference)])),
+      parse = TRUE,
+      hjust = -0.05, # slight offset so text doesn't hit border
+      vjust = -0.5 # , #slightly above line
+      # size = 5 # this is not foolproof
+    )
+  )
 }
 
 #------------------------------------------------------------------------------
@@ -1062,3 +1026,77 @@ plot_obsvpred <- function(
   }
   final
 }
+
+#------------------------------------------------------------------------------
+
+add_reference_line <- function(
+    dat,
+    ref_line,
+    label,
+    lbs = FALSE,
+    scale_amount = 1
+) {
+  # Add reference line
+  # Conditions for ref line
+  # 1. all are comparing same value = msy, target, unfished...
+  # 2. custom input vector = c("sable23_msy"=30, "sable25_msy"=40) -- user must indicate which model in label otherwise it will assign in order of dat
+  # 3. input vector of labels = c("msy", "target")
+  # getting data set - an ifelse statement in the fxn wasn't working
+  ref_lines_list <- list()
+  if (!is.null(ref_line)) {
+    # Check if length of ref_line = dat -- replicate value if not
+    if (!is.data.frame(dat) & length(ref_line) != length(dat)) ref_line <- rep(ref_line, length(dat))
+    # Put into for loop and add lines sequentially to plt
+    for (i in 1:length(ref_line)) {
+      # find the reference point value
+      if (is.null(names(ref_line[i]))) {
+        ref_line_x <- calculate_reference_point(
+          dat = if (is.data.frame(dat)) { dat } else { dat[[i]] },
+          reference_name = glue::glue("{label}_{ref_line[i]}"),
+          lbs = lbs
+        ) / scale_amount
+        if (length(ref_line_x) == 0) {
+          cli::cli_alert_warning("{label}_{ref_line[i]} not found for {ifelse(is.data.frame(dat), 'data', names(dat[i]))}")
+          next
+        }
+        ref_line_x <- stats::setNames(ref_line_x, ref_line[i])
+      } else {
+        ref_line_x <- ref_line[i] / scale_amount
+      }
+      
+      if ("unfished" %in% names(ref_line_x)) {
+        if (is.data.frame(dat)) {
+          sel_dat <- dat
+        } else {
+          sel_dat <- dat[[i]]
+        }
+        plt_lab <- label
+        # find the minimum x axis value from the plot
+        min_year <- sel_dat |> 
+          dplyr::filter(grepl(plt_lab, label), year != 1) |>
+          dplyr::pull(year) |>
+          min(na.rm = TRUE) |>
+          round(digits = 2)
+        # add point to plot and add theme
+        # plt2 <- plt2 +
+        # TODO: set color for each point to match that of the line for the model
+        ref_lines_list <- append(ref_lines_list,
+          ggplot2::geom_point(ggplot2::aes(x = min_year - 1, y = ref_line_x, color = model)) # should I keep -1 or set as first year?
+        )
+      } else {
+        # add apply/purrr/or for loop for reference lines -- not just the first anymore
+        # plt2 <- plt2 +
+        ref_lines_list <- append(ref_lines_list,
+          reference_line(
+            # conditionally add label name
+            label_name = ifelse(length(names(dat)[i]) == 1, label, names(dat)[i]), #"spawning_biomass",
+            ref_line = ref_line_x,
+            scale_amount = scale_amount,
+            model_name = ifelse(is.data.frame(dat), "1", names(dat)[i])
+          )
+        )
+      }
+    } # close ref_line for loop
+  } # close if ref_line NULL
+  ref_lines_list
+}
diff --git a/R/utils_rda.R b/R/utils_rda.R
index a5b84d2c..27721710 100644
--- a/R/utils_rda.R
+++ b/R/utils_rda.R
@@ -177,7 +177,7 @@ calc_kqs <- function(returned_kq,
   if (returned_kq == "F.min") {
     return(
       prepared_data$estimate |>
-        na.omit() |>
+        data.table::na.omit() |>
         min() |>
         round(digits = 3)
     )
@@ -186,7 +186,7 @@ calc_kqs <- function(returned_kq,
   if (returned_kq == "F.max") {
     return(
       prepared_data$estimate |>
-        na.omit() |>
+        data.table::na.omit() |>
         max() |>
         round(digits = 3)
     )
diff --git a/man/reference_line.Rd b/man/reference_line.Rd
index fbd71ec4..72ac2bd7 100644
--- a/man/reference_line.Rd
+++ b/man/reference_line.Rd
@@ -4,29 +4,19 @@
 \alias{reference_line}
 \title{Preformatted reference line}
 \usage{
-reference_line(plot, dat, label_name, reference, scale_amount = 1, lbs = FALSE)
+reference_line(label_name, ref_line, scale_amount = 1)
 }
 \arguments{
-\item{plot}{Plot object. A ggplot2 object where the reference line will be added}
-
-\item{dat}{Data frame. Standard data frame where reference point should be extracted}
-
 \item{label_name}{String. Name of the quantity that users want to
 extract the reference point from}
 
-\item{reference}{String of the reference point
+\item{ref_line}{String of the reference point
 
 Options: Including, but not limited to: "msy", "unfished", "target"}
 
 \item{scale_amount}{Number. A number to scale the y-axis values.
 
 Default: 1}
-
-\item{lbs}{Logical. TRUE/FALSE; indicate whether to convert the y-axis values from
-kilograms to pounds. The default units match the default in the
-unit_label argument - 'mt'.
-
-Default: `FALSE`}
 }
 \value{
 A ggplot2 geom_hline object for a reference point that can be added