3 Results and discussion

3.1 Tests of the model

Our spin axis model includes all instantaneous effects, e.g. precession and nutation. This is demonstrated by a number of tests. Figure 3 displays the nutational wobble of the Earth. The curve shows the familiar 18.6 year period with an 9.2 arcsec amplitude (as measured w.r.t present ecliptic), in agreement with observation. Faster and smaller oscillations, with periods of one half year and one half sidereal month are also seen in the inset figure.

In Fig. 4 we plot the instantaneous spin precession rate (see Sect. 2.9). A short period corresponding to one half of a sidereal month is easily seen, as is a half year period. A longer period corresponding to the nutation period (18.6 years) is also present (not seen in Fig. 4 due to the limited time span). Averaging over an integer number of nutation periods (93 years) yields an annual average precession rate of 50.42 arcsec/yr, again in agreement with observation.

Figure 5 shows the variation of the Earth's precession cycle over the last 2 million years. The cycle is rather stable with a maximum difference of only about 800 years, or 3%. The mean precession cycle length is 25709 yr. This corresponds to a mean annual precession of 50.41 arcsec/yr, in excellent agreement with the aforementioned result of 50.42. These results are a confirmation of the stability of the model.

Figure 1: Instantaneous torques (divided by $\sin2\theta$) on Mars from Phobos in x-, y- and z-direction. The bottom subfigure displays the orbit inclination of Phobos w.r.t. the Mars equator.

Figure 2: Instantaneous torques (divided by $\sin2\theta$) on Mars from Deimos in x-, y- and z-direction. The bottom subfigure displays the orbit inclination of Deimos w.r.t. the Mars equator.

Figure 3: Nutation of the Earth. The main figure shows the 18.6 yr period caused by the precessing nodes of the lunar orbit. Inset displays the fine structure.

Figure 4: Instantaneous spin precession rate of the Earth over one yr.

Figure 5: Length of the precession cycles of Earth's spin axis during the last 2 Myr. The dashed line shows the mean value of 25 709 yr.

Figure 6 gives the inclination of the instantaneous ecliptic relative to the present ecliptic for the time interval -1 Myr to +1 Myr. The result is in good agreement with that of Muller & MacDonald (1997). From the above tests we conclude that the model yields trustable results to a high degree of accuracy.

Figure 6: Inclination of the Earth orbit w.r.t. the present ecliptic.

3.2 Earth

Figure 7 shows the variations of the Earth's obliquity, i.e. the angle between the Earth's spin axis and the normal of the instantaneous ecliptic, for the last 2 million years. The main period is around 41000 years, which is close to the expected value (cf. Berger & Loutre 1991). The maximum amplitude is around $\pm$1.3 degrees. We have compared our results with the classical findings of Berger & Loutre (1991) and Berger (1992). The curve centered on 22 degrees shows the differences between Berger's and our results. The differences increase considerably for times earlier than 1 Myr before present. A closer analysis reveals that the differences are mainly caused by a temporal phase shift between the two sets of results, while the amplitudes remain quite similar. In a comparison with Fig. 7 of Quinn et al. (1991) the above mentioned discrepancy is much smaller despite the fact that Quinn et al. also included relativistic effects and tides. The curve centered at 21 degrees is the difference between the La93 program of Laskar et al. (1993c) and our findings. In this case relativity is included but not tides (nominal solution). The agreement with our data is also here better than the agreement with

Figure 7: (Top) Obliquity of the Earth over the last 2 Myr. The curve centered at 22 degrees shows the difference between the results of Berger & Loutre (1991) and ours while the curve centered at 21 degrees shows the difference between the relativistic results of Laskar et al. (1993c) and ours. (Bottom) Obliquity of the Earth with and without the Moon.

Berger & Loutre. The differences between our data and those of Laskar et al. are due to relativity and differences in model intricacies and input parameters. The curves indicate that over at least a 1 Myr time-span, the influence of relativity and tides on obliquity can in practise be disregarded. For a more detailed discussion, see Appendix. The bottom part of Fig. 7 shows the variations of obliquity with and without the Moon. As has been pointed out by others (e.g. Laskar et al. 1993a) the Moon is very important in stabilizing the Earth. Without the Moon, the amplitude shows a tenfold increase compared to the case where the Moon is present. Without the Moon, however, the obliquity variations are much slower.

3.3 Mars

We have computed the spin precession rate of Mars, and the results for the last 10 years are presented in Fig. 8. This integration included both Mars moons and was carried out at the extreme time step $\triangle t =3.6525 \times 10^{-3}$ d. The mean precession rate over this period of time is 7.57 arcsec/yr, in excellent agreement with the results given by Folkner et al. (1997) based on the Mars Pathfinder mission.

Figure 8: Instantaneous Mars spin precession rate. A period of 1 marsian year is clearly seen. Variations in amplitude are due to the large eccentricity of the Mars orbit.

