In defense of Milankovitch

1. Introduction

[2] Milutin Milankovitch was among the first to highlight the role that periodic variations in the Earth's orbit might play in climate [e.g., Milankovitch, 1941]. Working from the orbital calculations of Joseph Adhémar [Adhémar, 1842], he computed time series of insolation as a function of latitude and season, and also undertook basic energy balance studies to estimate the temperature changes that might result. In collaboration with Alfred Wegner and Wladimer Köppen, [e.g., Köppen and Wegener, 1924], he further argued that periods of minima in summertime insolation in the northern high latitudes coincided with the half dozen then-known instances of glacial advances in Europe. The landmark paper of Hays et al. [1976] built upon this earlier work, reconstructing global ice volume using oxygen isotopes in two deep-sea sediment cores. Hays et al. showed that global ice volume in the late Pleistocene also reflected these orbital variations. Since this early work, there has been a profusion of data analysis and modeling aimed at fleshing out the links between climate and insolation variations [e.g., Roe and Allen, 1999; Paillard, 2001]. Much of this research has been viewed as in some way evaluating or verifying a ‘Milankovitch hypothesis’ (or Milankovitch theory) of climate. However, several challenges have recently emerged, with studies (1) arguing that variations of atmospheric CO₂ and tropical sea surface temperatures played an important role in ice-age cycles [e.g., Shackleton, 2000; Lea, 2004], (2) questioning the timing of temperature changes relative to insolation and the global extent of the ice age climates [Winograd et al., 1992; Gillespie and Molnar, 1995], and (3) highlighting the relatively small fraction of total variance in global ice volume associated with insolation variations [Wunsch, 2004]. Moreover, progress has been impeded by the lack of a well-formulated, specific, and generally-accepted hypothesis. The term ‘Milankovitch hypothesis’ is used in a variety of ways, ranging from the simple expectation that one ought to see orbital frequencies in time series of paleoclimate proxies, to the implication that all climate variability with time scales between 10³ and 10⁶ yr is fundamentally driven by orbital variations. Somewhere in the middle of this are the more vague statements found in some form in many textbooks, that orbital variations are the cause, or pacemaker, of the Pleistocene ice ages. Phrases like Milankovitch curves, Milankovitch insolation, Milankovitch frequencies, Milankovitch forcing, and Milankovitch cycles pervade the literature, adding to the somewhat nebulous picture.

[3] Three main orbital parameters induce insolation variations: eccentricity (∼100 and ∼400 kyr periods), obliquity (41 kyr), and climatic precession (∼19 and ∼23 kyr) [e.g., Imbrie and Imbrie, 1980]. Myriad paleoclimate proxies around the globe have been analyzed to evaluate climate variability at these timescales. Theoretical understanding of climate dynamics (and of many of the proxies themselves) is not yet at the stage where the physical causes of the variations the proxies record can be known in any detail. It is constructive therefore to try to separate the general question of how orbital variations are expressed in the climate of different regions and in different paleoclimate proxies, from the specific question of the causes of changes in the extent or volume of ice sheets. While Milankovitch's contributions are obviously relevant to both, any specific hypothesis or theory named after him should most appropriately be about the latter, since it is much closer to the original compass of his work. In this paper, a specific formulation of the Milankovitch hypothesis is suggested and defended: orbitally-induced variations in summertime insolation in the northern high latitudes are in antiphase with the time rate of change of ice sheet volume.

2. Results

[4] Reconstructions of global ice volume during the Pleistocene rely on the measurement of oxygen isotopes incorporated into the shells of foraminifera, records of which are recovered from deep-sea sediment cores (and which also in part reflect deep ocean temperatures [e.g., Shackleton, 2000]). We compare two such records in Figure 1. It has been estimated that at the last glacial maximum (∼21 kyrs ago) approximately 85% of the extra ice volume compared to the present was in the Northern Hemisphere [Peltier, 2004], so the records may be taken as predominantly reflecting Northern Hemisphere ice volume. This is the interpretation placed on the record in this paper although we note that Raymo et al. [2006] suggest the possibility that Southern Hemisphere ice fluctuations might play a disproportionate role in the oxygen isotope signal.

