R/simulate.R

In paleolimbot/pb210: Lead-210 dating utilities

#' Simulate sediment accumulation
#'
#' `pb210_simulate_accumulation()` simulates the accumulation of lead-210 over
#' time, returning a data frame with one row per `time_step`, whereas
#' `pb210_simulate_core()` returns a data frame with one row per `depth_step`
#' (like extruding a core). These functions were created to test the
#' [pb210_crs()] and [pb210_cic()] functions, and to create
#' the lead-210 theory vignette. Ages are in years, depths are in cm,
#' masses are in kg, densities are in kg / m^3^, and specific activities
#' are in Bq / kg. [pb210_simulate_counting()] simulates putting a subsample
#' of each slice on a counter for a specified amount of time. This allows
#' assigning an error to the known lead-210 activity, as well as provide
#' more realistic data for simulations. See [pb210_error_from_activity()]
#' for details.
#'
#' @param max_age The maximum age of the simulation (years)
#' @param time_step The time to consider in each step (years)
#' @param supply A function such as that generated by [pb210_supply_constant()].
#'   The flux of lead-210 to the surface, in Bq/ m^2^ / year
#' @param mass_accumulation A function such as that generated by [pb210_mass_accumulation_constant()].
#'   Values returned by the function should be in kg / m^2^ / year.
#' @param compressibility A function such as that generated by [pb210_compressibility_constant()].
#'   Values returned by the function should be in 1 / Pa.
#' @param initial_density A function such as that generated by [pb210_density_constant()].
#'   Values returned by the function should be in kg / m^3^.
#' @param decay_constant The decay constant to use (see [pb210_decay_constant()])
#' @param accumulation An age simulation, created by [pb210_simulate_accumulation()].
#' @param depth_step A vector of depth steps to consider, from the top to the bottom
#'   of the core. There should be one element in the vector for each slice considered.
#' @param core_area The internal area of the core in m^2^.
#' @inheritParams pb210_error_from_activity
#'
#' @return A tibble with columns:
#'
#'   - **age** (years): the mean age of the slice (weighted by mass)
#'   - **depth** (cm): the midpoint depth of the slice
#'   - **activity** (Bq/kg)
#'   - **age_top**, **age_bottom** (years)
#'   - **depth_top**, **depth_bottom** (cm)
#'   - **slice_mass** (kg)
#'   - **slice_density** (kg / m^3^)
#'
#'   The output of `pb210_simulate_counting()` also includes the columns:
#'
#'   - **activity_estimate** (Bq/kg): A plausible activity that might be
#'     obtained by a lead-210 counting method, sampled using [stats::rpois()].
#'   - **activity_se** (Bq/kg): The standard error (estimate of the standard
#'     deviation) of **activity_estimate**. For the expected error,
#'     use [pb210_error_from_activity()] on
#'
#'
#' @importFrom rlang .data
#' @export
#'
#' @examples
#' # 1 row per year
#' accumulation <- pb210_simulate_accumulation()
#' # 1 row per 0.5 cm
#' core <- pb210_simulate_core(accumulation)
#' counted_core <- pb210_simulate_counting(core)
#'
pb210_simulate_accumulation <- function(mass_accumulation = pb210_mass_accumulation_constant(),
                                        max_age = 500, time_step = 1,
                                        supply = pb210_supply_constant(),
                                        compressibility = pb210_compressibility_constant(),
                                        initial_density = pb210_density_constant(),
                                        decay_constant = pb210_decay_constant()) {
  # we do not propogate error in this simulation
  decay_constant <- without_errors(decay_constant)

  stopifnot(
    is.numeric(max_age), length(max_age) == 1,
    is.numeric(time_step), length(time_step) == 1, time_step > 0,
    is.numeric(decay_constant), length(decay_constant) == 1, decay_constant > 0
  )

  # all calculations here are for a 1 m^2 section.
  core_area <- 1 # m^2

  mass_accumulation <- rlang::as_function(mass_accumulation, env = rlang::caller_env())
  supply <- rlang::as_function(supply, env = rlang::caller_env())
  compressibility <- rlang::as_function(compressibility, env = rlang::caller_env())
  initial_density <- rlang::as_function(initial_density, env = rlang::caller_env())

  tbl <- tibble::tibble(
    # ages
    age_top = seq(max_age - time_step, 0, by = -time_step),
    age_bottom = seq(max_age, 0 + time_step, by = -time_step),
    age = (.data$age_top + .data$age_bottom) / 2,

    # initial values
    initial_pb210_content = supply(.data$age) * time_step, # Bq
    slice_mass = mass_accumulation(.data$age) * time_step * core_area, # kg
    compressibility = compressibility(.data$age), # 1/Pa
    initial_density = initial_density(.data$age), # kg / m^3

    # radioactive decay of pb210, normalized to mass
    activity = .data$initial_pb210_content * exp(-decay_constant * .data$age) /
      .data$slice_mass, # Bq / kg

    # compression of sediment
    cumulative_mass_above = c(rev(cumsum(rev(.data$slice_mass[-1]))), 0), # kg
    delta_pressure = 9.806 * .data$cumulative_mass_above, # Pa
    relative_volume = 1 - (.data$delta_pressure * .data$compressibility), # unitless (volume / volume at deposition)

    # use density and relative_volume to map slice mass to thickness
    # kg / (kg / m^3) * unitless = m^3
    slice_volume = .data$slice_mass / .data$initial_density * .data$relative_volume,
    slice_density = .data$slice_mass / .data$slice_volume,
    thickness = .data$slice_volume / core_area * 100, # cm
    depth_bottom = rev(cumsum(rev(.data$thickness))),
    depth_top = c(rev(cumsum(rev(.data$thickness[-1]))), 0),
    depth = (.data$depth_top + .data$depth_bottom) / 2
  )

  # sanity checks on outputs
  stopifnot(
    all(tbl$slice_mass >= 0),
    all(tbl$initial_pb210_content >= 0),
    all(tbl$initial_density > 0)
  )

  # return like most core data, backwards in time
  tbl[
    order(tbl$depth),
    c(
      "age", "depth", "activity", "age_top", "age_bottom",
      "depth_top", "depth_bottom", "slice_mass", "slice_density"
    )
  ]
}


