A reference time scale for Site U1385 (Shackleton Site) on the SW Iberian Margin

Abstract We produced a composite depth scale and chronology for Site U1385 on the SW Iberian Margin. Using log(Ca/Ti) measured by core scanning XRF at 1-cm resolution in all holes, a composite section was constructed to 166.5 meter composite depth (mcd) that corrects for stretching and squeezing in each core. Oxygen isotopes of benthic foraminifera were correlated to a stacked δ18O reference signal (LR04) to produce an oxygen isotope stratigraphy and age model. Variations in sediment color contain very strong precession signals at Site U1385, and the amplitude modulation of these cycles provides a powerful tool for developing an orbitally-tuned age model. We tuned the U1385 record by correlating peaks in L* to the local summer insolation maxima at 37°N. The benthic δ18O record of Site U1385, when placed on the tuned age model, generally agrees with other time scales within their respective chronologic uncertainties. The age model is transferred to down-core data to produce a continuous time series of log(Ca/Ti) that reflect relative changes of biogenic carbonate and detrital sediment. Biogenic carbonate increases during interglacial and interstadial climate states and decreases during glacial and stadial periods. Much of the variance in the log(Ca/Ti) is explained by a linear combination of orbital frequencies (precession, tilt and eccentricity), whereas the residual signal reflects suborbital climate variability. The strong correlation between suborbital log(Ca/Ti) variability and Greenland temperature over the last glacial cycle at Site U1385 suggests that this signal can be used as a proxy for millennial-scale climate variability over the past 1.5 Ma. Millennial climate variability, as expressed by log(Ca/Ti) at Site U1385, was a persistent feature of glacial climates over the past 1.5 Ma, including glacial periods of the early Pleistocene (‘41-kyr world’) when boundary conditions differed significantly from those of the late Pleistocene (‘100-kyr world’). Suborbital variability was suppressed during interglacial stages and enhanced during glacial periods, especially when benthic δ18O surpassed ~ 3.3–3.5‰. Each glacial inception was marked by appearance of strong millennial variability and each deglaciation was preceded by a terminal stadial event. Suborbital variability may be a symptomatic feature of glacial climate or, alternatively, may play a more active role in the inception and/or termination of glacial cycles.