[5] The SPECMAP record [Imbrie et al., 1984] assumes a priori that ice volume and orbital forcing are related, and tunes the depth profile of the ocean core to yield an age scale. Because of this tuning procedure, it is likely that variability at orbital frequencies is overestimated [e.g., Huybers and Wunsch, 2004]. Huybers and Wunsch have recently developed a record (HW04) that is independent of any such orbital-tuning assumptions. The overall variations and timing of the deglaciations still agree with the SPECMAP record to within the estimated errors. Remaining errors probably displace energy from the orbital frequencies, and thus HW04 likely underestimates the fraction of energy at those frequencies. Therefore these two records can be taken as providing guiding bounds on the strength of the association between insolation variations and global ice volume.

[6] Figure 1 shows a comparison of the two ice volume reconstructions with variations in daily-averaged summer solstice insolation at 65N (hereafter referred to as June 65N) for the last 750 kyr. Peak-to-peak amplitudes are close to 100 Wm⁻². The maximum lag correlations of the SPECMAP and HW04 records with the insolation variations are −0.4 and −0.2 respectively. As has been noted many times elsewhere, the maximum correlation occurs when the ice volume lags the June 65N insolation curve. The lag is 6 kyr in the SPECMAP record and 8 kyr in the HW04 record. This time lag has been variously attributed to the dynamical response time of ice sheets, to the role of ocean circulation in the global transmission of the climate signal, or to the role of CO₂ or tropical sea surface temperatures in driving ice age cycles [e.g., Shackleton, 2000; Lea, 2004; Imbrie and Imbrie, 1980; Imbrie et al., 1992, 1993]. The correlation between insolation and ice volume arises from the shared variance at the precessional and obliquity frequencies [e.g., Hays et al., 1976; Imbrie et al., 1992]. However as evident from Figure 1, there is much greater variability at lower frequencies in the ice volume than in the insolation. We present results in the time domain here. An alternative statistical analysis using cross-spectral estimates supports these results (see auxiliary materials).

[7] While most studies have focused on the connection between insolation and ice volume (V), there is a more direct physical connection between insolation and the rate of change of ice volume (dV/dt) [e.g., Roe and Allen, 1999; Wunsch, 2003]. This distinction is crucial. First, the mass balance of ice sheets is acutely sensitive to summertime temperature: the characteristically convex profile of ice sheets means that the area of ablation at land based-margins varies strongly with summertime temperature. This effect renders the total ablation rate proportional to approximately the third power of the summertime temperature above some reference value [e.g., Pollard et al., 2000; Roe and Lindzen, 2001; Ohmura et al., 1996]. Second, while the convergence of atmospheric and oceanic heat fluxes plays a large role in wintertime climates, summertime climates in continental interiors are much more strongly controlled by the local radiation balance [e.g., Peixoto and Oort, 1992]. Thus there are strong physical grounds, supported by model studies [e.g., Felzer et al., 1995], for expecting a direct response of summertime temperatures, and hence of ice-sheet ablation rates and dV/dt, to local summertime insolation variations. Figure 2 compares June 65N insolation to dV/dt from the SPECMAP and HW04 records. The maximum correlations are −0.8 and −0.4, and occur with no lag and a lag of 1 kyr, respectively. The strong lag-correlations relative to Figure 1, and the absence of a large lag, are striking and demonstrate essentially concurrent variations in dV/dt and summertime insolation in the northern high latitudes. Both reconstructions therefore support the Milankovitch hypothesis as formulated above.

[8] Were the characteristic time scale of ice volume adjustment short compared to the time scale of insolation variations, the ice would come into equilibrium with the forcing, and there would be an antiphased response between ice volume and insolation. Such would be the case for small mountain glaciers and ice caps. However at the scale of the great continental ice sheets, theory, dynamic models, and geological reconstructions all indicate long adjustment times of several tens of thousands of years [e.g., Weertman, 1961; Greve, 1997; Roe and Lindzen, 2001; Clark et al., 1993]. That an antiphase relationship between dV/dt and insolation is observed is very strong evidence of such a long time scale of adjustment. We note also that the nonlinearities of ice flow and mass balance preclude the application of a single adjustment time scale to an ice sheet. Ice sheet height, area, and volume all have different adjustment times [e.g., Roe and Lindzen, 2001], all of which are dependent on the ice sheet size, the climate state, and the magnitude of the forcing.

