Speleothem stable isotope records for east-central Europe: resampling sedimentary proxy records to obtain evenly spaced time series with spectral guidance

Uneven spacing is a common feature of sedimentary paleoclimate records, in many cases causing difficulties in the application of classical statistical and time series methods. Although special statistical tools do exist to assess unevenly spaced data directly, the transformation of such data into a temporally equidistant time series which may then be examined using commonly employed statistical tools remains, however, an unachieved goal. The present paper, therefore, introduces an approach to obtain evenly spaced time series (using cubic spline fitting) from unevenly spaced speleothem records with the application of a spectral guidance to avoid the spectral bias caused by interpolation and retain the original spectral characteristics of the data. The methodology was applied to stable carbon and oxygen isotope records derived from two stalagmites from the Baradla Cave (NE Hungary) dating back to the late 18th century. To show the benefit of the equally spaced records to climate studies, their coherence with climate parameters is explored using wavelet transform coherence and discussed. The obtained equally spaced time series are available at https://doi.org/10.1594/PANGAEA.875917.

1 Introduction

With more than a hundred speleothem-related 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 they (i) are formed in relatively protected cave environments that render late-stage alteration less frequent; (ii) can be precisely dated by radiometric methods (e.g., U–Th dating); (iii) can be sampled at high spatial and temporal resolution, depending on the growth rate of the speleothem; (iv) are well distributed worldwide; and (v) provide multiple independent proxies (e.g., textural characteristics, trace element and stable isotope composition; for further details, see the comprehensive overview by Fairchild and Baker (2012) suitable to the reconstruction of past climate conditions.

Global, regional, and local processes collectively determine the final geochemical signature archived in a speleothem. However, some proxies are more likely to be influenced by large-scale factors (e.g., stable oxygen isotopes; Lachniet, 2009), others are clearly locally controlled (e.g., trace elements; Fairchild and Treble, 2009). Any of the previously mentioned proxies can be used in paleoclimate evaluation, assuming that a clear relationship between climate conditions and measured proxy values can be determined and described. To assist in the interpretation of the geochemical and stable isotope composition of speleothems, the most frequently used and widely accepted method is cave monitoring, in which the physical and chemical parameters of the cave environment and carbonate precipitation are determined in a multi-year study (e.g., Breitenbach et al., 2015; Mattey et al., 2008, 2016; Surić et al., 2018; Riechelmann et al., 2011). These can be used in forward modeling to quantify and understand the role of the dominant factors regulating the geochemical signal in speleothems (Baker and Bradley, 2010). The advantage of this approach is the possibility of gaining direct information on cave dynamics and speleothem formation in the course of environmental changes; the drawback is the usually short timescale, which rarely extends as far as a decade (Fairchild and Baker, 2012; Mattey et al., 2016). The appearance or absence of seasonality (i.e., the seasonality of cave air temperature, ventilation, the chemical composition of drip water) changes in the studied cave can be determined, but the effects of stronger environmental and climate changes may not be observed due to the short periods covered.

Another approach is the comparison of meteorological parameters and speleothem data that cover a sufficient calibration period (e.g., Demény et al., 2017; Jex et al., 2010). If compared with climate parameters, such data can serve as benchmark records, whose comparison may elucidate regional proxy correlations, or can be used to test complex climate models (Wackerbarth et al., 2012). The comparison can be made visually, but detailed statistical data processing (e.g., Baker et al., 2015; von Gunten et al., 2012; Wassenburg et al., 2016) can provide much more objective results.

Unfortunately, despite all efforts, uneven spacing is a common feature of sedimentary paleoclimate records. This characteristic usually biases the results of classical statistical methods and prohibits the application of time series analysis tools in many cases (Cosford et al., 2008; Mudelsee, 2010). There are excellent toolsets that can solve certain problems directly in the case of unevenly spaced data, e.g., determining whether it has a first-order autoregressive characteristic (Schulz and Mudelsee, 2002), estimating correlation with uncertainty limits between unevenly spaced autocorrelated series (Roberts et al., 2017), or conducting a variety of spectral analyses (Schulz et al., 1999; Schulz and Stattegger, 1997). However, a problem remains: the transformation of unevenly spaced data to a temporally equidistant time series is required in order for it to be suitable for the application of, for example, wavelet spectrum, multitaper analysis, and simple smoothers (Hammer, 2017). In the case of unevenly spaced sedimentary records, the required preprocessing step is most commonly performed using linear (e.g., Holmgren et al., 2003, Lachniet et al., 2004), or nonlinear interpolation techniques (e.g., Ersek et al., 2012), rescaling of the data (e.g., Deininger et al., 2017), or by insufficiently documented methods (e.g., McCabe-Glynn et al., 2013; Duan et al., 2014). A more complex approach combining spline smoothing and linear interpolation has been presented in a multi-paleoproxy study, where in addition the spectral characteristics of the data were assessed using multiple methods as well (Cosford et al., 2008).

However, any transformation which adds data or removes data from the original record unavoidably changes its spectral characteristics. This is a factor which can easily be overlooked even in advanced studies, although it deserves a greater degree of attention.

The present paper therefore aims to introduce an approach applied to selected stable carbon and oxygen isotope records from the Baradla Cave (NE Hungary for details see Sect. 2.1). The motivation was to obtain evenly spaced time series from unevenly spaced speleothem records with the application of a spectral guidance to avoid spectral bias and retain the genuine spectral characteristics, which can then be used in further proxy-climate evaluation studies. Section 2 describes the data sources and the proposed methodology for interpolation and resampling, combining the abilities of two available software packages (Björg Ólafsdóttir et al., 2016; Schulz and Mudelsee, 2002) and the development of the spectral control. Section 3 presents the evenly spaced data generated by the application of the methodology.

2 Materials and methods

2.1 Data sources

The stable carbon and oxygen isotope compositions (expressed as δ¹³C and δ¹⁸O values in ‰ relative to V-PDB) presented in this paper were derived from two stalagmites (VK1 and NU2; Fig. 1a) from the Baradla Cave (NE Hungary; 48^∘28^′ N, 20^∘30^′ E). Despite being 500 m apart in the cave, the two slow-growing speleothems show a clear similarity in growth rates (VK1: 0.37 and NU2: 0.32 mm yr⁻¹; Fig. 9 in Demény et al., 2017) based on lamina counting and pattern in stable isotope variability. The values, for VK1 and NU2, respectively (Demény et al., 2017), are as follows:

• δ¹³C: median, −10.47 and −9.97; range, 2.16 and 3.15; SD, 0.53 and 0.64.

• δ¹⁸O: median, −6.86 and −7.41; range, 1.69 and 2.39; SD, 0.4 and 0.43.

This is not surprising, since the rock cover above the two sample locations is comparable (100–150 m). Both deposits were seasonally laminated, providing a relative chronology of annual accuracy. Besides age estimation based on lamina counting, ¹⁴C activity was determined in the selected stalagmites. Radiocarbon measurements clearly reflected the post-1950 elevated atmospheric signal. Thus, the best match was obtained by shifting the NU2 and VK1 records by 3 and 5 years, respectively (Demény et al., 2017), supporting the synchronization of the individual isotope records. The good match of the isotopic peaks and the similar δ¹³C (Fig. 1a) and δ¹⁸O (Fig. 1b) amplitudes both argue for comparable water storage time in the feeding karst water system. Based on lamina counting, the chronological error at the base of VK1 is estimated to be between 3 and 5 years. However, the perfect matches found for δ¹³C (see Fig. 1a) and δ¹⁸O (see Fig. 1b) suggest that the chronological difference is much smaller. Moreover, the stable isotope records of the two stalagmites complement one another by filling the hiatuses found (Fig. 1a). To construct composite isotope records from the two stalagmites' δ¹³C and δ¹⁸O records, the original raw data were normalized using mean and standard deviation for the 1950–2000 period (an interval of simultaneous growth for the VK1 and NU2 stalagmites) and merged to a common timescale to provide a regional reference (Demény et al., 2017). Hereinafter, the composite records are expressed as Δ¹³C_Baradla and Δ¹⁸O_Baradla.

To test the skill of the methodology, the relationship of climate, and the processed evenly spaced Δ¹³C_Baradla and Δ¹⁸O_Baradla products was assessed. Thus, monthly averages of temperature and monthly precipitation totals corresponding to the cave location were retrieved a global gridded climate dataset (CRU TS3.23; Harris et al., 2014) with a time span covering 1901–2014 and a grid space of 0.5^∘.

2.2 Methodology

2.2.1 Interpolation and resampling with spectral guidance

To achieve a regular time axis, the gaps in the Δ¹³C_Baradla and Δ¹⁸O_Baradla 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 advantages of cubic spline fitting are that it is considered to be highly effective in the preservation of the spectral characteristics of the original record (Horowitz, 1974); it is robust and computationally efficient, thus being suitable to process large data sets (Musial et al., 2011); it outperforms linear interpolation – especially in the higher frequency domain (Schulz and Stattegger, 1997); and it has a smaller chance of overfitting as a higher-order spline.

Furthermore, spline interpolation does not introduce outliers as the Lomb–Scargle reconstruction (Hocke and Kämpfer, 2009) or the Kondrashov and Ghil technique (Kondrashov and Ghil, 2006) do when missing values are – for example – seasonally clustered, thus rendering the estimation of ranges of frequencies and the power spectrum unreliable (Musial et al., 2011). However, because of its “local” estimation characteristic (i.e., lack of predictive capacity), with the increase in a prolonged gap the accuracy of spline interpolation decreases. With the use of rigorous tests it has been documented that if the length of the gap is, for example, slightly above 10 % of the total data length, the mean absolute error of the estimation will be below 0.1. In the case of the Baradla composite stable isotope records, the maximum gap length was 4.85 years, corresponding to 2.24 %, meaning that the use of spline interpolation could be suggested in opposition to the previously mentioned methods. For further details on the benefits and drawbacks of the previously mentioned methods and the detailed documentation on the effects of gap length on mean absolute error the reader is referred to Musial et al. (2011).

It should also be noted that in the case of interpolation, regardless of the method applied, the high-frequency components in the spectrum will be underestimated; thus the interpolated time series will become “reddened” (Schulz and Stattegger, 1997).

Figure 1 Raw δ¹³C (a) and δ¹⁸O (b) records of the NU2 and VK1 speleothems (hiatuses indicated by gaps in the line), cubic spline fitted to the (c) original unevenly spaced Δ¹⁸O_Baradla time series and the (d) annual averages derived from the interpolated time series. The blue and orange rectangles in panels (a, b) mark the hiatuses of the NU2 and VK1 speleothems, respectively, as discussed in Demény et al. (2017).

Therefore, the spectral bias caused by interpolation has to be objectively quantified and checked. The spectral characteristics of the original and the interpolated data have to be compared to detect potential spectral artifacts. Thus, a threshold frequency is needed to be quantified, beyond which (i) the interpolation (cubic spline in the present case) leaves the original spectrum mainly unchanged, and (ii) the significant powers of the original time series are in coherence with the interpolated spectrum.

Figure 2 Illustration of the performance and robustness of the resampling protocol. (a) Three examples of the randomly resampled annual mean sunspot numbers (left panels), along with the histograms representing the distribution of time differences between the resampled values, and the rLSPs of the original and the corresponding resampled annual mean sunspot number data (right panels). (b) Histograms of the time differences between the Baradla stable isotope records' samples. (c) The rLSPs of the original and the 100 randomly resampled annual mean sunspot number data. The rLSP of the annual mean sunspot numbers is in yellow (right column in a and c). The rLSPs of the randomly resampled un-evenly spaced artificial (sunspot) time series, spline interpolated and annualized from the monthly mean sunspot numbers are in orange, purple and red for the three examples respectively (a), and in black for all the 100 randomly resampled time series (c). The bias-corrected 95 % chi-squared limit of a fitted AR(1) process for the rLSPs is shown in grey (a, c).

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 (LSPs) (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, as in the works of Hocke and Kämpfer (2009) and Press et al. (1996), and hifac = 1, the number of Monte Carlo simulations to obtain the significance levels was nsim = 1000, in line with other studies dealing with speleothem time series (Holzkämper et al., 2004; Neff et al., 2001), and the red-noise 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 (rLSPs). The obtained rLSPs of the original and the interpolated time series were then compared visually.

A record with pronounced periodic signals was selected to test the performance and robustness of the presented approach. Seasonal averages were computed from the monthly mean total sunspot numbers (SSN; WDC SILSO, 2017) and randomly resampled. Seven data points were taken randomly from each block of 10 years to replicate the resolution of the Baradla speleothem (avg. 0.7 stable isotope data per year). This was repeated a 100 times, so an ensemble of 100 randomly resampled time series were obtained, replicating the nonuniform resolution of the Baradla δ¹⁸O record. The maximum difference between the consecutive samples in the artificially resampled dataset was practically the same as in that in the Baradla speleothem record (∼ 7; Fig. 2a, b). Moreover, the relatively large gaps (≥ 4 year) represented about 5–10 % of the artificial time series, which was in harmony with the original speleothem (raw and composite) record (Fig. 2b). The only slight difference between the speleothem record and the artificial time series was seen in the most frequent gaps. In the case of the VK1 record, the most frequent gap between two samples was the same as in the artificial time series – 1 year – while in one of the raw speleothems and the composite record, this was 2 years.

As a final step in line with the proposed protocol, the 100 randomly resampled time series were spline-interpolated and annual averages were calculated, and these were then processed using REDFIT to assess their spectral characteristics and compare them to the spectra of the actual annual averages of the original (sunspot) record.

All the 100 rLSPs replicated the well-known ∼ 11-year sunspot cycle (α=0.95). Moreover, most of the rLSPs of the unequally spaced artificial time series mirrored an even smaller peak (∼ 8 years). However, the noise caused by spline interpolation can be clearly observed at the higher frequencies (Fig. 2c), highlighting the necessity of controlling/checking the spectral bias caused by interpolation.

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 the cross-spectral analysis of unevenly spaced paleoclimate time series. The run parameters were the same as those 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 then 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).

2.2.2 Wavelet transform coherence

Data preprocessing was necessary to ensure a time series is evenly spaced in time (Sect. 2.2.1) in order 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 (WTC) (Torrence and Webster, 1999) was used to assess and visualize the coherence and phase-angle difference of the time series on so-called power spectrum density (PSD) graphs (e.g., Fig. 4). The phase-angle differences are indicated by arrows. It should be noted that, depending on the chronological error, the phase information can be biased. For example, an error of 4 years could already reverse the phase angle of an 8-year period signal. Thus, keeping the chronological error of the input data (see Sect. 2.1) in mind, to avoid over-interpretation, only the dominantly in-phase relationships (indicated by arrows pointing to the right) and dominantly anti-phase relationships (arrows pointing to the left) are considered in the study.

WTC 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 (X_n, $n=\mathrm{1},\mathrm{\dots },N$) with a Δt degree of uniform resolution Eq. (1):

$$\begin{array}{}\text{(1)}& {\displaystyle }{W}_n^X\left(s\right)=\sqrt{{\displaystyle \frac{\mathrm{\Delta }t}{S}}}{\sum }_{n^{\prime }=\mathrm{1}}^N{X}_{n^{\prime }}{\mathrm{\Psi }}^{\ast }\left[\left(n{}^{\prime }-n\right){\displaystyle \frac{\mathrm{\Delta }t}{S}}\right],\end{array}$$

where N is the length of the time series, ψ^* is the wavelet function, and s is the scale. In this 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).

