Speleothem stable isotope reference records for East-Central 1 Europe-Resampling sedimentary proxy records to get evenly 2 spaced time-series with spectral control 3 4

4 István Gábor Hatvani, Zoltán Kern, Szabolcs Leél-Őssy, Attila Demény 5 Institute for Geological and Geochemical Research, Research Centre for Astronomy and Earth Sciences, Hungarian 6 Academy of Sciences, Budapest, Budaörsi út 45, H-1112, Hungary 7 Eötvös Loránd University, Department of Physical and Applied Geology, Budapest, Pázmány Péter stny. 1/C, H-1117, 8 Hungary 9 10 *Correspondence to: István Gábor Hatvani (hatvaniig@gmail.com) 11 12 13 Abstract. Uneven spacing is a common feature of sedimentary paleoclimate records, in many cases causing difficulties 14 in the application of classical statistical and time series methods. Although special statistical tools do exist to assess 15 unevenly spaced data directly, the transformation of such data to a temporarily equidistant time series applicable to 16 commonly used statistical tools remains, however, an unachieved goal. The present paper, therefore, introduces an 17 approach to obtain evenly spaced time series (using cubic spline fitting) from unevenly spaced speleothem records with 18 the application of a spectral control to avoid spectral bias caused by interpolation and retain the original spectral 19 characteristics of the data. The methodology was applied to stable carbon and oxygen isotope records derived from two 20 stalagmites of the Baradla Cave (NE Hungary) dating back to the late 18 century; it was also applied to additional well21 dated and high-resolution stable isotope records from the Han-sur-Lesse Cave (Belgium). To show the benefit of these 22 equally spaced records to climate studies, their coherence with primary and complex climate indices is explored using 23 wavelet transform coherence and discussed. The results shed light on clear relationships with climate and NAO indices, 24 lending support to the approach utilized in this study. Moreover, these suggest that the Baradla composite stable isotope 25 data can serve as regional reference records for the past ~200 years. The equally spaced time series obtained, are available 26 at doi: 10.1594/PANGAEA.875917. 27 28


Introduction
With more than a hundred speleothem studies published every year, it is trivial to state that speleothems are one of the most important objects of paleoclimate research.Compared to other continental carbonate deposits, they are especially valuable since (i) they are formed in relatively protected cave environments that render late-stage alteration less frequent, (ii) they can be dated using numerical dating methods (e.g.U-Th series and radiocarbon dating), (iii) they can be sampled at high spatial and temporal resolution, and (iv) they provide a number of data that serve as proxies for climate conditions (e.g.textural characteristics, trace element and stable isotope compositions, for further details see the comprehensive review by (Fairchild and Baker, 2012)).To achieve a regular time axis, the gaps in the Δ 13 CBaradla and Δ 18 OBaradla time series were filled by cubic-spline fitting, and resampled to an annual resolution by averaging (Fig. 1) (De Boor, 1978) in an R statistical environment (R Core Team, 2016) using the stringr package (Wickham, 2015) and the spline function of the stats package (R Core Team, 2016).The advantage of cubic spline fitting is that it is considered to be highly effective in preserving the spectral characteristics of the original record (Horowitz, 1974), and it outperforms linear interpolation -especially in the higher frequency domain (Schulz and Stattegger, 1997)-, and has a smaller chance of overfitting as a higher order spline.In the case of interpolation, regardless of the applied method, the high frequency components in a spectrum will be underestimated, thus the interpolated time series will become "reddened" (Schulz and Stattegger, 1997).Thus, the spectral bias caused by interpolation has to be objectively quantified and controlled.The spectral characteristics of the original and the interpolated data have to be compared to detect potential spectral artifacts.Thus, the threshold frequency was quantified beyond which (i) the interpolation (cubic spline in the present case) left the original spectrum mainly unchanged, and (ii) the significant powers of the original time series were in coherence with the interpolated spectrum.
To determine this threshold frequency objectively, first the spectral characteristics of the original-, and the interpolated time series were explored to see if the powers of the time series ' Lomb-Scargle Fourier Transform periodograms (LSP) (Lomb, 1976;Scargle, 1982) are significant against red-noise background from a first-order autoregressive (AR(1)) process using REDFIT (Schulz and Mudelsee, 2002).In the calculations, a Welch I window with two overlapping (50%) segments was used, the oversampling parameter was set as ofac=4 according to (Hocke and Kämpfer, 2009;Press et al., 1996), and hifac=1, the number of Monte-Carlo simulations to obtain the significance levels was nsim=1000 in line with studies dealing with speleothem time series (Holzkämper et al., 2004;Neff et al., 2001) boundary was estimated as the bias-corrected 80% chi-squared limit of a fitted AR(1) process.For the computations the redfit function of the dplR package was used (Bunn, 2008).The obtained periodograms will be referred to hereinafter as redfit Lomb-Scargle Fourier Transform periodograms (rLSP).The obtained rLSPs of the original and the interpolated time series were than visually compared.
After the visual comparison of the rLSPs, the coherence of the unevenly spaced original and the interpolated time series was also computed, using REDFIT-X (Björg Ólafsdóttir et al., 2016) developed for cross-spectral analysis of unevenly spaced paleoclimate time series.The run parameters were the same as indicated above.
Combining the visual comparison of the original and interpolated time series' rLSPs with their quantified coherence spectrum, the smallest significant period of the original data was determined, which was also present, and in coherence with the spline interpolated one.This could be set as the threshold frequency above which the original spectra could be taken as unbiased.Finally, the variance below this threshold frequency was removed/filtered from the spline interpolated time series using the bandpass function of the astrochron package (Meyers, 2014).