Introduction
The Iberian Margin is well known as a source of rapidly accumulating sediment that contains a high-fidelity record of millennial-scale climate variability for the late Pleistocene Shackleton et al., 2004;Sanchez Goñi et al., 2002;Tzedakis et al., 2004;Vautravers and Shackleton, 2006;Skinner & Elderfield, 2007;Martrat et al., 2007;Margari et al., 2010Margari et al., , 2014Hodell et al., 2013a). A unique aspect of the record is the potential to correlate millennial-scale climate events from the marine environment to polar ice cores and terrestrial records. Few places exist in the world where such detailed and unambiguous marine-iceterrestrial linkages are possible. The continuity, high sedimentation rates, and fidelity of the sediments on the SW Iberian Margin make this region a prime target for ocean drilling. In 2011, an Ancillary Program Letter (APL 763) was submitted to the Integrated Ocean Drilling Program (IODP) requesting four days of ship time to drill multiple holes at a single site to 150 m. The program was approved and scheduled as part of the IODP Expedition 339 (Mediterranean Outflow).
During IODP Expedition 339, four holes were drilled at Site U1385 ( Fig. 1) with each reaching~150 m below sea floor (Hodell et al., 2013b). We refer to Site U1385 as the "Shackleton site" in honor of Nick Shackleton's pioneering work on SW Iberian Margin cores, which has played a pivotal role in our understanding of millennial-scale climate variability over the last glacial cycle Shackleton et al., 2004). The prime objective of IODP Site U1385 was to extend the piston core records from the SW Iberian Margin and provide a continuous marine record of Pleistocene millennial climate variability that can be correlated to polar ice cores and European terrestrial sequences.
A necessary prerequisite for realizing the full potential of Site U1385 is a continuous composite section and an accurate age model. Without a robust stratigraphical framework, paleoenvironmental observations and interpretations remain ambiguous and limit our ability to infer mechanisms of past climate change. Here we provide a revised composite depth scale and chronology for Site U1385, which will serve as the reference stratigraphic and chronologic framework for future paleoceanographic studies of the region.
We scanned cores in all holes of Site U1385 at 1-cm resolution using two Avaatech XRF core scanners in Texel (Netherlands) and one at Cambridge (U.K.) to acquire relative changes in sediment composition. We used log(Ca/Ti) to correlate among holes and revised the shipboard composite spliced section, including correction for distortion within individual cores by stretching and squeezing of cored intervals to match the composite splice. All holes were aligned to a composite reference splice, thereby permitting the integration of data among all holes.
Because of the need to correlate Site U1385 with marine sediment, ice core and speleothem records, we present four age models that can be applied in accordance with the chronostratigraphic needs of particular researchers and projects: (1) "Age_Depth_Radiocarbon" applies to the upper 5 rmcd of Site U1385 only and is based on correlating the log(Ca/Ti) of Site U1385 to the carbonate record of core MD99-2334, which has a detailed radiocarbon chronology . (2) "Age_Depth_GreenSyn" covers the last 400 kyr and is derived by correlating the log(Ca/Ti) record from Site U1385 to the log(Ca/Ti) record of nearby Cores MD01-2444 and -2443 (Hodell et al., 2013a). The age model is a modified version of the one derived by correlating millennial-scale cold events in planktonic δ 18 O and alkenone sea surface temperature (SST) to the synthetic Greenland temperature record using the Speleo-Age time scale (Barker et al., 2011).
2. Site U1385 location and core recovery IODP Site U1385 (37°34.285′ N, 10°7.562′ W, 2578 m below sea level) was drilled in December 2011 during Expedition 339. The site is located on a spur on the continental slope of the SW Iberian margin (Promonotorio dos Prinicipes de Avis), which is elevated above the abyssal plain and influence of turbidites (Fig. 1). It is very near the  (Hodell et al., 2013a). Five holes (A-E) were cored at Site U1385 using the advanced piston corer (APC) system to a maximum depth of~155.9 mbsf. A total of 67 cores were recovered from Holes A (17 cores), B (16 cores), C (1 core), D (16 cores), and E (17 cores), totaling 621.8 m of sediment with a nominal recovery of 103.2% (N 100% due to post-recovery core expansion). The offset distance among holes is about 20 m (Expedition 339 Scientists, 2013). The sediment lithology consists of uniform, nannofossil muds and clays, with varying proportions of biogenic carbonate and terrigenous sediment (Expedition 339 Scientists, 2013).

XRF core scanning
The split archive halves of Site U1385 were analyzed using three Avaatech X-ray fluorescence (XRF) core scanners to obtain semi-quantitative elemental data at 1-cm spatial resolution. Cores from Holes U1385A and U1385B were measured at the Royal Netherlands Institute for Sea Research (NIOZ) and Holes U1385D and U1385E at the Godwin Laboratory in Cambridge. The core surface was carefully scraped cleaned and covered with a 4-μm thin SPEXCertiPrep Ultralene foil to avoid contamination and reduce desiccation (Richter et al., 2006). All XRF scanners irradiated a surface of 10-mm high by 12-mm wide every 1 cm with identical instrument settings used at both NIOZ and Cambridge (Table 1). Two additional replicates were measured every 25 measurements to provide confidence limits on the elemental data obtained. Element intensities were obtained by post-processing of the XRF spectra using the Canberra WinAxil software with standard software settings and spectrum-fit models. The element intensities depend on the element concentration but also on matrix effects, physical properties, the sample geometry, and hardware settings of the scanner (Tjallingii et al., 2007). However, log-ratios of element intensities are simple linear functions of log-ratios of concentrations, minimize sample geometry effects, and allow calculation of statistically meaningful confidence limits (Weltje and Tjallingii, 2008). Therefore, log-ratios provide an objective and convenient way to compare different XRF records measured on different instruments in terms of relative chemical variations. The Ca/Ti data are archived on the PANGAEA and NOAA Paleoclimate data systems and the complete XRF data set will be presented in a companion paper.