Figure 3 rLSPs (upper panels) and coherency spectra (lower panels) of the original, and the spline-interpolated (a) Δ¹⁸O_Baradla and (b) Δ¹³C_Baradla 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.

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 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. Then the monthly data which indicated the strongest coherence were averaged (e.g., December–January–February–March) and their WTC with the proxy time series was calculated again. Even if there was a month in between those indicating a strong relationship with the climate data (e.g., strong relationship with December and February–March; weak relationship with January), the multi-monthly averages (December–March and February–March) were calculated and assessed using WTC. 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.coherency function of the Wavelet-comp package (Roesch and Schmidbauer, 2014) in R.

Figure 4 Time series of the 7.5-year low-passed composite Δ¹³C_Baradla records (on a reversed axis) with the gridded December–March average (a) temperature and (b) 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, where an in-phase relationship is indicated by arrows pointing to the right, whereas anti-phase relationship is shown by arrows pointing to the left. The dotted area marks the period band below the cut-off threshold.

Figure 5 Time series of the 4.5-year low-passed composite Δ¹⁸O records with the gridded climate data (Harris et al., 2014) December–May average (a) precipitation and (b) temperature (upper panels) and their time–frequency coherency image (lower panels). For further details see the caption of Fig. 4.

3 Results and discussions

