diff --git a/NAMESPACE b/NAMESPACE
index 3836bbb8c..2a6d95c3c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -48,6 +48,7 @@ export(extract_CrIs)
 export(extract_inits)
 export(extract_stan_param)
 export(fix_dist)
+export(forecast_infections)
 export(forecast_secondary)
 export(gamma_dist_def)
 export(generation_time_opts)
diff --git a/NEWS.md b/NEWS.md
index d6c62f148..c006184ca 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -6,6 +6,7 @@
 * The utility function `update_list()` has been deprecated in favour of `utils::modifyList()` because it comes with an installation of R. By @jamesmbaazam in #491 and reviewed by @seabbs.
 * The `fixed` argument to `dist_spec` has been deprecated and replaced by a `fix_dist()` function. By @sbfnk in #503 and reviewed by @seabbs.
 * Updated `estimate_infections()` so that rather than imputing missing data, it now skips these data points in the likelihood. This is a breaking change as it alters the behaviour of the model when dates are missing from a time series but are known to be zero. We recommend that users check their results when updating to this version but expect this to in most cases improve performance. By @seabbs in #528 and reviewed by @sbfnk.
+* `simulate_infections` has been renamed to `forecast_infections` in line with `simulate_secondary` and `forecast_secondary`. The terminology is: a forecast is done from a fit to existing data, a simulation from first principles. By @sbfnk in #544 and reviewed by @seabbs.
 
 ## Documentation
 
@@ -168,8 +169,8 @@ reporting delay and the generation time. These are based on an implementation in
  
 ## Deprecated features
 
-* `simulate_cases()` and `forecast_infections()` have been deprecated and have been removed. These functions depend on `EpiSoon` which itself is archived and near equivalent functionality is available within `EpiNow2` and in other packages (@seabbs).
-* Functions supporting secondary forecasting using `forecast_infections()` (i.e in `epinow())
+* `simulate_cases()` and `simulate_infections()` have been deprecated and have been removed. These functions depend on `EpiSoon` which itself is archived and near equivalent functionality is available within `EpiNow2` and in other packages (@seabbs).
+* Functions supporting secondary forecasting using `simulate_infections()` (i.e in `epinow())
 have been removed along with the arguments that supported them (@seabbs).
 * `global_map()`, `country_map()`, and `theme_map()` have all been deprecated and have been removed. These functions were used to support reporting of reproduction number
 estimates and are considered out of scope for `EpiNow2`. If finding useful contacting the
diff --git a/R/epinow.R b/R/epinow.R
index 5eefe62e9..cd2156db1 100644
--- a/R/epinow.R
+++ b/R/epinow.R
@@ -24,7 +24,7 @@
 #'   summarising results and reporting errors if they have occurred.
 #' @author Sam Abbott
 #' @export
-#' @seealso [estimate_infections()] [simulate_infections()] [regional_epinow()]
+#' @seealso [estimate_infections()] [forecast_infections()] [regional_epinow()]
 #' @inheritParams setup_target_folder
 #' @inheritParams estimate_infections
 #' @inheritParams setup_default_logging
diff --git a/R/estimate_infections.R b/R/estimate_infections.R
index 63a945919..8ec103d19 100644
--- a/R/estimate_infections.R
+++ b/R/estimate_infections.R
@@ -51,7 +51,7 @@
 #' samples, data used to fit the model, and the fit object itself.
 #'
 #' @author Sam Abbott
-#' @seealso [epinow()] [regional_epinow()] [simulate_infections()]
+#' @seealso [epinow()] [regional_epinow()] [forecast_infections()]
 #' [estimate_truncation()]
 #' @inheritParams create_stan_args
 #' @inheritParams create_stan_data
