Low-latitude climate variability in the Heinrich frequency band of the Late Cretaceous Greenhouse world

Introduction Conclusions References


Introduction
The origin of sub-Milankovitch climate variability with periods between 5 and 15 kyr (e.g., Hagelberg et al., 1994) remains elusive and is far less understood than cycles directly related to orbital climate forcing. So-called Heinrich events of the last glacial cycle from the North Atlantic (Heinrich, 1988;Bond et al., 1992;Hemming, 2004) represent the best known and most widely studied case of sub-Milankovitch variability, although variability on this timescale has been described from older glacials of the last 2.5 Myr as well (e.g., Becker et al., 2006 and references therein). Heinrich events denote rapid and massive discharges of icebergs into the North Atlantic during the last glacial cycle, suggesting that ice sheets are either directly responsible for or amplify an initial climate signal. They were first described by Heinrich (1988), who linked these to low-latitude climate to explain their recurrence time of ∼ 7-8 kyr. On the contrary, MacAyeal (1993) held internal ice-sheet oscillations responsible for their formation, while other autogenic models favor amplified jökulhlaupt events (Johnson and Lauritzen, 1995) or ice shelf instability (Hulbe, 1997; see also Hemming, 2004). They were related to white-noise stochastic forcing with a magnitude similar to random changes in insolation of ≥ 0.5 Wm −2 by Hyde and Crowley (2002).
However, it cannot be excluded that such internal or stochastic events are phase locked with initial external climate changes that act on the same timescale. Indeed, other studies involved an external climatic forcing to explain instability of fringing ice shelves (Hulbe et al., 2004) linked to precession paced variations in El Niño-Southern Oscillation (ENSO) intensity in the Pacific (Clement et al., 1999). Following Heinrich (1988), a low-latitude origin of Heinrich events was favored by McIntyre and Molfino (1996), who studied late Pleistocene abundance variations in the coccolith species Florisphaera profunda in the equatorial Atlantic and related these to precession-driven changes in zonal wind-driven divergence in the sub-Milankovitch frequency band at times of minimum eccentricity. Ziegler (2009) and Turner (2004) relate Heinrich events to more frequent El Niño events associated with precession-induced climate variability that originates between the tropics (to explain the sub-Milankovitch periods) and is exported to higher latitudes. Moreover, several modeling studies predict sub-Milankovitch variability through non-linear climate response (e.g., Braun et al., 2005, see also Hagelberg et al., 1994 and references therein), and numerous records (e.g., Steenbrink et al., 2003;Elrick and Hinnov, 1995;Shackleton, 1998, 1999;Zhao et al., 2006) indeed show such a pattern. Thus sub-Milankovitch variability has in particular been found in glacials of the last 2.5 Myr (Heinrich, 1988), although its origin may lie within the tropics (e.g., Hagelberg et al., 1994), suggesting that ice sheets amplify an initial (sub)orbital climate signal. Rutherford and D'Hondt (2000) suggest that ∼ 1.5 myr ago semi-precession cycles propagated from tropical to higher latitudes, and argue for a causal relation between semi-precession cycles, precession and eccentricity.
In case Heinrich events are related to low-latitude climate change and paced by precessional forcing and nonlinear feedback in the sub-Milankovitch band, cycles with the same period may be expected from pre-Quaternary times. Such cycles have indeed been described from Pliocene lignite-bearing lacustrine successions of northern Greece (Steenbrink et al., 2003) and platform carbonates of the Great Bahama Bank (Reuning et al., 2006). Sub-Milankovitch cycles have further been described in Upper Jurassic marls (Rodríguez-Tovar and Pardo-Igúzquiza, 2003;Boulila et al., 2010), Triassic marine sections (Zühlke et al., 2003;Wu et al., 2012) and Devonian carbonates (De Vleeschouwer et al., 2012). They have further been found in fluvial floodplain successions in the Eocene Willwood Formation in North America (Abdul Aziz et al., 2008), although the latter may equally well be associated with autogenic processes acting on the floodplain (Bown and Kraus, 1992).
In the present case, a logical choice to look for the persistence of sub-Milankovitch cycles is in marine successions of greenhouse periods of Earth history, such as the Cretaceous and Jurassic, when ice sheets were absent or played a minor role, and ice-driven responses are not expected. Evidence for the presence of such climate variability has been found in, for example, the Late Jurassic (Rodríguez-Tovar and Pardo-Igúzquiza, 2003;Boulila et al., 2010), while a prime example comes from the marine Cretaceous (Campanian) of the South Atlantic (Park et al., 1993). The latter report semi-precession cycles based on time series analysis of high-resolution magnetic susceptibility records from DSDP (Deep Sea Drilling Project) Site 516F located at a paleolatitude of ∼ 30 • S. However, using MS as a proxy for carbonate content does not exclude the possibility that the precession-related cycles represent sedimentary quadruplets. In these quadruplets, the two carbonate minima, or MS maxima per precessionrelated cycle, have a different climatic origin. Up to now, such quadruplets have only been described from the Mediterranean (e.g., de Visser et al., 1989;Hilgen et al., 2003), and the early South Atlantic was located at similar -but southern Hemisphere -paleolatitudes. In addition, visual inspection of the MS records of Park et al. (1993) reveals three rather than the expected two MS maxima per precession-related cycle, casting doubt on the inferred semi-precession nature of the cycles. Here, XRF elemental analysis is used to test the potential quadruplet structure of the precession-related cycles. Time-series analysis is applied to examine whether these variations represent true semi-precession cycles or have a shorter period similar to that of Heinrich events.