3.1 Data preprocessing before WTC analysis

After obtaining the rLSPs and the coherency spectra of the original and spline-interpolated Δ¹³C_Baradla and Δ¹⁸O_Baradla time series, these were compared as discussed in Sect. 2.2.1.

As low as a period of ∼ 5 years, the significant powers (α=0.8) of the original Δ¹⁸O_Baradla record were all mirrored in the spline-interpolated one, and to a high degree of coherency, and the original, and the interpolated record indicated a similar pattern on their rLSPs (Fig. 3a). However, around and below the 4.5-year period, the significant powers of the original Δ¹⁸O_Baradla record were no longer reflected in the spline-interpolated one, and the two time series' level of coherency became generally low as well. Therefore, to be consistent and to take a conservative (that is, careful) approach, the period domain below 4.5 years was omitted with a low-pass filter for the Δ¹⁸O_Baradla time series. It should be noted, however, that in cases when the storage component in the aquifer is larger than the numerically chosen cut-off period, the water storage time in the karst may amplify the auto-correlation of the isotopic time series.

The same steps were then performed on the Δ¹³C_Baradla record as well. Thus, around and below the ∼ 7.5-year period, the significant powers of the original Δ¹³C_Baradla record were no longer reflected in the spline-interpolated one, and the level of the two time series' coherency became generally low as well. Above the threshold, the avg. squared coherency is 0.9 and only drops below 0.8 in exceptional cases, while below it, its average is 0.57 and frequently drops to zero (Fig. 3b). Thus, the domain below 7.5 years was cut off from the Δ¹³C_Baradla spectrum to avoid the spectral bias caused by spline interpolation (Fig. 3b).

