New Long-Term Climate Oscillations

Joseph J. Smulsky

Institute of Earth’s Cryosphere, Tyum. SC of SB RAS, Federal Research Center, Tyumen, Russia.

DOI: 10.4236/ns.2021.138028                

Abstract

The astronomical theory of climate change is based on the solution of differential equations describing Earth’s orbital and rotational motions. The equations are used to calculate the change in insolation over the Earth’s surface. As a result of the author’s solution of the orbital problem, the periods and amplitudes of Earth-orbit variations and their evolution have been refined. Unlike previous studies, the equations of Earth’s rotational motion are solved completely. The Earth’s rotational axis precesses relative to a direction different from the direction of the orbit’s axial precession, and oscillates with periods of half a month, half a year and 18.6 years. Also, its oscillations occur with irregular periods of several tens of thousands of years and more. All these motions lead to oscillations of the obliquity in the range of 14.7° to 32.1°, which prove to be 7 - 8 times larger than obtained by a previous theory. In the same proportion, the Earth’s insolation oscillations increase in amplitude, with insolation extremes occurring in other epochs than those in the previous theory. The amplitudes and the onset times of the extremes correlate with known paleoclimate changes. Thirteen insolation periods of paleoclimate variation over an interval of 200 thousand years are identified. From the insolation evolution calculated over a time interval of 1 million years, 6 climate gradations from very cold to very warm are identified.

Keywords

Earth’s Orbit and Axis Evolution, Insolation, Climate Change, Ice Ages

1. INTRODUCTION

Long-term climate oscillations are analyzed in the Astronomical theory of climate change (or alternately named, Astronomical theory of ice ages). The theory is based on the solution of the following three problems: 1) what are the changes in the Earth’s orbit? 2) what are the changes in the Earth’s axis of rotation? 3) what are the changes in the amount of solar radiation over the Earth’s latitude based on the first two changes? The original version of the theory was developed by M. Milankovitch [1] in the first quarter of the 20^th century. Subsequently, the proposed approach was improved by other researchers [2 - 6]. At the end of the 20^th century, activities aimed at revisiting the above problems were initiated [7]. As the result of a more precise solution to the problem of Earth’s rotational motion, a second version of the climate change theory succeeded in explaining the factors that cause long-term climate oscillations.

This paper presents the main results of the new Astronomical theory of climate change, how they were obtained and their reliability.

2. METHOD

This paper presents the results of decades of theoretical studies. The foundations of mechanics were analyzed, unsubstantiated hypotheses were eliminated and mechanical operations were formulated with regard to their actual content [8]. The differential equations of the orbital motion of planets in the Solar system were redefined, a new method for solving them was developed, and the Galactica program was created for their numerical solution [9 - 11]. Differential equations for orbital and rotational motion were derived in a new way and a method for analyzing, solving and integrating them was developed. This evolved into a program for their numerical solution [11 - 13]. Orbital motion was considered in the coordinates associated with the fixed equator of Earth, and rotational motion was investigated in coordinates associated with a fixed plane of the Earth’s orbit. Transformations of the parameters of orbital motion relative to the coordinates associated with the plane of the moving equator are derived [14]. The equations for the irradiation of the Earth by the Sun, depending on the parameters of its orbital and rotational motions, are derived in a new way. A program for calculating the Earth insolation has been developed [15]. The numerical integration of these equations on supercomputers has been performed. These tasks are complex and laborious. For example, the integration of the equations of orbital motion over 100 million years took two years of machine time. The analysis of the solutions of these equations was carried out over several years.

These studies are unique, so all their stages were tested. Algorithms for checking the accuracy of solutions were introduced into the programs for integrating equations. For example, there are about two dozen in the Galactica program. All stages of the solution were checked, compared with the observations and conclusions of other researchers. Checking took longer than the solutions themselves. For example, verification of solutions to a problem on the rotation of the Earth was carried out over a period of three years.

3. RESULTS

3.1. Earth’s Motions and Their Variations

The Earth moves in an elliptical orbit around the Sun, which is located at one focus of an ellipse (Figure 1). The shortest Earth-Sun distance in perihelion is denoted by R_p, and the largest distance in aphelion by R_a. The period of Earth’s motion with respect to motionless space connected with the Solar system is P_sd = 365.25636042 days [16 , 17]. The quantity P_sd is called the sidereal period of the Earth’s rotation around the Sun. The Earth’s orbital motion proceeds in an anticlockwise direction, based on the orbit being viewed from the Earth’s North Pole N. The normal to the orbital plane is denoted as $\stackrel{\to }{S}$ , and is called the orbital axis.

With respect to motionless space, the Earth rotates around its axis, $\stackrel{\to }{N}$ at an angular velocity of ω_E = 7.292115∙10⁻⁵ 1/sec in an anticlockwise direction coincident with the direction of the Earth’s orbital motion. The value of ω_E corresponds to a full revolution performed by the Earth in 0.99726968 days. The Earth’s rotational axis $\stackrel{\to }{N}$ is inclined to the orbital axis $\stackrel{\to }{S}$ at an angle equal in the contemporary epoch, to ε = 23.43˚. This inclination is called the obliquity. During the orbital motion of Earth, the orientation of its rotational axis $\stackrel{\to }{N}$ remains unchanged in space (Figure 1). That is why at two points of the orbit, at times March 20, (20.03), and September 22, (22.09), the axis $\stackrel{\to }{N}$ is normal to the Earth-Sun direction. With respect to Earth, the Sun is in its equatorial plane. That is why the southern and northern hemispheres receive identical amounts of solar radiation and the day appears to be equal in duration to night. These points are called the day of vernal equinox, (20.03), and the day of autumnal equinox, (22.09). At the time, June 21, (21.06), the axis $\stackrel{\to }{N}$ is least inclined to the Earth-Sun line, and the northern hemisphere is therefore illuminated with maximum solar radiation. At the time of December 21, (21.12), the axis, $\stackrel{\to }{N}$

Figure 1. The Earth’s position in its orbit in 2025 during the days of vernal equinox (20.03), summer solstice (21.06), autumnal equinox (22.09), and winter solstice (21.12), and the time expressed in Earth motion during the days in spring (92.7 d), in summer (93.7 d), in autumn (89.9 d), and in winter (89.0 d): $\stackrel{\to }{N}$ is the Earth’s axis of rotation; ${\stackrel{\to }{M}}_2$ is a vector relative to which the axis $\stackrel{\to }{N}$ precesses in a period of 25.74 thousand years; $\stackrel{\to }{S}$ is the Earth’s rotational axis, and $\stackrel{\to }{M}$ is a vector relative to which the $\stackrel{\to }{S}$ axis precesses over a period of 68.7 thousand years [14 , 15 , 18 , 19].

is at maximum inclination to the Earth-Sun line. This is the reason the southern hemisphere is at maximum illumination at that time, and polar night approaches at high latitudes in the northern hemisphere. Since the time spent on reaching and leaving the extreme angles lasts for several days, these points are called respectively, the summer solstice day (21.06), and the winter solstice day (21.12).

