Stratigraphic continuity and fragmentary sedimentation: the success of cyclostratigraphy as part of integrated stratigraphy

Abstract The Milankovitch theory of climate change is widely accepted, but the registration of the climate changes in the stratigraphic record and their use in building high-resolution astronomically tuned timescales has been disputed due to the complex and fragmentary nature of the stratigraphic record. However, results of time series analysis and consistency with independent magnetobiostratigraphic and/or radio-isotopic age models show that Milankovitch cycles are recorded not only in deep marine and lacustrine successions, but also in ice cores and speleothems, and in eolian and fluvial successions. Integrated stratigraphic studies further provide evidence for continuous sedimentation at Milankovitch time scales (104 years up to 106 years). This combined approach also shows that strict application of statistical confidence limits in spectral analysis to verify astronomical forcing in climate proxy records is not fully justified and may lead to false negatives. This is in contrast to recent claims that failure to apply strict statistical standards can lead to false positives in the search for periodic signals. Finally, and contrary to the argument that changes in insolation are too small to effect significant climate change, seasonal insolation variations resulting from orbital extremes can be significant (20% and more) and, as shown by climate modelling, generate large climate changes that can be expected to leave a marked imprint in the stratigraphic record. The tuning of long and continuous cyclic successions now underlies the standard geological time scale for much of the Cenozoic and also for extended intervals of the Mesozoic. Such successions have to be taken into account to fully comprehend the (cyclic) nature of the stratigraphic record.

changes in insolation are too small to effect significant climate change, seasonal insolation variations resulting from orbital extremes can be significant (20% and more) and, as shown by climate modelling, generate large climate changes that can be expected to leave a marked imprint in the stratigraphic record. The tuning of long and continuous cyclic successions now underlies the standard geological time scale for much of the Cenozoic and also for extended intervals of the Mesozoic. Such successions have to be taken into account to fully comprehend the (cyclic) nature of the stratigraphic record.
Gold Open Access: This article is published under the terms of the CC-BY 3.0 license.
Astronomically induced climate forcing and its expression as cycles in the stratigraphic record play an important role in the construction of highresolution time scales and in understanding past sedimentation and climate change on Milankovitch (10 4 years up to 10 6 years) time scales. The tuning of stratigraphic cycles to astronomical target curves calculated with the help of astronomical solutions for the Solar System now underlies the age calibration of the Geological Time Scale (GTS) for most of the Cenozoic Era Vandenberghe et al. 2012). A significant part of the Mesozoic Era has been astronomically scaled for the GTS as well (Ogg & Hinnov 2012a, b;Ogg 2012). Moreover, all 40 Ar/ 39 Ar ages in the new GTS are calculated relative to an astronomically calibrated age for the Fish Canyon sanidine (FCs) dating standard (Kuiper et al. 2008;Schmitz 2012). Despite this progress, critical papers on cyclostratigraphy have been published (e.g. Miall & Miall 2004;Bailey 2009;Vaughan et al. 2011) focusing on the following points: (1) Stratigraphic successions are punctuated by hiatuses and changes in sedimentation rate and are thus by definition discontinuous and unsuitable for astronomical calibration (Miall & Miall 2004;Bailey 2009;Miall 2014). The fragmentary character of the stratigraphic record will not allow the study of Milankovitch cycles in detail and, in particular, the use of these cycles to build high-resolution time scales. It has been asserted that the stratigraphic record is 'more gap than record' (Ager 1973), but this statement depends on the time scale of reference, as sedimentation rate v. duration follows an inverse power law with an increasing percentage of time missing in hiatuses at longer time scales (Sadler 1981). Hence, it is concluded that stratigraphic continuity and constant sedimentation rates are myths that require balance between subsidence and sedimentation that in practice does not exist.
(2) The statistical basis for Milankovitch cyclicity is weak as low significance levels are often employed in combination with improper statistical treatment, which may lead to a situation in which false positives might become the norm (Vaughan et al. 2011). This point is also related to the question of whether the discrimination of a stratal (i.e. cycle) hierarchy in the stratigraphic record is real or merely an arbitrary subdivision of an uninterrupted continuum (Bailey & Smith 2008).
(3) Other critical points concern the weakness of the astronomical forcing, uncertainty in the periods of the astronomical cycles in the geological past, and insufficient precision of independent time calibration of cyclostratigraphic data to establish models based on tuning, filtering and other kinds of statistical treatment.
Here, we address these critical points first by reviewing a number of cases in which Milankovitch cycles occur in long sedimentary successions that are shown to be continuous on Milankovitch time scales. We present these cases for various time intervals, focussing in particular on the youngest, Cenozoic, part of the record because it includes sediments deposited under widely varying environmental conditions for which good time control is available. This approach is key for demonstrating that sedimentary cycles were controlled by astronomical climate forcing, that cyclic successions can be stratigraphically continuous up to million year time scales, and that these successions can be used to build high-resolution tuned time scales. We then address critical issues concerning the nature and continuity of the stratigraphic record, and the statistical analysis of cyclostratigraphy, and discuss additional concerns regarding the inferred weakness of the forcing, independent age constraints for testing the astronomical theory of climate change, the stability of astronomical frequencies in the past and the primary character of limestonemarl alternations used for tuning.

Pleistocene
The study of the influence of Earth's orbitalrotational cycles on climate, including the construction of astronomically tuned time scales, goes back to the nineteenth century (e.g. Imbrie & Imbrie 1979;Hilgen 2010 and references therein). Research focused initially on the astronomical theory of the Ice Ages, as first formulated by Adhémar (1842) and Croll (1864). Progress, however, was slow as only discontinuous records of river terraces and moraine deposits were available for study, while independent age control was essentially lacking (Imbrie & Imbrie 1979).
With the recovery of the first piston cores of deep marine sediments from the ocean floor during the first Swedish Deep-Sea Expedition in 1947 (Kullenberg 1947;Pettersson et al. 1951), continuous records of Ice Age history became available for the first time (Emiliani 1955). This recovery heralded the beginning of the revival -and the general acceptance -of the astronomical theory of the Ice Ages and the construction of astronomical time scales (Hays et al. 1976;Imbrie et al. 1984). The development of new climate proxies, statistical techniques and dating methods also played a crucial role in the acceptance of the theory (see Imbrie & Imbrie 1979).
The detailed marine records of the Ice Ages of the last 800 kyr (kiloyears) are dominated by a c. 100 kyr cycle and were tuned to boreal summer insolation or, alternatively, an ice volume model that used the astronomical parameters as input (Imbrie et al. 1984). Extension beyond 800 ka (kiloyears ago) using normal piston cores proved problematic as low sedimentation rate areas had to be targeted, which lacked adequate temporal resolution. This difficulty was solved by the adoption of a multiple hole drilling strategy in deep-sea drilling, which was developed to overcome problems of stratigraphic completeness at core breaks (e.g. Ruddiman et al. 1987). This strategy was used to extend the tuned marine oxygen isotope record back to 2.6 Ma (million years ago) Ruddiman et al. 1989), the time that marks the onset of major Northern Hemisphere glaciations. In contrast to the late Pleistocene, these older records are dominated by the 41 kyr axial obliquity cycle predicted by Milankovitch (1941) for all glacials.
The Marine Isotope Stages (MIS), based on the standardization of the characteristic pattern of benthic oxygen isotope variations as stages, provided a critical tool for the study of the Ice Ages in the marine realm. The MIS that were introduced with the publication of Emiliani (1955) allow deep-sea records to be synchronized along a common highresolution astronomically tuned age model, with the latest version being the LR04 stack of Lisiecki & Raymo (2005) based on a stack of 57 globally distributed isotope records.
The expression of the Quaternary ice ages is not limited to the deep marine realm but is also found in shallow marine and continental successions, as well as in climate archives of ice cores and speleothems. High-resolution records of sea-level change were identified in stratigraphic sequences of shallowmarine successions in the Wanganui Basin, New Zealand (Naish et al. 1998), linking sequence stratigraphy to cyclostratigraphy, and confirmed by magnetobiostratigraphic and radio-isotopic dating. Integrated stratigraphic correlations were further used to identify ice ages in continental successions, most notably in the lacustrine succession of Lake Baikal (e.g. Prokopenko et al. 2006) and Lake El'Gygytgyn (Melles et al. 2012), and the eolian record of the Chinese loess plateau. Astrochronologies have long been built for these loess deposits, which have been based on independent tuning (Ding et al. 1994) and correlation of the loess-soil alternations to the MIS (e.g. Hovan et al. 1989). This approach has led to a detailed understanding of the evolution of the East Asian monsoon .
Fluvial successions have been considered least suitable for cyclostratigraphic studies because it is assumed that tectonics and autocyclic processes dominate in the dynamic fluvial system (e.g. Beerbower 1964). Nevertheless, fluvial successions of Quaternary age have revealed the imprint of cyclic climate change (e.g. Blum et al. 1994;Törnqvist 1998), and a further example of astronomical control on a fluvial system comes from Pleistocene sequences in the Pannonian Basin (Nador et al. 2003).
Other archives of climate change in ice cores and speleothems are of great importance as well. The Antarctic ice core record extends back to 0.8 Ma, covering the last eight full glacial -interglacial cycles (e.g. Epica community members 2004). In particular, deuterium (dD), as a proxy for mid-to high-latitude temperature, reveals an excellent fit with the MIS (EPICA community members 2004; Loulergue et al. 2008). The ice cores also provide records of atmospheric CO 2 and CH 4 (Siegenthaler et al. 2005;Loulergue et al. 2008;Lüthi et al. 2008) that match the patterns observed in the marine record. The precession signal is amplified in the CH 4 record, as it probably picks up additional low latitude sources that operate independently from the ice ages (Ruddiman & Raymo 2003). A low latitude monsoon-related and precession-dominated signal is well documented in the speleothem d 18 O record of the Sanbao/Hulu caves (Wang et al. 2008).

Neogene
The Astronomical Time Scale (ATS) for the last 800 kyr was extended to the base of the Pliocene (5.33 Ma) following the introduction of the multiple hole drilling strategy in deep-sea drilling and the incorporation of land-based marine sections (Shackleton et al. 1990;Hilgen 1991a, b;Fig. 1). The ATS was extended into the Miocene with deep marine sections now exposed on land in the Mediterranean (e.g. Hilgen et al. 1995Hilgen et al. , 2003Krijgsman et al. 1999;Hüsing et al. 2007Hüsing et al. , 2009Figs 1 & 2). As in the Pliocene, integrated stratigraphic correlations to parallel sections, using magnetostratigraphy combined with planktonic foraminiferal biostratigraphy, were used to verify the inferred continuity of the successions at the Milankovitch scale (Fig. 2).
Astronomical tuning of the Mediterranean sections underlies the age calibration of the standard GTS for the Neogene Period (Lourens et al. 2004;Hilgen et al. 2012). Hilgen et al. (2006) proposed that some of the sections be formally designated as unit stratotypes for stages, with the Rossello Composite section as unit stratotype for the Zanclean and Piacenzian stages of the Pliocene, and Monte dei Corvi as unit stratotype for the Tortonian stage of the Miocene; the Global Stratotype Section and Points (GSSPs) of these stages had already been defined in these sections (Castradori et al. 1998;Van Couvering et al. 2000;Hilgen et al. 2005). Moreover, units defined by the astronomically controlled cyclicity can be designated as chronozones, that is, formal chronostratigraphic units of minor rank (Hilgen et al. 2006).
Cyclic terrestrial successions were also incorporated in the high-resolution astrochronologic framework of the Mediterranean Neogene. Independently from one another, magnetostratigraphic and 40 Ar/ 39 Ar dating revealed a precession origin for lignite -marl alternations of Pliocene age in the Ptolemais Basin of Greece (Steenbrink et al. 1999); these alternations have been correlated in detail to the marine Capo Rossello section (van Vugt et al. 1998). Miocene lacustrine successions in Spain also proved suitable for the study of astronomically forced climate change (Abdul Aziz et al. , 2004Abels et al. 2009a, b). Cycles in these lacustrine successions are dominated by precession and their number fits that of the deep marine Monte dei Corvi section for each magnetic polarity interval (Abels et al. 2009a, b;Fig. 3). The study of these continental successions contributes to a better understanding of astronomical climate forcing in the circum-Mediterranean area and may shed light on the climate system responsible for sapropel formation.
In addition to the Mediterranean, Neogene marine cyclic successions are encountered in numerous deep-sea cores. Sediments from Ceara Rise in the eastern Equatorial Atlantic are a prime example, as they reveal the expression of all astronomical parameters in sediment colour and magnetic susceptibility variations, and have been used for tuning (Shackleton & Crowhurst 1997;Zeeden et al. 2012;Fig. 4). High-resolution stable isotope records have also been generated that reveal the clear imprint of astronomical climate forcing (e.g. Holbourn et al. 2007Holbourn et al. , 2013Liebrand et al. 2011).

