Title: Simulation of Sediment Archived Climate Proxy Records
Version: 0.7.6
Description: Proxy forward modelling for sediment archived climate proxies such as Mg/Ca, d18O or Alkenones. The user provides a hypothesised "true" past climate, such as output from a climate model, and details of the sedimentation rate and sampling scheme of a sediment core. Sedproxy returns simulated proxy records. Implements the methods described in Dolman and Laepple (2018) <doi:10.5194/cp-14-1851-2018>.
License: MIT + file LICENSE
Encoding: UTF-8
Language: en-gb
LazyData: true
RoxygenNote: 7.3.2
Depends: R (≥ 3.5.0)
Imports: stats, mvtnorm, dplyr, tidyr, ggplot2, rlang
Suggests: testthat, knitr, rmarkdown
VignetteBuilder: knitr
URL: https://earthsystemdiagnostics.github.io/sedproxy/, https://github.com/EarthSystemDiagnostics/sedproxy
NeedsCompilation: no
Packaged: 2025-06-15 17:33:02 UTC; andrew-dolman
Author: Thomas Laepple [aut], Andrew Dolman [aut, cre]
Maintainer: Andrew Dolman <andrew.dolman@awi.de>
Repository: CRAN
Date/Publication: 2025-06-16 19:40:02 UTC

sedproxy: Simulation of Sediment Archived Climate Proxy Records

Description

The sedproxy package provides functions to simulate sediment archived proxies

Author(s)

Maintainer: Andrew Dolman andrew.dolman@awi.de

Authors:

See Also

Useful links:

Bioturbation weights

Description

For a given focal depth (or time), this function returns the probability that material collected from that depth was originally deposited at depth(s) z. In other words, that the material would have been found at depth z if there had been no bioturbation. It is the convolution of the depth solution from Berger and Heath (1968) with a uniform distribution to account for the width of the sediment layer from which samples were picked/extracted. It is a probability density function.

Usage

BioturbationWeights(
  z,
  focal.z,
  layer.width = 1,
  sed.acc.rate,
  bio.depth,
  scale = c("time", "depth")
)

Arguments

z

A vector of times or depths at which to evaluate the bioturbation weights

focal.z

The depth (or time) for which source dates are wanted

layer.width

the width of the sediment layer from which samples were taken, e.g. foraminifera were picked or alkenones were extracted, in cm. Defaults to 1 cm. If bio.depth and layer.width are both set to zero, each timepoint samples from a single year of the clim.signal, equivalent to sampling an annually laminated sediment core.

sed.acc.rate

Sediment accumulation rate in cm per 1000 years. Defaults to 50 cm per ka. Either a single value, or vector of same length as "timepoints"

bio.depth

Depth of the bioturbated layer in cm, defaults to 10 cm.

scale

Whether to scale depths by sediment accumulation rate to give positions in terms of time. Defaults to time.

Value

a numerical vector of weights.

References

Berger, W. H., & Heath, G. R. (1968). Vertical mixing in pelagic sediments. Journal of Marine Research, 26(2), 134–143.

Examples

z <- 0:10000
w <- BioturbationWeights(z, focal.z = 4000, layer.width = 1, sed.acc.rate = 5, bio.depth = 10)
plot(z, w, "l")

Simulate sediment archived proxy records from an input climate signal.

Description

ClimToProxyClim simulates the creation of a proxy climate record from a climate signal that is assumed to be true.

The following aspects of proxy creation are currently modelled.

  1. Seasonal bias in the encoding of a proxy due to the interaction between climate seasonality and any seasonality in the life cycle of the organism encoding the climate signal (e.g. Foraminifera for Mg/Ca ratios, or phytoplankton for Alkenone unsaturation indices).

  2. Bioturbation of the sediment archived proxy. For each requested timepoint, the simulated proxy consists of a weighted mean of the climate signal over a time window that is determined by the sediment accumulation rate sed.acc.rate and the bioturbation depth bio.depth which defaults to 10 cm. The weights are given by the depth solution to an impulse response function (Berger and Heath, 1968).

  3. Aliasing of seasonal and inter-annual climate variation onto to bioturbated (smoothed) signal. For proxies measured on a small number of discrete particles both seasonal and inter-annual climate variation is aliased into the proxy record. For example, Foraminifera have a life-cycle of approximately 1 month, so they record something like the mean temperature from a single month. If Mg/Ca is measured on e.g. n.samples = 30 individuals, the measured proxy signal is a mean of 30 distinct monthly mean temperatures and will thus be a stochastic sample of the true mean climate.

  4. Measurement noise/error is added as a pure Gaussian white noise process with mean = 0, standard deviation = sqrt(sigma.meas^2 + sigma.ind^2/n.samples).

  5. Additionally, a random bias can be added to each realisation of a proxy record. Bias is simulated as a Gaussian random variable with mean = 0, standard deviation = meas.bias. The same randomly generated bias value is applied to all timepoints in a simulated proxy record, when multiple replicate proxies are generated (n.replicates > 1) each replicate has a different bias applied.

ClimToProxyClim returns one or more replicates of the final simulated proxy as well as several intermediate stages (see section Value below).

Usage

ClimToProxyClim(
  clim.signal,
  timepoints,
  calibration.type = c("identity", "Uk37", "MgCa"),
  calibration = switch(calibration.type, identity = NA, Uk37 = "Mueller global", MgCa =
    "Ten planktonic species_350-500"),
  slp.int.means = NULL,
  slp.int.vcov = NULL,
  noise.type = switch(calibration.type, identity = "additive", Uk37 = "additive", MgCa =
    "multiplicative"),
  plot.sig.res = 100,
  habitat.weights = rep(1/ncol(clim.signal), ncol(clim.signal)),
  habitat.wt.args = NULL,
  bio.depth = 10,
  sed.acc.rate = 50,
  layer.width = 1,
  sigma.meas = 0,
  sigma.ind = 0,
  meas.bias = 0,
  scale.noise = switch(calibration.type, identity = FALSE, Uk37 = TRUE, MgCa = TRUE),
  n.samples = Inf,
  n.replicates = 1,
  top.of.core = NULL,
  n.bd = 3
)

Arguments

clim.signal

The "assumed true" climate signal, e.g. climate model output or instrumental record. A ts object consisting of a years x 12 (months) x n habitats (e.g. depths) matrix of temperatures. The time series should be at annual resolution and in reverse, i.e. "most recent timepoint first" order.

timepoints

The timepoints for which the proxy record is to be modelled

calibration.type

Type of proxy, e.g. Uk'37 or MgCa, to which the clim.signal is converted before the archiving and measurement of the proxy is simulated. Defaults to "identity" which means no conversion takes place.

calibration

The name of a specific calibration for which calibration parameters are provided by sedproxy. Currently applies only to calibration.type MgCa.

slp.int.means

Optional user supplied vector of values for the slope and intercept of the calibration function. Overides the defaults.

slp.int.vcov

Optional user supplied variance covariance matrix calibration parameters. Overides the defaults.

noise.type

Determines whether additive or multiplicative measurement noise is added. The appropriate type depends on the units of the proxy. Defaults to multiplicative for MgCa, additive for Uk'37 and identity (none) calibration types. Can be overidden with a string, "additive" or "multiplicative" in the case that pre-converted climate signal and measurement noise values are used in combination with an "identity" calibration type.

plot.sig.res

The resolution, in years, of the smoothed (block averaged) version of the input climate signal returned for plotting. This does not affect what the proxy model uses as input. If set to NA, no smoothed climate output is generated, this can speed up some simulations.

habitat.weights

Production weights for the proxy / proxy-carrier either as a vector of values with length = ncol(clim.signal), i.e. 1 weight for each month x habitat combination, a matrix of the same dimensions as the input climate signal matrix, or a function that produces an index of productivity as a function of temperature. Defaults to a vector of length = ncol(clim.signal) of equal weights.

habitat.wt.args

A named list of parameter values to be passed to a function named in habitat.weights.

bio.depth

Depth of the bioturbated layer in cm, defaults to 10 cm.

sed.acc.rate

Sediment accumulation rate in cm per 1000 years. Defaults to 50 cm per ka. Either a single value, or vector of same length as "timepoints"

layer.width

the width of the sediment layer from which samples were taken, e.g. foraminifera were picked or alkenones were extracted, in cm. Defaults to 1 cm. If bio.depth and layer.width are both set to zero, each timepoint samples from a single year of the clim.signal, equivalent to sampling an annually laminated sediment core.

sigma.meas

The standard deviation of the measurement error added to each simulated proxy value.

sigma.ind

The standard deviation of error between individuals (e.g. Forams) not otherwise modelled. This could included "vital effects" or aliasing of depth habitat variation not modelled via a depth resolved input climate signal and habitat weights. sigma.ind is scaled by n.samples before being combined with sigma.meas.

meas.bias

The amount of bias to add to each simulated proxy time-series. Each replicate proxy time-series has a constant bias added, drawn from a normal distribution with mean = 0, sd = meas.bias. Bias defaults to zero.

scale.noise

Scale noise to proxy units. Defaults to TRUE if calibration.type is not "identity"

n.samples

Number of e.g. Foraminifera sampled per timepoint, this can be either a single number, or a vector of length = timepoints. Can be set to Inf for non-discrete proxies, e.g. for Uk’37.

n.replicates

Number of replicate proxy time-series to simulate from the climate signal

top.of.core

The theoretical minimum age at the top of the core, ie. the year the core was sampled, defaults to the start of clim.in

n.bd

Number of multiples of the bioturbation width at which to truncate the bioturbation filter

Value

ClimToProxyClim returns an object of class "sedproxy.pfm", a list with three elements:

  1. a dataframe simulated.proxy

  2. a dataframe smoothed.signal

  3. a dataframe everything

The dataframe simulated.proxy contains a single realisation of the final forward modelled proxy, as well as the intermediate stages and the original climate signal at the requested timepoints.

The dataframe smoothed.signal contains a block averaged version the input climate signal, defaults to 100 year means but this is set by the parameter plot.sig.res. This is useful for plotting against the resulting simulated proxy.

The dataframe everything contains all of the above but with multiple replicates of the pseudo-proxy records if requested. The data are in "long form", with the column "stage" inidcating the proxy stage or input climate resolution and column "value" giving the values.

Named elements of the returned proxy record:

timepoints

Requested timepoints

clim.signal.ann

Input climate signal at requested timepoints at annual resolution

clim.signal.smoothed

Input climate signal at regular time intervals and resolution = plot.sig.res

clim.timepoints.ssr

Input climate signal at requested timepoints, smoothed to resolution = plot.sig.res

proxy.bt

Climate signal after bioturbation

proxy.bt.sb

Climate signal after bioturbation and habitat bias

proxy.bt.sb.inf.b

Climate signal after bioturbation, habitat bias, and calibration bias

proxy.bt.sb.inf.b.n

Climate signal after bioturbation, habitat bias, and measurement error

proxy.bt.sb.sampY

Climate signal after bioturbation, habitat bias, and aliasing of inter-annual variation

proxy.bt.sb.sampYM

Climate signal after bioturbation, habitat bias, and aliasing of inter-annual and intra-annual variation such as monthly temperatures or depth habitats

proxy.bt.sb.sampYM.b

Climate signal after bioturbation, habitat bias, and aliasing of inter-annual and intra-annual variation such as monthly temperatures or depth habitats, and calibration bias

proxy.bt.sb.sampYM.b.n

Climate signal after bioturbation, habitat bias, aliasing, and measurement error

simulated.proxy

Final simulated pseudo-proxy, this will be same as proxy.bt.sb.inf.b.n when n.samples = Inf, and proxy.bt.sb.sampYM.b.n when n.samples is finite

observed.proxy

True observed proxy (when supplied)

Examples

library(ggplot2)
set.seed(26052017)
clim.in <- ts(N41.t21k.climate[nrow(N41.t21k.climate):1,] - 273.15, start = -39)

PFM <- ClimToProxyClim(clim.signal = clim.in,
                       timepoints = round(N41.proxy$Published.age),
                       calibration.type = "identity",
                       habitat.weights = N41.G.ruber.seasonality,
                       sed.acc.rate = N41.proxy$Sed.acc.rate.cm.ka,
                       layer.width = 1,
                       sigma.meas = 0.46,
                       sigma.ind = 0, n.samples = Inf,
                       plot.sig.res = 10, meas.bias = 1,
                       n.replicates = 10)

PlotPFMs(PFM$everything, max.replicates = 1, stage.order = "seq") +
  facet_wrap(~stage)

PlotPFMs(PFM$everything, max.replicates = 1, stage.order = "var")

PlotPFMs(PFM$everything, stage.order = "var", plot.stages = "all")

Foraminifer Growth Rate Function from Lombard et al. (2009)

Description

Implements the function for foraminifer growth rate as a function of temperature from Lombard et al. (2009) with parametrization from FAME 1.0 (Roche et al, 2018).

Usage

ForamGrowthfT(
  foram = c("sacculifer", "bulloides", "pachy_d", "siphonifera", "universa", "pachy_s",
    "dutertrei", "ruber"),
  temperature_K,
  norm = FALSE,
  min.growth.thresh = 0
)

Arguments

foram

Name of foram.

temperature_K

Temperature in Kelvin

norm

Optional normalizing factor

min.growth.thresh

Sets a lower cutoff for growth as a proportion of the maximum growth rate for that taxon. For example in Roche et al (2018) a cutoff of 0.1 was used, meaning all growth rates less than 0.1*max were set to zero.

Value

A numerical vector or matrix with the same dimensions as the object passed to temperature_K. Units are daily growth rate, unless norm == TRUE.

References

Lombard, F., Labeyrie, L., Michel, E., Spero, H. J., and Lea, D. W.: Modelling the temperature dependent growth rates of planktic foraminifera, Marine Micropaleontology, 70, 1–7, https://doi.org/10.1016/j.marmicro.2008.09.004, 2009.

Roche, D. M., Waelbroeck, C., Metcalfe, B. and Caley, T.: FAME (v1.0): a simple module to simulate the effect of planktonic foraminifer species-specific habitat on their oxygen isotopic content, Geosci. Model Dev. Discuss., 2017, 1–22, doi:10.5194/gmd-2017-251, 2017.

Examples

ForamGrowthfT(foram = 'ruber', temperature_K = (c(280, 290)), norm = 1)

Convert "everything" part of output from ClimToProxyClim to dataframe. Used internally.

Description

Convert "everything" part of output from ClimToProxyClim to dataframe. Used internally.

Usage

MakePFMDataframe(PFM)

Arguments

PFM

output from ClimToProxyClim

Value

a dataframe

Seasonality of Globigerinoides ruber at core MD97-2141

Description

Seasonality of Globigerinoides ruber at core MD97-2141 predicted by the PLAFOM model (Fraile et al. 2008).

Usage

N41.G.ruber.seasonality

Format

A vector of 12 values

Source

Fraile, I., Schulz, M., Mulitza, S., & Kucera, M. (2008). Predicting the global distribution of planktonic foraminifera using a dynamic ecosystem model. Biogeosciences, 5(3), 891–911.

Mg/Ca proxy based temperature reconstruction for core MD97-2141

Description

Mg/Ca proxy based temperature reconstruction for core MD97-2141, originally published in Rosenthal et al. (2003), extracted from Shakun et al. (2012).

Usage

N41.proxy

Format

A data frame with 216 rows and 4 variables:

Published.age Age in yr BP
Published.temperature Reconstructed temperature in degrees C
Sed.acc.rate.cm.ka Sediment accumulation rate in cm per kyr

Details

Published age and published temperature were extracted from Shakun, J. D., Clark, P. U., He, F., Marcott, S. A., Mix, A. C., Liu, Z., … Bard, E. (2012). Global warming preceded by increasing carbon dioxide concentrations during the last deglaciation. Nature, 484(7392), 49–54. doi:10.1038/nature10915.

Sediment accumulation rates were estimated by fitting splines to the published age and depth estimates.

Source

Original reference: Rosenthal, Y., Oppo, D. W., & Linsley, B. K. (2003). The amplitude and phasing of climate change during the last deglaciation in the Sulu Sea, western equatorial Pacific. Geophysical Research Letters, 30(8), 1428. doi:10.1029/2002GL016612.

Metadata for datset N41.proxy

Description

Metadata for core MD97-2141 taken from Shakun et al. (2012)

Usage

N41.proxy.details

Format

A data frame with 1 rows and 17 variables:

Number Proxy ID number from Shakun et al. (2012)
ID.no Proxy ID number from Shakun et al. (2012) with prefix "N"
Core ID code of sediment core
Location Location of core
Proxy Proxy type
Lat Latutide of core
Lon Longitude of core
Elevation Depth of core-top below sea level in metres
Reference Original reference for the proxy record
Resolution Average time resolution in years
Calibration.ref Reference for Mg/Ca calibration
Calibration Mg/Ca calibration formula
Foram.sp Foram species analyses for proxy
Ref.14C Reference for carbon dating
Notes
Geo.cluster Coarse geographic location
Archive.type Type of proxy archive

Source

Shakun, J. D., Clark, P. U., He, F., Marcott, S. A., Mix, A. C., Liu, Z., … Bard, E. (2012). Global warming preceded by increasing carbon dioxide concentrations during the last deglaciation. Nature, 484(7392), 49–54. doi:10.1038/nature10915

Climate (surface temperature) at core MD97-2141 from TraCE-21ka

Description

Modelled surface temperature at core MD97-2141. Model output from TraCE-21ka simulations.

Usage

N41.t21k.climate

Format

A matrix with 22040 rows and 12 columns

Source

Liu, Z., Otto-Bliesner, B. L., He, F., Brady, E. C., Tomas, R., Clark, P. U., … Cheng, J. (2009). Transient Simulation of Last Deglaciation with a New Mechanism for Bølling-Allerød Warming. Science, 325(5938), 310–314. https://doi.org/10.1126/science.1171041

Plot forward modelled sedimentary proxies

Description

Plot forward modelled sedimentary proxies

Usage

PlotPFMs(
  PFMs,
  stage.order = c("var", "seq"),
  plot.stages = c("default"),
  max.replicates = 5,
  colr.palette = "default",
  alpha.palette = "default",
  levl.labels = "default"
)

Arguments

PFMs

A dataframe of forward modelled proxies

stage.order

Controls the order in which proxy stages are plotted, either sequentially, "seq", or in order of variance, "var". Defaults to var.

plot.stages

Proxy stages to be plotted, "default", "all", or a custom character vector

max.replicates

Maximum number of replicates to plot at once

colr.palette

Colours for the proxy stages

alpha.palette

Alpha levels for the proxy stages

levl.labels

Labels for the proxy stages

Value

a ggplot object of class "gg" "ggplot"

Examples

library(ggplot2)
set.seed(26052017)
clim.in <- ts(N41.t21k.climate[nrow(N41.t21k.climate):1,] - 273.15)

PFM <- ClimToProxyClim(clim.signal = clim.in,
                       timepoints = round(N41.proxy$Published.age),
                       calibration.type = "identity",
                       habitat.weights = N41.G.ruber.seasonality,
                       sed.acc.rate = N41.proxy$Sed.acc.rate.cm.ka,
                       sigma.meas = 0.45,
                       sigma.ind = 0,
                       n.samples = Inf,
                       plot.sig.res = 10, meas.bias = 1,
                       n.replicates = 10)

PlotPFMs(PFM$everything, max.replicates = 1, stage.order = "seq") +
  facet_wrap(~stage)

PlotPFMs(PFM$everything, max.replicates = 1, stage.order = "var")

PlotPFMs(PFM$everything, stage.order = "var", plot.stages = "all")

Convert between Temperature in Degrees C and Proxy Units

Description

A wrapper function for accessing proxy - temperature conversion functions

Usage

ProxyConversion(
  temperature = NULL,
  proxy.value = NULL,
  calibration.type = "identity",
  slp.int.means = NULL,
  slp.int.vcov = NULL,
  calibration = switch(calibration.type, identity = NA, Uk37 = "Mueller global", MgCa =
    "Ten planktonic species_350-500"),
  point.or.sample = c("point", "sample"),
  n = 1
)

Arguments

temperature

Temperature in degrees C

proxy.value

Temperature in proxy units

calibration.type

Type of proxy, e.g. Uk37 or MgCa

slp.int.means

Optional user supplied vector of values for the slope and intercept of the calibration function. Overides the defaults.

slp.int.vcov

Optional user supplied variance covariance matrix calibration parameters. Overides the defaults.

calibration

The name of a specific calibration for which calibration parameters are provided by sedproxy. Currently applies only to calibration.type MgCa.

point.or.sample

Use the "best estimate" calibration parameters, or parameters sampled from the fitted calibration model

n

the number of replicate conversions to make in the case of sampled calibration parameters

Details

Valid entries for calibration are: "Ten planktonic species_350-500", "G. aequilateralis_350-500", "G. aequilateralis_500-1000", "G. conglobatus_350-500", "G. hirsuta_350-500", "G. inflata_350-500", "G. ruber pink_250-350", "G. ruber pink_350-500", "G. ruber white_250-350", "G. ruber white_350-500", "G. sacculifer with sac_350-500", "G. sacculifer without sac_350-500", "G. truncatulinoides_350-500", "G. truncatulinoides_500-1000", "N. dutertrei_350-500", "O. univesa_350-500", "P. obliquiloculata_350-500"

Value

a vector of temperatures or proxy values

Examples

# From temperature to UK'37
## With fixed calibration
ProxyConversion(temperature = c(10, 20), point.or.sample = "point",
                calibration.type = "Uk37")

## With random calibration, 5 replicates
ProxyConversion(temperature = c(1, 2), n = 5, point.or.sample = "sample",
                calibration.type = "Uk37")


## Back-transformation with same calibration
ProxyConversion(
  proxy.value = as.vector(
    ProxyConversion(
      temperature = c(21, 22),
      calibration.type = "Uk37",
      point.or.sample = "point"
    )
  ),
  point.or.sample = "point",
  calibration.type = "Uk37"
)

## Back-transformation with random calibration
ProxyConversion(
  proxy.value = as.vector(
    ProxyConversion(
      temperature = c(21, 22),
      calibration.type = "Uk37",
     point.or.sample = "point"
    )
  )
  ,
  n = 5,
  point.or.sample = "sample",
  calibration.type = "Uk37"
)

## Incompatible arguments
## Not run: 
ProxyConversion(temperature = 1, proxy.value = 1)

## End(Not run)

Data for running the examples in package sedproxy

Description

A set of data objects for running the examples in sedproxy

N41.proxy

N41.proxy.details

N41.G.ruber.seasonality

N41.t21k.climate

Details

Example data

gisp2 ice core data at annual resolution

Description

gisp2 ice core data interpolated to annual resolution

Usage

gisp2.ann

Format

A data frame with 49885 rows and 3 variables:

age.yr.bp

Age in years BP

temperature

temperature in deg C

temperature.rescaled

temperature rescaled to resemble a d18O record as in Löwemark et al 2008

Details

GISP2 data from Alley, R.B.. 2004. GISP2 Ice Core Temperature and Accumulation Data. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series 2004-013. NOAA/NGDC Paleoclimatology Program, Boulder CO, USA. Interpolated to annual resolution and additionally rescaled to resemble a d18O record. As in Fig.2 of Löwemark, L., Konstantinou, K. I. and Steinke, S.: Bias in foraminiferal multispecies reconstructions of paleohydrographic conditions caused by foraminiferal abundance variations and bioturbational mixing: A model approach, Marine Geology, 256(1–4), 101–106, doi:10.1016/j.margeo.2008.10.005, 2008.

sedproxy parameters

Description

Parameters and variables required to generate a pseudo-proxy with ClimToProxyClim

Usage

param.tab

Format

A data frame with 12 rows and 4 variables:

Function argument

character Argument name in ClimToProxyClim

Description

character Description of argument and corresponding variable/parameter

Possible sources

character Source or possible sources of values for arguments

Default

character Default values for arguments

Scussolini et al. (2013) Table 1

Description

Data from table 1 in Scussolini et al. (2013)

Usage

scussolini.tab1

Format

A dataframe with 22 rows and 6 columns

Source

Scussolini, P., van Sebille, E., & Durgadoo, J. V. (2013). Paleo Agulhas rings enter the subtropical gyre during the penultimate deglaciation. Climate of the Past, 9(6), 2631–2639. https://doi.org/10.5194/cp-9-2631-2013

Labels for proxy stages

Description

Labels for proxy stages. For plotting.

Usage

stage.labels

Format

A named character vector

Description of proxy stages

Description

A description of the proxy stages in the output of ClimToProxyClim and default labels, colours and order for plotting

Usage

stages.key

Format

A data frame with 13 rows and 6 variables:

stage proxy stage
label label for proxy stage
description description of proxy stage
plot.order default plotting order of stages
plotting.colour default colour for plotting
plotting.alpha default alpha level for plotting