A spectral approach to estimating the timescale-dependent uncertainty of paleoclimate records – Part 1: Theoretical concept

Proxy records represent an invaluable source of information for reconstructing past climatic variations, but they are associated with considerable uncertainties. For a systematic quantification of these reconstruction errors, however, knowledge is required not only of their individual sources but also of their auto-correlation structure as this determines the timescale dependence of their magnitude, an issue that has been often ignored until now. Here a spectral approach to uncertainty analysis is provided for paleoclimate reconstructions obtained from single sediment proxy records. The formulation in the spectral domain rather than the time domain allows for an explicit demonstration and quantification of the timescale dependence that is inherent in any proxy-based reconstruction uncertainty. This study is published in two parts.

In this first part, the theoretical concept is presented, and analytic expressions are derived for the power spectral density of the reconstruction error of sediment proxy records. The underlying model takes into account the spectral structure of the climate signal, seasonal and orbital variations, bioturbation, sampling of a finite number of signal carriers, and uncorrelated measurement noise, and it includes the effects of spectral aliasing and leakage. The uncertainty estimation method, based upon this model, is illustrated by simple examples. In the second part of this study, published separately, the method is implemented in an application-oriented context, and more detailed examples are presented.

1 Introduction

The central issues of climate sciences include the estimation, understanding, and prediction of climatic variations across ranges of space and timescales that are relevant to the specific field of study. From an inductive perspective, such studies are necessarily based on observational data which the variability may be estimated from, whereas from a deductive perspective observational data are needed in the course of the validation of theories and models. For certain fields of study, instrumental or satellite data may provide a useful data source. Nonetheless, once processes are studied that involve climate states or variations at times before the instrumental era or that involve timescales longer than this, reconstructions obtained from paleoclimate proxies become indispensable. Such proxy records reveal imprints of past climatic conditions created by, for example, impacts on the calcification of the shells of marine organisms (Nürnberg et al., 1996), now preserved in sea sediments, on terrestrial pollen assemblages archived in lake sediments (Birks and Seppä, 2004), or on stable water isotopes that can be recovered from ice cores (Jouzel et al., 1997). Proxy-based reconstructions, however, are associated with notable uncertainties that are often much larger than those of instrumental data (Münch and Laepple, 2018; Reschke et al., 2019) and that can emerge from a variety of sources – they are essentially highly noisy and distorted observations of selected climate variables. Hence, an important task of the paleoclimate research field is to provide thorough quantitative estimates of these reconstruction uncertainties.

Reconstruction uncertainties may arise, for example, during the sampling and measurement procedure associated with measurement errors occurring in the laboratory (Rosell-Melé et al., 2001; Greaves et al., 2008) and with aliasing of variability from frequencies higher than those resolved (e.g., from El Niño–Southern Oscillation, ENSO, or the seasonal cycle; see, for example, Thirumalai et al., 2013; Laepple et al., 2018). Uncertainties may also arise from archive formation processes associated with errors induced by smoothing processes like bioturbation affecting sediment archives (Berger and Heath, 1968; Goreau, 1980) or diffusion within ice cores (Johnsen, 1977; Whillans and Grootes, 1985), with proxy seasonality (Jonkers and Kučera, 2015) potentially interacting with modulations of the seasonal cycle amplitude caused by slow orbital variations (Huybers and Wunsch, 2003; Laepple et al., 2011), and with uncertainties in the understanding of the climate–proxy relationship (including calibration errors; Tierney and Tingley, 2014). Further uncertainty sources may exist depending on the type of proxy used.

It turns out that a careful and systematic investigation of these reconstruction uncertainties is indispensable if we are to properly exploit the source of information contained in proxy archives, for such important issues like the estimation of the future evolution of natural and forced climate variability. Until now, however, reconstruction uncertainty estimates have often lacked the required accuracy (Lohmann et al., 2013; Reschke et al., 2019). In particular, one issue that deserves more detailed consideration is the timescale dependence of the reconstruction uncertainties (Amrhein, 2020). Although some of their sources like measurement errors will often be independent and, thus, uncorrelated between individual measurements, others like smoothing processes and orbital variations, in conjunction with proxy seasonality, have the potential to create serially correlated uncertainties (i.e., they are auto-correlated in the time domain). Thus, some uncertainty components may be described by white noise, while others may have the properties of red noise or an even more complex auto-correlation structure. The direct, and practically relevant, implication of this is the fact that, when averaging the proxy-based climate reconstruction over some time interval (e.g., by applying a moving average filter), the uncertainties may shrink at a different rate than if they were purely white noise.

One possibility to estimate the auto-correlation structure of reconstruction uncertainties is in the application of proxy forward models that generate proxy time series from climate (model) time series (see, for example, Evans et al., 2013; Dee et al., 2015; Dolman and Laepple, 2018). Specifically, the auto-correlation structure may then be inferred from ensembles of such simulated proxy time series. This approach is flexible regarding the complexity of the uncertainty-generating processes included in the model, but the insights gained from its application are limited by the fact that it represents a trial-and-error strategy. Moreover, the involved numerical simulations easily become computationally expensive. Therefore, it is useful and desirable to complement this with an alternative approach that allows for a systematic understanding of the auto-correlation structure of the reconstruction error components from an analytic point of view.

Accordingly, the aim of this paper is to provide a conceptual approach and, based thereon, an analytically derived method to estimate timescale-dependent reconstruction uncertainties for the example of sediment archives. Specifically, the method yields uncertainty estimates having been given a set of parameters that specify (i) the spectral structure of a supposed true climate signal, (ii) seasonal and orbital variations, (iii) proxy seasonality, (iv) bioturbation, (v) archive sampling parameters, (vi) sampling of a finite number of signal carriers, and (vii) uncorrelated measurement noise, and it takes into account the effects of spectral aliasing and leakage. The fact that archive smoothing is represented by bioturbation limits the validity of the method in its current form to proxy archives from sea and lake sediments. However, it has the potential to be generalized to other sedimentary archives such as ice cores by modifying the smoothing operator to represent isotopic diffusion.

The pivotal idea of our approach to address the timescale dependence of the uncertainty is in the derivation of its power spectrum since the spectrum is directly related (by the Wiener–Khintchine theorem; see Priestley, 1981) to the auto-correlation structure which, in turn, determines how the uncertainty scales with timescale (e.g., the length of an averaging interval). The convenience of the obtained mathematical expressions for the uncertainty power spectrum is twofold: (a) they can be used to acquire a qualitative understanding of the effects and relative importance of the various sources of uncertainty, and (b) they can serve to obtain quantitative uncertainty estimates for specific practical applications in paleoclimate science.

Part 1 of this study provides the theoretical basis of the uncertainty estimation method. In Sect. 2, the underlying reconstruction uncertainty model is defined in the time domain. Section 3 translates the model into the spectral domain by deriving the corresponding uncertainty power spectrum. Section 4 summarizes the results and demonstrates how timescale-dependent uncertainties can be obtained from the spectrum. The method and its limitations are discussed in Sect. 5, followed by the final conclusions in Sect. 6. Part 2 of this study, published separately (see Dolman et al., 2020), demonstrates the practical applicability of the method and also provides a software implementation for practical uncertainty estimation purposes, the so-called Proxy Spectral Error Model (PSEM).

2 Reconstruction uncertainty model

Before we can formulate our timescale-dependent uncertainty estimation method, we have to provide a precise definition of the underlying reconstruction uncertainty model, including our assumptions and simplifications that allow for an analytic treatment of the problem. Specifically, in order to define the uncertainty model, we need to

• suppose a structure of the true climate signal, which the final uncertainty estimates will be based upon, because some uncertainty components and their timescale dependence are subject to that structure;

• make simplifying assumptions regarding the archive formation, concerning proxy seasonality, the climate–proxy relationship, the sediment accumulation rate, and the effects of bioturbation mixing;

• specify the archive sampling and measurement procedure;

• define the reconstruction error as the difference between the obtained climate reconstruction and a suitable reference climate;

• define the reconstruction uncertainty in terms of the expected value of the squared reconstruction error.

Accordingly, the reconstruction uncertainty model can be thought of, conceptually, as an operator that takes as its arguments the supposed structure of the true climate signal and a set of parameters that appear in the mathematical formulation of the abovementioned assumptions. The remainder of this section is concerned with the details of the abovementioned five steps, including an explanation of the involved parameters. A complete list of the model parameters is provided in Appendix B. For possible sources and specific choices of parameter values, see Part 2 of this study (Dolman et al., 2020) and, in particular, its Table 1. Note that the reconstruction uncertainty model defined in this section is closely related to the proxy forward model of Dolman and Laepple (2018).

2.1 Climate signal

We assume that the supposed true climate signal consists of two components: a stochastic signal X(t) that represents the signal to be reconstructed from the proxy record and a deterministic signal Y(t) that represents the seasonal cycle, the amplitude of which is modulated by slow orbital variations. In addition, we make the simplifying assumption that X(t) and Y(t) are stochastically independent.

The stochastic signal X(t) is modeled as a zero-mean stochastically continuous stationary random process with infinite and continuous time parameter t, which has a purely continuous power spectrum (i.e., the spectrum has no discrete components). The actual structure of X(t) is to be specified in the spectral domain (see Sect. 3), and the uncertainty estimation method is constructed such that any spectral structure can be specified as long as it is consistent with the abovementioned properties of X(t). An illustration of one realization of such a random process is given by the gray–red line in Fig. 1b (obtained from a surrogate time series that obeys a simple power-law frequency scaling). It can be thought of as a time series of anomalies of a climate signal after the removal of the climatological seasonal cycle.

The deterministic signal Y(t) is modeled as a single harmonic oscillation that represents a simplified seasonal cycle, which is amplitude modulated by another single harmonic oscillation with a much longer period. Thus, Y(t) has a purely discrete power spectrum. Such a deterministic signal can be written as

$$\begin{array}{}\text{(1)}& Y\left(t\right)={\mathit{\sigma }}_{\mathrm{c}}\sqrt{\mathrm{2}}\cos \left({\mathit{\varphi }}_{\mathrm{c}}+\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{c}}t\right)\left[\mathrm{1}+{\mathit{\sigma }}_{\mathrm{a}}\sqrt{\mathrm{2}}\cos \left({\mathit{\varphi }}_{\mathrm{a}}+\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}t\right)\right],\end{array}$$

where ${\mathit{\nu }}_{\mathrm{c}}=(\mathrm{1}\mathrm{yr}{)}^{-\mathrm{1}}$ and ${\mathit{\nu }}_{\mathrm{a}}\ll (\mathrm{1}\mathrm{yr}{)}^{-\mathrm{1}}$ are the frequencies of the seasonal cycle and of its slow amplitude modulation, respectively, ${\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}$ are the corresponding variances of those oscillations, and ϕ_c and ϕ_a are their phases. The square bracket term represents the amplitude modulation factor that specifies the time-varying envelope of the seasonal cycle oscillation. Note that σ_c has the same units as X(t) and Y(t), whereas σ_a is dimensionless as it determines the fraction by which the amplitude of the seasonal cycle varies. In particular, ${\mathit{\sigma }}_{\mathrm{c}}\sqrt{\mathrm{2}}$ is the half-amplitude of the unmodulated seasonal cycle, and ${\mathit{\sigma }}_{\mathrm{a}}\sqrt{\mathrm{2}}$ is the fraction by which the seasonal cycle amplitude changes over an orbital modulation cycle. Furthermore, it is required that ${\mathit{\sigma }}_{\mathrm{a}}\sqrt{\mathrm{2}}<\mathrm{1}$, or equivalently ${\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}<\mathrm{1}/\mathrm{2}$, to avoid flipping seasons by a negative amplitude modulation factor (which would correspond to unrealistically strong effects of orbital variations). The deterministic signal is illustrated by the gray–red line in Fig. 1c. Note that only for the purpose of illustration has the modulation frequency been set to ${\mathit{\nu }}_{\mathrm{a}}=(\mathrm{130}\mathrm{yr}{)}^{-\mathrm{1}}$ in this figure, although a realistic value would be ${\mathit{\nu }}_{\mathrm{a}}=(\mathrm{23}\mathrm{kyr}{)}^{-\mathrm{1}}$, for example, if it were to represent an idealized planetary precession cycle.

Figure 1 Schematic illustration of the reconstruction uncertainty model. (a) Probability density function (as a function of time lag) that describes the combined effects of bioturbation (with timescale τ_b), of sediment sample thickness, and of proxy seasonality. (b) Stochastic component X(t) of the climate signal (gray line) with dates highlighted in red that fall into the proxy abundance window (of length τ_p) to reflect proxy seasonality. Blue rectangles indicate time intervals of length τ_s covered by the sediment slices. From each slice a finite number of signal carriers (N=3 in this example) are retrieved from random positions within the slice indicated by the blue dots, each of which carries the signal from the time at which it settled down on the surface of the sediment before it was mixed to its current position by bioturbation (indicated by black arrows for the central slice). Green squares indicate the reference climate signal obtained by averaging X(t) over intervals of length τ_r as indicated by the green line for the central point. (c) Same as (b) but for the deterministic component Y(t), which represents the amplitude modulated seasonal cycle. (d) Total reconstructed signal (blue) obtained by averaging over the N signal carriers from each slice, reference climate signal (green), and the difference between them (magenta) that represents the reconstruction error at a sampling interval Δt. Measurement errors are neglected in this illustration. Timescale parameters are set to τ_b=10 years, τ_s=5 years, τ_r=9 years, Δt=9 years, and ${\mathit{\tau }}_{\mathrm{p}}=\mathrm{1}/\mathrm{3}$ year.

2.2 Archive formation

To reflect proxy seasonality, we assume a seasonally confined time window during which the proxy is abundant. Thus, the climate signal and, in particular, the seasonal cycle are recorded only during those seasons. The length of this proxy abundance window is specified by the parameter τ_p, and the timing of the center of this window with respect to the seasonal cycle is specified by ϕ_c, as it appears in Eq. (1). Accordingly, if ϕ_c=0, then the abundance window is centered at the maximum of the seasonal cycle (i.e., the summer season if the climate signal X(t)+Y(t) represents temperature, for example), setting ${\mathit{\varphi }}_{\mathrm{c}}=\pm \mathit{\pi }/\mathrm{2}$ centers the window at either of the zero crossings (spring or autumn), ϕ_c=π at the minimum (winter), and likewise for all other phases. The seasonality parameters are required to fulfill the relations τ_p≤1 year and $-\mathit{\pi }<{\mathit{\varphi }}_{\mathrm{c}}\le \mathit{\pi }$. If τ_p=1 year, there is no seasonality and the parameter ϕ_c has no effect. Since in this formulation τ_p and ϕ_c are fixed, the abovementioned assumptions imply that we are neglecting any changes of proxy seasonality caused by, for example, habitat tracking. Specifically, there is no adaptation of proxy seasonality to changes in the seasonal cycle amplitude nor to variations of the stochastic component of the climate signal at any timescale. The effect of proxy seasonality defined in this way is illustrated in Fig. 1b and c by the red line segments which highlight that part of the signal that is recorded by the proxy. In this example, the proxy abundance window is set to cover the seasons around the maximum of the seasonal cycle.

In the following, we also neglect any uncertainties regarding the climate–proxy relationship, including calibration errors. Furthermore, we assume a known and constant sediment accumulation rate. Thus, we will treat all signals simply as a function of time and assume that the constant time–depth relationship is given as independent information.

Signal smoothing by sediment mixing caused by bioturbation is assumed to occur instantly and uniformly within the uppermost layer of the sediment. The thickness of this layer, the bioturbation depth, can be divided by the sediment accumulation rate to obtain the corresponding bioturbation timescale τ_b. Under the aforementioned assumptions, the effects of bioturbation can be described by a probability density function (PDF) of the form (Berger and Heath, 1968)