The inclination of Earth’s axis, $\stackrel{\to }{N}$ to the orbital axis, $\stackrel{\to }{S}$ leads to the variation in duration of sunlight, both during the year and on the same day at different latitudes. On summer solstice day, (Figure 1, 21.06), we have a polar day in the whole region between the North Pole and the Arctic Circle. Then, as the latitude decreases, the day gets shorter, reaching a 12-hour duration at the equator, and we have a polar night established below the Antarctic Circle. Contrary to this, on the day of winter solstice, (21.12), in the territory between the North Pole and the Arctic Circle we have a polar night. Then the day starts increasing in duration. At the equator, the day lasts for 12 hours, and a polar day sets in below the Antarctic Circle. As we approach the equinoctial points of 20.03, and 22.09, the difference between the days in latitude decreases in value. The day’s duration along all latitudes becomes identical at 12 hours.

As the Earth moves along its orbit, the alteration of seasons occurs. The duration of the seasons is defined by the Earth’s motion over certain orbital segments. From the vernal equinox day, 20.03, until the summer solstice day, 21.06, the duration of spring is 92.7 days. Over the summer segment, the duration is 93.7 days. Over the autumnal segment, the duration is 89.9 days. Over the winter segment, the duration is 89.0 days.

The Earth’s orbital and rotational motions define the variation in climate in the current epoch. However, the motions vary in time, and the climate undergoes change. The position of Earth’s orbit precesses in space. The Earth’s orbital axis $\stackrel{\to }{S}$ (Figure 1), rotates, or in other words, it precesses about the direction of $\stackrel{\to }{M}$ , which is motionless in space. The precession proceeds clockwise with a period of 68.7 thousand years. In a clockwise direction, the Earth’s axis $\stackrel{\to }{N}$ precesses about the direction of ${\stackrel{\to }{M}}_2$ , also motionless in space. The precession period is 25.74 thousand years. Besides this, the axes $\stackrel{\to }{S}$ and $\stackrel{\to }{N}$ execute oscillations, each with respect to its own precessional axis $\stackrel{\to }{M}$ and ${\stackrel{\to }{M}}_2$ , respectively. In addition, the shape of the orbit (its eccentricity, whose value varies from 0 to 0.064; the current value being equal to 0.016) and the perihelion position, both undergo variations. Today, the perihelion is over the winter segment (Figure 1), when winter sets in the northern hemisphere. Since the Earth’s orbital perihelion rotates in anticlockwise direction at a mean period of 147 thousand years, its position in other epochs can be at any point in the Earth’s orbit.

The changes of the Earth’s orbital and rotational motions lead to changes in its climate. Below, the latter changes are considered in more detail.

3.2. Evolution of the Earth’s Orbital Motion

The evolution of the Earth’s orbital motion is considered in a motionless frame, xyz whose origin is located at the center, O of celestial sphere 1 (Figure 2). Note that, depending on the particular problem of interest, the point O can be located either at the center of mass of the Solar system, at the center of the Sun, or at the center of the Earth. As a result of the interaction of Solar-system bodies, the Earth’s equatorial plane 2 and its orbital plane 3 both alter their positions, which are denoted with digits 2 and 3 in the epoch Т₀, and with digits 4 and 5 in the epoch Т, respectively. At epoch Т₀, we consider the epoch of the year 2000.0.

Since the annual motion of the Sun over celestial sphere 1 with respect to the Earth, proceeds along circles 5 or 3, the planes of the circles are additionally called, the planes of moving and motionless ecliptic, respectively. The frame xyz is related to the plane of the motionless Equator 2. The moving plane 5 of the Earth’s orbit is defined by the angle φ_Ω = γ₀γ₂ of the ascending-node position γ₂ and by the inclination of the plane i.

The Earth moves around the Sun along an open trajectory which is close to the shape of an ellipse. At one point in the trajectory, the perihelion, the Earth approaches the Sun to the shortest distance R_p, and at the opposite point, the aphelion, it moves away from the Sun to the largest distance R_a. In Figure 2, the projection of the perihelion onto the celestial sphere 1 is denoted with the character B and its position is coordinated with the angle φ_p = γ₂B. The shape of the orbit is defined by its eccentricity $e=\left(R_a-R_p\right)/\left(R_a+R_p\right)$ .

The oscillations of the angles φ_Ω and I reflect the rotation process of the orbit’s axis $\stackrel{\to }{S}$ (Figure 2), that is, the normal to the orbit plane, with a period of 68.7 thousand years about the motionless vector $\stackrel{\to }{M}$ [20]. The latter vector is the sum of the angular momentum vectors of all bodies in the Solar system. The plane normal to $\stackrel{\to }{M}$ is called the Laplace plane. The angles that define the position of the Laplace plane with respect to the motionless equatorial plane xOy in radians are equal to i_M = 0.401834 andφ_M = 0.0680946. The rotation or, in other words, the precession of the axis $\stackrel{\to }{S}$ about the vector $\stackrel{\to }{M}$ proceeds

Figure 2. Parameters of the Earth’s orbit and axis in the motionless equatorial xyz and motionless ecliptic x_ey_ez_e coordinate systems. 1—the celestial sphere; the planes in the epoch of Т₀: 2—that of the Earth’s equator, 3—that of the Earth’s orbit; in the epoch of Т: 4—that of the Earth’s equator, 5—that of the Earth’s orbit; unit vectors: $\stackrel{\to }{N}$ —Earth’s axis, $\stackrel{\to }{S}$ —Earth’s orbit axis, $\stackrel{\to }{M}$ —the angular momentum of the Solar system; γ₀—the vernal equinoctial point in the epoch of Т₀; B—the perihelion’s position on the celestial sphere; φ_Ω = γ₀γ₂—the angle of the orbit’s ascending node; φ_р = γ₂B—the perihelion angle; i—the orbit’s inclination.

in clockwise direction. In addition to the precession, the orbit’s axis $\stackrel{\to }{S}$ executes oscillations at periods of 97.35 thousand years, 1.164 million years, and 2.32 million years, relative to the vector $\stackrel{\to }{M}$ . During those oscillations, the angle between the vectors $\stackrel{\to }{S}$ and $\stackrel{\to }{M}$ never exceed a value of 2.94˚. All those motions are manifested in the behavior of angles φ_Ω and i in Figure 3.

Thus, the evolution of the Earth’s orbit proceeds as a result of the following four motions: 1) precession of the orbit’s axis $\stackrel{\to }{S}$ ; 2) oscillations of the orbit’s axis $\stackrel{\to }{S}$ ; 3) oscillations of the orbit’s eccentricitye; and 4) rotation of the orbit in its own plane (perihelion rotation).

Figure 3 shows the evolution of the Earth’s orbital parameters e, φ_Ω, i and φ_p over a span of one million years. The shortest eccentricity oscillation period is T_e1 = 94.5 thousand years, and the longer ones are the periods T_e2 = 413 thousand years andT_e3 = 2.31 million years [19 , 20]. Both the angle, φ_Ω of the orbit’s ascending node and the orbit’s inclination, i oscillate with a period T_Ω = T_i = 68.7 thousand years. The increasing behavior of the angle, φ_p toward future (Figure 3) reflects the non-uniform rotation of the perihelion in anticlockwise direction with a mean period of T_p = 152 thousand years (over a time interval of one million years). As it is seen from the graph, there exists an epoch of a reverse, or clockwise, perihelion’s motion.