diff --git a/R/simulate_infections.R b/R/simulate_infections.R
index 6ac8b1cb5..af917fd66 100644
--- a/R/simulate_infections.R
+++ b/R/simulate_infections.R
@@ -1,4 +1,22 @@
-#' Simulate infections using a given trajectory of the time-varying
+#' Deprecated; use [forecast_infections()] instead
+#'
+#' Calling this function passes all arguments to [forecast_infections()]
+#' @description `r lifecycle::badge("deprecated")`
+#' @param ... Arguments to be passed to [forecast_infections()]
+#' @return the result of [forecast_infections()]
+#' @export
+simulate_infections <- function(...) {
+    deprecate_warn(
+      "2.0.0",
+      "simulate_infections()",
+      "forecast_infections()",
+      "A new [simulate_infections()] function for simulating from given ",
+      "parameters is planned for implementation in the future."
+    )
+   forecast_infections(...)
+}
+
+#' Forecast infections from a given fit and trajectory of the time-varying
 #' reproduction number
 #'
 #' @description `r lifecycle::badge("stable")`
@@ -63,7 +81,7 @@
 #'
 #' # update Rt trajectory and simulate new infections using it
 #' R <- c(rep(NA_real_, 26), rep(0.5, 10), rep(0.8, 7))
-#' sims <- simulate_infections(est, R)
+#' sims <- forecast_infections(est, R)
 #' plot(sims)
 #'
 #' # with a data.frame input of samples
@@ -74,7 +92,7 @@
 #'   ),
 #'   value = R
 #' )
-#' sims <- simulate_infections(est, R_dt)
+#' sims <- forecast_infections(est, R_dt)
 #' plot(sims)
 #'
 #' #' # with a data.frame input of samples
@@ -83,12 +101,12 @@
 #'  .(date, sample, value)][sample <= 1000][date <= "2020-04-10"
 #' ]
 #' R_samples <- R_samples[date >= "2020-04-01", value := 1.1]
-#' sims <- simulate_infections(est, R_samples)
+#' sims <- forecast_infections(est, R_samples)
 #' plot(sims)
 #'
 #' options(old_opts)
 #' }
