A spectral approach to estimating the timescale-dependent uncertainty of paleoclimate records – Part 2: Application and interpretation

Proxy climate records are an invaluable source of information about the earth’s climate prior to the instrumental record. The temporal and spatial coverage of records continues to increase; however, these records of past climate are associated with significant uncertainties due to non-climate processes that influence the recorded and measured proxy values. Generally, these uncertainties are timescale dependent and correlated in time. Accounting for structure in the errors is essential for providing realistic error estimates for smoothed or stacked records, detecting anomalies, and identifying trends, but this structure is seldom accounted for. In the first of these companion articles, we outlined a theoretical framework for handling proxy uncertainties by deriving the power spectrum of proxy error components from which it is possible to obtain timescale-dependent error estimates. Here in Part 2, we demonstrate the practical application of this theoretical framework using the example of marine sediment cores. We consider how to obtain estimates for the required parameters and give examples of the application of this approach for typical marine sediment proxy records. Our new approach of estimating and providing timescale-dependent proxy errors overcomes the limitations of simplistic singlevalue error estimates. We aim to provide the conceptual basis for a more quantitative use of paleo-records for applications such as model–data comparison, regional and global synthesis of past climate states, and data assimilation.

all the information required to derive timescale dependent uncertainties, and working in the frequency domain further simplifies the estimation of the different error components. A number of additional useful quantities such as the uncertainty in a timeslice mean, the uncertainty in the difference between two time-slices and the expected timescale-dependent correlation between replicates of proxy records and between proxy records and the true climate, can easily be derived from the error spectrum.
In part I we discuss the theoretical basis for and give a full mathematical derivation of the Proxy Spectral Error Model 30 (PSEM). Here in part II we aim to facilitate the use of PSEM in paleo-climate applications. Thus, we 1) sketch the concepts behind the different error components in a more applied way, 2) provide heuristic approaches to parametrise the climate spectra and other parameters of the error model, 3) provide examples using virtual and actual sediment cores, and 4) provide an Rpackage implementing the spectral error estimation method. 2 Error spectra as a framework for timescale-dependent proxy uncertainty In this section we illustrate the error spectrum framework for the specific example of temperature related proxies in the shells of planktic foraminifera recovered from marine sediment cores. The major processes contributing uncertainty to these proxy records have been explored using physically motivated proxy forward models that simulate pseudo-proxies from an assumed true input climate (Laepple and Huybers, 2013;Dolman and Laepple, 2018). These processes include seasonality in the creation 5 of proxy signal carriers (e.g. foraminiferal tests), aliasing due to under-sampling of the seasonal climate cycle, mixing and smoothing of the signal due to bioturbation, and independent measurement and processing error. One approach to estimating the uncertainty for a given metric and proxy record is to use such a forward model to simulate ensembles of pseudo-proxy records, calculate the metric for each and then examine their statistical properties. In the approach we propose here, we do not simulate pseudo-proxies for a specific climate time-series, rather we make some simplifying assumptions about the properties Frequency [year −1 ] Power spectral density [K 2 year] Figure 2. An example spliced empirical and power-law power spectrum of ocean temperature (0-120 m water depth) for 20°N latitude. A small discontinuity at ν = 0.03 is visible, as to increase the robustness of the intercept estimate, the splicing is implemented by matching the integrated spectra between ν = 0.03 and ν = 0.1. S(ν) (see sections 2.1 and 3.1 of part I). While the developed approach allows for any choice of a climate spectrum, e.g. from complex climate models or estimated from observations, we here outline a heuristic method suitable for marine sea-surfacetemperature (SST) and surface δ 18 O calcite records. This method is implemented in the PSEM R package and requires only the core position and habitat depth range as input parameters.
A detailed, site-specific power spectrum of the climate can only easily be estimated empirically for timescales up to the 5 length of the instrumental record. At longer time scales the climate spectrum can be approximated by a power-law type scaling, S(ν) = αν −β , where the exponent β characterises the scaling behaviour and is thought to lie between 0 and 2 (Lovejoy, 2015;Schmitt et al., 1995) and α scales the amplitude of climate variation. Here we take the pragmatic approach of splicing zonally averaged empirical climate spectra for frequencies above 1/33 years with theoretical power law spectra at lower frequencies (Laepple and Huybers, 2014b) (Fig. 2). As the empirical spectra were estimated from annual resolution ocean temperature 10 records, we set power to zero at frequencies above 1/2 years.
For the power-law section of the climate spectra α was chosen so that the low frequency power law spectra are continuous with the empirically estimated high frequency regions of the spectra. We typically assume a value of 1 for β as this was found to be a good description of Holocene SST variability (Laepple and Huybers, 2014a) but this parameter can be freely specified and we illustrate the affect of varying this parameter in Appendix A. To allow these spectra to also be used for δ 18 O records 15 we re-calibrate them to δ 18 O calcite units using a standard calibration (Bemis et al., 1998), which in terms of variance is effectively just a division by a factor of 4.8 2 , assuming that the δ 18 O calcite is generally dominated by temperature variations at these timescales. A function to generate these spliced empirical and power-law spectra is supplied as part of the PSEM R package.
To this stochastic climate we add a deterministic seasonal cycle modeled as the power spectrum of a sine wave with frequency 1/1 year (see section 3.2 in part I). For a given location, we parametrise the amplitude of the seasonal cycle using a gridded data set of assimilated physically consistent δ 18 O sw and temperature (Breitkreuz et al., 2018). δ 18 O calcite was calculated on the PDB scale from δ 18 O sw and temperature using the equations of Shackleton (1974).
Additionally, the amplitude of earth's seasonal temperature cycle varies over a precessionary orbital cycle with approximate 5 frequency 1/23 kyr. The magnitude of this amplitude modulation depends primarily on latitude and we assume this is equal to the amplitude modulation of incoming solar radiation at a given latitude (Berger and Loutre, 1991). This introduces a low frequency (1/23 kyr) deterministic signal to the climate.