3.3. Evolution of the Earth’s Rotational Motion

Evolution of the Earth’s rotational motion is treated in the motionless frame x_ey_ez_e (Figure 2) connected

Figure 3. Evolution of Earth’s orbital parameters over a time interval of 1 million years: е—eccentricity; φ_Ω—angle of the orbit’s ascending node; i—orbit’s inclination; φ_р—perihelion angle; Т—time in million years reckoned from December 30, 1949; Т_е1, Т_Ω, and Т_i—the least oscillation periods for the eccentricity, the ascending node, and the orbit inclination expressed in a thousand years; Т_p—perihelion rotation period averaged over a time interval of 1 million years.

with the Earth orbit’s motionless plane, 3. The inclination angle and the precession angle ψ = γ₀γ₁ both define the position of equator, 4 moving with respect to the orbit’s motionless plane, 3. The precession angle ψ decreases toward future in an oscillatory manner at a mean rate of ${\stackrel{\dot{}}{\psi }}_m=-\text{2}\cdot \pi /P_{pr}$ according to a linear law. Here P_pr = 25.738 thousand years is the period of the Earth axisial precession averaged over a time interval of 1 million years [7 , 13 , 14]. The negative sign of ${\stackrel{\dot{}}{\psi }}_m$ implies that the unit vector $\stackrel{\to }{N}$ of the Earth’s rotation axis precesses in a clockwise direction (Figure 2). Precession of the Earth’s rotational axis, $\stackrel{\to }{N}$ proceeds around the second motionless vector ${\stackrel{\to }{M}}_2$ ; the angles defining the orientation of the plane normal to this vector with respect to the plane xOy are i_M2 = 0.417728 andφ_M2 = −0.0664662. The angle between the motionless vectors ${\stackrel{\to }{M}}_1$ and ${\stackrel{\to }{M}}_2$ is 3.201402˚.

The variation of the difference, $\Delta \psi =\psi -{\psi }_a$ for the precession angle, ψ (here, ${\psi }_a={\psi }_0+{\stackrel{\dot{}}{\psi }}_m\cdot T$ is the change of the precession angle in the process proceeding at the mean rate ${\stackrel{\dot{}}{\psi }}_m$ ) over a time interval of 1 million years, is shown in Figure 4. Over the latter time interval, the quantity, Δψ varies from −0.184 to 0.233 radian, so that the full oscillation swing amounts to 0.417 radian.

The inclination difference, $\Delta \theta =\theta -{\theta }_0$ in Figure 4, is given with respect to the initial angle θ₀ = 0.40904645. The quantity Δθ oscillates similar to, Δψ, yet in a narrower range from −0.0845 to 0.0855, so that the oscillation range here amounts to 0.17 rad. Thus, the oscillation amplitude for angle θ is 2.45 times smaller than that for angle ψ. Besides, the oscillations of Δθ do not coincide in phase with the oscillations of Δψ; they are shifted along the time axis T by −7.5 thousand years.

3.4. Evolution of the Earth’s Orbital Motion Relative to Its Rotational Motion

The orbital-motion parameters i, φ_Ω and φ_p and the rotational-motion parameters ψ and θ are defined by the obliquity, ε and perihelion angle φ_pγ of the orbit’s moving plane 5 with respect to the moving equator 4 (Figure 2). The oscillation spectrum of φ_pγ is rather broad since the angles φ_p, φ_Ω, i, ψ andθ contribute to the oscillations of this angle. The average variation of the angle, φ_pγ proceeds according to the law φ_pγm = φ_p −(2π·T/P_pr).

Figure 5 shows the variation of ε over five different time intervals In [13]. Over short time intervals, the oscillations of θ and ε are roughly identical. Indicated in the graphs are the main oscillation periods T_ni and amplitudes (θ_ai and ε_a4) of the inclination angle: half-month period T_n2 and half-year periods T_n3 and T_n4 = 18.6 years. Those oscillations are called nutation oscillations. The precession angle ψ exhibits similar oscillation periods, the amplitudes being two or three times greater.

Over the time interval In = 0.1 year, half-monthly oscillations are observed and diurnal oscillations with period T_n1 = 0.9973 day can be traced; over the interval In = 1 year, half-year oscillations emerge; over the interval In = 10 years, an oscillatory trend with a T_n4 = 18.6-year period is observed, with oscillations at this frequency prevailing over the time interval of In = 100 years.

Over the time interval In = 100 years, it is seen that the calculated obliquity ε 1) oscillates about its mean value, 2) received by S. Newcomb [16] and J. Simon et al. [17]. The oscillation amplitude ε_a4 = 9.2" of the period, T_n4 = 18.6 years also coincides with the observations. In astronomy, this quantity is called the nutation constant.

Figure 4. Evolution of the Earth’s rotational motion over a time interval of 1 million years. The differences; the precession Δψ and inclination Δθ angles are given in radians.

Figure 5. The dynamic of the obliquity ε (in radians) over five time intervals In: yr—year; Δθ ≈ ε − ε₀; ε₀—the obliquity at the initial epoch of December 30, 1949. T_n2, T_n3, T_n4 and θ_a2, θ_a3, ε_a4—the oscillation periods and amplitudes of the inclination angles; 1—according to the solutions of [7 , 12 , 13]; 2—approximation of observation data according to S. Newcomb [16] and J. Simon et al. [17]; 3—according to the solution by J. Laskar et al. [6]; 4—according to the solution by Sh. G. Sharaf and N. A. Budnikova [3].

3.5. Evolution of the Earth’s Obliquity and Insolation over a Span of 1 Million Years

As it is evident from Figure 5, over the time interval of In = 10 thousand years, a coincidence of the new obliquity ε 1 with the data 2 and 3 yielded by the first version of the Astronomical theory of climate change [1 - 6] is observed over a span of 2000 years.

It should be noted that these authors solved the problem in question over a large time interval: Sharaf and Budnikova [3] for 30 million years, and Laskar et al. [6] for an even longer period. We showed [18] that their results do not fundamentally differ from the results of Milankovitch [1], Berger and Loutre [4], Edvardsson et al. [5]. Therefore, we compare our results with these authors who are typical representatives of the former Astronomical theory.

As can be seen from Figure 5, after 2000 years, obliquity, ε calculated within the new version 1 of the theory shows clear deviation. As it is seen from the graphs of Figure 6, over the time interval of 1 million years the oscillations of ε as yielded by the second version of the theory proceed in the range of 14.7˚ to 32.1˚, whereas the same range in the previous theory was from 22.08˚ to 24.45˚; in other words, the range of oscillations in the second version of the theory proves to be seven times greater.

This difference is due to the fact that in the second version of Astronomical theory, the Earth’s rotation problem was treated in full, without simplifications. The solution of this problem and various checks of obtained data were analyzed at length in the publications [7 , 13 , 21] and will be covered in the subsection 4.3.