Stable isotopes
Samples were wet sieved at 63 μm and dried in an oven at 50°C. Stable isotopes were measured on the planktic foraminifera Globigerina bulloides selected from the 250 to 355 um size fraction and the benthic foraminifer Cibicidoides wuellerstorfi from the N212 μm fraction. Foraminifer tests were crushed and soaked in a solution of 1% hydrogen peroxide for 30 min in individual vials. Acetone was added and the samples placed in an ultra-sonic bath for 10 s, after which the liquid was carefully decanted to remove any contaminants. The samples were dried in an oven at 50°C overnight. Isotopic analyses of the samples were performed using a VG SIRA mass spectrometer with a Multicarb system for samples with a mass exceeding 80 μg. Analytical precision is estimated to be ±0.08‰ for both δ 18 O and δ 13 C. For smaller samples (b 80 μg), measurements were performed on a Thermo Finnigan MAT253 mass spectrometer fitted with a Kiel device. Analytical precision is estimated to be ±0.08‰ for δ 18 O and ±0.06‰ for δ 13 C, respectively. Results are reported relative to V-PDB. All isotope measurements were made in the Godwin Laboratory, University of Cambridge and are archived on the PANGAEA and NOAA Paleoclimate data systems.

Color scanning
Two independent sets of sediment color data were measured. During Expedition 339, color reflectance was measured on the archive halves of sediment cores using an Ocean Optics USB4000 spectrophotometer mounted on the automated "Section Half Multisensor Logger" (Expedition 339 Scientists, Methods, 2013). Freshly split cores were covered with clear plastic wrap and each measurement consisted of sequential 2-nm wide spectral bands from 400 to 900 nm.
Color data were also collected using the Avaatech line-scan camera system. Split cores surfaces were prepared by fine scraping with a stainless steel or glass blade. The cores were placed in the Avaatech XRF scanner and the shade over the window closed to prevent interference from external light. Standard image color cards were placed alongside the core for calibration. The white fluorescent lights were allowed to warm up and stabilize for 30 min before beginning each camera scan. All cores were scanned with the camera aperture set to 8+ with an exposure time of 5 msecs. The monochromatic filters in the camera split the light into RGB channels at 630 nm, 535 nm and 450 nm with good separation among channels. The RGB data were converted to the CIEl-a*-b* color system for comparison with the shipboard measurements. The data were subsequently averaged with a running mean of 10 data points for subsequent data compression. Color data are archived on the PANGAEA and NOAA Paleoclimate data systems.

Construction of the composite section
Sediment missing at core breaks is a common problem of core recovery using the hydraulic piston core (Ruddiman et al., 1987;Hagelberg et al., 1995). Constructing a composite section by splicing stratigraphic intervals from multiple holes has become a routine procedure for the IODP and is performed onboard ship (Expedition 339 Scientists, Methods, Stratigraphic Correlation, 2013). Each core is shifted by a constant offset using a single tie point, which provides a crude linear offset that does not account for non-linear drilling distortion. For example, this method does not correct for stretching (expansion) or squeezing (contraction) within a particular core and, consequently, cores are only aligned at the single splice tie point. Other intervals may not be perfectly correlated and the same sedimentary feature may have a different composite depth in each hole, causing misalignments when comparing data from different holes.
Here we used the log(Ca/Ti) measured at 1-cm intervals to correlate among holes and define a composite spliced section (Supplementary Tables 1 & 2; Fig. 2 and Supplementary Fig. 1), containing no gaps or disturbed intervals, to 166.5 m composite depth (mcd). We follow the terminology of Pälike et al. (2005) who suggested the use of revised meter composite depth (rmcd) to refer to revision of the shipboard mcd, and corrected revised meter composite depth (crmcd) for the correction that accounts for distortion in each core. We re-evaluated the shipboard splice and established new core offsets (Supplementary Table 1). To the extent possible, we attempted to maintain the same tie points as defined by the shipboard composite. A spliced stratigraphic section was constructed by selecting the best and most representative sections from overlapping holes ( Fig. 2 and Supplementary Fig. 1). Some minor deformation features (microfaults, contorted strata) were noted in some holes (see Table 5 of Expedition 339 Scientists, 2013), but often these did not extend to neighboring holes. Disturbed sections were avoided when constructing the revised composite splice.
Because of the large sampling demand for Site U1385, it was necessary to sample all holes; thus, we correlated all holes to the spliced composite section. The crmcd was defined by aligning features in the log(Ca/ Ti) record in each hole to the revised spliced composite section using Analyseries (Paillard et al., 1996). The differences between rmcd and crmcd reflect the relative degree of stretching and squeezing that was required for each core. These differences vary from −0.86 to +1.55 m ( Fig. 2 and Supplementary Fig. 3). This method corrects for distortion within individual cores, permitting the integration of data from all holes. Because all holes are aligned to the splice, the same crmcd should correspond to the same stratigraphic horizon in all holes. In practice, the accuracy of the alignment among holes is dependent upon the scale of the correlative features and variability of the log(Ca/Ti) record. One measure of this accuracy is the standard deviation of the stacked log(Ca/Ti) records averaged among the four holes ( Fig. 2 and Supplementary Fig. 2). We estimate that log(Ca/Ti) features are correlated to the decimeter level or better.