-simulate_infections <- function(estimates,
+forecast_infections <- function(estimates,
                                 R = NULL,
                                 model = NULL,
                                 samples = NULL,
diff --git a/R/utilities.R b/R/utilities.R
index 8dfde59a0..c253a5704 100644
--- a/R/utilities.R
+++ b/R/utilities.R
@@ -163,7 +163,7 @@ allocate_delays <- function(delay_var, no_delays) {
 #'
 #' @description `r lifecycle::badge("stable")`
 #' Allocate missing parameters to be empty two dimensional arrays. Used
-#' internally by [simulate_infections()].
+#' internally by [forecast_infections()].
 #' @param data A list of parameters
 #' @param params A character vector of parameters to allocate to
 #' empty if missing.
@@ -423,7 +423,7 @@ set_dt_single_thread <- function() {
 
 #' @importFrom stats glm median na.omit pexp pgamma plnorm quasipoisson rexp
 #' @importFrom lifecycle deprecate_warn
-#' @importFrom stats rlnorm rnorm rpois runif sd var rgamma
+#' @importFrom stats rlnorm rnorm rpois runif sd var rgamma pnorm
 globalVariables(
   c(
     "bottom", "cases", "confidence", "confirm", "country_code", "crps",
diff --git a/_pkgdown.yml b/_pkgdown.yml
index cc4ff0a49..a324e2be9 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -68,7 +68,7 @@ reference:
     desc: Function to estimate, simulate and forecast parameters of interest.
     contents:
       - estimate_infections
-      - simulate_infections
+      - forecast_infections
       - estimate_secondary
       - forecast_secondary
       - estimate_delay
@@ -126,7 +126,7 @@ reference:
   - title: Simulate
     desc: Functions to help with simulating data or mapping to reported cases
     contents:
-      - simulate_infections
+      - forecast_infections
       - simulate_secondary
       - adjust_infection_to_report
   - title: Data
diff --git a/inst/dev/recover-synthetic/rt.R b/inst/dev/recover-synthetic/rt.R
index a447debf0..b6dce1aa1 100644
--- a/inst/dev/recover-synthetic/rt.R
+++ b/inst/dev/recover-synthetic/rt.R
@@ -26,7 +26,7 @@ R <- c(
 )
 noisy_R <- R * rnorm(length(R), 1, 0.05)
 # update Rt trajectory and simulate new infections using it
-sims <- simulate_infections(init, R = noisy_R, samples = 10)
+sims <- forecast_infections(init, R = noisy_R, samples = 10)
 
 sim_R <- sims$summarised[variable == "R"]$median
 
diff --git a/man/allocate_empty.Rd b/man/allocate_empty.Rd
index 595d48f62..e269e2aac 100644
--- a/man/allocate_empty.Rd
+++ b/man/allocate_empty.Rd
@@ -20,5 +20,5 @@ A list of parameters some allocated to be empty
 \description{
 \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}}
 Allocate missing parameters to be empty two dimensional arrays. Used
-internally by \code{\link[=simulate_infections]{simulate_infections()}}.
+internally by \code{\link[=forecast_infections]{forecast_infections()}}.
 }
diff --git a/man/epinow.Rd b/man/epinow.Rd
index a4a302f8a..76cc6676f 100644
--- a/man/epinow.Rd
+++ b/man/epinow.Rd
@@ -176,7 +176,7 @@ options(old_opts)
 }
 }
 \seealso{
-\code{\link[=estimate_infections]{estimate_infections()}} \code{\link[=simulate_infections]{simulate_infections()}} \code{\link[=regional_epinow]{regional_epinow()}}
+\code{\link[=estimate_infections]{estimate_infections()}} \code{\link[=forecast_infections]{forecast_infections()}} \code{\link[=regional_epinow]{regional_epinow()}}
 }
 \author{
 Sam Abbott
diff --git a/man/estimate_infections.Rd b/man/estimate_infections.Rd
index 851e3a295..156a733d2 100644
--- a/man/estimate_infections.Rd
+++ b/man/estimate_infections.Rd
@@ -159,7 +159,7 @@ options(old_opts)
 }
 }
 \seealso{
-\code{\link[=epinow]{epinow()}} \code{\link[=regional_epinow]{regional_epinow()}} \code{\link[=simulate_infections]{simulate_infections()}}
+\code{\link[=epinow]{epinow()}} \code{\link[=regional_epinow]{regional_epinow()}} \code{\link[=forecast_infections]{forecast_infections()}}
 \code{\link[=estimate_truncation]{estimate_truncation()}}
 }
 \author{
diff --git a/man/forecast_infections.Rd b/man/forecast_infections.Rd
new file mode 100644
index 000000000..ec90e946c
--- /dev/null
+++ b/man/forecast_infections.Rd
@@ -0,0 +1,102 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/simulate_infections.R
+\name{forecast_infections}
+\alias{forecast_infections}
+\title{Forecast infections from a given fit and trajectory of the time-varying
+reproduction number}
+\usage{
+forecast_infections(
+  estimates,
+  R = NULL,
+  model = NULL,
+  samples = NULL,
+  batch_size = 10,
+  verbose = interactive()
+)
+}
+\arguments{
+\item{estimates}{The \code{estimates} element of an \code{\link[=epinow]{epinow()}} run that
+has been done with output = "fit", or the result of
+\code{\link[=estimate_infections]{estimate_infections()}} with \code{return_fit} set to TRUE.}
+
+\item{R}{A numeric vector of reproduction numbers; these will overwrite the
+reproduction numbers contained in \code{estimates}, except elements set to
+NA. Alternatively accepts a \verb{<data.frame>} containing at least \code{date} and
+\code{value} (integer) variables and optionally \code{sample}. More (or fewer) days
+than in the original fit can be simulated.}
+
+\item{model}{A compiled stan model as returned by \code{\link[rstan:stan_model]{rstan::stan_model()}}.}
+
+\item{samples}{Numeric, number of posterior samples to simulate from. The
+default is to use all samples in the \code{estimates} input.}
+
+\item{batch_size}{Numeric, defaults to 10. Size of batches in which to
+simulate. May decrease run times due to reduced IO costs but this is still
+being evaluated. If set to NULL then all simulations are done at once.}
+
+\item{verbose}{Logical defaults to \code{\link[=interactive]{interactive()}}. Should a progress bar
+(from \code{progressr}) be shown.}
+}
+\value{
+A list of output as returned by \code{\link[=estimate_infections]{estimate_infections()}} but based on
+results from the specified scenario rather than fitting.
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}}
+This function simulates infections using an existing fit to observed cases
+but with a modified time-varying reproduction number. This can be used to
+explore forecast models or past counterfactuals. Simulations can be run in
+parallel using \code{\link[future:plan]{future::plan()}}.
+}
+\examples{
+\donttest{
+# set number of cores to use
+old_opts <- options()
+options(mc.cores = ifelse(interactive(), 4, 1))
+
+# get example case counts
+reported_cases <- example_confirmed[1:50]
+
+# fit model to data to recover Rt estimates
+est <- estimate_infections(reported_cases,
+  generation_time = generation_time_opts(example_generation_time),
+  delays = delay_opts(example_incubation_period + example_reporting_delay),
+  rt = rt_opts(prior = list(mean = 2, sd = 0.1), rw = 7),
+  stan = stan_opts(control = list(adapt_delta = 0.9)),
+  obs = obs_opts(scale = list(mean = 0.1, sd = 0.01)),
+  gp = NULL, horizon = 0
+)
+
+# update Rt trajectory and simulate new infections using it
+R <- c(rep(NA_real_, 26), rep(0.5, 10), rep(0.8, 7))
+sims <- forecast_infections(est, R)
+plot(sims)
+
+# with a data.frame input of samples
+R_dt <- data.frame(
+  date = seq(
+    min(summary(est, type = "parameters", param = "R")$date),
+    by = "day", length.out = length(R)
+  ),
+  value = R
+)
+sims <- forecast_infections(est, R_dt)
+plot(sims)
+
+#' # with a data.frame input of samples
+R_samples <- summary(est, type = "samples", param = "R")
+R_samples <- R_samples[,
+ .(date, sample, value)][sample <= 1000][date <= "2020-04-10"
+]
+R_samples <- R_samples[date >= "2020-04-01", value := 1.1]
+sims <- forecast_infections(est, R_samples)
+plot(sims)
+
+options(old_opts)
+}
+}
+\seealso{
+\code{\link[=dist_spec]{dist_spec()}} \code{\link[=generation_time_opts]{generation_time_opts()}} \code{\link[=delay_opts]{delay_opts()}} \code{\link[=rt_opts]{rt_opts()}}
+\code{\link[=estimate_infections]{estimate_infections()}} \code{\link[=trunc_opts]{trunc_opts()}} \code{\link[=stan_opts]{stan_opts()}} \code{\link[=obs_opts]{obs_opts()}}
+\code{\link[=gp_opts]{gp_opts()}}
+}
diff --git a/man/simulate_infections.Rd b/man/simulate_infections.Rd
index 4b3046c6e..f1f3f577d 100644
--- a/man/simulate_infections.Rd
+++ b/man/simulate_infections.Rd
@@ -2,101 +2,19 @@
 % Please edit documentation in R/simulate_infections.R
 \name{simulate_infections}
 \alias{simulate_infections}