The amount of solar radiation reaching the Earth’s surface, also called the Earth’s insolation, is defined by the parameters e, ε and φ_pγ. Figure 6 gives a comparison of the insolation changes $Q_s^{\text{65}N}$ occurring during the summer caloric half-year at the 65-deg northern latitude in the second version of the theory, (curve 1) [21] with the changes as calculated by the first version of the theory [6]. Here, the amplitude of insolation oscillations is also seven times greater than that in the previous theory. Besides, the insolation extremes occur at other times, and the oscillation periods are different. Note that the astronomical

Figure 6. Evolution of the obliquity ε and that of summer half-year insolations $Q_s^{\text{65}N}$ and I over a time interval of 1 million years. Comparison of results yielded by the second version of the Astronomical theory of climate change (curve 1), with the results of the first version of this theory (curve 2), demonstrated using, as an example, the work by J. Laskar et al. [6]. $Q_s^{\text{65}N}$ —insolation in GJ/m² over the summer caloric half-year at 65-deg northern latitude; I—insolation at the equivalent latitudes over the summer caloric half-year at 65-deg northern latitude. Indicated in degrees are the maximum and minimum values of ε.

summer and winter half-years measured from the vernal equinox day to the autumnal equinox day and visa versa differ in duration for different epochs. That is why it is caloric half-year, equal in duration, that are considered here.

In order to compare climates in other epochs with the current climate, we consider the insolations at equivalent latitudes I. For calculating of I, we consider the Earth’s latitude φ characterized by receiving same amount of summer solar radiation, Q_s as in the current epoch. Figure 6 shows the insolation oscillations at equivalent latitudes I over a time interval of 1 million years. The lowest values I ≈ 90˚ indicate that at latitude 65˚N in summertime, there was less solar radiation on the pole than now. The highest values, such as I ≈ 23˚ at the time −0.031 million years denote epochs in which, in summertime, the amount of solar radiation having reached the Earth at latitude 65˚N exceeds the amount of solar radiation having fallen onto it presently in the tropics, i.e. in the equatorial area. Such profound insolation oscillations lead to substantial climate oscillations. As it is seen from curve 2, the oscillations of I in the previous theory were less significant.

3.6. Variation of Insolation over the Earth’s Latitude

Figure 7 shows the variation over the latitude φ of the year Q_T, and during caloric half-years summer Q_s, and winter Q_w insolations in three epochs: in contemporary epoch T = 0, in the warmest epoch T = −31.28 kyr and in the coldest epoch T = −46.44 kyr (over a time interval of 200 thousand years). Those epochs are characterized by the following values of 65˚N summer insolation: $Q_s^{\text{65}N}=5.9$ , 7.4, and 4.7 GJ/m², respectively. In those epochs, the obliquities were ε = 23.44˚, 32.10˚ and 14.8˚, respectively.

The summer insolation Q_s (Figure 7, dashed lines) in contemporary epoch, 1 exhibits a minimum value on the poles, and it reaches a maximum value at the tropics φ = ε, and in the vicinity of the equator, this insolation exhibits a minimum. As we move from the cold (line 3) to the warm epoch 2, the summer insolation Q_s on the poles increases by a factor of 2.07. At latitude 65˚N this insolation increases by a factor of 1.57. Since, on the average, this latitude well represents the variation of insolation at high latitudes, the above insolation was adopted by Milankovitch [1] as reference one in the climate characterization procedure. In warm epoch 2, the summer insolation Q_s has an equatorial minimum in the southern hemisphere, and in cold epoch 3, in the northern hemisphere.

The winter insolation Q_w (Figure 7) on the poles is zero, and it monotonically increases in the equatorial region. In the equatorial region, the insolation Q_w exhibits a maximum at a latitude, φ at which the summer insolation Q_s shows a minimum. Over the period from cold epoch 3 to warm epoch 2, the winter insolation Q_w exhibits most pronounced variations at middle latitudes. In the latter situation, for epochs 2 and 3 under consideration, e.g. at the latitude φ = 40˚, the change of the winter insolation is 1.38 times greater in the northern hemisphere than in the southern hemisphere. In cold epoch 3, the winter insolation at all latitudes is greater than that in warm epoch 2. In other words, during the cold epochs the winter seasons are warmer than those during the warm epochs.

The annular insolation Q_T (Figure 7) monotonically increases from the poles toward the equator. At the equator, the annular insolation exhibits a maximum, with the annular insolation being symmetrical with respect to the equator. In other words, the amounts of heat per year are identical in both hemispheres. From cold epoch 3, to warm epoch 2, the annular insolation Q_T on the poles increases by the same factor as the summer insolation Q_s. With decreasing latitude, the difference between the annular insolations decreases,

Figure 7. The distribution over latitude φ of insolations over three epochs. (Q_s, is summer half-year insolation; Q_w is winter half-year insolation; Q_T, is annular insolations; 1 is contemporary epoch; 2 is the warmest epoch and 3 is the coldest epoch over a time interval of 200 thousand years; $Q_s^{\text{65}N}$ is insolation in GJ/m² over the summer caloric half-year at 65-deg northern latitude; T, kyr—time in thousand years).

and at the latitude, φ = 45˚ the annular insolation experiences no changes. In the equatorial region, the changes of Q_T are reciprocal to its changes at the high latitudes: in cold epoch 3, the amount of heat per year exceed that in the warm epoch. In the latter situation, the change of insolation Q_T is four times smaller than that in the high-latitude region. That is why the main changes of the annular insolation occur at high latitudes.

3.7. Periods and Gradations of Earth’s Climate Changes

Over the previous interval of 200 thousand years (see Figure 8), 13 climatic periods O_I, 1_I, 2_I,12_I were identified [19 , 22]. As a result of the comparison of these periods with paleoclimate data for Western Siberia over 50 thousand years, it was found that the periods 3_I, 2_I, 1_I, O_I refer respectively to the Ermakov ice age, Karginsky warming, Sartan glaciation, and Holocene optimum. Those events also correspond to ice ages and interglacial periods in Europe and North America.

Also, the following gradations of the warm and cold climate were introduced (Figure 8): moderately warm, warm, and extremely warm climate levels, and moderately cold, cold, and extremely cold climate levels. During the past period of 1 million years (see Figure 9), the Earth has experienced six extremely cold (e.c.) periods and four extremely warm (e.w.) periods. The total number of cold (c.) and warm (w.) periods was 16 warm periods and 16 cold periods. Other periods were moderately cold (m.c.) ones and

Figure 8. Insolation periods O_I, 1_I, 2_I, … 12_I over a time interval of 200 thousand years and their boundaries: 1—mean insolation Q_sm; 1t and 2t—the first or second boundaries of warm levels; 1c and 2c—the first and second boundaries of cold levels; m.w., w., e.w—moderately warm, warm, and extremely warm levels; m.c., c., e.c. —moderately cold, cold, and extremely cold levels.