As is seen in Fig. 9, however, the precession rate of Mars is not very stable. The precession cycle has varied with around 26000 years, corresponding to about 16%, during the last 1.6 million years. The mean precession cycle length is 167000 yr, corresponding to 7.76 arcsec/yr.

Figure 9: Length of the precession cycles of Mars during the last 1.6 Myr. The dashed line shows the mean value of 167 000 yr.

The inclination of the Mars orbit from 2 Myr before present up to 1 Myr after present is plotted in Fig. 10. There are two main periods - one shorter of around 65000 years and one longer of about 1.3 Myr.

Figure 10: Inclination of the Mars orbit w.r.t. the present orbital plane.

Figure 11 displays the length of successive summer half years of Mars 500 kyr before present (Sect. 2.8). The oscillations show a period of about 50 kyr and the variation of amplitudes is $\pm$15%; since the eccentricity of the Earth orbit is small, this effect for Earth is less significant. The lengths of summers is further discussed in the climate section below.

Figure 11: Length of summers in the northern hemisphere of Mars.

The obliquity variations of Mars are much greater than those for the Earth. The small moons cannot stabilize the spin axis as the Moon does for the Earth. The influence of Phobos and Deimos are still important, however, as can be seen in Fig. 12, which gives the obliquity for the last 500 kyr with and without the moons.

Figure 12: Mars obliquity with and without its moons.

When the moons are included, the obliquity is stretched to the left and the difference between the two curves is seen to become significant. Our curve without the moons is very similar to the curve given by Bouquillon & Souchay (1999, Fig. 5). It is interesting to note that their computation included the moons. One would therefore expect that their curve should match our curve with the moons better. The reason for this discrepancy is that our mean torques are taken over a considerably longer period resulting in much greater mean torques, see Sect. 2.10. This is also demonstrated in Fig. 13. Here we show the tiny difference between the curves with and without the moons for the case where the mean torques are based on a too short time period (10 years). This short period gives an underestimate of the influence of the moons due to the fact that they are at present very close to the equatorial plane of Mars thus giving very small torques. For long simulations it is therefore important to treat the influence of the moons appropriately.

Figure 13: Obliquity difference for Mars. The underestimated mean treatment (10 yr) - isolated Mars. This difference is clearly very small; the correct treatment yields a difference that is 10-15 times larger cf. Fig. 12.

3.4 A climatic application

The usual way for computing astronomical influence on climate is to use the eccentricity, the obliquity, and the longitude of the perihelion to compute the climatic precession parameter $e\sin (\omega )$ and then to compute the mean insolation for various latitudes. Although similar, we have adopted a slightly different approach. The overall climate of Earth is assumed to mainly be governed by events at high northern latitudes. The present land mass distribution dictates that it is for northern latitudes that the continental ice sheets develop. Due to the ocean, the size of Antarctica is assumed to be rather inert in comparison. Below it shall also be argued for the vital role that is played by the summer radiation. Since the planet's spin axis is computed at each time step, the instants for millions of northern summer solstices can easily be identified (see Sect. 2.6). When the time of summer solstice is identified, the obliquity $\epsilon$ and sun-planet distance r are simultaneously recorded. The quantity $P=\sin\epsilon /r^{2}$ is proportional to the summer radiation power received at high northern latitudes. We shall now argue for a model where

Figure 14: (Top) Mean summer solar radiation power (insolation) and differentiated ice volume ( $\left\langle P\right\rangle \propto {\rm d}V_{\rm ice}/{\rm d}t$). (Bottom) Ice volume (Imbrie et al. 1990).

the knowledge of this radiation power (P) will be very useful. According to thermodynamics, the energy dQ flowing into the glacial system causes melting of ice into water of mass dm:

$${\rm d}Q=L_{\rm f} {\rm d}m$$ (33)

where the proportionality constant $L_{\rm f}$ is the heat of fusion. The power is the incoming energy per unit of time:

$${\rm power}=\frac{{\rm d}Q}{{\rm d}t}=L_{\rm f}\frac{{\rm d}m}{{\rm d}t}=L_{\rm f}\rho \frac{{\rm d}V}{{\rm d}t}$$ (34)