-\title{Simulate infections using a given trajectory of the time-varying
-reproduction number}
+\title{Deprecated; use \code{\link[=forecast_infections]{forecast_infections()}} instead}
 \usage{
-simulate_infections(
-  estimates,
-  R = NULL,
-  model = NULL,
-  samples = NULL,
-  batch_size = 10,
-  verbose = interactive()
-)
+simulate_infections(...)
 }
 \arguments{
-\item{estimates}{The \code{estimates} element of an \code{\link[=epinow]{epinow()}} run that
-has been done with output = "fit", or the result of
-\code{\link[=estimate_infections]{estimate_infections()}} with \code{return_fit} set to TRUE.}
-
-\item{R}{A numeric vector of reproduction numbers; these will overwrite the
-reproduction numbers contained in \code{estimates}, except elements set to
-NA. Alternatively accepts a \verb{<data.frame>} containing at least \code{date} and
-\code{value} (integer) variables and optionally \code{sample}. More (or fewer) days
-than in the original fit can be simulated.}
-
-\item{model}{A compiled stan model as returned by \code{\link[rstan:stan_model]{rstan::stan_model()}}.}
-
-\item{samples}{Numeric, number of posterior samples to simulate from. The
-default is to use all samples in the \code{estimates} input.}
-
-\item{batch_size}{Numeric, defaults to 10. Size of batches in which to
-simulate. May decrease run times due to reduced IO costs but this is still
-being evaluated. If set to NULL then all simulations are done at once.}
-
-\item{verbose}{Logical defaults to \code{\link[=interactive]{interactive()}}. Should a progress bar
-(from \code{progressr}) be shown.}
+\item{...}{Arguments to be passed to \code{\link[=forecast_infections]{forecast_infections()}}}
 }
 \value{
-A list of output as returned by \code{\link[=estimate_infections]{estimate_infections()}} but based on
-results from the specified scenario rather than fitting.
+the result of \code{\link[=forecast_infections]{forecast_infections()}}
 }
 \description{
-\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}}
-This function simulates infections using an existing fit to observed cases
-but with a modified time-varying reproduction number. This can be used to
-explore forecast models or past counterfactuals. Simulations can be run in
-parallel using \code{\link[future:plan]{future::plan()}}.
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}}
 }