Figure 9. Climate levels over a time interval of 1 million years: 1—mean insolation Q_sm; 1t and 2t—the first and second boundaries of the warm levels; 1c and 2c—the first and second boundaries of the cold levels; m.w., w., e.w.—warm climate levels; m.c., c., e.c.—cold climate levels. m. cl.— moderate-climate periods.

moderately warm (m.w.) ones. Besides, there were nine moderate-climate (m.cl.) periods, which included both cooling and warming phases.

4. DISCUSSION

4.1. Statements of the Problems and the Differences among Them

In the previous Astronomical theory of climate change, variations of Earth’s orbital elements were received on the basis of the theory of secular perturbations. This is an approximate analytical method for solving the interaction problem for Solar-system bodies. In that theory, the changes of the equator plane 4 (Figure 2) were analyzed approximately. With respect to this plane, the angle ε of the inclination of the orbital plane 5 and the perihelion position (being, in our notation, the angle φ_pγ in Figure 2) were determined.

In the new Astronomical theory the interaction problem for Solar-system bodies (for simplicity, we will call this the orbital problem) was solved, using no simplifications, with a high-precision numerical method using a specially developed Galactica system [9 - 11]. In the latter case, we consider a change in the orbital plane (item 5 in Figure 2) relative to the fixed space represented by the equatorial plane 2 on a certain epoch. Here, the inclination, i of the orbital plane 5 and the ascending-node angle φ_Ω differ from the angles figuring in the theory of secular perturbations. As it was noted above, in this theory the angle ε between the moving equatorial plane 4 and the moving orbital plane 5 was used.

Unlike in the theory of secular perturbations, we also numerically solved a second problem; the one about the Earth’s rotation, governed by its own differential equations. As a result, we obtain the laws of variation of the inclination angle θ and the precession angle ψ of the moving equatorial plane 4 with respect to the motionless orbital plane 3 (Figure 2).

A third complex geometric problem was then solved analytically for determining both the obliquity ε of the moving equatorial plane 4 with respect to the moving orbital plane 5 and the perihelion angle φ_pγ. The solution of the problem over short time intervals of the order of one thousand years is available in celestial mechanics. Over time intervals of several million years, during which those two planes and the orbit perihelion executed many irregular rotations in both directions, no solution was known.

A fourth problem, which yields the variation of the insolation versus the change of Earth’s orbital and its axis parameters, thus offering an insolation theory, was presented in its complete form by M. Milankovitch in the first quarter of the twentieth century. We have also solved this problem in a new way [15]. Here, a new mathematical algorithm for elliptic motion, more understandable to non-specialists and useful for computer calculations, was employed.

All equations, including the differential equations for the orbital and rotational motions, were derived using a novel approach. Since the resulting data were other quantities in different representations, new methods for their analysis were developed. All that activity was accompanied with the development of computer codes written in various programming languages.

4.2. Adequacy of the Solution of the Orbital Problem

