Milankovitch, the father of paleoclimate modelling

The origin of the long-term variations of the astronomical elements used by Milankovitch are first described, followed by the value of the astronomical periods. The detailed calculations by Milankovitch of the incoming solar radiation during the astronomical and caloric half-years are summarized, stressing the originality of the caloric ones. The second original contribution of Milankovitch to paleoclimate research was without any doubt his mathematical climate. How this model 10 allowed him to give the caloric summer and winter insolation a climatological meaning is illustrated.


Introduction
Paleoclimatology is primarily a reconstruction of past climatic variations on the basis of proxy records. It aims also to explain these variations from principles of climatic behaviour using climate models. Milankovitch has contributed significantly to this second objective by using the astronomical parameters 15 to compute the long-term variations of his caloric insolation which he used in a climate model (although very simple) to reconstruct the past climates. This paper intends to underline the fundamental and original contributions of Milankovitch in the understanding of the long-term climatic variations over the last one million years. 20 The two remarkable books of Milutin Milankovitch, his 1920 Théorie Mathématique written in French and -his 1941 Kanon der Erdbestrahlung written in German and translated in English in 1969, have largely contributed to his reputation. The celebration of the 100th anniversary of his 1920 French book is a good opportunity to stress what were his main contributions and to "rendre à César ce qui appartient à César (give back to Caesar what belongs to Caesar)". There is indeed a tendency to over-inflate one's 25 work for reasons that have nothing to do with the scientist, but for reasons that have to do with corporative or national politics of history of science. It is important to stress here that Milankovitch was always very careful through all his publications referring properly to the publications of others when he was using their results.
The first to calculate these astronomical elements was the French astronomer Joseph-Louis Lagrange (1736Lagrange ( -1813. He made the calculation for the six great planets (1781)(1782). At that time, Uranus was 35 only going to be discovered and the mass of the planets could only be roughly estimated. Aware of this uncertainty on the masses, Lagrange investigated its possible influence on his calculations, a formulation that was going to be used more than one century later by Vojislav Miskovitch (1892-1976, a colleague whom Milankovitch solicited his collaboration. 40 During the early 19th century, Pierre-Simon Laplace (1749-1827) wrote his 5-volume Celestial Mechanics between 1799 and 1825. Philippe Gustave le Doulcet, Comte de Pontécoulant (1795de Pontécoulant ( -1874 carried out the computation of the long-term variations of the elements of the great planets but with a few decimals only. 45 It is during the second part of the nineteen century that Urbain Le Verrier (1811-1877) came with a new theory of the planetary motion (1855) and the calculation of the secular perturbations (1856). He published the numerical values of eccentricity (with a precision of 10-4), longitude of the perihelion (in arc minutes), inclination (in arc seconds) and longitude of the node over 100,000 years before and after 1800 A.D. each tenth millennium. His calculations were carried out before he discovered the planet 50 Neptune. As this planet could not therefore be included in the Le Verrier calculations, John Nelson Stockwell (1832-1920 computed the secular perturbations by considering all the eight planets known at his time (1873): Mercury, Venus, the Earth, Mars, Jupiter, Saturn, Uranus and Neptune. Stockwell, as Le Verrier, had the possibility of correcting its computations by using better values of masses. This work completed the calculations of the secular perturbations of the great planets, but an error was discovered 55 by Harzer (1895) twenty-two years later. Finally it is the German Ludwig Pilgrim (1879Pilgrim ( -1935 (known mostly as a pioneer in colorimetry) who computed, at a time when Milankovitch was completing his doctoral degree thesis in 1904, the astronomical elements requested for the computation of insolation. Pilgrim extended (Pilgrim, 1904) the numerical computations of eccentricity, obliquity and longitude of the moving perihelion, using the Stockwell integrals for every fifth millennium over 1,010,000 before 60 and 40,000 years after 1850 A.D. and also for dates where the longitude of the perihelion was either 90° or 270° (Northern Hemisphere summer at perihelion or at aphelion). In Milankovitch's own words, Pilgrim was the first to compute adequately the elements affecting the long-term variations of insolation, so that Milankovitch could use them for his research. Pilgrim tried also to treat the Ice Age mathematically, but according to Milankovitch, this treatment of the climatological part by Pilgrim was a failure.
Milankovitch clearly indicates that for the calculation of his incoming solar radiation (insolation in short), he used first Stockwell-Pilgrim's values of the eccentricity, obliquity and precession for the last 1,000,000 years before 1850 A.D.. The numerical values of these astronomical parameters were reproduced in his 70 1920 book for the last 500,000 years. The insolation values for 55°, 60° and 65°N were published in Köppen-Wegener (1924) before the calculation was extended to other geographical latitudes and the new values published in his Mathematische Klimalehre in 1930.
Because of some errors in Stockwell (already detected by Harzer) and because he wanted to use the 75 astronomical parameters based on the most reliable values of the planetary masses, Milankovitch decided to use the Le Verrier calculation including his corrections for the masses. For completing this work, he asked the collaboration of his colleague, Prof V. Miskovitch (1892Miskovitch ( -1976. Miskovitch made the necessary corrections of the masses following Le Verrier procedure and computed the long-term variations of eccentricity e, obliquity ε, and climatic precession esinᴨᵞ , (ᴨᵞ the longitude of the perihelion ) for the past 80 600,000 years before 1800 A.D. with the following initial values: e0=0.0168 ᴨᵞ0=99°30´ ε0=23°27´55" All these values were published in Mathematische Klimalehre, 1930 andin Astronomische Mittel in 1938 85 and reproduced in Table IX  The comparison of the insolation values that Milankovitch calculated from these solutions show a good agreement. Milankovitch concluded that a further improvement of the planetary masses using the formulation by Le Verrier would not change "the essential features of the secular course of insolation as I have calculated". 95 Milankovitch was however well aware that the solution by Le Verrier could not be extended over millions of years because of the limited accuracy of the perturbation calculation that was based upon classical mechanics, missing the Einstein relativistic displacement of the perihelion of the planets. 100 Milankovitch calculated with great details the incoming solar radiation on the Earth, but seemed to have not been much interested in the astronomical periodicities themselves. From his table (based on Stockwell-Pilgrim), he simply deduced the average period of the oscillations of eccentricity as being 92 kyr, varying between 77 and 103 kyr. For precession, he found an average period of 21 kyr, varying between 16.3 kyr and 25.8 kyr. For the longitude of the moving perihelion, he explains that its 105 irregularities are due to the longitude of the fixed perihelion, but whether the perihelion has a mean motion remains an open question. For obliquity, he noticed that it "oscillates between extremely narrow limits" with a relatively stable period of 40 kyr varying between 38 and 45 kyr. It is actually the French mathematician Joseph Alphonse Adhémar (1797-1862) who was the first in 1842 to deduce the value of 21,000 years for precession by combining the astronomical precession calculated from the value of 50.1" The harmonics of precession, in particular those with a period of 19,000 and of 23,000 yearsalso found 115 by Hays et al. (1976) in their geological data -, and the 400,000-year of eccentricity were discovered by Berger (1973Berger ( , 1978a who calculated all the periods in the expansion of the long-term variations of the astronomical variables used in the calculation of insolation.

Periods of astronomical parameters
It is worth to note that Milankovitch was very much interested by obliquity. One of the reasons is probably 120 the strong obliquity signal in his caloric half-year insolation (see Figure 1). Milankovitch mentioned many times that the authors who were mainly stressing precession, were not sufficiently or not properly taking into account obliquity. Tables XII and XIII of his Canon give the change of the radiation of his Table XI for an increase of the obliquity by one degree (respectively in canonic units and in percent). These three tables were first published in his "The Problem of the Astronomical theories of the Ice Ages" in 1914 125 showing in details for the first time the influence of obliquity upon insolation. To put it briefly, an increase of obliquity reduces the latitudinal contrast between the equator and the poles mainly in the annual irradiation and increases the seasonal one (with an augmentation of the summer irradiation and a reduction of the winter one), with a similar effect in both hemispheres.  . It was actually Murphy who was the first in 1869 to put forward the idea that a long, cool summer and a short, mild winter are the most favourable conditions for glaciations (a hypothesis totally opposite to Croll, 1864). This idea was taken up by the Austrian climatologist Rudolf Spitaler (1859Spitaler ( -1946 half a century later (in 1921). It does therefore follow that it is not Milankovitch who originated 140 this principle as some authors have claimed and still claim calling it the "Milankovitch model". Milankovitch has actually popularized and spread the idea under the advice of Köppen (1941) who claimed : "the diminution of heat during the summer half-year is the decisive factor in glaciation" and also following the comments made earlier by Penck and Brückner (1909) and Bruckner et al. 1925: "From the climatological point of view, glaciers are not favoured by severe winter…but by a mild winter and a 145 cool summer".

Caloric Seasons
When dealing with the astronomical seasons, the long-term variations of both their total irradiation and of their length must be taken into account. To avoid this difficulty, Milankovitch introduced the caloric seasons. This concept of caloric seasons is one of the two most important and original contributions of 150 Milankovitch. These divide the year into two equally long seasons, one of whichthe caloric summercomprises all days during which the irradiation at the given latitude is stronger than on any day of the other half-yearthe caloric winter. Because the semi-major axis of the Earth's orbit, the sidereal period of revolution of the Earth around the Sun and, to an excellent approximation, the tropical year do not change with time, the length of these caloric seasons are exactly 182.6211 mean solar days when the 155 tropical year is used. This, however, does not solve the problem completely because the start and end of these half-year seasons change with time and because the double maximum and minimum characterizing the insolation in the intertropical regions.
Milankovitch noticed that he discovered these caloric seasons after his 1920 book on "Théorie In his Canon (1941), Milankovitch devoted twenty pages to the "quantities of heat received by a latitude during a caloric summer and winter half-year". From the formulas that he developed, it is clear that, 165 during their local season, the impact of the variations of obliquity is the same in both hemispheres and maximum in the high latitudes, whereas the impact of climatic precession is opposite in the two hemispheres and maximum at the low latitudes.
In Chap XX of his Canon, Milankovitch gave the numerical values of the caloric Northern and Southern 170 hemisphere summer half-years for 1800 (Table XXIII) and over the last 600,000 years (Table XXV) in canonic units (the canonic units introduced by Milankovitch are the units obtained if the solar constant is the unit of solar radiation and if the unit of time is 100,000 instead of seconds). Since no hypotheses were introduced for these calculations, Milankovitch, who was convinced of the perennity of his work, decided to call his results "Kanon der Erdbestrahlung" (Canon of Insolation). 175 With Köppen approval, Milankovitch preferred not to reproduce anymore the numerical values of insolation themselves, but rather to transform them in fictitious latitudes, called the 65°N equivalent latitudes; this was first published in his Mathematische Klimalehre. These latitudes are actually the present-day latitudes which received during the Northern Hemisphere caloric summer half-year the same 180 irradiation than 65°N in the past. A fictitious motion of these latitudes to the south corresponds therefore to an increase of the summer irradiation.

The mathematical climate
It must be stressed that the main contributions of Milankovitch were not only based on his insolation and radiation curves, but also on his mathematical computation of the thermal effects of the secular march of 185 insolation, his so-called mathematical climate. The direct effect calculated if insolation only was varying was published in his 1930 Mathematische Klimalehere. From the Sefan law and a grey body model with the reflective power of the surface and the absorption coefficient in IR kept invariable, Milankovitch calculated the long-term variations of the mean temperatures of the caloric summer and winter half-years (ΔT=ΔQ/150, ΔQ in canonic units).
If the ice cover and other feedbacks are taken into account, the indirect effects could be estimated. This was published in 1938 Astronomische Mittel where Milankovitch calculated first the altitude of the snow limit, Hi, as a function of the caloric summer insolation. This calculation was based on the correlation between these variables according to Köppen snow limit data for different latitudes. According to his 195 relationship, any variation of the summer irradiation by 1 canonic unit produces a shift of the snow limit altitude by 1m (ΔHi=1.09ΔQS).
These relationships allowed Milankovitch to give the caloric summer and winter insolation a climatological meaning. This shows how much he was concerned by climate and its variations. As these 200 relationships are simple and straightforward, Milankovitch did not published any additional tables and referred only to his tables providing the long-term variations of the caloric summer and winter half-year insolation. From his Table XXV, we can see that the deficit in summer radiation can reach 573 canonic units at 75°N, 22100 years ago, which according to his formula means a drop of the altitude of the snow line by more than 500m. Following the Köppen table of the altitude of the snow limit for different 205 latitudes, the polar cap can then have extended from 75°N up to 65°N, it means covering an area 2.75 times greater. This kind of deficit can also be reached in the tropical latitudes with an accompanying lowering of the snow limit altitude which, as noted by Milankovitch "refutes the opinion expressed by some geologists that insolation cannot explain such displacement". Milankovitch also pointed out that "owing to such variations of the summer irradiation, the mean summer temperature dropped from time to 210 time by more than 5° in the high and temperate latitudes of both hemispheres and even in the tropical latitudes".

Irradiation over the polar caps
Most important is that such an increase of the size of the polar snow cover changes the reflective power of the Earth. This is why to complete his Kanon, Milankovitch decided to compute the long-term 215 variations of the mean summer and winter insolation per unit surface area of the Northern and Southern polar snow-caps over the last 600,000 years (Table XXVIII). The extent of these polar caps was deduced from the treatise by Wundt (1933). In this treatise, the Northern cap extends presently to 75°N and reached 55°N at the maximum of the Ice Age. From these values, Milankovitch could compute the long-term variations of the insolation over such polar caps delimited by the parallel 55° assuming that the extension 220 of the snow cap was always proportional to the corresponding deficit in summer radiation. If the albedo is kept constant, it can be seen that the minimum summer radiation over the Northern cap reaching 55°N occurred 230,000 years ago, with a radiation deficit Δ1QS=363Δε-7900Δ(esinᴨᵞ), compared to the present, amounting 660 canonic units (ε is the obliquity, e the eccentricity and ᴨᵞ the longitude of the perihelion). As this deficit caused a southward extension of the cap of 20°, Milankovitch concluded that a change by 225 one canonic unit corresponds to a meridional change in the extent of the Northern polar cap by about 1´82 or 3.37 km. This means also that the Northern snow cap totally disappears for an increase of the summer insolation by 495 canonic units. This occurred quite a few times over the last 600,000 years, as for example 10,000 and 127,000 years ago.
Taking into account the reflective power of snow at the Earth' surface, for the cap reaching 55°N, Δ2QS=0.802(13520+Δ1QS)[sin(75°+1´82Δ1QS)-sin75°] shows that large negative amplitudes occurred several times over the last 600,000 years. At 230 kyr BP, the total deficit ΔQS= Δ1QS+Δ2QS amounts now 2180 canonic units. This is far more than the 660 deficit calculated if the reflective power of snow is not taken into account. It corresponds to a displacement downward of the snow limit of 2180 m which is 235 about the present altitude of the snow limit at 55°N. This implies that the polar cap must have reached this latitude at that time, what according to Milankovitch was actually observed in the geological reconstruction. It is also interesting to notice with Milankovitch that the deficit of the annual radiation at 230 kyr BP amounted 1920 caloric units which means a decrease in the annual temperature by 6.4 °C ( Using the data for the cap reaching 45°N, nine large deficits can be observed at 590.3, 550, 475.6, 435, 230, 187.5, 115, 71.9 and 25 kyr BP. These can be assembled in groups corresponding to the four glacial periods of the Penck-Bruckner scheme recognized by Köppen in the Milankovitch 65°N equivalent 245 latitude.
These new results taking into account the reflective power of the polar caps in addition to the long-term variations of insolation were published in "Neue Ergebnissen" (Milankovitch, 1937(Milankovitch, -1938 and, according to Milankovitch are "absolutely sufficient to explain the full extent of even the greatest climatic events of 250 the Quaternary and to clearly show their causes". All these calculations show clearly that Milankovitch can be named the "father of paleoclimate modelling", certainly more specifically than the father of the astronomical theory in general (the first to propose the variations of the Earth's orbit as the causes of climate changes was Jens Esmark (1763-1839), Thanks are due to Prof Vladimir Jankovic from the Center for the History of Science, Technology and Medicine at University of Manchester for his thoughtful comments and his proof-reading suggestions. Thanks also to Dr. Qiuzhen Yin from Université catholique de Louvain for reading the paper and for providing Figure 1.