$$\begin{array}{}\text{(2)}& {\displaystyle }{f}_{\mathrm{b}}\left(\mathit{ϵ}\right)=\left\{\begin{array}{ll}{\mathit{\tau }}_{\mathrm{b}}^{-\mathrm{1}}\exp \left[\left(\mathit{ϵ}-{\mathit{\tau }}_{\mathrm{b}}\right)/{\mathit{\tau }}_{\mathrm{b}}\right]& \text{if}\mathit{ϵ}\le {\mathit{\tau }}_{\mathrm{b}}\\ \mathrm{0}& \text{otherwise}\end{array}\right.,\end{array}$$

where ϵ has units of time. This PDF specifies at which probability a single signal carrier, retrieved from the archive at position t=t₀, has settled down on the surface of the sediment and, thus, has recorded the climate signal at a given time $t=t_{\mathrm{0}}+\mathit{ϵ}$. Essentially, it states that a signal carrier retrieved at t=t₀ cannot have its origin at times later than $t=t_{\mathrm{0}}+{\mathit{\tau }}_{\mathrm{b}}$ but that it can have its origin arbitrarily far in the past relative to t₀, although with exponentially decreasing probability. Thus, ϵ can be interpreted as the timing error caused by bioturbation that is associated with the signal recorded by an individual signal carrier. In Fig. 1b and c, the effect of bioturbation is illustrated by the black arrows, indicating the net mixing paths (i.e., the timing errors ϵ) of three selected signal carriers (blue dots).

2.3 Sampling and measurement procedure

We assume that the archive is sampled by taking slices of sediment, the thickness of which corresponds to time intervals of length τ_s and which are taken at distances (measured from center to center) corresponding to a sampling interval Δt. This sampling procedure is illustrated in Fig. 1b and c by the blue rectangles, indicating individual sediment slices. The total length of the record is denoted by T. For mathematical reasons that become apparent in Sect. 3, it is required that T is a multiple of Δt, and that Δt and τ_s are multiples of 1 year. Setting τ_s<Δt corresponds to discontinuous sampling, τ_s=Δt to continuous sampling, and τ_s>Δt to sampling with overlap. The effects of these cases are discussed by Amrhein (2020).

Because signal carriers are retrieved from arbitrary positions within each slice, the effect of the sediment sample timescale τ_s can be described by convolving the bioturbation PDF f_b(ϵ) with a slice PDF that has the shape of a moving average window,

$$\begin{array}{}\text{(3)}& f_{\mathrm{s}}\left(\mathit{ϵ}\right)={\mathit{\tau }}_{\mathrm{s}}^{-\mathrm{1}}\mathrm{\Pi }(\mathit{ϵ};{\mathit{\tau }}_{\mathrm{s}}),\end{array}$$

and which is essentially blurring the edge of f_b(ϵ). The symbol Π(t;τ) denotes the rectangle function

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{\Pi }(t;\mathit{\tau })=\left\{\begin{array}{ll}\mathrm{1}& \text{if}-\mathit{\tau }/\mathrm{2}<t\le \mathit{\tau }/\mathrm{2}\\ \mathrm{0}& \text{otherwise}\end{array}\right..\end{array}$$

Thus, the PDF of the timing errors, which describes the combined effects of bioturbation and of sampling slices of sediment, can be written as

$$\begin{array}{}\text{(5)}& f_{\mathrm{bs}}\left(\mathit{ϵ}\right)=f_{\mathrm{s}}\left(\mathit{ϵ}\right)\ast f_{\mathrm{b}}\left(\mathit{ϵ}\right).\end{array}$$

Hence, if there were no bioturbation (τ_b→0) and if single signal carriers were retrieved from infinitesimally thin slices (τ_s→0), then this PDF would reduce to a Dirac delta function, f_bs(ϵ)→δ(ϵ), in which case the abovementioned sampling procedure would yield the discrete climate signal $X_n+Y_n=X\left(t_n\right)+Y\left(t_n\right)$, where t_n=nΔt (with $n=\mathrm{0},\pm \mathrm{1},\pm \mathrm{2},\mathrm{\dots }$). In the general case with τ_b>0 and τ_s>0, we can express the result of the sampling procedure as a discrete signal with jittered sampling,

$$\begin{array}{}\text{(6)}& X_n^{\left(j\right)}+Y_n^{\left(j\right)}=X\left(t_n+{\mathit{ϵ}}_n^{\left(j\right)}\right)+Y\left(t_n+{\mathit{ϵ}}_n^{\left(j\right)}\right),\end{array}$$

with ${\mathit{ϵ}}_n^{\left(j\right)}\sim f_{\mathrm{bs}}\left(\mathit{ϵ}\right)$ and ${\mathit{ϵ}}_n^{\left(j\right)}$ representing the sampling jitter and f_bs(ϵ) the jitter PDF. In the abovementioned terminology, ${\mathit{ϵ}}_n^{\left(j\right)}$ represents the timing error of a single signal carrier (labeled j) retrieved from a slice centered at t=t_n.

Finally, we need to include the effect of proxy seasonality as defined in the previous subsection. This is accomplished through multiplying f_bs(ϵ) by a proxy seasonality function p(ϵ) which is given by the convolution of the Dirac comb function $\text{III}(\mathit{ϵ};{\mathit{\nu }}_{\mathrm{c}}^{-\mathrm{1}})$ with the rectangle function $\left({\mathit{\tau }}_{\mathrm{p}}{\mathit{\nu }}_{\mathrm{c}}{)}^{-\mathrm{1}}\mathrm{\Pi }\right(\mathit{ϵ};{\mathit{\tau }}_{\mathrm{p}})$,

$$\begin{array}{}\text{(7)}& \begin{array}{rl}p\left(\mathit{ϵ}\right)& ={\left({\mathit{\tau }}_{\mathrm{p}}{\mathit{\nu }}_{\mathrm{c}}\right)}^{-\mathrm{1}}\mathrm{\Pi }\left(\mathit{ϵ};{\mathit{\tau }}_{\mathrm{p}}\right)\ast \text{III}\left(\mathit{ϵ};{\mathit{\nu }}_{\mathrm{c}}^{-\mathrm{1}}\right)\\ & ={\left({\mathit{\tau }}_{\mathrm{p}}{\mathit{\nu }}_{\mathrm{c}}\right)}^{-\mathrm{1}}\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{\Pi }\left(\mathit{ϵ}-k{\mathit{\nu }}_{\mathrm{c}}^{-\mathrm{1}};{\mathit{\tau }}_{\mathrm{p}}\right),\end{array}\end{array}$$

where the Dirac comb function III(t;τ) is defined as a series of Dirac delta functions δ(t):

$$\begin{array}{}\text{(8)}& \text{III}(t;\mathit{\tau })=\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}\mathit{\delta }(t-k\mathit{\tau }).\end{array}$$

It turns out that in the limit of vanishing proxy seasonality (τ_p→1 year), the proxy seasonality function becomes constantly one, p(ϵ)→1, whereas in the limit of maximum proxy seasonality (τ_p→0), the proxy seasonality function reduces to the Dirac comb function, $p\left(\mathit{ϵ}\right)\to \text{III}(\mathit{ϵ};{\mathit{\nu }}_{\mathrm{c}}^{-\mathrm{1}})$. From the above, the discrete climate signal, obtained from the sampling procedure, may still be written as in Eq. (6) but with the sampling jitter,

$$\begin{array}{}\text{(9)}& {\mathit{ϵ}}_n^{\left(j\right)}\sim p\left(\mathit{ϵ}\right)f_{\mathrm{bs}}\left(\mathit{ϵ}\right),\end{array}$$

now being drawn from the full jitter PDF, p(ϵ)f_bs(ϵ), which describes the combined effects of bioturbation, of sampling slices of sediment, and of proxy seasonality. The proof that the full jitter PDF defined in this way integrates to unity follows in Sect. 3. The structure of the full jitter PDF is illustrated in Fig. 1a.

In practice, a finite number, N≥1, of signal carriers is retrieved from each sediment slice rather than just a single signal carrier, and subsequently a single proxy measurement is performed in the laboratory on the collection of those N signal carriers, which represent an average proxy value. This can be expressed as

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{‾}}{X}}_n+{\stackrel{\mathrm{‾}}{Y}}_n={\displaystyle \frac{\mathrm{1}}{N}}\sum _{j=\mathrm{1}}^N[X_n^{\left(j\right)}+Y_n^{\left(j\right)}].\end{array}$$

In addition, we assume that the involved sampling jitter ${\mathit{ϵ}}_n^{\left(j\right)}$ is uncorrelated (i.e., white) in terms of both n and j:

$$\begin{array}{}\text{(11)}& \text{Cor}\left({\mathit{ϵ}}_n^{\left(j\right)},{\mathit{ϵ}}_{n^{\prime }}^{\left(j^{\prime }\right)}\right)=\mathrm{0}\text{if}n\ne n{}^{\prime }\text{or}j\ne j^{\prime },\end{array}$$

which reflects our assumption that the bioturbation mixing paths of the individual signal carriers within the sediment do not affect each other. Furthermore, it is required that ${\mathit{ϵ}}_n^{\left(j\right)}$ is a stationary process to reflect our assumptions of a fixed bioturbation depth and a constant sediment accumulation rate. Finally, we take ${\mathit{ϵ}}_n^{\left(j\right)}$ as independent of X(t) and Y(t), which corresponds to the assumption that bioturbation does not depend on the climate.

In general, each laboratory measurement is associated with a measurement error μ_n, the magnitude of which may be characterized in terms of its variance ${\mathit{\sigma }}_{\mathit{\mu }}^{\mathrm{2}}$. Thus, the final reconstruction time series is given by ${\stackrel{\mathrm{‾}}{X}}_n+{\stackrel{\mathrm{‾}}{Y}}_n+{\mathit{\mu }}_n$, although we will omit μ_n in the following as it is assumed to be white noise and can thus easily be added at the very end of the entire uncertainty estimation procedure.

2.4 Definition of reconstruction error

The reconstruction error can now be defined as the difference between the obtained climate reconstruction (Eq. 10) and a suitable reference climate:

$$\begin{array}{}\text{(12)}& {\stackrel{\mathrm{̃}}{X}}_n+{\stackrel{\mathrm{̃}}{Y}}_n=\stackrel{\mathrm{̃}}{X}\left(t_n\right)+\stackrel{\mathrm{̃}}{Y}\left(t_n\right),\end{array}$$

where

$$\begin{array}{}\text{(13)}& \stackrel{\mathrm{̃}}{X}\left(t\right)+\stackrel{\mathrm{̃}}{Y}\left(t\right)={\mathit{\tau }}_{\mathrm{r}}^{-\mathrm{1}}\mathrm{\Pi }(t;{\mathit{\tau }}_{\mathrm{r}})\ast \left[X\right(t)+Y(t\left)\right]\end{array}$$

is the supposed true climate signal smoothed with a moving average filter with timescale τ_r, which is then subsampled at the same discrete times t_n. Here we require that τ_r is a multiple of 1 year such that ${\stackrel{\mathrm{̃}}{Y}}_n=\mathrm{0}$ because it is then an average over a number of complete seasonal cycles. Thus, we obtain the reconstruction error time series

$$\begin{array}{}\text{(14)}& E_n=E_{X,n}+E_{Y,n}\end{array}$$

with the error components

$$\begin{array}{}\text{(15)}& E_{X,n}={\stackrel{\mathrm{‾}}{X}}_n-{\stackrel{\mathrm{̃}}{X}}_n\text{and}{E}_{Y,n}={\stackrel{\mathrm{‾}}{Y}}_n.\end{array}$$

An example of one realization of the discrete climate reconstruction (Eq. 10), reference climate (Eq. 12), and reconstruction error time series (Eq. 14) is illustrated in Fig. 1d.

2.5 Definition of reconstruction uncertainty

The reconstruction error E_n refers only to a single realization of the stochastic processes X and ϵ, which are specified in terms of their power spectral density and PDF, respectively. Thus, to obtain a suitable measure of the reconstruction uncertainty that characterizes the magnitude of possible errors under the specified stochastic properties of X and ϵ, we define the root mean square (rms) reconstruction error ℰ_n by

$$\begin{array}{}\text{(16)}& {\mathcal{E}}_n^{\mathrm{2}}=\langle \langle E_n^{\mathrm{2}}{\rangle }_X{\rangle }_{\mathit{ϵ}},\end{array}$$

where 〈⋅〉_X and 〈⋅〉_ϵ denote the expected value operators with respect to X and ϵ, respectively. Then substitution from Eq. (14) yields

$$\begin{array}{}\text{(17)}& {\mathcal{E}}_n^{\mathrm{2}}=\langle \langle (E_{X,n}+E_{Y,n}{)}^{\mathrm{2}}{\rangle }_X{\rangle }_{\mathit{ϵ}}=\langle \langle E_{X,n}^{\mathrm{2}}{\rangle }_X{\rangle }_{\mathit{ϵ}}+\langle E_{Y,n}^{\mathrm{2}}{\rangle }_{\mathit{ϵ}},\end{array}$$

because X and Y are assumed to be independent and because X is a zero-mean process, which is to say 〈X〉_X=0 and, thus, $\langle \langle E_{X,n}{\rangle }_X{\rangle }_{\mathit{ϵ}}=\mathrm{0}$.

As will be shown in Sect. 3, E_X,n can be decomposed into two uncorrelated zero-mean stationary components: $E_{X,n}=F_{X,n}+W_{X,n}$ such that F_X,n can be expressed as the result obtained by bandpass filtering the signal X(t) in time and then subsampling it at the discrete times t_n, and W_X,n is a white noise process. Furthermore, it will be shown that E_Y,n can be decomposed into two uncorrelated and generally nonstationary components: $E_{Y,n}=F_{Y,n}+W_{Y,n}$ such that F_Y,n can be expressed as the result obtained by filtering and then subsampling the signal Y(t), and W_Y,n is a zero-mean white noise process. Thus, we can write

$$\begin{array}{}\text{(18)}& {\mathcal{E}}_n^{\mathrm{2}}=\langle \langle F_{X,n}^{\mathrm{2}}{\rangle }_X{\rangle }_{\mathit{ϵ}}+\langle \langle W_{X,n}^{\mathrm{2}}{\rangle }_X{\rangle }_{\mathit{ϵ}}+F_{Y,n}^{\mathrm{2}}+\langle W_{Y,n}^{\mathrm{2}}{\rangle }_{\mathit{ϵ}},\end{array}$$

where F_Y,n is a deterministic signal.

In addition to the uncertainty caused by the stochasticity of X and ϵ, we can in principle equip any of the model parameters with its own uncertainty and investigate how this contributes to the obtained reconstruction uncertainty. In the following, we apply this procedure to the seasonal phase ϕ_c as the seasonal timing of the proxy abundance is often a poorly constrained parameter. For this purpose, we need to specify a corresponding PDF of the seasonal phase. For simplicity, we choose the wrapped uniform PDF

$$\begin{array}{}\text{(19)}& f_{{\mathit{\varphi }}_{\mathrm{c}}}\left({\mathit{\varphi }}_{\mathrm{c}}\right)=\sum _{k=-\mathrm{1}}^{\mathrm{1}}{\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}^{-\mathrm{1}}\mathrm{\Pi }\left({\mathit{\varphi }}_{\mathrm{c}}-\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}+\mathrm{2}\mathit{\pi }k;{\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}\right),\end{array}$$

with $-\mathit{\pi }<{\mathit{\varphi }}_{\mathrm{c}}\le \mathit{\pi }$, the expected seasonal phase $-\mathit{\pi }<\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}\le \mathit{\pi }$, and the seasonal phase uncertainty $\mathrm{0}\le {\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}<\mathrm{2}\mathit{\pi }$. Note that setting ${\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}=\mathrm{2}\mathit{\pi }$ does not imply vanishing proxy seasonality (as this is expressed by setting τ_p=1 year), but it merely means that the seasonal timing of the proxy abundance window is completely unknown. The model parameters to be specified are now $\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}$ and ${\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}$ (see also Appendix B) rather than the single parameter ϕ_c, which is treated as unknown according to $f_{{\mathit{\varphi }}_{\mathrm{c}}}\left({\mathit{\varphi }}_{\mathrm{c}}\right)$. Now, to include the effect on the reconstruction uncertainty, we redefine the rms reconstruction error ℰ_n by applying the additional expected value operator $\langle \cdot {\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}$ with respect to ϕ_c to the right-hand side of Eq. (16) or likewise Eq. (18). Since X does not depend on ϕ_c, we obtain

$$\begin{array}{}\text{(20)}& {\displaystyle }{\mathcal{E}}_n^{\mathrm{2}}& {\displaystyle }=\langle \langle \langle E_n^{\mathrm{2}}{\rangle }_X{\rangle }_{\mathit{ϵ}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}\text{(21)}& {\displaystyle }& {\displaystyle }=\langle \langle F_{X,n}^{\mathrm{2}}{\rangle }_X{\rangle }_{\mathit{ϵ}}+\langle \langle W_{X,n}^{\mathrm{2}}{\rangle }_X{\rangle }_{\mathit{ϵ}}+\langle F_{Y,n}^{\mathrm{2}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}+\langle \langle W_{Y,n}^{\mathrm{2}}{\rangle }_{\mathit{ϵ}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}.\end{array}$$

Hence, by noting that $\langle F_{Y,n}^{\mathrm{2}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}=\langle F_{Y,n}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}^{\mathrm{2}}+{\text{Var}}_{\left({\mathit{\varphi }}_{\mathrm{c}}\right)}\left(F_{Y,n}\right)$, we can finally write

$$\begin{array}{}\text{(22)}& {\mathcal{E}}_n^{\mathrm{2}}={\mathcal{B}}_n^{\mathrm{2}}+{\mathcal{U}}_n^{\mathrm{2}},\end{array}$$

with the squared reconstruction bias,

$$\begin{array}{}\text{(23)}& {\mathcal{B}}_n^{\mathrm{2}}=\langle F_{Y,n}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}^{\mathrm{2}},\end{array}$$

and the squared reconstruction uncertainty,

$$\begin{array}{}\text{(24)}& {\mathcal{U}}_n^{\mathrm{2}}={\mathcal{U}}_{\left(\mathrm{1}\right)}^{\mathrm{2}}+{\mathcal{U}}_{\left(\mathrm{2}\right)}^{\mathrm{2}}+{\mathcal{U}}_{\left(\mathrm{3}\right),n}^{\mathrm{2}}+{\mathcal{U}}_{\left(\mathrm{4}\right),n}^{\mathrm{2}},\end{array}$$

the components of which are given by

$$\begin{array}{}\text{(25)}& \begin{array}{rl}& {\mathcal{U}}_{\left(\mathrm{1}\right)}^{\mathrm{2}}={\text{Var}}_{(X,\mathit{ϵ})}\left(F_{X,n}\right),{\mathcal{U}}_{\left(\mathrm{2}\right)}^{\mathrm{2}}={\text{Var}}_{(X,\mathit{ϵ})}\left(W_{X,n}\right),\\ & {\mathcal{U}}_{\left(\mathrm{3}\right),n}^{\mathrm{2}}={\text{Var}}_{\left({\mathit{\varphi }}_{\mathrm{c}}\right)}\left(F_{Y,n}\right),{\mathcal{U}}_{\left(\mathrm{4}\right),n}^{\mathrm{2}}=\langle {\text{Var}}_{\left(\mathit{ϵ}\right)}\left(W_{Y,n}\right){\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}.\end{array}\end{array}$$

Note that from ${\mathcal{U}}_{\left(\mathrm{1}\right)}^{\mathrm{2}}$ and ${\mathcal{U}}_{\left(\mathrm{2}\right)}^{\mathrm{2}}$ the time index n has been dropped to indicate the stationarity of these uncertainty components. It turns out that ${\mathcal{E}}_n^{\mathrm{2}}$ represents the expected power of the reconstruction error at a given time t_n, which, according to Eq. (22), is decomposed into the power ${\mathcal{B}}_n^{\mathrm{2}}$ contained in the reconstruction bias and the variance ${\mathcal{U}}_n^{\mathrm{2}}$ that quantifies the scatter around the bias.