In connection with the new solutions of all four problems, at due stages, all the problems were subjected to checking. For checking the adequacy in the solution of the orbital problem, nine trustworthiness criteria were developed. Some of those criteria were included into the Galactica program; that is why control was exercised right in the course of solving the problem. For all bodies having an observation base (these are the planets from Mercury to Neptune, and also the Moon, the solutions over the time interval of several thousand years were compared with the secular variations of orbital parameters. The coincidence was found to be excellent [20 , 23].

Over intervals of hundred thousand and million years, the orbital parameters were compared to the results obtained by the previous researchers [3 , 4 , 6]. The data were also found to be coincidental. Each of subsequent authors took into account the experience gained by previous researchers and made the theory of secular perturbations more refined. The theory was confirmed by making necessary comparisons. The later the works were published, the longer was the time interval over which the corresponding results were coincident with our data [20].

As it was noted previously, the perturbation theory is an approximate method for solving the orbital problem. After 20 million years, the solution obviously started departing from the actual data: the orbits of individual planets started increasing in size and, later, the planets could leave the Solar system [24]. We have solved the orbital problem over a time interval of 100 million years. All the orbital parameters of the planets and Moon executed steady oscillations, and there existed no tendency toward the change of those oscillations [20].

4.3. Adequacy of the Solution of the Earth’s Rotation Problem

Points concerning the reliability of the solution of the Earth’s rotation problem were analyzed in detail in publications [7 , 13 , 14]. Within the adopted solution method, all necessary checks were performed. For instance, the problem was solved in succession, implying the action of one of ten bodies (the Sun, the eight planets, and the Moon) [12]. The obtained oscillation periods of the Earth’s axis were confirmed with general conclusions made on the basis of the theorem of change in angular momentum and also, with the results of other authors [25]. With the actions due to all bodies, the problem was solved over different time intervals; as it was shown previously in sec. 4 and sec. 5, the obtained data proved to be coincident with observations. Integration of the equations over a time interval of 200 thousand years was performed with different initial conditions and integration steps. The latter resulted in no changes of the oscillation periods, oscillation amplitudes, and the onset times of the extremes.

From the graph with the interval I_n = 100 yr in Figure 5, one can see that the middle of the oscillations of the calculated 1 obliquity ε coincides with the observed average angle 2. That average angle 2 complied with all ancient observations made over a period of 2.5 thousand years of their history. The obtained oscillations with 9.2’’ amplitude and 18.6-year period were found to be precisely coincident with the observed oscillations. In the graph of Figure 5, with the interval I_n = 10,000 yr., it is seen that the deviation of the calculated angle, ε from the linear trend established for observational data, starts manifesting itself after 2.5 thousand years. In Figure 5, solutions for future time intervals are shown. The solutions for previous time intervals have a similar general appearance.

Following 2.5 thousand years, solutions obtained by previous researchers also started departing from the linear trend. From this time on, differences between solutions 2 and 3 by those authors and our solution 1 was observed. Over a longer time interval, namely, 200 thousand years, the Earth’s rotation problem was solved initially towards the future [26]. The solution for, ε was found to feature a different oscillation structure and other onset times for extremes, and what was most important, the amplitude of new oscillations exceeding the amplitude of the previous oscillations by a factor of 7 - 8 was revealed. From this time on, a solution check for the Earth’s rotation problem was initiated. This check lasted for three years.

The Earth’s rotation problem is one of the most complex problems in mechanics. Its solution can depend on the fundamental assumptions made while deriving the equations, on the choice of initial conditions, and on the procedure of solution reduction to the Earth’s moving orbit. That is why a cardinal check of obtained results would be their obtaining, without solving differential equations for rotational motion.

While studying the orbits of the bodies, we found that the evolution of the Moon’s orbital axis was similar to that of Earth’s axis of rotation. That result had led us to a compound model for Earth rotation in which, part of the Earth’s mass was uniformly distributed among peripheral bodies rotating around a central body, along a circular orbit. Under the action of the Moon, the Sun, and the planets, the orbits of the peripheral bodies started changing. It should be noted that under the axis of the orbit is meant a perpendicular to its plane. The evolution of the orbital axis of one of the bodies, modeled the evolution of the Earth’s axis. Such modeling of the Earth’s rotational motion included performing several solution stages of the orbital problem being treated with the Galactica program. In the initial series of our studies [20 , 27], three models were studied, and the possibility of modeling the evolution of the Earth’s axis was confirmed. In those models, the precession periods of the orbit axes were 170 and 2604 years, whereas the average period of the Earth’s precession axis is known to be 25,740 years. Was subsequently yet developed 11 models, while has not been reached the required period of precession [7 , 13]. In the 13^th model, the orbital radius of the peripheral bodies was equal to the Earth’s radius, the rotation period of the bodies was 0.142 hour, and the interaction between the model’s bodies was amplified by a factor of 9.6 in comparison with the gravitational interaction. Thus, the bodies of the 13^th model rotated 170 times faster than the Earth does. For studying the evolution of such models, the Galactica program was developed; in this model, the possibility of changing the interaction between the model’s bodies was provided.

The solution of the problem for the 13^th compound model of the Earth, over a time interval of 300 years [7 , 13], has yielded all characteristics of the dynamics of the Earth’s axis, including the half-monthly and half-year oscillations of the angles ε and ψ and the oscillations of those angles with a 18.6-year period. The amplitudes of those oscillations have also proved to be coincident with the solution results of the Earth’s rotation problem. Such a coincidence in the results of the model problem with the results of the direct problem occurs once over a period of 3 thousand years. Further solution of the problem is hampered by the necessity to reduce the integration step to values for which the computing time turns out to be too long. Thus, over a time interval of 3000 years the compound model of the Earth has confirmed the obtained results of integration of the differential equation for Earth’s rotation. The latter indicates that the assumptions and simplifications adopted in the derivation of the equations, the derivation of the equations itself, the solution method, and the transforming of the integrated data into the final form have also been confirmed.

A second check consisted in using an alternative integration method. In the program DfEqAl1.for solving the Earth’s rotation problem, a fourth-order Runge-Kutta integration method as implemented by Krut’ko et al. [28] was used. Over a time interval of 200 thousand years, a growth of daily oscillations of the derivatives $\stackrel{\dot{}}{\psi }$ and $\stackrel{\dot{}}{\theta }$ was identified. Next, a code for solving the DfEqADP8.for program was developed, which uses the Dormand-Prince method, i.e. the eight-order Runge-Kutta integration method [29]. In integrating the equations of rotational motion over a period of 200 thousand years all previously obtained results have found confirmation. In the latter case, the amplitude of the daily oscillations of $\stackrel{\dot{}}{\psi }$ and $\stackrel{\dot{}}{\theta }$ showed no increase and remained at one and the same level. Thus, the particular method for integrating the equations does not affect the obtained results, and the application of a more accurate method confirms them.

A third check consisted in employing another method for solving the problem. The differential equations for rotational motion contained the coordinates of ten bodies acting on the Earth (eight planets, the Sun, and the Moon). Those coordinates were determined while solving the orbital problem with the Galactica program. However, the data array over large time intervals obtained with an integration step equal to 10⁻⁵ - 10⁻⁴ years took an overly large memory space. That is why we have developed a mathematical model for the Solar system [30] that yielded the coordinates of the bodies at a desired time, calculated by the formulas for elliptic motion in which the ellipse parameters at each time, were determined from the data initially calculated by the Galactica program. In the course of the solution process of the problem, the mathematical model for the Solar system was subjected to a thorough check. Nonetheless, there still existed a probability that, over large time intervals, the insignificant differences between the results of the mathematical model for the Solar system and the coordinate values obtained by the Galactica problem, could affect the evolution of the rotational parameters ε and ψ. We have developed a new program, glc3rte2.for, for joint solution of the orbital problem and the Earth’s rotation problem. In the latter program, in one step over time, the Galactica program solved the orbital problem and, then, the Dormand-Prince method was used to solve, in the same step, the Earth’s rotation problem. With the help of the new program, we have obtained solutions of those problems over different time intervals, including the interval of 200 thousand years. All previously obtained results have found confirmation. The latter check also confirms the mathematical model for the Solar system over large time intervals.

Table 1 gives a quantitative comparison of precession periods P_prN, and the minimum and maximum obliquities (ε_min and ε_max, respectively), accurate to five significant digits. Within this accuracy level, the first two methods have yielded identical results. As it is evident from Table 1, the data calculated by the

Table 1. Comparison of results obtained by three methods for integration of rotational-motion equations over a period of 200 thousand years: RK-4—Runge-Kutta method of the fourth order; DP-8— Runge-Kutta method of the eighth order in Dormand-Prince realization; Gal—the bodies’ coordinates are determined by the Galactica program, and the rotational-motion equations are solved by the DP-8 method.

third method, differed in terms of precession period P_prN in the fourth digit, and in terms of ε in the fifth digit. Since the latter method was a most accurate one, the latter values provided refined results in comparison with the data obtained by the first two methods.

Thus, the various tests and verification of the initial solution method of the Earth’s rotation problem and also, independent solutions of the same problem by three other methods, have confirmed the fact that the Earth’s rotational axis executes oscillations at an amplitude of 7 - 8 times greater than that obtained in the previous solutions.

4.4. Adequacy of the Solution of the Third and Forth Problems

As it was noted previously, the third problem on the determination of the obliquity ε between the moving equatorial plane 4 and the moving orbital plane 5, and the perihelion angle φ_pγ (Figure 2), was a geometrically complex problem because a multitude of revolutions on the axes $\stackrel{\to }{N}$ and $\stackrel{\to }{S}$ , and the perihelion B, were executing about different axes in different directions. For instance, the Earth axis, $\stackrel{\to }{N}$ executed 777 revolutions per 20 million years of problem solution. In the problem, it was required to determine, from the obtained solution of the orbital problem for i,φ_Ω, and φ_р (Figure 2) and the Earth’s rotation problem for θ and ψ, the relative-motion characteristics, ε and φ_pγ.

Those transformations were derived; yet, they involved inverse trigonometric functions, known to be many-valued. The expressions themselves are cumbersome, with some their parts presenting imaginary or infinitely growing expressions. Those singularities needed to be identified for the revealing factors underlying their behavior and for developing algorithms for elimination of these factors. Initially, this problem was solved by means of spherical geometry. However, because of the complexity of involved logical concepts, there was no firm belief in adequacy of the obtained solution. Fortunately, an idea suggesting a second method had emerged. The axial vectors, $\stackrel{\to }{N}$ and $\stackrel{\to }{S}$ were replaced with their projections onto the axes of a Cartesian coordinate system. Then, trigonometric means were employed to derive the required transformations. The two transformation systems have allowed us to reveal errors in each of the systems and fix them, unless both systems started yielding the same results in the examined 20-million-year time interval.

The new algorithm developed for calculating the insolation in the fourth problem was checked [19] by performing insolation calculations using the orbital data by Milankovitch [1], J. Laskar et al. [6]. The new algorithm has yielded results, the same as those yielded by the Milankovitch algorithm.

4.5. Physical Cause of the Difference between the New and Previous Theory

As it was noted previously, according to the new solutions for the angles θ and ψ, the Earth’s rotation axis $\stackrel{\to }{N}$ rotates around the vector ${\stackrel{\to }{M}}_2$ (Figure 2), and it additionally exerts oscillations about that vector. We have calculated the angle θ_M2 between the vector ${\stackrel{\to }{M}}_2$ and the axis $\stackrel{\to }{N}$ and also, the angle ψ_M2 of the equator-plane precession, with respect to the plane normal to ${\stackrel{\to }{M}}_2$ . The main oscillation period of the angles θ_M2 and ψ_M2 proved to be equal to 41.1 thousand years. As was noted previously, the Earth’s orbital axis $\stackrel{\to }{S}$ precesses in a clockwise direction around the vector $\stackrel{\to }{M}$ at a period P_prS = −68.7 kyr, and the Earth’s axis $\stackrel{\to }{N}$ rotates around the vector ${\stackrel{\to }{M}}_2$ as well, at a period P_prN = −25.74 kyr. Since both axes precess in one direction, the effect of the orbital motion on the rotational motion will depend on the relative precession velocity. The period of that action P_rl, is defined by the reciprocal ratio of angular velocities

$P_{rl}=\frac{P_{prS}\cdot P_{prN}}{P_{prS}-P_{prN}}$ . After substitution of actual values, we have obtained a value P_rl = −41.1 kyr, that is, the

very oscillation period of the angles θ_M2 and ψ_M2. In the previous Astronomical theory of climate change, the main oscillation period of the obliquity, ε was equal to 41 thousand years. As was noted previously, that theory was based on a simplified solution of the Earth’s rotation problem. The obtained simplified solution has resulted in identical positions of the vectors ${\stackrel{\to }{M}}_2$ and $\stackrel{\to }{M}$ (Figure 2) and finally, in the obtained period of 41 thousand years.

In the new theory (Figure 2), the vectors ${\stackrel{\to }{M}}_2$ and $\stackrel{\to }{M}$ of the precessional axes have different orientations. That is why the moments of forces by which all bodies act on the Earth exhibit a wide range of oscillations. As a result, the oscillation range of θ with respect to the motionless orbit, (3) also increases. In addition, the oscillation range of the angle, ε between the moving orbital axis $\stackrel{\to }{S}$ and the moving axis $\stackrel{\to }{N}$ increases as well. All in all, in the new theory, the oscillation amplitude of angle, ε turns out to be exceeding the same amplitude in the previous theory by a factor of 7 - 8.

Thus, the oscillation period of obliquity ε in the previous theory of Earth-axis dynamics, was the period P_rl of the relative precession of Earth’s rotational axis $\stackrel{\to }{N}$ and its orbit $\stackrel{\to }{S}$ . That is why it is a fact that the precessional axes $\stackrel{\to }{M}$ and ${\stackrel{\to }{M}}_2$ of the orbital axis $\stackrel{\to }{S}$ and the Earth axis $\stackrel{\to }{N}$ were assumed coincident, was the physical cause behind obtaining erroneous results in the previous theory.

4.6. Final Verification

The final check of the Astronomical theory of climate change consisted in the comparison of its results with paleoclimate. While analyzing data gained by geologists, geographers, and other specialists in the field of paleoclimate, we have found [22 , 31 , 32] four extremes of summer insolation $Q_s^{\text{65}N}$ having occurred over the past 50 thousand years; namely, 4.16, 15.88, 31.28, and 46.44 thousand years ago; those dates correspond to the middle of Holocene, to the middle of the last ice age, to the middle of a warm period, and to the middle of the penultimate ice maximum, respectively. Those events are called differently in different regions of the world; yet, all of them have left their traces in Siberia, Europe, and Northern America.

The whole complex of performed studies and their tests outlined here give us grounds to assert that, in the present article, results of Astronomical theory of climate change obtained by taking into account all studies performed during past centuries, are reported.

Changes of obliquity, ε and angle of the perihelion, φ_pγ, as well as the eccentricity of the orbit, e are available at http://www.ikz.ru/~smulski//Data/Insol/ in the file OrAl1c_8.prn for ±100 years, in the file OrAl-5kyr.prn for 5 thousand years ago (ka), in the file OrAl-200ky.prn for 200 ka, in the file OrAl0-5My.prn for 5 Ma. The program, Ins12bdEn.mcd for analyzing these data, calculating the insolation of the Earth and the building of graphs is available on this site.

5. CONCLUSION AND FURTHER DEVELOPMENT OF THE ASTRONOMICAL THEORY OF CLIMATE CHANGE

As a result of the interaction of Solar-system bodies, evolution of the Earth’s orbital and rotational motions proceeds; this evolution in turn gives rise to insolation oscillations being the cause of climate changes observed over time intervals of tens of thousands of years. The same interactions lead to the evolution of the Sun’s motion about the center of mass of the Solar system [20 , 33] and also, to a change in the rotational motion of the Sun. Studies show [34], that the change in motions presents the cause of variations in solar activity. The radiation fluxes due to the Sun and its substance act on the Earth’s upper shells, thus leading to circulation changes arising in its atmosphere and ocean. Very likely, these factors act as the cause of short-term climate variations occurring over periods of several ten and hundred years. Further development of the Astronomical theory of climate change will be related to the determination of those oscillations.

ACKNOWLEDGEMENTS

The materials of this work were obtained as a result of research at the Institute of Earth’s Cryosphere, Tyum. SC of SB RAS, Federal Research Center for two decades, and in recent years, the research has been carried out under the theme 121041600047-2.

The results of the new Astronomical theory of climate change are based on the solution of the problems about the Earth’s orbital and rotational motion, obtained from supercomputers in the shared use center at the Siberian Supercomputer Center, Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk, Russia.

English language within the paper has been edited by Walter Babin, researcher and founder of the General Science Journal, Canada.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Milankovitch, M. (1930) Mathematische Klimalhre und Astronomische Theorie der Klimaschwankungen (Mathematical Climatology and the Astronomical Theory of Climate Change). Gebruder Borntraeger, Berlin.
[2]	Brouwer, D. and Van Woerkom, A.J.J. (1950) The Secular Variation of the Orbital Elements of the Principal Planets. Astronomical Papers for American Ephemeris and Nautical Almanac, 13, 81-107.
[3]	Sharaf, Sh.G. and Budnikova, N.A. (1969) The Secular Variation of the Orbital Elements of the Earth and the Astronomic Theory of Climate Change. Transactions, No. 14, 48-109. (In Russian)
[4]	Berger, A. and Loutre, M.F. (1991) Insolation Values for the Climate of the Last 10 Million Years. Quaternary Science Reviews, 1, 297-317. https://doi.org/10.1016/0277-3791(91)90033-Q
[5]	Edvardsson, S., Karlsson, K.G. and Engholm, M. (2002) Accurate Spin Axes and Solar System Dynamics: Climatic Variations for the Earth and Mars. Astronomy & Astrophysics, 384, 689-701. https://doi.org/10.1051/0004-6361:20020029
[6]	Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A.C.M. and Levrard, B. (2004) A Long-Term Numerical Solution for the Insolation Quantities of the Earth. Astronomy & Astrophysics, 428, 261-285. https://doi.org/10.1051/0004-6361:20041335
[7]	Smulsky, J.J. (2014) Basic Positions and New Results of the Astronomical Theory of Climate Change. VINITI, Tyumen, No. 258-V2014, 30 p. (In Russian) http://www.ikz.ru/~smulski/Papers/OsPoATLP3.pdf
[8]	Smulsky, J.J. (2004) The Theory of Interaction. Cultural Information Bank, Ekaterinburg, 304 p. http://www.ikz.ru/~smulski/TVEnA5_2.pdf
[9]	Smulsky, J.J. (2012) Galactica Software for Solving Gravitational Interaction Problems. Applied Physics Research, 4, 110-123. https://doi.org/10.5539/apr.v4n2p110
[10]	Smulsky, J.J. (2012) The System of Free Access Galactica to Compute Interactions of N-Bodies. International Journal of Modern Education and Computer Science, 11, 1-20. https://doi.org/10.5815/ijmecs.2012.11.01
[11]	Smulsky, J.J. (2018) Future Space Problems and Their Solutions. Nova Science Publishers, New York, 269 p. http://www.ikz.ru/~smulski/Papers/InfFSPS.pdf
[12]	Smulsky, J.J. (2011) The Influence of the Planets, Sun and Moon on the Evolution of the Earth’s Axis. International Journal of Astronomy and Astrophysics, 1, 117-134. https://doi.org/10.4236/ijaa.2011.13017
[13]	Smulsky, J.J. (2016) Fundamental Principles and Results of a New Astronomic Theory of Climate Change. Advances in Astrophysics, 1, 1-21. https://doi.org/10.22606/adap.2016.11001
[14]	Smulsky, J.J. (2020) The Evolution of the Earth’s Rotational Movement for Millions Years. The Complex Systems, 1, 3-42. http://www.ikz.ru/~smulski/Papers/EVDZ03_1EnJc.pdf
[15]	Smulsky, J.J. and Krotov, O.I. (2014) New Computing Algorithm of the Earth’s Insolation. Applied Physics Research, 6, 56-82. https://doi.org/10.5539/apr.v6n4p56
[16]	Newcomb, S. (1895) The Elements of the Four Inner Planets and the Fundamental Constants of Astronomy. Government Printing Office, Washington DC, 202 p. https://doi.org/10.1086/102176
[17]	Simon, J.L., Bretagnon, P., Chapront, J., Chapront-Touze, M., Francou, G. and Laskar, J. (1994) Numerical Expression for Precession Formulae and Mean Elements for the Moon and the Planets. Astronomy & Astrophysics, 282, 663-683.
[18]	Smulsky, J.J. (2018) The Sun’s Movement in the Sky Now and in the Past. Open Access Library Journal, 5, e4250. https://doi.org/10.4236/oalib.1104250
[19]	Smulsky, J.J. (2018) New Astronomical Theory of Ice Ages. LAP Lambert Academic Publishing, Riga, 132 p. (In Russian) http://www.ikz.ru/~smulski/Papers/InfNwATLPEn.pdf
[20]	Melnikov, V.P. and Smulsky, J.J. (2009) Astronomical Theory of Ice Ages: New Approximations. Solutions and Challenges. Academic Publishing House, Novosibirsk. http://www.ikz.ru/~smulski/Papers/AsThAnE.pdf
[21]	Smulsky, J.J. (2016) Evolution of the Earth’s Axis and Paleoclimate for 200,000 Years. LAP Lambert Academic Publishing, Saarbrucken, 228 p. (In Russian) http://www.ikz.ru/~smulski/Papers/InfEvEAPC02MEn.pdf
[22]	Smulsky, J.J. (2016) New Results on the Earth Insolation and Their Correlation with the Late Pleistocene Paleoclimate of West Siberia. Russian Geology and Geophysics, 57, 1099-1110. https://doi.org/10.1016/j.rgg.2016.06.009
[23]	Grebenikov, E.A. and Smulsky, J.J. (2007) Evolution of the Mars Orbit on Time Span in Hundred Millions Years. Reports on Applied Mathematics. Russian Academy of Sciences: A. A. Dorodnitsyn Computing Center, Moscow, 63 p. (In Russian) http://www.ikz.ru/~smulski/Papers/EvMa100m4t2.pdf
[24]	Laskar, J. (1996) Marginal Stability and Chaos in the Solar System. In: Ferraz Mello, S., et al., Eds., Dynamics, Ephemerides and Astrometry of the Solar System, IAU, Amsterdam, 75-88. https://doi.org/10.1017/S0074180900127160
[25]	Bretagnon, P., Rocher, P. and Simon, J.L. (1997) Theory of the Rotation of the Rigid Earth. Astronomy and Astrophysics, 319, 305-317.
[26]	Smul’skii, I.I. (2013) Analyzing the Lessons of the Development of the Orbital Theory of the Paleoclimate. Herald of the Russian Academy of Sciences, 83, 46-54. https://doi.org/10.1134/S1019331613010073
[27]	Mel’nikov, V.P., Smul’skii, I.I. and Smul’skii, Ya.I. (2008) Compound Modeling of Earth Rotation and Possible Implications for Interaction of Continents. Russian Geology and Geophysics, 49, 851-858. https://doi.org/10.1016/j.rgg.2008.04.003
[28]	Krutko, P.D., Maksimov, A.I. and Skvortsov, L.M. (1988) Automated System Engineering Algorithms and Programs. Radio and Communication, Moscow, 298 p. (In Russian).
[29]	Hairer, E., Nersett, S. and Vanner, G. (1990) Solving Ordinary Differential Equations. Mir, Moscow, 513 p. (In Russian) https://doi.org/10.1007/978-3-662-09947-6
[30]	Smulsky, J.J. (2007) A Mathematical Model of the Solar System. In: Theoretical and Applied Problems in Nonlinear Analysis, Russian Academy of Sciences: A.A. Dorodnitsyn Computing Center, Moscow, 119-138. (In Russian) http://www.ikz.ru/~smulski/Papers/MatMdSS5.pdf
[31]	Smulsky, J.J. and Ivanova, A.A. (2018) Experience of Paleoclimate Reconstruction on Insolation Change on Example of Western Siberia in the Late Pleistocene. Climate and Nature, 1, 3-21. (In Russian) http://www.ikz.ru/~smulski/Papers/OpRcnPClmt6J.pdf
[32]	Smulsky, J.J. (2019) The Upcoming Tasks of Fundamental Science. Sputnik + Publishing House, Moscow, 134 p. (In Russian) http://www.ikz.ru/~smulski/Papers/InfPrZaFN.pdf
[33]	Smulsky, J.J. (2017) Cosmic Impacts on the Earth and Their Influence on the Arctic. The Complex Systems, 4, 27-42. (In Russian) http://www.ikz.ru/~smulski/Papers/CsmAcEIA.pdf
[34]	Mörner, N.-A. (2016) Planetary Influence on the Sun and the Earth, and a Modern Book-Burning. Nova Science Publishers, New York, 196 p.