Wavelet transform coherence (WTC)
The data preprocessing to ensure a time series is evenly spaced in time (Section 2.2.1) was necessary to find the areas with common powers between the speleothem stable isotope time series and the climatic data in the time-frequency space.
Wavelet transform coherence (Torrence and Webster, 1999) was used to assess and visualize the coherence of the time series on so-called power spectrum density (PSD) graphs (e.g.Fig. 3).This method may be considered similar to a correlation coefficient, but with the difference that here we are dealing with a localized time-frequency space (Grinsted et al., 2004).It is based on wavelet spectrum analysis, a function localized in both frequency and time, with a mean of zero (Grinsted et al., 2004), and may be taken as the convolution of the data and the wavelet function (Kovács et al., 2010) for a time series (Χn, n=1,. . ., N) with a '' degree of uniform resolution Eq. ( 1): (1) Whereby N is the length of the time series, 'ψ*' is the wavelet function and 's' is the scale.In the study the Morlet mother wavelet (Morlet et al., 1982) was used to generate daughter wavelets, because it establishes a clear distinction between random fluctuations and periodic regions (Andreo et al., 2006) and has also been used in other speleothem studies (e.g.Ersek et al., 2012;Holmgren et al., 2003).
If two time series correlate, it does not mean that their WTC will indicate a strong common periodic behavior, because the periodic component has to be present in order to find a meaningful WTC.In the course of the evaluation only those positive signals were considered which were significant (α=0.1)against an AR(1) process, for details see (Roesch and Schmidbauer, 2014).Since, the wavelet functions are normalized to have unit energy, the obtained wavelet transforms may even be compared with other time series (Torrence and Compo, 1998).
From a practical point of view, the PSD graphs of the WTC analysis between the stable isotope time series and the monthly climate data were calculated to find the months with the highest response.These chosen months were averaged and their WTC with the proxy time series was calculated again.Thus, consecutive multi-monthly averages (seasons) were obtained, indicating the maximum response forming the final output of the analysis.The WTC PSD graphs were generated with the analyze.coherencyfunction of the Wavelet-comp package (Roesch and Schmidbauer, 2014) in R.

Data preprocessing before WTC analysis
After obtaining the rLSPs and the coherency spectra of the original-and spline interpolated Δ 13 CBaradla and Δ 18 OBaradla time series, these were compared as discussed in Section 2.2.1.
As low as the ~5yr period, the significant powers (α=0.8) of the original Δ 18 OBaradla record were all mirrored in the spline interpolated one with high coherency and the original-and the interpolated record indicated a similar pattern on their rLSPs (Fig. 2a).However, below the ~5yr period, the significant powers of the original Δ 18 OBaradla record were no longer reflected any more in the spline interpolated one, and the two time series' coherency became generally low, as well.
Therefore, to be consistent and conservative, the period domain below 4.5 yrs was omitted with a lowpass filter for the Δ 18 OBaradla time series.The same steps were then performed for the Δ 13 CBaradla record, as well, thus the domain below 7.5yrs was cut off from the Δ 13 CBaradla spectrum to avoid the spectral bias caused by spline interpolation (Fig. 2b).