The individual components are to be interpreted as follows: the component 𝒰₍₁₎ quantifies the reconstruction uncertainty that arises from the difference between (i) the smoothing effect on X(t) caused by bioturbation and by sampling from slices of sediment and (ii) the smoothing effect on X(t) caused by the moving average window used to obtain the reference climate. Since the two smoothing effects represent low-pass filters with different cut-off frequencies, they act together as a bandpass filter on X(t) (as shown by Amrhein, 2020). This uncertainty component represents the total smoothing effect in the limit of infinitely many signal carriers being retrieved from each slice of sediment (N→∞). If only a finite number of signal carriers is retrieved from each slice, there is an additional residual that is not averaged out in this case. This residual is quantified by the component 𝒰₍₂₎. Likewise, the component 𝒰_(4),n quantifies the additional residual that arises from sampling only a finite number of signal carriers but now pertains to the deterministic signal Y(t). This residual component also depends on the timing of the seasonal proxy abundance and its uncertainty, as specified by Eq. (19), because of the non-linear relation between the variance, aliased from the seasonal cycle, and the seasonal timing. Because the seasonal cycle amplitude is modulated over time by orbital variations, this uncertainty component is nonstationary. On the other hand, in the limit of infinitely many signal carriers, the smoothing effects on Y(t) leave nothing but a deterministic bias that obtains its only uncertainty, quantified by 𝒰_(3),n, from the seasonal timing uncertainty. Finally, when averaging this bias across all possible seasonal timings that are allowed according to Eq. (19), a purely deterministic error component is obtained which is quantified by the reconstruction bias ℬ_n. Further clarification of these interpretations will emerge in Sect. 3.

Now, in order to formulate our timescale-dependent uncertainty estimation method, we need a spectral representation of the expected power ${\mathcal{E}}_n^{\mathrm{2}}$. We will achieve this by deriving the power spectral density of the reconstruction error E_n separately for its individual components to obtain spectral representations of the squared reconstruction uncertainty components ${\mathcal{U}}_{\left(\mathrm{1}\right)}^{\mathrm{2}}$, ${\mathcal{U}}_{\left(\mathrm{2}\right)}^{\mathrm{2}}$, ${\mathcal{U}}_{\left(\mathrm{3}\right),n}^{\mathrm{2}}$, ${\mathcal{U}}_{\left(\mathrm{4}\right),n}^{\mathrm{2}}$, and of the squared reconstruction bias ${\mathcal{B}}_n^{\mathrm{2}}$. This task is addressed in the following section. The reader who does not intend to follow the entire derivation may proceed directly to Sect. 4, which summarizes the main results of Sect. 3 and then illustrates the method based thereon for estimating timescale-dependent reconstruction uncertainties.

3 Spectral representation of reconstruction uncertainty

The reconstruction uncertainty model, defined in the previous section, is now translated into the spectral domain. Since the two components of the supposed true climate signal have different properties – in the sense that X(t) is a stationary random process with a continuous power spectrum, whereas Y(t) is an nonstationary deterministic signal with a discrete power spectrum – the two components require separate mathematical treatment. Accordingly, the derivation of the power spectral density of the reconstruction error E_n is accomplished separately in the following two subsections for the components of E_X,n that are based on the stochastic signal X(t) and for those of E_Y,n that are based on the deterministic signal Y(t).

3.1 Stochastic signal: continuous climate spectrum

A spectral representation of the stochastic signal component X(t), with infinite and continuous time parameter t, is given by the Riemann–Stieltjes integral (see Priestley, 1981, Sect. 4.11):

$$\begin{array}{}\text{(26)}& X\left(t\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }t}\mathrm{d}Z\left(\mathit{\nu }\right),\end{array}$$

where Z(ν) is a complex-valued stochastic process such that the power spectral density of X(t) is given by

$$\begin{array}{}\text{(27)}& S_X\left(\mathit{\nu }\right)=\langle \left|\mathrm{d}Z\right(\mathit{\nu }){|}^{\mathrm{2}}{\rangle }_X/\mathrm{d}\mathit{\nu },\end{array}$$

and where the dZ(ν) are zero-mean, orthogonal increments. Note that a conventional Fourier representation of X(t) does not exist because of the stochastic nature of the signal and that it is dZ(ν)∕dν rather than Z(ν) which formally plays the role of the Fourier transform in the abovementioned representation (Priestley, 1981). Likewise, the signal X_n=X(t_n), sampled at the discrete times t_n=nΔt, has the spectral representation

$$\begin{array}{}\text{(28)}& X_n=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{t}_n}\mathrm{d}Z\left(\mathit{\nu }\right),\end{array}$$

and the signal with jittered sampling, $X_n^{\left(j\right)}=X(t_n+{\mathit{ϵ}}_n^{\left(j\right)})$, as defined by (6), can be expressed as (Moore and Thomson, 1991)

$$\begin{array}{}\text{(29)}& X_n^{\left(j\right)}=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }(t_n+{\mathit{ϵ}}_n^{\left(j\right)})}\mathrm{d}Z\left(\mathit{\nu }\right).\end{array}$$

Following the approach of Balakrishnan (1962), we consider its auto-covariance function $\langle \langle {X_n^{\left(j\right)}}^{\star }{X}_{n^{\prime }}^{\left(j\right)}{\rangle }_X{\rangle }_{\mathit{ϵ}}$, where $(\cdot {)}^{\star }$ denotes the complex conjugate. By substitution from Eq. (29) and expressing the product of integrals as a double integral, we obtain (see Priestley, 1981, 249–250, where the same is shown for the case without sampling jitter)

$$\begin{array}{}\text{(30)}& {\displaystyle }\begin{array}{rl}& \langle \langle {X_n^{\left(j\right)}}^{\star }{X}_{n^{\prime }}^{\left(j\right)}{\rangle }_X{\rangle }_{\mathit{ϵ}}\\ & ={〈{〈\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{-i\mathrm{2}\mathit{\pi }\mathit{\nu }(t_n+{\mathit{ϵ}}_n^{\left(j\right)})}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{}^{\prime }(t_{n{}^{\prime }}+{\mathit{ϵ}}_{n^{\prime }}^{\left(j\right)})}\mathrm{d}{Z}^{\star }\left(\mathit{\nu }\right)\mathrm{d}Z\left(\mathit{\nu }{}^{\prime }\right)〉}_X〉}_{\mathit{ϵ}}\end{array}\text{(31)}& {\displaystyle }\begin{array}{rl}& =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }{t}_n({\mathit{\nu }}^{\prime }-\mathit{\nu })}{e}^{i\mathrm{2}\mathit{\pi }{\mathit{\nu }}^{\prime }(t_{n{}^{\prime }}-t_n)}\langle e^{i\mathrm{2}\mathit{\pi }({\mathit{\nu }}^{\prime }{\mathit{ϵ}}_{n^{\prime }}^{\left(j\right)}-\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)})}{\rangle }_{\mathit{ϵ}}\\ & \langle \mathrm{d}{Z}^{\star }\left(\mathit{\nu }\right)\mathrm{d}Z\left({\mathit{\nu }}^{\prime }\right){\rangle }_X,\end{array}\end{array}$$

where we have used the independence of ${\mathit{ϵ}}_n^{\left(j\right)}$ and X(t). Now, from the orthogonality of the dZ(ν), it follows that $\langle \mathrm{d}{Z}^{\star }\left(\mathit{\nu }\right)\mathrm{d}Z\left({\mathit{\nu }}^{\prime }\right){\rangle }_X=\mathrm{0}$ whenever $\mathit{\nu }\ne {\mathit{\nu }}^{\prime }$. Thus, the contribution to the integral (Eq. 31) is nonzero only for $\mathit{\nu }={\mathit{\nu }}^{\prime }$, and the auto-covariance function can then be expressed by the single integral (also using Eq. 27)

$$\begin{array}{}\text{(32)}& \langle \langle {X_n^{\left(j\right)}}^{\star }{X}_{n^{\prime }}^{\left(j\right)}{\rangle }_X{\rangle }_{\mathit{ϵ}}=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }(t_{n{}^{\prime }}-t_n)}{C}_{n,n{}^{\prime }}(-\mathit{\nu },\mathit{\nu })S_X\left(\mathit{\nu }\right)\mathrm{d}\mathit{\nu }\end{array}$$

with the characteristic function

$$\begin{array}{}\text{(33)}& C_{n,n{}^{\prime }}({\mathit{\nu }}_{\mathrm{1}},{\mathit{\nu }}_{\mathrm{2}})=\langle e^{i\mathrm{2}\mathit{\pi }({\mathit{\nu }}_{\mathrm{1}}{\mathit{ϵ}}_n^{\left(j\right)}+{\mathit{\nu }}_{\mathrm{2}}{\mathit{ϵ}}_{n^{\prime }}^{\left(j\right)})}{\rangle }_{\mathit{ϵ}}.\end{array}$$

Note that without sampling jitter (i.e., with ${\mathit{ϵ}}_n^{\left(j\right)}=\mathrm{0}$), expression (32) reduces to the Wiener–Khintchine theorem (see Priestley, 1981, for example), which states that the auto-covariance function of a signal and its power spectral density are a Fourier transform pair. Because ${\mathit{ϵ}}_n^{\left(j\right)}$ is white, we have $\langle e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }({\mathit{ϵ}}_{n^{\prime }}^{\left(j\right)}-{\mathit{ϵ}}_n^{\left(j\right)})}{\rangle }_{\mathit{ϵ}}=\langle e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_{n^{\prime }}^{\left(j\right)}}{\rangle }_{\mathit{ϵ}}\langle e^{-i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)}}{\rangle }_{\mathit{ϵ}}$ if $n\ne n{}^{\prime }$ and, thus,

$$\begin{array}{}\text{(34)}& {\displaystyle }{C}_{n,n{}^{\prime }}(-\mathit{\nu },\mathit{\nu })=\left\{\begin{array}{ll}\mathrm{1}& \text{if}n=n{}^{\prime }\\ \left|C\right(\mathit{\nu }){|}^{\mathrm{2}}& \text{if}n\ne n{}^{\prime }\end{array}\right.,\end{array}$$

where

$$\begin{array}{}\text{(35)}& C\left(\mathit{\nu }\right)=\langle e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)}}{\rangle }_{\mathit{ϵ}}\end{array}$$

is the characteristic function (or the complex conjugate of the Fourier transform) of the jitter PDF, p(ϵ)f_bs(ϵ), since using the definition of the expected value yields

$$\begin{array}{}\text{(36)}& {\displaystyle }C\left(\mathit{\nu }\right)& {\displaystyle }=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }\mathit{ϵ}}p\left(\mathit{ϵ}\right)f_{\mathrm{bs}}\left(\mathit{ϵ}\right)\mathrm{d}\mathit{ϵ}\text{(37)}& {\displaystyle }& {\displaystyle }=\widehat{p}\left(\mathit{\nu }\right)\ast {\widehat{f}}_{\mathrm{bs}}^{\star }\left(\mathit{\nu }\right)\text{(38)}& {\displaystyle }& {\displaystyle }=\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}\text{sinc}\left(k{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right){\widehat{f}}_{\mathrm{bs}}^{\star }(\mathit{\nu }+k{\mathit{\nu }}_{\mathrm{c}}).\end{array}$$

Here $\widehat{x}\left(\mathit{\nu }\right)$ denotes the Fourier transform of a function x(t), and we are using the fact that the Fourier transform of Π(t;τ) is given by τsinc(ντ) and the Fourier transform of III(t;τ) by ${\mathit{\tau }}^{-\mathrm{1}}\text{III}(\mathit{\nu };{\mathit{\tau }}^{-\mathrm{1}})$, and f_bs(ϵ) and p(ϵ) are defined by Eqs. (5) and (7), respectively. The cardinal sine function is defined as

$$\begin{array}{}\text{(39)}& {\displaystyle }\text{sinc}\left(x\right)=\left\{\begin{array}{ll}\mathrm{1}& \text{if}x=\mathrm{0}\\ \sin \left(\mathit{\pi }x\right)/\left(\mathit{\pi }x\right)& \text{if}x\ne \mathrm{0}\end{array}\right.,\end{array}$$

and, for the steps from Eq. (36) to Eq. (38), the convolution theorem is used. Expression (38) represents a series of amplitude-modulated and shifted versions of the function ${\widehat{f}}_{\mathrm{bs}}^{\star }\left(\mathit{\nu }\right)$, which is obtained by taking the Fourier transform of Eq. (5):

$$\begin{array}{}\text{(40)}& {\widehat{f}}_{\mathrm{bs}}^{\star }\left(\mathit{\nu }\right)=\text{sinc}\left(\mathit{\nu }{\mathit{\tau }}_{\mathrm{s}}\right)[\mathrm{1}+i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{\tau }}_{\mathrm{b}}{]}^{-\mathrm{1}}\exp (i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{\tau }}_{\mathrm{b}}),\end{array}$$

and, thus,

$$\begin{array}{}\text{(41)}& \left|{\widehat{f}}_{\mathrm{bs}}^{\star }\right(\mathit{\nu }){|}^{\mathrm{2}}={\text{sinc}}^{\mathrm{2}}(\mathit{\nu }{\mathit{\tau }}_{\mathrm{s}}\left)\right[\mathrm{1}+(\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{\tau }}_{\mathrm{b}}{)}^{\mathrm{2}}{]}^{-\mathrm{1}},\end{array}$$

which is the product of a squared sinc function and a Lorentzian function. In the following, we assume that $[max({\mathit{\tau }}_{\mathrm{b}},{\mathit{\tau }}_{\mathrm{s}}){]}^{-\mathrm{1}}\ll {\mathit{\nu }}_{\mathrm{c}}$ such that the characteristic width of the functions ${\widehat{f}}_{\mathrm{bs}}^{\star }(\mathit{\nu }+k{\mathit{\nu }}_{\mathrm{c}})$ is much less than the shift increment ν_c. These functions then have negligible overlap, and we can write

$$\begin{array}{}\text{(42)}& {\displaystyle }\left|C\right(\mathit{\nu }){|}^{\mathrm{2}}=\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}{\text{sinc}}^{\mathrm{2}}(k{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\left)\right|{\widehat{f}}_{\mathrm{bs}}^{\star }(\mathit{\nu }+k{\mathit{\nu }}_{\mathrm{c}}){|}^{\mathrm{2}}.\end{array}$$

The structure of $\left|C\right(\mathit{\nu }){|}^{\mathrm{2}}$ is illustrated in Fig. 2a, and it represents the squared modulus of the Fourier transform of the jitter PDF shown in Fig. 1a. Since it is shown for τ_b=10 years and τ_s=5 years, we have $[max({\mathit{\tau }}_{\mathrm{b}},{\mathit{\tau }}_{\mathrm{s}}){]}^{-\mathrm{1}}=\mathrm{1}/(\mathrm{10}\mathrm{years})\ll {\mathit{\nu }}_{\mathrm{c}}$, and, thus, the $\left|{\widehat{f}}_{\mathrm{bs}}^{\star }\right(\mathit{\nu }+k{\mathit{\nu }}_{\mathrm{c}}){|}^{\mathrm{2}}$ peaks are well separated along the frequency axis.

The proof that the jitter PDF p(ϵ)f_bs(ϵ), as defined by Eq. (9), does indeed integrate to unity is equivalent to showing that C(0)=1, as can be seen from Eq. (36) with ν=0. To demonstrate this, we evaluate C(0) using Eq. (38), noting (i) that the term with k=0 is equal to one at ν=0 because ${\widehat{f}}_{\mathrm{bs}}^{\star }\left(\mathrm{0}\right)=\mathrm{1}$ according to Eq. (40) and (ii) that the remaining terms with k≠0 are all equal to 0 at ν=0 because ${\widehat{f}}_{\mathrm{bs}}^{\star }\left(k{\mathit{\nu }}_{\mathrm{c}}\right)=\mathrm{0}$ since sinc(kν_cτ_s)=0 according to the requirement that τ_s be a multiple of 1 year (see Sect. 2.3).

Figure 2 (a) Squared modulus of the characteristic function (or of the Fourier transform) of the jitter PDF, $\left|C\right(\mathit{\nu }){|}^{\mathrm{2}}$, given by Eq. (42) (solid line), for the same parameters as used in Fig. 1 (i.e., τ_b=10 years, τ_s=5 years, ${\mathit{\tau }}_{\mathrm{p}}=\mathrm{1}/\mathrm{3}$ year) and the envelope function sinc²(ντ_p) (dashed line). (b) Squared modulus of the error transfer function, $\left|C\right(\mathit{\nu })-\text{sinc}(\mathit{\nu }{\mathit{\tau }}_{\mathrm{r}}){|}^{\mathrm{2}}$, as it appears in Eq. (58) and with τ_r=9 years (dashed line as in a). The frequency axis is normalized by the seasonal cycle frequency ${\mathit{\nu }}_{\mathrm{c}}=(\mathrm{1}\mathrm{yr}{)}^{-\mathrm{1}}$.

To obtain the power spectral density of the reconstruction error components of E_X,n, we rewrite the integrand of Eq. (29) as $e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{t}_n}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)}}$ and split the jitter factor $e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)}}$ into its expected value C(ν) and the deviation thereof, $e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)}}-C\left(\mathit{\nu }\right)$, as in Moore and Thomson (1991). Then we can decompose $X_n^{\left(j\right)}$ as

$$\begin{array}{}\text{(43)}& X_n^{\left(j\right)}=U_n+V_n^{\left(j\right)},\end{array}$$

with the components

$$\begin{array}{}\text{(44)}& U_n=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{t}_n}C\left(\mathit{\nu }\right)\mathrm{d}Z\left(\mathit{\nu }\right)\end{array}$$

and

$$\begin{array}{}\text{(45)}& V_n^{\left(j\right)}=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{t}_n}[e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)}}-C(\mathit{\nu }\left)\right]\mathrm{d}Z\left(\mathit{\nu }\right).\end{array}$$

From this we obtain, by analogy with the steps from Eq. (30) to Eq. (32), the auto-covariance functions of U_n and $V_n^{\left(j\right)}$, as well as their cross-covariance function:

$$\begin{array}{}\text{(46)}& {\displaystyle }\langle U_n^{\star }{U}_{n{}^{\prime }}{\rangle }_X=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }(t_{n{}^{\prime }}-t_n)}\left|C\right(\mathit{\nu }\left){|}^{\mathrm{2}}{S}_X\right(\mathit{\nu })\mathrm{d}\mathit{\nu },\text{(47)}& {\displaystyle }\begin{array}{rl}& \langle \langle {V_n^{\left(j\right)}}^{\star }{V}_{n^{\prime }}^{\left(j\right)}{\rangle }_X{\rangle }_{\mathit{ϵ}}\\ & =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }(t_{n{}^{\prime }}-t_n)}\left[C_{n,n{}^{\prime }}\right(-\mathit{\nu },\mathit{\nu })-|C\left(\mathit{\nu }\right){|}^{\mathrm{2}}\left]S_X\right(\mathit{\nu })\mathrm{d}\mathit{\nu },\end{array}\text{(48)}& {\displaystyle }\begin{array}{rl}& \langle \langle U_n^{\star }{V}_{n^{\prime }}^{\left(j\right)}{\rangle }_X{\rangle }_{\mathit{ϵ}}\\ & =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }(t_{n{}^{\prime }}-t_n)}{C}^{\star }\left(\mathit{\nu }\right)\langle e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)}}-C\left(\mathit{\nu }\right){\rangle }_{\mathit{ϵ}}{S}_X\left(\mathit{\nu }\right)\mathrm{d}\mathit{\nu }.\end{array}\end{array}$$

Since the term $\langle e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)}}-C\left(\mathit{\nu }\right){\rangle }_{\mathit{ϵ}}$ in Eq. (48) is 0, the cross power spectral density of U_n and $V_n^{\left(j\right)}$ vanishes at all frequencies, and, thus, the two processes are uncorrelated. Accordingly, the sum of their auto-covariance functions, Eqs. (46) and (47), equals the auto-covariance function of $X_n^{\left(j\right)}$, given by Eq. (32). Furthermore, Eq. (32) with $n=n{}^{\prime }$ shows that ${\text{Var}}_{(X,\mathit{ϵ})}\left(X_n^{\left(j\right)}\right)={\text{Var}}_{\left(X\right)}\left(X\right)$.

Note that the square bracket term in Eq. (47) represents the auto-covariance function of the jitter factor $e^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{\mathit{ϵ}}_n^{\left(j\right)}}$, and Eq. (34) implies that

$$\begin{array}{}\text{(49)}& \begin{array}{rl}& \left[C_{n,n{}^{\prime }}\right(-\mathit{\nu },\mathit{\nu })-|C\left(\mathit{\nu }\right){|}^{\mathrm{2}}]\\ & =\left\{\begin{array}{ll}\mathrm{1}-\left|C\right(\mathit{\nu }){|}^{\mathrm{2}}& \text{if}n=n{}^{\prime }\\ \mathrm{0}& \text{if}n\ne n{}^{\prime }\end{array}\right..\end{array}\end{array}$$

Thus, the auto-covariance function of $V_n^{\left(j\right)}$ is nonzero only at lag zero ($n=n{}^{\prime }$) and zero at all other lags ($n\ne n{}^{\prime }$), and, hence, $V_n^{\left(j\right)}$ is a white noise process. On the other hand, U_n can be seen as the result of linearly filtering the signal X(t) with the jitter PDF, p(ϵ)f_bs(ϵ), and then subsampling it at the discrete times t_n, with $\left|C\right(\mathit{\nu }){|}^{\mathrm{2}}$ in Eq. (46) being interpreted as the squared modulus of the spectral transfer function. Since neither the linear filter nor the subsampling alters the expected value, we have $\langle U_n{\rangle }_X=\langle X_n^{\left(j\right)}{\rangle }_X=\langle X{\rangle }_X=\mathrm{0}$, and from Eq. (43) it follows that also $\langle V_n^{\left(j\right)}{\rangle }_X=\langle X_n^{\left(j\right)}-U_n{\rangle }_X=\mathrm{0}$. Finally, the stationarity of X and ϵ implies that U_n and $V_n^{\left(j\right)}$ are stationary. The structure of the auto-covariance function of $X_n^{\left(j\right)}$, given by Eq. (32), is illustrated schematically in Fig. 3, highlighting its decomposition into the respective contributions from U_n and $V_n^{\left(j\right)}$. In particular, it turns out that the magnitude of the variance of the white noise component $V_n^{\left(j\right)}$ is obtained by extrapolating the auto-covariance function from nonzero lags towards lag zero. This separates the full variance into two components (indicated in the figure by the transition in color at lag zero) such that, by setting $n=n{}^{\prime }$ in Eqs. (46), (47), and (49) and substitution from Eq. (49) into Eq. (47),

$$\begin{array}{}\text{(50)}& {\displaystyle }{\text{Var}}_{(X,\mathit{ϵ})}\left(V_n^{\left(j\right)}\right)={\text{Var}}_{\left(X\right)}\left(X\right)-{\text{Var}}_{\left(X\right)}\left(U_n\right).\end{array}$$

This is the key idea of the approach of Balakrishnan (1962), and we will return to this idea in Sect. 3.2 in the context of the deterministic signal component Y(t).

Figure 3 Schematic illustration of the auto-covariance function of the discrete process $X_n^{\left(j\right)}$ (black circles), as given by Eq. (32), normalized by the variance of X(t) and the auto-covariance contribution from U_n (red lines), as given by Eq. (46), and from $V_n^{\left(j\right)}$(green lines), as given by Eq. (47), with green dots indicating zero contribution.

With these properties of the abovementioned components U_n and $V_n^{\left(j\right)}$, we can now rewrite the error component E_X,n, defined by Eq. (15), also using the X component of Eq. (10):

$$\begin{array}{}\text{(51)}& {\displaystyle }{E}_{X,n}=F_{X,n}+W_{X,n},\end{array}$$

with

$$\begin{array}{}\text{(52)}& {\displaystyle }{F}_{X,n}=U_n-{\stackrel{\mathrm{̃}}{X}}_n\end{array}$$

and

$$\begin{array}{}\text{(53)}& {\displaystyle }{W}_{X,n}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{j=\mathrm{1}}^N{V}_n^{\left(j\right)},\end{array}$$

where F_X,n and W_X,n represent the components of E_X,n explained in Sect. 2.5, which is to say a component obtained by filtering the signal X and a white noise component.

By analogy with Eq. (28), a spectral representation of the X component ${\stackrel{\mathrm{̃}}{X}}_n$ of the reference climate signal, defined by Eqs. (12) and (13), is given by

$$\begin{array}{}\text{(54)}& {\displaystyle }{\stackrel{\mathrm{̃}}{X}}_n=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }{t}_n}\text{sinc}\left(\mathit{\nu }{\mathit{\tau }}_{\mathrm{r}}\right)\mathrm{d}Z\left(\mathit{\nu }\right),\end{array}$$

which also uses the convolution theorem and where the sinc function represents the Fourier transform of the moving average window in Eq. (13). Then the auto-covariance function of F_X,n is obtained from Eqs. (44), (52), and (54) as

$$\begin{array}{}\text{(55)}& \begin{array}{rl}& \langle F_{X,n}^{\star }{F}_{X,n{}^{\prime }}{\rangle }_X\\ & =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathrm{2}\mathit{\pi }\mathit{\nu }(t_{n{}^{\prime }}-t_n)}\left|C\right(\mathit{\nu })-\text{sinc}(\mathit{\nu }{\mathit{\tau }}_{\mathrm{r}}\left){|}^{\mathrm{2}}{S}_X\right(\mathit{\nu })\mathrm{d}\mathit{\nu }.\end{array}\end{array}$$

Since, by analogy with Eq. (48), it can be shown that the cross power spectral density of ${\stackrel{\mathrm{̃}}{X}}_n$ and $V_n^{\left(j\right)}$ vanishes at all frequencies, the same holds for the error components F_X,n and W_X,n in Eq. (51), and, thus, the power spectral density of their sum equals the sum of their spectral densities. Finally, because ${\mathit{ϵ}}_n^{\left(j\right)}$ is also white in terms of j, we have, from Eqs. (50) and (53),

$$\begin{array}{}\text{(56)}& {\displaystyle }{\text{Var}}_{(X,\mathit{ϵ})}\left(W_{X,n}\right)& {\displaystyle }=\left[{\text{Var}}_{\left(X\right)}\right(X)-{\text{Var}}_{\left(X\right)}(U_n\left)\right]/N\text{(57)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{N}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[\mathrm{1}-|C\left(\mathit{\nu }\right){|}^{\mathrm{2}}\left]S_X\right(\mathit{\nu })\mathrm{d}\mathit{\nu },\end{array}$$

where the second step may be obtained directly from Eq. (53) by substitution from Eq. (49) into Eq. (47) with $n=n{}^{\prime }$.

From this we can now obtain a spectral representation of the squared reconstruction uncertainty components ${\mathcal{U}}_{\left(\mathrm{1}\right)}^{\mathrm{2}}$ and ${\mathcal{U}}_{\left(\mathrm{2}\right)}^{\mathrm{2}}$ as defined by Eq. (25) by writing the power spectral density of F_X,n, denoted by $S_{{\mathcal{U}}_{\left(\mathrm{1}\right)}}\left(\mathit{\nu }\right)$, and the power spectral density of W_X,n, denoted by $S_{{\mathcal{U}}_{\left(\mathrm{2}\right)}}\left(\mathit{\nu }\right)$. Specifically, from Eq. (55) we obtain (also taking into account spectral aliasing and leakage; see Priestley, 1981, for example)

$$\begin{array}{}\text{(58)}& \begin{array}{rl}{S}_{{\mathcal{U}}_{\left(\mathrm{1}\right)}}\left(\mathit{\nu }\right)& =\text{III}\left(\mathit{\nu };\mathrm{\Delta }{t}^{-\mathrm{1}}\right)\ast \mathit{\left\{}T{\text{sinc}}^{\mathrm{2}}\right(\mathit{\nu }T)\\ & \ast \left[{\left|C\left(\mathit{\nu }\right)-\text{sinc}\left(\mathit{\nu }{\mathit{\tau }}_{\mathrm{r}}\right)\right|}^{\mathrm{2}}{S}_X\right(\mathit{\nu }\left)\right]\mathit{\}},\end{array}\end{array}$$

with $-{\mathit{\nu }}_{*}<\mathit{\nu }\le {\mathit{\nu }}_{*}$. Here, ${\mathit{\nu }}_{*}=(\mathrm{2}\mathrm{\Delta }t{)}^{-\mathrm{1}}$ denotes the Nyquist frequency, Δt the sampling interval between the discrete sampling times t_n=nΔt, and T the length of the proxy record (being a multiple of Δt). Likewise, by confining the variance of the white noise process W_X,n, given by Eq. (57), to the same frequency interval, we obtain the constant spectral density

$$\begin{array}{}\text{(59)}& {\displaystyle }{S}_{{\mathcal{U}}_{\left(\mathrm{2}\right)}}\left(\mathit{\nu }\right)={\displaystyle \frac{\mathrm{\Delta }t}{N}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[\mathrm{1}-|C\left({\mathit{\nu }}^{\prime }\right){|}^{\mathrm{2}}\left]S_X\right({\mathit{\nu }}^{\prime })\mathrm{d}{\mathit{\nu }}^{\prime },\end{array}$$

with $-{\mathit{\nu }}_{*}<\mathit{\nu }\le {\mathit{\nu }}_{*}$. To understand the structure of $S_{{\mathcal{U}}_{\left(\mathrm{1}\right)}}\left(\mathit{\nu }\right)$, note that the term $\left|C\right(\mathit{\nu })-\text{sinc}(\mathit{\nu }{\mathit{\tau }}_{\mathrm{r}}){|}^{\mathrm{2}}$ in Eq. (58), referred to as the squared modulus of the error transfer function, acts as a linear filter on the stochastic component X(t) of the supposed true climate signal. Its structure is illustrated in Fig. 2b under the additional assumption that ${\mathit{\tau }}_{\mathrm{r}}^{-\mathrm{1}}\ll {\mathit{\nu }}_{\mathrm{c}}$. It turns out that it represents a multi-bandpass filter with the low-frequency band being confined between $[max({\mathit{\tau }}_{\mathrm{b}},{\mathit{\tau }}_{\mathrm{s}}){]}^{-\mathrm{1}}$ and ${\mathit{\tau }}_{\mathrm{r}}^{-\mathrm{1}}$ (this corresponds to the frequency band of the transfer function discussed by Amrhein, 2020; see his Fig. 2), whereas each high-frequency band is confined to an interval bounded by $k{\mathit{\nu }}_{\mathrm{c}}\pm [max({\mathit{\tau }}_{\mathrm{b}},{\mathit{\tau }}_{\mathrm{s}}){]}^{-\mathrm{1}}$ with $k=\pm \mathrm{1},\pm \mathrm{2},\mathrm{\dots }$, according to Eq. (42). The consequences of this particular filter structure, in conjunction with the effects of spectral aliasing, are discussed in Sect. 4. Finally, according to the finite length of the proxy record, we need to subsample the abovementioned power spectral densities at the discrete frequencies ν_m=mΔν (with $m=\mathrm{0},\pm \mathrm{1},\pm \mathrm{2},\mathrm{\dots }$ and $\mathrm{\Delta }\mathit{\nu }=\mathrm{1}/T$), which yields

$$\begin{array}{}\text{(60)}& {\displaystyle }{S}_{{\mathcal{U}}_{\left(\mathrm{1}\right)},m}=S_{{\mathcal{U}}_{\left(\mathrm{1}\right)}}\left({\mathit{\nu }}_m\right)\end{array}$$

and

$$\begin{array}{}\text{(61)}& {\displaystyle }{S}_{{\mathcal{U}}_{\left(\mathrm{2}\right)},m}=S_{{\mathcal{U}}_{\left(\mathrm{2}\right)}}\left({\mathit{\nu }}_m\right).\end{array}$$

Since F_X,n and W_X,n have zero cross power spectral densities, the power spectral density of E_X,n is then given by

$$\begin{array}{}\text{(62)}& {\displaystyle }{S}_{{\mathcal{U}}_{(\mathrm{1},\mathrm{2})},m}=S_{{\mathcal{U}}_{\left(\mathrm{1}\right)},m}+S_{{\mathcal{U}}_{\left(\mathrm{2}\right)},m}.\end{array}$$

3.2 Deterministic signal: discrete orbital spectrum

The deterministic signal Y(t), defined by Eq. (1), can be expressed as

$$\begin{array}{}\text{(63)}& {\displaystyle }Y\left(t\right)=Y_{\mathrm{c}}\left(t\right)[\mathrm{1}+Y_{\mathrm{a}}(t\left)\right],\end{array}$$

with the seasonal cycle oscillation,

$$\begin{array}{}\text{(64)}& {\displaystyle }{Y}_{\mathrm{c}}\left(t\right)={\mathit{\alpha }}_{\mathrm{c}}\left[Y_{\mathrm{c}}^{-}\right(t)+Y_{\mathrm{c}}^{+}(t\left)\right],\end{array}$$

and the amplitude modulating orbital oscillation,

$$\begin{array}{}\text{(65)}& {\displaystyle }{Y}_{\mathrm{a}}\left(t\right)={\mathit{\alpha }}_{\mathrm{a}}\left[Y_{\mathrm{a}}^{-}\right(t)+Y_{\mathrm{a}}^{+}(t\left)\right],\end{array}$$

where

$$\begin{array}{}\text{(66)}& {\displaystyle }{Y}_{\mathrm{c}}^{\pm }\left(t\right)=e^{\pm i(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{c}t+{\mathit{\varphi }}_{\mathrm{c}})},Y_{\mathrm{a}}^{\pm }\left(t\right)=e^{\pm i(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}t+{\mathit{\varphi }}_{\mathrm{a}})},\end{array}$$

and

$$\begin{array}{}\text{(67)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{c}}={\mathit{\sigma }}_{\mathrm{c}}/\sqrt{\mathrm{2}},{\mathit{\alpha }}_{\mathrm{a}}={\mathit{\sigma }}_{\mathrm{a}}/\sqrt{\mathrm{2}}.\end{array}$$

Then we can rewrite the signal (Eq. 63) as a complex Fourier series,

$$\begin{array}{}\text{(68)}& {\displaystyle }Y\left(t\right)={\mathcal{Y}}^{-}\left(t\right)+{\mathcal{Y}}^{+}\left(t\right),\end{array}$$

with

$$\begin{array}{}\text{(69)}& \begin{array}{rl}& {\mathcal{Y}}^{\pm }\left(t\right)\\ & ={\mathit{\alpha }}_{\mathrm{c}}{\mathit{\alpha }}_{\mathrm{a}}{Y}_{\mathrm{c}}^{\pm }\left(t\right)Y_{\mathrm{a}}^{\mp }\left(t\right)+{\mathit{\alpha }}_c{Y}_{\mathrm{c}}^{\pm }\left(t\right)+{\mathit{\alpha }}_{\mathrm{c}}{\mathit{\alpha }}_a{Y}_{\mathrm{c}}^{\pm }\left(t\right)Y_{\mathrm{a}}^{\pm }\left(t\right)\end{array}\end{array}$$

such that ${\mathcal{Y}}^{+}\left(t\right)=\left[{\mathcal{Y}}^{-}\right(t){]}^{\star }$. Thus, the right-hand side of Eq. (68) is the sum of six Fourier modes, two of which occur at the frequencies ±ν_c (the carrier wave in amplitude modulation terminology) and four of which occur at the frequencies $\pm {\mathit{\nu }}_{\mathrm{c}}\pm {\mathit{\nu }}_{\mathrm{a}}$ (representing the sidebands of ±ν_c).

Again following the approach of Balakrishnan (1962) and by analogy with Sect. 3.1, we evaluate $\langle {Y_n^{\left(j\right)}}^{\star }{Y}_{n^{\prime }}^{\left(j\right)}{\rangle }_{\mathit{ϵ}}$, where $Y_n^{\left(j\right)}=Y(t_n+{\mathit{ϵ}}_n^{\left(j\right)})$, as defined by Eq. (6). However, this is not the auto-covariance function in this case because sampling the seasonal cycle oscillation Y_c(t) at an average interval of Δt may leave a nonzero bias as Δt is a multiple of 1 year ($={\mathit{\nu }}_{\mathrm{c}}^{-\mathrm{1}}$). Thus, in general, we may have $\langle Y_n^{\left(j\right)}{\rangle }_{\mathit{ϵ}}\ne \mathrm{0}$, and it turns out that

$$\begin{array}{}\text{(70)}& {\displaystyle }\langle {Y_n^{\left(j\right)}}^{\star }{Y}_{n^{\prime }}^{\left(j\right)}{\rangle }_{\mathit{ϵ}}=\langle Y_n^{\left(j\right)}{\rangle }_{\mathit{ϵ}}^{\star }\langle Y_{n^{\prime }}^{\left(j\right)}{\rangle }_{\mathit{ϵ}}+{\text{Cov}}_{\left(\mathit{ϵ}\right)}\left({Y_n^{\left(j\right)}}^{\star },Y_{n^{\prime }}^{\left(j\right)}\right),\end{array}$$

where ${\text{Cov}}_{\left(\mathit{ϵ}\right)}({Y_n^{\left(j\right)}}^{\star },Y_{n^{\prime }}^{\left(j\right)})$ is the nonstationary auto-covariance function of $Y_n^{\left(j\right)}$. We now decompose this signal as

$$\begin{array}{}\text{(71)}& {\displaystyle }{Y}_n^{\left(j\right)}=A_n+B_n^{\left(j\right)},\end{array}$$

with the components

$$\begin{array}{}\text{(72)}& {\displaystyle }{A}_n=\langle Y_n^{\left(j\right)}{\rangle }_{\mathit{ϵ}},\end{array}$$

$$\begin{array}{}\text{(73)}& {\displaystyle }{B}_n^{\left(j\right)}=Y_n^{\left(j\right)}-\langle Y_n^{\left(j\right)}{\rangle }_{\mathit{ϵ}},\end{array}$$

and by analogy with the components U_n and $V_n^{\left(j\right)}$, respectively, in Sect. 3.1 such that we can rewrite Eq. (70) as

$$\begin{array}{}\text{(74)}& {\displaystyle }\langle {Y_n^{\left(j\right)}}^{\star }{Y}_{n^{\prime }}^{\left(j\right)}{\rangle }_{\mathit{ϵ}}=A_n^{\star }{A}_{n{}^{\prime }}+{\text{Cov}}_{\left(\mathit{ϵ}\right)}\left({B_n^{\left(j\right)}}^{\star },B_{n^{\prime }}^{\left(j\right)}\right).\end{array}$$

This also implies that A_n and $B_n^{\left(j\right)}$ are uncorrelated.

The structure of A_n is obtained from Eq. (68) by replacing t in the exponential terms in Eq. (66) by $n\mathrm{\Delta }t+{\mathit{ϵ}}_n^{\left(j\right)}$ and then applying the expected value operator. Note that because Δt is a multiple of 1 year, the modes at ±ν_c become aliases of ν=0, the modes at $\pm {\mathit{\nu }}_{\mathrm{c}}-{\mathit{\nu }}_{\mathrm{a}}$ become aliases of $\mathit{\nu }=-{\mathit{\nu }}_{\mathrm{a}}$, and those at $\pm {\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}}$ become aliases of ν=ν_a. Then considering phase interference caused by the aliasing, using Eqs. (35), (38), and (40), noting that $e^{\pm i\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{c}n\mathrm{\Delta }t}=\mathrm{1}$, and exploiting the symmetry property $C(-k{\mathit{\nu }}_{\mathrm{c}}+\mathit{\nu })=C(k{\mathit{\nu }}_{\mathrm{c}}+\mathit{\nu })$, we obtain

$$\begin{array}{}\text{(75)}& {\displaystyle }{A}_n=\mathrm{2}{\mathit{\alpha }}_{\mathrm{c}}\cos \left({\mathit{\varphi }}_{\mathrm{c}}\right)\text{sinc}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)A_n{}^{\prime },\end{array}$$

with

$$\begin{array}{}\text{(76)}& \begin{array}{rl}& A_n{}^{\prime }=\mathrm{1}+{\mathit{\alpha }}_{\mathrm{a}}\\ & \left({\widehat{f}}_{\mathrm{bs}}^{\star }(-{\mathit{\nu }}_{\mathrm{a}})e^{-i(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}})}+{\widehat{f}}_{\mathrm{bs}}^{\star }\left({\mathit{\nu }}_{\mathrm{a}}\right)e^{i(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}})}\right).\end{array}\end{array}$$

If we explicitly express the argument and the modulus of ${\widehat{f}}_{\mathrm{b}}^{\star }\left({\mathit{\nu }}_{\mathrm{a}}\right)$ as

$$\begin{array}{}\text{(77)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{b}\mathrm{1}}=\arg \left[{\widehat{f}}_{\mathrm{b}}^{\star }\right({\mathit{\nu }}_{\mathrm{a}}\left)\right]=\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{b}}-\arctan \left(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{b}}\right)\end{array}$$

and

$$\begin{array}{}\text{(78)}& {\displaystyle }{M}_{\mathrm{b}\mathrm{1}}=\left|{\widehat{f}}_{\mathrm{b}}^{\star }\right({\mathit{\nu }}_{\mathrm{a}}\left)\right|={\left[\mathrm{1}+(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{b}}{)}^{\mathrm{2}}\right]}^{-\mathrm{1}/\mathrm{2}},\end{array}$$

respectively, then we can rewrite Eqs. (75) and (76) as

$$\begin{array}{}\text{(79)}& {\displaystyle }{A}_n={\mathit{\sigma }}_{\mathrm{c}}\sqrt{\mathrm{2}}\cos \left({\mathit{\varphi }}_{\mathrm{c}}\right)\text{sinc}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)A_n{}^{\prime },\end{array}$$

with

$$\begin{array}{}\text{(80)}& {\displaystyle }{A}_n{}^{\prime }=\mathrm{1}+{\mathit{\sigma }}_{\mathrm{a}}\sqrt{\mathrm{2}}{M}_{\mathrm{b}\mathrm{1}}\text{sinc}\left({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right)\cos (\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}}+{\mathit{\varphi }}_{\mathrm{b}\mathrm{1}}).\end{array}$$

It turns out that, as long as we take the seasonal phase ϕ_c as fixed, A_n represents a deterministic bias caused by uneven sampling of the seasonal cycle due to proxy seasonality and that this bias varies in time because the amplitude of the seasonal cycle is modulated by orbital variations. Since ${\mathit{\sigma }}_{\mathrm{a}}\sqrt{\mathrm{2}}<\mathrm{1}$, the term A_n^′ is always positive, and, thus, the sign of the bias A_n is determined only by the seasonal phase ϕ_c. Note that the phase component ϕ_b1 of the oscillation results from the asymmetry of the bioturbation PDF, f_b(ϵ), defined by Eq. (2), which creates a time lag caused by bioturbation. However, if ${\mathit{\tau }}_{\mathrm{b}}\ll {\mathit{\nu }}_{\mathrm{a}}^{-\mathrm{1}}$, we have ϕ_b1≈0 and the time lag vanishes.

To understand the structure of $B_n^{\left(j\right)}$, we consider its variance given by the auto-covariance function at lag zero, ${\text{Cov}}_{\left(\mathit{ϵ}\right)}({B_n^{\left(j\right)}}^{\star },B_{n^{\prime }}^{\left(j\right)}){|}_{n=n{}^{\prime }}$. From this we find, as is shown in Appendix A, that $B_n^{\left(j\right)}$ is a nonstationary zero-mean white noise process. In particular, the variance has a stationary component and two time-varying components which oscillate at the frequencies ν_a and 2ν_a (see Eq. A12). To illustrate this behavior, we consider the following three simplified cases.

First, if there is no amplitude modulation of the seasonal cycle (${\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}=\mathrm{0}$), the variance of $B_n^{\left(j\right)}$ is stationary and is given by

$$\begin{array}{}\text{(81)}& \begin{array}{rl}& {\text{Var}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right)={\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}\\ & \mathit{\{}\mathrm{1}-{\text{sinc}}^{\mathrm{2}}({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}})+\cos (\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\left)\right[\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)-{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)\left]\mathit{\right\}}.\end{array}\end{array}$$

The dependence of this variance on the width of the proxy abundance window τ_p and its seasonal timing ϕ_c is illustrated in Fig. 4. If τ_p=0, the white noise variance vanishes because each year the same value is sampled from the seasonal cycle. If τ_p=1 year, the white noise variance equals the seasonal cycle variance ${\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}$. For intermediate values of τ_p the white noise variance depends on the seasonal phase. Note that for phases around $\left|{\mathit{\varphi }}_{\mathrm{c}}\right|=\mathit{\pi }/\mathrm{2}$ the white noise variance can exceed the seasonal cycle variance by up to 22 %.

Figure 4 The white noise variance ${\text{Var}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right)$ that is sampled from the seasonal cycle (normalized by the seasonal cycle variance ${\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}$) for the simplified case without amplitude modulation by orbital variations. The variance is shown as a function of the width of the proxy abundance window τ_p and of the seasonal phase ϕ_c, which together characterize proxy seasonality. Values greater than 1 are indicated by gray shading.

Second, if the seasonal cycle is modulated by orbital variations ($\mathrm{0}<{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}<\mathrm{1}/\mathrm{2}$) and we set τ_p=0 (and, for simplicity, we choose (i) τ_b=0 so as to avoid any additional phase lags, ${\mathit{\varphi }}_{\mathrm{b}\mathrm{1}}={\mathit{\varphi }}_{\mathrm{b}\mathrm{2}}=\mathrm{0}$, and (ii) ϕ_c=0 for a maximum effect of proxy seasonality), then the variance of $B_n^{\left(j\right)}$ is given by

$$\begin{array}{}\text{(82)}& \begin{array}{rl}& {\text{Var}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right)\\ & =\mathrm{2}{\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}\mathit{\{}\mathrm{1}-{\text{sinc}}^{\mathrm{2}}({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}})+\cos (\mathrm{4}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+\mathrm{2}{\mathit{\varphi }}_{\mathrm{a}})\\ & \left[\text{sinc}\right(\mathrm{2}{\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}})-{\text{sinc}}^{\mathrm{2}}({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\left)\right]\mathit{\}}.\end{array}\end{array}$$

Note the analogy with Eq. (81) but with ν_cτ_p and ϕ_c being replaced by ν_aτ_s and ϕ_a, respectively. This is because the white noise variance is now sampled from the orbital oscillation. Also it now has a time-varying component with frequency 2ν_a.

Third, if we consider the same case but with τ_p=1 year, we obtain

$$\begin{array}{}\text{(83)}& \begin{array}{rl}& {\text{Var}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right)\\ & ={\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}\mathit{\left\{}\right[\mathrm{1}+\mathrm{2}{\mathit{\sigma }}_{\mathrm{a}}\sqrt{\mathrm{2}}\text{sinc}\left({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right)\cos (\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}})]\\ & +{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}[\mathrm{1}+\text{sinc}(\mathrm{2}{\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\left)\cos \right(\mathrm{4}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}}\left)\right]\mathit{\}}.\end{array}\end{array}$$

In this case, the white noise variance has two time-varying components with frequencies ν_a and 2ν_a because the amplitude modulation factor has the basic structure 1+cos (2πν_at), the square of which, as it appears in the variance, is $\mathrm{1}+\mathrm{2}\cos \left(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}t\right)+{\cos }^{\mathrm{2}}\left(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}t\right)$. Note that the seasonal phase ϕ_c has no effect in this case with τ_p covering the full seasonal cycle.

With these properties of the abovementioned components A_n and $B_n^{\left(j\right)}$, we can now rewrite the error component E_Y,n, defined by Eq. (15), also using the Y component of Eq. (10) as

$$\begin{array}{}\text{(84)}& {\displaystyle }{E}_{Y,n}=F_{Y,n}+W_{Y,n},\end{array}$$

with

$$\begin{array}{}\text{(85)}& {\displaystyle }{F}_{Y,n}=A_n\end{array}$$

because ${\stackrel{\mathrm{̃}}{Y}}_n=\mathrm{0}$, and

$$\begin{array}{}\text{(86)}& {\displaystyle }{W}_{Y,n}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{j=\mathrm{1}}^N{B}_n^{\left(j\right)}.\end{array}$$

Then we obtain the reconstruction bias, defined by Eq. (23), as

$$\begin{array}{}\text{(87)}& {\displaystyle }{\mathcal{B}}_n=\langle A_n{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}={\mathit{\sigma }}_{\mathrm{c}}\sqrt{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{1}}{A}_n{}^{\prime },\end{array}$$

with

$$\begin{array}{}\text{(88)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{1}}=\cos (\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}})\text{sinc}({\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}/\mathrm{2}\mathit{\pi })\text{sinc}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right),\end{array}$$

where the seasonal phase uncertainty ${\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}$ and the expected seasonal cycle phase $\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}$ are defined by Eq. (19) and A_n^′ by Eq. (80). With this, the third squared reconstruction uncertainty component, defined by Eq. (25), can be expressed as

$$\begin{array}{}\text{(89)}& {\displaystyle }{\mathcal{U}}_{\left(\mathrm{3}\right),n}^{\mathrm{2}}=\langle A_n^{\mathrm{2}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}-\langle A_n{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{2}}{A}_n{{}^{\prime }}^{\mathrm{2}},\end{array}$$

with

$$\begin{array}{}\text{(90)}& \begin{array}{rl}& {\mathit{\gamma }}_{\mathrm{2}}=\mathit{\{}\mathrm{1}-{\text{sinc}}^{\mathrm{2}}({\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}/\mathrm{2}\mathit{\pi })+\cos (\mathrm{2}\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}})\\ & \left[\text{sinc}\right(\mathrm{2}{\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}/\mathrm{2}\mathit{\pi })-{\text{sinc}}^{\mathrm{2}}({\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}/\mathrm{2}\mathit{\pi }\left)\right]\mathit{\left\}}{\text{sinc}}^{\mathrm{2}}\right({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}).\end{array}\end{array}$$

Finally, from Eq. (86) and because ${\mathit{ϵ}}_n^{\left(j\right)}$ is white in terms of j, we have

$$\begin{array}{}\text{(91)}& {\displaystyle }{\text{Var}}_{\left(\mathit{ϵ}\right)}\left(W_{Y,n}\right)={\text{Var}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right)/N,\end{array}$$

and, thus, the fourth squared reconstruction uncertainty component is obtained as

$$\begin{array}{}\text{(92)}& {\displaystyle }{\mathcal{U}}_{\left(\mathrm{4}\right),n}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N}}\langle {\text{Var}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right){\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}},\end{array}$$

with ${\text{Var}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right)$ given by Eq. (A12) and where applying the expected value operator $\langle \cdot {\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}$ amounts to replacing each instance of cos (2ϕ_c) as it appears multiple times in the components of Eq. (A12), according to

$$\begin{array}{}\text{(93)}& {\displaystyle }\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\right)\stackrel{\text{repl.}}{⟶}\cos (\mathrm{2}\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}})\text{sinc}(\mathrm{2}{\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}/\mathrm{2}\mathit{\pi }).\end{array}$$

Note that $\text{sinc}(\mathrm{2}{\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}/\mathrm{2}\mathit{\pi })=\mathrm{0}$ if ${\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}=\mathit{\pi }$ or 2π, in which case the expected value in Eq. (92) simplifies because all terms in the components of Eq. (A12) that are multiplied by cos (2ϕ_c) vanish.

To obtain spectral representations of ${\mathcal{B}}_n^{\mathrm{2}}$ and ${\mathcal{U}}_{\left(\mathrm{3}\right),n}^{\mathrm{2}}$, we consider first the power spectral density of the signal A_n^′, limited to a finite time interval of length T (centered at t=0) and interpreted as the length of the proxy record. Specifically, we can express this discretized power spectral density as

$$\begin{array}{}\text{(94)}& {\displaystyle }{S}_{\mathrm{0},m}=\left|\mathcal{F}\right[\mathrm{\Pi }(t_n;T/\mathrm{\Delta }t)A_n{}^{\prime }\left]\right({\mathit{\nu }}_m){|}^{\mathrm{2}}/T,\end{array}$$

given at the discrete frequencies ν_m=mΔν (with $m=\mathrm{0},\pm \mathrm{1},\pm \mathrm{2},\mathrm{\dots }$ and $\mathrm{\Delta }\mathit{\nu }=\mathrm{1}/T$). In the above,

$$\mathcal{F}\left[x_n\right]\left(\mathit{\nu }\right)=\sum _{n=-n_{\mathrm{h}}}^{n_{\mathrm{h}}}{e}^{-i\mathrm{2}\mathit{\pi }\mathit{\nu }n\mathrm{\Delta }t}{x}_n\mathrm{\Delta }t$$

denotes the discrete time Fourier transform of a sequence x_n (with $-{\mathit{\nu }}_{*}<\mathit{\nu }\le {\mathit{\nu }}_{*}$), the discrete rectangle function acts as a window to confine A_n^′ to the finite time interval, and T∕Δt is the number of sampling times t_n=nΔt with $n=\mathrm{0},\pm \mathrm{1},\pm \mathrm{2},\mathrm{\dots },\pm n_{\mathrm{h}}$, where $n_{\mathrm{h}}=(T/\mathrm{\Delta }t-\mathrm{1})/\mathrm{2}$, for odd numbers of sampling times.

The discrete time Fourier transform of the rectangle function is given by the Dirichlet kernel (see Priestley, 1981, p. 437) which can be expressed as T multiplied by the aliased (or periodic) sinc function, defined by

$$\begin{array}{}\text{(95)}& \begin{array}{rl}& \text{asinc}(\mathit{\nu };T,\mathrm{\Delta }t)\\ & =\left\{\begin{array}{ll}\mathrm{1}& \text{if}\mathit{\nu }=\mathrm{0}\\ \sin \left(\mathit{\pi }\mathit{\nu }T\right)/\left[\sin \right(\mathit{\pi }\mathit{\nu }\mathrm{\Delta }t)T/\mathrm{\Delta }t]& \text{if}\mathit{\nu }\ne \mathrm{0}\end{array}\right.,\end{array}\end{array}$$

within the interval $-\mathrm{2}{\mathit{\nu }}_{*}<\mathit{\nu }<\mathrm{2}{\mathit{\nu }}_{*}$. With this, and if we express the discrete time Fourier transform of A_n^′ as a series of Dirac delta functions, we obtain from Eq. (76)

$$\begin{array}{}\text{(96)}& \begin{array}{rl}& \mathcal{F}\left[\mathrm{\Pi }\right(t_n;T/\mathrm{\Delta }t\left)A_n{}^{\prime }\right]\left(\mathit{\nu }\right)=T\text{asinc}(\mathit{\nu };T,\mathrm{\Delta }t)\\ & \ast \mathit{\left\{}\mathit{\delta }\right(\mathit{\nu })+{\mathit{\alpha }}_{\mathrm{a}}[\mathit{\delta }(\mathit{\nu }+{\mathit{\nu }}_{\mathrm{a}}){\widehat{f}}_{\mathrm{bs}}^{\star }(-{\mathit{\nu }}_{\mathrm{a}})e^{-i{\mathit{\varphi }}_{\mathrm{a}}}\\ & +\mathit{\delta }(\mathit{\nu }-{\mathit{\nu }}_{\mathrm{a}}){\widehat{f}}_{\mathrm{bs}}^{\star }\left({\mathit{\nu }}_{\mathrm{a}}\right)e^{i{\mathit{\varphi }}_{\mathrm{a}}}\left]\mathit{\right\}}.\end{array}\end{array}$$

Then from Eq. (94), also considering phase interference between the asinc functions, using Eqs. (77) and (78), and noting that the central asinc function centered at ν=0 has its zeros at the discrete frequencies ν_m, we have

$$\begin{array}{}\text{(97)}& S_{\mathrm{0},m}=T(S_{\mathrm{c},m}+S_{\mathrm{ca},m}+S_{\mathrm{a},m}),\end{array}$$

with

$$\begin{array}{}\text{(98)}& \begin{array}{rl}& S_{\mathrm{c},m}={\mathit{\delta }}_m,\\ & S_{\mathrm{ca},m}={\mathit{\delta }}_m{\mathit{\sigma }}_{\mathrm{a}}\sqrt{\mathrm{2}}{M}_{\mathrm{b}\mathrm{1}}\text{sinc}\left({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right)\mathrm{2}\mathit{\xi }\left(\mathrm{0}\right)\cos ({\mathit{\varphi }}_{\mathrm{a}}+{\mathit{\varphi }}_{\mathrm{b}\mathrm{1}}),\end{array}\end{array}$$

and

$$\begin{array}{}\text{(99)}& \begin{array}{rl}& S_{\mathrm{a},m}=({\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}/\mathrm{2})M_{\mathrm{b}\mathrm{1}}^{\mathrm{2}}{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right)\mathit{\left\{}{\mathit{\xi }}_{+}^{\mathrm{2}}\right({\mathit{\nu }}_m)+{\mathit{\xi }}_{-}^{\mathrm{2}}({\mathit{\nu }}_m)\\ & +\mathrm{2}{\mathit{\xi }}_{+}\left({\mathit{\nu }}_m\right){\mathit{\xi }}_{-}\left({\mathit{\nu }}_m\right)\cos \left[\mathrm{2}\right({\mathit{\varphi }}_{\mathrm{a}}+{\mathit{\varphi }}_{\mathrm{b}\mathrm{1}}\left)\right]\mathit{\}},\end{array}\end{array}$$

and where

$$\begin{array}{}\text{(100)}& {\displaystyle }{\mathit{\xi }}_{\pm }\left(\mathit{\nu }\right)=\text{asinc}(\mathit{\nu }\pm {\mathit{\nu }}_{\mathrm{a}};T,\mathrm{\Delta }t),\end{array}$$

with δ_m denoting the single-argument Kronecker delta with ${\mathit{\delta }}_{m=\mathrm{0}}=\mathrm{1}$ and ${\mathit{\delta }}_{m\ne \mathrm{0}}=\mathrm{0}$.

With Eq. (97) we obtain, by analogy with Eqs. (87) and (89), the spectral representation of ${\mathcal{B}}_n^{\mathrm{2}}$ given by

$$\begin{array}{}\text{(101)}& S_{\mathcal{B},m}=\mathrm{2}{\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}{S}_{\mathrm{0},m},\end{array}$$

with $-{\mathit{\nu }}_{*}<{\mathit{\nu }}_m\le {\mathit{\nu }}_{*}$ and the spectral representation of ${\mathcal{U}}_{\left(\mathrm{3}\right),n}^{\mathrm{2}}$ given by

$$\begin{array}{}\text{(102)}& S_{{\mathcal{U}}_{\left(\mathrm{3}\right)},m}={\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{2}}{S}_{\mathrm{0},m},\end{array}$$

with $-{\mathit{\nu }}_{*}<{\mathit{\nu }}_m\le {\mathit{\nu }}_{*}$. Since W_Y,n, given by Eq. (86), is white noise, we can express the spectral representation of ${\mathcal{U}}_{\left(\mathrm{4}\right),n}^{\mathrm{2}}$ as

$$\begin{array}{}\text{(103)}& S_{{\mathcal{U}}_{\left(\mathrm{4}\right)},m}={\displaystyle \frac{\mathrm{\Delta }t}{N}}\langle {\stackrel{\mathrm{‾}}{\text{Var}}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right){\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}},\end{array}$$

with $-{\mathit{\nu }}_{*}<{\mathit{\nu }}_m\le {\mathit{\nu }}_{*}$, and where ${\stackrel{\mathrm{‾}}{\text{Var}}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right)$ is given by Eq. (A19). Note that, as in Eq. (92), applying the expected value operator $\langle \cdot {\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}$ amounts to applying the replacement (Eq. 93) to the components of Eq. (A19). From the above, the power spectral density of E_Y,n is then given by

$$\begin{array}{}\text{(104)}& S_{\mathcal{B},{\mathcal{U}}_{(\mathrm{3},\mathrm{4})},m}=S_{\mathcal{B},m}+S_{{\mathcal{U}}_{\left(\mathrm{3}\right)},m}+S_{{\mathcal{U}}_{\left(\mathrm{4}\right)},m}.\end{array}$$

Hence, we obtain the spectral representation of Eq. (22) that relates the rms reconstruction error ℰ_n, defined by Eq. (20), to the reconstruction bias ℬ_n and the reconstruction uncertainty 𝒰_n:

$$\begin{array}{}\text{(105)}& S_{\mathcal{E},m}=S_{\mathcal{B},m}+S_{\mathcal{U},m},\end{array}$$

where $S_{\mathcal{U},m}=S_{{\mathcal{U}}_{\left(\mathrm{1}\right)},m}+S_{{\mathcal{U}}_{\left(\mathrm{2}\right)},m}+S_{{\mathcal{U}}_{\left(\mathrm{3}\right)},m}+S_{{\mathcal{U}}_{\left(\mathrm{4}\right)},m}$, the first two components of which are given by Eqs. (60) and (61), respectively, at the end of Sect. 3.1. Thus, S_ℰ,m is the power spectral density (given at the discrete frequencies ν_m) of the reconstruction error E_n (given at the discrete times t_n).

4 Timescale-dependent reconstruction uncertainty

The reconstruction uncertainty components 𝒰₍₁₎, 𝒰₍₂₎, 𝒰_(3),n, and 𝒰_(4),n and the reconstruction bias ℬ_n, defined in Sect. 2.5, can now be quantified using the expressions derived in Sect. 3. Specifically, given the set of parameters of the reconstruction uncertainty model (see Appendix B), including the specifications of the deterministic component of the supposed true climate signal, and given the power spectral density of the stochastic signal component, we obtain the following.

• First is the uncertainty component 𝒰₍₁₎ that arises from the various smoothing processes which affect the stochastic signal component X(t) in the limit of infinitely many signal carriers retrieved from each slice of sediment (N→∞). From Eq. (60) or from Eq. (58), noting that for stochastic signals spectral aliasing and leakage neither generate nor destroy but only redistribute power spectral density, we have

  $$\begin{array}{}\text{(106)}& {\displaystyle }{\mathcal{U}}_{\left(\mathrm{1}\right)}^{\mathrm{2}}& {\displaystyle }=\sum _{m=-m_{\mathrm{h}}}^{m_{\mathrm{h}}}{S}_{{\mathcal{U}}_{\left(\mathrm{1}\right)},m}\mathrm{\Delta }\mathit{\nu }\text{(107)}& {\displaystyle }& {\displaystyle }=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\left|C\right(\mathit{\nu })-\text{sinc}(\mathit{\nu }{\mathit{\tau }}_{\mathrm{r}}\left){|}^{\mathrm{2}}{S}_X\right(\mathit{\nu })\mathrm{d}\mathit{\nu },\end{array}$$

  with $m_{\mathrm{h}}=(T/\mathrm{\Delta }t-\mathrm{1})/\mathrm{2}$ for odd numbers of sampling times, and C(ν) is given by Eq. (38). This uncertainty component depends on the timescale parameters τ_b, τ_s, τ_r, and τ_p and on S_X(ν).

• Second is the white noise uncertainty component 𝒰₍₂₎ that arises from sampling only a finite number N of signal carriers from each slice of sediment. From Eq. (61), or likewise from Eq. (59), we obtain

  $$\begin{array}{}\text{(108)}& {\displaystyle }{\mathcal{U}}_{\left(\mathrm{2}\right)}^{\mathrm{2}}& {\displaystyle }=\sum _{m=-m_{\mathrm{h}}}^{m_{\mathrm{h}}}{S}_{{\mathcal{U}}_{\left(\mathrm{2}\right)},m}\mathrm{\Delta }\mathit{\nu }\text{(109)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{N}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[\mathrm{1}-|C\left(\mathit{\nu }\right){|}^{\mathrm{2}}\left]S_X\right(\mathit{\nu })\mathrm{d}\mathit{\nu }.\end{array}$$

  This uncertainty component depends on the timescale parameters τ_b, τ_s, and τ_p and on S_X(ν) and N.

• Third are the reconstruction bias ℬ_n, its uncertainty 𝒰_(3),n caused by the imperfectly known seasonal timing of the proxy abundance, and the white noise uncertainty component 𝒰_(4),n arising from sampling only a finite number of signal carriers, which are all readily given in the time domain by Eqs. (87), (89), and (92), respectively. These components depend on the timescale parameters τ_b, τ_s, and τ_p, on the seasonal phase parameters $\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}$ and ${\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}$, and on the specifications of the deterministic signal components σ_c, σ_a, ν_c, ν_a, and ϕ_a. The white noise component 𝒰_(4),n also depends on N. The time averages of the squares of these components over the length of the proxy record T can be obtained directly from their spectral representations, given by Eqs. (101) to (103) as

  $$\begin{array}{}\text{(110)}& \begin{array}{rl}& \stackrel{\mathrm{‾}}{{\mathcal{B}}_n^{\mathrm{2}}}=\sum _{m=-m_{\mathrm{h}}}^{m_{\mathrm{h}}}{S}_{\mathcal{B},m}\mathrm{\Delta }\mathit{\nu },\\ & \stackrel{\mathrm{‾}}{{\mathcal{U}}_{\left(\mathrm{3}\right),n}^{\mathrm{2}}}=\sum _{m=-m_{\mathrm{h}}}^{m_{\mathrm{h}}}{S}_{{\mathcal{U}}_{\left(\mathrm{3}\right)},m}\mathrm{\Delta }\mathit{\nu },\\ & \stackrel{\mathrm{‾}}{{\mathcal{U}}_{\left(\mathrm{4}\right),n}^{\mathrm{2}}}=\sum _{m=-m_{\mathrm{h}}}^{m_{\mathrm{h}}}{S}_{{\mathcal{U}}_{\left(\mathrm{4}\right)},m}\mathrm{\Delta }\mathit{\nu },\end{array}\end{array}$$

  which then also depend on T.

Since all of the abovementioned uncertainty components and the bias depend on a number of timescale parameters, the rms reconstruction error ℰ_n, in this sense, already represents a timescale-dependent uncertainty measure. However, we may extend the concept of uncertainty timescale dependence as follows.

In practice, during the process of data analysis, climate reconstructions are often smoothed by some linear filter either because one is explicitly interested in time averages of the reconstructed climate variable or because one may hope to reduce the reconstruction uncertainty by averaging out short-timescale noise. However, the extent to which the uncertainty actually shrinks depends on the auto-correlation structure of the reconstruction error, which, by the Wiener–Khintchine theorem, is directly related to the power spectral density of the error. Thus, from the expressions of the error power spectral densities, derived in Sect. 3, we can directly quantify the uncertainty reduction that is achieved by applying a linear filter, as is shown in the following for the uncertainty components 𝒰₍₁₎ and 𝒰₍₂₎ as an example.

If the reconstruction error time series is smoothed by, for simplicity, a discrete moving average filter of width τ₀ (being a multiple of Δt), then the squared uncertainty obtained after smoothing is given by

$$\begin{array}{}\text{(111)}& {\mathcal{U}}_{(\mathrm{1},\mathrm{2})}^{\mathrm{2}}\left({\mathit{\tau }}_{\mathrm{0}}\right)=\sum _{m=-m_{\mathrm{h}}}^{m_{\mathrm{h}}}{\text{asinc}}^{\mathrm{2}}({\mathit{\nu }}_m;{\mathit{\tau }}_{\mathrm{0}},\mathrm{\Delta }t)S_{{\mathcal{U}}_{(\mathrm{1},\mathrm{2})},m}\mathrm{\Delta }\mathit{\nu },\end{array}$$

where the squared asinc function represents the squared modulus of the discrete time Fourier transform of the filter window acting as a spectral transfer function, and the asinc function is defined by Eq. (95). Note that if τ₀=Δt (i.e., no smoothing), then this transfer function is constantly 1 across all frequencies, and that if τ₀=T, then it is equal to 1 at frequency zero, and it is 0 at all other frequencies. Thus, in the latter case the uncertainty of the time average over the entire proxy record is obtained. Figure 5 illustrates the above for some choices of parameters designed to exemplify as many aspects as possible of the uncertainty estimation procedure in a single example at the expense of using somewhat unrealistic parameter values. More realistic application examples of the method follow in Part 2 of this study (Dolman et al., 2020). Specifically, we set τ_b=10 years, ${\mathit{\tau }}_{\mathrm{s}}={\mathit{\tau }}_{\mathrm{r}}=\mathrm{6}$ years, ${\mathit{\tau }}_{\mathrm{p}}=\mathrm{1}/\mathrm{3}$ year, Δt=6 years, $T=\mathrm{23}\mathrm{\Delta }t=\mathrm{138}$ years, ${\mathit{\tau }}_{\mathrm{0}}=\mathrm{3}\mathrm{\Delta }t=\mathrm{18}$ years, and N=100. For the power spectral density of X(t) we assume a Lorentzian shaped AR(1) red noise spectrum given by $S_X\left(\mathit{\nu }\right)=\mathrm{2}\mathit{\alpha }/\left[\right(\mathrm{2}\mathit{\pi }\mathit{\nu }{)}^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}]$ with the characteristic timescale ${\mathit{\alpha }}^{-\mathrm{1}}=(\mathrm{1}/\mathrm{10})$ years such that the process X is only weakly red.

Figure 5 Illustration of the method for estimating timescale-dependent reconstruction uncertainties in the spectral domain, shown for the uncertainty components 𝒰₍₁₎ and 𝒰₍₂₎ which are based on the stochastic component X(t) of the supposed true climate signal. (a) Power spectral density S_X(ν) (gray line) of this signal component, defined on a continuous and infinite frequency axis (normalized by the seasonal cycle frequency ${\mathit{\nu }}_{\mathrm{c}}=(\mathrm{1}\mathrm{yr}{)}^{-\mathrm{1}}$), together with the product $\left|C\right(\mathit{\nu }\left){|}^{\mathrm{2}}{S}_X\right(\mathit{\nu })$ (red line) and the product $\left|C\right(\mathit{\nu })-\text{sinc}(\mathit{\nu }{\mathit{\tau }}_{\mathrm{r}}\left){|}^{\mathrm{2}}{S}_X\right(\mathit{\nu })$ (blue line). The green area equals the integral ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\left[\mathrm{1}-{\left|C\left(\mathit{\nu }\right)\right|}^{\mathrm{2}}\right]S_X\left(\mathit{\nu }\right)\mathrm{d}\mathit{\nu }$, as it appears in Eq. (59), which measures the variance of the white noise error component caused by sampling only a finite number of signal carriers. (b) The same white noise variance (green area) but divided by N (the number of signal carriers) after spectral aliasing and leakage have been applied to obtain the power spectral density $S_{{\mathcal{U}}_{\left(\mathrm{2}\right)},m}$ (green dots), defined by Eq. (61), on a finite and discrete frequency axis now normalized by the Nyquist frequency ${\mathit{\nu }}_{*}=(\mathrm{2}\mathrm{\Delta }t{)}^{-\mathrm{1}}$, together with $S_{{\mathcal{U}}_{\left(\mathrm{1}\right)},m}$ (blue dots), defined by Eq. (60), and $S_{{\mathcal{U}}_{(\mathrm{1},\mathrm{2})},m}$ (magenta dots), defined by Eq. (62). Cyan dots indicate the same as blue dots but neglect the effect of spectral leakage for comparison. (c) The product ${\text{asinc}}^{\mathrm{2}}({\mathit{\nu }}_m;{\mathit{\tau }}_{\mathrm{0}},\mathrm{\Delta }t)S_{{\mathcal{U}}_{(\mathrm{1},\mathrm{2})},m}$ (magenta dots), the integral of which (magenta area) equals the squared reconstruction uncertainty ${\mathcal{U}}_{(\mathrm{1},\mathrm{2})}^{\mathrm{2}}\left({\mathit{\tau }}_{\mathrm{0}}\right)$, defined by Eq. (111).

This power spectral density is shown in Fig. 5a by the gray line. According to the reconstruction uncertainty model defined in Sect. 2, S_X(ν) is decomposed into two components: (i) |C(ν)|²S_X(ν), shown by the red line, the integral of which equals the variance of U_n defined by Eq. (44) and where |C(ν)|² (shown in Fig. 2a) acts as a spectral transfer function on S_X(ν); (ii) [1−|C(ν)|²]S_X(ν), the integral of which, indicated by the green area, equals the variance of the white noise component $V_n^{\left(j\right)}$ defined by Eq. (45). If S_X(ν) is multiplied by the squared modulus of the error transfer function (shown in Fig. 2b), the component ${\left|C\left(\mathit{\nu }\right)-\text{sinc}\left(\mathit{\nu }{\mathit{\tau }}_{\mathrm{r}}\right)\right|}^{\mathrm{2}}{S}_X\left(\mathit{\nu }\right)$, shown by the blue line, is obtained, the integral of which equals the variance of F_X,n defined by Eq. (52).

This component (blue dots), as well as the white noise component (green dots), is shown again in Fig. 5b but after spectral aliasing and leakage have been applied according to the sampling and measurement procedure described in Sect. 2.3. These components represent the discretized power spectral densities $S_{{\mathcal{U}}_{\left(\mathrm{1}\right)},m}$ and $S_{{\mathcal{U}}_{\left(\mathrm{2}\right)},m}$. Note that the broad peaks at nonzero frequencies in Fig. 5b are direct images of the low-frequency peaks in Fig. 5a (blue line), whereas the bump centered at ν=0 represents the summed aliases of the high-frequency peaks in Fig. 5a (blue line) at ±kν_c. Without proxy seasonality (τ_p=1 year), those peaks do not exist, and, thus, the $S_{{\mathcal{U}}_{\left(\mathrm{1}\right)},m}$ power spectrum in Fig. 5b falls off to near zero at ν=0. Only spectral leakage may then lead to nonzero $S_{{\mathcal{U}}_{\left(\mathrm{1}\right)},m}$ at ν=0, although in the example shown here the effect of the leakage is small (cyan dots). However, in cases with small T (implying large Δν), spectral leakage can provide a relevant contribution of power at ν=0 as the power from the neighboring broad spectral peaks is then effectively redistributed to the center of the frequency domain. This contribution is particularly important if the uncertainty of the time average over the entire proxy record is computed by setting τ₀=T in Eq. (111) as in this case it is the only contribution to ${\mathcal{U}}_{\left(\mathrm{1}\right)}^{\mathrm{2}}\left({\mathit{\tau }}_{\mathrm{0}}\right)$.

Finally, the sum of the abovementioned components, given by $S_{{\mathcal{U}}_{(\mathrm{1},\mathrm{2})},m}$, is shown by the magenta dots in Fig. 5c. By multiplying this summed power spectral density by the spectral transfer function of the discrete moving average window mentioned above and then integrating it, we obtain the squared reconstruction uncertainty after smoothing, ${\mathcal{U}}_{(\mathrm{1},\mathrm{2})}^{\mathrm{2}}\left({\mathit{\tau }}_{\mathrm{0}}\right)$, defined by Eq. (111) and shown by the magenta area in Fig. 5c.

Likewise, one can define other timescale-dependent uncertainty metrics. For example, one might be interested in the uncertainty of the difference between the time averages over two periods of length T₁ and T₂, which are separated in time by the interval $\mathit{\delta }t=(n{}^{\prime }-n)\mathrm{\Delta }t$, measured from center to center. This can be expressed by the variance

$$\begin{array}{}\text{(112)}& \begin{array}{rl}& \text{Var}\left(\right[T_{\mathrm{1}}^{-\mathrm{1}}\mathrm{\Pi }(t_n;T_{\mathrm{1}}/\mathrm{\Delta }t)\ast E_{X,n}]\\ & -[T_{\mathrm{2}}^{-\mathrm{1}}\mathrm{\Pi }\left(t_n;T_{\mathrm{2}}/\mathrm{\Delta }t\right)\ast E_{X,n+n{}^{\prime }}])\\ & =\text{Var}\left(T_{\mathrm{1}}^{-\mathrm{1}}\mathrm{\Pi }\left(t_n;T_{\mathrm{1}}/\mathrm{\Delta }t\right)\ast E_{X,n}\right)\\ & +\text{Var}(T_{\mathrm{2}}^{-\mathrm{1}}\mathrm{\Pi }\left(t_n;T_{\mathrm{2}}/\mathrm{\Delta }t\right)\ast E_{X,n+n{}^{\prime }})-\mathrm{2}\text{Cov}\\ & \left(T_{\mathrm{1}}^{-\mathrm{1}}\mathrm{\Pi }\left(t_n;T_{\mathrm{1}}/\mathrm{\Delta }t\right)\ast E_{X,n},T_{\mathrm{2}}^{-\mathrm{1}}\mathrm{\Pi }\left(t_n;T_{\mathrm{2}}/\mathrm{\Delta }t\right)\ast E_{X,n+n{}^{\prime }}\right)\end{array}\end{array}$$

if it were to be computed for the uncertainty components based on X(t) and where T₁, T₂ and δt are multiples of Δt. Then, using the Wiener–Khintchine theorem, we obtain the difference uncertainty metric

$$\begin{array}{}\text{(113)}& \begin{array}{rl}& \mathit{\delta }{\mathcal{U}}_{(\mathrm{1},\mathrm{2})}^{\mathrm{2}}(T_{\mathrm{1}},T_{\mathrm{2}},\mathit{\delta }t)\\ & =\sum _{m=-m_{\mathrm{h}}}^{m_{\mathrm{h}}}\mathit{\left\{}\right[{\text{asinc}}^{\mathrm{2}}({\mathit{\nu }}_m;T_{\mathrm{1}},\mathrm{\Delta }t)+{\text{asinc}}^{\mathrm{2}}({\mathit{\nu }}_m;T_{\mathrm{2}},\mathrm{\Delta }t)]\\ & S_{{\mathcal{U}}_{(\mathrm{1},\mathrm{2})},m}\mathrm{\Delta }\mathit{\nu }\mathit{\}}-\mathrm{2}{\mathcal{F}}^{-\mathrm{1}}\\ & \left[\text{asinc}\right({\mathit{\nu }}_m;T_{\mathrm{1}},\mathrm{\Delta }t\left)\text{asinc}\right({\mathit{\nu }}_m;T_{\mathrm{2}},\mathrm{\Delta }t\left)S_{{\mathcal{U}}_{(\mathrm{1},\mathrm{2})},m}\right]\left(\mathit{\delta }t\right),\end{array}\end{array}$$

where ${\mathcal{F}}^{-\mathrm{1}}\left[x_m\right]\left(k\mathrm{\Delta }t\right)={\sum }_{m=-m_{\mathrm{h}}}^{m_{\mathrm{h}}}{e}^{i\mathrm{2}\mathit{\pi }m\mathrm{\Delta }\mathit{\nu }k\mathrm{\Delta }t}{x}_m\mathrm{\Delta }\mathit{\nu }$ denotes the inverse discrete Fourier transform of a sequence x_m. Note that in this form the abovementioned difference uncertainty metric is valid only for stationary uncertainty components. If the seasonal cycle amplitude is constant over time (i.e., no orbital variations), then all uncertainty components are stationary. If orbital variations are taken into account, however, only 𝒰₍₁₎ and 𝒰₍₂₎ are stationary, as shown in Sect. 3. For the components 𝒰_(3),n, 𝒰_(4),n, and ℬ_n, the difference metric, in this case, had to be computed directly from their time domain expressions; see Eqs. (87), (89), and (92), respectively.

To conclude this section, we briefly present an example of the time series and power spectra of the reconstruction bias ℬ_n and of the reconstruction uncertainty components 𝒰_(3),n and 𝒰_(4),n. Specifically, we set ${\mathit{\tau }}_{\mathrm{b}}={\mathit{\tau }}_{\mathrm{s}}=\mathrm{100}$ years, ${\mathit{\tau }}_{\mathrm{p}}=\mathrm{1}/\mathrm{3}$ year, Δt=100 years, $T=\mathrm{101}\mathrm{\Delta }t=\mathrm{10}\mathrm{100}$ years, and N=5. Furthermore, the deterministic signal is specified in this example by the parameters ${\mathit{\sigma }}_{\mathrm{c}}=\sqrt{\mathrm{1}/\mathrm{2}}$, $\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}=\mathit{\pi }/\mathrm{4}$, ${\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}=\mathit{\pi }/\mathrm{2}$, ${\mathit{\sigma }}_{\mathrm{a}}=\sqrt{\mathrm{1}/\mathrm{8}}$, ${\mathit{\nu }}_{\mathrm{a}}=(\mathrm{23}\mathrm{kyrs}{)}^{-\mathrm{1}}$, and ${\mathit{\varphi }}_{\mathrm{a}}=\mathit{\pi }/\mathrm{2}$, which implies that the amplitude of the seasonal cycle decreases during the 10 100 years, as is the case during the Holocene. According to the proxy seasonality parameter values chosen here, the bias ℬ_n is positive and exhibits a negative trend, as shown in Fig. 6a (cyan line). Likewise the uncertainty components 𝒰_(3),n and 𝒰_(4),n also decrease over time (blue and green lines). Since the orbital modulation frequency ν_a is located on the discrete frequency axis between ν₀ and ν₁, its power spectral density is distributed by spectral leakage across all frequencies. This yields a highly red power spectrum of ℬ_n and 𝒰_(3),n, shown in Fig. 6b. Thus, at high frequencies the summed power spectral densities $S_{{\mathcal{U}}_{(\mathrm{3},\mathrm{4})},m}$ (dashed black line) and $S_{\mathcal{B},{\mathcal{U}}_{(\mathrm{3},\mathrm{4})},m}$ (dashed gray line) are dominated by the white noise component, whereas at low frequencies they are dominated by the effect of orbital variations. Hence, if we were to compute the timescale-dependent uncertainty metric (Eq. 111) but for 𝒰_(3),n and 𝒰_(4),n, denoted by ${\mathcal{U}}_{(\mathrm{3},\mathrm{4}),n}^{\mathrm{2}}\left({\mathit{\tau }}_{\mathrm{0}}\right)$, the uncertainty would shrink only slowly for increasing values of τ₀ because the orbital variations are associated with a highly correlated error in time at long timescales.

Figure 6 Example of the reconstruction uncertainty based on the deterministic component Y(t) of the supposed true climate signal. (a) Time series of the reconstruction bias ℬ_n (cyan line), of the uncertainty components 𝒰_(3),n (blue line) and 𝒰_(4),n (green line), of $[{\mathcal{U}}_{(\mathrm{3},\mathrm{4}),n}^{\mathrm{2}}{]}^{\mathrm{1}/\mathrm{2}}$ (dashed black line), and $[{\mathcal{B}}_n^{\mathrm{2}}+{\mathcal{U}}_{(\mathrm{3},\mathrm{4}),n}^{\mathrm{2}}{]}^{\mathrm{1}/\mathrm{2}}$ (dashed gray line). (b) Power spectral densities of the corresponding error components: S_ℬ,m (cyan line), $S_{{\mathcal{U}}_{\left(\mathrm{3}\right)},m}$ (blue line), $S_{{\mathcal{U}}_{\left(\mathrm{4}\right)},m}$ (green line), $S_{{\mathcal{U}}_{(\mathrm{3},\mathrm{4})},m}$ (dashed black line), and $S_{\mathcal{B},{\mathcal{U}}_{(\mathrm{3},\mathrm{4})},m}$ (gray line).

5 Discussion

To allow for an analytic treatment of the problem, the method for estimating timescale-dependent reconstruction uncertainties, presented in Sects. 2 to 4, is necessarily based on a number of simplifying assumptions.

• We assume a fixed proxy seasonality in the sense of applying every year the same seasonal timing of a prescribed proxy abundance period, characterized by the parameters τ_p and ϕ_c. For this reason we have to separate the supposed true climate signal into a stochastic component X(t) and a deterministic component Y(t) that represents the seasonal cycle because proxy seasonality then implies an in-phase subsampling from Y(t) which, in turn, affects the amount of variance aliased from the seasonal cycle 𝒰_(4),n and which may also lead to a reconstruction bias ℬ_n and associated uncertainty 𝒰_(3),n. This scenario represents the extreme case where a seasonal abundance period is completely imposed on the proxy by an external process (see, for example, Leduc et al., 2010), such as, for example, seasonally determined nutrient supply possibly controlled by the seasonality of solar irradiance or oceanic upwelling. By contrast, in the opposite extreme case where no seasonality is imposed at all, we do not need to separate the climate signal into X(t) and Y(t). In this case the total climate signal is fully recorded by the proxy, but its total variance is reduced by some factor because of habitat tracking if the habitat PDF of the proxy is narrower than the PDF of the climate signal. According to the idea of Mix (1987), this factor can be obtained by multiplying the two PDFs, and it may possibly also be expressed as a frequency-dependent spectral transfer function. This scenario corresponds to setting τ_p=1 year in our method and subsequently multiplying the obtained error power spectrum by the aforementioned transfer function. Hence, if we introduce some parameter, $\mathrm{0}\le s\le \mathrm{1}$, that measures the extent to which seasonality is imposed for a specific proxy record (with s=0 indicating no imposition of seasonality), then we may express the actual uncertainty as a linear combination of the uncertainties obtained from the abovementioned two extreme scenarios, weighted by s and by 1−s, respectively. Note, however, that the effects of seasonality can be rather complex (see, for example, Jonkers and Kučera, 2017), depending on the type of proxy used, and, thus, the optimal strategy for modeling the associated uncertainties depends on the specific application.

• We neglect calibration errors representing uncertainties regarding the climate–proxy relationship. Assuming this relationship is linear and is obtained by linear regression, errors of this type may have two effects. Uncertainties in the intercept parameter will introduce a reconstruction uncertainty that is constant in time like the bias uncertainty 𝒰_(3),n (in the case without orbital variations). Uncertainties in the slope parameter, on the other hand, will introduce a frequency-independent uncertainty in the error variance. The mean of the possible error variances, however, might be close to the variance obtained from our method unless the PDF of the obtained error variances is strongly skewed. If the climate–proxy relationship is non-linear, or if there are uncertainties regarding the linearity itself, modeling of the implied uncertainties might be more complex, although it should still be possible to decompose those errors into a bias and a variance component.

• We assume a constant sediment accumulation rate and a constant bioturbation depth, and we also assume regular sampling from the sediment core and neglect dating uncertainties, although relaxing these assumptions may generate additional uncertainties of noticeable magnitude. For example, the relevance of dating uncertainties is demonstrated by Goswami et al. (2014) and Boers et al. (2017). If these sources of uncertainty are treated in a stochastic sense, they could, in principle, be included in our approach by allowing for correlated sampling jitter ϵ, the mathematical basis of which is given by Balakrishnan (1962) (see also Moore and Thomson, 1991). More generally, these uncertainties could be modeled by allowing for a variable time–depth relationship and perhaps by also allowing for the nonstationarity of the uncertainty components 𝒰₁ and 𝒰₂ to represent variations of the smoothing timescales τ_b and τ_s.

From the above, it turns out that, in its current form, the method is neither complete in terms of processes affecting the reconstruction uncertainty nor does it cover all possible reconstruction scenarios in terms of proxy type and application context. However, our formulation of the method outlines a conceptually and mathematically well-founded approach of how timescale-dependent reconstruction uncertainties could, and probably should, be estimated – in particular, when systematic and exact quantification is required. This latter point is highly relevant, for example, in the context of comparisons between circulation models and paleo-observations (e.g., Lohmann et al., 2013; Laepple and Huybers, 2014; Matsikaris et al., 2016) or likewise for any reanalysis efforts (e.g., Hakim et al., 2016) if data obtained from proxy records are involved. Thus, the fact that some of the neglected sources of uncertainty might be large compared to what is gained by our exact mathematical treatment does not qualify our approach as overly precise. The approach rather demonstrates the directions for future efforts in quantitative uncertainty estimation. As discussed above, our current formulation of the method may indeed be extended beyond the simplifications made. But as mathematical complexity increases in such case, extended formulations should be tailored to specific applications. In this sense, our formulation provides a minimal basis for the development of future uncertainty estimation methods.

Furthermore, the timescale-dependent uncertainties obtained from our method depend explicitly on assumptions regarding the structure of the supposed true climate signal X(t)+Y(t), although this climate signal is the unknown quantity to be reconstructed from the proxy record. However, it is an inevitable fact that the timescale-dependent reconstruction uncertainties do actually depend on this structure, a fact that is made obvious by our method and likewise by Amrhein (2020). One possible approach towards solving this problem would be an iterative procedure: (i) assume a specific structure for the supposed true climate signal; (ii) apply our method to obtain reconstruction uncertainties for a given proxy record; (iii) check whether the reconstructed signal is consistent under the obtained uncertainties with the assumed structure given its spectral or auto-correlation properties; and (iv) if this is not the case, update the assumptions and repeat these steps.

Finally, although our method provides an advancement in the quantification of reconstruction uncertainties, it also introduces a number of model parameters which are associated with their own uncertainty. However, if we are to improve quantitative uncertainty estimates, our reconstruction uncertainty model helps to identify those parameters which are most important and, therefore, need to be determined at higher precision. For example, how much seasonality is imposed on a certain proxy at a given geographical location within a specific local ecological system? On the other hand, it is possible to investigate how parameter uncertainties translate into reconstruction uncertainties, as was shown for the seasonal phase parameter ϕ_c. Nonetheless, the eventual benefit of uncertainty estimation methods like the one presented in this study and of extensions based thereon has still to be worked out in the future by systematically applying such methods to real data.

6 Conclusions

The present study introduces a method, the so-called Proxy Spectral Error Model (PSEM; see also Part 2 of this study by Dolman et al., 2020), for estimating timescale-dependent uncertainties of paleoclimate reconstructions obtained from single sediment proxy records. The method is based on an uncertainty model that takes into account proxy seasonality (together with orbital variations of the seasonal cycle amplitude), bioturbation, archive sampling parameters, and the effects of measuring only a finite number of signal carriers. For this model, analytic expressions are derived for the power spectrum of the reconstruction error from which timescale-dependent reconstruction uncertainties can be obtained. This approach is motivated by the fact that the spectral structure of the error is equivalent to its auto-correlation structure which, in turn, determines how archive smoothing, sampling, and averaging timescales affect the uncertainties. Various timescale-dependent uncertainty metrics can be defined and then computed from the error power spectrum by multiplying the spectrum by specific transfer functions and then integrating. This corresponds, in the time domain, to additional postprocessing steps performed on the reconstructed time series. For example, it is possible to investigate the uncertainty reduction achieved by a low pass filter with a given cut-off timescale or to quantify the uncertainty of the difference between two time averages with given averaging timescales.