3.2 Climate signals in the stable isotope composition of carbonates of the Baradla speleothems

To present the applicability of the now equally spaced paleoproxy data in climate studies, the Δ¹³C_Baradla and Δ¹⁸O_Baradla records were compared to precipitation and temperature (primary 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.5-year low-passed speleothem Δ¹³C_Baradla time series and similarly low-passed primary climate parameters (precipitation, temperature; Fig. 4). In the case of the composite Δ¹³C_Baradla record, the best coherence was found with the December–March averages of temperature and precipitation amount, although a more persistent relationship was detected with precipitation than with temperature (Fig. 4).

The generally anti-phase coherence/relationship between precipitation amount and the Δ¹³C_Baradla record is in accordance with the former interpretation (Demény et al., 2017), namely that a higher amount of precipitation results in enhanced biological soil activity, more biogenic CO₂ in the soil, and consequently more negative Δ¹³C values in the speleothems of the Baradla Cave. Although the relationship between the Δ¹³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–sedimentary proxy relationship (e.g., kinetic fractionation, vegetation change, prior calcite precipitation; see Fairchild and Baker, 2012, for the compilation of governing factors). Nevertheless, with WTC 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 anti-phase difference was revealed between Δ¹³C and precipitation amount (Fig. 4a).

As for the relationship between Δ¹³C and temperature, significant coherences were observed in the 8–16-year period band, but these, however, weakened between ∼ 1960 and ∼ 1975. Furthermore, the phase-angle differences were not persistent afterwards (Fig. 4b). The varying coherence with temperature – regarding the period and phase angle (Fig. 4b) – might be related to biological activity in the soil driven primarily by the interplay between temperature and water saturation conditions mutually affecting the decomposition of the organic matter that is the dominant source of speleothem carbon (Noronha et al., 2015).

Significant coherences were found between the Baradla speleothem's Δ¹⁸O_Baradla records and the multi-monthly (December–May) averages of precipitation and temperature. The coherence between precipitation and Δ¹⁸O_Baradla was patchily developed, and was significant in the intervals 1900 to  ∼ 1915 and ∼ 1920 to  ∼ 1955 for the period range of 4–16 years, and between ∼ 1970 and 2010 for around the ∼ 5-year period (Fig. 5a). In the case of temperature and Δ¹⁸O_Baradla a more developed, less patchy pattern was seen, with one longer period lacking coherence between ∼ 1960 and the early 1980s (Fig. 5b). Coherence is primarily present in the 4–16-year period range, as in the previous case. A mutual characteristic of the relationship of both precipitation and temperature to Δ¹⁸O_Baradla is the common discontinuity around the 1960s (Fig. 5). Such coherences were to be expected, since the oxygen isotope composition of the speleothem is driven by temperature and drip water composition (Lachniet, 2009). Drip water composition is directly related to meteoric water composition, governed by atmospheric temperature and moisture origin in the mid-latitudes (Dansgaard, 1964), and has been documented in the surroundings of the study area (Bottyán et al., 2017; Kaiser et al., 2001). In relation to the seasonal characteristic of the response, a strong contrast exists in the seasonal water balance at the Baradla Cave. Precipitation exceeds evapotranspiration from October to April (Demény et al., 2017), which is well reflected in the isotopic character of the drip water, arguing for the overwhelming dominance (∼ 90%) of winter precipitation (Czuppon et al., 2018). Large-scale tropospheric circulation controls precipitation amount and water stable isotope composition over Europe (Comas-Bru et al., 2016), which invites future investigations of coherence with atmospheric circulation indices.

However, with regard to the phase differences between the Δ¹⁸O_Baradla records and the climate parameters, a somewhat confusing picture can be observed in the investigated ∼ 110 years. The phase-angle differences changed multiple times, and a dominant direction could hardly be assigned. This may be a result of the complex interplay of governing factors.

4 Data availability

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

5 Conclusions

A cubic-spline-based universally applicable methodology was developed to handle the uneven spacing of sedimentary proxy records. Additionally, the bias that interpolation may cause was evaluated with spectral control. The methodology was successfully applied to the composite stable carbon and oxygen isotope records of the speleothems of the Baradla Cave, NE Hungary. The composite stable isotope records of the Baradla Cave speleothem were compared with temperature and precipitation-amount data at a monthly resolution 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. With temperature, however, a less persistent relationship was observed, and this is suspected to be related to biological activity in the soil. In the case of oxygen isotopic composition, its relationship with temperature and precipitation amount was less clear, probably due to the competing and/or superimposed factors determining the carbonate's oxygen isotopic composition. These observations provide a firm basis for the interpretation of stable isotope data obtained for the Baradla Cave system. The now evenly resampled and low-pass-filtered composite records may also serve as regional benchmarks in future proxy paleoclimate evaluations.