Age models
A robust time scale is essential for Site U1385 because the "Shackleton site" is likely to become an important reference section for marine-iceterrestrial core correlations and the study of orbital-and millennialscale climate variability. For this reason, we present four age models for Site U1385 that can be applied according to the requirements of particular researchers and projects. A comparison of these age models also provides an indication of the variability of age estimates associated with differing age modeling techniques. We closely monitored changes in sedimentation rate for each age model when defining age-depth points. If substantial changes in sedimentation rates were implied by an age model, we evaluated whether the implied changes in the flux of biogenic and/or detrital sediment were plausible and justified.

Radiocarbon-based age model
Kasten core MD99-2334K, which is near Site U1385 (Fig. 1), has a detailed age model based on high-resolution radiocarbon dating . We correlated the log(Ca/Ti) in the upper 5 rmcd of Site U1385 to the weight %CaCO 3 record of MD99-2334K and transferred the radiocarbon age model to U1385 (Fig. 3). Linear sedimentation rates for the last 25 kyr average~20 cm kyr −1 (Fig. 4A; Supplementary Table 3).

Synthetic Greenland
The log(Ca/Ti) record from the upper 56 rmcd of Site U1385 was correlated to the log(Ca/Ti) record from nearby piston cores MD01-2444 and -2443 using 55 tie points ( Fig. 5; Supplementary Table 4) (Hodell et al., 2013a). MD01-2444 is essentially at the same location as Site U1385 whereas MD01-2443 is~26 km to the north (Fig. 1). The age model is a modified version of the one derived by Barker et al. (2011) based on correlation of millennial-scale cold events expressed in planktic δ 18 O and SST (Tzedakis et al., 2004;Vautravers and Shackleton, 2006;Martrat et al., 2007;Margari et al., 2010) to the synthetic Greenland temperature record.
The log(Ca/Ti) records of MD01-2444 and U1385 are very similar for the last 190 ka, but there are some differences between U1385 and MD01-2443 (Fig. 5). The correlation is ambiguous during MIS 10 between 350 and 360 kyr and during early MIS 11. Between 400 and 420 ka, the log(Ca/Ti) at Site U1385 is flat whereas MD01-2443 shows considerably more variability including a peak at~408 kyr that is not seen in the U1385 record. It is uncertain whether these differences are the result of problems with MD01-2443, U1385, or both.
Interval sedimentation rates vary between~5 and 30 cm kyr −1 , increasing towards the top of the core owing to low compaction of sediment (Fig. 4B).

Oxygen isotopes
The benthic δ 18 O record of Site U1385 was correlated to the oxygen isotope stack of Lisiecki and Raymo (2005) (Fig. 6; Supplementary Table 5). The correlation between U1385 and LR04 demonstrates that Site U1385 contains a continuous record from the Holocene tõ 1.45 Ma (Marine Isotope Stage 47). Linear sedimentation rates average 10.9 cm kyr −1 for the entire section with interval sedimentation rates varying between~1 and 25 cm kyr −1 (Fig. 4C). The lowest sedimentation rate is confined to the interval between 490 and 430 ka (Termination V and early MIS 11) when a condensed section or brief hiatus is suspected.