Geological setting and stratigraphic framework
DSDP Site 516 was drilled during DSDP Leg 72 at 30 • 16.59 S latitude and 35 • 17.10 W longitude, in a water depth of more than 1300 m on the Rio Grande rise in the South Atlantic . This site was the shallowest site drilled during DSDP Leg 72 and was designated for the study of upper water column characteristics of the southwestern Atlantic during the Neogene-Quaternary and Late Cretaceous; sediments deposited in deeper water were recovered at other sites. Hydraulic piston coring provided the base for detailed studies of astronomically induced climate forcing during the Quaternary and Miocene, while the Eocene-Oligocene and Cretaceous-Paleogene boundaries were targeted in rotary drilled carbonate-rich successions.
Part of the Upper Cretaceous at Site 516F was subjected to a detailed cyclostratigraphic study by Herbert and d'Hondt (1990), which revealed the influence of precession and eccentricity. To study sub-Milankovitch climate variability, Park et al. (1993) selected the interval between 1145 and 1166 m b.s.f. (meters below seafloor). This interval was placed by Weiss (1983) in the Globotruncana arca and G. ventricosa zones, based on the successive first occurrences of the nominate species. He positioned the Campanian base at the G. arca FO, but this boundary is now provisionally placed at the bottom of C33r in GTS2012 (Ogg and Hinnov, 2012). In the meantime the planktonic foraminiferal biozonal scheme has been significantly modified (Petrizzo et al., 2011;Ogg and Hinnov, 2012), largely because of the marked diachronous nature of events such as the G. ventricosa FO (Petrizzo et al., 2011). The magnetostratigraphy of the Upper Cretaceous at Site 516F was studied by Hamilton and Suzyumov (1983). Based on the common identification of the Cretaceous-Paleogene boundary, they showed that the magnetic polarity sequence identified at Site 516F could be correlated to the magnetostratigraphy of the Gubbio section of Alvarez et al. (1977). This correlation shows that the succession ranges from Chron 34 to 29 and that the interval studied by Park et al. (1993) falls within C33r and belongs to the early Campanian.
At that time, Site 516 was located at a paleolatitude of ∼ 30 • S, and the South Atlantic Ocean was a smaller, much more confined ocean basin than at present. However, Ndisotope studies indicate that a major change in deep ocean circulation occurred during the Campanian, attributed to the opening of a pathway through the Southern Ocean (Robinson et al., 2010;Martin et al., 2012;Robinson and Vance, 2012).
Though it is argued that the Rio Grande Rise formed a restriction for deep water flow (Frank and Arthur, 1999), Ndisotopes suggest that a pathway must have existed from the mid-Campanian onwards (Robinson et al., 2010;Robinson and Vance, 2012).
Since the exact date of the opening of the central Atlantic gateway (CAG) is not well constrained, the possibility of deep water circulation in the Southern Atlantic during the early Campanian cannot be excluded. However, because Nd-isotope studies estimate the opening of the CAG to be in the middle to late Campanian (Robinson et al., 2010;Robinson and Vance, 2012), it seems reasonable to assume that deep water circulation during the studied time interval was minimal.