#' @rdname pb210_simulate_accumulation
#' @export
pb210_simulate_core <- function(accumulation = pb210_simulate_accumulation(),
                                depth_step = rep(0.5, 60), core_area = pb210_core_area()) {
  stopifnot(
    is.data.frame(accumulation),
    all(
      c(
        "age", "depth_top", "depth_bottom", "age_top", "age_bottom",
        "activity", "slice_mass",
        "slice_density"
      ) %in% colnames(accumulation)
    ),
    is.numeric(depth_step), all(depth_step > 0),
    is.numeric(core_area), length(core_area) == 1, core_area > 0
  )

  # all calculations for `accumulation` were for 1 m^2^. slice_mass is the only parameter
  # that depends on this. For the calculations in this resampling to be correct,
  # we need to rescale slice_mass here.
  accumulation$slice_mass <- accumulation$slice_mass * core_area / 1.

  # calculations on accumulation that apply to all of the new depth slices
  accumulation$thickness <- accumulation$depth_bottom - accumulation$depth_top # cm
  accumulation$slice_volume <- accumulation$slice_mass / accumulation$slice_density # m^3

  accumulation$activity_total <- accumulation$activity * accumulation$slice_mass # Bq
  accumulation$age_total <- accumulation$age_bottom - accumulation$age_top # yr
  accumulation$sediment_accumulation_rate <- accumulation$age_total / accumulation$thickness # yr

  depths <- tibble::tibble(
    # make new depth slices
    depth_top = c(0, cumsum(depth_step[-1])), # cm
    depth_bottom = cumsum(depth_step), # cm
    depth = (.data$depth_top + .data$depth_bottom) / 2, # cm

    # these attributes are calculated for each new depth slice
    interpolated_attributes = mapply(
      function(top, bottom) {
        # these calculations only apply to this new depth slice
        # they are all of length nrow(accumulation)
        depth_bottom_match <- pmin(bottom, accumulation$depth_bottom) # cm depth
        depth_top_match <- pmax(top, accumulation$depth_top) # cm depth
        depth_in_slice <- pmax(0, depth_bottom_match - depth_top_match) # cm
        proportion_in_slice <- depth_in_slice / accumulation$thickness # proportion
        mass_in_slice <- accumulation$slice_mass * proportion_in_slice # kg

        # return a 1-row tibble of summary values
        tibble::tibble(
          slice_mass = sum(mass_in_slice),
          activity = sum(accumulation$activity_total * proportion_in_slice) /
            .data$slice_mass,
          age = stats::weighted.mean(accumulation$age, mass_in_slice),
          slice_volume = (bottom - top) / 100 * core_area,
          slice_density = .data$slice_mass / .data$slice_volume
        )
      },
      .data$depth_top, .data$depth_bottom,
      SIMPLIFY = FALSE
    )
  )

  # interpolate age_top and age_bottom using stats::approx
  known_depths <- c(0, accumulation$depth_bottom)
  known_ages <- c(0, accumulation$age_bottom)
  depths$age_top <- stats::approx(known_depths, known_ages, depths$depth_top)$y
  depths$age_bottom <- stats::approx(known_depths, known_ages, depths$depth_bottom)$y

  tbl <- tibble::as_tibble(
    cbind(
      depths[c("depth_top", "depth_bottom", "depth", "age_top", "age_bottom")],
      do.call(rbind, depths$interpolated_attributes)
    )
  )

  # the same format as the pb210_simulate_accumulation()
  tbl[
    order(tbl$depth),
    c(
      "age", "depth", "activity", "age_top", "age_bottom",
      "depth_top", "depth_bottom", "slice_mass", "slice_density"
    )
  ]
}