Palaeogene
The tuning of climate proxy records of Ocean Drilling Program (ODP) Site 1218 in the Pacific with its detailed magneto-and biostratigraphy resulted in an ATS for the entire Oligocene, providing insight into the astronomical climate forcing and the functioning of the global carbon cycle, as expressed in the carbon isotope record (Pälike et al. 2006a). Deep marine cyclic sections in the northern Apennines, including the Massignano section that houses the formally defined Eocene -Oligocene boundary, have been analysed cyclostratigraphically as well, resulting in a further extension of the ATS (Jovane et al. 2006;Brown et al. 2009;Hyland et al. 2009).
Following an Eocene gap, resulting from a relatively shallow position of the Carbonate Compensation Depth (CCD) in the ocean (Pälike & Hilgen 2008), attempts have also been made to tune the older part of the Palaeogene down to the Cretaceous -Palaeogene (K/Pg) boundary (e.g. Kuiper et al. 2008;Westerhold et al. 2008). The Paleocene ATS currently has uncertainties in both the exact number and tuning of 405 kyr orbital eccentricity cycles Westerhold et al. 2012;Renne et al. 2013). Nevertheless, sedimentary cycle patterns in the Zumaia section in Spain and Atlantic ODP Leg 208 sites exhibit the well known c. 1:5:20 ratio characteristic of c. 20 precession cycles proportionate to c. 4 short eccentricity (c. 100 kyr) and one long eccentricity (405 kyr) cycle; this interpretation fits with results of integrated stratigraphy ( Fig. 3 in Kuiper et al. 2008;Westerhold et al. 2007). Detailed cyclostratigraphies down to the precession scale have been established in particular for the interval between the Paleocene -Eocene Thermal Maximum (PETM) and Eocene Thermal Maximum (ETM) 2 in Leg 208 sites Westerhold et al. 2007) and in the classical succession of northern Italy (Galeotti et al. 2010).
Cycles in the Zumaia section have been correlated at the precession scale to the Bjala section in Bulgaria for the interval of the Danian-Selandian transition; this interval covers magnetochronozone C27n and includes the so-called top C27n carbon isotope excursion (Dinarès-Turell et al. 2010. Paleocene and early Eocene hyperthermals are further recorded as distinct carbon isotope excursions in the fluvial succession in the Bighorn Basin of the USA (Bowen et al. 2001;Abels et al. 2012). This floodplain succession also exhibits distinct cyclicity, which on the basis of magneto-and , and (c) shallow lacustrine-floodplain successions exposed in the Orera and (d) Cascante sections (Miocene, Spain). Punta di Maiata is a partial section of the Rossello Composite proposed as unit stratotype for the Zanclean and Piacenzian Stages. Sapropels at Monte Gibliscemi show characteristic cycle hierarchy with sapropels grouped into bundles, reflecting the amplitude modulation of precession by eccentricity. Fig. 2. Astronomical tuning of sapropels and associated grey marls in land-based deep marine sections in the Mediterranean for the interval between 10 Ma and 7 Ma. Colours in the lithological columns indicate sapropels (black), associated grey marls (grey) and homogeneous marls (yellow). Colours in the magnetostratigraphic columns indicate normal polarities (black), reversed polarities (white) and uncertain polarities (grey). Sapropels and associated grey marls have been numbered per section and lumped into large-scale groups (roman numerals) and small scale groups (after Krijgsman et al. 1995;Hilgen et al. 1995). The initial age model is based on magnetobiostratigraphy. Phase relations between sapropel cycles and the orbital parameters/insolation used for the tuning are based on the comparison of the sapropel chronology for the last 0.5 myr with astronomical target curves. Tuning was carried out in successive steps starting with matching the large scale sapropel bundles to long period eccentricity and ending with matching the individual sapropels to precession minima and insolation maxima. The astronomical solution used is La93 (Laskar et al. 1993). biostratigraphic constraints has been related to precession (Abdul Aziz et al. 2008a;Abels et al. 2013). This interpretation is further validated by the consistency between astrochronologic age models developed independently for the ETM2 and H2 hyperthermals in the continental and marine realm (Stap et al. 2009;Abels et al. 2012). An important outcome of these studies is that floodplain sedimentation and avulsion frequency is regionally controlled by astronomically induced climate change rather than by autogenic processes alone (Abels et al. 2013). The sections studied thus far form part of a fluvial succession that seems essentially continuous over at least one million years. In this respect, Blum & Törnqvist (2000) have stated that 'the response of pre-Quaternary fluvial  (Abels et al. 2009a, b) and the marine reference section of Monte dei Corvi in Spain Hüsing et al. 2009). The correlations and tuning are tightly constrained by the excellent magnetostratigraphy in all sections. Note the similar number of cycles per magnetochronozone/subchron, despite differences in lock-in depth and -potentially -delayed acquisition of magnetization. Cycles are numbered per section. systems to climate change is one of the most challenging and potentially rewarding research topics in fluvial sedimentology for the new millennium'.
The interpretation of lacustrine cycles in terms of precession forcing goes back to Bradley (1929) with his classic study of the Eocene Green River Formation in North America. Later work revealed the additional influence of the c. 100 kyr and 405 kyr eccentricity cycles through bundling of the precession-related cycles (Fischer & Roberts 1991;Machlus et al. 2008). Recently, cycles in the Wilkins Peak Member have been placed in a basinscale cyclostratigraphic framework (Aswasereelert et al. 2012). The detailed stratigraphic framework The initial age model for constraining the tuning is based on calcareous plankton biostratigraphy. Tuning was established from the top downwards after establishing a spliced record. Cycle patterns in core 25 are dominantly controlled by precession/eccentricity while the lower part of core 26 is dominated by obliquity. The excellent fit between the complex cycle patterns in the cores and the insolation target partly depends on values for tidal dissipation and dynamical ellipticity introduced into the astronomical solution. The abrupt switch to obliquity in core 26 is also seen in the insolation target and is related to a minimum in the very long 2.4 myr period eccentricity cycle. reveals that hiatuses are present in marginal settings, as expected, while basinal successions are continuous at Milankovitch time scales.

Mesozoic
Recently, the Maastrichtian part of the Zumaia section has been investigated in detail, using an integrated stratigraphic approach (Batenburg et al. 2012;Dinarès-Turell et al. 2013). The study was directed at establishing a carbon isotope stratigraphy and an astronomically tuned age model based on the 405 kyr cycle (this cycle is stable in the astronomical solution beyond 50 Ma) (see Laskar et al. 2011a;Westerhold et al. 2012). Combined with the nearby Sopelana section, the record covers almost the entire Maastrichtian in an essentially continuous succession (Batenburg et al. 2012 Upper Cretaceous cyclostratigraphy from the Western Interior Basin (WIB), USA, was assessed by Gilbert (1895) for an early astronomically based estimate of a 20 myr duration (with an uncertainty of 'either twice or only one-half') for Late Cenomanian-Coniacian time (Fischer 1980;Hilgen 2010). Recently, a detailed cyclostratigraphy and floating astrochronology was developed for the entire Niobrara Formation by Locklair & Sageman (2008) based on the 405 kyr cycle, using geophysical well logs and covering the entire Coniacian and Santonian stages with astrochronologic durations of 3.40 + 0.13 myr and 2.39 + 0.15 myr, respectively. Further cyclostratigraphic studies of the WIB focused on the Cenomanian-Turonian boundary interval ; this interval was extended both upwards and downwards by Ma et al. (2014) (see also the section on Radio-isotopic ages consistent with Milankovitch forcing). Evidence of Milankovitch forcing in the terrestrial record of the Late Cretaceous (Turonian-Santonian) has also been reported in lacustrine sediments of the Songliao Basin in China (Wu et al. 2009(Wu et al. , 2013a. The Lower Cretaceous rhythmic pelagic succession exposed in the northern Apennines, Italy provides another classic example of Milankovitch cyclicity. The succession is well exposed in the Contessa and Bottaccione river valleys near Gubbio (e.g. Lowrie et al. 1982). Following integrated stratigraphic studies, including cyclostratigraphy, the succession may well be continuous over tens of millions of years . Classic studies further come from the 77-m-long Piobbico core (Herbert et al. 1995;Grippo et al. 2004). A recent detailed cyclostratigraphic study of this core produced a floating astrochronology of 405 kyr cycles indicating a duration of 25.85 myr for the combined Albian -Aptian stages in an apparently continuous succession (Huang et al. 2010a).
Magneto-and biostratigraphic boundary ages provide the main independent time controls for Jurassic cyclostratigraphy. The exception is the basal Jurassic, which is highly precisely radioisotope dated and intercalibrated with cyclostratigraphy (Blackburn et al. 2013). Multi-million year long cyclic marine sequences from the Kimmeridgian -Tithonian (Weedon et al. 2004;Boulila et al. 2008;Huang et al. 2010b), Oxfordian (Boulila et al. 2010), Toarcian (Boulila et al. 2014), and Sinemurian -Hettangian (Ruhl et al. 2010;Hüsing et al. 2014) with high-resolution records of total organic carbon, carbon isotopes, carbonate content and magnetic susceptibility show evidence for distinct 405 kyr orbital eccentricity cycles. Most of these sequences have been tuned to interpreted 405 kyr cycles, resulting in a sharpening of higherfrequency power preferentially in the obliquity and precession bands.
The Triassic provides excellent examples of Milankovitch forcing in the continental successions of the Newark Basin . The Newark series consists of cyclic lacustrine deposits that are supposedly continuous over c. 25 myr. The identification of 405 kyr, c. 100 kyr and c. 20 kyr cyclicity resulted in a floating astrochronology that has been anchored to an age of 201.464 Ma for the Triassic -Jurassic boundary based on radio-isotopic age constraints from basalts overlying the main Triassic portion of the lacustrine sediments (Kent & Olsen 2008;Olsen et al. 2011;Blackburn et al. 2013). Evidence for Milankovitch forcing also comes from Triassic fluvio-lacustrine and playa deposits in the North German Basin in Germany (Reinhardt & Ricken 2000;Szurlies 2007;Vollmer et al. 2008). Finally, the Triassic hemi-pelagic rhythmically bedded chert succession of Japan, covering some 30 myr, reveals a full hierarchy of precession and eccentricity cycles, including very low frequency components (Ikeda et al. 2010;Ikeda & Tada 2013).

Palaeozoic
Investigations are underway to seek evidence for Milankovitch forcing in the Palaeozoic. A prime example comes from the upper Permian marine sections of Meishan, the stratotype for the Changhsingian Stage, and Shangsi in China (Wu et al. 2013b). These sections were used to estimate an astronomical duration of 7.793 myr for the Lopingian Epoch. Combined with multiple radioisotopic ages, this signifies an important first step towards extending the ATS into the Palaeozoic. Anderson (1982Anderson ( , 2011 used annual laminae thickness counts to interpret Milankovitch forcing of more than 260 000 marine evaporite varves in the Late Permian (Ochoan) Castile Formation. Classical shallow marine cyclothems of the Carboniferous have been related to the c. 100 kyr and especially 405 kyr eccentricity cycles (Heckel 1986(Heckel , 1994. They have been correlated from the Donets Basin in the Ukraine to their North American counterparts, suggesting a global forcing mechanism of sea level at Milankovitch timescales (Davydov et al. 2010;Martin et al. 2012). Further back in the Palaeozoic, the Devonian has produced examples of astronomical climate forcing in marine successions by mainly the precession and eccentricity. The evidence indicates an astrochronologic duration of 6.5 + 0.4 myr for the Frasnian stage (House 1985;de Vleeschouwer et al. 2012ade Vleeschouwer et al. , b, 2013. The Early Palaeozoic has an extensive Milankovitch-band cyclostratigraphy (e.g. Read 1995) that is in need of high-quality geochronologic control and re-analysis (Hinnov 2013a).

Discussion
Milankovitch and the nature of the stratigraphic record Integrated stratigraphy and completeness. The examples given in the first part of this paper provide evidence that astronomical climate forcing is recorded, that both marine and continental cyclic successions can be continuous over multi-millionyear-long time scales, and that these successions can be used to build high-resolution time scales. The evidence mainly comes from applying an integrated stratigraphy approach. Integrated stratigraphy is the combined application of multiple stratigraphic subdisciplines, including biostratigraphy, magnetostratigraphy, chemostratigraphy, cyclostratigraphy and geochronology, to solve stratigraphic issues often related to geological time (e.g. Montanari et al. 1997;Abdul Aziz et al. 2008b). In the study of Milankovitch cycles, integrated stratigraphy is used to independently test whether sedimentary cycles are related to astronomical climate forcing by precession, obliquity and eccentricity, and whether successions are continuous at the Milankovitch time scale Hüsing et al. 2009). It remains difficult to demonstrate such continuity by showing that all cycles with the shortest orbital period are recorded. In fact, this is at present only possible for the Neogene where initial magnetobiostratigraphic and radio-isotopic age models were used as the starting point for a stepwise tuning; large(r)-scale cycles were first tuned to eccentricity followed by the tuning of small-scale cycles to precession and insolation. The astronomical target curves show that all precession-and/or obliquity-related cycles are recorded (e.g. Lourens et al. 1996;Hüsing et al. 2009). The development and application of the marine isotope stratigraphy (Lisiecki & Raymo 2005), fully integrated with magnetobiostratigraphy, tells the same story for the Plio-Pleistocene. Such a continuity does not only hold for cyclic deep marine successions but also for continental successions of Neogene age (Abdul ) (see also Fig. 3). For older time intervals, the integrated stratigraphic approach combined with an exact match in number of cycles between cyclic successions and target curves on the shortest Milankovitch time scales is not yet possible. Here, integrated stratigraphy is used to show that Milankovitch cycles are present. This approach further reveals that no major gaps are present, and there is no reason to assume that these successions might not be continuous, as we consider it unlikely that continuous successions are restricted to the Neogene. In an increasing number of cases, high-resolution precise radio-isotopic age determinations are fully consistent with, and thus confirm, the initial Milankovitch interpretation of the cyclicity. Issues related to the Milankovitch interpretation of cyclic successions, such as nature and continuity, independent testing and strength of forcing, are discussed in more detail below.
Chaos, continuity and sedimentation rate. Sadler (1981) used large compilations of accumulation rates and their dependence upon the measured time span to show that sedimentation rates follow an inverse linear relationship when plotted against time on a log-log scale, with proportionally more time missing at longer time scales. Such a negative power law is considered to be characteristic of fractal behaviour. Plotnick (1986) used the fractal 'Cantor bar' model of Mandelbrot (1983) to explain the log-linear relationship between sedimentation rates and the time span of Sadler (1981). Plotnick's hypothetical section showed an ever-increasing number of hiatuses at shorter time scales and that more time is missing at longer time scales. However, it is difficult to distinguish such a model in the real world from one where hiatuses are controlled by Milankovitch forcing (Kemp 2012). This is especially the case because shallow marine successions are notoriously difficult to date accurately, a problem that has troubled sequence stratigraphy from the beginning (e.g. Miall 1992). In that sense, shallow marine successions are more likely to follow the inverse power law of Sadler (1981) between sedimentation rates and time than deep marine archives and to a lesser extent (deep) lacustrine successions.
Alleged fractal attributes of the stratigraphic record, such as the increase in cumulative length of hiatuses or its self-similarity and non-scale dependent nature, suggest that sedimentary processes are governed by non-linear dynamics and chaotic behaviour (Bailey 1998). Chaos theory predicts complex non-random responses from systems in which feedback mechanisms and thresholds are important, potentially competing with the Milankovitch hypothesis as an explanation for sedimentary cyclicity (Smith 1994). Complex dynamic systems may involve pseudo-periodic repetition, but lack predictability. Smith (1994) does not envisage chaos theory as an alternative explanation for Milankovitch cyclicity, but rather that external forcing controlled by the astronomical parameters may interact with a complex dynamical system (i.e. the climate and depositional system), for example through phase locking with the external oscillator and reinforcing or dampening the original forcing. In that case, the initial astronomical forcing will still be preserved as cyclic variations in the stratigraphic record. Indeed, the response of the climate and depositional system to astronomical forcing is expected to include non-linearity and thresholds (see section on Rectification and distortion), but this does not preclude that the initial forcing is recorded as cycles in the stratigraphic record. Bailey (1998) argues that the dynamic systems that govern the stratigraphic record are so complex and chaotic, and their output so repetitive, that it is tenuous to assume that any recorded cyclicity may reflect the initial cyclic forcing by a deterministic periodic system. Algeo & Wilkinson (1988) argue that the recurrence of regular alternations in the Milankovitch frequency band is coincidental and related to sedimentary processes constrained by subsidence. They cite fluvial channel migration and deltaic lobe switching as examples of sedimentary processes governed by internal autogenic processes and thus by a non-linear dynamic system. However, there is evidence that river avulsion may in some cases be dictated by astronomically forced climate change rather than by autogenic control (see sections on Pleistocene and Palaeogene); this evidence is supported by independent time control (Abels et al. 2012). Thus, an integrated stratigraphic approach and, in particular, independent radioisotopic age control is critically important to distinguish among different working hypotheses for stratigraphic cyclicity, as the Milankovitch theory has well-defined expectations in terms of cycle thickness and scale. It is this approach that is advocated in the present paper. Kemp (2012) modelled the behaviour of the shallow marine depositional system, starting from a cyclic model of sedimentation in combination with stochastic variability. Accordingly, hiatuses pervade successions at shorter time scales potentially associated with Milankovitch control, but these successions may be continuous on time scales equal to and longer than the forcing period. Kemp further observed a step in the Milankovitch frequency band in the data of Sadler (1981) from shallow marine settings, and linked that step to the prevalence of hiatuses related to the c. 100 kyr cycle (Kemp 2012, his Fig. 4). Sadler (1999) arrived at the same Milankovitch interpretation of hiatuses to explain the observed slope steepening in a sedimentation rate v. time plot for shallow marine successions (see also Kemp & Sadler 2014).
Deep-sea data were included in the analysis of the continuity of the stratigraphic record in a similar way (Anders et al. 1987;Sadler & Strauss 1990). Anders et al. (1987) examined pelagic sediments from the DSDP database and concluded that sediment is preserved in at least 65% of 100 kyr scale intervals. Thus they supported the conclusion that pelagic environments produce sequences of calcareous oozes with only rare hiatuses and generally higher completeness. The deep marine archive is most suitable for demonstrating the registration and preservation of astronomical climate forcing in the form of Milankovitch cycles in the stratigraphic record and for using the cycles to build astronomically tuned time scales with unprecedented accuracy, precision and resolution. Before the recovery of the first deep-sea piston cores, it was generally assumed that deep marine records would be continuous. However, the deep marine record proved more fragmentary than anticipated, from disturbances caused by, for example, deepsea currents, basin starvation and slumping as a consequence of earthquakes and slope oversteepening (e.g. Keller & Barron 1983;Aubry 1991). Nevertheless, various areas remained relatively undisturbed during prolonged time intervals. These areas are often targeted for palaeoclimate-oriented legs in deep-sea drilling as they are located away from the continental margins on submarine highs and suitable for recovering the pelagic signal above the CCD. This approach makes carbonate-rich successions with higher sedimentation rates available that are excellent archives for palaeoclimatic studies and astronomical tuning. Today, much is known about the seafloor from previous drilling and seismic surveys, and sites can be carefully selected to ensure stratigraphic continuity for particular intervals. For instance, temporal reconstruction of the CCD played an important role in the selection of IODP (International Ocean Discovery Program) Leg 320/321 drilling sites (Pälike et al. 2009(Pälike et al. , 2012. As a consequence, numerous long records are now available from the ocean that are continuous often for several millions to tens of millions of years, as proven by integrated stratigraphic studies. Hiatuses are detected using integrated stratigraphy (e.g. Gale et al. 2011) or spectral analysis (Meyers & Sageman 2004).
The analysis of Anders et al. (1987) does not take into account the deep-sea sections and cores that were subsequently drilled and used to construct the Cenozoic ATS. The age control of these sections and cores, which is based on astronomical tuning and independently confirmed by magnetobiostratigraphic and radio-isotopic dating, is excellent. At the Milankovitch time scale, these successions are continuous over millions of years and have sedimentation rates that are near constant or only slightly varying. The variations may also be related to the cyclicity itself (Herbert 1994; Van der Laan et al. 2005) or represent longer-term tectonic or climatic trends. The sediment accumulation rates of these successions will likely have fluctuated on short time scales below 10 4 years. Sedimentation in the pelagic realm, for example, is partly dictated by seasonal changes in carbonate -biogenic opal production and terrigenous input (e.g. Turner 2002). As a consequence, El Niño and, on longertime scales, centennial-and millennial-scale cycles will have had their impact. Nevertheless, it has been shown that sedimentation rate can be near constant at 10 4 -yr to occasionally 10 6 -yr time scales. Consequently, these successions will plot as a horizontal line in sedimentation rate v. time on a log-log scale (see Fig. 3 of Anders et al. 1987), with sedimentation rate being essentially the same over five to seven temporal orders of magnitude. Prime examples of such successions are the Capo Rossello composite and Monte dei Corvi-La Vedova sections of the Mediterranean Neogene (Hilgen 1991a, b;Lourens et al. 1996;Hüsing et al. 2009), Ceara Rise and Walvis Ridge sites in the Atlantic (e.g. Lourens et al. 2005;Westerhold et al. 2007;Liebrand et al. 2011;Zeeden et al. 2012) and ODP Site 1218 for the Oligocene in the Pacific (Pälike et al. 2006a). As a consequence, these sections/cores do not follow the inverse power law between sedimentation rate and time of Sadler (1981) below 10 6 -yr time scales. However, the treatment of, for example, Sadler (1981), is a statistical one and does not exclude that successions are continuous on the Milankovitch time scale over 10 5 -10 6 years. Another example of a continuous marine record comes from the evaporite cycles of the Permian Castile Formation where Anderson (1982) demonstrated Milankovitch control of sedimentary cycles on annual lamina thickness counts, comprising in total more than 260 000 years. Annual laminae were shown to be continuous for distances up to 113 km in the basin (Anderson et al. 1972). These laminae further revealed the presence of sub-Milankovitch periodicities in addition to the annual cycle and precession and eccentricity control (Anderson 2011). It suggests an essentially continuous and only slightly variable sedimentation rate over five temporal orders of magnitude.
Such examples are not restricted to the marine record but include the continental record as well. Bradley (1929) used laminae thicknesses in the oil shales of the Eocene Green River Basin to underpin his precessional interpretation of sedimentary cycles, which are now known to be part of c. 100 kyr and 405 kyr eccentricity related bundles (Roehler 1993;Machlus et al. 2008;Meyers 2008). As with the Permian evaporites, this implies near constant and continuous sedimentation rates over at least six temporal orders. Other examples come from lacustrine successions of the Miocene in Spain and the Triassic Newark Basin succession in North America, while in these cases continuity and near constant sedimentation rates are not shown to start at the annual scale.
Fluvial successions may also be continuous over millions of years at Milankovitch time scales. An important example is the Eocene Bighorn Basin in North America (see section on Palaeogene). In this case, the formation of long and continuous fluvial successions occurs in settings favoured by relatively high subsidence rates. Despite continuity in the fluvial succession in the Bighorn Basin at the Milankovitch scale, sedimentation within the shortest precession-related cycles is likely discontinuous and sporadic (Abels et al. 2013). As a consequence, sediment accumulation rates remain constant from the precession cycle upwards over two to three temporal orders of magnitude, again not following the inverse relation between sediment accumulation rate and time. At the same time, successions deposited in a more marginal setting of the same basin might be less continuous and follow the inverse rule. This is likely also the case for the lacustrine successions of the Eocene Green River Basin where marginal successions do not record all the forcing cycles and are therefore less complete than the basinal successions (e.g. Aswasereelert et al. 2012).
Shallow marine and continental successions are vulnerable to erosion and are thus particularly susceptible to hiatuses. This opens the possibility that the inverse relation observed by Sadler (1981) is partly biased towards stratigraphic successions affected by (global) sea level located in marginal marine settings. Such a bias may also be related to the vast literature and interest in eustasy and sequence stratigraphy, while the cyclostratigraphic community has focused mainly on deep marine archives in the search for continuous records of climate change. The deep marine archives represent a significant portion of the total archive and should not be overlooked when exploring the nature of the stratigraphic record. Ideally, these complementary, and not opposing, views of the stratigraphic record should be reconciled before a true understanding of the nature of the stratigraphic record can be achieved. Steps in this direction have been made by Schlager (2010) and Kemp (2012).
The significance of cyclostratigraphic spectra Hypothesis testing. A fundamental problem in cyclostratigraphy is whether or not the hypothesis of astronomical forcing (H1) is supported by representative data. H1 is presented as a time series with astronomical frequencies -for example, insolation -presumed to be consistent with the data. The null hypothesis (H0) is taken as the case for no Milankovitch forcing, for which a noise time series (null model) is assumed that is also consistent with the data. The goal is to attempt to reject H0 by comparing the data with the null model within the statistical constraints of the data, and to accept H1. Two errors accompany this procedure: Type 1 errors, or 'false positives' or 'false alarms,' that is, rejecting H0 when H0 is true, and Type 2 errors, or 'false negatives,' that is, accepting H0 when H1 is true.
In the basic hypothesis test, the data spectrum ('spectrum' short for 'power spectrum') is compared to a noise (the 'null') spectrum. Spectrum estimators have statistical properties that allow the construction of confidence intervals defined by a (x 2 -distributed) probability level a. The lower confidence limit (CL) of the data spectrum is compared with the noise spectrum; if data power at the lower CL exceeds noise power at a given frequency (or frequencies), H0 may be rejected. The practice has developed to graph the noise spectrum, using the data spectrum CL factor to reposition the noise spectrum at an equivalent 'significance level' relative to the data spectrum. This allows convenient assessment of multiple CLs in terms of noise spectrum significance levels plotted together with the data spectrum.
Recently, Vaughan et al. (2011) identified shortcomings in this approach. First, there is usually no accounting for multiple tests in the case of many independent frequencies. This means that the probability level a used to set the threshold confidence interval should be adjusted downwards: this is known as the Bonferroni correction. A discussion of multiple tests in climate spectral analysis is given in Mudelsee (2010). Second, the first-order autoregressive (AR(1)) model commonly adopted as the null model for hypothesis testing often does not describe the data well. Alternative approaches include a simple or bending power law null spectrum (Vaughan et al. 2011;Kodama & Hinnov 2014).
The sedimentation rate problem. Sedimentary cyclicity implies changes in sedimentation rate and is inherent to cyclostratigraphy (Herbert 1994). The result is an uncertain, distorted time scale and frequency dispersal of spectral power (Meyers et al. 2001;Westphal et al. 2004). The problem leads to elevated false negatives if a is too small. Therefore, a balance must be found between the issues of false positives v. statistical power (or false negatives). The overriding challenge is to estimate the sedimentation rate variations.
Geologists have long recognized the role of the insolation 'canon' as a built-in time scale for cyclostratigraphy. Thus arose the practice of 'astronomical tuning' to estimate and reduce the distorting effects of variable sedimentation rates. Astronomical tuning can be as all-encompassing as matching a cyclic data sequence to an assumed insolationbased model -a long-used technique (e.g. Hays et al. 1976;Imbrie et al. 1984;Lourens et al. 1996Lourens et al. , 2004Lisiecki & Raymo 2005; and countless others, including results presented in this paper, Figs 2-5 & 7) -or as simple as tuning to a single frequency, for example, the 405 kyr eccentricity cycle (e.g. Grippo et al. 2004;Huang et al. 2010a, b;Wu et al. 2013b; and many others). An independent time scale is needed for initial calibration to an astronomical target, for example, radio-isotope dating of ashes in the section that is to be tuned, or dated ashes from elsewhere and projected into the section by bio-chemo-magnetostratigraphic correlation. Astronomical tuning also has disadvantages and must be applied with caution (see section on Tuning-induced Milankovitch spectra).
Rectification and distortion. Stratigraphic distortion leads to dispersal of spectral power with the emergence of side bands and harmonics as spectral artefacts (Herbert 1994;Meyers et al. 2001;Westphal et al. 2004). The cycle with the highest frequency (usually related to precession) undergoes the most intense distortion and peak broadening. In case of multiple distortions from eccentricitymodulated changes in sedimentation rates, it will become hard to recognize the higher frequency periodicities with conventional spectral methods even in the case of a purely deterministic sedimentary series (see Fig. 8 in Herbert 1994). Importantly for our discussion, this distortion is accompanied by an apparent lowering of the significance level of spectral power.  ) modelled a sinusoidal precession index forcing function, which is subsequently distorted by a non-linear response of the climate system, followed by non-linear recording in the stratigraphic domain, and lastly, subjected to effects of bioturbation. With each step, more power is transferred from the precession into the eccentricity band: the final spectrum is dominated by eccentricity, while precession-related peaks are almost eliminated. This modelled spectrum compares favourably with the carbonate spectrum of the Albian Fucoid Marls in the Piobbico core (Italy). The dominance of eccentricity in the spectrum cannot be explained by the direct effect of c. 0.25% of eccentricity on annual global insolation, but results from the eccentricity modulation of the precessional amplitude (e.g. Fig. 7 in Grippo et al. 2004). This implies that the precession cycle is present but masked by the disturbing effects of non-linear responses, thresholds and mixing. As a consequence of these complications, the power in the precession band is significantly reduced. Nevertheless, precessionrelated cycles can still be visually detected as thin black shale layers bundled in clusters that reflect the eccentricity modulation of the precession amplitude (see Fig. 11 in Hinnov 2013b). It is important to emphasize that eccentricity does not register in the spectrum of the precession index (or insolation forcing), as it only modulates precession amplitude. Thus, the presence of eccentricity in palaeoclimatic spectra is explained by non-linear or rectifying responses of the climate and/or depositional system to astronomical forcing (Weedon 2003;Huybers & Wunsch 2003). The eccentricity modulation of the precession signal can be investigated in these records (e.g. Hinnov 2000) although with caution (Huybers & Aharonson 2010;Zeeden et al. 2013). Such distortions also point to the problem of false negatives rather than that of false positives in the spectral analysis of palaeoclimatic records. The precession-related spectral peaks in  hardly rise above background, despite the fact that precession completely dominated the original forcing.
Examples. We demonstrate spectral analysis and hypothesis testing of two sedimentary records in which the presence of Milankovitch cycles is generally accepted: (1) the marine carbonate record of the Pliocene Capo Rossello composite section (RCS), Sicily (5.3-2.7 Ma) (Hilgen & Langereis 1989;Langereis & Hilgen 1991), and (2) the marine benthic oxygen isotope and weight percent aluminium records of the Late Miocene Ain el Beida (AEB) section, Morocco (6.5-5.5 Ma) Van der Laan et al. 2005. The astrochronology of both records is supported by independent magnetobiostratigraphic age models, which are based on magnetostratigraphic calibration to the geomagnetic polarity time scale (GPTS). The purpose of the analysis is to assess the (in-) adequacy of the AR(1) null model and to evaluate false negatives at the 99% CL. To highlight typical issues arising in these analyses, we present two approaches: (1) Lomb-Scargle (L-S) spectral analysis for unevenly spaced time series using REDFIT (Schulz & Mudelsee 2002); and (2) prolate multi-taper spectral analysis for evenly spaced time series using the Matlab procedure of RedNoise_ConfidenceLevels (Husson 2013). The confidence levels in the examples were estimated without considering a multiple test correction as advocated by Vaughan et al. (2011), although REDFIT provides information for such a correction ('Critical false-alarm level'; see Figs 6 & 8). While we acknowledge the strict statistical view taken by Vaughan et al. (2011), we also have to take account of competing problems such as 'false negative' assessments resulting from stratigraphic distortions, coupled with the low spectral power for the short eccentricity (c. 100 kyr) terms.
Pliocene Capo Rossello, Sicily. The RCS covers the entire Pliocene in a rhythmic deep marine succession (Hilgen & Langereis 1989). The tuning of the RSC underlies the standard GTS for this time interval. The origin of the carbonate cycles is complicated as carbonate dilution, dissolution and productivity all play a role (Van Os et al. 1994). A special characteristic of the basic precessionrelated carbonate cycles in the RCS is their quadripartite build-up with two carbonate maxima and minima per cycle (Fig. 5). The minima occur in the grey and beige marl beds of the basic greywhite-beige-white colour cycles. The grey marl beds have been tuned to summer insolation maxima    (Fig. 5), as they are equivalent to sapropels, which are known to correspond to insolation maxima (Lourens et al. 1996(Lourens et al. , 2004. The RCS carbonate record was originally analysed using the Blackman-Tukey correlogram method, applying an 80% CL (Hilgen & Langereis 1989). Here we analyse the RCS carbonate record, untuned in the stratigraphic domain and tuned in time to the La2004 astronomical solution (Lourens et al. 2004;Fig. 5). The astronomical tuning assigns grey marl and sapropel midpoints to maxima of the 658N summer insolation curve. The sedimentation rate curve that results from this tuning is shown in Figure 5.
The untuned carbonate L-S spectrum (Fig. 6a) reveals peaks of 15-20 m and of 0.95 m and 0.85 m that are significant at 99% CL: these correspond roughly to long eccentricity (405 kyr) and to 23 kyr and 19 kyr precession. Another peak above 99% CL occurs at 0.5 m, which is close to the Nyquist frequency of the original sample set. This results from the quadripartite structure of the precession-related cycles. The c. 5 m peak associated with the short eccentricity cycle is significant at 95% CL, while an obliquity-related peak at c. 2 m falls far below these CLs. This obliquity influence is weak compared to the precessioneccentricity combination, but it is consistently found in the Mediterranean Neogene, and precession-obliquity interference patterns in the RCS reveal a close fit with the astronomical target curve (Lourens et al. 1996;Hilgen et al. 2003).
The tuned carbonate L -S spectrum (Fig. 6b) reveals 405 kyr, 19 kyr and 11.5 kyr peaks significant at 99% CL, 24 kyr and 22.5 kyr peaks significant at 95%, and 124 kyr, 95 kyr and 41 kyr peaks significant at 90%. Enhancement and sharpening of spectral peaks at the obliquity and precession frequencies is expected, due to tuning, to a mix of obliquity and precession in the insolation target curve. However, the low (90%) CL of the obliquity peak is unexpected, because it is a tuned frequency. The appearance of eccentricity terms is also not expected and can be interpreted as independent evidence for astronomical forcing. The eccentricity is present due to signal rectification of the precession forcing by deposition (see discussion above) . Moreover, the observed CLs of the long (405 kyr) and short (124 kyr and 95 kyr) eccentricity terms follow expectation: over multi-million year-long time intervals, the long eccentricity cycle advances as a single 405 kyr term and therefore registers at a high spectral CL, but short eccentricity continually fluctuates in periodicity between 132 kyr and 95 kyr and cannot achieve a high spectral CL (e.g. Table 3 in Meyers 2012). This is reflected in the .99% CL of the 405 kyr peak and the c. 90% CL of the 124 kyr and 95 kyr peaks in the tuned spectrum. The presence of two carbonate minima per precession-related cycle is reflected by the dominant 11.8 kyr and 11.5 kyr peaks.
The multi-tapered spectra ( Fig. 6c & d) bear out similar results as the L -S spectra, but there are also significant differences. The resolution of the spectra is similar: the multi-tapered spectra have eight degrees of freedom (dofs) compared with seven dofs in the L -S spectra. The estimated AR(1) null spectra for the L-S spectra are very 'white', that is, there is little decline to lower power toward the Nyquist frequency. In the multi-tapered spectra, the AR(1) null spectra, which are computed with the data linearly interpolated to the average sample spacing, have a classic 'red' structure, tapering to low power toward the Nyquist frequency. In the multi-tapered spectra, not one but two high-power spectral peaks are measured in the lowest frequencies, at 42 m and 16 m in the untuned spectrum, corresponding to 1425 kyr and 405 kyr periods in the tuned spectrum. Possibly the original non-uniform sampling has a systematic bias that enhances low frequencies when the data are linearly interpolated. In the untuned precession band, one peak at 0.95 m exceeds the 99% CL, although in the tuned multitapered spectrum all three precession terms (24 kyr, 22.5 kyr and 19 kyr) are present at the 99% CL. As with the tuned L-S spectrum, two short eccentricity terms are well separated and visible at 124 kyr and 95 kyr in the tuned multi-tapered spectrum but do not even attain 90% CL. As with the tuned L-S spectrum, there are dominant peaks at 11.8 kyr and 11.5 kyr, related to the double carbonate minima between sapropels. Other significant terms possibly related to the precession index are present at 14.8 kyr and 12.7 kyr (small theoretical terms in the La2004 solution occur at 14.9 kyr, 14.4 kyr and 13.0 kyr, see also Stability of astronomical frequencies in the past).
Miocene Ain el Beida (AEB), Morocco. An example of distortion in the stratigraphic domain is provided by the Miocene marine AEB section, where sedimentary cyclicity is dominantly related to carbonate dilution by clastic input ( Van der Laan et al. 2005). The AEB section has a reliable magnetobiostratigraphy and records the onset of the Messinian Salinity Crisis (MSC) with strong cyclic sedimentation . Changes in sedimentation rate with increases up to five times background values occur in the thickest and most prominent reddish marl layers. These have been interpreted as controlled by eccentricity-related changes in precession amplitude at times of precession minima (strong c. 100 kyr cycles in Fig. 7f ). This causes frequency displacements of the precession-related spectral peaks and broadening of the obliquity-related peak. These distortions disappear from the spectrum after applying astronomical tuning to the section, as will be demonstrated below. For the astronomical tuning, midpoints of the reddish and beige layers (Fig. 7i) were used as calibration points: reddish layers are correlated to La2004 (1,1) 658N summer insolation maxima, and beige layers to insolation minima (Van der Laan et al. 2012). Here we analyse the benthic marine oxygen isotope (d 18 O) and weight percent aluminium (%Al) records of the AEB section.
The untuned d 18 O L -S spectrum (Fig. 8A a) reveals 1.20 m, 1.82 m, 2.43 m and 3.46 m spectral peaks at the 99% CL. Following tuning (Fig. 8A b), 12.4 kyr, 19 kyr, 23 kyr and 41 kyr spectral peaks appear, all significant at the 99% CL, as well as sub-Milankovitch peaks at 6.6 kyr and 6.0 kyr (these latter, however, are at power levels that are an order of magnitude lower than the Milankovitch terms). The tuned and untuned d 18 O multi-tapered spectra have the same four spectral peaks that calibrate to precession and obliquity, but no significant sub-Milankovitch power.
The untuned %Al L-S spectrum (Fig. 8B a) reveals many significant peaks in the spectrum, except for f , 0.3 cycles/m (wavelengths greater than 3 m). It has a high-power 1.85 m peak at the 99% CL in common with the untuned d 18 O spectrum, as well as a 2.43 m peak at the 95% CL, but no c. 3.4 m peak: and there is an 8.8 m peak at the 99% CL in the %Al spectrum that does not appear in the untuned d 18 O L -S spectrum. The tuned %Al L-S spectrum (Fig. 8B b) shows how the high-power untuned peaks have now shifted into the precession band (19 kyr and 23 kyr) at the 99% CL: there is no obliquity peak, but a weakly significant 100 kyr peak at the 90% CL. Additionally, there are multiple sub-Milankovitch peaks exceeding the 99% CL at very low power levels. The tuned and untuned %Al multi-tapered spectra each have five spectral peaks that correspond to precession, obliquity and short eccentricity, but register no significant sub-Milankovitch power.
Implications. Astronomical tuning produces palaeoclimatic time series that are vulnerable to circular reasoning (see section on Tuning-induced Milankovitch spectra). The RCS and AEB records were tuned to insolation dominated by the obliquity and precession; consequently, the obliquity and precession bands of the tuned records are expected to acquire spectral peaks with high significance. The eccentricity, however, is not present in the insolation at a measurable level, and so tuned spectral terms that are sharpened in the eccentricity band provide impartial evidence for astronomical forcing.
The following conclusions can be drawn for RCS: (1) Significant sedimentation rate variations distort the time scale and hence the spectrum of the carbonate record.
(2) Correcting the RCS time scale by astronomical tuning sharpens precession and obliquity frequencies, and importantly, aligns lowfrequency variations to the three main eccentricity terms: 405 kyr, 124 kyr and 95 kyr. (3) The short eccentricity terms at 124 kyr and 95 kyr do not achieve a high CL, but this is an expected outcome (Meyers 2012) and an example of false negatives. (4) The low CL of the tuned obliquity term is evidence of another false negative. (5) The non-uniform sampling of the RCS carbonate record significantly affects the spectrum estimation and AR(1) null modelling, although in this case the applied L-S and multi-taper algorithms resulted in consistent outcomes (see Appendix for discussion on Differences in the L -S and multi-taper spectra). To summarize the AEB example: (1) The section shows evidence for significant changes in sedimentation rate related to precession amplitude that distort the spectrum (see also Van der Laan et al. 2005).
(2) The d 18 O spectrum shows significant obliquity and precession spectral peaks, whereas the %Al spectrum has significant precession and eccentricity spectral peaks and no obliquity power. (3) Eccentricity related components do not always reach a high CL and provide another example of false negatives. (4) Finally, the d 18 O data were not tuned directly, but the %Al data track the sedimentary cycles that were used for the tuning, that is, %Al is higher in the reddish marls and is lower in the beige layers, and so %Al is directly connected with the tuning. Thus, despite the tuning to insolation, which has substantial obliquity power, the tuned %Al spectrum does not have a statistically significant obliquity peak.
In both cases (RCS and AEB), statistically significant spectral peaks exceeding the 99% CL with respect to the AR(1) null model occur in the Milankovitch band. The 1.8 m and 5 m cycles in the untuned RCS carbonate record are most likely false negatives: these calibrate, respectively, to the obliquity and short eccentricity in the tuned RCS spectra (Fig. 6b & d). In the case of AEB, other proxies collected from this section have been     shown to have statistically significant short eccentricity components: records of Globigerinoides total %, nannofossils SST (sea surface temperature), and counts of planktonic foraminifera per gram correlated with the La2004 (1,1) astronomical solution all register a statistically significant spectral coherency peak in the short eccentricity band (Fig. 4 in Van der Laan et al. 2012). Thus, the short eccentricity peak in the tuned %Al spectrum (Fig. 8B b) is likely a false negative: it exceeds the 99% CL in the interpolation-based AR(1) model of the Multi-Taper method (MTM) spectrum (Fig. 8B d). The 29 kyr peak in the tuned d 18 O AEB spectrum that almost reaches the 99% CL may be evidence for a minor obliquity term. The prospects for detecting false negatives in the RCS and AEB spectra stem from the detailed integrated stratigraphy that accompanies these two records, in particular, the independent magnetobiostratigraphic age models (Langereis & Hilgen 1991;Krijgsman et al. 2004). The age models provide clear predictions (before tuning) at which frequencies spectral peaks are to be expected that are related to astronomical climate forcing. If such peaks are found close to or at the expected frequencies, then a logical conclusion is that these are related to the astronomical parameters, even if they do not reach high CLs.
In closing, a '99% CL' can be unrealistically high when analysing cyclostratigraphic records. The relevant CL will vary depending on the assumed null model, and with the degree of distortion of the recorded signal and the accompanying reduction in spectral power, and the expectation of what the CL should be. For example, short eccentricity terms are not expected to achieve a high CL. Unless these countervailing problems are taken into account, estimated statistical confidence will not provide a realistic basis for judging whether spectral peaks reflect astronomically controlled cyclicity or not. Certainly none of these problems are anticipated by a strict statistical approach. Instead, an integrated stratigraphic approach must be applied to test whether cyclostratigraphic variations are related to astronomically forced climate change. This should ideally be combined with climate modelling of astronomical extremes for the time and location under consideration (example in section on Weakness of forcing). Such an integrated approach will help to reveal which spectral peaks should be considered false positives (Type 1 error) or false negatives (Type 2 error).
Tuning-induced Milankovitch spectra. The potential for Milankovitch frequencies to be introduced into records through astronomical tuning is well known and is a serious drawback (e.g. Shackleton et al. 1995;Hinnov & Park 1998;Rial 1999;Rial & Anaclerio 2000;Proistosescu et al. 2012). Random noise series tuned to astronomical target curves will entrain periods of the target curve; however, tuning a random noise series will also reveal erratic changes in sedimentation rates, while sedimentation rate changes resulting from tuning an astronomically controlled data set can be expected to be related to the amplitude of the forcing (Herbert 1994;Van der Laan et al. 2005).
Furthermore, proxy records often reveal clear amplitude variations, which may reflect the eccentricity modulation of the precession amplitude. These real amplitude variations are not introduced by the tuning process but are an attribute of the climate proxy data themselves. Such variations can thus be used to test whether a tuned time series is consistent with an astronomical target, as was performed by, for example, Pälike & Shackleton (2000), Pälike et al. (2004Pälike et al. ( , 2006b, Shackleton & Crowhurst (1997), Shackleton et al. (1999), Westerhold et al. (2007) and Zeeden et al. (2013). Complex amplitude demodulation provides the statistical means to use amplitude for testing an astronomical tuning in the time domain (Shackleton & Crowhurst 1997). However, this method may also be affected by frequencies introduced by the tuning process, when no real amplitudes are present in the data (Huybers & Aharonson 2010). Nevertheless, given that real data amplitudes (and not frequencies introduced by the tuning process) are investigated, amplitudes of precession-related signals filtered from a geological record can be compared to orbital eccentricity, which modulates the precession amplitude. A reasonable fit between the amplitude of a precession-related component from a geological record and eccentricity is a very convincing argument for the orbital forcing of a geological time series (Zeeden et al. 2013). Datasets showing real amplitudes corresponding to amplitudes of Milankovitch parameters can be interpreted as astronomically forced. Much care should be taken when no clear amplitude variations are present, or when tuned time series are investigated.
Drawbacks resulting from tuning can be overcome by conducting time series analysis in the stratigraphic domain, or in the time domain with an independent (e.g. magnetobiostratigraphic) age model. For instance Hays et al. (1976) tested several age models using a very limited number of calibration points ultimately derived from radioisotopic dating. Huybers & Wunsch (2004) and Huybers (2007) developed a new chronology for the past two million years that is based on a depthderived age estimate: weak non-linearities involving short eccentricity and obliquity are observed in the response, which disappear if a tuned age model is assumed for the time series. Independent (i.e. non-astronomical) age models have been developed for non-sedimentary successions such as ice and speleothems, using for example, ice flow modelling and U/Th dating, respectively (Wang et al. 2008). These are not discussed here, but the results provide strong support for the presence of Milankovitch cyclicity.
A new class of 'objective' tuning techniques has recently appeared in the literature that approaches the problem of estimating sedimentation rates with statistical models. These harken back to the methods developed by Martinson et al. (1982), Schiffelbein &Dorman (1986) andBrüggemann (1992). Meyers & Sageman (2007) developed an 'average spectral misfit' (ASM) metric to evaluate the fit of a stratigraphic spectrum to an astronomical spectrum, given a range of likely sedimentation rates. The applied sedimentation rate resulting in the lowest number of misfits to a large number of realizations of randomized spectra with the same resolution constraints as the stratigraphic spectrum is taken as the most likely solution. A related method based on Bayesian inversion to evaluate the fit of test sedimentation rates to an astronomical model is described in Malinverno et al. (2010). These transformational methodologies provide statistics on the fitted sedimentation rates, and allow for failure, that is, lack of a unique solution for a data fit to an astronomical model. Miall & Miall (2004) raise the issue of whether independent age constraints are precise enough to validate orbital tuning and the interpretation of recurrent variations in the stratigraphic record in terms of orbital climate forcing. Radio-isotopic dating is often used for this purpose, either by direct dating of sections or indirectly by employing the GPTS. Here we discuss uncertainties in radioisotopic ages, radio-isotopic ages that are consistent with Milankovitch cyclicity, radio-isotopic age constraints for astronomical tuning and the intercalibration of independent numerical dating methods ( 40 Ar/ 39 Ar, U/Pb, astronomical).

Independent age constraints
Uncertainties in radio-isotopic ages. Radio-isotopic dating is often used as an independent test for the presence of Milankovitch cycles. Radio-isotopic dates can provide absolute numerical tie points while pairs or sets of radio-isotopic dates can also be used to calculate durations. When calculating durations from a set of dates derived from the same decay scheme (e.g. 40 K/ 40 Ar or 238 U/ 206 Pb), employing the same reference material (natural standard or tracer solution) and the same decay constant(s), systematic uncertainties associated with the calibration of the reference material and the decay constant(s) can be ignored. In the case of the 40 Ar/ 39 Ar system, precise and accurate knowledge of the neutron fluence monitor age and the 40 K branching decay constants is not required to precisely determine durations, given that all experiments employ the same values for these constants. Similarly, for any of the U -Pb systems, tracer solution calibration and the decay constant uncertainties can be ignored when calculating the duration between two U -Pb dates, provided that the dates are calculated using the same tracer solution. In contrast, systematic uncertainties must be taken into account when calculating durations based on dates derived from different decay schemes (e.g. between a U-Pb and a 40 Ar/ 39 Ar date) or when dates based on a single decay scheme are calculated relative to different reference materials (e.g. different U -Pb tracer solutions or 40 Ar/ 39 Ar neutron fluence monitors). These systematic uncertainties also have to be considered when absolute radio-isotopic dates are compared with ages derived from orbital tuning. This implies that for any radio-isotopic system, durations can be constrained more accurately and more precisely than absolute ages. However, in both 40 Ar/ 39 Ar and U -Pb dating, geological uncertainties may play a critical role in estimating both the duration and age (e.g. residence time, or loss of radiogenic parent or daughter). Schoene et al. (2013) provide a detailed description of uncertainties in radio-isotopic dating. The 40 Ar/ 39 Ar dating of sanidine and U -Pb dating of zircon are presently the 'gold standards' of geochronology. Researchers now report full (analytical and systematic) 2s level age uncertainties that are on the order of 0.1% for 40 Ar/ 39 Ar and U -Pb dating in inter-comparison studies (Blackburn et al. 2013;Renne et al. 2013). This is at the scale of the orbital eccentricity cycle for Mesozoic rocks.
Radio-isotopic ages consistent with Milankovitch forcing. While systematic and geological uncertainties limit the accuracy of radio-isotopic dates, testing the hypothesis that astronomical climate forcing is recorded in the stratigraphic record requires accurate and precise quantification of durations rather than absolute numerical ages. As outlined before, systematic uncertainties (e.g. the uncertainty associated with the age of the FCs dating standard) can be ignored when calculating the age difference between two dates that were measured with the same method relative to the same reference material. This allows durations to be measured accurately and precisely enough to test the potential Milankovitch origin of cyclicities in the sedimentary record.
Precession induced monsoonal variations recorded in speleothem d 18 O records of the Sanbao/ Hulu caves (China) are consistent with precise 230 Th dating (Wang et al. 2008). Single crystal 40 Ar/ 39 Ar sanidine ages for a number of ash beds in a Pliocene lignite-bearing lacustrine succession (Ptolemaïs, northern Greece) are in excellent agreement with the magnetostratigraphic age model for the same succession and with a precession origin of the lignite -marl alternations (Van Vugt et al. 1998;Steenbrink et al. 1999). The astrochronology for the Mediterranean Neogene was independently confirmed using 40 Ar/ 39 Ar and U/Pb dating (e.g. Kuiper 2004;Wotzlaw et al. 2014).
The Milankovitch origin of lacustrine cyclicity in the Eocene Green River Formation of North America is also validated through consistency between cyclostratigraphy and radio-isotopic ages (Machlus et al. 2008;Smith et al. 2010). Application of the astronomically calibrated FCs age of 28.201 Ma of Kuiper et al. (2008) to the sanidine ages of the volcanic ash beds allows a direct comparison of the Green River cycles to the astronomical solution for the Early Eocene for the first time Both studies revealed consistency between cyclostratigraphic interpretations going back more than one century (Gilbert 1895) and state-of-the-art radio-isotopic ages. The 40 Ar/ 39 Ar and U/Pb age pairs further support the astronomically calibrated FCs age of 28.201 Ma of Kuiper et al. (2008) (see under Age of Fish Canyon sanidine (FCs) standard ). The cyclostratigraphic interpretation and associated astrochronology of the lacustrine succession of Early Turonian-Early Campanian age in the Songliao Basin, northeastern China, is also consistent within error with published U -Pb dates of four bentonites (Wu et al. 2013a).
Furthermore, the late Triassic astrochronology of Olsen et al. (2011) developed on the basis of the sedimentary cyclicity in the Newark Basin has recently been corroborated in detail by highresolution U -Pb zircon ages (Blackburn et al. 2013). Consistency is further established between highprecision U -Pb dating and Milankovitch interpretation of the upper Permian sections of Meishan, the stratotype for the Changhsingian Stage, and Shangsi in China (Wu et al. 2013b). These sections were used to obtain an astrochronologic duration of 7.793 myr for the Lopingian epoch. Similarly, new high precision U -Pb zircon ages of numerous ash layers in the Donets basin (Davydov et al. 2010) are consistent with the classic interpretation of shallow marine to continental sequences of Carboniferous age in terms of astronomical forcing and glacio-eustatic sea-level change (e.g. Heckel 1986;Heckel et al. 2007).
Intercalibration of radio-isotopic and astronomical time, and age of the Fish Canyon sanidine (FCs) standard. Intercalibration of radio-isotopic and astronomical dating methods has been instrumental to developments in both fields and has revealed some fundamental limitations in the accuracy and precision of the different methods. The extension of the ATS from the late Pleistocene to the Miocene -Pliocene boundary resulted in chron boundary ages that were 3-12% older than in the existing GPTS (Shackleton et al. 1990;Hilgen 1991a, b). These discrepancies were explained by incomplete degassing of (bulk) basalt samples dated with the K/Ar method when new 40 Ar/ 39 Ar ages became available (e.g. Baksi et al. 1992;Spell & McDougall 1992;Renne et al. 1993). A decade later, the 40 Ar/ 39 Ar community realized that the full uncertainty of c. 2.5% in 40 Ar/ 39 Ar seriously compromised numerical dating by this method (Renne et al. 1998;Min et al. 2000). Following early attempts of intercalibrating 40 Ar/ 39 Ar and astronomical dating methods (Renne et al. 1994;Hilgen et al. 1997), the most widely used FCs fluence monitor was determined to have an astronomically calibrated 40 Ar/ 39 Ar age of 28.201 + 0.046 Ma (Kuiper et al. 2008). This age, adopted in GTS2012 (Schmitz 2012), has been confirmed by Rivera et al. (2011: 28.172 + 0.028 Ma), using the same method as Kuiper et al. (2008), but a different astronomically dated standard, and by Wotzlaw et al. (2013: 28.196 + 0.038 Ma), who applied U -Pb zircon dating on the Fish Canyon Tuff itself using state-of-the-art techniques. The consistency between U -Pb and astronomical dating has been substantiated by single crystal U -Pb zircon dating of several ash layers in the tuned Monte dei Corvi and La Vedova sections of northern Italy (Wotzlaw et al. 2014). Single crystal 40 Ar/ 39 Ar dating of Pleistocene silicic volcanics that incorporate an astronomically calibrated age of the FCs neutron fluence monitor yield eruption ages that achieve concordance or near concordance with CA-ID-TIMS 238 U/ 206 Pb zircon latecrystallization dates (Rivera et al. 2013(Rivera et al. , 2014Singer 2014). The use of these (or similar) radio-isotopic dating methods also nears or achieves concordance with independently determined astronomical ages for the same units, extending from the Pleistocene to the Cretaceous (Smith et al. 2010;Meyers et al. 2012;Wu et al. 2013a;Ma et al. 2014;Sageman et al. 2014). In particular, new 40 Ar/ 39 Ar -U/Pb age pairs of ash beds from the cyclic succession of the Eocene Green River Formation (Smith et al. 2010) and the Upper Cretaceous of the WIB support the FCs age of 28.201 Ma Sageman et al. 2014) and cast doubt on other ages (see below).
An astronomically calibrated FCs age has been questioned by Renne et al. ( , 2011, who used a statistical optimization model with input of both 40 Ar/ 39 Ar and U -Pb data to obtain an FCs age of 28.294 + 0.036 Ma. Channell et al. (2010) Kuiper et al. (2008) is used for computation, although others are in good agreement.
Finally, Phillips & Matchan (2013), using ultrahigh precision 40 Ar/ 39 Ar step-heating analyses of multi-crystal aliquots of FCs, calculated a revised eruption age of 28.008 + 0.040 Ma for the Fish Canyon Tuff relative to the astronomically calibrated age of the A1 ashbed of Rivera et al. (2011). However, it is not clear whether the sanidine of the Fish Canyon tuff used by Kuiper et al. (2008) and Rivera et al. (2011) also contained extraneous argon, while other factors such as isotopic fractionation during step heating may play a critical role as well. Clearly, this interpretation of a young eruption age is inconsistent with the youngest U -Pb zircon age of 28.196 + 0.038 Ma of Wotzlaw et al. (2013), although it cannot be excluded that even the youngest dated zircons crystallized some tens of thousands of years before eruption. Thus, the problem of the FCs age is still unresolved, although the numerical age seems to converge to approximately 28.2 Ma.
Radio-isotopic age constraints for astronomical tuning. The tuning of the Neogene has now been extended to the K -Pg boundary and beyond. In this older time interval, only the c. 405 kyr eccentricity cycle can be reliably used for tuning as the chaotic behaviour of the Solar System limits the accuracy of the most recent La2011 solution to c. 50 -55 Ma as far as full eccentricity (i.e. including the c. 100 kyr and 2.4 myr cycles) is concerned . As a consequence of the Eocene gap in the ATS (Pälike & Hilgen 2008), independent age constraints are needed for the tuning of cyclic marine successions to the correct 405 kyr cycle in order to establish an ATS for the early Palaeogene (Kuiper et al. 2008). The absolute accuracy of such 'floating astrochronologies' thus depends on the accuracy of the radioisotopic dates that provide the numerical tie points for the cyclostratigraphic record. In such a case, systematic uncertainties associated with the different radioisotopic dating methods have to be taken into account. Renne et al. (2013) provide 40 Ar/ 39 Ar age constraints for the tuning of the K -Pg boundary interval to the beginning of the 405 kyr minimum around c. 66.0 Ma, using either the Kuiper et al. (2008) or Renne et al. (2011) ages for the FCs dating standard (see above). The tuning of the K -Pg boundary proposed by Kuiper et al. (2008) stood at the base of the astronomical age model that underlies the early Palaeogene time scale in GTS2012 (Vandenberghe et al. 2012). However, an alternative ATS gives much younger ages for the boundary (65.25 Ma) and the FCs (27.89 Ma) ). This tuning differs by two 405 kyr cycles from the one proposed by Kuiper et al. (2008) and Hilgen et al. (2010), the latter being consistent with an age of 28.201 Ma for the FCs. The (early) Palaeogene time scale of Vandenberghe et al. (2012) is further consistent with single crystal U -Pb zircon ages of 55.785 + 0.034 Ma for an ash layer in the Paleocene -Eocene Thermal Maximum (PETM) on Spitsbergen (Charles et al. 2011) and with U -Pb dates between 56.0 Ma and 55.8 Ma for intrusive rocks associated with the East Greenland flood basalts of the North Atlantic Magmatic Province, one of the proposed triggers for the PETM (Wotzlaw et al. 2012). Ideally, a combination of different radio-isotopic dating methods (U -Pb, 40 Ar/ 39 Ar) should yield identical results to constrain the tuning.

Weakness of forcing
An argument that has been made against Milankovitch theory is a supposed weakness of the astronomical forcing. Bailey (2009) states that 'it seems implausible that minor, largely precession-induced, variations in terrestrial insolation, averaging + 3.5% about the mean value during the last 20 Ma , should routinely have elicited a linear (direct, proportionate and, thus, correspondingly cyclic) sedimentary response'. This view is echoed in an analysis of observed Quaternary climate change by Wunsch (2004) (although see Meyers et al. 2008) and in theoretical arguments made by Rubincam (1994Rubincam ( , 2004 (although see comment by Berger 1996 and reply by Rubincam 1996). Indeed the insolation has a weak obliquity component in a spectrum that is dominated by annual and daily cycles (Huybers & Wunsch 2003;Huybers & Curry 2006). The evidence that the palaeoclimate system has responded strongly and persistently to the astronomical parameter variations demands an explanation that is best provided by climate modelling.
One of the earliest palaeoclimate modelling studies (Kutzbach & Otto-Bliesner 1982) suggested that astronomically forced insolation changes are at least partly responsible for the different climate of the early Holocene. In particular, the monsoon system was found to respond strongly to insolation changes: the c. 7% higher July insolation 9000 years B.P. (compared to present-day values) results in a precipitation difference of 26% in the summer months over the North-African and Asian monsoon regions (Kutzbach & Otto-Bliesner 1982). Many more modelling studies have followed, showing that monsoonal precipitation can have a large response to insolation changes, likely due to the non-linear relationship between temperature and saturation vapour pressure (e.g. Kutzbach & Street-Perrott 1985;Kutzbach & Guetter 1986;Prell & Kutzbach 1987Kutzbach et al. 2008). For the Mid-Holocene, the astronomically induced changes are rather consistent among the models used in the Paleoclimate Model Intercomparison Project (Joussaume et al. 1999;Braconnot et al. 2007).
A more recent study using a high-resolution ocean-atmosphere general circulation model shows that the increased Northern Hemisphere summer insolation during the Mid-Holocene (c. 5% higher than present) results in precipitation changes of up to 46% for the North African summer monsoon (Bosmans et al. 2012). These changes are stronger than in other models (e.g. Braconnot et al. 2007), likely due to higher resolution and a more sophisticated parameterization (Fig. 9a). In particular, a modelled shift in precipitation from ocean to land, induced by precession, may be sufficient to explain precession signals in sedimentary records such as the Mediterranean sapropels (e.g. Rossignol-Strick 1985). The annual mean temperature change is non-zero due to climatic feedbacks, showing that, despite a near-zero change in annual insolation, climate responds strongly to interannual changes caused by the astronomical parameters. Figure 9b shows that annual mean temperature differences between the Mid-Holocene and preindustrial time are mostly negative. This is due to reduced insolation in boreal winter and to feedbacks that cool the monsoon-affected areas in response to strengthened monsoons induced by increased summer insolation (Bosmans et al. 2012). The cooler monsoon areas are in agreement with previous modelling (Braconnot et al. 2007).
Thus, climate and depositional systems appear to have been sensitive to small changes in insolation: various thresholds and (positive) feedbacks are at play, such as the temperature response described above and vegetation feedbacks (Kutzbach et al. 1996;Claussen 2009). For example, a closed lake system will be very sensitive to any change in net evaporation: the evidence for dominantly precession-related pluvial phases in the large lakes of North Africa and adjacent Asia is overwhelming (e.g. Crombie et al. 1997;Vaks et al. 2010).
During ice ages, the effect of changes in insolation can be strongly amplified by feedback mechanisms such as greenhouse gases and ice-albedo. This amplification was recently demonstrated in transient simulations using CLIMBER 2, a low-resolution Fig. 9. (a) Differences in annual mean precipitation between the Mid-Holocene and the pre-industrial, from EC-Earth model experiments (Bosmans et al. 2012). Annual global mean precipitation over land is increased by 2.9% (21 mm a 21 ) and by 4.5% (48 mm a 21 ) over the tropical land areas (238S-238N). Precipitation over the ocean is decreased by 0.8% globally, and by 1.2% over the tropics. The precipitation changes over the tropics are dominated by changes in the summer monsoons (see Bosmans et al. 2012). (b) Differences in annual mean surface air temperature between the Mid-Holocene and the pre-industrial, from EC-Earth model experiments (Bosmans et al. 2012). Annual global mean temperature is decreased by 0.118, and by 0.358 over the tropical land areas (238S-238N) despite a near-zero change in annual mean insolation. coupled model of intermediate complexity (Petoukhov et al. 2000), which was used in transient simulations with adequate results (Claussen et al. 1999;Tuenter et al. 2005;Weber & Tuenter 2011). This version of CLIMBER-2 (i.e. version 2.3) includes a dynamic vegetation model and a thermodynamic sea-ice model, but it does not include a carbon-cycle model and a dynamical ice-sheet model. The reader is referred to the Appendix for a more detailed description of the model and input data.
Using CLIMBER-2, two transient simulations for the past five million years were performed. For the first simulation (Run O) only the astronomical parameters from the La2004 solution  varied, while the ice sheets were kept fixed at their present day size, and the atmospheric CO 2 concentration was set to the pre-industrial value of 280 ppmv. For the second simulation (Run OIG) the same astronomical parameters as for Run O were used but the height and area of the ice sheets (Greenland, Eurasia, North America and Antarctica), as well as the CO 2 concentration were allowed to vary. The fraction of the surface area covered by ice in the spatial grid of CLIMBER-2 was obtained by summing the areas of all icecovered grid points of the ice-sheet model located within one grid cell. The change in height due to ice thickness variations was obtained by averaging all height changes in the same grid points of the ice-sheet model.
For Run O the annual surface air temperature for the grid point in which Lake Baikal is located shows an astronomically forced amplitude of about 18C while the seasonal amplitude is much larger, reaching up to 12 8C in summer (Fig. 10). For the period 3 -5 Ma, the amplitude for Run OIG is comparable to Run O, but for the last 3 million years, the CO 2 and the ice-sheet forcings strongly amplify the temperature change up to 10 8C (annual and winter) and locally a difference of about 20 8C between glacials and interglacials in summer (Fig. 10). Furthermore, in Run OIG a strong cooling trend beginning at 3 Ma is observed that is absent in Run O.
Spectral analyses of the modelled temperature series (not shown) reveal that the annual surface air temperatures respond to obliquity and precession with comparable amplitude in Run O. In the seasonal air temperature data, precession strongly dominates over obliquity. In Run OIG, the obliquity signal in the surface air temperatures is strongly amplified. This results in comparable strengths of precession and obliquity in the seasonal temperatures and a much larger obliquity signal compared to precession in the annual temperature response. Furthermore, atmospheric CO 2 and ice sheets introduce a c. 95 kyr signal in the temperature for Run OIG which is absent in Run O.
From these modelling results it can be concluded that astronomically forced insolation can cause large seasonal climate variations that can also lead to large annual variations. Furthermore, the direct astronomically forced insolation can be strongly amplified by feedback mechanisms such as greenhouse gases and ice sheets (Clark et al. 1999;Loutre et al. 2007;Yin & Berger 2012). Nevertheless, driving an Atmosphere-Ocean General Circulation Model (AOGCM) with only astronomical forcing and without additional constraints will not lead to ice ages, indicating that the feedbacks in the system are still not understood.
Stability of astronomical frequencies in the geological past Miall & Miall (2004) touch on the issue of the longterm stability of the astronomical cycles. In the geologic past, frequencies were different for precession and obliquity due to Earth-Moon dynamics. Figure 11 summarizes the spectrum of the astronomical parameters over the past ten million years. The Earth's rotation rate has been decelerating through geologic time due to the phenomenon of tidal dissipation; consequently, the Earth-Moon distance was shorter in the past, and the Moon has been receding over time (recent review in Coughenour et al. 2012). The outcome is an everlengthening precession rate for the Earth, leading to a lengthening of the precession and obliquity cycle periodicities. This effect has been calculated for different time intervals in the past (e.g. Berger et al. 1989); a decreasing rotation rate based on empirical evidence is included in the La2004 astronomical solution . The exact values of the dissipation effect are not known and depend on the geologic evidence. Independent evidence from corals, bivalves and tidalites appears to confirm a faster rotation rate in the past (e.g. Kvale et al. 1999;Williams 2000). Spectral peak ratios in a Permian red bed succession in France are also consistent with the inferred shortening (Kruiver et al. 2000).
Precession and obliquity are also affected by Earth's dynamical ellipticity. In principle, this effect can be determined through modelling mantle viscosity and convection (Forte & Mitrovica 1997): in practice, it can also be assessed through detailed comparison with cyclostratigraphy (Pälike & Shackleton 2000;Lourens et al. 2001;Zeeden et al. 2014). The latter studies reveal that present day values of tidal dissipation and dynamical ellipticity can be assumed over the last ten million years: slight adjustments may be necessary for older geologic intervals Hüsing et al. 2007). The uncertainty associated with these effects increases back in time and implies a risk that precession-and obliquity-forced stratigraphic cycles may end up tuned to the wrong astronomical cycle (Lourens et al. 2004). However, the amplitude modulation of the precession by eccentricity provides a means to correct for such errors as far back as 50-55 Ma.
Eccentricity remains stable back to c. 50-55 Ma, according to the most recent astronomical solutions La2010 and La2011 (Laskar et al. 2011a;Westerhold et al. 2012). Various terms such as the 2.4 myr eccentricity modulation cannot be calculated reliably further back in time due to the chaotic behaviour of the Solar System (caused by the asteroids Vesta and Ceres: Laskar et al. 2011b). In fact, the 2.4 myr cycle could shift to a 1.2 myr period (and back) as a consequence of transient resonances in the Solar System (Laskar 1999). Such transitions should leave a marked imprint in the cyclostratigraphic record, but the timings cannot be precisely predicted from the astronomical models. Results thus far suggest that the period of this cycle was indeed shorter during part of the Mesozoic. For instance, a 1.7 myr period was found in the Triassic Newark Basin ) and a period between 1.6 and 1.8 myr in a Triassic to Jurassic chert succession in Japan (Ikeda & Tada 2013).
However, the 405 kyr eccentricity cycle remains stable over much longer time intervals (at least back to 250 Ma), although the uncertainty in its computation increases with time. It is the 405 kyr cycle that is used for tuning early Cenozoic, Mesozoic and older cyclostratigraphies (Hinnov & Hilgen 2012;Hinnov 2013a, b).

Primary origin of limestone-marl alternations
A final issue to be discussed is the inferred primary climatic origin of the limestone -marl alternations used for tuning. This primary character has been questioned (e.g. Westphal et al. 2008) with diagenetic unmixing providing an alternative mechanism to explain such alternations (e.g. Hallam 1986). It is widely accepted that diagenesis often played an important role in enhancing the sometimes subtle primary variations in carbonate content, but in the absence of primary criteria, rhythmic successions have also been explained as resulting merely from differential diagenesis (e.g. Munnecke & Samtleben 1996). For that purpose, criteria for a primary origin of calcareous rhythmites were presented as a means to assess their suitability for high-resolution stratigraphic applications and orbital dating (Westphal et al. 2010). The problem encountered in the Monte dei Corvi section (Westphal et al. 2008) is Fig. 11. Power spectrum of Earth's astronomical parameters, valid for the past ten million years. The nominal La2004 astronomical solution ) was processed at Dt ¼ 1 kyr over 0 -10 Ma as an ETP time series (Imbrie et al. 1984). Its spectrum is displayed above as the squared modulus of the Fast Fourier Transform, that is, spectral density v. frequency. The orbital eccentricity, obliquity (tilt) and precession index frequency bands are indicated. of special importance as it hosts the Tortonian GSSP and its tuning underlies the GTS for the interval between 8 Ma and 12.5 Ma. The primary character of the cycles at Monte dei Corvi stems from the characteristic cycle patterns, which reflect precession-obliquity interference in addition to the precession-eccentricity hierarchy Hüsing et al. 2007). It is further substantiated by the high-resolution proxy (e.g. Ti/Al) records from the slightly older cycles in the La Vedova High Cliff section, which clearly indicate a primary origin (Mourik et al. 2010). The lack of indications for a primary origin at Monte dei Corvi may come from problems in recognizing the complex quadruplet structure of the basic precession-dominated cyclicity in this section.

Conclusions
From the available literature and collective experience of the community of scientists involved in cyclostratigraphy, we conclude that: † Sedimentary successions can be continuous at Milankovitch time scales over millions of years. This is shown for deep marine and lacustrine successions, and probably even some fluvial successions. In Cenozoic, Mesozoic and Palaeozoic examples, this continuity and the presence of Milankovitch cyclicity has been independently confirmed by magnetobiostratigraphic and/or radio-isotopic dating. † There are many examples of continuous successions at Milankovitch time scales with relatively constant sedimentation rates over multiple temporal orders of magnitude that do not follow the empirical inverse power law between sedimentation rate and time. This needs to be taken into account for the nature of the stratigraphic record to be understood in all its facets. † The authenticity of Milankovitch cycles in the stratigraphic record is substantiated by climate modelling, in both transient experiments and time-snapshot experiments for astronomical extremes, although much work still needs to be done to explain Milankovitch variability in stratigraphy (e.g. the Mid-Pleistocene Transition, high-magnitude 100-kyr cyclicity, sea level forcing, oceanographic forcing, etc.). † The use of high significance levels in spectral analysis is not realistic in view of the statistics of the astronomical forcing signal and the nature of the (cyclo-) stratigraphic record with timedepth distortions stemming from bioturbation and non-linear responses in both climate and depositional systems and resulting variations in sediment accumulation rate.
The two referees R. Bailey and G. Weedon are thanked for their thorough review of the original version of the manuscript. D. Smith is thanked as editor for his comments and for his patience when the revision of the manuscript was delayed.

Appendix
Differences in the L-S and multi-taper spectra There are striking differences between the L -S and MTM spectra and their null models, especially for the RCS carbonate series (compare Fig. 6a-d). In the RCS spectra, the most noticeable difference is in the very low frequencies. Computation of the L -S spectrum requires dividing the input series into overlapping segments (in this case, five 50%-overlapping segments), which means that for the tuned RCS spectrum periodicities lower than 2.623 cycles per 3 myr ¼ 0.874 myr per cycle cannot be measured. MTM spectrum estimation, while requiring interpolation of the tuned RCS series to a uniform sample spacing (here taken as the average sample spacing of 5.569 kyr) does permit estimation of periodicities as low as 2.623 myr (i.e. the Rayleigh frequency resolution is Df ¼ 1/(2.623 myr)). Thus the tuned RCS MTM spectrum (Fig. 6d) registers a low-frequency term at 1/(1425 kyr); in contrast, the tuned RCS L-S spectrum (Fig. 6b) indicates no power at this frequency. For the AR(1) model shown with the L -S spectra, the estimated r (lag-1 autocorrelation coefficient) is calculated by regression on the nonuniformly sampled data, plus a modelled bias correction (Mudelsee 2002), which results in r ¼ 0.03. For the AR(1) model shown with the MTM spectra, r is estimated directly from the autocorrelation operation on the original interpolated data set and is sensitive to the interpolated sample spacing, here taken as the average sample spacing. For the RCS interpolated series, r is 0.37-0.42, that is, more than one order of magnitude larger than the regression estimate. The difference between the values of r arises mainly from methodology. The regression technique (for non-uniformly sampled series) is sensitive to high-power signal components; it is therefore recommended that the series be pre-processed to remove such high frequency components (Usage.pdf in Schulz & Mudelsee 2002). In the RCS spectra, the high-power, highfrequency spectral peaks close to the Nyquist frequency clearly have an influence on regression-based r that is not reflected in the interpolation-based r. The robust red noise modelling of Mann & Lees (1996), not shown here, addresses the signal contamination problem with robust smoothing of the spectrum to attenuate high-power signal effects and regression on the AR(1) model. However, Meyers (2012) has recently described a lowfrequency bias inherent to the Mann & Lees (1996) model, proposing an alternative LOWSPEC procedure to reduce this bias. In summary, defining an appropriate null model for cyclostratigraphic spectra is a 'grand challenge' that has yet to be solved. Further discussion of problems with the Mann & Lees (1996) method can be found in Vaughan et al. (2011).

CLIMBER-2, a coupled model of intermediate complexity
The prescribed ice-sheet areas and heights in CLIMBER-2 were obtained from a five million year simulation using the 3D coupled ice sheet-ice shelf-bedrock model ANICE (De Boer et al. 2014). The model reconstructs ice volume of Antarctica, Eurasia, North America (using a grid distance of 40 km for these ice sheets) and Greenland (grid distance of 20 km). Basically, ANICE follows the principles outlined in Huybrechts (1990) and Van de Wal (1999). The ice velocities are determined using a combination of two stress-balance approximations. For the velocities on the ice sheets, the shallow ice approximation (Hutter 1983) is used while the computations of ice-shelf and sliding velocities are based on the shallow shelf approximation (Morland 1987). These two velocities are combined following the approach used by the Parallel Ice Sheet Model (PISM) (Winkelmann et al. 2011). The mass balance model of ANICE is forced by monthly temperature and precipitation rate (De Boer et al. 2013). Over Antarctica and Greenland, precipitation is a function of the condensation temperature above the surface inversion layer (Huybrechts 2002) while for the Eurasian and North American ice sheets a precipitation model is used that includes orographic forcing of precipitation and changes in the moisture content (Roe & Lindzen 2001;Roe 2002;De Boer et al. 2013). The monthly temperature forcing consists of the continental mean (40-808N) surface air temperature and is computed by an inverse forward modelling routine (Bintanja et al. 2005;Bintanja & Van de Wal 2008;De Boer et al. 2010). This routine linearly relates the continental temperature to the difference between the modelled benthic d 18 O and the observed d 18 O 100 years later. The input (i.e. the observed benthic d 18 O) is the LR04 benthic stack (Lisiecki & Raymo 2005). The advantage of the forward modelling routine is that a self-consistent record is constructed of surface air temperature, ice volume and benthic d 18 O.
The CO 2 concentration as used by CLIMBER-2 was reconstructed using the method as described in Van de Wal et al. (2011). They derived an empirical relationship between atmospheric temperatures and CO 2 . With this relation the 5 myr CO 2 record is constructed using atmospheric temperatures as simulated by the 3D ice-sheet model (De Boer et al. 2014). As a result, the prescribed CO 2 concentration is mutually consistent with the prescribed ice sheets.