Orbital tuning
Astronomical tuning offers the potential of a more accurate time scale than can be achieved by solely correlating benthic δ 18 O to the LR04 stack if more precise orbitally-paced signals are available than those used to construct the LR04 record. Variations in sediment color of SW Iberian Margin sediments contain a strongly modulated precession signal over the past 400 kyr (Hodell et al., 2013a). The modulation of precession by eccentricity constitutes a powerful tool for both developing the orbitally-tuned age model and for testing its robustness . Hodell et al. (2013a) suggested sediment redness (a*) at this site is an almost pure precession signal with no apparent phase lag. However, we could not use a* for Site U1385 because the signal is greatly attenuated downcore because of the reductive dissolution of Fe minerals in the zone of sulfate reduction (Expedition 339 Scientists, 2013). Instead, we use the lightness (L*) of the sediment for tuning.
Because sedimentation rates are fairly constant throughout Site U1385, we initially conducted time series analysis in the depth domain. The L* power spectrum of the detrended data reveals cycles at~2.5 and 2 m that correspond to precession cycles of 23 and 19 kyr, respectively, assuming a mean sedimentation rate of~11 cm kyr −1 . The L* signal was bandpass filtered at a central frequency of~0.45 +/− 0.1 cycles m −1 to capture both of the precession peaks (Fig. 7).
The filtered signal shows an amplitude modulation in the depth domain that resembles the eccentricity modulation of precession, especially for the interval below 100 rmcd (Fig. 7). Whereas amplitude modulation has been criticized as a test of orbital tuning in the time domain (Huybers and Aharonson, 2010), the modulation in depth cannot be an artifact of the tuning process.
It is a relatively straightforward process to bring the L* signal into line with precession. However, the lag between precession forcing and the response of sediment lightness is unknown. To estimate this lag, we rely on the radiocarbon dating of the CaCO 3 record in kasten core MD99-2334K . The precession minimum occurs at 11.3 ka and the peak in carbonate and lightness occurs at 8.5 kyr (Fig. 3). This 2.8-kyr lag is similar to the estimated~3-kyr lag between the insolation maximum and sapropel deposition in the Mediterranean and attributed to a delayed increase in monsoon intensity caused by precession-paced North Atlantic cold events (Ziegler et al., 2010). At 37°N, the mean summer insolation maximum also lags the precession index by about 2.8 kyr and is roughly in phase with L* (Fig. 3). Thus, we correlated peaks in the L* signal to local summer insolation at 37°N (Fig. 8; Supplementary Table 6). Peaks in L* are also well correlated with lows in Ti/Al and dark color at Sites 967/968, which mark the occurrence of sapropels associated with precession minima (Konijnendijk et al., 2014). Shackleton et al. (1995) suggested several criteria for assessing the robustness of a tuned time scale, including: (i.) it should not conflict with other geochronologic information; (ii.) it should not imply changes in sedimentation rate that are implausible; (iii.) it should concentrate variance at the Milankovitch frequencies; and (iv.) the coherence between the proxy used and astronomical target should be high and,  (Lisiecki and Raymo, 2005); and (D) astronomically tuned age model derived by correlating L* to local summer insolation at 37°N. more importantly, the amplitude modulation of the tuned record should be closely aligned with the astronomical forcing.

Evaluating the age models
The tuned time scale of Site U1385 meets all of these criteria: (i.) There is generally good agreement among the four independently determined age models and with the revised positions of polarity reversal boundaries ( Fig. 9) (C. Xuan, pers. comm.). The U1385 age model is also consistent with other astronomicallytuned time scales (see Section 5.1.1). (ii.) Sedimentation rates over the entire section average 9.2 cm kyr −1 and vary from a low of 1 cm kyr −1 during Termination V to a high of 33 cm kyr − 1 . The lowest rate is attributed to a condensed section or short hiatus associated with Termination V and early MIS 11. Variability in sedimentation rates is greatest with the Greenland synthetic age model because of a greater number of age control points over a shorter period of time (Sadler, 1981). (iii.) The tuned age model, which is based on correlating L* to precession, results in the transfer of a large amount of variance to the obliquity band for both L* and log(Ca/Ti).
(iv.) The coherence is high at all orbital frequencies but, more importantly, the eccentricity modulation of the precession component is evident in both depth and time domains. Shackleton et al. (1995) suggested that complex demodulation is superior to coherence for evaluating the reliability of a time scale. Fig. 10 shows the result of complex demodulation of the orbitally tuned L* signal at a central frequency of 0.048 cycles/kyr (21-kyr period) relative to the Earth's eccentricity cycle. The amplitude modulation of the L* is clearly aligned with orbital eccentricity. This amplitude modulation of precession by eccentricity is also evident in the depth domain (Fig. 7), which is independent of orbital tuning.