Climate-composition relationship for the Baradla speleothems
To present the applicability of the now equally spaced paleoproxy data in regional climate models, the Δ 13 CBaradla and Δ 18 OBaradla records were compared to local (precipitation and temperature) and regional (NAO indices) climate variables.
The first step was to explore whether meaningful and significant coherences could be found in the time-frequency space between the equally spaced and 7.5yr lowpassed speleothem Δ 13 CBaradla time series and similarly lowpassed primary climate parameters (precipitation, temperature).In the case of the composite Δ 13 CBaradla record, the best coherence was found with the December-March average precipitation amount, and this was mostly in anti-phase.
The generally anti-phase coherence/relationship between precipitation amount and the composite carbon isotope record is in accordance with the general theory that a higher amount of precipitation results in enhanced biological soil activity, more biogenic CO2 in the soil, and consequently more negative Δ 13 C values in the speleothems of the Baradla Cave (Demény et al., 2017).Although the relationship between the Δ 13 C and precipitation amount had previously been observed for individual records using visual comparisons and regression analysis (Demény et al., 2017), it was found to be weak due to dating uncertainties and the varying effect of additional factors contributing to the precipitation-composition relationship (e.g.kinetic fractionation, vegetation change, prior calcite precipitation; see (Fairchild and Baker, 2012) for the compilation of governing factors).Nevertheless, with wavelet coherence analysis, it was not the whole spectra which was taken into account all at once (as in the case of linear regression analysis), but the relationship was mapped for each frequency band and observed over time.The phase of the coherence was observed to vary on the decadal-scale.A generally negative phase difference was revealed between Δ 13 C and precipitation amount (Fig. 3).Significant coherences were found between the Baradla speleothem's Δ 18 OBaradla records and the monthly-and multimonthly averages of the primary climate parameters with a common discontinuity around the 1960s (Fig. 4).Such coherences were to be expected, since the oxygen isotope composition of the speleothem is driven by temperature and drip water composition (Lachniet, 2009), with the latter directly related to meteoric water composition governed by atmospheric temperature and moisture origin (Dansgaard, 1964;Kaiser et al., 2001).
However, with regard to the phase differences between the Δ 18 OBaradla records and the climate parameters, a somewhat confusing picture can be observed in the investigated ~110 years.The phase differences changed multiple times, and a dominant direction could hardly be assigned.This may be a result of the complex interplay of governing factors.Eastern Alps (Kaiser et al., 2001), and also to the Baradla Cave region, in accordance with the anti-phase relationship observed in this study.This observation can be used to interpret further the isotope data of older speleothems in the light of the large scale atmospheric and oceanic variations of the North Atlantic realm.

Speleothem stable isotope records in a strongly NAO-influenced area
The methodology proposed and utilized on the Baradla speleothem was tested on an additional high resolution and accurately dated speleothem dataset.The Proserpine stalagmite's δ 13 C and δ 18 O records (Han-sur-Lesse Cave, Belgium; Van Rampelbergh et al., 2015) covering the period 1479-2000 AD were processed following the same procedure.
Although in the text only the section comparable with the previously discussed NAOi spanning 1857-2000 is processed and compared using WTC with the NAOi, the processed datasets for the entire time can be found in the supplementary material.After the necessary processing to achieve the even spacing of the records, the spectral assessment of the rLSPs of the original and the interpolated data indicated that a 3.5yr and 3.3yr lowpass filtering is necessary for the full δ 18 O and δ 13 C records of the Proserpine speleothem, respectively, to avoid bias caused by the interpolation.As a next step the season with strongest response was selected, by forming multi-monthly averages of the climate data, as in the case of the Baradla speleothem stable isotope records.Note here that the climate data was lowpassed, just as the speleothem data assessed together with it.
The relationship with the precipitation reconstruction was explored, and it proved to be mostly in-phase (where it gave any relationship at all, only from ~1960 onwards), most probably due to the superimposed effects of local factors and regional climate variations.With the NAOi, however, it indicated changing (Fig. 6a) and an anti-phase coherency (Fig. 6b) at the ~8yr period band.
The Proserpine δ 18 O record indicated a strengthened response with the November-March averages of the NAOst and the November-April averages of the NAOPC indices at the ~7yr period, predominantly in anti-phase (Fig. 7).The anti-phase relationship was most explicit in the case of the NAOPC index (Fig. 7b).The strength of the coherence was somewhere lower for the shorter comparison with the NAOPC than for the NAOst (Fig. 7) record, which had an extra ~50yrs overlap with the speleothem stable isotope record.
The reason why the NAO and the Proserpine δ 18O records are in anti-phase is because when the NAO index is positive, the cave receives more winter (low-δ 18 O) precipitation, while the negative NAO state induces cold and dry periods around the cave environment, resulting in increased δ 18 O values in the speleothem (Van Rampelbergh et al., 2015).The   The composite stable isotope records of the Baradla Cave speleothem were compared with monthly-resolved temperature and precipitation-amount data using wavelet transform coherence analyses.A generally negative relationship between the carbon isotope record and cold season precipitation amount was revealed, in accordance with earlier assumptions.In the case of the oxygen isotopic composition, its relationship with primary climate variables (temperature and precipitation amount) was less clear, probably due to the competing and/or superimposed factors that determine the carbonate oxygen isotopic composition.Nevertheless, the Δ 18 O record indicated an anti-phase relationship with North Atlantic Oscillation indices reflecting moisture origin.Specifically, in a negative NAO mode the moisture transport trajectory is shifted to the Mediterranean from where the Baradla Cave receives δ 18 O-enriched precipitation.These observations provide a firm base for the interpretation of stable isotope data obtained for the Baradla Cave system, NE Hungary.The now evenly resampled and lowpass filtered composite records can serve as regional benchmarks in future proxy paleoclimate evaluations.
Moreover, the methodology was successfully applied using published data for a stalagmite from Belgium (δ 13 C and δ 18 O records of the Proserpine stalagmite, Han-sur-Lesse Cave).Expressed relationships with the NAOi, especially with the δ 18 O records, were found in accordance with the earlier assumptions, lending credibility to the methodology applied and developed in this study.