-\examples{
-\donttest{
-# set number of cores to use
-old_opts <- options()
-options(mc.cores = ifelse(interactive(), 4, 1))
-
-# get example case counts
-reported_cases <- example_confirmed[1:50]
-
-# fit model to data to recover Rt estimates
-est <- estimate_infections(reported_cases,
-  generation_time = generation_time_opts(example_generation_time),
-  delays = delay_opts(example_incubation_period + example_reporting_delay),
-  rt = rt_opts(prior = list(mean = 2, sd = 0.1), rw = 7),
-  stan = stan_opts(control = list(adapt_delta = 0.9)),
-  obs = obs_opts(scale = list(mean = 0.1, sd = 0.01)),
-  gp = NULL, horizon = 0
-)
-
-# update Rt trajectory and simulate new infections using it
-R <- c(rep(NA_real_, 26), rep(0.5, 10), rep(0.8, 7))
-sims <- simulate_infections(est, R)
-plot(sims)
-
-# with a data.frame input of samples
-R_dt <- data.frame(
-  date = seq(
-    min(summary(est, type = "parameters", param = "R")$date),
-    by = "day", length.out = length(R)
-  ),
-  value = R
-)
-sims <- simulate_infections(est, R_dt)
-plot(sims)
-
-#' # with a data.frame input of samples
-R_samples <- summary(est, type = "samples", param = "R")
-R_samples <- R_samples[,
- .(date, sample, value)][sample <= 1000][date <= "2020-04-10"
-]
-R_samples <- R_samples[date >= "2020-04-01", value := 1.1]
-sims <- simulate_infections(est, R_samples)
-plot(sims)
-
-options(old_opts)
-}
-}
-\seealso{
-\code{\link[=dist_spec]{dist_spec()}} \code{\link[=generation_time_opts]{generation_time_opts()}} \code{\link[=delay_opts]{delay_opts()}} \code{\link[=rt_opts]{rt_opts()}}
-\code{\link[=estimate_infections]{estimate_infections()}} \code{\link[=trunc_opts]{trunc_opts()}} \code{\link[=stan_opts]{stan_opts()}} \code{\link[=obs_opts]{obs_opts()}}
-\code{\link[=gp_opts]{gp_opts()}}
+\details{
+Calling this function passes all arguments to \code{\link[=forecast_infections]{forecast_infections()}}
 }
diff --git a/tests/testthat/test-simulate_infections.R b/tests/testthat/test-simulate_infections.R
index dbe55e677..4d626fa8a 100644
--- a/tests/testthat/test-simulate_infections.R
+++ b/tests/testthat/test-simulate_infections.R
@@ -15,57 +15,63 @@ out <- suppressWarnings(estimate_infections(reported_cases,
 ))
 
 
-test_that("simulate_infections works to simulate a passed in estimate_infections object", {
-  sims <- simulate_infections(out)
+test_that("forecast_infections works to simulate a passed in estimate_infections object", {
+  sims <- forecast_infections(out)
   expect_equal(names(sims), c("samples", "summarised", "observations"))
 })
 
-test_that("simulate_infections works to simulate a passed in estimate_infections object with an adjusted Rt", {
+test_that("forecast_infections works to simulate a passed in estimate_infections object with an adjusted Rt", {
   R <- c(rep(NA_real_, 40), rep(0.5, 17))
-  sims <- simulate_infections(out, R)
+  sims <- forecast_infections(out, R)
   expect_equal(names(sims), c("samples", "summarised", "observations"))
   expect_equal(tail(sims$summarised[variable == "R"]$median, 9), rep(0.5, 9))
 })
 
-test_that("simulate_infections works to simulate a passed in estimate_infections object with a short adjusted Rt", {
+test_that("forecast_infections works to simulate a passed in estimate_infections object with a short adjusted Rt", {
   R <- c(rep(NA_real_, 40), rep(0.5, 17))
-  sims <- simulate_infections(out, R)
+  sims <- forecast_infections(out, R)
   expect_equal(names(sims), c("samples", "summarised", "observations"))
   expect_equal(tail(sims$summarised[variable == "R"]$median, 9), rep(0.5, 9))
 })
 
-test_that("simulate_infections works to simulate a passed in estimate_infections object with a long adjusted Rt", {
+test_that("forecast_infections works to simulate a passed in estimate_infections object with a long adjusted Rt", {
   R <- c(rep(NA_real_, 40), rep(1.2, 15), rep(0.8, 15))
-  sims <- simulate_infections(out, R)
-  sims10 <- simulate_infections(out, R, samples = 10)
+  sims <- forecast_infections(out, R)
+  sims10 <- forecast_infections(out, R, samples = 10)
   expect_equal(names(sims), c("samples", "summarised", "observations"))
   expect_equal(tail(sims$summarised[variable == "R"]$median, 30), R[41:70])
 })
 
 test_that("simulate infections can be run with a limited number of samples", {
   R <- c(rep(NA_real_, 40), rep(1.2, 15), rep(0.8, 15))
-  sims <- simulate_infections(out, R, samples = 10)
+  sims <- forecast_infections(out, R, samples = 10)
   expect_equal(names(sims), c("samples", "summarised", "observations"))
   expect_equal(tail(sims$summarised[variable == "R"]$median, 30), R[41:70])
   expect_equal(max(sims$samples$sample), 10)
 })
 
 test_that("simulate infections fails as expected", {
-  expect_error(simulate_infections())
-  expect_error(simulate_infections(out[-"fit"]))
+  expect_error(forecast_infections())
+  expect_error(forecast_infections(out[-"fit"]))
 })
 
-test_that("simulate_infections works to simulate a passed in estimate_infections object with an adjusted Rt in data frame", {
+test_that("forecast_infections works to simulate a passed in estimate_infections object with an adjusted Rt in data frame", {
   R <- c(rep(1.4, 40), rep(0.5, 17))
   R_dt <- data.frame(date = summary(out, type = "parameters", param = "R")$date, value = R)
-  sims_dt <- simulate_infections(out, R_dt)
+  sims_dt <- forecast_infections(out, R_dt)
   expect_equal(names(sims_dt), c("samples", "summarised", "observations"))
 })
 
-test_that("simulate_infections works to simulate a passed in estimate_infections object with samples of Rt in a data frame", {
+test_that("forecast_infections works to simulate a passed in estimate_infections object with samples of Rt in a data frame", {
   R_samples <- summary(out, type = "samples", param = "R")
   R_samples <- R_samples[, .(date, sample, value)][sample <= 1000]
   R_samples <- R_samples[date >= "2020-04-01", value := 1.1]
-  sims_sample <- simulate_infections(out, R_samples)
+  sims_sample <- forecast_infections(out, R_samples)
   expect_equal(names(sims_sample), c("samples", "summarised", "observations"))
 })
+
+test_that("simulate_infections is deprecated", {
+  expect_deprecated(
+    sims <- simulate_infections(out)
+  )
+})