Comparison to other time scales
The prediction, and subsequent verification, that the accepted ages of recent polarity reversal boundaries were too young by 6-7% was a triumph for validating the astronomical tuning method and refinement of the geologic time scale (Shackleton et al., 1990). The reason ODP Site   Lisiecki and Raymo (2005). Red crosses represent the tie points. 677 was well suited for astronomical tuning was the combined effects of obliquity in benthic δ 18 O and a strong precession cycle (modulated by eccentricity) in the planktonic δ 18 O record. Shackleton et al. (1990) used the output of the ice volume model of Imbrie and Imbrie (1980) as the tuning target for the planktonic and benthic δ 18 O records of Site 677. The forcing was boreal summer insolation at 65°N with the assumption of a constant lag of ice sheet response to orbital forcing (T M = 15 kyr) and the ratio of the response time for ice growth and decay (b = 0.6). Shackleton et al. (1990) used the same constants for the entire record acknowledging that these parameters likely changed as the size of ice sheets increased during the Quaternary. Although there have been minor modifications to the time scale of Shackleton et al. (1990), confirmation of the Site 677 age model has been provided by many subsequent studies (Hilgen, 1991;Bassinot et al., 1994;Lisiecki and Raymo, 2005).
Benthic δ 18 O on the independently tuned time scale of Site U1385 agrees well with the Site 677 benthic δ 18 O record (Fig. 11A). There is a slight lead of the Site U1385 benthic δ 18 O record over the Pacific especially at glacial terminations after 650 ka. This phasing is consistent with a demonstrated lead of Atlantic δ 18 O over the Pacific by~3000-5000 years for terminations I and II (Skinner and Shackleton, 2005).
The most commonly used Quaternary time scale today is LR04 (Lisiecki and Raymo, 2005). Construction of this age model involved stacking benthic δ 18 O records and correlating the stack to an ice volume model driven by insolation with the assumption of a fixed lag between insolation forcing and ice sheet response at particular periodicities for the last 1.5 Ma. However, benthic δ 18 O is not an unambiguous proxy for ice volume (Shackleton, 2000;Elderfield et al., 2010Elderfield et al., , 2012, and the temperature and ice volume components of δ 18 O have different lag times with respect to insolation forcing (Hodell et al., 2013a).
In general, the tuned time scale of Site U1385 compares favorably with LR04 within the estimated error of the chronology (Fig. 11B; Lisiecki and Raymo, 2005), which is ±4 kyr for the past million years and ± 6 kyr for the interval from 1.0 to 1.5 Ma. The tuning of Site U1385 disagrees with LR04 near the base of the record between MIS 44 and 47 where the ages of the stage boundaries are in closer agreement with those of Site 677 (Fig. 11B). Our positioning of some isotope stage boundaries also differs from Lisiecki and Raymo (2005) for the  interval from 1 to 1.5 Ma (Fig. 14F & G). This difference results partly because we emphasized precession in the tuning process, whereas Lisiecki and Raymo (2005) emphasized obliquity and because we have used millennial variability, in addition to benthic δ 18 O, to recognize the end of interglacial stages or substages (Tzedakis et al., 2012a,b).
The tuned time scale of Site U1385 agrees well with the EDC3 and AICC2012 chronologies for the EPICA-DomeC Antarctic ice core for the past 800 ka (Fig. 11C). EDC3 and AICC2012 were derived by inversion of a snow accumulation and mechanical ice flow model using a set of independently dated horizons (e.g., age markers) along the core (Parrenin et al., 2007;Bazin et al., 2013). The differences between EDC3 and AICC2012 are small and within the uncertainty of the EDC3 age model, which is estimated to be ± 3 kyr at 100 ka and increasing to ± 6 kyr below 130 ka to the base of the core.
The strong precession cycle in L* at Site U1385 offers the opportunity to compare the chronostratigraphy to the Mediterranean cyclostratigraphy derived by tuning sapropel layers to precession (Hilgen, 1991;Lourens et al., 1996;Lourens, 2004). A recent revision to the sapropel chronology was proposed by Konijnendijk et al. (2014) who correlated lows in Ti/Al at Sites 967 and 968 to peaks in boreal summer insolation at 65°N assuming a 2.7-kyr time lag. A comparison of the independently tuned L* record of Site U1385 and the Ti/Al cyclostratigraphy reveals reasonably good agreement (Fig. 12), especially for the interval below 800 kyr. The correlation can be improved further by fine-tuning the peaks in L* with lows in Ti/Al at Sites 967/968 (Supplementary  Table 7), which directly synchronizes the Site U1385 with Mediterranean cyclostratigraphy. Hodell et al. (2013a) demonstrated the relative proportion of biogenic carbonate and detrital input to the sediment, as reflected by log(Ca/Ti) in core MD01-2444, is highly correlated with carbonate content. On orbital time scales, the primary control on CaCO 3 content is variations in clay flux rather than changes in carbonate productivity (Thomson et al., 1999). The log(Ca/Ti) record shows strong cyclicity at each of the Milankovitch frequencies with the residual reflecting suborbital variability (Fig. 13).

Millennial variability over the last 1.5 Ma
Millennial-scale variations in log(Ca/Ti) on the SW Iberian Margin are correlated with planktonic δ 18 O and alkenone SST, resembling the Greenland ice core δ 18 O record for the last glacial period (Hodell et al., 2013a). Cold stadial events are marked by decreases in log(Ca/Ti) whereas warmer interstadials coincide with increases in log(Ca/Ti). Assuming this relationship held for older glacial periods, log(Ca/Ti) offers a means to assess changes in millennial variability on the Iberian Margin during glacial periods for the past 1.5 Ma.
The log(Ca/Ti) of the upper 5 m of Site U1385 shows the classic deglaciation structure observed in many SW Iberian Margin cores with low values during H2, H1, and the Younger Dryas and relatively high values during the last glacial maximum (LGM), Bolling-Allerod (BA), and early Holocene (8 kyr) (Fig. 3). Each of the older Heinrich events is also marked by minima in log(Ca/Ti), whereas prominent Greenland insterstadials are represented by peaks in log(Ca/Ti) (Fig. 14A). The log(Ca/Ti) records suggests that millennial-scale variability was a persistent feature of the SW Iberian Margin record during each of the glacial periods of the last 1.5 Ma, including those of the early Pleistocene. Glacial periods during the 41-kyr world displayed significant millennialscale variability despite the fact that glacial boundary conditions differed from the latest Pleistocene. This observation is consistent with Raymo et al. (1998) who also demonstrated the persistence of millennial-scale in the early Pleistocene on the basis of proxies of iceberg discharge and deep-water chemistry at ODP Site 983. Contrary to other studies (Weirauch et al., 2008), we do not observe a significant change in millennial variability of log(Ca/Ti) across the Middle Pleistocene transition at Site U1385.
Suborbital climate variability was enhanced during glacial periods and suppressed during full interglacial periods throughout the Pleistocene. McManus et al. (1999) suggested that millennial-scale variability in ice-rafted debris increases when benthic δ 18 O surpasses a critical threshold of 3.5‰ (unadjusted δ 18 O for C. wuellerstorfi). Although δ 18 O is affected by both temperature and ice volume, the threshold presumably reflects the time when ice sheets became large enough to reach the coast and interact with the ocean. Raymo et al. (1998) found a slightly lower threshold value of 3.3‰ for the middle Pleistocene (MIS 40-44), whereas Bailey et al. (2010) proposed an even lower value (3‰) for the late Pliocene and earliest Pleistocene. Our log(Ca/Ti) and δ 18 O results are generally consistent with a threshold for increased millennial variability of~3.3 to 3.5‰ (Fig. 14A-G). McManus et al. (1999) further suggested that millennial variability was suppressed when δ 18 O surpassed a second threshold of 4.5 to 4.6‰ in benthic δ 18 O, creating a "window" in which millennial variability is enhanced between 3.5 and 4.5‰. The 4.5‰ threshold is not exceeded in glacial periods beyond 650 ka (MIS 16).
Each of the glacial terminations is marked by a terminal stadial event that has been labeled x.1 where x is the number of the preceding glacial marine isotope stage (Fig. 14A-G). It's unclear whether these terminal stadial events are actively involved in, or simply a consequence of, the deglacial process ( 10. Amplitude of the 21-kyr signal of the L* record at Site U1385 (red line), as determined by complex demodulation, compared with the signal of orbital eccentricity (Laskar et al., 1993;black line). For precession we demodulated the L* record at a central frequency of 0.048 cycles/kyr and used a low-pass Tukey filter with 100 weights and a 50% cutoff at 0.015 cycles/kyr. Fig. 11. Comparison of the tuned age model of Site U1385 (black) with other timescales including Site 677 (Shackleton et al., 1990), LR04 (Lisiecki and Raymo, 2005), and EDC3 (Parrenin et al., 2007).

2011
). For example, strong terminal millennial events may lead to CO 2 degassing in the Southern Ocean through a bipolar seesaw mechanism , thereby hastening deglaciation.
Millennial variability in log(Ca/Ti) is particularly strong during the "early glacial" period associated with glacial inception. Raymo et al. (1998) made a similar observation for the early glacial intervals of MIS 40 and 44 in the middle Pleistocene. The climate system appears to be particularly sensitive to millennial variability at the onset of glaciations. This may reflect an enhanced sensitivity of overturning circulation to fresh water forcing when AMOC is still strong following an interglacial period. Additionally, the location of iceberg discharge relative to the source areas of deep-water formation may be important for determining the sensitivity of the climate system to freshwater forcing. Shackleton (2006) highlighted the importance of stratigraphy and time scale development for all aspects of paleoceanography. We have produced a robust composite section and time scale for Site U1385 that will underpin future paleoceanographic studies. The age models back to 1.5 Ma are in good agreement with each other and the tuned age model compares favorably with existing Fig. 13. The log(Ca/Ti) record (black) at Site U1385 with bandpass filters for eccentricity (brown), obliquity (green), and precession (blue). The sum (orange) of the orbital components explains a significant amount of the variance in the log(Ca/Ti) signal. time scales (e.g., 677, LR04, EDC3, sapropel cyclostratigraphy). The tuned age model is applied down-core to produce a continuous time series of log(Ca/Ti) variations. This record offers a rapid means of assessing past changes in millennial-scale climate variability on the Iberian Margin for the last 1.5 Ma.