Data availability:
The data sets produced with the methodology presented in the paper are available at doi: 10.1594/PANGAEA.875917.

Competing interests:
The authors declare that they have no conflict of interest.

Figure 1 :
Figure 1: Cubic spline fitted to the (a) original unevenly spaced Δ 18 OBaradla time series and the (b) annual averages derived from the interpolated time series.

Figure 2 :
Figure 2: rLSPs (upper panels) and coherency spectra (lower panels) of the original, and the spline interpolated (a) Δ 18 OBaradla and (b) Δ 13 CBaradla records of the Baradla speleothem.The bias-corrected 80% chi-squared limit of a fitted AR(1) process for the rLSPs is shown in grey.The coherency spectra were produced with a 95% Monte Carlo false-alarm levels (lower panels).The vertical black line indicates the determined cut-off period.

Figure 3 :
Figure 3: Time series of the 7.5yr lowpassed composite Δ 13 CBaradla records (on a reversed axis) with the gridded December-March average precipitation amount (Harris et al., 2014) (upper panel) and their time-frequency coherency image (lower panel).The white contours in the lower panel show the 90% confidence levels calculated based on 1000 AR (1) surrogates.The black arrows indicate the phase-angle difference.

Figure 4 :
Figure 4: Time series of the 4.5yr lowpassed composite Δ 18 O records with the gridded climate data (Harris et al., 2014) December-May average (a) temperature and (b) precipitation (upper panels) and their time-frequency coherency image (lower panels).For further details see the caption of Fig. 3.

Figure 5 :
Figure 5: Time series of the 4.5yr lowpassed composite Δ 18 O records with the December-May average (a) NAOst and (b) NAOPC reconstructions (upper panels) and their time-frequency coherency images (lower panels).For further details see the caption of Fig. 3.
anti-phase δ 18 O-NAOi relationship is especially significant for the PC-based NAOi.The detection of the relationship between the processed Proserpine speleothem record(Van Rampelbergh et al., 2015) and the NAOi further validates the procedure elaborated in this study and demonstrates the use of the established algorithm.Earth Syst.Sci.Data Discuss., https://doi.org/10.5194/essd-2017forjournal Earth Syst.Sci.Data Discussion started: 16 June 2017 c Author(s) 2017.CC BY 3.0 License.

Figure 6 :Fig
Figure 6: Time series of the 3.3yr lowpassed composite Δ 13 C records with the December-May average (a) NAOst

Figure 7 :
Figure 7: Time series of the 3.5yr lowpassed Proserpine Δ 18 O records with the (a) November-March average . Hatvani developed the model code and performed the simulations.Z. Kern preprocessed the data and provided the composite records.A. Demény conceived the study and provided the geochemical interpretation.Sz.Leél-Őssy was responsible for the cave studies.All authors took part in the manuscript preparation.