The method proves useful in different ways. First, it can serve to obtain quantitative uncertainty estimates for practical applications in paleoclimate science. This is demonstrated in Part 2 of this study (Dolman et al., 2020) in which a number of application examples are presented. Second, the derived analytic expressions can be used to acquire a better qualitative understanding of the structure of the uncertainties. In particular, we can conclude the following.

• The reconstruction uncertainties can be decomposed into two components. The first is a component whose variance is obtained by multiplying the power spectrum of the supposed true climate signal by a transfer function and then integrating. This so-called error transfer function has a structure corresponding to a bandpass filter with its cut-off timescales given by the longest applied archive smoothing timescale and by a suitably chosen reference smoothing timescale (by analogy with the transfer function discussed by Amrhein, 2020). Thus, multiplying the spectrum by the error transfer function corresponds to applying that bandpass filter to the supposed true climate signal. The second is a white noise component that scales inversely with the number of signal carriers retrieved from each slice of sediment (and being subject to the same single laboratory measurement). Thus, in the asymptotic limit of infinitely many signal carriers, this component vanishes. In the opposite limit, with only a single signal carrier being measured from each slice, the variance of this component equals the variance that is contained in the supposed true climate signal at timescales shorter than the longest applied archive smoothing timescale. This component corresponds to what is referred to, according to Dolman and Laepple (2018), as the noise created by the aliasing of variability from interannual and intra-annual timescales. Depending on geographical location and climatic conditions, this white noise uncertainty component may be dominated by ENSO variability or by the seasonal cycle, for example.

• In the presence of proxy seasonality such that the climate signal is recorded by the proxy only during a limited seasonal window each year, the abovementioned error transfer function has additional high-frequency peaks at the seasonal cycle frequency and at its higher harmonics, and, thus, it corresponds to a multi-bandpass filter in this case. In consequence of this, a certain amount of variance is reallocated from the abovementioned white noise uncertainty component to the first component, although it appears there at the lowest frequencies because of spectral aliasing. Thus, proxy seasonality may generate uncertainties that are highly correlated in the time domain. In most cases, this low-frequency uncertainty will be dominated by the seasonal cycle and its amplitude modulation caused by orbital variations (as demonstrated by Huybers and Wunsch, 2003, for example). Nonetheless, if the stochastic climate variability is only weakly red such that it is associated with notable power near the seasonal cycle frequency, it may also give rise to low-frequency uncertainties, in particular if the seasonal cycle is weak by comparison.

• If, in addition, the proxy abundance window is known to have a preferred seasonal timing throughout the year, then the contribution that the seasonal cycle signal (with its deterministic phase) makes to both of the abovementioned two uncertainty components is further modified. The white noise component can be larger or smaller than for random seasonal timing, and, in particular, the first uncertainty component may include a (potentially time-varying) deterministic bias in this case. Moreover, the sum of their variances may change because of the in-phase subsampling from a deterministic signal.

• Uncertainties caused by laboratory measurement errors are independent of the abovementioned components, and, thus, the associated power spectral density can simply be added to the error power spectrum obtained from our method. In practice, this uncertainty component is assumed to be white noise such that it scales inversely with any averaging timescale.

Another interesting and future application of the derived analytic expressions would be the inference of the power spectrum of the true climate signal. Specifically, by setting the reference climate in our method to zero and then repeating the entire derivation, one obtains the analytic expressions for the power spectrum of the climate reconstruction itself rather than of its error. Thus, one obtains an operator that transforms the power spectrum of the supposed true climate signal into a spectrum subject to the distortions caused by the processes included in our reconstruction uncertainty model. Then, given all of the parameters of the uncertainty model, and assuming a parametric form for the true climate signal, it might be possible to estimate its parameters by means of a maximum likelihood approach (that investigates the likelihood, under a given set of parameters, of the power spectrum estimated from a specific proxy record). This essentially amounts to inverting the aforementioned operator, which is similar to the correction technique used by Laepple and Huybers (2013) and motivated by the anti-aliasing approach of Kirchner (2005).

Appendix A: Nonstationary variance of the white noise component $B_n^{\left(j\right)}$

The variance of $B_n^{\left(j\right)}$ is given by its auto-covariance function at lag zero, ${\text{Cov}}_{\left(\mathit{ϵ}\right)}\left({B_n^{\left(j\right)}}^{\star },B_{n^{\prime }}^{\left(j\right)}\right){|}_{n=n{}^{\prime }}$. By substitution from Eq. (73), it can be shown, after some algebraic transformations, that

$$\begin{array}{}\text{(A1)}& {\left.{\text{Cov}}_{\left(\mathit{ϵ}\right)}\left({B_n^{\left(j\right)}}^{\star },B_{n^{\prime }}^{\left(j\right)}\right)\right|}_{n=n{}^{\prime }}=\mathcal{R}+{\mathcal{R}}^{\star },\end{array}$$

with

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& \mathcal{R}={\mathit{\alpha }}_{\mathrm{c}}^{\mathrm{2}}\left(R_{\mathrm{1}}+R_{\mathrm{2}}\right)+\mathrm{4}{\mathit{\alpha }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{a}}\left(R_{\mathrm{3}}+R_{\mathrm{4}}\right)\\ & +\mathrm{2}{\mathit{\alpha }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{a}}^{\mathrm{2}}\left(R_{\mathrm{5}}+R_{\mathrm{6}}+R_{\mathrm{7}}+R_{\mathrm{8}}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& R_{\mathrm{1}}=D_{\mathrm{1}},R_{\mathrm{2}}=e^{-i\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}}{D}_{\mathrm{2}},R_{\mathrm{3}}=e^{i\left(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}}\right)}{D}_{\mathrm{3}},\\ & R_{\mathrm{4}}=\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\right)e^{i\left(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}}\right)}{D}_{\mathrm{4}},R_{\mathrm{5}}=D_{\mathrm{5}},\\ & R_{\mathrm{6}}=\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\right)D_{\mathrm{6}},R_{\mathrm{7}}=e^{i\mathrm{2}\left(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}}\right)}{D}_{\mathrm{7}},\\ & R_{\mathrm{8}}=\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\right)e^{i\mathrm{2}\left(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}}\right)}{D}_{\mathrm{8}},\end{array}\end{array}$$

and the characteristic function differences D₁ to D₈, also using definition (33), are given by

$$\begin{array}{}\text{(A4)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{1}}& =C_{n,n{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}},{\mathit{\nu }}_{\mathrm{c}}\right)\right|}_{n=n{}^{\prime }}-C_{l,l^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}},{\mathit{\nu }}_{\mathrm{c}}\right)\right|}_{l\ne l^{\prime }}\\ & =\mathrm{1}-{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right),\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{2}}& =C_{n,n{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}},-{\mathit{\nu }}_{\mathrm{c}}\right)\right|}_{n=n{}^{\prime }}-C_{l,l{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}},-{\mathit{\nu }}_{\mathrm{c}}\right)\right|}_{l\ne l{}^{\prime }}\\ & =\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)-{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right),\end{array}\text{(A6)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{3}}& =C_{n,n{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}},{\mathit{\nu }}_{\mathrm{c}}\right)\right|}_{n=n{}^{\prime }}\\ & -C_{l,l{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}},{\mathit{\nu }}_{\mathrm{c}}\right)\right|}_{l\ne l{}^{\prime }}\\ & ={\widehat{f}}_{\mathrm{bs}}^{\star }\left({\mathit{\nu }}_{\mathrm{a}}\right)[\mathrm{1}-{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)],\end{array}\text{(A7)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{4}}& =C_{n,n{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}},-{\mathit{\nu }}_{\mathrm{c}}\right)\right|}_{n=n{}^{\prime }}\\ & -C_{l,l{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}},-{\mathit{\nu }}_{\mathrm{c}}\right)\right|}_{l\ne l{}^{\prime }}\\ & ={\widehat{f}}_{\mathrm{bs}}^{\star }\left({\mathit{\nu }}_{\mathrm{a}}\right)[\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)-{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)],\end{array}\text{(A8)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{5}}& =C_{n,n{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}-{\mathit{\nu }}_{\mathrm{a}},{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}}\right)\right|}_{n=n{}^{\prime }}\\ & -C_{l,l{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}-{\mathit{\nu }}_{\mathrm{a}},{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}}\right)\right|}_{l\ne l{}^{\prime }}\\ & =\mathrm{1}-{\left|{\widehat{f}}_{\mathrm{bs}}^{\star }\left({\mathit{\nu }}_{\mathrm{a}}\right)\right|}^{\mathrm{2}}{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right),\end{array}\text{(A9)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{6}}& =C_{n,n{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}-{\mathit{\nu }}_{\mathrm{a}},-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}}\right)\right|}_{n=n{}^{\prime }}\\ & -C_{l,l{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}-{\mathit{\nu }}_{\mathrm{a}},-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}}\right)\right|}_{l\ne l{}^{\prime }}\\ & =\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)-{\left|{\widehat{f}}_{\mathrm{bs}}^{\star }\left({\mathit{\nu }}_{\mathrm{a}}\right)\right|}^{\mathrm{2}}{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right),\end{array}\text{(A10)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{7}}& =C_{n,n{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}},{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}}\right)\right|}_{n=n{}^{\prime }}\\ & -C_{l,l{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}},{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}}\right)\right|}_{l\ne l{}^{\prime }}\\ & ={\widehat{f}}_{\mathrm{bs}}^{\star }\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{a}}\right)-{\widehat{f}}_{\mathrm{bs}}^{\star \mathrm{2}}\left({\mathit{\nu }}_{\mathrm{a}}\right){\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right),\end{array}\text{(A11)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{8}}& =C_{n,n{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}},-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}}\right)\right|}_{n=n{}^{\prime }}\\ & -C_{l,l{}^{\prime }}{\left.\left(-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}},-{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{a}}\right)\right|}_{l\ne l{}^{\prime }}\\ & ={\widehat{f}}_{\mathrm{bs}}^{\star }\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{a}}\right)\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)-{\widehat{f}}_{\mathrm{bs}}^{\star \mathrm{2}}\left({\mathit{\nu }}_{\mathrm{a}}\right){\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right).\end{array}\end{array}$$

Since the auto-covariance contributions R₃, R₄, R₇, and R₈ depend on t_n, the variance of $B_n^{\left(j\right)}$ is nonstationary. Furthermore, it turns out that with $n\ne n{}^{\prime }$ the characteristic function differences D₁ to D₈ are all zero, and, thus, the auto-covariance contributions R₁ to R₈ are all zero. This implies that the auto-covariance function of $B_n^{\left(j\right)}$ is nonzero only at lag zero ($n=n{}^{\prime }$) and zero at all other lags ($n\ne n{}^{\prime }$). Hence, $B_n^{\left(j\right)}$ is a white noise process, and from its definition (Eq. 73) it follows that it has zero mean. Note that Eq. (A4) is identical to Eq. (49) in Sect. 3.1 but with ν=ν_c, and so the abovementioned procedure follows the same key idea (according to the approach of Balakrishnan, 1962) of extrapolating the auto-covariance function from nonzero lags towards lag zero.

From the abovementioned expressions, we can write the variance of $B_n^{\left(j\right)}$ as

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{\text{Var}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right)& ={\mathcal{V}}_B^{\left(\mathrm{0}\right)}+{\mathcal{V}}_{B,n}^{\left(\mathrm{1}\right)}\cos \left(\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+{\mathit{\varphi }}_{\mathrm{a}}+{\mathit{\varphi }}_{\mathrm{b}\mathrm{1}}\right)\\ & +{\mathcal{V}}_{B,n}^{\left({\mathrm{2}}^{\prime }\right)}\cos \left(\mathrm{4}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+\mathrm{2}{\mathit{\varphi }}_{\mathrm{a}}+{\mathit{\varphi }}_{\mathrm{b}\mathrm{2}}\right)\\ & +{\mathcal{V}}_{B,n}^{\left(\mathrm{2}{}^{\prime \prime }\right)}\cos \left(\mathrm{4}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{t}_n+\mathrm{2}{\mathit{\varphi }}_{\mathrm{a}}+{\mathit{\varphi }}_{\mathrm{b}\mathrm{1}}\right),\end{array}\end{array}$$

with the amplitude of the stationary variance component,

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{\mathcal{V}}_B^{\left(\mathrm{0}\right)}& ={\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}\mathit{\{}\mathrm{1}-{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)+\cos (\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}})\\ & [\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)-{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)]\mathit{\}}+{\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}\\ & \mathit{\{}\mathrm{1}-M_{\mathrm{b}\mathrm{1}}^{\mathrm{2}}{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right){\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)+\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\right)\\ & [\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)-M_{\mathrm{b}\mathrm{1}}^{\mathrm{2}}{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right){\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)]\mathit{\}},\end{array}\end{array}$$

the amplitude of the variance component oscillating at frequency ν_a,

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}{\mathcal{V}}_{B,n}^{\left(\mathrm{1}\right)}& =\mathrm{2}{\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{a}}\sqrt{\mathrm{2}}\mathit{\left\{}{M}_{\mathrm{b}\mathrm{1}}\text{sinc}\left({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right)\right[\mathrm{1}-{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)]\\ & +\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\right)M_{\mathrm{b}\mathrm{1}}\text{sinc}\left({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right)\\ & [\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)-{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)]\mathit{\}},\end{array}\end{array}$$

and the amplitudes of the variance components oscillating at frequency 2ν_a,

$$\begin{array}{}\text{(A15)}& {\displaystyle }& {\displaystyle }{\mathcal{V}}_{B,n}^{\left(\mathrm{2}{}^{\prime }\right)}={\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}\mathit{\left\{}{M}_{\mathrm{b}\mathrm{2}}\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right)\right[\mathrm{1}+\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\right)\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)\left]\mathit{\right\}},\text{(A16)}& {\displaystyle }& {\displaystyle }{\mathcal{V}}_{B,n}^{\left({\mathrm{2}}^{\prime \prime }\right)}=-{\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}\mathit{\left\{}{M}_{\mathrm{b}\mathrm{1}}^{\mathrm{2}}{\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{s}}\right){\text{sinc}}^{\mathrm{2}}\left({\mathit{\nu }}_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{p}}\right)\right[\mathrm{1}+\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\right)\left]\mathit{\right\}},\end{array}$$

and where

$$\begin{array}{}\text{(A17)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{b}\mathrm{2}}=\arg \left[{\widehat{f}}_{\mathrm{b}}^{\star }\right(\mathrm{2}{\mathit{\nu }}_{\mathrm{a}}\left)\right]=\mathrm{4}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{b}}-\arctan \left(\mathrm{4}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{b}}\right)\end{array}$$

and

$$\begin{array}{}\text{(A18)}& {\displaystyle }{M}_{\mathrm{b}\mathrm{2}}=\left|{\widehat{f}}_{\mathrm{b}}^{\star }\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{a}}\right)\right|={\left[\mathrm{1}+{\left(\mathrm{4}\mathit{\pi }{\mathit{\nu }}_{\mathrm{a}}{\mathit{\tau }}_{\mathrm{b}}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}/\mathrm{2}},\end{array}$$

and ϕ_b1 and M_b1 are defined by Eqs. (77) and (78), respectively. The time average of this variance over an infinitely long time interval is then given by ${\mathcal{V}}_B^{\left(\mathrm{0}\right)}$, provided that Δt is not a multiple of ${\mathit{\nu }}_{\mathrm{a}}^{-\mathrm{1}}$. If the time average is taken over a finite time interval of length T, centered at t=0, the time mean variance is given by

$$\begin{array}{}\text{(A19)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{\text{Var}}}_{\left(\mathit{ϵ}\right)}\left(B_n^{\left(j\right)}\right)& ={\mathcal{V}}_B^{\left(\mathrm{0}\right)}+{\mathcal{V}}_{B,n}^{\left(\mathrm{1}\right)}\cos \left({\mathit{\varphi }}_{\mathrm{a}}+{\mathit{\varphi }}_{\mathrm{b}\mathrm{1}}\right)\text{sinc}\left({\mathit{\nu }}_{\mathrm{a}}T\right)\\ & +{\mathcal{V}}_{B,n}^{\left({\mathrm{2}}^{\prime }\right)}\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{a}}+{\mathit{\varphi }}_{\mathrm{b}\mathrm{2}}\right)\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{a}}T\right)\\ & +{\mathcal{V}}_{B,n}^{\left(\mathrm{2}{}^{\prime \prime }\right)}\cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{a}}+{\mathit{\varphi }}_{\mathrm{b}\mathrm{1}}\right)\text{sinc}\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{a}}T\right),\end{array}\end{array}$$

provided that $\mathrm{\Delta }t\ll {\mathit{\nu }}_{\mathrm{a}}^{-\mathrm{1}}$.

Appendix B: Parameters of the reconstruction uncertainty model as defined in Sect. 2

Parameter	Symbol
Seasonal cycle variance	${\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}$
Seasonal cycle frequency	νc
Expected seasonal cycle phase	$\langle {\mathit{\varphi }}_{\mathrm{c}}{\rangle }_{{\mathit{\varphi }}_{\mathrm{c}}}$
Seasonal phase uncertainty	${\mathrm{\Delta }}_{{\mathit{\varphi }}_{\mathrm{c}}}$
Amplitude modulation variance	${\mathit{\sigma }}_{\mathrm{a}}^{\mathrm{2}}$
Amplitude modulation frequency	νa
Amplitude modulation phase	ϕa
Proxy abundance timescale	τp
Bioturbation timescale	τb
Sediment sampling timescale	τs
Sampling interval	Δt
Length of proxy record	T
Number of signal carriers	N
Measurement error variance	${\mathit{\sigma }}_{\mathit{\mu }}^{\mathrm{2}}$
Reference climate averaging timescale	τr