Conclusion
Millennial-scale variability is an intrinsic characteristic of glacial climate states of the Pleistocene. Suborbital climate variability at Site U1385 was enhanced during glacial periods when benthic δ 18 O exceeded 3.3-3.5‰ (McManus et al., 1999), and strongly suppressed during full interglacial periods. No two glacial cycles were exactly alike with respect to either the magnitude or pacing of millennial climate variability during the past 1.5 Ma. Each glacial onset was marked by a strong increase in millennial variability and, in contrast, each deglaciation was marked by a 'terminal stadial event' followed by an interglacial period when millennial variability was suppressed. An important outstanding problem in paleoclimatology is to understand how millennial-and orbital-scale climate variability interact to shape the observed patterns of Pleistocene climate change. In this regard, Site U1385 is an important reference section for reconstructing and studying orbital and suborbital climate variability over the past 1.5 Ma.  (Railsback et al., 2015). We have marked the start of the interglacials at the terminal stadial event indicated by salmon shading. The end of the interglacial stages or substages is placed where benthic δ 18 O decreases and millennial variability begins in the Ca/Ti record. The position of stage boundaries along top axis are those identified by Lisiecki and Raymo (2005). Dashed lines on the benthic δ 18 O record correspond to values of 3, 3.5, 4, and 4.5‰ to aid recognition of δ 18 O thresholds in millennial variability as expressed in Ca/Ti.