diff --git a/DESCRIPTION b/DESCRIPTION
index aca0384e..192f2448 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: ss3sim
 Title: Fisheries Stock Assessment Simulation Testing with Stock Synthesis
-Version: 1.20.2
+Version: 1.21.0
 Authors@R: c(
     person(c("Kelli", "F."), "Johnson", , "kelli.johnson@noaa.gov", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-5149-451X")),
diff --git a/NEWS.md b/NEWS.md
index 4d93a089..401fc062 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,14 @@
 # News
 
+## ss3sim 1.21.0
+
+* Added season to the time series results and stopped hard-coding to
+  season 1.
+* Added month to the time series results.
+* Use better semantic versioning protocols where the minor number
+  is reserved for features and the patch number (terminal) is for
+  fixes.
+
 ## ss3sim 1.20.2
 
 * Fixed bug in Introduction vignette regarding self-test estimation method to
diff --git a/R/get-results.r b/R/get-results.r
index f518680c..ed4609a3 100644
--- a/R/get-results.r
+++ b/R/get-results.r
@@ -382,15 +382,40 @@ get_results_mod <- function(dir = getwd(), is_EM = NULL, is_OM = NULL) {
   )
 }
 
-#' Extract time series from a model run.
+#' Return the time series information from an iteration
 #'
-#' Extract time series from an [r4ss::SS_output()] list from a model run.
-#' Returns a data.frame of the results for SSB, recruitment and effort by year.
+#' Extract and return time series from an [r4ss::SS_output()] list, that is
+#' read in from the estimation method of a single iteration. The main time
+#' series information is included but no information about the uncertainty of
+#' those measurements is available. See the derived quantities for uncertainty.
+#'
+#' @details
+#' Information about both season and area are included in the data frame. For
+#' values that have no associated season or area, i.e., are summary values over
+#' all areas and seasons, the values are repeated for each area/season
+#' combination within a given year. For example, the recruitment deviation is
+#' for all areas and is thus repeated in each row across areas for a given year.
 #'
 #' @template report.file
 #' @export
 #' @family get-results
 #' @author Cole Monnahan
+#' @return
+#' A data frame with the following columns:
+#' * year
+#' * Area
+#' * Seas
+#' * Bio_smry
+#' * SpawnBio
+#' * Recruit_0
+#' * retainB_[0-9]+
+#' * retainN_[0-9]+
+#' * deadB_[0-9]+
+#' * deadN_[0-9]+
+#' * F_[0-9]+
+#' * SPRratio
+#' * rec_dev
+#' * raw_rec_dev
 get_results_timeseries <- function(report.file) {
   years <- report.file$startyr:(report.file$endyr +
     ifelse(is.na(report.file$nforecastyears),
@@ -401,7 +426,7 @@ get_results_timeseries <- function(report.file) {
   catch_cols <- grep("^retain\\([B|N]\\):_", colnames(report.file$timeseries))
   dead_cols <- grep("^dead\\([B|N]\\):_", colnames(report.file$timeseries))
   other_cols <- which(colnames(report.file$timeseries) %in%
-    c("Yr", "Seas", "Bio_smry", "SpawnBio", "Recruit_0"))
+    c("Yr", "Area", "Seas", "Bio_smry", "SpawnBio", "Recruit_0"))
   xx <- report.file$timeseries[, c(other_cols, catch_cols, dead_cols, F_cols)]
   # remove parentheses and colons from column names
   colnames(xx) <- gsub("\\(|\\)|\\:", "", colnames(xx))
@@ -421,8 +446,13 @@ get_results_timeseries <- function(report.file) {
       spr[, grep("label", colnames(spr), ignore.case = TRUE)], "_"
     ), "[", 2))
     colnames(spr)[which(colnames(spr) == "Value")] <- "SPRratio"
-    spr[["Seas"]] <- 1 # need to add seasonal column; just assign to first? Or should be NA?
-    df <- merge(xx, spr[, c("SPRratio", "Yr", "Seas")], by = c("Yr", "Seas"), all.x = TRUE)
+    df <- merge(
+      xx,
+      spr[, c("SPRratio", "Yr")],
+      by = c("Yr"),
+      all.x = TRUE,
+      all.y = FALSE
+    )
     df$SPRratio[is.na(df$SPRratio)] <- 0
   } else {
     df <- xx
@@ -436,17 +466,16 @@ get_results_timeseries <- function(report.file) {
   )
   dev <- dev[dev[, getcols[1]] %in% years, getcols]
   colnames(dev) <- gsub("dev", "rec_dev", colnames(dev), ignore.case = TRUE)
-  dev[["Seas"]] <- 1 # Add Seas; just assign to 1? or should be NA?
-  ## create final data.frame
-  df <- merge(df, dev,
-    by.x = c("Yr", "Seas"),
-    by.y = c(colnames(dev)[getcols[1]], "Seas"), all.x = TRUE, all.y = TRUE
-  )
-  rownames(df) <- NULL
-  # change year name
-  df$year <- df$Yr
-  df$Yr <- NULL
-  df
+  dev[["Seas"]] <- report.file[["spawnseas"]]
+  # create final data.frame
+  out <- dplyr::full_join(
+    x = df |> dplyr::mutate(Yr = as.integer(Yr)),
+    y = dev,
+    by = c(Yr = colnames(dev)[getcols[1]], "Seas")
+  ) |>
+    dplyr::arrange(Yr, Area, Seas) |>
+    dplyr::rename(year = Yr)
+  return(out)
 }
 
 #' Extract time series from a model run with the associated standard deviation.
diff --git a/man/get_results_timeseries.Rd b/man/get_results_timeseries.Rd
index c7287ec5..eb931baf 100644
--- a/man/get_results_timeseries.Rd
+++ b/man/get_results_timeseries.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/get-results.r
 \name{get_results_timeseries}
 \alias{get_results_timeseries}
-\title{Extract time series from a model run.}
+\title{Return the time series information from an iteration}
 \usage{
 get_results_timeseries(report.file)
 }
@@ -10,9 +10,37 @@ get_results_timeseries(report.file)
 \item{report.file}{An \code{\link[r4ss]{SS_output}} list for a model
 (operating model or estimation model).}
 }
+\value{
+A data frame with the following columns:
+\itemize{
+\item year
+\item Area
+\item Seas
+\item Bio_smry
+\item SpawnBio
+\item Recruit_0
+\item retainB_\link{0-9}+
+\item retainN_\link{0-9}+
+\item deadB_\link{0-9}+
+\item deadN_\link{0-9}+
+\item F_\link{0-9}+
+\item SPRratio
+\item rec_dev
+\item raw_rec_dev
+}
+}
 \description{
-Extract time series from an \code{\link[r4ss:SS_output]{r4ss::SS_output()}} list from a model run.
-Returns a data.frame of the results for SSB, recruitment and effort by year.
+Extract and return time series from an \code{\link[r4ss:SS_output]{r4ss::SS_output()}} list, that is
+read in from the estimation method of a single iteration. The main time
+series information is included but no information about the uncertainty of
+those measurements is available. See the derived quantities for uncertainty.
+}
+\details{
+Information about both season and area are included in the data frame. For
+values that have no associated season or area, i.e., are summary values over
+all areas and seasons, the values are repeated for each area/season
+combination within a given year. For example, the recruitment deviation is
+for all areas and is thus repeated in each row across areas for a given year.
 }
 \seealso{
 Other get-results: