In EarthSystemDiagnostics/proxysnr: Separate Signal and Noise in Climate Proxy Records

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Overview

This vignette documents the calculation of the spectral transfer functions for the effects of vapour diffusion and time uncertainty. The transfer functions describe the loss in spectral power as a function of frequency that is expected due to these two processes: for time uncertainty, the loss applies to the spectrum of the common signal recorded by the spatial average of an array of proxy records (the "stacked" record); for diffusion, it applies to the loss in average spectral power of the individual records. Proxy records refer here, more specifically, to any proxy time series that are based on layer counting in the case of time uncertainty and to stable isotope records from firn/ice cores in the case of diffusion.

The documentation utilises as examples the code for calculating the transfer functions for the DML and WAIS isotope records studied in Münch and Laepple (2018). In there, you will also find a detailed description of the basic method (see Appendix B).

Dependencies

library(proxysnr)

For the modelling of time uncertainty we need the package simproxyage:

# check if package is available
pkg.flag <- requireNamespace("simproxyage", quietly = TRUE)

if (!pkg.flag) {
   message(paste("Package \"simproxyage\" not found.",
                 "No time uncertainty transfer functions will be calculated."))
} else {
  library(simproxyage)
}

If you have installed proxysnr without additionally installing the package simproxyage and you try to rebuild this vignette, you will receive the warning message here that simproxyage is not available; then, no time uncertainty transfer functions can be calculated. To do so, head over to the github page for simproxyage and install the package using remotes::install_github.

General settings

First, we set general parameters that are required for the calculations.

We set a seed for reproducible results:

set.seed(4122018)

For convenience, list objects are used for saving the transfer functions:

# for diffusion transfer function
dtf <- list()
# for time uncertainty transfer function
ttf <- list()

A specific resolution of the created surrogate white noise time series is used for the diffusion calculations, for which, in general, a higher resolution compared to the target resolution of the isotope records is adopted. This is necessary for numerical stability of the diffusion results at the high-frequency end.

Here, the target resolution of the isotope records is 1 yr (annual resolution) and a simulation resolution of 0.5 yr has proven sufficient for the numerical stability of the results:

res <- 1/2

If the simulation resolution differs from the target resolution, also the simulated frequency axis of the diffusion transfer function will differ from the respective frequency axis of the proxy spectra. Therefore, the transfer functions are interpolated in frequency space on the target frequency axis of the proxy spectra, which hence needs to be supplied.

Here, for convenience, we load the required target frequency axes from the data supplied with the package.

ff <- list()
ff$dml1 <- diffusion.tf$dml1$freq
ff$dml2 <- diffusion.tf$dml2$freq
ff$wais <- diffusion.tf$wais$freq

Finally, we specify the number of Monte Carlo simulations used to obtain the transfer functions. In general, an average of a large number of simulations is needed to reduce the spectral uncertainties and obtain smooth estimates. However, since this takes considerable amount of computation time, for the vignette building only a small number of simulations is used here:

ns <- 10

Note that for the results of Münch and Laepple (2018), ns has been set to 100,000, which required a computation time of 1--3 hours depending on the data set (using a 2.4 GHz Intel Core i5 processor).

DML1 isotope records

In the following, we document all data-set-specific parameters needed for the calculations along the example of the DML1 data set of isotope records.

Specify the number of proxy records ("cores") in the studied array:

nc <- 15

Specify the time axis common to all proxy records in the array in units of the original (target) resolution of the records (here, years), which gives the number of proxy data points:

# set time axis
t <- seq(1994, 1801)
# number of proxy data points
nt <- length(t)

Specify the time points of age control points where the time uncertainty is assumed to be zero:

acp <- c(1994, 1884, 1816, 1810)

Set the rate of the random process which perturbs the proxy record chronologies, so here the rate of missing or double-counting of a layer. We use symmetric rates here:

rate <- 0.013

Specify a numeric vector with two entries to use asymmetric rates.

For the diffusion simulations, estimates of the diffusion length (in units of the record resolution) are needed as a function of time since deposition (see the parameter sigma in the call of the function DiffusionTF below, and ?DiffusionTF). Along with proxysnr, the diffusion length estimates for the DML and WAIS cores are provided in diffusion.length; for details see ?diffusion.length. Calculation of these estimates is documented in the package source under ./data-raw/calculate-diffusion-length.R, where . is the root of the package source directory. The relevant site parameters needed for the diffusion length calculations are provided in clim.site.par; for details see ?clim.site.par. For further information see the publication Münch and Laepple (2018).

Finally, run the calculations and store only the ratio of average output over input spectra, i.e. the transfer function.

Run the diffusion model:

tmp <- DiffusionTF(nt = nt / res, nc = nc, ns = ns,
                   sigma = diffusion.length$dml1[, -1],
                   res = res)$ratio

Since we used a certain simulation resolution, we have to interpolate onto the target resolution of our proxy data. Additionally, we store meta data on the simulations:

dtf$dml1 <- tmp
dtf$dml1$freq <- ff$dml1
dtf$dml1$spec <- approx(tmp$freq, tmp$spec, dtf$dml1$freq)$y

attr(dtf$dml1, "version") <- sprintf("Creation date: %s.", Sys.time())
attr(dtf$dml1, "N.sim") <- sprintf("Number of simulations used: N = %s",
                       formatC(ns, big.mark = ",", format = "d"))
attr(dtf$dml1, "log-smooth") <- "Log-smooth applied: No."

class(dtf$dml1) <- "spec"