[9] It is not possible to unequivocally attribute a climate response to an insolation forcing at a particular latitude and season because, to a good approximation, any such forcing can be constructed from a linear combination of climatic precession and obliquity indices [Imbrie and Imbrie, 1980]. In theory then, an infinite set of insolation curves (or their meridional gradients) can be matched with a given climate signal, leaving the physical mechanisms producing the climate response ambiguous. To address this, regression analyses were performed to find the best-fit linear combinations of the obliquity and climatic precession indices for dV/dt from the SPECMAP and HW04 records (see auxiliary material for methods). These best-fit combinations of the orbital parameters were then compared with characteristic insolation curves. Figure 3 shows that, for the SPECMAP record, the best-fit combination is almost identical to the June 65N insolation. In the case of the HW04 record, the best-fit combination closely matches the summer half-year (April to September average) 65N insolation. Summer might arguably be better defined by the radiation half-year [Huybers, 2006]. We have not tried to distinguish here which should be preferred - the nonlinearity of total ablation (i.e., not just local ablation rate) to summertime temperatures makes it likely that peak summertime temperatures ought to be weighted more than average summertime temperatures. The key point for the purpose of this paper is that, while these analyses do not prove a causal link absolutely, both records confirm the a priori hypothesis based on physical arguments of a direct connection between some sensible measure of summertime insolation in the northern high latitudes and the rate of change of ice volume there.

[10] Atmospheric CO₂ has also been suggested as driving changes in global ice volume [e.g., Shackleton, 2000; Lea, 2004]. The concentration of CO₂ varied between about 200 and 280 ppmv over the last several ice age cycles, and caused approximately 2 Wm⁻² variations in surface longwave radiation forcing [e.g., Ramaswamy et al., 2001]. Comparisons of the impacts of shortwave and longwave radiative forcing appropriate over the ice sheets are not straightforward, but taking summer half-year insolation variations in shortwave (Figure 3), and assuming an albedo of 0.5 for melting ice, variations in summertime shortwave forcing exceed the direct CO₂ radiative forcing by about a factor of five. It has also been reported that the ice volume lags behind CO₂, and this has led to the suggestion that CO₂ variations drive ice age cycles [Shackleton, 2000; Lea, 2004; Ruddiman and Raymo, 2003]. However, cross-spectral analyses in Figure 4 (and lag correlations, auxiliary materials) show that, at frequencies where there is significant coherence between the records, atmospheric CO₂ lags, or is at most synchronous with, dV/dt. In other words, variations in melting precede variations in CO₂. Thus, the relatively small amplitude of the CO₂ radiative forcing and the absence of a lead over dV/dt both suggest that CO₂ variations play a relatively weak role in driving changes in global ice volume compared to insolation variations. This certainly does not rule out CO₂ as a primary cause of tropical or other climate variations, or of the apparent synchronization of the ice-age signal between hemispheres [Lea, 2004]. We note though that dV/dt contributes both to the rate of freshwater input into high latitude oceans (changes in which have been argued to cause global climate signals), and also to the rate of sea level rise (with a world-wide destabilizing tendency on ice shelves and coastal ice masses). The Milankovitch hypothesis as formulated here does not explain the large rapid deglaciations that occurred at the end of some of the ice age cycles (i.e., the several large negative excursions in Figure 2): many studies point to the need to invoke internal dynamics of ice sheets as a mechanism for occasional rapid collapses if a threshold size is exceeded [e.g., Imbrie and Imbrie, 1980; Paillard, 2001]. Nor do the results explain the mid-Pleistocene transition between an earlier interval characterized by ∼40 kyr durations of ice ages and a later interval with ∼80 kyr to ∼120 kyr durations: it has been suggested that this transition may have been due to changes in basal conditions under the ice sheet or perhaps chaotic ice sheet dynamics [Clark and Pollard, 1998; Huybers and Wunsch, 2005]. The prevailing view to date has been that ice sheet volume is the most important variable to consider. While this is obviously the case for global sea level, it is ice sheet extent that matters most for albedo, and ice sheet height that matters for atmospheric circulation [e.g., Broccoli and Manabe, 1987; Shinn and Barron, 1989]. Ice sheets are dynamic systems and these properties can vary quite differently from each other. However, the results presented here demonstrate the critical physical importance of focusing on the rate of change of ice volume, as opposed to the ice volume itself. The available evidence supports the essence of the original idea of Köppen, Wegner, and Milankovitch as expressed in their classic papers [Milankovitch, 1941; Köppen and Wegener, 1924], and its consequence: (1) the strong expectation on physical grounds that summertime insolation is the key player in the mass balance of great Northern Hemisphere continental ice sheets of the ice ages; and (2) the rate of change of global ice volume is in antiphase with variations in summertime insolation in the northern high latitudes that, in turn, are due to the changing orbit of the Earth.