A&A 428, 261-285 (2004)
J. Laskar1 - P. Robutel1 - F. Joutel 1 - M. Gastineau1 - A. C. M. Correia1,2 - B. Levrard1
1 - Astronomie et Systèmes Dynamiques, IMCCE - CNRS UMR8028, 77 Av. Denfert-Rochereau, 75014 Paris, France
2 - Departamento de Física da Universidade de Aveiro, Campus Universitário de Santiago, 3810-193 Aveiro, Portugal
Received 23 May 2004 / Accepted 11 August 2004
We present here a new solution for the astronomical computation of the insolation quantities on Earth spanning from -250 Myr to 250 Myr. This solution has been improved with respect to La93 (Laskar et al. 1993) by using a direct integration of the gravitational equations for the orbital motion, and by improving the dissipative contributions, in particular in the evolution of the Earth-Moon System. The orbital solution has been used for the calibration of the Neogene period (Lourens et al. 2004), and is expected to be used for age calibrations of paleoclimatic data over 40 to 50 Myr, eventually over the full Palaeogene period (65 Myr) with caution. Beyond this time span, the chaotic evolution of the orbits prevents a precise determination of the Earth's motion. However, the most regular components of the orbital solution could still be used over a much longer time span, which is why we provide here the solution over 250 Myr. Over this time interval, the most striking feature of the obliquity solution, apart from a secular global increase due to tidal dissipation, is a strong decrease of about 0.38 degree in the next few millions of years, due to the crossing of the s6+g5-g6 resonance (Laskar et al. 1993). For the calibration of the Mesozoic time scale (about 65 to 250 Myr), we propose to use the term of largest amplitude in the eccentricity, related to g2-g5, with a fixed frequency of 3.200''/yr, corresponding to a period of 405 000 yr. The uncertainty of this time scale over 100 Myr should be about , and over the full Mesozoic era.
Key words: chaos - celestial mechanics - ephemerides - Earth
After the publication of a more precise solution of the Earth by Le Verrier (1856), that took into account the secular perturbations of all main planets except for Neptune, Croll (1890) proposed that the variation of the Earth's eccentricity was also an important parameter for the understanding of the past climates of the Earth.
The first computations of the variations of the obliquity (angle between the equator and orbital plane) due to the secular variations of the orbital plane of the Earth are due to Pilgrim (1904), and were later used by Milankovitch (1941) to establish his theory of the Earth's insolation parameters. Since then, the understanding of the climate response to the orbital forcing has evolved, but all the necessary ingredients for the insolation computations were present in Milankovitch's work.
The revival of the Milankovitch theory of paleoclimate can be related to the landmark work of Hays et al. (1976), that established a correlation between astronomical forcing and the O records over the past 500 kyr. The Milankovitch theory has since been confirmed overall with variations in the climate response to the insolation forcing (see Imbrie & Imbrie 1979; Imbrie 1982, for more historical details; and Zachos et al. 2001; Grastein et al. 2004, for a recent review on the astronomical calibration of geological data).
Since the work of Pilgrim (1904) and Milankovitch (1941), the orbital and precession quantities of the Earth have undergone several improvements. Le Verrier's solution (1856) consisted in the linearized equations for the secular evolution of the planetary orbits. Stockwell (1873), and Harzer (1895) added the planet Neptune to Le Verrier's computations. A significant improvement is due to Hill (1897) who discovered that the proximity of a resonance in Jupiter and Saturn's motion induces some important complements at the second order with respect to the masses. The solution of Brouwer & Van Woerkom (1950) is essentially the solution of Le Verrier and Stockwell, complemented by the higher order contributions computed by Hill. This solution was used for insolation computations by Sharav & Boudnikova (1967a,b) with updated values of the parameters, and later on by Vernekar (1972). The computations of Vernekar were actually used by Hays et al. (1976).
The next improvement in the computation of the orbital data is by Bretagnon (1974), who computed the terms of second order and degree 3 in eccentricity and inclination in the secular equations, but omitted the terms of degree 5 in the Jupiter-Saturn system from Hill (1897). This solution was then used by Berger (1978) for the computation of the precession and insolation quantities of the Earth, following Sharav & Boudnikova (1967a,b). All these works assumed implicitly that the motion of the Solar system was regular and that the solutions could be obtained as quasiperiodic series, using perturbation theory.
When Laskar (1984, 1985, 1986) computed in an extensive way the secular equations for the Solar system, including all terms up to order 2 in the masses and 5 in eccentricity and inclination, he realized that the traditional perturbative methods could not be used for the integration of the secular equations, due to strong divergences that became apparent in the system of the inner planets (Laskar 1984). This difficulty was overcome by switching to a numerical integration of the secular equations, which could be done in a very effective way, with a very large stepsize of 500 years. These computations provided a much more accurate solution for the orbital motion of the Solar system (Laskar 1986, 1988), which also included a full solution for the precession and obliquity of the Earth for insolation computations over 10 Myr (million of years). Extending his integration to 200 Myr, Laskar (1989, 1990) demonstrated that the orbital motion of the planets, and especially of the terrestrial planets, is chaotic, with an exponential divergence corresponding to an increase of the error by a factor of 10 every 10 Myr, thus destroying the hope to obtain a precise astronomical solution for paleoclimate studies over more than a few tens of million of years (Laskar 1999).
The first long term direct numerical integration (without averaging) of a realistic model of the Solar system, together with the precession and obliquity equations, was made by Quinn et al. (1991) over 3 Myr. Over its range, this solution presented very small differences with the updated secular solution of Laskar et al. (1993) that was computed over 20 Myr, and has since been extensively used for paleoclimate computations under the acronym La93. The orbital motion of the full Solar system has also been computed over 100 Myr by Sussman & Wisdom (1992), using a symplectic integrator with mixed variables (Wisdom & Holman 1991), confirming the chaotic behaviour found by Laskar (1989, 1990). Following the improvement of computer technology, long term integrations of realistic models of the Solar system become more easy to perform, but are still challenging when the time interval is the age of the Solar system or if the accuracy of the model is comparable to the precision of short term planetary ephemeris. The longest integration was made over several billions of years by Ito & Tanikawa (2002) with a Newtonian model that did not contain general relativity or lunar contributions, while a recent long term integration of the orbital motion of the Solar system including general relativity, and where the Moon is resolved was made by Varadi et al. (2003) over about 50 Myr.
Here we present the result of the efforts that we have conducted in the past years in our group in order to obtain a new solution for Earth paleoclimate studies over the Neogene period (23 Myr) and beyond. After a detailed analysis of the main limiting factors for long term integrations (Laskar 1999), and the design of new symplectic integrators (Laskar & Robutel 2001), we have obtained a new numerical solution for both orbital and rotational motion of the Earth that can be used over about 50 Myr for paleoclimate studies, and even over longer periods of time if only the most stable features of the solution are used. Apart from the use of a very complete dynamical model for the orbital motion of the planets, and of the new symplectic integrator of Laskar & Robutel (2001) (see next section), a major change with respect to La93 (Laskar et al. 1993) is the consideration, in the precession solution, of a more complete model for the Earth-Moon tidal interactions (Sect. 4.1).
Table 1: Main constants used in La2004. IAU76 refers to the resolutions of the International Astronmical Union of 1976, IERS1992, and IERS2000, refers to the IERS conventions (McCarthy 1992; McCarthy & Petit 2004).The present solution (La2004) or some of its previous variants have already been distributed (Laskar 2001) and used for calibration of sedimentary records over extended time interval (Pälike 2002), and in particular for the construction of a new astronomically calibrated geological time scale for the Neogene period (Lourens et al. 2004). The full solution is available on the Web site www.imcce.fr/Equipes/ASD/insola/earth/earth.html, together with a set of routines for the computation of the insolation quantities following (Laskar et al. 1993).
In addition to the numerical output of the orbital and rotational parameter of the Earth that is available on the Web site, we have made a special effort to provide in Sect. 8 very compact analytical approximations for the orbital and rotational quantities of the Earth, that can be used in many cases for a better analytical understanding of the insolation variations. Finally, in the last sections, the stability of the solution La2004 is discussed, as well as the possible chaotic transitions of the arguments present in the main secular resonances.
The orbital solutions La90-93 (Laskar 1990; Laskar et al. 1993) were obtained by a numerical integration of the averaged equations of the Solar system, including the main general relativity and Lunar perturbations. The averaging process was performed using dedicated computer algebra routines. The resulting equations were huge, with about 150 000 polynomial terms, but as the short period terms were no longer present, these equations could be integrated with a step size of 200 to 500 years, allowing very extensive long term orbital computations for the Solar system.
Although these averaged equations solutions could be improved by some new adjustment of the initial conditions and parameters (Laskar et al. 2004), it appears that because of the improvements of computer technology, it becomes now possible to obtain more precise results over a few tens of million of years using a more direct numerical integration of the gravitational equations, and this is how the present new solution La2004 is obtained. A detailed account of the dynamical equations and numerical integrator is not in the scope of the present paper, but will be made in a forthcoming publication. Nevertheless, we will report here on the main issues concerning the numerical integrator, and will more extensively discuss on the dissipative terms in the integration (see Sect. 4.1).
The orbital model differs from La93, as it comprises now all 9 planets of the Solar System, including Pluto. The post-Newtonian general relativity corrections of order 1/c2 due to the Sun are included following Saha & Tremaine (1994).
The Moon is treated as a separate object. In order to obtain a realistic evolution of the Earth-Moon system, we also take into account the most important coefficient ( ) in the gravitational potential of the Earth and of the Moon (Table 1), and the tidal dissipation in the Earth-Moon System (see Sect. 4.1). We also integrate at the same time the precession and obliquity equations for the Earth and the evolution of its rotation period in a comprehensive and coherent way, following the lines of Néron de Surgy & Laskar (1997), Correia et al. (2003) (see Sect. 3).
In order to minimize the accumulation of roundoff errors,
the numerical integration was performed with the new
symplectic integrator scheme
of (Laskar & Robutel 2001),
with a correction step for the integration of the Moon.
This integrator is particularly adapted to perturbed systems where the
Hamiltonian governing the equations of motion can be
written in the form
as the sum of an integrable part A
(the Keplerian equations of the planets orbiting the
Sun and of the Moon around the Earth), and a small perturbation potential
(here the small parameter
is of the order of the
Using this integrator with step size is then equivalent to integrate exactly a close by Hamiltonian
the error of method
is of the order of
when the correction step
is added, while the same quantity is of the order
in the widely used symplectic integrator of (Wisdom & Holman 1991).
More precisely, the integration step
In Fig. 1a is plotted the evolution of the total energy of the system from -250 Myr to +250 Myr. The energy presents a secular trend that corresponds to the dissipation in the Earth-Moon system. Indeed, after removing the computed energy change due to the secular evolution of the Earth-Moon semi-major axis, the secular trend in the energy evolution disappears and we are left with a residual that is smaller than after 250 Myr, and seems to behave as a random walk (Fig. 1b). The normal component of the angular momentum is conserved over the same time with a relative error of less than (Fig. 1c).
|Figure 1: Conservation of integrals. a) Relative variation of the total energy of the system versus time (in Myr) from -250 Myr to +250 Myr. b) The same after correction of the secular trend due to the tidal dissipation in the Earth-Moon system. c) Relative variation of the normal component of the total angular momentum.|
|Open with DEXTER|
|Figure 2: Normalized repartition of the energy numerical error for 1000 yr steps, over the whole integration, from -250 Myr to +250 Myr. Practically superposed to this curve is the computed normal distribution given by Eqs. (4) and (5).|
|Open with DEXTER|
In order to test the randomness of the numerical error in the integration,
we have plotted the distribution of the difference of total energy
from one 1000 year step to the next. After normalization to 1 of the area
of the distribution curve (Fig. 2), we obtain a curve that is almost
identical to the normalized Gaussian function
|Figure 3: Fundamental planes for the definition of precession and obliquity. Eqt and Ect are the mean equator and ecliptic at date t. Ec0 is the fixed ecliptic at Julian date J2000, with equinox . The general precession in longitude is defined by . is the longitude of the node, and i the inclination. The angle between Eqt and Ect is the obliquity.|
|Open with DEXTER|
Following (Darwin 1880; Mignard 1979), we assume that the torque resulting from tidal friction is proportional to the time lag needed for the deformation to reach the equilibrium. This time lag is supposed to be constant, and the angle between the direction of the tide-raising body and the direction of the high tide (which is carried out of the former by the rotation of the Earth) is proportional to the speed of rotation. Such a model is called "viscous'', and corresponds to the case for which 1/Q is proportional to the tidal frequency.
After averaging over the mean anomaly of the Moon and Earth (indices M and ),
and over the longitude of node and perigee of the Moon,
we have at second order in eccentricity for the contributions of the solar tides
At second order in eccentricity, we
have for the orbital evolution of the Moon (Mignard 1979, 1980, 1981; Néron de Surgy & Laskar 1997)
We have also taken into account the tidal effect raised by the Earth on the Moon,
following Néron de Surgy & Laskar (1997), with the following
simplifying assumptions that the Moon is locked in synchronous spin-orbit resonance
with the Earth (i.e.
), and expanding only at first order in the Moon
). We thus obtain the additional contributions
The core and the mantle have different dynamical ellipticities, so they tend to have different precession rates. This trend produces a viscous friction at the core-mantle boundary (Poincaré 1910; Rochester 1976; Lumb & Aldridge 1991). We will consider here that an effective viscosity can account for a weak laminar friction, as well as a strong turbulent one which thickens the boundary layer (e.g. Correia et al. 2003).
The core-mantle friction tends to slow down the rotation and to bring the obliquity down to if and up to otherwise, which contrasts with the effect of the tides. Furthermore, one can see that, despite the strong coupling by pressure forces, there can be a substantial contribution to the variation of the spin for high viscosities and moderate speeds of rotation.
The Earth's atmosphere also undergoes some torques which can be transmitted to the surface by friction: a torque caused by the gravitational tides raised by the Moon and the Sun, a magnetic one generated by interactions between the magnetosphere and the solar wind. Both effects are negligible; see respectively Chapman & Lindzen (1970) and Volland (1978).
Finally, a torque is produced by the daily solar heating which induces a redistribution of the air pressure, mainly driven by a semidiurnal wave, hence the so-called thermal atmospheric tides (Chapman & Lindzen 1970). The axis of symmetry of the resulting bulge of mass is permanently shifted out of the direction of the Sun by the Earth's rotation. As for the body tides, this loss of symmetry is responsible for the torque which, in the present conditions, tends to accelerate the spin (Correia & Laskar 2003).
Volland (1978) showed that this effect is not negligible, reducing the Earth's despinning by about 7.5. But, although the estimate of their long term contributions deserves careful attention, we have not taken the atmospheric tides into account in the computations presented in the next sections, assuming that, as the uncertainty on some of the other factors is still large, the global results obtained here will not differ much when taking this additional effect into consideration.
Several internal geophysical processes are able to affect the moments of
inertia of the Earth and thereby its precession constant value.
Redistribution of mass within the Earth can
occur as a consequence of plate subduction, upwelling plumes, mantle
avalanches or density anomalies driven by mantle convection. Typical
time-scales of these phenomena range from 10 to 100 Myr. However, their impact is still poorly known.
Forte & Mitrovica (1997) computed convection-induced perturbations in
the dynamical ellipticity of the Earth over the last 20 Myr, using
different sets of 3D seismic heterogeneity models and mantle viscosity
profiles. Most of the models indicate that mantle convection led to a
mean relative increase of the precession constant close to since the beginning of the Neogene period (23 Myr ago).
Although this value is similar to our estimate of the tidal dissipation
effect (but in the opposite direction, see Sect. 9), we did not take
account of its impact in our computations.
Because of its still large uncertainty, we feel that it is more
appropriate to include its possible contribution as a perturbation of the
tidal dissipation parameter.
|Figure 4: Differences La2003-DE406 over the full range of DE406 (-5000 yr to +1000 from J2000). The units for semi-major axis (a) are AU, arcsec for mean longitude (l) and longitude of perihelion (perihelion). The eccentricity (e) and the inclinations variables ( ), where i and are the inclination and node from the ecliptic and equinox J2000.|
|Open with DEXTER|
Climate friction is a positive and dissipative feedback between obliquity variations and climate which may cause a secular drift of the spin axis. In response to quasi-periodic variations in the obliquity, glacial and interglacial conditions drive transport of water into and out of the polar regions, affecting the dynamical ellipticity of the Earth. A significant fraction of the surface loading is compensated by viscous flow within the Earth, but delayed responses in both climatic and viscous relaxation processes may introduce a secular term in the obliquity evolution (Rubincam 1990, 1995). Because both the ice load history and the visco-elastic structure of the Earth are not strongly constrained, it is difficult to produce accurate simulations and predictions of the coupled response of the entire system. Nevertheless, simplifying assumptions have been used to estimate the magnitude and the direction of the secular drift. Several analyses have examined this phenomenon, suggesting that climate friction could have changed the Earth's obliquity by more that 20 over its whole geological history (Bills 1994; Ito et al. 1995; Williams et al. 1998). A more detailed theoretical and numerical treatment of climate friction has been developped by Levrard & Laskar (2003). Using available constraints on the ice volume response to obliquity forcing based on O oxygen-isotope records, they showed that climate friction impact is likely negligible. Over the last 3 Ma, corresponding to the intensification of the Northern Hemisphere glaciation, a maximal absolute drift of only /Myr has been estimated for a realistic lower mantle viscosity of 1022 Pa. s. Before 3 Ma, climate friction impact is expected to be much more negligible because of the lack of massive ice caps. In this context, this effect was not taken into account in our long-term computations. Indeed, the uncertainty on the tidal dissipation is still important and it is expected that most of the possible effect of climate friction would be absorbed by a small change of the main tidal dissipation parameters.
Nevertheless, glaciation-induced perturbations in the dynamical ellipticity of the Earth change the precession frequency and thereby the obliquity frequencies leading to a temporal shift between a perturbated and a nominal solution (see Sect. 9). Estimates of this offset have been performed by Mitrovica & Forte (1995) and Levrard & Laskar (2003) for a large set of mantle viscosity profiles. It may reach between 0.7 and 8 kyr after 3 Myr, corresponding respectively to extreme lower mantle viscosities of and 1023 Pa. s. In the last case, the Earth is nearly rigid on orbital timescales and the corresponding offset can be thus reasonably considered as a maximal value.
Table 2: Maximum difference between La2004 and DE406 over the whole time interval of DE406 (-5000 yr to +1000 yr with origin at J2000); Col. 1: -100 to + 100 yr; Col. 2:-1000 to + 1000 yr; Col. 3: -5000 to +1000 yr. EMB is the Earth-Moon barycenter.Using a direct numerical integrator, our goal is to provide a long term solution for the orbital and precessional elements of the Earth with a precision that is comparable with the usual accuracy of a short time ephemeris. We have thus compared our solution with the most advanced present numerical integration, DE406, that was itself adjusted to the observations (Standish 1998). Over the full range of DE406, that is from -5000 yr to + 1000 yr from the present date, the maximum difference in the position of the Earth-Moon barycenter is less than 0.09 arcsec in longitude, and the difference in eccentricity less than 10-8. The difference in the longitude of the Moon are more important, as the dissipative models that are used in the two integrations are slightly different. They amount to 240 arcsec after 5000 years, and the eccentricity difference is of the order of , but it should be noted that only the averaged motion of the Moon will have some significant influence on the precession and obliquity of the Earth. Table 2 summarizes the maximum differences between La2004 and DE406 for the main orbital elements, over the whole interval of DE406, for different time intervals.
|Figure 5: Eccentricity of the Earth over 25 Myr in negative time from J2000. The solid line stands for the present solution La2004, while the dotted line is the eccentricity in the La93 solution (Laskar et al. 1993). The differences of the two solutions becomes noticeable after 10 Myr, and significantly different after 15-20 Myr.|
|Open with DEXTER|
|Figure 6: Inclination (in degrees) of the Earth with respect to the fixed ecliptic J2000 over 25 Myr in negative time from J2000. The solid line stands for the present solution La2004, while the dotted line is the eccentricity in the La93 solution (Laskar et al. 1993). The differences of the two solutions become noticeable after 15 Myr, and significantly different after 20 Myr.|
|Open with DEXTER|
Table 3: Main secular frequencies gi and si of La2004 determined over 20 Ma for the four inner planets, and over 50 Ma for the 5 outer planets (in arcsec yr-1). and are the observed variations of the frequencies over respectively 100 and 250 Myr. In the last column, the period of the secular term are given.We have compared in Figs. 5 and 6 the eccentricity and inclination of the secular solution La93 and the new present solution La2004. Over 10 Myr, the two orbital solutions La93 and La2004 are nearly identical, and do not differ significantly over 20 Myr. The differences in eccentricity increase regularly with time, amounting to about 0.02 in the eccentricity after 20 Myr, about 1/3 of the total amplitude 0.063, while the difference in the orbital inclination reaches 1 degree, to compare to a maximum variation of 4.3 degrees. These differences result mostly from a small difference in the main secular frequency g6 from Jupiter and Saturn that is now g6=28.2450 arcsec/year (Table 3) instead of 28.2207 arcsec/year in the previous solution. Beyond 20 million years, the differences between the two solutions become more noticeable. The solutions in obliquity and climatic precession are plotted in Figs. 7 and 8. The difference in obliquity over 20 Myr amounts to about 2 degrees, which means that the obliquity cycles of the two solutions become nearly out of phase after this date. The differences La2004-La93 are plotted in (a), while in (b) only the precession model has been changed and in (c) only the orbital model has been changed. Comparing 7b and 7c shows that most of the differences between the La93 and La2004 obliquity solutions result from the change in the dissipative model of the Earth-Moon system, while the only change of the orbital solution would lead to a much smaller change of only 0.6 degree in obliquity after 20 Myr. Similar conclusions can be made for the climatic precession (Figs. 8a-c).
|Figure 7: Differences in obliquity (in degrees) La2004-La93(1, 1) versus time (in Myr) over 20 Myr from present a). In La2004 both orbital and precession models have been changed. In b) only the precession model has been changed, while in c) only the orbital motion has been changed. It is thus clear that the main change in La2004 arise from the change of precession model.|
|Open with DEXTER|
|Figure 8: Differences in climatic precession ( , where is the perihelion from the moving equinox) La2004-La93(1, 1) versus time (in Myr) over 20 Myr from present a). In La2004 both orbital and precession models have been changed. In b) only the precession model has been changed, while in c) only the orbital motion has been changed. It is thus clear that the main change in La2004 arise from the change of precession model.|
|Open with DEXTER|
Over long time scales, the determination of longitude of the planet on its orbit is not very important, and we are more concerned by the slow evolution of the Earth orbit under secular planetary perturbations. The semi-major axis of the Earth has only very small variations (Fig. 11). Indeed, as shown by Laplace (1773) and Lagrange (1776), there are no secular variations of the semi-major axis of the planets at first order with respect to the masses, while some terms of higher order can be present (Haretu 1885). These terms are very small for the inner planets, but more visible in the solutions of the outer planets where the proximity of some mean motion resonances increases their contribution at second order with respect to the masses (Milani et al. 1987; Bretagnon & Simon 1990). For the solution of the Earth, they are so small that they will not induce any noticeable change of the mean Earth-Sun distance, or of the duration of the year in the geological past, at least over the past 250 millions of years.
It was first shown by Lagrange (1774, 1777) that the inclination and nodes of the planets
suffer long term quasiperiodic variations. Shortly after, Laplace (1776) demonstrated that
this was the same for the eccentricity and longitudes.
Both computations were made using the linearization of the first order averaged equations of
motion. In this approximation, if we use the complex notations
Over a finite time of a few millions of years, one can still numerically integrate the equation of motion, and then search through frequency analysis (Laskar 1990, 1999, 2003) for a quasiperiodic approximation of the solution of the form (23). The non regular behavior of the solution will induce a drift of the frequencies in time, related to the chaotic diffusion of the trajectories (Laskar 1990, 1993, 1999).
In Table 3, the fundamental frequencies of the solution La2004 have been computed by frequency analysis over 20 Myr of the proper modes related to the inner planets (Mercury, Venus, Earth, Mars), and over 50 Myr for the proper modes associated to the outer planets (Jupiter, Saturn, Uranus, Neptune, Pluto) (Table 3). This difference reflects the fact that the main source of chaotic behavior comes from overlap of secular resonances in the inner planet system, and that the outer planets secular systems is more regular than the inner one. Actually, in Cols. 2 and 3 of Table 3, we give the maximum observed in the variation of the secular frequencies over 100 and 250 Myr, in both positive and negative time, for our nominal solution La2004. These figures are in very good agreement with the equivalent quantities given in Table X of (Laskar 1990) for the integration of the secular equations. The number of digits of the fundamental frequencies in Table 3 is different for each frequency and reflects the stability of the considered frequency with time. Additionally, we have given in Figs. 9 and 10 the actual evolution of the secular frequencies obtained with a sliding window with a step size of 1 Myr (compare to Laskar 1990, Figs. 8 and 9).
It is interesting for practical use to have an analytical expression for
the main orbital quantities of the Earth.
From the numerical values of La2004, we have performed a frequency analysis
(Laskar 1990, 1999, 2003) in order to obtain a quasiperiodic approximation of
the solutions over a few Myr.
In order to be consistent with the remaining part of the paper,
we chose a time interval covering 20 Myr.
As we are mostly interested over negative time, we made the analysis
from -15 Myr to +5 Myr, as usually, the precision of the approximation decreases
at the edges of the time interval.
The analysis of the eccentricity variables
is given in Table 4.
z is obtained as
|Figure 9: Variation of the secular frequencies g1-9 from -250 to +250 Myr. The frequencies are computed over 20 Myr for g1-4 and over 50 Myr for g5-9, after transformation of elliptical elements to proper modes (Laskar 1990).|
|Open with DEXTER|
|Figure 10: Variation of the secular frequencies s1-9 from -250 to +250 Myr. The frequencies are computed over 20 Myr for s1-4 and over 50 Myr for s6-9, after transformation of elliptical elements to proper modes (Laskar 1990).|
|Open with DEXTER|
Table 4: Frequency decomposition of for the Earth on the time interval [-15,+5] Myr (Eq. (25)).
Nevertheless, for information, we will give also here the leading terms in the expansion of the eccentricity as they may be useful in paleoclimate studies (Table 6). The three leading terms in eccentricity are well known in paleoclimate studies g2-g5 (405 kyr period), g4-g5 (95 kyr period), and g4-g2 (124 kyr period).
The solution for the inclination variables
are the Earth
inclination and longitude of node with respect to the fixed ecliptic and equinox J2000, limited to 24 quasiperiodic terms
The solutions for precession and obliquity are obtained through the precession Eq. (6).
In absence of planetary perturbations, the obliquity is constant (
Another difficulty, arising with the obliquity solution, is due to the dissipation in the Earth-Moon system,
which induces a significant variation of the precession frequency with time, as the Earth rotation
slows down, and the Earth-Moon distance
|Figure 11: Variation of the semi-major axis of the Earth-Moon barycenter (in AU) from -250 to +250 Myr.|
|Open with DEXTER|
Table 5: Frequency decomposition of for the Earth on the time interval [-15,+5] Myr (Eq. (26)).
|Figure 12: ( Top) Eccentricity of the Earth from -11 to +1 Myr. The solid curve is the solution La2004. Almost completely hidden behind the solid curve is a dashed curve representing the quasiperiodic approximation of Table 4 with 26 periodic terms. ( Bottom) Inclination of the Earth from -11 to +1 Myr. The solid curve is the solution La2004. The dashed line (almost identical to the solid curve) is the quasiperiodic approximation of Table 5 with 24 periodic terms. On both plots, as the two curves are nearly identical, the difference of the two solutions is also plotted in dotted curve.|
|Open with DEXTER|
We have plotted in Fig. 14 the evolution of the obliquity of the Earth from -250 to + 250 Myr.
The effect of the tidal dissipation is clearly visible on this timescale, and there is a general increase
of the obliquity from -250 Myr to present time. In positive time, there is an obvious singularity
in the obliquity, with a decrease of about 0.4 degree. This results from the crossing of the resonance
with the term
s6+ g5-g6 that was already described in a different context in Laskar et al. (1993).
After this sudden decrease, the obliquity increases again, with roughly the same trend.
An additional (but smaller singularity is observable around +150 Myr, that corresponds
to the crossing of the resonance
with the term
(Table 5). The averaged obliquity can be well approximated in negative time
with the expression
Table 6: Frequency decomposition of the eccentricity of the Earth on the time interval [-15,+5] Myr. For all computations, the data of Table 4 should be preferred. Here, with e0=0.0275579. For each term, in the first column is the corresponding combination of frequencies where gi are the fundamental frequencies (Table 3), and the frequencies of the terms in z from Table 4.
The comparison with the full solution is given in Fig. 15, over selected periods of time, from -11 Myr to the present. In Table 7, the terms of frequency are recognized as combinations of the form , where p is the precession frequency, and a secular frequency from the secular motion of the orbital plane (from Table 5). The approximation is made over the interval with p = p0 +p1 T, where p0 and p1 are taken from Eq. (35).
Table 7: Approximation for the obliquity of the Earth, following Eq. (33). This expression is not strictly quasiperiodic, because of the presence of the dissipative term p1 in the evolution of the precesion frequency (33).
|Figure 13: a) Evolution of length of the day for the Earth, in hours, from -250 to +250 Myr. b) Residuals with the fit of the averaged solution with the polynomial expression (41).|
|Open with DEXTER|
The full approximate solution for the precession is more complicated because of the presence of
the secular resonance
s6+g5-g6. In fact, the proximity of this resonance is responsible for the largest periodic term in the expression of the precessionangle (Fig. 16).
The precession angle can thus be approximate by (
in arcsec, t in years)
|Figure 14: Evolution of the obliquity of the Earth in degrees, from -250 to +250 Myr. The grey zone is the actual obliquity, while the black curve is the averaged value of the obliquity over 0.5 Myr time intervals. The dotted line is a straight line fitted to the average obliquity in the past.|
|Open with DEXTER|
|Figure 15: Comparison of the solution of the obliquity La2004 (solid line) with its approximation using Table 7 (dotted line) with 26 periodic terms. The difference of two solutions (+22 degrees) is also plotted.|
|Open with DEXTER|
In Fig. 18, the evolution of the Earth-Moon semi-major axis is plotted over 250 Myr in positive and negative time. The grey area corresponds to the short period variations of the Earth-Moon distance, obtained in the full integration of the system, while the black curve is the integration of the averaged equations that is used for the precession computations. The agreement between the two solutions is very satisfying.
We have searched for a polynomial approximation for the secular evolution of ,
in order to provide
a useful formula for simple analytical computations. In fact, the Earth-Moon
evolution is driven by the tidal dissipation equations (Sect. 4.1),
that involve the obliquity of the planet.
The important singularity of the obliquity in the future (see previous section)
prevents from using a single approximation
over the whole time interval, and we have used a different polynomial in positive and negative time.
In the same way, it is possible to approximate the evolution of the precession frequency of the Earth p
as the Moon goes away and the rotation of the Earth slows down with time. The polynomial approximation
of the precession frequency is then in arcsec yr-1, with T in Gyr
|Figure 16: Comparison of the solution of the precession La2004 with its approximation from -15 to +2 Myr. The grey curve is obtained after removing uniquely the secular trend . The dark curve are the residual (in radians) after removing the resonant term (Eq. (34)), which is also displayed in solid line.|
|Open with DEXTER|
|Figure 17: Comparison of the solution of the climatic precession of La2004 with its approximation from -11 to +1 Myr. The grey curve is the full climatic precession . The dark curve are the residual after removing the approximation given by Eq. (37).|
|Open with DEXTER|
|Figure 18: a) Evolution of the Earth-Moon semi major axis (in Earth radii) from -250 to +250 Myr. The grey zone is the result of the integration of the full equations, while the black curve is the integration of the averaged equations, as used in the precession computations; b) residuals with the fit of the averaged solution with the polynomial expression (39).|
|Open with DEXTER|
|Figure 19: a) Evolution of the precession frequency of the Earth p (in arcsec/yr) -250 to +250 Myr; b) residuals with the fit of the averaged solution with the polynomial expression (40).|
|Open with DEXTER|
In such long term numerical computations, the sources of errors are the error of method due to the limitation of the numerical integration algorithm, the numerical roundoff error resulting from the limitation of the representation of the real numbers in the computer, and the errors of the initial conditions and of the model which are due to our imperfect knowledge of all the physical parameters in the Solar system, and to the necessary limitations of the model. Some of the main sources of uncertainty in the model for long term integrations were reviewed in (Laskar 1999). The effect of all these errors is amplified by the chaotic behaviour of the system, with an exponential increase of the difference between two solutions with different settings, until saturation due to the limited range of the considered variables. A summary of these limitations, extracted from (Laskar 1999) is given in Table 8, but it should be noted that these times of validity are probably pessimistic, as they correspond to analytical estimates based on the "worst case'' situation. Our requirement for the precision of a long term solution is also not the same as for precise short term ephemeris. For the long term solutions, aimed at paleoclimate or qualitative studies, two solutions can be considered as similar, as long as the pattern of the two solutions (eccentricity for example), resulting (at first order) from the combination of various proper modes remain similar. This will last until some of the main proper modes get completely out of phase (see Sect. 7). It will be the same for the solutions in obliquity and precession.
Table 8: Main sources of uncertainty in the orbital solution (from Laskar 1999). For each limiting factor, an analytical estimate of the time of validity of the solution is given (in Myr), taking into account the exponential growth of the error.
|Figure 20: Stability of the solution for eccentricity of the Earth. a) Difference of the nominal solution La2004 with stepsize years, and La2004*, obtained with years; b) difference of the nominal solution with the solution obtained while setting for the Sun (instead of in the nominal solution).|
|Open with DEXTER|
Practically, we have tested the stability of our solution by comparison of the eccentricity (Fig. 20) and inclination (Fig. 21) with different solutions. In Figs. 20a and 21a, the nominal solution La2004 (with stepsize years) is compared to an alternate solution, La2004*, with the same dynamical model, and a very close stepsize years (Table 9). This special value was chosen in order that our output time span h=1000 years corresponds to an integer number (204 800) of steps, in order to avoid any interpolation problems in the check of the numerical accuracy. This is thus a test of the time of validity for obtaining a precise numerical solution, resulting from method and roundoff errors of the integrator. It is thus limited here to about 60 Myr.
This limitation of 60 Myr is a limitation for the time of validity of an orbital solution, independently of the precision of the dynamical model. In order to go beyond this limit, the only way will be to increase the numerical accuracy of our computations, by improving the numerical algorithm, or with an extended precision for the number representation in the computer.
A second test is made on the uncertainty of the model. Several sources of
uncertainty were listed in (Laskar 1999), and it was found that
one of the main sources of uncertainty was due to the imprecise knowledge
of the Solar oblateness (
Its most recent determination, obtained with
the SOHO and GONG helioseismic data give
(Pijpers 1998), with a very similar value adopted in DE406 (
Standish 1998), while in DE200, it was not taken into account (Newhall et al. 1983).
We will thus consider that comparing the nominal solution La2004 (with
with an alternate solution (La20040) with
is representative of the uncertainty
of the dynamical model for our long term integrations
|Figure 21: Stability of the solution for inclination of the Earth. a) Difference of for the nominal solution La2004 with stepsize years, and for La2004*, obtained with years; b) difference of for the nominal solution with the solution obtained while setting for the Sun (instead of in the nominal solution).|
|Open with DEXTER|
Table 9: Variants of the La2004 solutions. The nominal solution La2004 is obtained for a stepsize , and a Solar oblateness . The other solutions, with different setting have been computed to test the stability of the nominal solution. In the last column is the tidal dissipation factor.
The results La2004 - La20040 for the Earth's eccentricity and inclination are displayed in (Figs. 20b and 21b) over 100 Myr. As for Mars (Laskar et al. 2004), the effect of becomes noticeable after about 30 Myr (26 Myr were predicted with an analytical estimate in (Laskar 1999), and the solution remains very similar over 40 Myr, and totally out of phase after 45 Myr. We will thus consider here that 40 Myr is about the time of validity of our present orbital solution for the Earth. Other integrations made with and (Table 9) gave practically the same estimates. It should be noted that with our present algorithms, we are much more limited by the precision of the model (Figs. 20b and 21b) than by the numerical accuracy (Figs. 20a and 21a).
For the obliquity and precession, as was already mentioned in (Laskar et al. 1993), the main source of uncertainty is the evolution of the tidal dissipation and other related factors (see Sect. 4.1). Among these factors, one can mention the variations of the dynamical ellipticity of the Earth resulting from changes of the ice caps, or mantle convections. The maximal contribution of these effects can be estimated (see Levrard & Laskar 2003), but it is very difficult to obtain at present a precise value of their contributions. Moreover, the largest dissipative term, resulting from body tides can evolve substantially as a result of the variation of the distribution of the continents, as most of the tidal dissipation occurs in shallow seas (Jeffreys 1920, 1976; Egbert & Ray 2000).
In Fig. 22, we have explicited the result of
a difference of 5% (Figs. 22a and c) and 10% (Figs. 22b and d)
in the tidal dissipation term of the Earth and the Moon (
in Table 1).
With 5% error (La2004(0.95)), the solution of
obliquity is valid over about 20 Myr, but with a 10% error (La2004(0.90)),
the solution is out of phase after 20 Myr.
Nevertheless, the uncertainty in the tidal dissipation appears mostly
as a small change of the precession frequency p, which
reveals in the obliquity solution as
a time offset that does not change much the obliquity pattern.
Indeed, in Figs. 22e and f, we have made a readjustment of the
time scale for the solutions La2004(0.95) and La2004(0.90)of respectively
|Figure 22: Difference in obliquity for different tidal dissipation factor (Table 9): a) ; b) . In c) and d), the nominal solution of the obliquity is plotted in solid line, while La2004(0.95) c) and La2004(0.90) d) are in dotted line. After an adjustment of the time scale, the differences of the solutions are much smaller over 40 Myr: e) ; f) .|
|Open with DEXTER|
|Figure 23: Resonant arguments (in radians versus time in Myr) ( top) and ( bottom) for different solutions: a) is the secular solution La93 (Laskar 1990; Laskar et al. 1993) while b) is the present nominal solution La2004. The characteristics of the other variants La20041.5 c), La20040 d), La2004** e) are listed in Table 9. The arrows indicate the portions of the solutions when the argument (s4-s3)-(g4-g3) is in libration.|
|Open with DEXTER|
In all variants of La2004 (Fig. 23 (top)), the resonant argument corresponding to was in libration over the first 40 Myr, while the transition to circulation occurred at about 25 Myr in La93. All the solutions La2004, La2004**, La20040, La20041.5 differ significantly beyond 40 to 50 Myr, but they all enter the (s4-s3)-(g4-g3) resonance within 100 Myr, comforting the conclusions of (Laskar 1990, 1992). In Fig. 23 (top), the portions of the solutions where the argument of (s4-s3)-(g4-g3) is in libration are pointed by arrows. It should be noted that at least over 250 Myr, the diffusion seems to remain confined between the resonances (s4-s3)-2(g4-g3) and (s4-s3)-(g4-g3). The other important resonant argument ( (g1-g5) -(s1-s2)) that was identified by Laskar (1990, 1992) as the origin of the chaotic behavior of the inner planets presents in all the variations of the numerical solutions La2004, La2004**, La20040, La20041.5 a similar behavior as for the secular solution La93 (Fig. 23 (bottom)).
The date of the first transition from libration to circulation of the resonant argument related to is a macroscopic feature displaying the chaotic diffusion of the planetary orbits. It is thus very tempting to search whether this feature could be observed in the geological data that extend now in continuous records beyond 30 Myr.
The direct observation of the individual arguments related to g3,g4,s3,s4 is certainly out of reach. But it may be possible to follow the evolution in time of the differences g3-g4 and s3-s4. Indeed, s3-s4 appears as a beat of about 1.2 million of years in the solution of obliquity, as the result of the beat between the p+s3 and p+s4 components of the obliquity, where p is the precession frequency of the axis (see Laskar 2004, Fig. 7). In a similar way, g3-g4 appears as a beat with period of about 2.4 Myr in the climatic precession. One understands that because of the occurrence of these beats, the detection of the resonant state in the geological data can be possible, as one has to search now for phenomena of large amplitude in the geological signal. Indeed, in the newly collected data from Ocean Drilling Program Site 926, the modulation of 1.2 Myr of the obliquity appears clearly in the spectral analysis of the paleoclimate record (Zachos et al. 2001). Moreover, using the ODP legs 154 and 199, Pälike et al. (2004) could find some evidence that the critical argument of did not show a transition to circulation at 25 Myr, as in La93, but remained in libration over 40 Myr, as in the present new solution La2004. This result is very impressive, although the search for the determination of these resonant angles in the geological data is just starting. We may expect that a careful analysis, and new data spanning the most recent 60 to 100 Myr, will allow to determine the possible succession of resonant states in the past, allowing, for example, to discriminate between the solutions displayed in Fig. 23. Of particular importance would be the detection of the first transition from libration to circulation of the resonant argument (s4-s3)-2(g4-g3). This program, if completed, will provide some extreme constraint for the gravitational model of the Solar System. Indeed, the observation of a characteristic feature of the solution at 40 to 100 Myr in the past, because of the exponential divergence of the solutions, will provide a constraint of 10-4 to 10-10 on the dynamical model of the Solar System.
|Figure 24: Stability of the g5 ( top) and g2-g5 ( bottom) arguments. The difference (in radians) of the angles related to g5 and g2-g5 from the nominal solution La2004 and an alternate solution La20041.5 a), La20040 b), and La2004* c) (Table 9).|
|Open with DEXTER|
An astronomically calibrated time scale for geological data has been achieved successfully over the Neogene period (0-23.03 Myr) (Lourens et al. 2004), and the next challenge will be to improve the astronomical precision of the solution, in order to extend this over the Paleogene (23 to 65 Myr). Beyond this time length, it is hopeless to expect to obtain a precise solution for the orbital elements of the Earth. On the other hand, some features of the solution have a longer time of stability, if they are related to the outer planets secular frequencies (Laskar 1990).
A very important period present in geological records is the 405 kyr period related to the g2-g5 argument, that is the leading term in the solution of eccentricity (Table 6). In Fig. 24, we have tested the stability of this argument over the full period of our integrations, that is over 250 Myr, that extents over the Mesozoic geological era (about 65 to 250 Myr). First, we have compared the argument related to g5 over 3 variants (La2004*, La20040, La20041.5) of La2004, that reflect both the numerical errors and the errors of the model. This argument, related to Jupiter and Saturn is extremely stable, and does not show variations of more than 0.05 radians over the full 250 Myr period (Fig. 24 (top)). In Fig. 24 (bottom), we have plotted the same quantities for g2-g5. The dispersions over 50, 100, and 250 Myr are about 0.3, 1.0, 6.2 radians, that correspond to absolute error in time respectively of about 20, 70, and 400 kyr over the same periods.
It is thus possible to use this argument, that is visible in many sedimentary records
of the Early Mesozoic (Olsen & Kent 1999, and references therein) as an absolute geological time scale
over the whole Mesozoic period. In order to facilitate the establishment of this timescale,
we have searched for a best fit or the period of this argument
over all the solutions of Table 9.
Within the uncertainty of this procedure, the best value for the
frequency g2-g5 is
|Figure 25: 405 kyr term in eccentricity. Maximum difference (in radians) of the argument of g2-g5 in all 6 solutions of Table 9 (La2004, La2004*, La2004**, La20040, La20041.0, La20041.5) with respect to the linear approximation where t is in yr.|
|Open with DEXTER|
|Figure 26: Eccentricity of the Earth. Nominal solution La2004 filtered in the interval (solid line), and approximation with a single periodic term (47) (dotted line).|
|Open with DEXTER|
When taking into account the dissipation in the Earth-Moon system, it is remarkable that in the time interval , the only major event in the obliquity of the Earth occurs in the vicinity of the present time (Fig. 14). We have discarded the possibility of some numerical artefact by integrating different models. In particular, we find the same slow decrease of about 0.38 degrees in the obliquity over about 20 Myr either when we integrate the obliquity and precession equations of Sect. 4 at the same time as the orbital equations, or when the precession and obliquity are computed afterwards, with a forced orbital motion. Moreover, the behavior of the solution is perfectly explained by the crossing of the s6+g5-g6 resonance, as it is explicited in full details in the appendix.
The crossing of this resonance will occur in the next 20 Myr, but one could wonder whether this crossing could have also occured in the recent past and thus had a possible impact on the past Earth climate.
In (Laskar et al. 1993), we searched for the possibility to encounter this resonance in the past by changing the momentum of inertia of the Earth as a result of the ice load changes occuring during an ice age, but it seems that our values for the change of momentum of inertia during an ice age was overestimated (Mitrovica & Forte 1995; Levrard & Laskar 2003). In a similar way, Forte & Mitrovica (1997) have advocated the mantle convection as a mean to change the values of the momentum of inertia of the Earth in the past, and thus to eventually cross the s6+g5-g6 resonance in the past. Their results are very sensitive to their models, and it is not clear that the induce change of momentum of inertia will be able to counteract the effect of tidal dissipation that drives the obliquity away from the resonance in the past. Moreover, the analysis of geological data in the past few Myr (Lourens et al. 2001; Pälike & Shackleton 2001) tends to show that the precession frequency has actually increased in the past and thus, that the system was driven away from the resonance. Additional studies on this topic are certainly welcome, and we should still search for the possibility of a crossing of the s6+g5-g6 resonance in the past, as it may have impacted the climatic history of the Earth, but in absence of new evidence, we should consider that the proximity of this resonance is pure coincidence.
The improvement of the model and parameters is much more difficult, and will require additional work. First, we will have to include in our model supplementary contributions in order to compare much more closely to DE406. In particular, in the present work, to save computer time, we had to simplify the interactions in the Earth-Moon system. In the next integration, we will take into account additional contributions in the body body interactions of the Earth-Moon system, and also the libration of the Moon. We will also have to adapt the tidal model in the Earth-Moon system. We will need to take into account the effect of the 5 main asteroids that are present in DE406, and probably derive a small model for the contribution of the 295 additional asteroids that are considered in DE406. Indeed, the consideration of the main asteroids appears as the next most important source of uncertainty for the long term solution (Table 8). Some of the other factors listed in Table 8 can easily be taken into account, as the consideration of the main Galilean satellites, or the mass loss of the Sun.
Finally, it will also be important to reestablish the precision of the parameter of the model by a direct comparison to the observations. For this part, we expect to benefit from the effort that is developed at present in our laboratory for the improvement of short time ephemeris devoted to the reduction of Earth based and space missions observations (Fienga & Simon 2004).
The fulfilment of this program will allow to improve on the main limiting factors listed in Table 8. It should be stressed again here that the limiting times listed in Table 8 are somewhat pessimistic. We can thus hope that after completing this effort, the orbital solution may be valid over 50 to 65 Myr.
Nevertheless, although this program may lead to a solution valid over about 60 Myr, obtaining a precise solution for the precession and obliquity of the Earth's axis over 65 Myr is much more difficult, as it will require a modelization of the evolution of the various dissipative contributions (see Sect. 4.1) over this period. An alternate way will be to fit the evolution of the main precession frequency to the geological data over the Paleogene period. Although some results have been already obtained in this direction over some shorter time span (Lourens et al. 2001; Pälike & Shackleton 2001), to extract the variation of the precession frequency from paleoclimate stratigraphic data is still a challenge, but we hope that the present paper will provide a good basis for this future work.
In this section, we show that Fresnel integrals are very good approximations
for the crossing of the resonance.
Let us define the Fresnel integrals (we have dropped here the usual
|Figure A.2: Obliquity variation (in degree) versus time resulting from the crossing of the resonance in the term (solid line). In the same plot, the dependence of the precession phase is illustrated. for the obliquity offset is null (dashed line) while it is maximal for . The nominal solution (solid line) corresponds nearly to a maximal decrease of the obliquity.|
J.L. thanks J. Landwehr for providing a copy of Pilgrim's work, A. Vernekar for a copy of his monograph, and H. Pälike for his comments. This work benefited from support from INSU-CNRS, PNP-CNRS, and CS, Paris Observatory.