Material and methods
A Konica Minolta CM-600d spectrophotometer was used to measure color reflectance with a resolution of ∼ 1.5 m on core surfaces of DSDP Site 516F at the BCR (Bremen Core Repository) at MARUM (Bremen). This was done for the interval between 1145 and 1166 m b.s.f. studied by Park et al. (1993). In addition, samples of ∼ 1.5 cm length (∼ 12 cm 3 , quarter cores) were taken at the same spacing from selected shorter intervals within this long interval (Sections 5 (81-150 cm) and 6 (0-93 cm) of core 113 and Section 1 (0-135 cm) of core 114). This resulted in a set of 198 core samples. No continuous splice is available for DSDP Site 516. A composite of cores 113 and 114 was created by connecting these cores directly. A new depth scale (meters composite depth, or mcd) was made for the composite, which was 53 cm longer compared to the record by Park et al. (1993) due to the inclusion of core section 113-7.
Samples were weighted using a Mettler P1210. A Niton XLi XRF (X-ray fluorescence) analyzer at Utrecht University was used to measure elemental abundances of the samples. Sawed down-core surfaces of samples were used for XRF measurements. The MS of the samples was measured using a KLY-2 Kappabridge. Finally, color reflectance of the selected shorter interval was measured on the core samples and core surface directly at the BCR. The interval for color measurements is longer than that for other measurement techniques, because some of the color data were generated directly from the core surfaces where no samples were taken. Due to cracks in the cores no equal spacing was possible for the color measurements. All the unprocessed results of measurements described above can be found in a Pangaea online database at http://doi.pangaea.de/10.1594/PANGAEA.828372.
Analyseries version 1.1 (Paillard et al., 1996) is used for Blackman-Tukey spectral analysis and Gaussian bandpass filtering. In addition, point-moving averages (pma's) of the L * , a * , MS and XRF records were calculated to reduce noise; this results in a smoothed data set. Pma's are chosen so that the cycles discussed are clearly revealed. A stronger (less strong) smoothing is applied to illustrate longer (shorter) period cyclicity.
Wavelet analysis including a significance test was applied using the software provided by C. Torrence and G. Compo. This software is available at URL: http://paos.colorado.edu/ research/wavelets/, and the reader is referred to Torrence and Compo (1998) for a description of the method. Prior to wavelet analysis, data were detrended and normalized using a Morlet mother wavelet. Power spectra were generated using the Redfit program (Schulz and Mudelsee, 2002, version 3.8e) based on the (bias-corrected) Lomb-Scargle Fourier transform; a Welch window is used. Results for the confidence levels are based on 1000 simulations.

Results
Results will first be presented for the color records of the long interval between 1145 and 1166 mcd (Fig. 1) and then for the color, MS and XRF records of the selected shorter intervals (Fig. 2).

Long-color record
The lightness (L * ) and redness (a * ) records of the long interval are shown in Fig. 1, including pma's and bandpass filtered components. The L * record and its 10 pma (Fig. 1a) reveal marked shifts to low values, with minima located approximately 50 cm apart. These well-defined minima occur in bundles of 3-5, with approximately 2.5 m between bundles. Moreover, three intervals (1145-1147, 1152-1156 and 1159-1163 mcd) can be recognized in which the averaged L * values reach most extreme minima. These intervals can also easily be recognized when a 30 pma is applied. Higher frequency variations can be observed throughout the L * record as well, with one or two additional, usually less distinct, minima in between the well-defined minima of the ∼ 50 cm cycle.
The a * record follows a different pattern with less marked (bundling of) peaks, and a trend towards higher values before returning to lower values in the top 4-5 m (Fig. 1b). The interval from 1145 to 1150 mcd reveals three distinct shifts to lower values. These follow a sawtooth pattern, marked by an abrupt change to high values followed by a gradual return to normal values interrupted by higher frequency variations (see 10 pma, Fig. 1b). This is succeeded by an interval marked by minima that are ∼ 90 cm apart (∼ 1149-1153 mcd), and an interval (up to 1156 mcd) in which the maxima are separated by ∼ 50 cm and, thus, more closely spaced (see Fig. 1b, 10 pma's of a * ). The uppermost interval is marked by a return to negative values and shows a less regular pattern and less pronounced cyclic variations.
The L * Blackman-Tukey power spectrum (Fig. 3) reveals a marked double peak (of ∼ 43 and 53 cm) that corresponds to the ∼ 50 cm cycle, which reflects the prominent shifts to minimum values observed in the long record. The bandpass filtered component of this cycle, which encompasses both spectral peaks, is shown above the L * record in Fig. 1a. This filtered component reveals marked amplitude changes that track the grouping of this cycle in bundles described above. The L * Blackman-Tukey spectrum, in addition, reveals less distinct peaks centered around ∼ 1 and ∼ 2.5 m. The latter is depicted as filtered component next to the L * record in Fig. 1a. This cycle follows in part the bundling of the ∼ 50 cm cycle and the associated amplitude changes of the extracted ∼ 50 cm cycle (Fig. 1a). Finally, spectral peaks occur with periods of around half and one third of that of the ∼ 50 cm cycle.
The L * Lomb-Scargle periodogram (Fig. 4) shows that the 50, 21 and 16 cm spectral peaks are significant at the 80 % confidence level. The peak in the ∼ 2.5 m band is not statistically significant at the 80 % confidence level and is not present in the Lomb-Scargle spectrum. However, the wavelet analysis of L * (Fig. 5) does show a significant influence of The a * spectrum (Fig. 3) also shows a double peak around 50 cm (about 45 and 49 cm). In contrast to L * , the a * spectrum does not reveal a 2.5 m cycle, but rather well-defined ∼ 90 cm and ∼ 1.6 m cycles. The extracted components of these cycles are shown above the a * record in Fig. 1 as well. Finally, the a * spectrum contains peaks that correspond to approximately half the thickness of the ∼ 50 cm cycle. The Lomb-Scargle periodogram of the a * record (Fig. 4) shows that the ∼ 1.6 m, ∼ 90 cm and ∼ 50 cm period peaks are significant at 90 % confidence level. The wavelet plot of the a * record (Fig. 5) shows the expression of the ∼ 50 cm period and some higher frequency variations.

Short-color, elemental abundance and MS records
Color, MS and elemental records of the two selected intervals are shown in Fig. 2 with pma's. The MS record of both intervals reveals well-defined maxima associated with the 40-50 cm cycle (of the long interval) with higher frequency variations in between (Fig. 2). This 40-50 cm cycle follows a sawtooth pattern with an abrupt shift to the prominent maximum, in several cases followed upcore by two less distinct peaks of decreasing amplitude associated with the higher frequency variability. The MS spectrum (Fig. 6) reveals a distinct ∼ 45 cm peak as well as peaks consistent with ∼ 21 and 14 cm cycles that reflect the higher frequency variations in MS. The Al spectrum also shows peaks at 21 and 14 cm. The Ti/Al spectrum reveals a different pattern, showing peaks at 11 and 16 cm (Fig. 7c).
The bandpass filtered 50 cm cycle from the L * data picks up the prominent MS maxima, while the filtered 16 cm cycle (also from L * data, Fig. 7b) traces most of the observed higher frequency variability. The filtered 22 cm cycle recognizes the distinct MS maxima of the 50 cm cycle and places an extra cycle in between some of these maxima (see Fig. 7b).
The same cycles are also recognized in the long-L * record and spectrum (Figs. 1, 3 and 7a), with L * minima corresponding to MS maxima, but the sawtooth pattern in MS is less clearly seen in L * . The ∼ 50 cm cycle is clearer in the 10 pma's. The a * record and spectrum also show the ∼ 50 cm cycle (see Fig. 3), but the shorter period cycles are not readily observed in the a * record and spectrum. Elemental abundance records follow the distinct pattern observed in the MS record, presumably as a consequence of the closed sum effect of Al (and Ca; Fig. A1 in the Appendix). We examined the elemental data as ratio over Al to eliminate this effect and reveal changes relative to terrestrial input, represented by, for example, Al (Figs. 2 and A2 in the Appendix). Elemental data of Zr, Ti and Fe (over Al) are positively correlated with Al and MS variations associated with the 45, 21 and 16 cm cycles. This implies that they show a relative increase compared to Al (and MS) maxima in these cycles. On the other hand, Sr, Si and Ba (over Al) show a negative correlation with Al and MS variations in these cycles (Figs. 2 and A2 in the Appendix).
Elemental ratio spectra all reveal apparent periodicity with a cycle thickness between ∼ 45 and 55 cm, in correspondence to the a * , L * and MS spectra. Like the L * and MS data, elemental ratio spectra often reveal the ∼ 22 and ∼ 16 cm cycles as well, but less distinct and sometimes at slightly different frequencies.

Eccentricity, obliquity and precession-related variations
The ∼ 2.5 m and ∼ 50 cm cycles, which are evident from visual inspection and spectral analysis of the records, explain a large part of the variation in the color data of the long record. The filtered ∼ 2.5 m cycle in addition traces changes in the amplitude of the filtered ∼ 50 cm cycle in the lower part (1145-1153 mcd) of the record. This amplitude modulation of the ∼ 50 cm cycle by the ∼ 2.5 m cycle and the characteristic ∼ 1 : 5 ratio between these cycles suggest that the ∼ 2.5 m cycle represents the short ∼ 100 kyr eccentricity cycle, while the 40-55 cm cycle corresponds to the climatic precession cycle with a ∼ 19-23 kyr period. This interpretation results in an estimated average sedimentation rate of ∼ 2.3 cm kyr −1 and follows the interpretation of Park et al. (1993). The ∼ 50 cm cycle is thicker than the climatic precession cycles documented from Site 516F by Herbert and D'Hondt (1990), who report a lower sedimentation rate of 1.3 cm kyr −1 . In addition, a long-period cycle with a thickness of 8-10 m can be recognized in the L * record of the long interval. This cycle likely corresponds to long ∼ 405 kyr eccentricity, which is stable through geologic time (Laskar et al., 2004). These ∼ 100 and 405 kyr cycles are well expressed by the 30 pma of L * (Fig. 1a). The studied interval belongs to the middle part of Chron C33r according to the calibration of the Site 516F magnetostratigraphy to the polarity timescale . Such an early Campanian age is consistent with the planktonic foraminiferal biostratigraphy with first occurrences of Globotruncana arca and G. ventricosa recorded in the middle of core 114 and in core 112, respectively (Weiss, 1983), although the biostratigraphic significance of these events has been questioned as a consequence of diachroneity (Petrizzo et al., 2011). C33r with the Santonian-Campanian boundary provisionally placed at its base has a duration of 3.74 myr in GTS2012 (Ogg and Hinnov, 2012) while the correlative reversed polarity interval at Site 516 is 73 m thick, resulting in an average sedimentation rate of ∼ 2 cm kyr −1 . This rate is in very good agreement with our assumption that the ∼ 50 cm cycle is related to climatic precession if we start our calculation from an average precession period of 20.8 kyr for the Campanian (Berger et al., 1989; see also d'Hondt, 1990 andPark et al., 1993). The slight mismatch may be explained by the presence of slump scars, which indicate that some unknown amount of sediments is missing in cores 108 and 110 (Herbert and d'Hondt, 1990).
The cyclostratigraphic interpretation is further consistent with Herbert and d'Hondt (1990) who carried out a remarkable cyclostratigraphic study of deep-sea cores of Late Cretaceous to early Paleogene age from the South Atlantic. They  sediment color as proxy for carbonate content. All their records reveal a dominant influence of precession, while the ∼ 100 kyr eccentricity cycle is often expressed by bundles of 4 to 6 precession-related cycles, as is the case in our L * record and the MS record of Park et al. (1993). This bundling can also be observed in the core photographs, which in addition reveal the presence of long intervals without distinct precession-related cyclicity that may reflect 405 kyr eccentricity minima (as eccentricity determines precession amplitude). Unfortunately, individual cores are too short to fully portray this cycle. Such extended intervals without distinct precession-related cycles are also evident around 1149 and 1158 mcd in the L * and MS records of the studied interval (see Fig. 1).
The ∼ 2.5 m cycle, which is supposedly related to ∼ 100 kyr eccentricity, is extracted from the color records, even though it is not significant at the 80 % CL. Another support for this interpretation is that the expression of eccentricity stems from modulating the precession amplitude and not from its small direct effect on the mean global annual insolation (of < 0.25 %). A pure amplitude modulator is not picked up by classical spectral analysis, as shown by the lack of eccentricity-related peaks in climatic precession spectra; eccentricity-related peaks enter proxy spectra through nonlinear responses in the climate-depositional system, as already suggested by Park et al. (1993). This also explains why the influence of eccentricity is recorded as power spectrum maxima in the precession frequency band of the L * wavelet spectrum. The ∼ 90 cm cycle in the a * spectrum may correspond to double precession or, more likely, to obliquity. In that case, the ∼ 1.6 m a * cycle may well reflect double obliquity, as also suggested by the results of bandpass filtering between 1148 and 1153 mcd (Fig. 1b). However, a detailed comparison with the L * record suggests that the ∼ 1.6 m cycle in a * also partly traces the inferred ∼ 100 kyr cycle in L * , with a * minima corresponding to L * minima. The contrast in the L * and a * spectral characteristics suggest that these parameters reflect different parts of the climate system. The L * record correlates inversely with the MS and most elemental abundance records ( Fig. 2; see also Giosan et al., 2002). Therefore, lightness (L * ) data can in many cases be interpreted as a proxy for the concentration of non-carbonate and nonsilica matter (exported from the water column), representing terrestrial input. This is also suggested for the data presented here.
The a * record reveals different spectral characteristics, as mentioned before, and lacks the high-frequency sub-Milankovitch variability. It therefore has a different and less straightforward origin as L * . The hue (shown in part by the a * index) of the sediment is determined predominantly by the hematite content of the sediment and the relative abundance of Fe 2+ and Fe 3+ in clay minerals (e.g., Giosan et al., 2002). The inferred presence of obliquity in a * is in first instance difficult to explain as the influence of obliquity is much stronger at high latitudes. However, it is also found in low-latitude climate records, even at times when obliquity controlled glacials are absent and ice-driven responses can be excluded (e.g., Lourens et al., 1996). An alternative explanation is a link to inter-hemispheric low-latitude insolation gradients that control monsoonal activity (Rossignol-Strick, 1985;Lourens and Reichart, 1996;Leuschner and Sirocko, 2003). On the other hand, results of climate modeling of orbital extremes indicate that obliquity may exert a noticeable influence on low-latitude climate systems, such as the African monsoon, through teleconnections with high-latitude insolation forcing (Tuenter et al., 2003), although new simulations using a full GCM suggests that the obliquity signal originates from low-latitude (cross-)equatorial moisture transport (Bosmans, 2014).

Quadruplets and semi-precession
The color and MS spectra of the long interval in addition reveal intervals of increased spectral power in the sub-Milankovitch band of the spectrum. To understand this cyclicity, the selected short intervals were studied in detail using geochemical elements as additional proxy data. Sub-Milankovitch variations are particularly evident in the MS record, but are also recognized in color and elemental over Al ratios (see Fig. 2). The thicknesses of these sub-Milankovitch cycles (14 and 22 cm) can be translated into periods of ∼ 7 and 10 kyr, using the inferred sedimentation rate based on a thickness of ∼ 45-55 cm for the precession-related cycles combined with an average duration of ∼ 20 kyr for precession at 90 Ma (Berger et al., 1992).
Elemental ratios such as Ti/Al and Zr/Al follow the pattern of the sub-Milankovitch variations observed in MS within a precession-related cycle (Fig. 2). The Al and Ti/Al records in Fig. 2 suggest that the semi-precession cycle of ∼ 22 cm inferred by Park et al. (1993) is present in both records (especially in 1152 till 1153.2 mcd and 1154.9 till 1155.7 m). This points to a similar source of the terrigenous material. This observation excludes a quadripartite structure of the precession-related cycles as a possible explanation for the inferred semi-precession cyclicity in MS of Park et al. (1993). Up to now, quadripartite precession-related cycles with two carbonate minima/maxima per cycle having a different climatic origin have only been described from the Mediterranean Neogene (e.g., de Visser et al., 1989;Hilgen et al., 2003). In this case, the basic precession-related cycle reveals two carbonate (and MS) maxima per cycle, but terrigenous elements over Al (e.g., Ti/Al) reveal a single cyclic variation associated with the basic cycle itself (de Visser et al., 1989;Van Os et al., 1994). This points to two different sources of the terrigenous material (fluvial vs. eolian), and a precession-related variation in their relative contribution associated with the relative humid and arid phase of the cycle (Foucault and Mélières, 2000), while diagenetic unmixing of carbonate plays a role as well (van Os et al., 1994). Such a quadripartite explanation is also inconsistent with the presence of three MS maxima in well-defined precessionrelated cycles at DSDP Site 516F and their associated 16 cm (∼ 7 kyr) peak in the power spectra.
However, the latter is also in contradiction with a semiprecession origin of the sub-Milankovitch variability inferred by Park et al. (1993), as this would imply two MS maxima per cycle separated by 10-11 kyr (i.e., half a precession cycle). Such a semi-precession cycle was suggested by the results of their spectral analysis. However, the filtered semi-precession component alternately picks up successive prominent precession-related maxima in MS, but intermediate semi-precession cycles do not pick up sub-Milankovitch variability in between the prominent MS maxima. The less good fit between the filtered half-precession cycle and the data shows that the semi-precession cycle can at least partly be explained as an artifact of the spectral analysis, reflecting the first harmonic of the precession-related cyclicity rather than a cycle on its own (see also Herbert, 1994). In fact, the observed sub-Milankovitch variations, especially in MS and Al, are better approximated by the ∼ 7 kyr cycle (Figs. 2 and 7). As a consequence, this ∼ 7 kyr (16 cm) cycle better approximates the sub-Milankovitch variability observed in the proxy records. It is this cycle that is responsible for the observed triple-peak signature of the most prominent precession-related cycles in the short intervals. Similar ∼ 7 kyr sub-Milankovitch cycles have been observed in deeptime successions (e.g., De Vleeschouwer et al., 2012;Boulila et al., 2010).

Sub-Milankovitch variations in the Heinrich frequency band
As the ∼ 7 kyr cycle describes the actual sub-Milankovitch cycle in the DSDP Site 516F proxy records best, a climatic interpretation for this cycle is sought. The most logical explanation is to link this cycle to the influence of climatic precession between the tropics, as this would most easily explain the period shorter than that of precession in the sub-Milankovitch frequency band. At the equator, the two overhead passages of the sun per precession cycle (during the vernal and autumnal equinoxes) lead to two insolation maxima separated by half a precession cycle (e.g., Berger et al., 2006). This would explain a semi-precession cycle, but not the ∼ 7 kyr cycle. Berger and Loutre (1997) further investigated the amplitude of the seasonal cycle at the equator; the resulting spectrum reveals peaks associated with the semiprecession and quarter-precession components in addition to peaks associated with eccentricity (dominant) and obliquity (but not precession). In this way, a 5-6 kyr cycle may be explained, but not the present 7-8 kyr cycle. The influence of semi-precession and quarter-precession has been detected in paleoclimatic records both of the Pleistocene as well as older time intervals (e.g., Ferretti et al., 2010;Hernandez-Almeida et al., 2012;Steenbrink et al., 2003;Anderson, 2011). A possible clue to the origin of the 7-8 kyr cycle may come from simulations of sub-Milankovitch climate variability associated with dynamic vegetation, using transient climate model runs (Tuenter et al., 2007). In the model of Tuenter et al. (2007), these variations result from the dynamic vegetation response to precession forcing. For instance, the monsoonal runoff reveals semi-precession and quarter-precession periods in runs with interactive vegetation. However, in several instances, also periods in between 10-12 and 5-6 kyr are found. For instance, the July runoff originating from the modeled Asian Monsoon shows three peaks in a precession cycle that are ∼ 7 kyr apart ( Fig. 9c in Tuenter et al., 2007). This fits the climatic interpretation of the ∼ 7 kyr cycles at Site 516F as they are also recorded in elemental abundance ratios such as Ti/Al and Zr/Al, suggesting that they are likely to result from terrestrial sediment input via fluvial runoff into the paleo-South Atlantic Ocean (see also Park et al., 1993). Tuenter et al. (2007) also report changes in salinity as a result of the changes in runoff, particularly in semi-enclosed basins that may be considered as similar to the paleo-Atlantic Ocean. As the paleolatitude of DSDP Site 516 is at latitudes most likely influenced by trade winds, large drylands in the African hinterland may interact with vegetation in a similar way as described by Tuenter et al. (2007). However, the transient climate simulations should be repeated using more sophisticated climate-vegetation models, as an artifact related to the limited number of vegetation components cannot be excluded. The ∼ 7 kyr period of the sub-Milankovitch cycle at Site 516F fits the recurrence time of Heinrich events (massive input of ice-rafted debris in the North Atlantic) of the last glacial cycle remarkably well. Cycles with similar periods are not restricted to the last glacial, but are also found in δ 18 O and other proxy records of marine isotope stages (MIS) 96-100 at ∼ 2.4 Ma (Becker et al., 2005(Becker et al., , 2006, though these lack the typical Heinrich event signature (Bond et al., 1992; results ODP Leg 162/171 IODP 303). Our observation is in agreement with the interpretation of Heinrich events being controlled by low-latitude climate changes induced by a non-linear response to the precession cycle between the tropics. Changes in monsoonal controlled runoff on sub-Milankovitch timescales will alter salinities and temperatures in the (paleo)Atlantic and, hence, influence the meridional mass transport to higher latitudes. Sub-Milankovitch cycles with a similar period were also detected in laminae thickness counts of Permian evaporites (Anderson, 1983(Anderson, , 2011. For these cycles (as well the recorded semi-and quarterprecession cycles), a similar monsoonal origin is inferred, supporting that these originate at low latitudes.

Conclusions
Based on results of the generated data and time series analysis, it can be concluded that a dominant 10-11 kyr semi-precession cycle does not solely describe the sub-Milankovitch variability observed in high-resolution MS, XRF and color records of DSDP Site 516F cores 113 and 114. Also, a quadruplet structure of precession-related cycles as recorded in the Mediterranean Noegene can be excluded as the cycles in the sub-Milankovitch frequency band have a similar chemical signature. The spacing between the maxima in magnetic susceptibility, color and terrestrial element concentration rather suggests a dominance of a shorter 7-8 kyr sub-Milankovitch cycle instead of the semi-precession cyclicity observed by Park et al. (1993). The observation of the 7 kyr cycle in a low-latitude Cretaceous greenhouse setting supports previous claims that the Heinrich events of the late Pleistocene are likely related to sub-Milankovitch climate variability that is exported from equatorial areas to higher latitudes. This sub-Milankovitch cyclicity may result from a non-linear response to precession forcing at low latitudes, possibly related to a slow response of vegetation to changes in insolation (as modeled by Tuenter et al., 2007).