North Atlantic forcing of moisture delivery to Europe throughout the Holocene

Century-to-millennial scale fluctuations in precipitation and temperature are an established feature of European Holocene climates. Changes in moisture delivery are driven by complex interactions between ocean moisture sources and atmospheric circulation modes, making it difficult to resolve the drivers behind millennial scale variability in European precipitation. Here, we present two overlapping decadal resolution speleothem oxygen isotope (δ18O) records from a cave on the Atlantic coastline of northern Iberia, covering the period 12.1–0 ka. Speleothem δ18O reveals nine quasi-cyclical events of relatively wet-to-dry climatic conditions during the Holocene. Dynamic Harmonic Regression modelling indicates that changes in precipitation occurred with a ~1500 year frequency during the late Holocene and at a shorter length during the early Holocene. The timing of these cycles coincides with changes in North Atlantic Ocean conditions, indicating a connectivity between ocean conditions and Holocene moisture delivery. Early Holocene climate is potentially dominated by freshwater outburst events, whilst ~1500 year cycles in the late Holocene are more likely driven by changes internal to the ocean system. This is the first continental record of its type that clearly demonstrates millennial scale connectivity between the pulse of the ocean and precipitation over Europe through the entirety of the Holocene.

and X skimmer cone combination was used during U analysis. The high sensitivity Jet sample and X skimmer cone combination was used for Th analysis as the samples had relatively small amounts of radiogenic 230 Th available for analysis. Sample delivery to the plasma was via an ESI PFA nebulizer tip (50 µl/min uptake rate) and a Cetac Aridus II desolvating nebulizer. Samples were introduced into the nebulizer in 0.2 M HCl -0.05 M HF and the system cleaned between samples using 0.25 M HCl -0.1 M HF following Potter et al. (2005). Instrument sensitivity was c. 500-600 V/ppm for U and c. 1100-1200 V/ppm for Th at the above sample uptake rates and with the specified cone combinations. Ar with trace N 2 was used as the transfer gas, with N 2 employed to minimize U and Th oxide production in the plasma. A standard-sample-standard bracketing protocol was employed for data collection. U was analysed in static multicollection mode with 234 U measured in the axial SEM and other U isotopes measured in Faraday cups. CRM 112a was used to determine the SEM/Faraday bias, and CRM 112a spiked with the IRMM 3636 233 U-236 U tracer was used to determine the instrumental mass bias during analytical sessions (primary correction based on 233 U/ 236 U and checked using the spike-corrected 238 U/ 235 U). Th was analysed in static multicollection mode with 230 Th measured in the axial SEM and 229 Th and 232 Th measured on Faraday cups. Mass bias and SEM/Faraday gain was monitored during Th analytical sessions by analysing an in-house reference solution prepared from Ames Th (mainly 232 Th) mixed with 229 Th and 230 Th high purity tracers, previously calibrated against CRM 112a U. All raw data were corrected for hydride (~1-2 ppm) and down mass peak tailing (~2-3 ppm at 1 AMU, < 0.4 ppm at 2 AMU).
Data were reduced using an in-house spreadsheet, where ages and their uncertainties were calculated using an uncertainty propagation protocol based on McLean et al. (2011). Correction for initial Th contributions to measured activity ratios and calculation of initial [ 234 U/ 238 U] employed U series functions in Isoplot version 4 Excel 2010 add-in (see Ludwig, 2003b). Activity ratios and ages were calculated using the decay constants of (Cheng et al., 2013). A U/Th composition of [ 238 U/ 230 Th] = 0.8 ± 0.4 was used to correct all data for initial ("detrital") Th.
Speleothem age models (Figures S1 and S2) were generated from the detritus-corrected U-Th ages (Table S1) using the StalAge algorithm (Scholz and Hoffmann, 2011). StalAge identifies and removes U-Th date outliers, followed by a Monte Carlo simulation to create the age model based upon the U-Th analyses and their associated uncertainties (see detailed discussion in Scholz and Hoffmann (2011). Table S1: Detritus corrected U-Th ages for speleothems ASR and ASM. Data pertaining to the sampling locations and relevant U-Th measurement data from the analysis of speleothem's ASR and ASM.  Figure S1: Speleothem ASR chronology (green line). Developed from 10 U/Th analyses on speleothem ASR using StalAge modelling; red lines denote age model error. U series dates and associated errors are presented (black diamonds) to show sample spacing. Speleothem is shown at base of figure, as are the approximate locations of U-Th analyses (black arrows).

Speleothem growth phases
Speleothem ASR grew between 12,500 and 500 BP with most rapid deposition occurring during the Younger Dryas and early Holocene (1 data point every 2-3 years). Carbonate deposition was interrupted at 8,600 ±400 BP for approximately four thousand years until ~4,000 BP. This hiatus was characterised by intermittent speleothem growth and changes to speleothem crystal structure and colour. δ 18 O data from this period is removed from further analysis due to potential isotope fractionation during carbonate deposition. Speleothem deposition then continued during the late Holocene until 500 BP. Within this data, episodes of drying are observed at 11,800, 10,600, 9200, 8600, 6400, 4400, 2700, 1500 and 230 BP (Fig. 2).
Speleothem ASM offered a shorter record growing between 7,850 and 0 BP. Importantly this speleothem grew relatively rapidly during the mid Holocene, almost completely covering the period of growth hiatus in ASR; a period of 650 ±300 years exists between the end of initial growth in ASR and the start of growth in ASM. One possible break in growth for ASM is identified by the StalAge algorithm between 3,800 and 2,150 BP, no obvious change in speleothem structure or colour is identified for this period. However, to ensure consistency, data related to this hiatus is also excluded from the final data set.

Stable isotope analysis
The speleothem δ 18 O profile consists of 1244 isotope measurements from ASR and 880 from ASM, sampled using an automated micro-mill with 0.3 mm diamond encrusted drill bit along the central portion of the growth axis, creating between 50-100 μg of carbonate powder per sample. Stable isotope analysis was then undertaken at the NERC Stable Isotope Facility, British Geological Survey, using an IsoPrime isotope ratio mass spectrometer with Multiprep device; average 2σ uncertainty is 0.07‰. Isotope values are reported relative to the international VPDB standard.

Estimation of cycle periodicity and time varying amplitudes
Cycle periodicity and amplitudes have been estimated using a Dynamic Harmonic Regression (DHR) model with cycle length selection through maximisation of variance explained by single frequency DHR analysis.
Dynamic Harmonic Regression is a special case of Unobserved Components Model often used in econometrics (see e.g. Harvey and Proietti (2005)) and in this application it takes the form of: Subscript t denotes temporal variability of time series or its parameters and in a time series is usually the sample number. The additive components of the model, which are not observed individually, form in this case a plausible decomposition of the series in the frequency domain, and thus in Equation 1 is a low frequency trend component, is what is usually termed irregular component, including noncyclic elements, is the observation disturbance normally assumed to be a Gaussian distributed serially uncorrelated series, and is the harmonic cyclical component.
In Equation 2 is the number of frequencies , = 1, … , used in the model, , and , , = 1, … , are stochastic time varying parameters (TVP's). Estimation of these time-varying parameters is fully described in Young et al., (1999) and implemented in the CAPTAIN Toolbox for Matlab™. Implementation of a multi-frequency model is through simple addition of cyclic components, and is numerically robust due to the orthogonality of the harmonic components as long as there is a sufficient spectral separation between them, and when they are present in the data. Such multifrequency DHR has been used in a different application in Tych et al., (2002).
One of the main stages of solving this class of problems is selection of the covariance parameters of the Kalman Filter and Fixed Interval Smoother. Their definitions and some approaches are provided and reviewed in Young et al., (1999) and Harvey and Proietti (2005). In the present context these parameters determine the time scales of variation of the estimated trend levels and amplitudes of the cyclic components. These need to be selected in such a way that the trend component does not cover the spectrum of the cycles (that would lead to multiple colinearity in the model and associated numerical and interpretation problems.) The cyclic components' amplitudes need to be allowed to vary but to make any interpretation feasible their variability needs to be slower than the cycle they define.
It has to be noted here that the aim of tuning the Unobserved Component Models is not to achieve a perfect fit to the data. While this is achievable through choosing the variance parameters, the estimated components would then lose their interpretation, for instance the trend becoming fast varying and following the data closely. The aim in this case is to achieve objective and interpretable estimates of trend and of cyclic components, while keeping the fit to the data at a 'sensible' level.
The other parameters of the model include the cycle periodicities. These too need to be estimated objectively based on the available data. In the present analysis both periodicities are estimated using the following data based optimisation process. The set of , = 1, … , in a single cycle analysis consists of harmonics of the hypothetical frequency 0 = 2 , so that = 0 . This leads to the method applied in our analysis of finding the 'best' periodicity as the one that produces the highest power of the cyclic component signal. We use a simple quadratic norm of the cycle component: where ( 0 ) is the cyclic component (2) with the proviso that the model fit expressed as R 2 does not drop significantly.
is a commonly used proportion of the variance of observations explained by the DHR model output ̂.
In other words, we seek the cycle frequency 0 which gives the best proportion of the data explained by the cyclic component of the DHR model with fixed variance parameters. Such frequency can be identified by a simple univariate search using objective function (3) within Matlab ™ environment. Because of the orthogonal nature of the harmonic functions such searches can be performed individually, given their sufficient spectral separation, as orthogonality excludes interaction between the components.
For the presented data the periodicity search generates two well-defined maxima of Equation (3) at 1290 and 1490 years. These values were then used in a dual frequency DHR cycle analysis.
There is a strong indication of two separate periodicities present in the signal, as shown in Fig.S3 where the periodic component norm as in Equation 3 above (effectively: variance contribution of the periodic component) shows strong and well defined peaks at both periodicities, with a well defined trough in between them. The spectral resolution of this type of analysis is robust thanks to the orthogonality of the harmonic functional base of the DHR method, so we are confident in this result.
In terms of sensitivity of the DHR estimates to the NVR (variance parameters) the observed relationships between the periodicities of the two cycles persist throughout a broad range of the variance parameters spanning several orders of magnitude, so this dual-peak result is highly unlikely to be an artefact. The full, dual frequency DHR model explains 73.8% of the data variance (Fig. S4), showing that it is not only providing a detection tool but also indicating that the two cyclic components dominate the data. The remaining series of model residuals does not show any significant structure at longer time-scales, indicating that the DHR model is effective in explaining the entire structure of the data at the time scale longer than 150 years (10 samples). Evaluation of these hypotheses would require additional data records not available at this time.
(a) (b) Figure S4: Combined and detrended speleothem data with the modelled harmonic components (thick black line).
Both Matlab ™ scripts used in this analysis and the CAPTAIN Toolbox for Matlab™ are available upon request from one the contributing authors (Tych), who is also one of the principal developers of the CAPTAIN Toolbox.

Record comparisons
Data pre-processing for comparisons of our data set with other climate records involved putting the compared series into a common time-base, which was regularly sampled to allow the use of time series analysis methods, including cross-correlation.
Nudge procedure: The records have been transformed from the original sampling time base into regular 15-year sampled time series. In this nudge procedure the series was first interpolated at the required 15-year grid (dt=15), then the interpolated samples which did not have a 'real' data point within dt/2 from the interpolated sample time were replaced by missing value flags (NaN in Matlab and other systems), so that the procedure did not create artefacts in the form of 'data' where there had been none originally.
Subsampling procedure: The uniformly sampled series could, if this was required, be subsampled (decimated every k-samples) to the required lower resolution by first applying an anti-aliasing lowpass no-lag filter (IRW: see e.g. Young et al, 1999, Young 2011) followed by taking every k-th sample of the smoothed (filtered) series. This procedure was used to bring the compared time series to the same temporal resolution and to remove higher frequency components in the frequency range above that of interest (smoothing). This allowed the comparison of like with like in terms of the temporal patterns of the series.
Detrending: Trend is understood here as a slow and smooth change in the series covering the lowest part of the spectrum, below the periodicities or patterns of interest (IRW trend: see e.g. Young et al, 1999, Young 2011. This amounts to low-pass filtering of the series with no static or dynamic lags involved.

Standardisation:
The detrended series are then simply scaled to have standard deviation of one. This is a standard procedure in time series and generally data analysis, used when different time series need to be presented 'scale-free' to compare their patterns, including cycles.

Cross-correlation analysis:
This standard time series analysis method requires the time series to be provided on the same time-base. In the classical definition, standard correlation between the samples of the time series is calculated for lags 0 (simultaneous), +-1, +-2, … so providing the level of linear relationship between the series shifted by 0, +-1 etc. (Box et al 2011).

Uncertainty estimate of Asiul record:
We did not assume any initial error, the vertical uncertainty band in any comparison at reduced resolution is a conservative estimate of the 95% confidence interval resulting from smoothing the combined and detrended Asiul record required for the resolution change with no aliasing. The horizontal uncertainty is derived from the StalAge model outputs based on the individual speleothems U/Th chronologies.
The cross-analysis of series and timing errors included: 1. Comparative analysis and visualisation of both Cueva de Asiul records with timing errors. Analysing the consistency of the two Cueva de Asiul records (ASM, ASR) was based on the overlap period (550-2100AD) where both records are available ( Fig S5); both Cueva de Asiul ASR and ASM series were 'nudged' into the 15-year uniform sampling and decimated with anti-aliasing; cross-correlation function was estimated with a peak close to zero lag (coincidental) of 0.2, which was significant. This correlation was calculated for the estimated values. Bearing in mind the timing uncertainty of both of these values ( Fig S5 shows the overlap period with timing uncertainty bands) the two series show similarity, justifying the combining of these records as has been done in the main analysis.

Figure S5
: Combined and detrended speleothem data shown during the growth overlap period. Speleothem time series (thick black line) with U/Th error envelopes (grey bands) between 550 and 2100AD, 0 lag correlation of 0.2 is statistically significant. Visual analysis between the records indicates that a stronger relationship maybe derived if the lower resolution ASR record was allowed to vary within its U/Th error envelope.
2. Comparative analysis of Bond et al., (2001) and Cueva de Asiul records. Bond et al., data was nudged and resampled into a 60 year sampling period to allow direct comparisons with the Cueva de Asiul data subsampled with anti-aliasing into 60 year samples, then standardised as above. During the first 4000 years of the records the correlation between the two series reaches 0.505 with standard error of 0.1250, and is therefore significant. At a 15 year subsampling this relationship is 0.482 with standard error 0.129. This is visible in Fig. 3 of the main text, which shows the Bond record alongside the Cueva de Asiul data set with U/Th error plotted. If the series from present up to the hiatus in Cueva de Asiul data at 7500 is taken, then the correlation is weaker, reaching 0.165 with standard error of 0.09, so not significant. For the Early Holocene, there is no significant correlation when using the same statistical constraints as above. This is to be expected given the high timing errors at this age. Repositioning of the oldest section of the Cueva de Asiul data set (9000 -12000BP, t= -390 years; within the boundaries of the StalAge errors) shows a high level of significant co-variation between the data sets (0.560, s.e = 0.106; Fig. S6), indicating statistically the visual similarities between the records (Fig. 3 main text) once temporal error is taken into consideration. 3. Comparative analysis of Thornalley et al., (2009) and Cueva de Asiul records Thornalley et al., (2009) data was nudged and resampled into 60 year sampling period, which was an average time between the samples near the beginning of the record (most recent); isolated gaps resulting from the nudging process were interpolated over using cubic splines, there were no larger gaps, where this process could introduce artefacts. The same process of detrending with variance intervention over the Cueva de Asiul hiatus period was applied as for the combined Cueva de Asiul record. This was due to the apparent step change in the Thornalley et al., (2009) data within the period of the Cueva de Asiul hiatus. The data was then standardised as above and a cross-correlation function was calculated for the two series for the entire period up to the Cueva de Asiul hiatus (0-7800BP). Correlation analysis after this period was not undertaken due to different processes acting on the ocean density record in the early Holocene (Thornalley et al., 2009). For the standardised Cueva de Asiul record the cross correlation peaked at a lag of -4 (Thornalley ocean series leading) with a crosscorrelation coefficient of 0.264 (s.e. 0.089). It is also apparent that there is a strong higher frequency component in the Cueva de Asiul data, not showing in the Thornalley series. To compare like with like, the Cueva de Asiul data was smoothed to cover a similar spectrum as that of Thornalley et al., (2009)