Run the time uncertainty model with Gaussian white noise as the surrogate data and document the meta data:

if (pkg.flag) {
   ttf$dml1 <- TimeUncertaintyTF(t = t, acp = acp, nt = nt, nc = nc, ns = ns,
                                 rate = rate, surrogate.fun = rnorm)$ratio
} else {
  ttf$dml1$freq <- dtf$dml1$freq
  ttf$dml1$spec <- rep(NA, length(dtf$dml1$freq))
}

attr(ttf$dml1, "version") <- sprintf("Creation date: %s.", Sys.time())
attr(ttf$dml1, "rate") <- sprintf("Process rate used: %1.3f.", rate)
attr(ttf$dml1, "N.sim") <- sprintf("Number of simulations used: N = %s",
                       formatC(ns, big.mark = ",", format = "d"))
attr(ttf$dml1, "log-smooth") <- "Log-smooth applied: No."

class(ttf$dml1) <- "spec"

The other isotope records

We run the calculations for the DML2 isotope records...

nc <- 3
t <- seq(1994, 1000)
nt <- length(t)
acp <- c(1994, 1884, 1816, 1810, 1459, 1259)
rate <- 0.013

tmp <- DiffusionTF(nt = nt / res, nc = nc, ns = ns,
                   sigma = diffusion.length$dml2[, -1],
                   res = res)$ratio

dtf$dml2 <- tmp
dtf$dml2$freq <- ff$dml2
dtf$dml2$spec <- approx(tmp$freq, tmp$spec, dtf$dml2$freq)$y

if (pkg.flag) {
   ttf$dml2 <- TimeUncertaintyTF(t = t, acp = acp, nt = nt, nc = nc, ns = ns,
                                 rate = rate, surrogate.fun = rnorm)$ratio
} else {
  ttf$dml2$freq <- dtf$dml2$freq
  ttf$dml2$spec <- rep(NA, length(dtf$dml2$freq))
}

attr(dtf$dml2, "version") <- sprintf("Creation date: %s.", Sys.time())
attr(dtf$dml2, "N.sim") <- sprintf("Number of simulations used: N = %s",
                       formatC(ns, big.mark = ",", format = "d"))
attr(dtf$dml2, "log-smooth") <- "Log-smooth applied: No."

attr(ttf$dml2, "version") <- sprintf("Creation date: %s.", Sys.time())
attr(ttf$dml2, "rate") <- sprintf("Process rate used: %1.3f.", rate)
attr(ttf$dml2, "N.sim") <- sprintf("Number of simulations used: N = %s",
                       formatC(ns, big.mark = ",", format = "d"))
attr(ttf$dml2, "log-smooth") <- "Log-smooth applied: No."

class(dtf$dml2) <- "spec"
class(ttf$dml2) <- "spec"

... and the WAIS isotope records:

nc <- 5
t <- seq(2000, 1800)
nt <- length(t)
acp <- c(2000, 1992, 1964, 1885, 1837, 1816, 1810)
rate <- 0.027

tmp <- DiffusionTF(nt = nt / res, nc = nc, ns = ns,
                   sigma = diffusion.length$wais[, -1],
                   res = res)$ratio

dtf$wais <- tmp
dtf$wais$freq <- ff$wais
dtf$wais$spec <- approx(tmp$freq, tmp$spec, dtf$wais$freq)$y

if (pkg.flag) {
   ttf$wais <- TimeUncertaintyTF(t = t, acp = acp, nt = nt, nc = nc, ns = ns,
                                 rate = rate, surrogate.fun = rnorm)$ratio
} else {
  ttf$wais$freq <- dtf$wais$freq
  ttf$wais$spec <- rep(NA, length(dtf$wais$freq))
}

attr(dtf$wais, "version") <- sprintf("Creation date: %s.", Sys.time())
attr(dtf$wais, "N.sim") <- sprintf("Number of simulations used: N = %s",
                       formatC(ns, big.mark = ",", format = "d"))
attr(dtf$wais, "log-smooth") <- "Log-smooth applied: No."

attr(ttf$wais, "version") <- sprintf("Creation date: %s.", Sys.time())
attr(ttf$wais, "rate") <- sprintf("Process rate used: %1.3f.", rate)
attr(ttf$wais, "N.sim") <- sprintf("Number of simulations used: N = %s",
                       formatC(ns, big.mark = ",", format = "d"))
attr(ttf$wais, "log-smooth") <- "Log-smooth applied: No."

class(dtf$wais) <- "spec"
class(ttf$wais) <- "spec"

Plot the transfer functions

The obtained transfer functions can be plotted with:

PlotTF(dtf = dtf, ttf = ttf,
       names = c("DML1", "DML2", "WAIS"),
       col = c("black", "firebrick", "dodgerblue"))

The transfer functions for simulations with ns = 10^5 as used and shown in Münch and Laepple (2018) are provided along proxysnr in the variables diffusion.tf and time.uncertainty.tf; see ?diffusion.tf, ?time.uncertainty.tf, str(diffusion.tf), str(time.uncertainty.tf) for details.

These can be plotted by simply calling

PlotTF(names = c("DML1", "DML2", "WAIS"),
       col = c("black", "firebrick", "dodgerblue"))

Literature cited

Münch, T. and Laepple, T.: What climate signal is contained in decadal- to centennial-scale isotope variations from Antarctic ice cores?, Clim. Past, 14, 2053-2070, doi: 10.5194/cp-14-2053-2018, 2018.

EarthSystemDiagnostics/proxysnr documentation built on Oct. 2, 2021, 3:03 p.m.