Borrower: YGM Call #: QE508 .G3956 2020 v.1 Lending String: CAU,*UUA,OBE,LLU,UUS Location: MAIN d > 'S D ■*-» D E ^—» D O C/3 CO 00 "O Patron: Journal Title: Geologic time scale 2020 / Volume: Issue: Month/Year: Pages: Chap. #4 Article Author: F M Gradstein; James G Ogg; Mark D Schmitz; Gabi Ogg Laskar Article Title: Geologic time scale 2020 / Chapter #4, Astrochronology Imprint: Amsterdam, Netherlands : Elsevier, 2020 ILL Number: 214058436 ODYSSEY ENABLED Charge Maxeost: 10.00IFM Shipping Address: Information Delivery Services Milne Library 1 College Circle Geneseo, New York 14454 United States Fax: Ariel: Email: ILL@GENESEO.EDU 214058436 98614 J. Laskar Chapter 4 Astrochronology Chapter outline 4.1 Introduction 139 4.1.1 Historical introduction 139 4.1.2 The astronomical solution 141 4.2 Eccentricity 141 4.2.1 Decomposition of the eccentricity 142 4.2.2 Mathematical intermezzo 143 4.2.3 Eccentricity modulations 144 4.3 Chaos in the solar system 145 4.3.1 Drifting frequencies 146 4.3.2 The 405-kyr g2 - g5 metronome 146 Abstract The long-term variations of the orbital and rotational parameters of the Earth are the key ingredients for the insolation forcing in the Milankovitch theory. This chapter describes the main aspects of these variations, concentrating on the aspects that are currently recovered in the stratigraphic record. A special emphasis is given to the very long periodic terms (> 1 Myr period) that modulate the astronomical solutions and that are essential for understanding the chaotic behavior of the solar system. 4.1 Introduction According to the Milankovitch theory (Milankovitch, 1941), some of the large climatic changes of the past originate from the variations of the Earth's orbit and of its spin axis resulting from the gravitational pull of the other planets and the Moon. These variations can be traced over many millions of years (Myr) in the geological sedimentary record, although the mechanisms that transfer the forcing insolation to the sedimentary variations are not precisely known. The recovery of astronomical signal in stratigraphic sequences has allowed local or global calibration of the stratigraphic records, and cyclostratigraphy is now a very active field of research. After the astronomical calibration of the Neogene Period (Lourens et al., 2004; Hilgen et al., 2012), focus turned toward the entire Paleogene Period Otologic Tinw Scale 2020. DOI: http»^/dol.on{/10.l016m978.0-12-824360.2.00004^ O 2020 Eltevier B.V. All righu reserved. 4.3.3 The g4-g3 2.4 Myr cycle 149 4.4 Inclination and obliquity 150 4.4.1 Simplified expressions 150 4.4.2 Tidal evolution 151 4.4.3 Obliquity solution 152 4.4.4 The 173-kyr s3 - s6 metronome 153 4.5 Chaotic diffusion and secular resonances 154 4.6 Discussion 155 Acknowledgments 155 References 155 (e.g., Kuiper et al., 2008; Westerhold et al., 2012, 2014, 2015; Boulila et al., 2018), covering the entire Cenozoic Era. Extending this procedure through the Mesozoic Era and beyond is difficult, as the solar system motion is chaotic (Laskar, 1989, 1990). It is thus not possible to retrieve the precise orbital motion of the planets beyond 60 Ma from their present state (Laskar et al., 201 lb). Nevertheless, the existence of a stable component in the astronomical forcing, the 405-kyr metronome (e.g., Laskar et al., 2004), has allowed the continuation of the astronomical calibration of geologic time deep into the Mesozoic Era and even into the Paleozoic Era and the Precambrian. Detailed compilations of currently available cyclostra-tigraphic records have been summarized recently (e.g., Hinnov and Hilgen, 2012; Hinnov, 2018b; Huang, 2018), and we refer to these. In this chapter, we will focus on the astronomical solution and especially on the long cycles of these solutions, with the aim to answer some of the common questions that arise in the analysis of long sequences of stratigraphic records. 4.1.1 Historical introduction During the 18th century the question of the stability of the solar system was of prime importance, as it was also necessary to decide whether Newton's law properly describes the motion of the celestial bodies (for details, refer to 139 140 PART I II Concepts and Methods 0.08 0.07 0.06 - 0.05 0.04 - 0.03 0.02 0.01 0 LAGR ■ LEV-STO BR74 • La93 La2004 • BR74 La93 La2004 time (Ma) before present FIGURE 4.1 Improvements in the computation of the past evolution of the Earth's eccentricity. LAGR is the first secular solution for the solar system, with only six planets (Lagrange. 1783.1784); Uranus was added in LEV (Le Verrier. 1840); and Neptune in STO (Stockwell, 1873); BR74 (Bretagnon. 1974) added some terms of order 2 with respect to the masses and degree four in eccentricity and inclination. La93 (Laskar. 1988. 1990; Laskar et al., 1993) is a numerical solution of the averaged equations that contains all terms up to order 2 with respect to the masses and degree six in eccentricity and inclination, the contribution of the Moon and general relativity. La04 (Laskar et al.. 2004) is a numerical integration of the full equations of motion. Laskar. 2013). A very important result of this quest was the derivation of the first long-term models for the solar system orbital evolution. A first result, of fundamental importance for cyclostratigraphy is the demonstration, at first order of the planetary masses, of the invariance of the semimajor axes of the orbits of the planets (Laplace, 1776). This result is also practically verified in the full, nonapproximated system of equations, with the major consequence that the orbital period of the Earth does not change over time. One can thus assume that the length of the year has practically not changed over the past billion years.1 By contrast, Lagrange and Laplace found that in the linear approximation of the averaged equations of motion, the eccentricity, inclination, and orientation of the orbits change significantly with time, in a quasiperiodic manner with frequencies of several tens of kyr to Myr, but in a way that does not allow for planetary collisions (Lagrange, 1778; Laplace, 1775)2; the first full computation of the long-term motion of the Earth's orbit is due to Lagrange at the end of the 18th century (Lagrange, 1783, 1784). Of course, the solution of Lagrange only includes the planets visible to naked eye (Mercury, Venus, Earth, Mars, Jupiter, and Saturn), but it already provides a very accurate representation of the Earth's orbital motion over the past million years. In Lagrange's solution, all the main features of the variation of the Earth's orbital elements are presenL but it was only after the work of Agassiz (1840), showing evidence of past Ice Ages, and the new solution of Le Verrier, including Uranus (Le Verrier, 1840, 1841, 1856) that it was advocated that the variations of the Earth's orbit could trigger the large climatic variations of the past (Croll, 1875) (see Hilgen (2010) for more historical details). The orbital solution was upgraded by Stockwell (1873) who added the contribution of Neptune (Fig. 4.1). This latest orbital solution was used by Pilgrim (1904) for the computation of the variation of the Earth spin axis evolution. 1. The relative loss of the mass of the Sun is of about 9 X 10"14/vear Iki™ • „« deduce that the mass loss of the Sun induces an increase o 9 x An TT™™™ °f ^ and Kepler'S *M period of the Earth over IGa. '° AU^ ^ ^™Jor axis and a decrease of only 1.5 hours in the orb.tal 2. Laplace's work was largely inspired by Ugrange's manuscript that was submitted in .774 (see Laskar. 2013). Astrochronology Chapter | 4 141 Nevertheless, in his theory of the insolation of the Earth, Milankovitch (1941) considered that the solution of Le Verrier (1856) was more reliable and asked his colleague Miskovic to update Le Verrier's solution for the new values of the planetary masses and to use it for the computation of the orientation of the spin axis of the Earth with respect to its orbit. After comparison to the solution of Stockwell (1873) and Pilgrim (1904), Milankovitch decided to limit his insolation computations to the most recent 600,000 years. With the use of computers, it was possible to extend these analytical computations significantly. The solution of Bretagnon (1974) for the solar system comprises 318 periodic terms, while the secular system of Laskar, 1988 (1990) and Laskar et al. (1993) contains 153,824 terms, including the averaged contribution of the Moon and general relativity. Nevertheless, these analytical perturbative methods always require some truncation in series expansions and thus have some limitations in precision. With the improvement in computer speed and numerical integration algorithms, it is now possible to directly integrate the equations of motion, as in the La2004 solution (Laskar et al., 2004). When comparing the various solutions that have been used in stratigraphic astrochronology (Fig. 4.1 A), it appears that although Lagrange solution is somewhat off in the first 500,000 years, it already provides a good measure of the qualitative behavior of the Earth's orbital solution. The other solutions are in quite good agreement over the first 600,000 years but begin to depart from one another after this date. On the contrary, the semianalytical solution La93 (Laskar, 1988, 1990; Laskar et al., 1993) is a perfect match to the full numerical solution La2004 (Laskar et al, 2004) over the most recent 2 Myr and even over the last 10 Myr (see Laskar et al., 2004). Starting with La93, the orbital solution can thus be considered as perfectly known over the past few Myr. The evolution of the precision of the solutions is particularly striking beyond 1.5 Ma (Fig. 4.IB). The difference is very large with respect to the solution of Bretagnon (1974) and Berger (1978) but insignificant with respect to the more recent La2004. 4.1.2 The astronomical solution Due to the gravitational interactions of the planets, the Earth's orbit and spin axis present significant variations in "me. The orbit precesses slowly on its plane in space (Fig. 4.2), and the equator precesses around the normal to toe orbit (Fig. 4.3). This slow precession motion of the Planetary orbits is described by a combination of periodic modes related to the precession of the perihelions with fundamental secular frequencies g, (i= 1, 2.....8) and precession of the orbital planes in space with fundamental secular frequencies st (i - 1, 2.....8) (Table 4.1). In addition, the eccentricity of the orbit and the inclination with respect to the fixed reference frame oscillates with the same FIGURE 4.2 The eccentricity e is the ratio of the distance between the two foci of the ellipse (SS) and the major axis of the ellipse (2a). At perihelion the Earth-Sun distance is a( 1 - e); at aphelion it becomes a(l + e). The horizontal line is the direction of the ascending node (Fig. 4.3). n FIGURE 4.3 Earth angular parameters. The instantaneous orbital plane of the Earth, the ecliptic of date, is referred with respect to a fixed reference frame [mean ecliptic J2000 in La2004 (Laskar et al., 2004) and invariant plane in La2010 (Laskar et al., 201 la)], with a fixed origin -i0 (equinox J2000 in U2004). The ecliptic of date is defined by the longitude of the ascending node ft and the inclination ;'. The argument of perihelion u is the angle from the line of node SN to the perihelion and the true anomaly v the angle from perihelion to the position of the Earth. The equinox of date 7 is the intersection of the equator with the ecliptic of date. The spin axis of the Earth is directed toward the North Pole (NP) and ^ is the spin angle. The obliquity e is the angle from the normal to the ecliptic of date (n) to the spin axis (NP). The precession angle v describes the motion of the spin axis of the Earth around n. The longitude of perihelion is the sum of the longitude of the node ft and the argument of perihelion ^ (x = ft + ui). It should be noted that the two angles, ft. u, are not on the same plane. frequencies. The processing of the motion of the spin axis enters an additional frequency, the precession frequency p. 4.2 Eccentricity The eccentricity of the Earth e is a measure of the shape of its orbit (Fig. 4.2). At perihelion the Sun-Earth distance (SE) is a(l _ e) and a(l + e) at aphelion. The insolation on 142 PART I II Concepts and Methods ( mniM^n-.....fc-r-rfP.p. and s, of La2004 and La2010a in arcsec 'year. La2004 La2010a •^100 Period (year) gi 5.59 5.59 0.13 231,843 82 7.452 7.453 0.019 173,913 83 17.368 17.368 0.20 74,620 g4 17.916 17.916 0.20 72,338 «5 4.257452 4.257482 0.000030 304,407 8b 28.2450 28.2449 0.0010 45,884 8? 3.087951 3.087946 0.000034 419,696 8b 0.673019 0.673019 0.000015 1,925,646 s\ -5.59 -5.61 0.15 231,843 S2 -7.05 -7.06 0.19 183,830 S3 -18.850 - 18.848 0.066 68,753 - 17.755 - 17.751 0.064 72,994 S5 0 0 -26.347855 -26.347841 0.000076 49,188 s? -2.9925259 -2.9925258 0.000025 433,079 Sa -0.691736 -0.691740 0.000010 1,873,547 A,00 are the observed variations, in arcsec/year, of the frequencies over 100 Myr (Laskar et al., 2011 a) Th column. v e periods of the secular term are g ven in the last the surface of the planet is / = yr2, where /0 is the insolation at 1 AU,3 and r = SE the Sun-Earth distance. When averaged over the year, that is, over the orbital period, we find the average annual insolation As the semimajor axis a is constant, 1M depends only on the eccentricity that varies from 0 to about 0.06 over 10 Ma. The relative variation of lM is thus 1.8 X 10~5 that is very small. By contrast, the ratio of insolation at perihelion versus aphelion is where 7" is the temperature expressed in Kelvin. Considering an average temperature of 7"= 285 K (14.85°C), we obtain t>7"=4.8K for the present eccentricity (*> = 0.0167), and 6T= 17.1 K for e = 0.06. These simple examples are quoted here to emphasize how the eccentricity can modulate the seasonal insolation. For more complete models, one can refer to Paillard (1998, 2001) and Bosnians et al. (2014). for l+e which amounts to 1.27 at maximum eccentricity e = 0 06 Averaging over the Earth surface and using a simple radU alive model, this relation translates into a relative varia tion of temperature of a planet 6T, from perihelion to aphelion, considered as a black body, as 6T — % e 4.2.1 Decomposition of the eccentricity The eccentricity signal is one of the major targets iui stratigraphic studies, especially for older times, before the Neogene Period. It is thus important to understand the main components of the eccentricity signal. The decomposition of this signal in terms of fundamental frequencies is given in Table 4.2. In this decomposition, all ot the terms are recognized as combinations of the fundamental secular frequencies (Table 4.1). More precisely, most terms are differences of two g, except ^ 82 gs- (g4 - g3). indeed, all combinations of frequences in the periodic decomposition of the eccentri- can 3he °f ? f0mi " = with * = » - °- ThiS oe easily understood when one realizes that the 3. Astronomical unit. Astrochronology Chapter | 4 143 TABLE 4.2 First 10 terms (in decreasing amplitude) of the frequency decomposition of the Earth's eccentricity over the time interval [ - 15, +5] Ma. k M(7year) P (kyr) bk X 104 i gl -g'' 3 200 405 107 2 g4~gs 13.652 95 81 3 g4-gj 10.456 124 62 4 gj-gí 13.110 99 53 5 gj-gi 9.910 131 45 6 g4-gl 0.546 2373 30 - gi -gs 1.326 978 28 6 g4-gl 12.325 105 21 9 g2-gS-(g4-gj) 2.665 486 20 10 g2-g> 1.884 688 18 The eccentricity e can be expressed as e = eo + Jil-i o»cos(/itf + r;4) with e0 = 0.0275579. Column two lists the corresponding combination of frequencies where gi are the fundamental frequencies (Table 4.1). Source: Adapted from Laskar et al. (2004). important variable in the dynamical evolution of the solar system is not the eccentricity (e), but the complex variable z = e exp(i:u), where w = ft + w. As shown in Figs. 4.2 and 4.3, = ELi bkexp(igkt + 6k). As this solution is composed of only five periodic terms, the frequency decomposition °f the eccentricity e(5) = lz(5)l is more straightforward and is provided in Table 4.4. In the first column, k is the index TABLE 4.3 The five leading terms in the frequency decomposition of the complex eccentricity variable z = e exp ix for the Earth over the time interval [ - 15, +5] Ma (Laskar et al., 2004) (j = EAexp(W + 0,)). n gk k gk ("/year) bk 8k (degree) 1 gs 5 4.257 0.0189 30.7 2 82 2 7.457 0.0163 - 157.8 3 8< 4 17.910 0.0130 140.6 4 g3 3 17.367 0.0088 -55.9 gl 1 5.579 0.0042 77., J of the term in the frequency decomposition (by decreasing amplitude) of e(5\ and K the rank of the same term in the decomposition of the full Earth eccentricity (Tables 2 and 6 from Laskar et al., 2004). It is important to note that all 10 leading terms of the Earth's eccentricity can be explained by only the first 5 terms of z3. In Section 4.2.3, we will discuss further the outcome of the decomposition of Table 4.4 in the observed aspects of the Earth eccentricity solution and their possible manifestations in the geological data. Before, it is instructive to understand the mathematical origin of the periodic terms involved in Table 4.4. 4.2.2 Mathematical intermezzo Let consider z = ^ake\p(\(gkt + ek)) t-i N Z = £\ztexp(i7i>) (4.1) the expression of z = e exp(i^), where the amplitudes ak are positive real numbers and nk = gkt + Qk. The eccentricity e is then e = Jzl, where z is the complex conjugate of z, that is, N z = ^ a*exp( - irrt). (4.2) *=i We have (4.3) zz = e2= £wa*0/exp(i(7rt - m)) = £* al + a*a/exp(i(7Ti - 717)). We see that the arguments that appear are differences nk-Tr, = (gk-gi)t + tfk-0i) with frequencies gk-g,. 144 PART I II Concepts and Methods (TABLE 4.4 H* ,3 periodic Wn» «n decree »mp« -ic limited to the first five linear terms of z (Table <^___---- sition of e(S) = |z<5)|, v /hen sr® * k Ú("/year)__ 3.200 P (kyr) 405 b'k x io4 109 1 2 1 2 «4-«5 13.653 95 82 3 3 g4~ft 10.453 124 66 4 4 gi~gs 13.110 99 53 5 5 gl-gl 9.910 131 44 6 6 g*-gl 0.543 2387 35 7 7 gl "gs 1.322 980 25 8 10 g2-gl 1.878 690 21 9 8 g4-gl 12.331 105 16 10 12 g2+g4-2g5 16.853 77 16 11 18 gs+g-t-g.'-gs 23.563 55 13 12 9 g2-g5-(g4-g3) 2.657 488 13 13 20 gi-gs + (g*-gi) 3.743 346 13 e*51 = 3> bicosOijf + 0,) with e0 - 0.0269. k is the rank of the term by decreasing amplitude in e'5', while K is the rank of the same term in e (Table 4.2). Column three is the corresponding combination of frequencies, g, are the fundamental frequencies (Table 4.1). involving all fundamental frequencies gt. But this is e2, and not e. With el = ^ka2k, we have e2=e2Q(l+X) (4.4) with x ~ H a*a'exP(K"* - n,)) (4.5) and thus assuming that X is small with respect to 1, and expanding up to second order in X, e = e0^Tx = e0(^+l-X-l-X2 + O(X^. (4.6) Thus X will involve differences of two frequencies ft - ft (Eq. 4.5) as terms 1-9 of Table 4.4, but the terms of the X part will be sums of two differences nk-n, These will be terms of order 4, involving four fundamental frequencies g„ as terms 10-13 of Table 4 4 It is also -nterestrng to note that the phase of these terms will £ opposite to the equivalent combination of 4« because of the minus sign in Eq (4 6) ^«"»»eiits Now we can return to the simple example of = lz I. z has only five periodic components A« H<»m strated earlier, e{5) will contain «„i„V em°n" order, sum of ^^C'^T^^T Indeed, this is well revealed bvT> *' ,(Table 4"4)' "» (Fig. 4.4). For £^J£F£* g** <* are easily identified and is thus remarkable that the most important features of the eccentricity solution of the Earth are provided by the simple model zi5) (Table 4.4). It should be noted that the largest periodic component of the eccentricity is the 405-kyr term g2 - g5. This term is fundamental in cyclostratigraphy as its period is very stable and can thus be used as a metronome for the establishment of local and global time scales (Olsen, 1986; Laskar, 1999; Laskar et al., 2004, 2011a; Boulila et al., 2008; Hinnov and Hilgen, 2012; Kent et al., 2018; Hinnov, 2018a; Huang, 2018). 4.2.3 Eccentricity modulations Due to the importance of the 405-kyr mode (g2 - ft). & is important to filter the data to retrieve its 405-kyr component. From Fig. 4.4, it is clear that the g2 - g5 mode does not occur in isolation but is surrounded by two nearby peaks, řr^f"8,10 82 ~ 85 ~ <* " ft) ■*» ft " ft + ^ ' 83) ^reps. 488 and 346 kyr period). These side terms produce a modulation of the 405-kyr component eh with frequency «4 ft (Fig. 4.5B). As the g4 - g3 term also appears in the eccentricity, the g4 - ft mode can also be directly retrieved oy filtering the eccentricity in the [0:l.l]7year interval (ea) ( 'g- 4.5A). By superposing ea with the envelope h of the ^>Kyr component eb> one can see that ea is almost identical, although in opposite phase to eb (Fig. 4.5) (see also Laskar et al., 201 la). (~imif478i modulation appears also in the high frequency eccentricity terms. These main terms appear Astrochronology Chapter | 4 145 a E aj, o a E o -4.5 10 15 20 25 freq ("/yr) FIGURE 4.4 Fast fourier transform (FFT) of the La2004 eccentricity solution over 33 Myr (top) and Fourier transform of the solution limited to the five main linear terms of ; (Table 4.3). For el5). all the terms can be easily identified, and a combination of their corresponding frequencies are reported in the figure (see also Table 4.4). The periods of the corresponding terms are displayed (in kyr) in the top figure. Frequencies are expressed in arcsec/year ("/year): 1 "/year = 0.7716 cycle/Myr. (Fig. 4.4) in two sets: (g3-g2, g4-&) and 0j3~ ft. g4 - g5). These components will both modulate with a g4 - g3 frequency (Fig. 4.5C and D). In this case, the modulation envelopes ec, ej are similar and in phase with g4 - g3 (ea). We can explain this using our simple 5 term model. Let us consider the filtered eccentricity in the [9.3, ll]7year band. We will have ec = a exp(i(7r3 - 7r2)) + a'exp(i(7r4 ~~ ^2))» (4-7) where a, a' are both positive (Eq. 4.6). Thus ec = [a + a'exp(i(7r4 - 773))] exp(i(7r3 - 7r2)), (4.8) and as a' is positive, the slow modulation a' exp(i (7r4 - 7r3) appears with the same phase as ea. This is the same for ed (Fig. 4.5D) and for all order 2 couples 8i~gp g*-gj, as for example, (g3~8i< 8*~8\< Fig- 4.4). Now let us consider the modulation of the 405-kyr term, eb. This term involves three components in its simple approximation (Eq. 4.6) Cfc = a exp (i(7r2 - tt5)) - b exp(i((7T2 - 7Ts) _ (*4 ~ K^ ~ f exp(i((7T2 - 7T5) + (7T4 - 7T3))) = exp(i(7r2-7r5)) X [a - b exp( - i(7r4 - 7r3)) - b' exp(i(^4 - 7r3)))J. (4.9) where a, b. and b' are positive. We have now a minus sign before b and b' (- \/&X2 in Eq. 4.6). This induces a modulation of eb with frequency g4 — g3, but due to these minus signs, it will be in opposite phase with respect to ea. Indeed, if instead of the eccentricity, expanded as 1 -1- 1/2X - 1/8X2 (Eq. 4.6), we consider a fictitious eccentricity like expression 1 + 1/ZX + 1/8X2. with the opposite sign in the terms X2 of fourth order,4 then the modulation in the 405-kyr band of this fictitious eccentricity is in phase with g4 - g3. 4.3 Chaos in the solar system Since the first semianalytical long-term solutions of Laskar (1988, 1990) and Laskar et al. (1993), it becomes possible to compute reliable orbital solutions starting from the present initial conditions (Section 4.1.1). This was confirmed later on by direct numerical integrations (Quinn et al., 1991; Laskar et al., 1992, 2004). It was previously thought that the progress of computers and of observational techniques would result in an astronomical solution with higher precision, so that time validity could be extended steadily both in the future and in the past as envisioned by Laplace (1812). But the discovery of the chaoticity of the orbital motion of the solar system put an end to this hope (Laskar, 1989, 1990). 4- X is of second order in the gt (Eq. 4.5). 146 PART | II Concepts and Methods 0.08 0.06 0.04 0.02 -0.02 -0.04 -0.06 0 10 15 20 25 time (Ma) FIGURE 4.5 Filtered eccentricity of the La2010a solution (Laskar et al.. 201 la). The filtered solutions are shifted in order to be plotted on the same graph. (A) ea is the filtered eccentricity in the [0. 1.1]"/year (period > 1.18 Myr) band (+0.03) (red). (B) eb is filtered in the [2.2. 4.3]7year ([301. 589] kyr period) band {purple). (C) ec is filtered in the [9.3, 11 ]7year ([ 139, 118] kyr period) band (-0.03) (green). (D) ed is filtered in the [12.6, 14.5]7year ([103, 89] kyr period) band (-0.06) (blue). The upper envelopes of eb. ec. ed. respectively, eb, ec. id are plotted in red. The thin black curve is the (ea) curvr shifted in order to compare to the envelopes eb, er. ed of fj„ er, ed. (e„) nearly coincide with er. ej and is phase opposite to the eb. See the text for discussion. See also Laskar a ul.. 2011a. Indeed, the uncertainty in the solutions grows exponentially, by a factor of 10 every 10 Myr (Laskar, 1989). More recently, it was shown that the motion of the minor planets Ceres and Vesta is itself chaotic, on much shorter time scales than the planets. Due to the perturbation of these celestial bodies on the planets, the possibility for constructing a precise orbital solution for the planets of the solar system from their present state is limited to about 60 Myr (Laskar et al., 201 lb). Thus the use of the Earth's eccentricity solution as a template for cyclostratigraphy will suffer the same limitation. In Fig. 4.6, five eccentricity solutions are plotted over 60 Myr into the past (only 2 Myr slices are plotted every 10 Myr). The La2004 model (Laskar et al 2004) is widely used and has been demonstrated to be precise for the past 40 Myr (Laskar et al., 201 la). The La2010a and La2010d solutions are improved versions, with a model including the five major asteroids (Ceres, Vesta, Pallas Iris' and Bamberga). Their initial conditions were obtained using a fit to a 1 Myr long high-precision planetary ephemeris rNPOP (Laskar et al.. 2011a). La2010a is adjusted to INPOP08 (Fienga et al.. 2009) and La2010d to INPOP06 (Fienga et al., 2008). As it was realized that INPOP06 is more accurate than INPOP08 (Fienga et al., 2011), La2010d should be preferred to U2010a that is in agreement with the comparison to the updated version La2011 (Laskar et al., 2011b) that is adjusted to INPOPlOa (Fienga et al., 2011). This was also confirmed through comparison with geological data (e.g., Boulila et al., 2012; Westerhold et al., 2012). 4.3.1 Drifting frequencies Another expression of this chaotic motion is the fact that the main frequencies of the system (Table 4.1) are not constant but can drift in a significant way (Laskar, 1990; Laskar et al., 2004), even if the system is largely conservative, with minor dissipation.5 These variations are summarized in col Aioo of Table 4.1 that represents the variation of the different fundamental frequencies observed over 100 Myr. As was already described in Laskar (1990), these variations depend largely on the involved planets. Indeed, the chaos is not evenly distributed among the planets. The frequencies related to the outer solar system (g5, g6, g7, gg, s6, s7, s8) are nearly constant over the age of the solar system and reflect the mostly regular behavior of the outer solar system (Jupiter, Saturn, Uranus, Neptune).6 By contrast, the frequencies related to the inner planets, (gi, g2, gj, g4, Si, $2. St) undergo significant variations, with some differences in their unstability. They can be put in three classes, depending their A,00 value (Table 4.1): 1. unstable frequencies: g(, g3, g4, jlt s2 2. moderately unstable frequencies: s3, s4 3. nearly stable frequencies: g2 This last frequency is of particular interest as it contributes to the g2 - g5 term with a 405-kyr period that is the largest term of the eccentricity signal (Table 4.2). Despite the chaotic motion of the solar system, this term can thus be used as a metronome for the time calibration of the stratigraphic record in the Mesozoic and beyond. 4.3.2 The 405-kyr g2 -g5 metronome The main periodic component of the Earth's eccentricity is the 405-kyr g2 - g5 term (Table 4.2). The value of gi is practically constant, and g2 presents only small chaotic diffusion (Table 4.1). This component can thus be approximated by a single periodic term that gives an approximate 5. This is not the case for the rotational motion of the Earth which ' w '- 6. As a rule of thumb, one can consider mat 17yr corresponds to a'rToTn T Eanh-M°°" system, of 2* after 1 Gyr. More precisely, 17yr = 0.77 X 10"6 cycles/year. (U% M>T MaCtly)- A variation of 0.0017yr will make an offset La2004 La2010d La2010 a La2011 PART I II Concepts and Methods eccentricity, including the constant term, expressed in a very simple form Laskar et al. (2004) ew5 = 0.027558 - 0.010739 cos(2434" + 3.2007). (4.10) This expression was established by a fit to La2004. but with the improved solutions of Laskar et al. (201 la), it appears that there was no need to change this formulation (Fig 4.7). One needs indeed to remember that beyond 60 Ma, as it is obvious from Fig. 4.7, the drift in frequencies becomes apparent and cannot be predicted only by the celestial mechanics computation. But this unknown drift is small and amounts to less than one period over 250 Myr (Fig. 4.7), which is about 405-kyr over 250 Myr (~ 1.6%.). This is better than most radioisotopic determinations (e.g.. Fig. 1.4 from Gradstein et al., 2012). Eq. (4.10) can thus be used for cyclostratigraphic tuning over the whole Mesozoic and beyond. The stability of this 405-kyr term was recently confirmed by precise U-Pb zircon dates at 210-215 Ma (Kent et al., 2018). In an equivalent way, one can use the following formula, expressed in radians t-405 = 0.027558 - 0.010739 cos 0.0118 + 2^ ^O.O 405.000; (4.11) where t is in years and counted negatively in the past. In Fig. 4.8, ^405 is plotted on selected time intervals over 250 Myr. It is compared with the filtered eccentricity in the [2.2, 4.3]7year ([301, 589] kyr period) band for four recent solutions La2004 (Laskar et al., 2004), La2010a, La2010d (Laskar et al., 2011a), and La2011 (Laskar et al., 2011b). It should be noted that these filtered solutions, as in Fig. 4.5B, include the side terms that induce a g4 - gi modulation of the g2 - g5 component, which is why, even in the most recent time, the amplitude of the filtered eccentricity does not strictly match the purely periodic 6405 solution (Eq. 4.10). Beyond 55 Ma, there is also some phase shift, but this is expected, due to the uncertainty of the behavior of the g2-gs mode beyond 60 Ma (Fig. 4.7). Even at 250 Ma, the phase shift is less than half a period, below the above-quoted ~ 1 ft. uncertainty. Warning: For stratigraphic calibration purposes, it is in general not recommended to use the filtered eccentricity inn5 L°nd 40Ma f°r 132004 "* 50M. for La2010 La2011, as beyond this age, their behavior is not consistent Moreover by tuning to this filtered eccentricity one S SaSS£5#=s introduced in the tunine S ^ conWnt a pure cosine IcZ^ ^TT"** l° USe *ere is no reason to ^ 3^ frequency (405,000 year period) that ^ radioisotopic daUng (Lnt e?^ 0t8 Z^"^ by «"•• -wis). 11 the improvement 100 150 time (Ma) FIGURE 4.7 Differences (in radians) of the argument 6^-^ 0f g2-gi in all solutions La2004 (Laskar et al.. 2004). La2010a,b,c,d (Laskar et al., 201 la) with respect to the pure single-frequency approxi-mation ftwsC) = 3.200"r. where r is in year. Adapted from Laskar et al. (2011a). La2004 ■ • La2010d La2010a La2011 247 248 time (Ma) time (Ma) ÜIC,UÜL4,8 405-kyr & - «5 exponent In black is plotted 4«. * angle-frequency periodic term provided by Eq (4.10). The colored curves correspond to four different eccentricity solutions U2004 (Laskar Bt aL »M£ La20,0a, U2010d (Laskar et al.. 2011a). and U20H (Laskar^ M lb). Each solution is filtered in the [2.2, 4.3]7year ([301. 589] kyr ^ band and thus mcludes the modulation by the g[ -g3 component (Fig. ^ Astrochronology Chapter | 4 149 of radioisotopic measures provides more precise time constraints in the future, then it will be possible to improve a tuning target by providing either a slightly different value for the fc-405 frequency or even a varying frequency for this term (Fu and Laskar, 2019). Such improvement is more than welcome, but meanwhile, one should stick to the constant 3.2007year frequency. By contrast, for ages that are within the validity time of the solution, that is 40 Ma for La2004 and 50 Ma for La2010 and La2011, one can use the full eccentricity solution, as well as the derived filtered eccentricity (Fig. 4.5B). 4.3.3 The g4 -g3 2.4 Myr cycle The g2 fundamental frequency is the most stable, not considering the outer planet ones. This led to the recognition of the g2 - gs metronome. By contrast, g3 and g4 are the most unstable frequencies (Table 4.1). Moreover, we have seen the important role of the g4 - g3 2.4 Myr term in the eccentricity (Section 4.2.3). g4-g3 is the sixth term in amplitude in the eccentricity (Table 4.2) but appears also as the main modulation of the g2 - g5 405-kyr term and also as the modulation of the ~ 100-kyr terms in the eccentricity (Fig. 4.5). However, this term cannot be used for time calibration, as its behavior is not stable, and its frequency, as for s4 - j3, will evolve because of the chaotic diffusion of the orbits (Fig. 4.9). This modulation has been recognized in sedimentary records of the Cenozoic and Mesozoic eras (Olsen and Kent, 1999; Palike et al., 2004; Boulila et al., 2014; Fang et al., 2015; Ma et al., 2017; Westerhold et al., 2017), although in Olsen and Kent (1999), the 405-kyr modulation was measured with a period of about 1.7 Myr, instead of the present 2.4 Myr value. The question arises as to whether this difference could be the expression of the chaotic diffusion of the solar system, and this was answered positively in Olsen et al. (2019). Indeed, in Fig. 4.10, extracted from Olsen et al. (2019), the period of the g4 - g3 argument is plotted versus time for 13 different orbital solutions. For the most recent 40 Myr, they all reveal the same ~ 2.4 Myr period, but then they depart from each other due to chaotic diffusion (Laskar, 1990; Laskar et al., 2004). The green horizontal line represents the 1.7 Myr value observed in the Newark-Hartford data (Olsen and Kent, 1999; Olsen et al., 2019). This value is attained by many of the solutions and in particular by La2010d (in black) at roughly the same 200 Ma age. It can also be observed that the excursion of the PgA-gi period is even larger and can evolve across the [1.4:2.6 Myr] period range during this time interval. The prediction of the evolution of the actual path of the Pgt_gi period in the past cannot be retrieved by only considering the present planetary positions and computing their past orbits using the laws of celestial mechanics. As 40 60 time (Ma) FIGURE 4.9 Top: differences (in radians) of the argument of g^-gi in solutions La20O4 (Laskar et al.. 2004), La2010a,b,c.d (Laskar et al.. 2011a) with respect to the linear evolution 2.6647", where T is in Myr. Bottom: differences (in radians) of the s4 - sy5 argument in La2004, La2010a.b,c,d with respect to the linear expression 2 x 2.6647", where 7" is in Myr. Adapted from Laskar et al., 201 la. 100 time (Ma) FIGURE 4.10 Evolution of the period of the g4- g} argument for 13 orbital solutions over 250 Myr in the past. The horizontal line is the 1.7 Myr value observed in the Newark-Hartford data. The red curve is La2004. and the black curve La2010d. Over the first 40 Myr, all values are of ~ 2.4 Myr, but they diverge after 50 Myr due to chaotic diffusion. La2010d (black) has nearly the same value as the one found in the Newark-Hartford data around the same age (200-220 Ma) (Olsen et al.. 2019). in Olsen et al. (2019), we will have to rely on geological data to retrieve this information. Recovering these long-period cycles in the geological data is in some sense recovering the planetary orbital motions through geological data beyond their horizon of predictability. 150 PART | II Concepts and Methods 4.4 Inclination and obliquity The shape of the Earth's orbit, regulated by the eccentricity, is not the only important parameter for the computation of the insolation on the Earth's surface. The other main ingredient is the orientation of the Earth's spin axis that is regulated by the obliquity e, the angle between the orbital plane of the Earth and its equator, and the precession angle. 0, that describes the orientation of the spin angle in its slow motion around the pole of the orbital plane, n (Fig. 4.3). Here we make the approximation that the spin axis is also the axis of inertia of the Earth. The precession ip and obliquity e (Fig. 4.3) equations for the rigid Earth in the presence of planetary perturbations are given by Kinoshita (1977), Laskar (1986), Laskar et al. (1993), Neron De Surgy and Laskar (1997), and Laskar et al. (2004) 1 - ^(2(/)sin^- A(t)cosip) ^(Mt) sin^' + 3(1) cos0) - 2C(t) dX -i *'\ chp = aX dt ~ L (4.12) with8 X = L cos e, L = C7, where 7 is the spin rate of the Earth, A< B < C are the principal momentum of inertia of the Earth, and s/\-p2-q2 [q+p{qp-pq)) g(0= , , = V*-P2~q2 Qt) = qp-pq \p-qiqp-pq)] (4'13) where q = sin(i/2) cos fi and p = sin(tf2) sin fi, and where q is the precession constant 3G 1 21 niQ Mm ( is constant) and no planetary perturbations. The elliptical elements are thus constant, and A = 2 = C = 0 in Eq. (4.13). Eq. (4.12) reduces then to = 0 dcose dt dip — = a cos eo i.e. e = £-0 = Cte. (4.15) The obliquity is then constant, and the precession angle ip evolves linearly with time at a constant angular speed of a cos e0. This is a zero order solution. We can go further by reducing (Eq. 4.12) to the first order terms. We obtain the solution of order one, de jt - lip sin ti> - q cos t/>) = 22te(C exp(iV)). (4- ,6> where C = sin(i2) exp(ifi) and %e denotes the real part of the complex number. With the quasiperiodic approximation (e.g., Table 4.5 of Laskar et al., 2004) N C- 2_,a*exp(i(^r + 04)), *=i (4.17) 7. The angle between the Earth's spin axis and its axis of inertia is less ,h„n I • 8. There ,s here a misprint in Uskar e, al. (2004, It f. and not X = cos, Astrochronology Chapter | 4 151 The first-order solution of the obliquity will be a similar quasiperiodic function e = e0 + 2}_^ —— cos((y* + p)t + k + ^ {4 ,8) k=\ k & The terms that appear in the obliquity have thus frequency uk + p, where p is the precession frequency, and uk die the frequencies of the inclination variables = sin (i/2) exp(/fi*) (Fig. 4.3). The amplitude of these terms is multiplied by Vi/{yk+p). High frequencies are thus favored (factor uk). Amplitudes are also divided by vk + p, and resonance will occur when i/k+p = 0. At present, p = 50.4758387year (Laskar et al., 2004), but due to tidal dissipation in the Earth-Moon system, p is not constant but evolves in time, as the spin rate of the Earth and the Earth-Moon distance evolves. 4.4.2 Tidal evolution The Lunar Laser ranging measurements have taken place since the Apollo and Lunokhod mission installed reflectors on the Moon nearly 50 years ago, with an accuracy that is now less than 2 cm (e.g., Viswanathan et al., 2018). This allows us to monitor the present recession of the Moon, at a rate of ~3.8 cm/year (Dickey et al., 1994; Laskar et al., 2004). Backward integration of the Earth-Moon system provides interpolation formulae for the Earth-Moon distance (aM, in Earth radius), the length of day (LOD. in hours), and the precession constant (p, in "/year) as provided in the La2004 solution (Laskar et al., 2004) aM = 60.142611 + 6.1008877 - 2.709407r2 + 1.366779F3 - 1.48406274 LOD = 23.934468 + 7.4321671 - 0.7270467"2 + 0.409572r' - 0.58969274 p = 50.475838 - 26.3685837 + 21.89086272 (4.19) where 7* is the time from the present (J2000), expressed in Gyr and counted negatively in the past (Fig. 4.12). These expressions have been established by a fit over 250 Myr but can be extrapolated over 500 Myr for a first estimate of the past evolution of these quantities. It should nevertheless be reminded that these expressions cannot be extrapolated over the age of the solar system, and the past evolution of the Earth-Moon system is still largely unknown. If one integrates back the evolution of the Earth-Moon system, owing to the present rheology parameters of the Earth, one finds that the Moon hits the Earth at about 1.5 Gyr ago, which is clearly not compati-b'e with our understanding of the origin of the Moon or J«tory of the Earth (Gerstenkorn, 1969; Walker and Zahnle, 1986). In order to reconcile this evolution with 35 freq 7yr FIGURE 4.11 Spectral analysis of the obliquity ;. The spectral analysis is performed over the interval (10:20) Ma. The main peaks are recognized as p + sh where s, is one of the fundamental frequencies of the inclination of the orbital plane (Table 4.1). On top, the periods are given in kyr. Two additional terms of higher order are given: p + s} - (gt - g3) and p + j4 - (g4 - g3) (see Table 4.5). Frequencies are expressed in arcsec/yr ("/year): 17 year = 0.7716 cycle/Myr. FIGURE 4.12 Past evolution of the Earth-Moon distance aM (top. in Earth radius RE). of the LOD (middle, in hours), and precession frequency p (bottom, in arcsec/year). These curves are obtained using Eq (4 19). which are extrapolated from the La2004 solution over 250 Myr (Laskar et al.. 2004). LOD. Length of day. PART | II Concepts and Methods 24 23 22 £ 21 Q O J 20 19 18 17 Bvalva ♦ Cart • Strom«»»« * CycWMttW1)' ' 500 1000 1500 time (Ma) 2000 2500 FIGURE 4.13 Length of day (LOD) evolution due to tidal dissipation in the Earth-Moon system. The dotted red line is the LOD provided by Eq. (4.19) (Laskar et al.. 2004). The dotted black line is an empirical fit using a simplified tidal model adjusted to the geological data (Walker and Zahnle, 1986; Lambeck. 1980; Berger and Loutre. 1994). top. Length of day. Compilation of various data from I Williams, 2000) and references therein. The cyclostratigraphic data are from (Meyers and Malinvemo. 2018). the age of the Moon, one needs to assume that the present tidal dissipation of the Earth is about three times its past averaged value. This is possible, as the present tidal quality factor 0 (~ 11) is largely due to the dissipation in the shallow seas and thus subject to change by a large amount with the repartition of the continents (by comparison, for Mars, 0~9O). Moreover, the tidal response of the oceans strongly depends on the rotation period of the Earth, and resonances may occur that increase the tidal dissipation (Webb, 1980, 1982; Auclair-Desrotour et al., 2018). But the precise past evolution of the Earth-Moon system will require some input from the geological record. There are numerous estimates of the past rotational state of the Earth, obtained from various indicators such as bivalves, corals, stromatolites, or tidal deposits. These records have been compiled in several publications (e.g., Lambeck, 1980; Williams, 2000) (Fig. 4.13). It should nevertheless be stressed that most of these data suffer from large uncertainties that are not always estimated. It is certainly needed that these data or other equivalent data are reanalyzed using clear, updated, methodologies, and procedures. All raw data should further be made publicly available. Tidal dissipation is also expressed in geological records by the shortening of the climatic precession and obliquity periods back in time (Eq. 4.19) (see also Berger et al, 1992; Berger and Loutre, 1994). These climate forcing terms have been recorded in sedimentary geological archives and associated datasets (e.g., Zeeden et al., 2014; Meyers and Malinvemo, 2018) (Fig 4 13) While this tidal dissipation effect can be seen as a phase shift of the precession/obliquity cycle relative to a solu tion assuming recent tidal dissipation rate in the Neogene (Lourens et al.. 2001; Zeeden et al.. 2014), a shortening of the precession and obliquity periods relative t0 £ stable eccentricity 405-kyr metronome is observed j„ Paleozoic and Mesozoic datasets (e.g., Wu et al., 2013a-Boulila et al., 2014. 2019). Such datasets from ^ Mesozoic and Cenozoic could be used to reconstruct the Earth's precession and obliquity periods in a quantitative manner, and it is desirable that the analysis of such records will be continued in order to improve the knowledge of the past evolution of the Earth-Moon system. In Fig. 4.13 is also plotted (in red dashed line) the computed variation of the LOD as obtained by Eq. (4.19) (Laskar et al., 2004). It should be stressed that this curve has not been fitted to the available geological data (Fig. 4.13) but is obtained through the sole use of the Lunar Laser ranging data over the past few decades. In addition to the variations expressed in Eq. (4.19), the tidal dissipation induces an average variation in the obliquity itself which can be written as £ = 23.270773 + 2.0112957" (4.20) where 7" is in billions of years (Laskar et al., 2004), counted negatively in the past. The obliquity was thus smaller going back in time (see Fig. 14 from Laskar et al., 2004). This formula, obtained through a fit over 250 Myr, could also be used over 500 Myr in the past, although as stated earlier, large uncertainties remain, which can only be improved by constraints provided by the geological record. In addition to the tidal dissipation in the Earth-Moon system, the variations of the Earth's spin rate and orientation can result from changes in the momentum of inertia of the Earth. These can result from change in the ice bold (e.g., Laskar et al., 1993; Levrard and Laskar, 2003) or plate tectonics (e.g., Mitrovica et al., 1997; Morrow etal., 2012). The problem with these effects is that their signature is not easy to disentangle from that of tidal dissipation, as they will also manifest themselves by a change in the precession rate (e.g., Palike and Shackleton. 2000; Lourens et al., 2001). Over Gyr time scales, it may further be necessary to take into account the mass loss of the sun that will affect also the orbital secular frequencies (e.g-Spalding etal., 2018). 4.4.3 Obliquity solution Due to the dissipation in the Earth-Moon system described earlier, the analysis of the obliquity solution is complex. It >s nevertheless interesting to look to the main features of the solution over a limited time of 20 Myr, where the dissipate aspect is moderate (Fig. 4.14). In Fig. 4.14 the obliquity t plotted, as well as various filtered expressions et, e* tered over respectively [28:38], [23:38], [42:47]'7year. These filtering intervals are dictated by the analysis of the Founf spectrum of the obliquity (Fig. 4.11). The envelopes Si, • Astrochronology Chapter | 4 153 0 5 10 15 20 time (Ma) FIGURE 4.14 Obliquity (e) evolution (in degrees) over 20 Myr from U2004 (Laskar et al.. 2004) (top). C] is the filtered obliquity in the window [28:38]7year ([34.1:46.3] kyr periods) (green). In red is plotted the envelope €\ of e\. £2 is the filtered obliquity in the wider window [23:38]7year ([34.1:56.3] kyr periods) (green). In red is plotted the envelope £2 of s2- £3 is the filtered obliquity in the window [42:47]7year ([27.6:30.9] kyr periods) (green). In red is plotted the envelope i\ of e\. The vertical scale is the same for e, e 1, e2 and five times larger for £3. and £3 of these filtered obliquity solutions allow to extract the most important components of the obliquity. In Fig. 4.15 are plotted the Fast Fourier Transform (FFT) analyses of these envelopes §1, e2, s3, of the filtered obliquities eu e2, £3. with the identification of the main terms. As expected, the main term in these envelopes is related to the s4 - 53 term, with a period of ~ 1.2 Myr. This term results from the beat of the p + s4 and p + s3 obliquity terms (Fig. 4.11 and Table 4.5). However, other terms appear as well. The term g4 - g3 is also present in the eccentricity solution with a period of ~2.4 Myr. This term results from both the beat of the p + s3 and p + s3 - (g4 - 83) terms and the P + «4 and p + 54 - (£4 - g3) terms (Fig. 4.11 and Table 4.5). Very important, is further the S3 - s6 term, appearing as the °eat of p + s6 with the main obliquity term p + s3. Finally, si - s2 appears as the beat of p + Ji and p + s2. The important feature of all these spectral terms is that they do not depend on the precession frequency p. but only on the orbital solution with secular main frequencies gi, Sj. These terms will thus not be affected by the strong variations in p (Eq. 4.19 and Fig. 4.12). Both s4 - s3 and s3 - s6 are of particular importance: the first one because it is at present in resonance with the modulation frequency of the eccentricity g4 - g3 (s4 -s3 = 2(g4 - g3)), and the second one because s6 is a stable frequency and s3 a moderately stable frequency (Section 4.3.1). It is thus possible to use the s3 - s6 inclination term as an additional chronometer for stratigraphic tuning, with a period of 173 kyr. 4.4.4 The 173-kyr s3 - s6 metronome The g2 ~ g$ 405-kyr metronome is a fundamental tool for establishing local or global time scales (see Section 4.3.2), but this signal is not always present. Recently, it has been demonstrated that in some cases the s3 - 56 173-kyr cycle can also be used as a metronome for the calibration of stratigraphic sequences (Boulila et al., 2018; Charbonnier et al., 2018). This cycle allows to calibrate obliquity dominated stratigraphic sequences. This 53 - 56 term, present in the modulation of the obliquity (Figs. 4.14 and 4.15B), does not depend on the precession frequency p and is quite stable in time (Fig. 4.16). Only the variation of the orbital plane of the Earth is involved. We can call this term the 173-kyr inclination metronome, analogous to the 405-kyr g2 - g5 eccentricity metronome. The time scale uncertainty associated with the inclination metronome is of the order of 400 kyr over 100 Myr that is about 0.49c. But contrary to the eccentricity metronome, the inclination metronome is not the largest term present in the obliquity and not even in the modulation of the obliquity. It is nevertheless quite isolated (Fig. 4.15B) which explains why it can be successfully used for stratigraphic calibration (Boulila et al., 2018; Charbonnier et al., 2018). A good approximation for this cycle can be given by the following expression ^3-^(0 = 0.144 cos(404,444" + 7.5"r) (4.21) where t is in years, counted negatively in the past. The angle is in arcseconds and should usually be converted to radians to compute the cosine. The frequency s3 - s6 has been rounded to 7.5 "/year, as it is meaningless to use the exact expression s3 - s6 = 7.4978557year obtained from Table 4.1 due to the variability of s3. Alternatively, one can use the same quantity expressed in radians and years (counted negatively in the past). £,3^(0 = 0.144 cos(l.961+2rrr^5). (4.22) 154 PART I II Concepts and Methods £0.01 4.5 Chaotic diffusion and secular resonances The present solar system is characterized by the presence of two main secular resonances (Laskar, 1990, 1992; Laskar et al., 1992, 2004, 201 la). This is expressed by a commen-surability relation among the secular main frequencies while the corresponding angular argument is oscillating (We ^ i( is in libration, like for the small oscillations of a pendulum) and not circulating (like a rigid pendulum with large initial velocity). These two resonances are 6 = 2(g4 ~gi) ~ (S4-S3) (4.23) and e = (g\ - ^5) - (*i -s2). "freq( /yr)6 FIGURE 4.15 FFT of the envelopes e,. h. h of the filtered obliquities e,, £2, E) from Fig. 4.14. Frequencies are expressed in arcsec/yr ("/year)-1 "/year = 0.7716 cycle/Myr. (4.24) Both are important in the dynamics of the system, but the first one draws particular attention as we have seen that the 2.4 Myr g4 - g3 term is the main long-term modulation of the eccentricity (Section 4.3.3). In the same way the 1.2 Myr s4 - s3 term is the largest modulation term of the obliquity (Fig. 4.14). These long-period cycles have been recognized in the geological record (e.g., Olsen and Kent, 1999; Shackleton et al., 2000; Zachos et al., 2001; Palikeetal., 2001,2004). The argument ip9 of 6 = 2{g4 - g3) - (s4 - s3) is in libration in all recent solutions up to nearly 50 Ma (Fig. 4.17), which seems to be consistent with the geological record (e.g., Palike et al., 2004). But over longer time intervals, it is most probable that departure from the 2(g4 - g3) - (s4 - sj) occurs, as what is observed in the numerical simulations (Fig. 4.18). It should be noted that observing a change in the pg4-g> period only is not sufficient to conclude that the system exit the 0 resonance, as the two Pgi-gi and PSt-S} periods can change, but stay in the same 2:1 ratio, corresponding to the black line of Fig. 4.18. In the recent years, there has been an increasing interest search of chaotic transition in the 9 = 2(g4 - g3) ~ ~ ft) 40 60 time (Ma) FIGURE 4.16 Stability of the t - , - of the argument i/co.pL d tc^H^T *» varia-soluuons U20O4 (Laskar ^ «™ for the four La2011 (Laskar « al., 2011b) and! ™ °10d,(Uskar « «L 2011a, made using initial conditions derived from?*, ^21m has ^ 2009; Laskar e. al.. 2004). ™ * & 10 DE421 (Folkner et al., 40 60 time (Ma) FIGURE 4.17 Evolution (in radians) of the argument FIGURE 4.18 Evolution of the period of the s4 - S3 argument versus the period Pgi-g, of g4 - g3 for 13 orbital solutions over 250 Myr in the past. The vertical line is the 1.7 Myr value observed in the Newark-Hartford data. The red curve is La2004 and the black curve is La2010d. The black line corresponds to the 2(g4 - g3) - Cs4 - s}) resonance. The red line corresponds to the (gi - gi) - (st - sj) resonance. The green dot is the origin of all solutions, corresponding to the present date, where all solutions start in the 2 (gi - g3) - (j4 - s3) resonance. Adapted from Olsen et al. (2019). secular resonance (e.g., Grippo et al., 2004; Huang et al.. 2010; Wu et al., 2013b; Ikeda and Tada, 2014; Fang et al., 2015; Ma et al., 2017; Gambacorta et al., 2018; Ma et al., 2019). This search is difficult, as it requires very long records of high quality that are not very numerous. Some convincing results are nevertheless obtained (e.g.. Ma et al., 2017), and we can expect that more will follow in the near future. 4.6 Discussion Since GTS20O4 (Gradstein et al., 2004) and the astronomical calibration of the Neogene (Lourens et al., 2004), huge progress has been made in the analysis of stratigraphic records, and the astronomical solutions are challenged to follow this evolution. Starting from the present initial conditions, despite a highly accurate fit to the most precise observational data, gathered from spacecraft orbiting around the planets, the astronomical solution is limited to 60 Ma (Laskar et al„ 2011b) because of its chaotic behavior. Meanwhile, recent solutions are valid over about 50 Ma (Laskar et al., 2011a). This is not sufficient to address the needs for stratigraphic studies that have covered the Cenozoic and are now being extended to cover the entire Mesozoic. This extension, beyond the 60 Ma limit, is made Possible by the use of both the 405-kyr g2 - gs eccentricity metronome and the 173-kyr s3-s6 inclination metronome (see Sections 4.3.2 and 4.4.4). In order to go beyond the use of these pure periodic terms, it will be necessary to extend the astronomical solutions, and this will only be made possible by using the geological record as an input for constraining the astronomical solution. Encouraging results have been obtained in this direction (Olsen et al., 2019; Zeebe and Lourens, 2019). In the same way the stratigraphic record can be used to constrain the past rotational evolution of the Earth (e.g., Meyers and Malinvemo, 2018), and it is most probable that similar studies will help to decipher the past tidal evolution of the Earth-Moon system in the near future. The search for chaotic transitions in the 2(g4 - g3) -(s4 - s3) secular resonance is a hunt that is shared by many, as well as analysis of other very long periodic components. But in order to obtain convincing results, the stratigraphic community needs to adopt rigorous methods with open shared data, processing techniques, and protocols. It will be the price to switch from qualitative analysis to quantitative results that can be cross compared and used as input for the next generation of astronomical solutions. Acknowledgments The author thanks L. Hinnov and F. Hilgen for helpful discussions and remarks that helped to improve this chapter. The support from ANR-AstroMeso is acknowledged. References Agassiz, L.. 1840, Etudes Sur Les Glaciers. Neuchatel. Auclair-Desrotour, P.. Malhis. S.. Laskar. J., and Leconte. J.. 2018, Oceanic tides from Earth-like to ocean planets. Astronomy & Astrophysics, 615: A23. Berger. A., 1978, Long-term variations of caloric insolation resulting from the Earth's orbital elements. Quaternary Research. 9 (2). 139-167. Berger. A., and Loutre. M., 1994, Astronomical forcing through geological time. Special Publication of the International Association of Sedimentologists, 19: 15-24. Berger. A.. Loutre. M.F.. and Laskar. J.. 1992. Stability of the astronomical frequencies over the Earth's history for paleoclimate studies. Science. 255 (5044), 560-566. Bosnians, J.H.C., Drijfhout, S.S., Tuenter. E., Hilgen. F.J.. and Lourens, L.J., 2014. Response of the North African summer monsoon to precession and obliquity forcings in the EC-Earth GCM. Climate Dynamics. 44 (1-2). 279-297. Boulila. S.. Galbrun. B., Hinnov. L.A.. and Collin, P.-Y., 2008. Orbital calibration of the Early Kimmeridgian (southeastern France): implications for geochronology and sequence stratigraphy. Terra Nova. 20 (6), 455-462. Boulila, S., Galbrun, B.. Laskar. J., and Paelike. H.. 2012. A similar to 9 myr cycle in Cenozoic delta C-13 record and long-term orbital eccentricity modulation: is there a link? Earth and Planetary Science Utters, 317: 273-281. WOS:000301616700027. Boulila, S.. Galbrun. B., Huret, E., Hinnov, L.A.. Rouget, I., Gardin, S., et al., 2014, Astronomical calibration of the Toarcian PART I II Concepts and Methods early Toarcian OAE 98-111 fi,rtft and Planetary Science Utters, .en. middle Eocene astronomical timescale. Earth and Planetary Science Utters. 486 (15). 94-107. . Boulila. S.. Galbru, B.. Sadki. D.. Gardi, S.. and BancW,0119. Constraints on the duration of the early Toarcan T-OAE and evidence for carbon-reservoir change from the High Atlas (Morocco). Global and Planetary Change. 175: 113-128. Bretagnon. P.. 1974. Termes a longues periodes dans le systeme solatre. Astronomy A Astrophysics. 30: 141-154. Cited By (since 1996): 42. Charbonnier. G.. Boulila. S.. Spangenberg. J.E.. Adatle. T.. Follnu. K.B.. and Laskar. J.. 2018. Obliquity pacing of the hydrological cycle during the Oceanic Anoxic Event 2. Earth and Planetary Science Utters. 499: 266-277. WOS:000444359800024. Croll. J.. 1875. Climate and Time in Their Geological Relations: A Theory of Secular Changes of the Earth's Climate. D. Appleton. Google-Books-ID: 4J8XAAAAYAAJ. Dickey. J.. Bender, P.. Faller. J.. Newhall. X.. Ricklefs. R.. Ries. J., et al.. 1994. Lunar laser ranging - a continuing legacy of the Apollo program. Science. 265 (5171). 482-490. Fang. Q.. Wu. R, Hinnov, L.A., Jing, X.. Wang. X., and Jiang, Q„ 2015, Geologic evidence for chaotic behavior of the planets and its constraints on the third-order eustatic sequences at the end of the Late Paleozoic Ice Age. Palaeogeography, Palaeoclimatology. Palaeoecology, 440: 848-859. Fienga. A.. Manche, H., Laskar. J., and Gastineau. M., 2008. INPOP06: a new numerical planetary ephemeris. Astronomy & Astrophysics. 477: 315-327. Fienga. A.. Laskar. J.. Morley. T.. Manche. H.. Kuchynka, P., Le Poncin-Lafitte. C, et al., 2009, INPOP08. a 4-D planetary ephemeris: from asteroid and time-scale computations to ESA Mars Express and Venus Express contributions. Astronomy & Astrophysics. 507: 1675-1686. Fienga. A.. Laskar. J.. Kuchynka, P., Manche, H., Desvignes G Gasuneau, M„ et al., 2011. The INPOPlOa planetary ephemeris and its apphcauons in fundamental physics. Celestial Mechanics and Dynamical Astronomy. Ill: 363-385. Folkner. W Williams. J., and Boggs, D.. 2009. The planetary and lunar ephemens de 421. In IPN Progress Report, pp 1 -34 Fu' rwivlfT1"- J;m Frequency analysis and re"°" * mm P ^ s°lutions- As,ronom> & A^ics. Gradstein. P.M.. Ogg. J.G., and Smith. A.G 2004 A r , • Scale 2004. Cambridge University Press '°g'C ^ Grafctein. P.. Ogg, J., Schmitz. K.. and Ogg G 2012 Th r 5.