#' @rdname pb210_simulate_accumulation
#' @export
pb210_simulate_counting <- function(accumulation = pb210_simulate_core(),
                                    count_mass = 0.5 / 1000,
                                    count_time = lubridate::ddays(1)) {
  stopifnot("activity" %in% colnames(accumulation))
  check_count_params(accumulation$activity, count_mass, count_time)
  count_time <- count_time_seconds(count_time)

  # model counts as a poisson distribution with lambda = the expected number of counts
  expected_decays <- pb210_counts_from_activity(
    accumulation$activity,
    count_mass,
    count_time
  )
  sampled_decays <- stats::rpois(nrow(accumulation), lambda = expected_decays) # counts

  # add estimate and standard error to the
  accumulation$activity_estimate <- pb210_activity_from_counts(
    sampled_decays,
    count_mass,
    count_time
  )
  accumulation$activity_se <- pb210_error_from_counts(
    sampled_decays,
    count_mass,
    count_time
  )

  accumulation
}


#' Parameter generators for lead-210 supply
#'
#' @param value An flux of lead 210, in Bq / m^2^ / year.
#'
#' @return A function of a single parameter, `age`, which is a (decreasing) vector of ages.
#' @export
#'
#' @references
#' Appleby, P.G. 2008. Three decades of dating recent sediments by fallout
#' radionuclides: a review. The Holocene, 18: 83–93. <https://10.1177/0959683607085598>
#'
#' @examples
#' ages <- 150:0
#'
pb210_supply_constant <- function(value = 100)  {
  force(value)
  function(age) {
    value
  }
}


#' Parameter generators for mass accumulation rates
#'
#' @param value,mean An accumulation rate, in kg / m2 / year.
#' @param sd For random (log)normal accumulation rates, the log standard deviation.
#'   The default is 1, which results in highly variable sedimentation rates. This is
#'   useful for testing [pb210_crs()].
#'
#' @return A function of a single parameter, `age`, which is a (decreasing) vector of ages.
#' @export
#'
#' @examples
#' age_compare <- tibble::tibble(
#'   ages = 150:0,
#'   constant = pb210_mass_accumulation_constant()(ages),
#'   rlnorm = pb210_mass_accumulation_rlnorm()(ages)
#' )
#'
pb210_mass_accumulation_constant <- function(value = 0.150)  {
  force(value)
  function(age) {
    value
  }
}

#' @rdname pb210_mass_accumulation_constant
#' @export
pb210_mass_accumulation_rlnorm <- function(mean = 0.150, sd = 1) {
  force(mean)
  force(sd)
  function(age) {
    stats::rlnorm(length(age), meanlog = log(mean), sdlog = sd)
  }
}

#' Parameter generators for density values
#'
#' Note that this refers to dry density at the time of deposition.
#'
#' @param value A density at the time of deposition, in kg / m^3^
#'
#' @return A function of a single parameter, `age`, which is a (decreasing) vector of ages.
#' @export
#'
#' @examples
#' ages <- 150:0
#' pb210_density_constant()(ages)
#'
pb210_density_constant <- function(value = 150)  {
  force(value)
  function(age) {
    value
  }
}

#' Parameter generators for compressibility
#'
#' @param value A compressibility at the time of deposition, in 1/Pa.
#'   Higher values (closer to 0) indicate a higher propensity for the material
#'   to be compressed.
#'
#' @return A function of a single parameter, `age`, which is a (decreasing) vector of ages.
#' @export
#'
#' @examples
#' ages <- 150:0
#' pb210_compressibility_constant()(ages)
#'
pb210_compressibility_constant <- function(value = 1e-3)  {
  force(value)
  function(age) {
    value
  }
}

#' @rdname pb210_compressibility_constant
#' @export
pb210_compressibility_none <- function()  {
  pb210_compressibility_constant(0)
}