where $\rho$ is the density for water and dV is the increase in water volume during a particular interval dt. There are certainly many effects that can melt ice such as the incoming solar radiation power, the oceans and various weather phenomena etc. Suppose that the solar radiation power plays an important role for the heating of the oceans, creation of the weather systems and (after absorption in the atmosphere) directly causing ice melting to take place. Our hypothesis here is that in average for long periods, the various complications smear out and the most important remaining effect that really changes is the solar radiation power. These changes are mainly caused by changes in the Earth's orbit and spin axis. As mentioned above, the value for the solar radiation power P is calculated at each summer solstice. This gives an indication how warm or cold a particular summer was. It is very interesting to look for proportionality by replacing "power'' above with the radiation power P so that $P=\sin(\varepsilon )/r^{2}\propto {\rm d}V/{\rm d}t$. Since dV was taken to be the amount of melted water it is replaced by $-(\rho _{\rm ice}/\rho _{\rm water}){\rm d}V_{\rm ice}$. When the water volume increases the ice volume decreases, thus the negative sign. The interesting proportionality to look for thus becomes: $P\propto {\rm d}V_{\rm ice}/{\rm d}t$. It could, of course, be argued that this is not the radiation power near the surface of the ice sheets. However, it is likely that a relative increase in power above the atmosphere in average leads to a corresponding relative increase in power near the ice sheets. The SPECMAP data (Imbrie et al. 1990) are based on sea sediment measurements in the Northern Atlantic. The resolution in these ice volume data is 1000 yrs. We assume that melting mostly takes place during the summer half years. In Sect. 2.8, it is shown how the length of each summer half year is computed. The total melting period dt is then obtained by summing all these 1000 summer half years. The power P is also averaged for 1000 summer solstices resulting in the appropriate $\left\langle P\right\rangle$. The lower part of Fig. 14 shows the ice volume as a function of time for the last 782000 years. The upper part of the figure shows a comparison between mean summer radiation power $\left\langle P\right\rangle$ and the differentiated ice volume. Here, the curves were brought together and "zoomed'' until they match each other. This simple procedure is presented in Fig. 15 (i.e. to move the curves together and determine the proportionality constant). If the proportionality assumption is perfect the match would be perfect. The top part of Fig. 14 indeed shows that the match is very good and that the proportionality hypothesis is sound. What needs to be further addressed is the perhaps surprising fact that also ice growth is very well predicted. It is evident that during one year we have two competing effects: namely ice growth during the winter half year and melting during the summer half year. If the summer radiation power is high the melting will be high. The corresponding winter (the complement) for that year is probably quite cold and dry meaning that the ice growth should be low. The net result for that year is high melting. On the other hand, if summer radiation power is low the melting should also be low. The corresponding winter will in this case be unusually warm and humid with a high probability of snowfall thus causing an increase in ice growth. The net for this year is ice growth. The true weather conditions are of course much more complicated. However, the resolution for the ice volume data is 1000 years. No matter how complicated a given year actually was, the average melting and growth for 1000 year intervals follows the proposed mechanism very well. One interesting interpretation is worthwhile to mention. If one takes the average of the whole power curve in Fig. 14 (top), one gets a line crossing through the whole power curve. All the crossing points between the power curve and this average power line correspond almost exactly to maxima or minima of the ice volume (bottom curve). In other words, the average summer radiation power corresponds to status quo. Ice melting equals ice growth during a complete year [or rather 1000 years since this is the resolution in ice volume] so the total change in ice volume becomes zero.

Figure 15: This least square fit of summer radiation power and differentiated ice volume determines the proportionality constant between the curves. r=0.76 is the correlation coefficient.

Figure 16: Fine structure of the mean summer radiation power (solid curve) and differentiated ice volume (dotted curve). No $\sim$5 kyr time lag is observed.

Figure 17: Comparison of relative summer radiation power for Mars (top) and Earth (bottom) for their northern latitudes.

Positions of peaks in the top two curves of Fig. 14 are seen to coincide almost perfectly over the entire time interval. The only major deviations are occurring at 10, 130, 340, 430, and 620 kyr before present. A comparison with the ice volume curve shows that these deviations always correspond to extremely rapid terminations of the ice ages. Perhaps this indicates an accelerated absorption of heat due to a reduction in continental ice cover. It is also possible that the initial melting is much more rapid for a fully developed continental ice sheet compared to a later stage when the ice-covered area is smaller. Other complicated reasons involve feedback mechanisms such as greenhouse gases (water vapor, carbon dioxide etc.). This calls for further investigations. It is interesting to observe that the usual frequency problem with the 100000 year period in the ice volume data disappears in our comparison. Neither are there any time lag problems between summer radiation power and ${\rm d}V_{\rm ice}/{\rm d}t$. This problem occurs only when comparing summer radiation power (or insolation) with ice volume ( $V_{\rm ice}$). Figure 16 displays the fine-structure of our Fig. 14. It is evident that there exists no time lag of about 4-5 kyr as reported e.g. by Berger et al. (1993).

Finally, we compare the variations in summer radiation power for Mars and the Earth. The results are displayed in Fig. 17. The huge variations ($\pm$50%) for Mars are obvious for the northern parts of Mars. The causes involve both large variations in obliquity and, due to the comparatively large eccentricity of the Mars orbit, large variations in sun-planet distance at summer solstice. The northern polar regions of the ancient Mars were probably subject to similarly large climatic variations that may have had a tremendous impact on glaciers, atmospheric pressure and the planet's ability to sustain life.