diff --git a/DESCRIPTION b/DESCRIPTION
index adb85ad..872f0c9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -9,7 +9,7 @@ URL: https://scfmutils.predictiveecology.org,
     https://github.com/PredictiveEcology/scfmutils,
     https://predictiveecology.github.io/scfmutils/
 Date: 2025-03-12
-Version: 2.0.9.9003
+Version: 2.0.9.9004
 Authors@R: c(
     person("Steve", "Cumming", email = "Steve.Cumming@sbf.ulaval.ca",
            role = c("aut")),
diff --git a/NEWS.md b/NEWS.md
index 3ceb22b..2ba8d1b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,7 @@
 - streamline `prepInputsFireRegimePolys()` for BEC;
 - `prepInputsFireRegimePolys()` gets a new argument `subsetType` to retrieve entire polygons
   bordering `studyArea` or only the spatial intersection (default, preserving previous behaviour) 
+- add a title parameter for `scfmDiagnostics` plots 
 
 # scfmutils 2.0.8
 
diff --git a/R/comparePredictions.R b/R/comparePredictions.R
index 163e4cc..e3b81c9 100644
--- a/R/comparePredictions.R
+++ b/R/comparePredictions.R
@@ -21,17 +21,19 @@ utils::globalVariables(c(
 #' @examples
 #' \dontrun{
 #' ## assumes user has run scfm to produce the simList `mySimOut`
-#' dt <- comparePredictions_summaryDT(fireRegimePoints = mySimOut$fireRegimePoints,
-#'                                    burnSummary = mySimOut$burnSummary,
-#'                                    fireRegimePolys = mySimOut$fireRegimePolys,
-#'                                    times = times(mySimOut))
+#' dt <- comparePredictions_summaryDT(
+#'   fireRegimePoints = mySimOut$fireRegimePoints,
+#'   burnSummary = mySimOut$burnSummary,
+#'   fireRegimePolys = mySimOut$fireRegimePolys,
+#'   times = times(mySimOut)
+#' )
 #'
 #' gg_mfs <- comparePredictions_meanFireSize(dt)
 #' gg_fri <- comparePredictions_fireReturnInterval(dt)
 #' gg_ign <- comparePredictions_annualIgnitions(dt)
 #' gg_frp <- plot_fireRegimePolys(mySimOut$fireRegimePolys)
 #'
-#' gridExtra::grid.arrange(fps,  gg_mfs,  gg_fri,  gg_ign, nrow = 2, ncol = 2)
+#' gridExtra::grid.arrange(fps, gg_mfs, gg_fri, gg_ign, nrow = 2, ncol = 2)
 #' }
 #'
 #' @author Ian Eddy
@@ -43,9 +45,11 @@ utils::globalVariables(c(
 #' @rdname comparePredictions
 comparePredictions_summaryDT <- function(fireRegimePoints = NULL, burnSummary = NULL,
                                          fireRegimePolys = NULL, times = NULL) {
-  if (any(is.null(fireRegimePolys$pSpread), is.null(fireRegimePolys$xBar),
-          is.null(fireRegimePolys$burnyArea), is.null(fireRegimePoints),
-          is.null(burnSummary), is.null(times))) {
+  if (any(
+    is.null(fireRegimePolys$pSpread), is.null(fireRegimePolys$xBar),
+    is.null(fireRegimePolys$burnyArea), is.null(fireRegimePoints),
+    is.null(burnSummary), is.null(times)
+  )) {
     stop("fireRegimePolys is missing columns or insufficient args provided")
   }
 
@@ -62,14 +66,16 @@ comparePredictions_summaryDT <- function(fireRegimePoints = NULL, burnSummary =
     ## median fire size is not used by scfm but is worth recording
     ## regimes where mean is much greater than median will be hard to recreate
     escaped <- fireRegimePoints[fireRegimePoints$SIZE_HA > fireRegimePoly$cellSize, ]
-    medianFireSize <- median(escaped$SIZE_HA) #should be no need for na.rm
+    medianFireSize <- median(escaped$SIZE_HA) # should be no need for na.rm
 
     pSpread <- fireRegimePoly$pSpread
     pIg <- fireRegimePoly$ignitionRate
 
     if (!"grp" %in% names(burnSummary)) {
-      stop("burnSummary data.table does not have a 'grp' column.\n",
-           "Are you running a recent version of scfmSpread (>= 2.0.0)?")
+      stop(
+        "burnSummary data.table does not have a 'grp' column.\n",
+        "Are you running a recent version of scfmSpread (>= 2.0.0)?"
+      )
     }
 
     ## burnSummary data.table from scfmSpread is coded as follows:
@@ -86,7 +92,7 @@ comparePredictions_summaryDT <- function(fireRegimePoints = NULL, burnSummary =
     achievedEscapes <- nrow(burnSum[grp %in% 1 & N > 1]) / simLength
 
     ## mean fire size: mean size of all fires ignited and escaped in SAR, regardless of where spread
-    burnSum1 <- burnSum[grp %in% c(1, 2), lapply(.SD, sum), by = c("igLoc", "year"), .SDcols = "areaBurned"]  # nolint
+    burnSum1 <- burnSum[grp %in% c(1, 2), lapply(.SD, sum), by = c("igLoc", "year"), .SDcols = "areaBurned"] # nolint
     burnSum1 <- burnSum1[areaBurned %>>% fireRegimePoly$cellSize, ]
     meanFireSize <- ifelse(nrow(burnSum1) == 0, 0, mean(burnSum1$areaBurned))
 
@@ -96,23 +102,25 @@ comparePredictions_summaryDT <- function(fireRegimePoints = NULL, burnSummary =
     MAAB <- sum(burnSum2$areaBurned) / simLength
 
     achievedFRI <- simLength / (sum(0, burnSum2$areaBurned) / fireRegimePoly$burnyArea)
-    targetFRI  <- 1 / fireRegimePoly$empiricalBurnRate
-
-    pred <- data.frame("PolyID" = x,
-                       "histMeanSize" = fireRegimePoly$xBar, ## predicted (empirical) mean fire size
-                       "histMedianSize" = medianFireSize,
-                       "modMeanSize" = meanFireSize,
-                       "achievedFRI" = achievedFRI,
-                       "targetFRI" = targetFRI,
-                       "burnableArea_ha" = fireRegimePoly$burnyArea,
-                       "targetIgnitions" = targetIgnitions,
-                       "achievedIgnitions" = achievedIgnitions,
-                       "targetEscapes" = targetEscapes,
-                       "achievedEscapes" = achievedEscapes,
-                       "pEscape" = fireRegimePoly$pEscape, ## escape prob (no. fires > cellSize / no. fires)
-                       "p0" = fireRegimePoly$p0,  ## p0 and pEscape may indicate something incorrect
-                       "pSpread" = pSpread, ## spread probability estimated from the SCAM model
-                       "pIgnition" = pIg) ## ignition probability of a single pixel
+    targetFRI <- 1 / fireRegimePoly$empiricalBurnRate
+
+    pred <- data.frame(
+      "PolyID" = x,
+      "histMeanSize" = fireRegimePoly$xBar, ## predicted (empirical) mean fire size
+      "histMedianSize" = medianFireSize,
+      "modMeanSize" = meanFireSize,
+      "achievedFRI" = achievedFRI,
+      "targetFRI" = targetFRI,
+      "burnableArea_ha" = fireRegimePoly$burnyArea,
+      "targetIgnitions" = targetIgnitions,
+      "achievedIgnitions" = achievedIgnitions,
+      "targetEscapes" = targetEscapes,
+      "achievedEscapes" = achievedEscapes,
+      "pEscape" = fireRegimePoly$pEscape, ## escape prob (no. fires > cellSize / no. fires)
+      "p0" = fireRegimePoly$p0, ## p0 and pEscape may indicate something incorrect
+      "pSpread" = pSpread, ## spread probability estimated from the SCAM model
+      "pIgnition" = pIg
+    ) ## ignition probability of a single pixel
     return(pred)
   })
   return(rbindlist(out))
@@ -120,16 +128,17 @@ comparePredictions_summaryDT <- function(fireRegimePoints = NULL, burnSummary =
 
 #' @param dt scfm summary `data.table` produced by `comparePredictions_summaryDT()`
 #'
+#' @param title character, the plot title
+#'
 #' @export
 #' @importFrom ggplot2 aes geom_abline geom_point geom_text ggplot labs
 #' @importFrom ggplot2 scale_x_continuous scale_y_continuous theme_bw xlab ylab
 #' @rdname comparePredictions
-comparePredictions_meanFireSize <- function(dt) {
+comparePredictions_meanFireSize <- function(dt, title = NULL) {
   if (any(is.null(dt))) {
     stop("all arguments must be provided and cannot be NULL.")
   }
-
-  ggplot(dt, aes(x = histMeanSize, y = modMeanSize)) +
+  gg_mfs <- ggplot(dt, aes(x = histMeanSize, y = modMeanSize)) +
     geom_point(aes(histMeanSize, modMeanSize)) +
     labs(x = "historical mean fire size (ha)", y = "modeled mean fire size (ha)") +
     theme_bw() +
@@ -137,13 +146,17 @@ comparePredictions_meanFireSize <- function(dt) {
     scale_x_continuous(limits = c(0, NA)) +
     geom_text(aes(label = PolyID), vjust = "inward", hjust = "inward") +
     geom_abline(slope = 1)
+  if (!is.null(title)) {
+    gg_mfs <- gg_mfs + ggtitle(title)
+  }
+  return(gg_mfs)
 }
 
 #' @export
 #' @rdname comparePredictions
-comparePredictions_fireReturnInterval <- function(dt, times) {
+comparePredictions_fireReturnInterval <- function(dt, times, title = NULL) {
   if (any(is.null(dt), is.null(times))) {
-    stop("all arguments must be provided and cannot be NULL.")
+    stop("dt and times must be provided and cannot be NULL.")
   }
 
   ## remove the infinite FRI caused by no simulated fires
@@ -153,10 +166,11 @@ comparePredictions_fireReturnInterval <- function(dt, times) {
     ## TODO: confirm wording of this, with comment above
     warning("achievedFRI may be off where targetFRI is less than 4x the simulated time.")
   }
-
   ## TODO: remove targetFRI filter below. plot those points differently to indicate poor estimates
-  ggplot(dt[!is.infinite(achievedFRI) & targetFRI < c(times$end - times$start) * 4],
-         aes(x = targetFRI, y = achievedFRI)) +
+  gg_fri <- ggplot(
+    dt[!is.infinite(achievedFRI) & targetFRI < c(times$end - times$start) * 4],
+    aes(x = targetFRI, y = achievedFRI)
+  ) +
     geom_point() +
     labs(y = "simulation FRI (years)", x = "estimated FRI (years)") +
     theme_bw() +
@@ -164,21 +178,24 @@ comparePredictions_fireReturnInterval <- function(dt, times) {
     scale_y_continuous(limits = c(0, NA)) +
     scale_x_continuous(limits = c(0, NA)) +
     geom_text(aes(label = PolyID, vjust = "inward", hjust = "inward"))
+  if (!is.null(title)) {
+    gg_fri <- gg_fri + ggtitle(title)
+  }
+  return(gg_fri)
 }
 
 #' @export
 #' @rdname comparePredictions
-comparePredictions_annualIgnitions <- function(dt) {
+comparePredictions_annualIgnitions <- function(dt, title = NULL) {
   if (any(is.null(dt))) {
     stop("all arguments must be provided and cannot be NULL.")
   }
 
-  dt <- copy(dt) #avoid adding per ha cols (or add them?)
+  dt <- copy(dt) # avoid adding per ha cols (or add them?)
   dt[, targetIgnitions_Mha := targetIgnitions / burnableArea_ha * 1e6]
   dt[, achievedIgnitions_Mha := achievedIgnitions / burnableArea_ha * 1e6]
 
-
-  ggplot(dt, aes(x = targetIgnitions_Mha, y = achievedIgnitions_Mha)) +
+  gg_ai <- ggplot(dt, aes(x = targetIgnitions_Mha, y = achievedIgnitions_Mha)) +
     geom_point() +
     labs(y = "simulation annual ignitions (per Mha)", x = "estimated annual ignitions (per Mha)") +
     theme_bw() +
@@ -186,21 +203,24 @@ comparePredictions_annualIgnitions <- function(dt) {
     scale_y_continuous(limits = c(0, NA)) +
     scale_x_continuous(limits = c(0, NA)) +
     geom_text(aes(label = PolyID, vjust = "inward", hjust = "inward"))
+  if (!is.null(title)) {
+    gg_ai <- gg_ai + ggtitle(title)
+  }
+  return(gg_ai)
 }
 
 #' @export
 #' @rdname comparePredictions
-comparePredictions_annualEscapes <- function(dt) {
+comparePredictions_annualEscapes <- function(dt, title = NULL) {
   if (any(is.null(dt))) {
     stop("all arguments must be provided and cannot be NULL.")
   }
 
-  dt <- copy(dt) #avoid adding per ha cols (or add them?)
+  dt <- copy(dt) # avoid adding per ha cols (or add them?)
   dt[, targetEscapes_Mha := targetEscapes / burnableArea_ha * 1e6]
   dt[, achievedEscapes_Mha := achievedEscapes / burnableArea_ha * 1e6]
 
-
-  ggplot(dt, aes(x = targetEscapes_Mha, y = achievedEscapes_Mha)) +
+  gg_ae <- ggplot(dt, aes(x = targetEscapes_Mha, y = achievedEscapes_Mha)) +
     geom_point() +
     labs(y = "simulation annual escapes (per Mha)", x = "estimated annual escapes (per Mha)") +
     theme_bw() +
@@ -208,18 +228,24 @@ comparePredictions_annualEscapes <- function(dt) {
     scale_y_continuous(limits = c(0, NA)) +
     scale_x_continuous(limits = c(0, NA)) +
     geom_text(aes(label = PolyID, vjust = "inward", hjust = "inward"))
+  if (!is.null(title)) {
+    gg_ae <- gg_ae + ggtitle(title)
+  }
+  return(gg_ae)
 }
 
 #' @param size minimum fire size (ha)
 #'
+#' @param title character, the plot title
+#'
 #' @export
 #' @importFrom data.table as.data.table setnames
 #' @importFrom fpCompare %>>%
 #' @importFrom ggplot2 facet_wrap geom_histogram
 #' @rdname comparePredictions
-comparePredictions_fireDistribution <- function(fireRegimePoints = NULL, burnSummary = NULL, size) {
+comparePredictions_fireDistribution <- function(fireRegimePoints = NULL, burnSummary = NULL, size, title = NULL) {
   histDistribution <- fireRegimePoints[fireRegimePoints$SIZE_HA %>>% size, ]
-  if (nrow(histDistribution) < 1){
+  if (nrow(histDistribution) < 1) {
     warning("no historical fires > escape size - showing all fires instead")
     histDistribution <- fireRegimePoints
   }
@@ -235,11 +261,15 @@ comparePredictions_fireDistribution <- function(fireRegimePoints = NULL, burnSum
   simDistribution[, source := "simulated"]
 
   allFires <- rbind(simDistribution, histDistribution, fill = TRUE)
-  allFires[ , PolyID := as.factor(PolyID)]
+  allFires[, PolyID := as.factor(PolyID)]
 
-  ggplot(allFires, aes(x = log(areaBurned), fill = PolyID)) +
+  gg_fd <- ggplot(allFires, aes(x = log(areaBurned), fill = PolyID)) +
     geom_histogram() +
     xlab("log of escaped fire size (ha)") +
     theme_bw() +
     facet_wrap(~source, scales = "free_y", nrow = 2)
+  if (!is.null(title)) {
+    gg_fd <- gg_fd + ggtitle(title)
+  }
+  return(gg_fd)
 }
diff --git a/R/plot_wrappers.R b/R/plot_wrappers.R
index 6372b2f..a6b77d0 100644
--- a/R/plot_wrappers.R
+++ b/R/plot_wrappers.R
@@ -6,19 +6,25 @@ utils::globalVariables(c(
 #'
 #' @template fireRegimePolys
 #'
+#' @param title character, the plot title
+#'
 #' @returns a `ggplot` object
 #'
 #' @export
 #' @importFrom ggplot2 aes geom_sf ggplot scale_fill_discrete theme_minimal
-plot_fireRegimePolys <- function(fireRegimePolys) {
+plot_fireRegimePolys <- function(fireRegimePolys, title = NULL) {
   if (!is.factor(fireRegimePolys$PolyID)) {
     fireRegimePolys$PolyID <- as.factor(fireRegimePolys$PolyID)
   }
-
-  ggplot(fireRegimePolys) +
+  gg_frp <- ggplot(fireRegimePolys) +
     geom_sf(aes(fill = PolyID)) +
     scale_fill_discrete() + ## TODO: use same palette as plot_fireRegimeRas ??
     theme_bw()
+
+  if (!is.null(title)) {
+    gg_frp <- gg_frp + ggtitle(title)
+  }
+  return(gg_frp)
 }
 
 #' Plot fire regime raster
@@ -33,11 +39,15 @@ plot_fireRegimePolys <- function(fireRegimePolys) {
 #' @importFrom ggplot2 ggplot ggtitle scale_fill_brewer theme_bw
 #' @importFrom tidyterra geom_spatraster
 plot_fireRegimeRas <- function(x, title) {
-  ggplot() +
+  gg_frr <- ggplot() +
     geom_spatraster(data = terra::as.factor(x)) +
     scale_fill_brewer(palette = "Paired", type = "qual", na.value = "transparent") +
-    ggtitle(title) +
     theme_bw()
+
+  if (!is.null(title)) {
+    gg_frr <- gg_frr + ggtitle(title)
+  }
+  return(gg_frr)
 }
 
 #' Plot age map
@@ -56,11 +66,15 @@ plot_fireRegimeRas <- function(x, title) {
 plot_ageMap <- function(x, title, maxAge) {
   x[x > maxAge] <- maxAge
 
-  ggplot() +
+  gg_am <- ggplot() +
     geom_spatraster(data = x) +
     scale_fill_distiller(palette = "Greens", direction = 1, na.value = "transparent") +
-    ggtitle(title) +
     theme_bw()
+
+  if (!is.null(title)) {
+    gg_am <- gg_am + ggtitle(title)
+  }
+  return(gg_am)
 }
 
 #' Plot burn maps
@@ -69,18 +83,23 @@ plot_ageMap <- function(x, title, maxAge) {
 #'
 #' @param title character, the plot title
 #'
+#' @param subtitle character, the plot subtitle
+#'
 #' @returns `ggplot` object
 #'
 #' @export
 #' @importFrom ggplot2 ggplot ggtitle theme_bw
 #' @importFrom tidyterra geom_spatraster
 #' @importFrom viridis scale_fill_viridis
-plot_burnMap <- function(x, title) {
-  ggplot() +
+plot_burnMap <- function(x, title = NULL, subtitle = NULL) {
+  gg_bm <- ggplot() +
     geom_spatraster(data = x) +
     scale_fill_viridis(na.value = "transparent") +
-    ggtitle(title) +
     theme_bw()
+
+    gg_bm <- gg_bm + ggtitle(title) + labs(subtitle = subtitle)
+
+  return(gg_bm)
 }
 
 #' Plot flammable map
@@ -95,9 +114,13 @@ plot_burnMap <- function(x, title) {
 #' @importFrom ggplot2 ggplot ggtitle scale_fill_distiller theme_bw
 #' @importFrom tidyterra geom_spatraster
 plot_flammableMap <- function(x, title) {
-  ggplot() +
+  gg_fm <- ggplot() +
     geom_spatraster(data = x) +
     scale_fill_distiller(palette = "RdBu", na.value = "transparent") +
-    ggtitle(title) +
     theme_bw()
-}
+  
+  if (!is.null(title)) {
+    gg_fm <- gg_fm + ggtitle(title)
+  }
+  return(gg_fm)
+}
\ No newline at end of file
diff --git a/man/comparePredictions.Rd b/man/comparePredictions.Rd
index 784f8fc..4dae0d1 100644
--- a/man/comparePredictions.Rd
+++ b/man/comparePredictions.Rd
@@ -16,18 +16,19 @@ comparePredictions_summaryDT(
   times = NULL
 )
 
-comparePredictions_meanFireSize(dt)
+comparePredictions_meanFireSize(dt, title = NULL)
 
-comparePredictions_fireReturnInterval(dt, times)
+comparePredictions_fireReturnInterval(dt, times, title = NULL)
 
-comparePredictions_annualIgnitions(dt)
+comparePredictions_annualIgnitions(dt, title = NULL)
 
-comparePredictions_annualEscapes(dt)
+comparePredictions_annualEscapes(dt, title = NULL)
 
 comparePredictions_fireDistribution(
   fireRegimePoints = NULL,
   burnSummary = NULL,
-  size
+  size,
+  title = NULL
 )
 }
 \arguments{
@@ -41,6 +42,8 @@ comparePredictions_fireDistribution(
 
 \item{dt}{scfm summary \code{data.table} produced by \code{comparePredictions_summaryDT()}}
 
+\item{title}{character, the plot title}
+
 \item{size}{minimum fire size (ha)}
 }
 \value{
@@ -53,17 +56,19 @@ Create \code{data.table} to compare scfm predictions with historical observation
 \examples{
 \dontrun{
 ## assumes user has run scfm to produce the simList `mySimOut`
-dt <- comparePredictions_summaryDT(fireRegimePoints = mySimOut$fireRegimePoints,
-                                   burnSummary = mySimOut$burnSummary,
-                                   fireRegimePolys = mySimOut$fireRegimePolys,
-                                   times = times(mySimOut))
+dt <- comparePredictions_summaryDT(
+  fireRegimePoints = mySimOut$fireRegimePoints,
+  burnSummary = mySimOut$burnSummary,
+  fireRegimePolys = mySimOut$fireRegimePolys,
+  times = times(mySimOut)
+)
 
 gg_mfs <- comparePredictions_meanFireSize(dt)
 gg_fri <- comparePredictions_fireReturnInterval(dt)
 gg_ign <- comparePredictions_annualIgnitions(dt)
 gg_frp <- plot_fireRegimePolys(mySimOut$fireRegimePolys)
 
-gridExtra::grid.arrange(fps,  gg_mfs,  gg_fri,  gg_ign, nrow = 2, ncol = 2)
+gridExtra::grid.arrange(fps, gg_mfs, gg_fri, gg_ign, nrow = 2, ncol = 2)
 }
 
 }
diff --git a/man/plot_burnMap.Rd b/man/plot_burnMap.Rd
index d07f7a8..2fd1744 100644
--- a/man/plot_burnMap.Rd
+++ b/man/plot_burnMap.Rd
@@ -4,12 +4,14 @@
 \alias{plot_burnMap}
 \title{Plot burn maps}
 \usage{
-plot_burnMap(x, title)
+plot_burnMap(x, title = NULL, subtitle = NULL)
 }
 \arguments{
 \item{x}{\code{SpatRaster} object corresponding to a current or cumulative burn map.}
 
 \item{title}{character, the plot title}
+
+\item{subtitle}{character, the plot subtitle}
 }
 \value{
 \code{ggplot} object
diff --git a/man/plot_fireRegimePolys.Rd b/man/plot_fireRegimePolys.Rd
index 1e53826..5392428 100644
--- a/man/plot_fireRegimePolys.Rd
+++ b/man/plot_fireRegimePolys.Rd
@@ -4,10 +4,12 @@
 \alias{plot_fireRegimePolys}
 \title{Plot fire regime polygons}
 \usage{
-plot_fireRegimePolys(fireRegimePolys)
+plot_fireRegimePolys(fireRegimePolys, title = NULL)
 }
 \arguments{
 \item{fireRegimePolys}{\code{sf} polygon or multipolygon object defining the fire regime polygons}
+
+\item{title}{character, the plot title}
 }
 \value{
 a \code{ggplot} object