Diffusion transfer function

This function implements an empirical Monte Carlo approach to estimate the spectral transfer function for the effect of firn diffusion on the spatial average of firn/ice-core stable isotope records.

Usage

CalculateDiffusionTF(
  nt,
  nc,
  ns,
  sigma,
  res = 1,
  window = NULL,
  coherent = FALSE,
  df.log = NULL,
  verbose.output = FALSE,
  ...
)

Arguments

nt: the length of the modelled isotope records (i.e. the number of data points in each record).
nc: the number of cores in the modelled core array.
ns: the number of Monte Carlo simulations for estimating the transfer function.
sigma: numeric vector of length nt or numeric array of dimension nt * nc providing diffusion length values. The nt diffusion length values are assumed to correspond to the respective nt isotope values. If only a numeric vector is provided, it is assumed that these diffusion lengths are valid for all nc cores. If an array is provided, each column provides the diffusion lengths for the respective core. Note that the units of the diffusion length must match the units of res.
res: the sampling (e.g., temporal) resolution of the isotope data; determines the frequency axis of the transfer function.
window: length-2 vector giving a start and an end time (within 1:nt) offering the possibility to only use a subset of the total length of the simulated records for the transfer function analysis, while the default of NULL means to use the records' entire lengths.
coherent: if TRUE, nc identical white noise time series are assumed to estimate the transfer function; else (the default) nc independent noise series.
df.log: width of the Gaussian kernel in logarithmic frequency units to smooth the spectral estimates; NULL (the default) suppresses smoothing. In general, smoothing should not be necessary when a sufficiently large number of Monte Carlo simulations (parameter ns) is used; nevertheless, it can be switched on here when the transfer function still appears too noisy.
verbose.output: logical controlling the size of the return object; per default, only the transfer function spectrum is returned, else also the spectra whose ratio determines the transfer function (see Details).
...: additional parameters which are passed to the spectral estimation function SpecMTM.

Value

either a spectral object (?spec.object) of the transfer function if verbose.output = FALSE (default), or a list of the spectral objects signal, diffused and ratio, providing the averages over the ns simulations of:

signal:: the undiffused spectrum;
diffused:: the diffused spectrum;
ratio: their ratio (diffused/undiffused), i.e. the transfer function.

Details

The approach is described in detail in Münch and Laepple (2018). In brief, nc Gaussian white noise time series are created and diffused and the average of these time series is calculated. The process is repeated ns times. For each of the ns realisations, spectra of the average diffused and undiffused records are calculated; subsequently, the ns spectra are averaged, and the ratio of the average diffused to the average undiffused spectrum yields the spectral transfer function.

Diffusion is modelled as the convolution of the undiffused record with a Gaussian with standard deviation given by the diffusion length sigma. The spectral estimates are calculated using Thomson’s multitaper method with three windows with linear detrending before analysis.

References

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, https://doi.org/10.5194/cp-14-2053-2018, 2018.

Author