Proxy error processes as filters.
During the creation of a climate proxy record, some of the processes that introduce errors, defined here as differences between 10 the proxy record and the true climate, can be thought of as acting like filters on the true climate signal. Here we illustrate the concept of PSEM by considering the smoothing effects of bioturbation and the width of sediment slices from which signal carriers, e.g. foraminiferal shells, are extracted.
Bioturbation at the water-sediment interface mixes the upper few centimeters of sediment thereby mixing together signal carriers of different ages. This acts like a smoothing filter on the climate signal, reducing the amplitude of climate variations in 15 a frequency dependent way. The magnitude of the reduction at each frequency depends on the filter characteristics; using the simple physical bioturbation model of Berger and Heath (1968), the filter width, τ b , is simply the bioturbation depth divided by the sedimentation rate. In the frequency domain, this is equivalent to multiplying the power spectrum of the climate with the transfer function of the filter and results in a power spectrum of the error as described in part I, sections 2.2 and 3.1.
Similarly, when a sediment core is sampled the proxy signal carriers are picked from a series of slices of sediment, each with 20 a finite width of typically 1-2 cm. This again acts as a filter, this time a running mean filter with width, τ s , equal to the ratio of the slice thickness and sedimentation rate. As for the bioturbation filter, in the frequency domain the effect on the original signal is obtained by multiplying the power spectrum with the transfer function of the filter (see sections 2.3 and 3.1 in part I).
2.2.1 Error relative to the reference climate.
So far we refer to a proxy error as a difference between the measured proxy value and the "true" climate signal. The bioturbation 25 and slice thickness filters smooth the climate signal so that the proxy differs from the true climate. However, in practice, values of proxy variables from marine sediments are rarely interpreted as representing the instantaneous climate state, it is understood implicitly that some smoothing has taken place. The error therefore depends on the interpreted timescale of the proxy record and so, in the spectral error model, we define error relative to the true climate at a specific timescale provided by the user (see section 2.4 in part I). The power-spectrum of this error is shown as the brown dashed line in Fig. 3a Figure 3. A conceptual representation of PSEM. a) The true climate signal is filtered (smoothed) by processes such as bioturbation. This modifies the power spectrum of the climate (red) in a frequency dependent way, producing the power spectrum of the climate signal after bioturbation (blue). Proxy records are assumed to represent the climate at a particular timescale, (e.g. centennial, millennial), the reference climate spectrum (purple) is the power spectrum of the true climate smoothed to this timescale. The error that bioturbation and other smoothing produces (dashed brown) is a function of the reference and bioturbated climate spectra. b) The power of the true climate signal that is removed by bioturbation and other smoothing processes (grey shaded region) can be redistributed as a white noise error component ( (Bijma et al., 1990;Spero, 1998). The bioturbation and slice thickness filters can be thought of as probability density functions (PDFs) that describe the portion of time from which the climate signal is sampled by these N signal 5 carriers. As this is a finite sample, the resulting proxy value is an estimate of the mean value and contains a stochastic noise component in addition to the deterministic error caused by the smoothing. The variance of this stochastic error term is equal to the integral of the difference between the climate spectrum and the spectrum of the smoothed climate signal, divided by the number of individual signal carriers in a sample, shown as the grey line and shaded area in Fig. 3b (see sections 2.3 and 3.1 in part I). As N tends to infinity, the error due to undersampling tends toward zero. This may be the case for proxies such as 10 Uk'37, where measurements consists of many thousands of organic molecules.

The seasonal cycle
The often large cyclical variation in climate variables associated with the seasonal cycle can also add noise to a proxy record if the seasonal cycle is not adequately sampled by individual signal carriers, each of which records a short portion of the cycle (Laepple and Huybers, 2013;Schiffelbein and Hills, 1984). Additionally, a bias can be introduced in the record if the signal 15 carriers are produced in greater numbers during a particular part of the year (Jonkers and Kučera, 2015;Leduc et al., 2010). In the case of orbital modulation of the amplitude of the seasonal cycle, this bias becomes a slowly varying error, and the variance of the noise process also varies over the course of an orbital cycle (see section 3.2 of part I).
In the spectral error model, the seasonal cycle is represented by the discrete power spectrum of a deterministic sine wave (Eq.1 in part I). The interaction between this signal and seasonality in the production of signal recorders determines the 20 magnitude of two errors -a white noise error generated by undersampling of the seasonal cycle and a bias due to only sampling a portion of the seasonal cycle.
We represent seasonality in the production of signal carriers here by saying that production occurs continuously over a set fraction of the year, τ p (see section 2.2 in part I). This can also be viewed as a kind of filter, this time on the discrete spectrum of the deterministic sine wave seasonal cycle. The transfer function of this filter can be constructed in an analogous way to 25 those for bioturbation and slice thickness, and the difference between the filtered and unfiltered spectra gives the error due to sampling only part of the seasonal cycle (see section 3.2 in part I). Again, the finite time-period each carrier records means that the portion of the seasonal cycle during which signal carriers are created is sampled by the individuals and this generates additional redistributed white noise in the proxy signal.
In the absence of orbital modulation, four parameters determine the errors generated by filtering and sampling the seasonal σ 2 c , the variance of the full seasonal cycle; τ p , the proportion of the seasonal cycle during which signal carriers are produced; φ c , the expected midpoint (phase) of the signal carrier production; and ∆φ c , which represents uncertainty in the phase of the carrier production (section 2.5 in part I).
If the signal carriers are produced all year round, τ p = 1, there is no bias and the white noise component has a variance equal to the variance of the full seasonal cycle divided by the number of signal carriers per sample, N . 5 If the signal carriers are produced for only part of the year, τ p < 1, but with completely unknown timing (we don't know which months, ∆φ c = 2π), then the expected variance of this white noise is equal to the difference between the variance of the full seasonal cycle and the variance of the seasonal cycle filtered with a running mean filter of width, τ p . In the spectral domain, this is analogous to the grey shaded area in Fig. 3(a), but this time for the discrete spectrum of the deterministic sine wave seasonal cycle. In this situation the sign of the seasonal bias is unknown. It appears in the error spectrum as power at 10 frequency zero.
If the timing of the production phase of the signal carriers is known precisely, ∆φ c = 0, then the white noise variance is equal to the variance of the piece of the sine wave that describes the portion of the year in which the carriers are produced.
In this case the sign and value of the bias are completely "known" given the parameters of the model, or put another way, the proxy record can be attributed to the correct season.

15
PSEM can handle intermediate situations where, for example, we can parametrise with a proxy season length, τ p ; an expected phase, φ c , which is the midpoint of the proxy production season; and an uncertainty in this phase, 0 < ∆φ c < 2π. In this situation there is both a bias with an expected value and sign, and an uncertainty around this expected value that comes from the phase uncertainty.
Finally, if we include amplitude modulation of the seasonal cycle over the course of a precessionary cycle, σ 2 a > 0, the 20 size of any seasonal bias and bias uncertainty will vary over time. In the spectral domain this manifests as leakage of power from frequencies lower than 1 / T to higher frequencies and creates additional timescale-dependent errors corresponding to uncertain changes in the magnitude of seasonal biases between time periods. Additionally, the magnitude of aliased seasonal cycle variation will vary over a precessionary period (see section 3.2 of part I).

Measurement error and individual variation 25
As, to a first order, the measurement error can be assumed to be independent between measurements, we simply add the power spectrum of a white noise error term σ meas . More complex measurement errors such as machine drift or memory between measurements could be integrated by adding a power spectrum characterising these machine characteristics. We add an additional error term, σ ind , to account for inter-individual variation in the encoded signal. This is a catch-all term to include things like differences in depth habitat occupied by individuals, variation in the encoding of the signal, i.e. "vital effects" (Haarmann et al., 30 2011;Schiffelbein and Hills, 1984;Sadekov et al., 2008;Duplessy et al., 1970). For a given proxy measurement, the variance of this term is scaled by the number of individual signal carriers in the sample, N .

Calibration error
Finally we add uncertainty in proxy's calibration to temperature as a constant error, applying to all values in a given record, for which we don't know the sign but do have some idea of the magnitude σ cal . This could for example be the standard error of the intercept term in a linear calibration model. This is implemented as additional power at frequency zero.
2.7 Power spectrum of the total error 5 As the individual error components are independent, their power spectra can be added together to get the spectrum of the total error. Once the spectral error model has been parametrised for a given core, proxy, and sampling scheme, the resulting empirical power spectrum provides the magnitude and full temporal correlation structure of the error components. From this a number of useful quantities can be obtained. These include the error on individual proxy measurements, the error after smoothing a record to a lower resolution, the error on a time-slice and on the difference between two time-slices, the expected correlation 10 between replicated proxy records. We illustrate these applications in the following two sections.
3 Illustration of the error spectrum approach for a hypothetical proxy record.
We first illustrate the error spectrum approach for a hypothetical 10 kyr foraminiferal Mg/Ca record. Parameter values have been chosen to be realistic while ensuring that all components of the error model are presented. We parametrise the climate spectrum as described in section 2.1, assuming a surface dwelling foraminifera at a virtual core position of 20°N, -18 E, 15 calcifying between the surface down to 120 m (Fig. 2). We further assume that this taxon forms tests for a 7 month period of the year, centered around the peak of the seasonal cycle but with an uncertainty in this phase of 2 months in either direction.
A bioturbation depth of 10 cm and sedimentation rate of 10 cm kyr −1 are assumed. Thirty foraminifera are picked from contiguous 1 cm thick sediment slices so that the resulting record has a sampling interval, ∆t, of 100 years. We assume a measurement error of 0.25°C and inter-individual variation of 1°C. All parameters are given in Table 1. 20 The power spectra of the individual error components, together with their sum and the assumed power spectrum of the climate, are shown for this example parameter set in Fig. 4. Power at ν = 0 corresponds to errors that are constant for a given proxy time-series and thus do not shrink as additional measurements are averaged together. Here it is composed of those parts of the seasonal bias and bias uncertainty that are constant over time (i.e. not orbitally modulated), and the calibration uncertainty.
The power spectral densities of the measurement error and individual variation components are horizontal lines, indicating 25 that their power is independent of frequency, i.e. that these error components have the property of white noise. This applies also to the two components of aliased / redistributed seasonal cycle and stochastic climate variation.
In contrast, the component due to bioturbation and sediment slice thickness smoothing shows strong frequency dependence.
The error due to smoothing is proportional to the variation in the climate signal, therefore, as variation in the climate is larger at lower frequencies, the error due to smoothing also increases towards lower frequencies. This is true up until the point at which 30  The individual values in record are most likely interpreted as a kind of mean of the time interval between observations and so we set the implicit reference timescale to be the same as the sampling resolution (∆t = 100). For this example, the timescale of the bioturbation smoothing, τ b , is 1000 years. The large difference between these timescales implies a large error due to bioturbation smoothing. At frequencies above about 2 x the bioturbation filter width, the power of the bioturbation error is equal to the power of the reference climate.

5
The timescale-dependent portions of the seasonal bias and seasonal bias uncertainty are due to orbital modulation of the amplitude of the seasonal cycle. As this is an approximately 23 kyr cycle, these errors only become large at timescales approaching 23 kyr.

Timescale-dependent proxy error
The error for the individual proxy values at their original sampled resolution can be obtained by integrating the error spectrum. 10 When the record is smoothed before its interpretation, the timescale represented by each point changes as does the error. The error for a given timescale can be obtained by integrating the error spectra after first multiplying with the transfer function of the smoothing filter (see section 4, Eq. 110 of part I).
Timescale-dependent error for the example parameter set is shown in Fig. 5 for timescales from 100 years (the original sampling resolution of the record) to 10000 years (a mean or time-slice of the entire length of the record), assuming a running

Error on a time-slice mean and the difference between two time-slices.
The information contained in the power spectra also allows us to estimate the uncertainty in the difference in the climate between two timepoints, or between two time-slices. A "time-slice" refers to an average taken over a set of proxy values  Figure 6. The uncertainty or error on the estimate of the mean climate over two time-slices covering the first and last 1000 years of the pseudo-proxy record, and the error in the estimate of the difference between these two time-slices. The realistic error estimate using PSEM (third column) is much smaller than the naive error estimate that one would obtain by just adding up the variances. within a certain time period of interest. For example, using the parametrisation from Table 1, if we wanted to compare the mean climate over the first and last 1000 years the proxy record, the error variance on each individual time-slice would be the value at timescale = 1000 years in Fig. 5. If all the errors were independent in time then the error variance on the difference between these two time-slices would simply be the sum of the two variances. However, as some of these error components are autocorrelated (or even constant over the entire time-series) the covariance in the errors for the two time-slices needs to be 5 accounted for. The information to do this is contained in the power spectrum of the error (see section 4, Eq. 111-112 of part I). The uncertainty, or error, on the estimate of the difference between two time-slices is much smaller than the errors on the time-slices themselves (Fig. 6).

A worked example: replicated Holocene δ 18 O records from South Java
Here we illustrate the use of the error spectrum model on a real proxy record by applying it to replicated foraminiferal δ 18 O 10 records taken from core GeoB 10054-4 off South Java in the Indian Ocean collected during the R/V SONNE 184 expedition (8°40.90'S, 112°40.10' E; 1076 m water depth; Hebbeln and cruise participants, 2006). Replicated δ 18 O records were created using two different sampling schemes. In record 1 (Rep1) measurements were made on samples consisting of 5 Globigerinoides ruber(s.s.) tests each, at a mean interval of 83 years; in Rep2, 30 tests of G. ruber(s.s.) were used per sample at a mean interval of 246 years. An age model was constructed using 9 AMS-14 C dates on mono-specific samples of Trilobatus sacculifer (see 15 supplement S1). The Marine13 radiocarbon calibration curve was used to calibrate the ages and construct a linear age model (Reimer et al., 2013). As both records come from the same core, more advanced age-depth modelling is not required here.
A replicated set of radiocarbon dates taken on ten samples each of 10 foraminifera indicated an inter-individual standard deviation in age of 720 years  (S1), which we use for τ b , corresponding to a bioturbation depth of about 14 cm.

Parametrisation
To assist potential users of PSEM here we describe the parameter choices step by step. 5 Formally, the spectral error model describes only regularly sampled time-series whose total length T is an odd multiple of the sampling interval ∆t. When applying the spectral error model to real proxy time-series we have to make some additional approximations to accommodate their inevitably irregular (in time) sampling intervals. Hence we approximate the sampling interval for each record ∆t by ∆t.
For the climate spectrum we again used a spliced empirical and theoretical spectrum and estimated the amplitude of the 10 seasonal cycle as described in section 2.1, this time for the surface down to 50 m, resulting in an amplitude of 0.53 permille.
At this location G. ruber is thought to produce tests at an approximately equal rate throughout the year and represent the annual mean surface temperature in this region . We therefore set τ p = 1, which implies year-round production of signal carriers. As τ p = 1, the parameters φ c and ∆φ c , which control the phase of signal carrier production relative to the seasonal cycle, have no effect on the error spectrum. Similarly, orbital modulation of the amplitude of the 15 seasonal cycle will have only a small effect and is ignored here.
The sediment accumulation rate estimated from the calibrated radiocarbon dates and retrieval depths is approximately 20 cm/kyr. The sediment slices were 1 cm thick so τ s was set to 50 years. A replicated set of radiocarbon dates taken on ten samples each of 10 foraminifera indicated an inter-individual standard deviation in age of 720 years (unpublished data), which we use for τ b , corresponding to a bioturbation depth of about 14 cm. 20 For σ meas we use the analytic replicability of 0.1 permille. σ ind is more difficult to parametrise, we use 0.32 permille which was estimated by Sadekov et al. (2008) as the contribution of "vital effects" to replicability estimates for G. ruber.

Timescale-dependent uncertainty
Total error variance at the highest frequency resolved (∆t) is much higher for Rep1 than Rep2 as there are fewer foraminifera per sample so that the individual variation, aliased seasonal cycle and aliased stochastic climate components are all larger for 25 Rep1 than for Rep2 (Fig. 7). The effect of this can be seen clearly in Fig. 8(a) which shows Rep1 and Rep2 at their original irregularly sampled time-points together with their PSEM estimated uncertainties. Rep1 shows much higher variance despite the fact that Rep1 and Rep2 come from the same sediment core and therefore both experienced the same climate signal and the same degree of bioturbation smoothing.
In Fig. 8(b), Rep1 and Rep2 have both been interpolated and smoothed to a regular 492 year resolution (492 = 2 x ∆t for 30 Rep2). As the original time-series were irregular, a different number of proxy measurement now contribute to each mean value.
To account for this, we evaluate PSEM separately for each point, setting (∆t) to the timescale τ smooth divided by the number of original proxy measurements. After this smoothing, the two series are in much closer agreement. The error for Rep1 has  shrunk much more than that for Rep2. In fact, smoothed to a timescale of 492 years, the error on Rep1 is now smaller than that for Rep2, due to the larger number of proxy measurements contributing on average to each point in the smoothed series (dotted vertical line in Fig. 7).

Expected correlation with the true climate
Finally, the property of a proxy record we are perhaps most interested in is its correlation with the true climate. The proxy and 5 climate can have low correlation due to a combination of non-climate variation (noise) in the proxy record and because variation in the climate has been smoothed out in the proxy record. As the noise and degree of smoothing are timescale dependent, so to is the proxy-climate correlation. Using the power spectra of the errors and the assumed power spectrum of the climate signal we can calculate the expected timescale dependent correlation between the proxy and climate. We do not in general know the true climate and so cannot test the accuracy of the climate-proxy correlation estimate, however, we can also calculate the 10 expected correlation between replicate proxy records and use this as a partial test of the model, under the assumption that only the processes considered in PSEM affect the proxy record. The results (Fig. 9) indicate an increasing expected correlation from around 0 at centennial timescales to around 0.5 at millennial timescales. This is an upper bound estimate as the chronological uncertainty and other effects not considered here will further decrease the climate to proxy relationship.
We estimated the timescale dependent correlation between Rep1 and Rep2 using the R package corit (Reschke et al., 2019).

15
The irregular time-series were first interpolated to high resolution regular timeseries and then smoothed with a set of increasingly wide running mean filters before calculating the correlation between them. The observed correlation between Rep1 and Rep2 is somewhat higher than that estimated from the error spectra ( Fig. 9), perhaps indicating that we have assumed for example slightly too high measurement error, although it is unclear if this difference is statistically significant. At timescales above 1000 years estimates of the observed correlation become very variable due to there being very few effective data points left after smoothing.

5
Understanding the errors associated with climate proxies is an important task for paleoclimate research. Proxy errors, defined here as differences between the inferred climate and the unknown true climate, can be large and can thus strongly influence our understanding of past climate history and the functioning of the climate system. Many components of proxy error have a complex temporal autocorrelation structure, making them timescale dependent and a challenge to properly quantify and account for. The model introduced here (PSEM) and in the companion paper (Kunz et al., 2020) offers a rigorous and compact way to calculate and express this structure as error spectra, specifically here in this first version for marine sediment cores. Once defined, the error spectra can be used to calculate many quantities that will be useful to paleoclimate research. In addition to the 5 error on individual measurements, these include the error after smoothing the record, the error on time-slices and differences between time-slices, the expected correlation between replicates of a proxy record and between a record and the true climate. As with every model, some challenges remain, in particular how to deal rigorously with the irregularity of real proxy timeseries, the climate dependency of the habitat of the recording organisms and the error associated with the age-uncertainty.
Nevertheless, we argue that the spectral error approach represents a significant advance towards obtaining reliable uncertainty 10 estimates for all quantities of interest rather than single estimates applied uniformly to the record.
Although we have restricted the examples given here to the Holocene, PSEM can be applied to other periods, albeit with greater uncertainty in its parametrisation. Many of the error components, such as the bioturbation smoothing and seasonal aliasing, should remain approximately correct; however, if we include glacial-interglacial cycles there will be larger variations in the sedimentation rate, and the seasonality of the signal carriers (e.g. foraminifera) may change in unexpected ways. The 15 amplitude of the seasonal temperature cycle and the precession driven modulation of the seasonal cycle will vary with the longer inclination and eccentricity orbital cycles -although these changes are proportionally relatively small and deterministic. For the assumed stochastic climate spectrum, the key issue is the assumption of stationarity. If multiple glacial cycles are included then one could argue that the spectrum is again stationary and still dominated by a power-law type variation. It becomes more difficult to justify if only one glacial-interglacial is included. Nonetheless the PSEM approach should be a significant improvement over assuming independent errors. 5 Simulation based forward modelling approaches such as sedproxy (Dolman and Laepple, 2018) and PRYSM (Dee et al., 2015) could also be used to estimate these quantities by generating and summarizing over many simulated pseudo-proxy records. The advantages of PSEM are that it provides an analytic understanding of the timescale dependence of error components while retaining the mechanistic understanding of the proxy generating process and can make these uncertainty estimates rapidly for large sets of parameters. For example, to directly model the error for a wide range of potential sediment core 10 characteristics (sedimentation rate and bioturbation depth), with different sampling schemes and at different locations with differing climate properties such as seasonal cycle amplitude. This allows us to both better interpret existing proxy record and to optimize future field work to answer specific questions.
Beyond modelling errors, PSEM also facilitates the use of the proxy variability itself to make inferences about the climate system. It allows prediction of the variability observed in individual foraminifera assemblages (IFA) (e.g., Groeneveld et al.,15 2019; Thirumalai et al., 2019) and thus to directly test the sensitivity of IFA statistics on the sedimentation rate, seasonality or the spectrum of climate variability. Finally, PSEM provides the basis to develop spectral correction approaches that infer the climate spectrum from the corrupted and distorted proxy spectrum, building on the approaches previously proposed for simpler sediment models (Laepple and Huybers, 2013) or for aliasing only (Kirchner, 2005).
The PSEM version proposed here includes the sediment proxy processes described earlier for the proxy forward model 20 sedproxy (Dolman and Laepple, 2018) and represents a trade-off between complexity and completeness. For example, the interaction of seasonality in the recording and climate signal is the only slowly varying process so far included. However, the general formulation of the PSEM allows other processes to be added. For example, depending on the timescale of interest and the proxy type, other slowly varying processes such as long-term changes in seawater Mg/Ca (Medina-Elizalde et al., 2008), or long-term instrumental drift and memory effects of the measurement process, could be included by specifying the power 25 spectra. When accounting for these processes, the use of PSEM vs. the classical single value uncertainty approach becomes even more important.
Here, and in part I, we have defined analytical expressions specifically for sediment archived climate proxies; however, the approach is applicable to other proxy types as most proxy types experience similar error generating and distorting processes.
For example, smoothing also affects water isotopes measured in ice-cores via water vapor diffusion (Johnsen et al., 2000) 30 and geochemical indices measured in coral records (Gagan et al., 2012) via successive incremental calcification in corals.
These processes can also be expressed as filters on the climate signal and the power spectra of the errors they produce derived in a similar way. Thus, we hope that PSEM presents an important step towards providing more realistic error estimates for paleoclimate research.
Code and data availability. The Proxy Spectral Error Model (PSEM) is implemented as an R package. Its source code is available from the public git repository https://github.com/EarthSystemDiagnostics/psem. A snapshot of the R package at the time these analyses were performed is archived at Zenodo (Dolman, 2020). Down-core radiocarbon dates and δ 18 O isotope data for GeoB 10054-4 are provided in Supplement S1.
Appendix A: The influence of the assumed climate spectrum on estimates of proxy error 5 An advantage of the spectral error method is that it does not require a specific realization of a climate model -we do however need to make some assumptions about the statistical properties of the true climate trajectory, which we encode via an assumed power spectrum of the climate. Our approach is to use a composite spectrum where we splice an empirical spectrum derived from observations for frequencies above 1/30 years with a theoretical power-law spectrum for lower frequencies. In the examples presented in the main body of the manuscript we have assumed a slope of 1 for the power-law portion of the spectrum. 10 Here we test the sensitivity of the method to the slope assumption by re-evaluating the error spectra for the example in section 3 using slopes of 0.5 and 1.5 in addition to 1 (Fig. A1). As changing the climate spectrum slope primarily influences the power at low frequencies, the effect of this on proxy error depends on the scale of the bioturbation -if bioturbation is low, so that it integrates only short timescale climate fluctuations, then the influence of the climate slope will be low. To illustrate this, we additionally use values of 20, 200, and 2000 for the parameter τ b -corresponding to sedimentation rates of 5, 50 and 500 15 cm kyr −1 for sediment with a 10 cm bioturbation depth.
Increasing the slope of the stochastic climate spectrum increases the error components due to the smoothing of the climate signal by bioturbation and the aliasing of this filtered climate variation as a white noise error term (see Fig. A2). However, the magnitude of these effects depends additionally on the parameter τ b , controlling the timescale of the smoothing filter. If τ b is small, e.g. 20 years (Fig. A2 a,b,c), the increase in error is minimal as the change in climate power only starts at frequencies 20 below 1/30 years, which are barely influenced by the smoothing filter. The result of this is that increasing the climate spectrum slope shifts the point at which the power of the climate signal exceeds the total error spectrum to higher frequencies. For τ b = 20 years, this point shifts from approximately 1/1000 years with beta = 1, to 1/333 years for beta = 1.5. When beta = 0.5, climate power remains below the total error at all shown frequencies, even with very low bioturbation.
For larger values of τ b , the bioturbation filter smooths the additional climate variation down to lower frequencies, such that 25 the increase in bioturbation and aliasing error is larger, offsetting the increase in climate power. For τ b = 200 (the same as in the main MS) shifting beta from 1 to 1.5 moves this from approximately 1/1000 to 1/500 years (Fig. A2 d,e,f). For very large values of τ b corresponding to a sedimentation rate of 5 cm ka-1 and bioturbation depth of 10 cm (τ b = 2000), increased climate variation is almost completely offset by increased bioturbation smoothing error down to frequencies of 1/5000 years (Fig. A2 g,h,i).  Figure A1. Spliced empirical and power-law power spectrum of ocean temperature (0-120 m water depth) for 20°N latitude, with slopes of 0.5, 1 and 1.5 for the low-frequency power-law region of the spectrum.   Figure A2. A comparison of error spectra using different slopes for the power-law portion of the stochastic climate spectrum. Comparisons are made for slopes, β, of 0.5, 1 and 1.5, and for different amounts of bioturbation smoothing, with values for τ b of 20, 200 and 2000 years.