A Paleoclimate Model of Ice-Albedo Feedback Forced by Variations in... Orbit Richard McGehee Clarence Lehman
c 2012 Society for Industrial and Applied Mathematics SIAM J. APPLIED DYNAMICAL SYSTEMS Vol. 11, No. 2, pp. 684–707 A Paleoclimate Model of Ice-Albedo Feedback Forced by Variations in Earth’s Orbit∗ Richard McGehee† and Clarence Lehman‡ Abstract. Earth undergoes long-term temperature cycles alternating between glacial and interglacial episodes. It is widely accepted that changes in Earth’s orbit and rotation axis cause variations in solar input which drive the glacial cycles. However, classic papers have clearly established that the response of Earth’s climate system to orbital forcing is not a simple linear phenomenon and must include nonlinear feedback mechanisms. One of these mechanisms is ice-albedo feedback, which can be modeled as a dynamical system. When combined with the cycles in the orbital elements and compared with the climate data, the model confirms that ice-albedo feedback is an important component of Earth’s climate. Key words. dynamical systems, paleoclimate, Milankovitch cycles, ice-albedo feedback AMS subject classifications. 37N99, 86A40 DOI. 10.1137/10079879X 1. Introduction. For millions of years, glaciers have been a permanent feature on Earth’s surface. The climate has cycled through relatively warm periods, such as we are experiencing now with glaciers conﬁned mostly to Antarctica and Greenland, to relatively cold periods with vast ice sheets extending over the continents. The sun provides the energy powering Earth’s weather and climate. The incoming solar radiation, or insolation, varies over several time scales. Changes in Earth’s orbit and rotation axis, called Milankovitch cycles, cause variations in the insolation at time scales comparable to the glacial cycles. These cycles contain three components: (1) the eccentricity of Earth’s elliptical orbit around the sun; (2) the tilt, or obliquity, of Earth’s axis of rotation; and (3) the precession of Earth’s rotation axis. These components will be discussed in detail below. Many studies, dating back to the mid–nineteenth century, have shown that the Milankovitch cycles aﬀect Earth’s climate . In a seminal paper, Hays, Imbrie, and Shackleton  ﬁrmly established that the Milankovitch cycles contribute substantially to the glacial cycles. Since the incoming radiation is a fairly complicated function of both time and space, it is often convenient to reduce it to a single quantity. A quantity traditionally used is the daily average insolation at a particular latitude on a particular day, typically at 65◦ north on the summer solstice [7, 11], a quantity we will denote by Q65 . This quantity changes slightly from ∗ Received by the editors June 15, 2010; accepted for publication (in revised form) by M. Zeeman January 31, 2012; published electronically May 31, 2012. http://www.siam.org/journals/siads/11-2/79879.html † School of Mathematics, University of Minnesota, Minneapolis, MN 55455 ([email protected]). This author was supported in part by NSF grant DMS-0940366. ‡ Department of Ecology, Evolution and Behavior, University of Minnesota, Minneapolis, MN 55455 ([email protected]). This author was supported in part by a Resident Fellowship Grant from the University of Minnesota’s Institute on the Environment. 684 Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 685 year to year, depending on the components of the Milankovitch cycles, and can be computed and analyzed as a time series. As measured by the power spectrum of Q65 , the major contribution to variations in insolation is caused by the precession, with changes in the obliquity contributing less, and changes in the eccentricity contributing the least . Using ocean sediment data from the last 468,000 years, Hays, Imbrie, and Shackleton noted that the major components of climate variation are in the frequencies associated with the eccentricity, followed by frequencies associated with obliquity, with frequencies associated with precession forming the smallest component, the exact opposite of the order found in Q65 . (See Figure 1.) They ascribed this reversal of the contributions to possible nonlinear eﬀects in climate dynamics . precession eccentricity obliquity obliquity obliquity obliquity eccentricity eccentricity eccentricity precession precession precession Q65 Hays climate data Zachos climate data Model Prediction Figure 1. Relative contributions of the Milankovitch components to insolation and climate. In a more recent study, using ocean sediment data from the last 4,500,000 years, Zachos et al.  reached a somewhat diﬀerent conclusion. They found the strongest contribution to climate response in frequencies associated with obliquity, followed by eccentricity, with negligible contribution from precession ; see Figure 1. One should note that the discrepancy between the Hays climate data and the Zachos climate data is not simply a question of better data collected between 1976 and 2001. The data show that the frequencies associated with eccentricity are more pronounced during the last million years. Indeed, an analysis of the modern data over the last million years still yields the result described by  and illustrated in Figure 1. Actual climate dynamics are very complex, involving much more than insolation and certainly much more than insolation distilled down to a single quantity, Q65 . Recent advances in modeling the relationship between Milankovitch cycles and glacial cycles have involved the introduction of triggering mechanisms, whereby changes in insolation caused by the Milankovitch cycles trigger glacial retreats, which then continue due to feedback eﬀects [3, 5, 10]. Other studies have replaced Q65 with the total insolation during the summer [4, 6]. One feedback mechanism is the ice-albedo eﬀect. The reﬂectivity of Earth’s surface, called the albedo, is diﬀerent for ice than for water or land. If more of Earth is covered with ice, then more of the solar energy is reﬂected back into space, leaving less to heat Earth. On the other hand, if less ice covers Earth, more solar energy is absorbed. This eﬀect constitutes positive feedback. As the ice advances, more solar energy is reﬂected, cooling Earth and causing the ice to advance further. Conversely, if the ice is receding, more energy is absorbed, heating Earth and melting more ice. The concept that ice-albedo feedback is a mechanism that ampliﬁes the eﬀect of Milankovitch cycles on climate originated in the mid–nineteenth century as well . The ice-albedo feedback can be modeled using ideas dating back to Budyko  and Sellers . The model is a dynamical system with a state space consisting of the annual average Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 686 R. McGEHEE AND C. LEHMAN temperature as a function of latitude. The energy input is the annual insolation as a function of latitude. Ice is assumed to form if the annual average temperature falls below a certain critical value and to melt if it rises above the critical value. The amount of insolation absorbed at a given latitude depends on whether there is ice at that latitude. Reradiation is a function of Earth surface temperature and hence depends on latitude. Heat transfer between latitudes is modeled simply as a linear transfer of energy from the warmer latitudes to the cooler ones. The details of this model are described in section 4. The annual average temperature as a function of latitude and time can be computed using the Milankovitch cycles to determine the input to the ice-albedo feedback model. The location of the ice boundary and the global mean surface temperature (GMT) as functions of time can also be computed. The model predicts that the ice boundary and the GMT depend mostly on the obliquity, followed by eccentricity, with no dependence on precession. Thus the model is consistent with the last 4,500,000 years of data, as analyzed by Zachos et al. (See Figure 1.) It is also consistent with recent studies showing that the glacial cycles are triggered by obliquity  and that they are driven by the summer average insolation . Although the model correctly predicts the relative contribution of the components of the Milankovitch cycles to Earth’s climate, there are large discrepancies between the predicted climate and the data, suggesting that other mechanisms play a larger role than that of the albedo or perhaps a synergistic role, amplifying the eﬀect of ice-albedo feedback. 2. The climate record. On the glacial time scale, the most complete climate record is provided by ocean sediment data. Small organisms called foraminifera form calcium carbonate shells which are fossilized in the sediments. The concentration of oxygen isotope 18 O in the shells was ﬁxed at the time the organism was living and depends on both the isotope concentration and the temperature of the sea water. Since the 18 O isotope is heavier than the more common 16 O isotope, water molecules evaporating from the ocean are less likely to contain 18 O than those left behind in the ocean. Therefore, the concentration of 18 O in sea water is higher during glacial periods, when a larger portion of 16 O is locked up in the glaciers. Also, the uptake of 18 O by the foraminifera is higher when the water temperature is colder. As a result, the 18 O concentration in the fossils is a proxy for climate: higher 18 O concentrations correspond to colder climates. Since the eﬀect due to temperature is smaller than the eﬀect due to concentration during the glacial cycles , we assume that the 18 O fossil record primarily indicates ice volume. The 18 O record for the last 5 million years is shown in Figure 2. The vertical axis labeled 18 δ O indicates a change from a ﬁxed standard. Note that the vertical scale is reversed; δ18 O decreases vertically, so that values higher on the vertical axis indicate higher temperatures. Data are taken from Lisecki and Raymo . It will be convenient to use the geologic time scale for various time periods represented in Figure 2. The Pliocene Epoch runs from about 5.3 million years ago to about 2.6 million years ago. The Pleistocene Epoch runs from the end of the Pliocene to about 12,000 years ago. We will use the term “early Pliocene” for the beginning of the Pliocene up until about 3.6 million years ago, and the term “late Pleistocene” for the last million years of the Pleistocene. These time periods are indicated in Figure 2. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 687 Late Pleistocene EarlyPliocene 2.5 3.0 ɷ18O 3.5 4.0 4.5 5.0 5.5 Ͳ5500 Ͳ5000 Ͳ4500 Ͳ4000 Ͳ3500 Ͳ3000 Ͳ2500 Ͳ2000 Ͳ1500 Ͳ1000 Ͳ500 0 Kyr Figure 2. Climate record for last 5.32 million years (data from ). 3. Milankovitch cycles. The incoming solar radiation at a point on Earth’s surface varies with time, most notably, daily cycles of night and day and yearly cycles of changing seasons. There are also longer-term cycles caused by variations in Earth’s orbit and rotation axis. These are the so-called Milankovitch cycles, whose three components are (1) eccentricity, (2) obliquity, and (3) precession. 3.1. Eccentricity. Earth moves around the sun on an approximately elliptical orbit. The orbit changes over geologic time due primarily to the inﬂuence of the other planets. As will be discussed below, the main component of Earth’s orbit aﬀecting the global annual insolation, averaged over the entire surface of Earth over an entire year, is the eccentricity of the ellipse. Over the last 5.5 million years, the eccentricity has ranged from almost zero (a nearly circular orbit) to about 0.06. The current value is about 0.017. Using the equations of celestial mechanics, Earth’s orbit can be computed both backward and forward in time. The eccentricity for the past 5.5 million years is shown in Figure 3. Here we are using computations due to Laskar et al. . 0.06 eccentricity 0.05 0.04 0.03 0.02 0.01 0 Ͳ5500 Ͳ5000 Ͳ4500 Ͳ4000 Ͳ3500 Ͳ3000 Ͳ2500 Ͳ2000 Ͳ1500 Ͳ1000 Ͳ500 0 Kyr Figure 3. Eccentricity over the last 5.5 million years. The graph shows a period of about 100 kyr superimposed on a period of about 400 kyr. The power spectrum of this time series is shown in Figure 4, which is computed using the discrete Fourier transform algorithm in MATLAB. Note that the apparent 100 kyr cycle breaks into a 95 kyr cycle and a smaller 125 kyr cycle. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 688 R. McGEHEE AND C. LEHMAN Figure 4. The power spectrum of the eccentricity over the last 5.5 million years. obliquity(degrees) 24.5 24.0 23.5 23.0 22.5 22.0 Ͳ5500 Ͳ5000 Ͳ4500 Ͳ4000 Ͳ3500 Ͳ3000 Ͳ2500 Ͳ2000 Ͳ1500 Ͳ1000 Ͳ500 0 Kyr Figure 5. Obliquity over the last 5.5 million years. 3.2. Obliquity. The obliquity is the angle between Earth’s rotational axis and the vector perpendicular to the plane of the ecliptic (the plane containing Earth’s orbit). The obliquity varies over geologic time between about 22◦ and about 24.5◦ . The current value is about 23.5◦ . Changes in the obliquity have an important eﬀect on Earth’s climate because the intensity of the solar input near Earth’s poles depends on the obliquity. Higher values of the obliquity produce more solar energy hitting the North Pole during the northern summer and the South Pole during the northern winter. The obliquity can also be computed using equations of classical mechanics. Figure 5 shows the obliquity of Earth’s rotation axis for the last 5.5 million years. Again we are using the computations of Laskar et al. . The graph shows a dominant frequency with a modulating amplitude. The power spectrum is shown in Figure 6, where we see that the dominant frequency has a period of about 41 kyr. The long-period modulation seen in Figure 5 is probably a beat between the two close frequencies seen in Figure 6. 3.3. Precession. The axis of Earth’s rotation changes due to small diﬀerences in the gravitational force on Earth’s equatorial bulge. Like a spinning top, Earth’s rotation axis precesses about the vector perpendicular to the plane of the ecliptic. The actual precession in a ﬁxed reference plane is not the important variable for Earth’s weather and climate. Instead, the appropriate precession variable is the longitudinal angle between the rotation vector and the major axis of Earth’s elliptical orbit. This angle determines Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 689 Figure 6. The power spectrum of the obliquity over the last 5.5 million years. precessionindex 0.06 0.04 0.02 0 Ͳ0.02 Ͳ0.04 Ͳ0.06 Ͳ5500 Ͳ5000 Ͳ4500 Ͳ4000 Ͳ3500 Ͳ3000 Ͳ2500 Ͳ2000 Ͳ1500 Ͳ1000 Ͳ500 0 Kyr Figure 7. Precession index over the last 5.5 million years. Figure 8. The power spectrum of the precession index over the last 5.5 million years. where the seasons occur along the ellipse and aﬀects the relative insolation during winter and summer months. In the computation of the insolation quantity Q65 , the precession shows up as the product of the eccentricity and the sine of the precession angle, called the precession index. This graph is shown in Figure 7. The amplitude modulation is the eﬀect of multiplying by the eccentricity. The underlying precession cycle has a period of about 23 kyr, although it is actually composed of three dominant periods, as shown in Figure 8. As is the case for Figures 5 and 6, the amplitude Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 690 R. McGEHEE AND C. LEHMAN 580 W/m2 540 500 460 420 Ͳ5500 Ͳ5000 Ͳ4500 Ͳ4000 Ͳ3500 Ͳ3000 Ͳ2500 Ͳ2000 Ͳ1500 Ͳ1000 Ͳ500 0 Kyr Figure 9. Daily average insolation at 65◦ N on summer solstice. Figure 10. The power spectrum of the daily insolation at 65◦ N on summer solstice for the last 5.5 million years. modulation shows up as a beat between the frequencies seen in Figure 8. 3.4. Daily insolation at 65◦ N. All of these variables together aﬀect Earth’s weather and climate. It is sometimes useful to measure their eﬀect on insolation by considering a single variable, the daily average insolation at 65◦ north latitude on the summer solstice, denoted Q65 . This quantity is a function of the three Milankovitch variables and hence varies with time. Its values for the last 5.5 million years are shown in Figure 9. A superﬁcial inspection of the graph leads to the speculation that this insolation closely tracks the precession index. This observation is conﬁrmed by the power spectrum, as can be seen in Figure 10. Note the prevalence of the frequencies clustered around a period of 23 kyr, in a pattern very similar to that of the precession index seen in Figure 8. Note also the frequencies around a period of 41 kyr, in the same pattern as those of the obliquity seen in Figure 6. Finally, note the lack of any discernible frequencies associated with eccentricity. 3.5. The climate record. In a highly cited work , Hays, Imbrie, and Shackleton related the Milankovitch cycles to the climate of the last 468,000 years, using data from ocean sediment core samples. Analyzing the power spectrum of the data, they concluded that “. . . climatic variance of these records is concentrated in three discrete spectral peaks at periods of 23,000, 42,000, and approximately 100,000 years. These peaks correspond to the dominant periods of the Earth’s solar orbit, and contain respectively about 10, 25, and 50 percent of the climatic variance.”  Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 691 On the other hand, analyzing the power spectrum of the solar forcing, they found the opposite order for the variance of the forcing, as illustrated in Figure 1. It should be noted that they analyzed as solar forcing both the insolation at 60◦ N at summer solstice and the insolation at 50◦ S at winter solstice, instead of the usual 65◦ N at summer solstice. However, as we saw in the previous section and illustrated in Figure 10, the conclusion holds as well for 65◦ N: the dominant variance occurs at periods associated with precession, followed by periods associated with obliquity. Periods associated with eccentricity contain a negligible portion of the variance. In a more recent work, Zachos et al.  repeated the previous approach using much more extensive data. In the record of the last 4.5 million years, they found the dominant periods in the power spectrum of the data to be those associated with obliquity, followed by those associated with eccentricity. This conclusion is illustrated in Figure 1. The climate data used by Zachos et al. is essentially the same as that shown above from Lisiecki and Raymo. The power spectrum is shown in Figure 11. Note the dominance of the frequency corresponding to a period of 41 kyr, and note the resemblance to the spectrum of the obliquity shown in Figure 6. Note also the suite of frequencies associated with periods around 100 kyr, corresponding to the spectrum associated with eccentricity and shown in Figure 4. We also see a small suite of frequencies around 23 kyr, presumably a small contribution from precession. A linear trend was subtracted from the data before the power spectrum was computed using MATLAB’s discrete Fourier transform function. 41Kyr 23Kyr power 100Kyr 0.00 0.01 0.02 0.03 0.04 0.05 0.06 frequency(1/Kyr) Figure 11. The power spectrum of the Lisiecki–Raymo stack. As mentioned in the introduction, several theories have been proposed to explain the dominance of the obliquity cycles in the climate data. The model analyzed below provides a possible explanation based on the eﬀect of ice-albedo feedback. The parameters in the model are all annual averages. As will be seen below, the contribution of precession cancels when the average annual insolation as a function of latitude is computed. A similar eﬀect occurs when computing the seasonal average insolation [4, 6]. Although theories based on seasonal averages probably produce more complete explanations of the process governing glacial retreats, it seems to us that a simpler model may also have some merit. 4. Ice-albedo feedback. Models exploring the eﬀect of ice-albedo feedback on Earth’s ice cover date back to Budyko  and Sellers . Many authors have explored similar models. Here we use a version due to Tung . Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 692 R. McGEHEE AND C. LEHMAN The basic variable is the annual average surface temperature T as a function of latitude. The dynamical equation can be written ∂T = Qs (y) (1 − α (y, η)) − (A + BT ) + C T̄ − T , ∂t where y is the sine of latitude and T is the annual average surface temperature as a function of y and time t. This equation is an example of an energy balance model, since it is of the form ∂T = heat imbalance. R ∂t The units on both sides of the equation are Watts per square meter (Wm−2 ). The quantity R is the speciﬁc heat of Earth’s surface, measured in units of Watts per square meter per degree centigrade. The actual value of R is irrelevant in this paper, because we consider only equilibrium solutions. The right-hand side of (1) breaks into (1) R heat imbalance = insolation − reﬂection − reradiation + transport. 4.1. Insolation. The annual average incoming solar radiation is given by the term insolation = Qs (y). The quantity Q is the global annual average insolation, while s (y) is the relative insolation normalized to satisfy 1 s (y) dy = 1. (2) 0 We can think of s (y) as the distribution of the incoming solar radiation as a function of latitude, averaged over an entire year. The variable y is chosen instead of the actual latitude so that, by Archimedes, the annual GMT can be written simply as 1 T (y, t) dy. T̄ (t) = 0 Note that we are assuming symmetry with respect to the equator, so the variable y takes values between 0 and 1. 4.2. Reflection. The incoming solar energy reﬂected back into space is given by the term reﬂection = Qs (y) α (y, η) . It is assumed that there is a single ice line at y = η, that ice covers the surface for y > η, and that the surface is ice-free for y < η. The albedo, α (y, η), has one value, α1 , where the surface is ice-free, and another value, α2 , where the surface is ice covered. Thus, α1 , y < η, (3) α (y, η) = α2 , y > η. Since the albedo is the proportion of energy reﬂected back into space, the annual rate at which solar energy is absorbed by Earth’s surface is Qs (y) (1 − α (y, η)). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 693 4.3. Reradiation. The energy reradiated into space at longer wavelengths is approximated linearly by the term reradiation = A + BT. This term includes many phenomena, all wrapped into a single linear term. The basic eﬀect is the emission of heat from Earth’s surface, which depends on the surface temperature. But before the heat escapes into space, some of it is absorbed by greenhouse gasses and returned to the surface. The reradiation term A + BT is the net loss of energy from the surface to space. The parameters A and B used below were determined empirically from satellite data . 4.4. Transport. The energy transported from warmer latitudes to cooler latitudes is approximated by the term transport = C T̄ − T . Once again, many phenomena are included in this single linear term, mostly atmospheric and oceanic circulation. Since all these phenomena are averaged over an entire year, it is perhaps not a terrible simpliﬁcation to assume that all parts of the surface are attempting to reach the GMT and that the annual heat transfer is proportional to the diﬀerence between the annual GMT and the annual mean temperature at a particular latitude. 4.5. Equilibrium solution. Following Tung , we use these values for the constants introduced above: Q = 343 Wm−2 , A = 202 Wm−2 , B = 1.9 Wm−2 (◦ C)−1 , C = 3.04 Wm−2 (◦ C)−1 , α1 = 0.32, α2 = 0.62. We now look for an equilibrium solution Tη∗ (y) having a single ice line at y = η. This equilibrium will satisfy (4) Qs (y) (1 − α (y, η)) − A + BTη∗ (y) + C T̄η∗ − Tη∗ (y) = 0. Integrating (4) yields 1 0 Qs (y) (1 − α (y, η)) − A + BTη∗ (y) + C T̄η∗ − Tη∗ (y) dy = 0, which simpliﬁes to Q (1 − ᾱ (η)) − A − B T̄η∗ = 0, (5) where ᾱ (η) = 1 0 α (y, η) s (y) dy = η 0 α1 s (y) dy + η 1 α2 s (y) dy. Using the normalization of s (y) given by (2), we can write ᾱ (η) = α2 − (α2 − α1 ) S (η), Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 694 R. McGEHEE AND C. LEHMAN where S(η) = 0 η s(y)dy. Thus, if we know the ice line η, we can solve (5) for the GMT, T̄η∗ = (6) 1 (Q (1 − ᾱ (η)) − A), B and, using this equation and (4), we can solve for the equilibrium temperature proﬁle, Tη∗ (y) = (7) 1 Qs (y) (1 − α (y, η)) − A + C T̄η∗ , B+C where α (y, η) is given by (3) and T̄η∗ is given by (6). Note that the discontinuity of the albedo produces an equilibrium temperature proﬁle that is discontinuous across the ice boundary, despite the assumption of heat transport between latitudes. The heat transport is not modeled as a diﬀusion process, which would produce a continuous temperature proﬁle, but is instead modeled as a relaxation toward the GMT. Continuing to follow Tung , we assume that, at equilibrium, the average temperature across the ice line is Tc = −10◦ C. Since s is continuous, we have 1 B+C 1 Tη∗ (η+) = B+C Tη∗ (η−) = Qs (η) (1 − α1 ) − A + C T̄η∗ , Qs (η) (1 − α2 ) − A + C T̄η∗ , so the ice line condition becomes (8) Tc = Tη∗ (η−) + Tη∗ (η+) 1 = Qs (η) (1 − α0 ) − A + C T̄η∗ , 2 B+C where α0 = α1 + α2 . 2 Combining (8) and (6), we have C 1 Qs (η) (1 − α0 ) − A + (Q (1 − ᾱ (η)) − A) = Tc , B+C B which reduces to (9) Q B+C A C s (η) (1 − α0 ) + (1 − ᾱ (η)) − − Tc = 0. B B This is the equation we will solve numerically for the location of the ice line η. Once we have computed the location of the ice line, the GMT follows from (6), and the temperature proﬁle follows from (7). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 695 5. Insolation. We turn now to the insolation Qs (y) as it appears in the ice-albedo feedback model. We will compute the global insolation Q and the insolation distribution function s (y) in terms of Earth’s eccentricity, obliquity, and precession. Let K be the solar output in Watts (Joules per second). At a distance r from the sun, the solar intensity is K Wm−2 . 4πr 2 The average annual insolation is diﬀerent for diﬀerent points on Earth’s surface. To compute these quantities, we use spherical coordinates to describe a point on the surface of Earth: ⎤ ⎡ ⎤ ⎡ cos ϕ cos γ u1 u = ⎣u2 ⎦ = ⎣ cos ϕ sin γ ⎦ , u3 sin ϕ where γ is the longitude and ϕ is the latitude. Note that u is a unit vector perpendicular to the surface of Earth. We introduce coordinates (x̂, ŷ, ẑ) with the sun at the origin. The (x̂, ŷ)-plane is the plane of Earth’s orbit (the ecliptic), the x̂-axis is the major axis of Earth’s elliptical orbit, and the positive x̂-axis points to the aphelion (the point on the ellipse furthest from the sun). We will refer to these coordinates as the “ecliptic coordinates.” We will also use polar coordinates in the (x̂, ŷ)-plane, assuming that the position of Earth in the ecliptic coordinates is (r (t), θ (t), 0). The following orthogonal matrix rotates the earth to an obliquity angle of β: ⎡ ⎤ cos β 0 sin β 1 0 ⎦. U1 (β) = ⎣ 0 − sin β 0 cos β It will be useful below to introduce the variables û, ϕ̂, and γ̂ by ⎡ ⎤ cos ϕ̂ cos γ̂ û = ⎣ cos ϕ̂ sin γ̂ ⎦ = U1 (β) u . sin ϕ̂ One can think of ϕ̂ and γ̂ as latitude and longitude on Earth’s surface measured not with respect to the axis of rotation but with respect to a vector perpendicular to the ecliptic plane. We also introduce the following matrix, which rotates the earth to a precession angle of ρ: ⎡ ⎤ cos ρ − sin ρ 0 U2 (ρ) = ⎣ sin ρ cos ρ 0⎦ . 0 0 1 In the ecliptic coordinates, the point on the surface of Earth becomes ⎡ ⎤ ⎡ ⎤ r cos θ r cos θ ⎣ r sin θ ⎦ + U2 (ρ) û = ⎣ r sin θ ⎦ + U2 (ρ) U1 (β) u . 0 0 Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 696 R. McGEHEE AND C. LEHMAN The insolation at this point is proportional to the cosine of the angle between the vector pointing from the sun to Earth and the vector perpendicular to Earth’s surface. Because of the way we have chosen the coordinates, the scalar product of these two vectors will be positive at night, when the insolation is zero, and negative during the day, when the insolation is positive. Therefore the insolation (in Wm−2 ) at the point on Earth’s surface is the positive part of K cos θ sin θ 0 U2 (ρ) U1 (β) u. (10) − 2 4πr If we multiply out the matrices, this quantity becomes K (cos ϕ (cos β cos (θ − ρ) cos γ + sin (θ − ρ) sin γ) + sin ϕ sin β cos (θ − ρ)) . 4πr 2 Therefore, the insolation at any given point on Earth’s surface (ϕ, γ) (latitude, longitude) at any given point (r, θ) along the orbit is K − cos ϕ cos β cos(θ − ρ) cos γ + sin(θ − ρ) sin γ I(β, ρ, r, θ, ϕ, γ) = 2 4πr + − sin ϕ sin β cos(θ − ρ) , − where β is the obliquity and ρ is the precession angle. We are interested in the annual average. Backing up to expression (10), we can write + K cos θ sin θ 0 U2 (ρ) û . I= − 4πr 2 The insolation then becomes the positive part of ⎡ ⎤T ⎡ ⎤ cos θ cos ρ + sin θ sin ρ cos ϕ̂ cos γ̂ K ⎣ − cos θ sin ρ + sin θ cos ρ⎦ ⎣ cos ϕ̂ sin γ̂ ⎦ − 4πr 2 0 sin ϕ̂ K cos ϕ̂ cos (θ − ρ − γ̂) . =− 4πr 2 Thus we can write the insolation at any given point on Earth’s surface as K cos ϕ̂[− cos (θ − ρ − γ̂)]+ , 4πr 2 where ϕ̂ is the latitude and γ̂ is the longitude of the point on Earth’s surface with respect to the ecliptic coordinates and where ρ, r, and θ are as before. Note that the explicit dependence on the obliquity has disappeared in these coordinates. We now think of eccentricity, obliquity, and precession as constant and integrate the insolation over a year. The only variables depending on time are r and θ, which move along an ellipse. The yearly average is K 1 P ¯ cos ϕ̂[− cos (θ (t) − ρ − γ̂)]+ dt I (ρ, ϕ̂, γ̂) = P 0 4πr(t)2 K cos ϕ̂ P [− cos (θ (t) − ρ − γ̂)]+ dt, = 4πP r(t)2 0 I (ρ, r, θ, ϕ̂, γ̂) = Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 697 where P is the period of Earth’s orbit. Recalling Kepler’s second law of planetary motion, we introduce the orbital speciﬁc angular momentum, Ω = r(t)2 dθ , dt and write the integral in terms of θ: K cos ϕ̂ I¯ (ρ, ϕ̂, γ̂) = 4πP Ω K cos ϕ̂ = 4πP Ω 0 2π [− cos (θ − ρ − γ̂)]+ dθ π/2 −π/2 cos θdθ = K cos ϕ̂ . 2πP Ω Note that the precession angle and the longitude both disappear completely from this equation. The annual average insolation depends only on the period and the angular momentum of Earth’s orbit and the latitude of the point on Earth’s surface, measured in the ecliptic coordinates. Returning to Earth’s coordinates, recall that û = U1 (β) u , which implies that sin ϕ̂ = − sin β cos ϕ cos γ + cos β sin ϕ and gives us an annual average insolation of K K cos ϕ̂ = 2πP Ω 1 − (− sin β cos ϕ cos γ + cos β sin ϕ)2 2πP Ω . We have so far completely ignored Earth’s rotation, focusing instead on a particular latitude and longitude on a nonrotating Earth. If we now average over the longitude γ, we ﬁnd the annual average insolation at a given latitude on a rotating Earth: 2π K 1 1 − (sin β cos ϕ cos γ − cos β sin ϕ)2 dγ 2πP Ω 2π 0 2π K 1 − (sin β cos ϕ cos γ − cos β sin ϕ)2 dγ . = 2 4π P Ω 0 I¯ = Note that this quantity is independent of the precession angle ρ. Note also that it is an even function of ϕ. In the ice-albedo model, this quantity I¯ that we just 1computed is the same as Qs (y), where y = sin ϕ and the function s is normalized so that 0 s (y) dy = 1. Therefore, K Qs (y) = 2 4π P Ω 2π 0 1− 1 − y 2 sin β cos γ − y cos β Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 2 dγ. 698 R. McGEHEE AND C. LEHMAN 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Figure 12. Insolation distribution as a function of y. Integrating over y, we have 1 s (y) dy 2Q = Q −1 K = 2 4π P Ω K = 2 4π P Ω 1 2π −1 0 π/2 2π −π/2 0 1− 1 − y 2 sin β cos γ − y cos β 2 dγdy 1 − (cos ϕ sin β cos γ − sin ϕ cos β)2 dγ cos ϕdϕ. Returning to the ecliptic coordinates (ϕ̂, γ̂) and noting that cos ϕ dϕdγ = cos ϕ̂ dϕ̂dγ̂, we have π/2 2π K cos ϕ̂ cos ϕ̂ dγ̂dϕ̂ 2Q = 2 4π P Ω −π/2 0 π/2 K K π K = , cos2 ϕ̂ dϕ̂ = = 2πP Ω −π/2 2πP Ω 2 4P Ω and so Q= K 8P Ω 2π 2 2 1− 1 − y 2 sin β cos γ − y cos β dγ. s (y) = 2 π 0 The graph of this function is shown in Figure 12 for the current value of obliquity β = 23.5◦ . Kepler’s third law tells us that P ∼ a−3/2 , and where a is the semimajor axis of Earth’s orbit. We also have that Ω ∼ a−1/2 1 − e2 , where e is the eccentricity of Earth’s orbit. According to Laskar et al. , the semimajor axis is essentially constant. Therefore, we can write Q = Q (e) = Q0 K =√ . 8P Ω 1 − e2 Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 699 Using a current eccentricity of 0.0167 and the current value of Q = 343 given by Tung , we see that Q0 = 343 1 − 0.01672 = 342.95, which is the value we use in the computations below. 6. Model simulation. 6.1. Equilibrium ice line and GMT. In section 4 we derived the ice-albedo model (9) which can be solved for the ice line y = η. For convenience, we repeat it here: Q (e) A C (11) h (η, e, β) = s (η, β) (1 − α0 ) + (1 − ᾱ (η, β)) − − Tc = 0, B+C B B where ᾱ (η, β) = α2 − (α2 − α1 ) η 0 s (y, β) dy. Note that we are explicitly displaying the dependence on Earth’s eccentricity e and obliquity β. We also derived (6), which expresses the GMT in terms of the location of the ice line: (12) T̄η∗ (e, β) = 1 (Q (e) (1 − ᾱ (η, β)) − A) . B In the previous section, we computed the global annual mean insolation Q and the insolation distribution function s (y) in terms of e and β: Q0 , Q (e) = √ 1 − e2 2π 2 2 1− 1 − y 2 sin β cos γ − y cos β dγ. s (y, β) = 2 π 0 These equations are inserted into (11) to get an implicit equation for the ice line η in terms of the eccentricity and obliquity. The current value for Earth’s eccentricity is e = 0.0167, while the current value of the obliquity is β = 23.5◦ = 0.410 radians. For these values, the graph of h (η, e, β) is shown in Figure 13 as the solid line. Note the two zeros, one at about η = 0.92 and one at about η = 0.26. The higher value corresponds to a stable ice line, while the lower value is unstable. For a discussion of the stability, see Tung  or Widiasih . Figure 13 also shows the graph of h for two other values of the parameters. The dotted line corresponds to eccentricity e = 0 and obliquity β = 22◦ = 0.384 radians, while the dashed line corresponds to eccentricity e = 0.06 and obliquity β = 24.5◦ = 0.428 radians. These values were chosen to illustrate the extremes over the glacial cycles. It is perhaps surprising that the variation of the stable equilibrium over glacial extremes is so small, given that, at the glacial maxima, the ice covered about 30% of Earth’s surface. This small variation indicates the existence of feedback mechanisms other than ice-albedo and is discussed below. Laskar et al.  provide tables for the eccentricity and obliquity for the past 5.5 million years. For each value of eccentricity and obliquity, we use Newton’s method to solve (11) for Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 700 R. McGEHEE AND C. LEHMAN 10 5 0 Ͳ5 Ͳ10 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Figure 13. The graph of h as a function of η. The solid line is for current values of eccentricity and obliquity, the dotted line is for e = 0 and β = 22◦ , and the dashed line is for e = 0.06 and β = 24.5◦ . The zeros correspond to equilibria of the ice line. IceLine 0.94 0.93 0.92 0.91 Ͳ5500 Ͳ5000 Ͳ4500 Ͳ4000 Ͳ3500 Ͳ3000 Ͳ2500 Ͳ2000 Ͳ1500 Ͳ1000 Ͳ500 0 time(Kyr) Figure 14. Ice line location from the ice-albedo feedback model. the ice line η, which we can plug into (12) to derive the GMT. We use a continuation method to follow the stable value of η through changes in e and β. We are assuming that the diﬀerential equation (1) comes to equilibrium quickly compared with the geologic time scale, so that it is necessary only to compute the equilibrium solution for each value of eccentricity and obliquity. 6.2. The model results. The results of the computation of the ice line for the past 5.32 million years are shown in Figure 14. We can compare the model output to the actual data shown in Figure 2. Under the assumption that the ice line is closely related to ice volume and therefore to the δ18 O record, Figure 14 should bear some resemblance to Figure 2. At ﬁrst glance, these graphs seem unrelated. The diﬀerences will be discussed below, but ﬁrst we consider the similarities by inspecting the power spectrum, shown in Figure 15. Note that the dominant period is at 41 kyr, exactly the dominant period for the obliquity. Indeed, one sees small spikes at frequencies of about 0.019 about 0.035, which also correspond to spikes seen in the obliquity spectrum (Figure 6). Although barely visible in Figure 15, there are also small spikes at the periods around 400 kyr and 100 kyr, reﬂecting the eﬀect of the eccentricity. A comparison of this ﬁgure with Figure 8 shows an absence of frequencies associated with precession—not a surprise, since the eﬀect due to precession cancels out of the model, which considers only annual averages. Thus we see the dominance of obliquity, with a very small contribution from eccentricity, and an absence of precession. Although the actual powers are substantially diﬀerent, these relative contributions by the Milankovitch cycles are Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 100Kyr 41Kyr 23Kyr power 400Kyr 701 0 0.01 0.02 0.03 0.04 0.05 0.06 frequency(1/Kyr) Figure 15. Power spectrum of the ice line computation. 14.8 GMT(°C) 14.6 14.4 14.2 14 13.8 Ͳ5500 Ͳ5000 Ͳ4500 Ͳ4000 Ͳ3500 Ͳ3000 Ͳ2500 Ͳ2000 Ͳ1500 Ͳ1000 Ͳ500 0 time(Kyr) Figure 16. GMT as computed in the ice-albedo feedback model. 100Kyr 41Kyr 23Kyr power 400Kyr 0 0.01 0.02 0.03 0.04 0.05 0.06 frequency(1/Kyr) Figure 17. Power spectrum of the GMT computation. in exactly the same order as they are in the climate data, as shown in the power spectrum illustrated in Figure 11. Recall from our discussion in section 2 that, although the biggest contribution to the δ18 O record comes from the ice volume, the record is also aﬀected by the deep ocean temperature. The ice-albedo feedback model readily allows for the computation of the GMT. Although the computed GMT is related to surface temperature rather than to deep ocean temperature, it still may be reﬂected in the δ18 O record. Figure 16 shows the GMT computed from the model, and Figure 17 shows its power spectrum. Note that, although the GMT is still dominated by the obliquity cycles, there is a substantial contribution from eccentricity. The model seems to indicate that, although Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 702 R. McGEHEE AND C. LEHMAN the eccentricity has a smaller eﬀect on climate than the obliquity, its eﬀect on the GMT is stronger than its eﬀect on ice volume. Note that both the computed ice boundary and the computed GMT show the same relative contributions from the orbital parameters. Thus we have established the main point of the paper: a simple model of ice-albedo feedback is suﬃcient to explain the relative contributions of the Milankovitch cycles in the climate data. However, as noted above, the model is far from providing a complete explanation of the data. 7. Comparison of the model with the data. 7.1. Amplitude comparison. Although the ice-albedo feedback model provides an explanation for the frequencies found in the climate record, it does not do very well at explaining the amplitude. To compare the amplitudes, we need to transform the model output of the ice boundary to the climate data of δ18 O. We will compare each to sea level. First recall that η corresponds to area, measured in units where 1 corresponds to the surface area of Earth. Assuming that the ice sheets have a constant thickness, then η corresponds to ice volume. Since δ18 O also corresponds to ice volume, we assume that there is a linear relation between η and δ18 O. Also, approximating with a constant surface area for the ocean, the sea level corresponds to volume of water and hence, inversely, to the volume of ice. We use the estimate that, if all the ice sheets melted, the sea level would rise about 65 meters, corresponding to an ice boundary of η = 1. According to the model, the current value of η is about 0.92. Therefore, the conversion factor is 65/0.08 = 810 meters per unit area (whole globe) of ice loss. The current value of δ18 O is about 3.2, while the value at the last glacial maximum was about 5.0 . If we use the estimate that the sea level at the last glacial maximum was about 125 meters below what it is today, we have a conversion factor of 125/1.8 = 70 meters of sea level drop per unit of δ18 O. Thus the conversion factor from η to δ18 O is 810/70 = 11.7 units of δ18 O per unit area. Since the current value of δ18 O is 3.2, while the current value of η is 0.92, we have the conversion formula δ18 O = 3.2 − 11.7(η − 0.92), which we used to produce Figure 18. 2.5 3 į18O 3.5 4 4.5 5 5.5 Ͳ5500 Ͳ5000 Ͳ4500 Ͳ4000 Ͳ3500 Ͳ3000 Ͳ2500 Ͳ2000 Ͳ1500 Ͳ1000 Ͳ500 Kyr Figure 18. The model result (red) superimposed on the climate record (blue). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 0 ICE-ALBEDO FEEDBACK 703 2.5 3 į18O 3.5 4 4.5 5 5.5 Ͳ1000 Ͳ900 Ͳ800 Ͳ700 Ͳ600 Ͳ500 Ͳ400 Ͳ300 Ͳ200 Ͳ100 0 Kyr Figure 19. Detail of Figure 18 for the late Pleistocene. į18O 2.5 3 3.5 Ͳ5500 Ͳ5000 Ͳ4500 Ͳ4000 Ͳ3500 Kyr Figure 20. Detail of Figure 18 for the early Pliocence. Figure 18 illustrates how far the model is from actually replicating the climate record, especially for the late Pleistocene (shown in detail in Figure 19), indicating that there are other phenomena at work much more powerful than ice-albedo feedback. However, as illustrated in Figure 20, the model more closely replicates the amplitude for the early Pliocene. 7.2. Delay. If one examines Figures 19 and 20 carefully, one can see a slight delay between the model output and the δ18 O data. We can quantify this delay by maximizing the correlation between the delayed model output and the data. For the early Pliocene, this delay is about 2.5 kyr, as shown in Figure 21, while for the late Pleistocene, this delay is about 7 kyr, shown in Figure 22. Lisiecki and Raymo noted a similar delay between the obliquity cycles and the δ18 O data, increasing from the early Pliocene to the late Pleistocene . They attributed the increased delay to the increased ice volume in the late Pleistocene, noting that it would take longer for the massive amounts of ice to respond to changes in the insolation. Indeed, they hypothesized such delays as part of the “tuning” of their age model, so it is not surprising that these delays show up again in the cross-correlation of their data with the ice-albedo feedback model, which, as noted above, responds primarily to the obliquity cycles. Our analysis tracks only the equilibrium solution of the Budyko model. Widiasih introduced a dynamic ice line to this model , which would cause a delay between the obliquity forcing and the model response. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 704 R. McGEHEE AND C. LEHMAN correlation 0.6 0.4 0.2 0 Ͳ0.2 Ͳ5 0 5 10 15 20 delay(Kyr) Figure 21. Cross-correlation between the model output and the climate data as a function of the delay for the early Pliocene. For this and the following ﬁgure, positive delay means the data lags the model. correlation 0.6 0.4 0.2 0 Ͳ0.2 Ͳ5 0 5 10 15 20 delay(Kyr) Figure 22. Cross-correlation between the model output and the climate data as a function of the delay for the late Pleistocene. 7.3. Skewness. Even with just a cursory look at the data during the Pleistocene, one notices that the ice retreats faster than it advances. Huybers has quantiﬁed this observation by computing the statistical quantity of “skewness” for the derivative of the data . Skewness is a measure of the departure of the distribution from symmetry about the mean. Huybers found that the skewness of the derivative of the δ18 O data increases during the Pleistocene, indicating an increasing asymmetry between the rate of glaciation and the rate of deglaciation. Figure 23 shows the cumulative distribution function (cdf) for the derivative of the δ18 O data for the late Pleistocene. The data are normalized so that the mean is zero and the standard deviation is one, so the horizontal scale is standard deviation. The vertical axis is the probability that the derivative of the δ18 O data is below the corresponding value on the horizontal scale. Using the standard deﬁnition of coeﬃcient of skewness as the third moment scaled by the standard deviation, we ﬁnd that the skewness of the derivative of the δ18 O data in the late Pleistocene is −0.46, meaning that the slopes are skewed toward negative values. Recall that the δ18 O is a measure of ice volume, so a negative slope corresponds to decreasing ice volume, or deglaciation. Thus a negative skewness indicates that the deglaciations proceed faster than the glaciations. A simpler measure of the asymmetry of the data is “Pearson’s second moment of skewness,” Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 705 cumulativeprobability 1 0.75 0.5 0.25 0 Ͳ3 Ͳ2 Ͳ1 0 1 2 3 slope(standarddeviationsfrommean) Figure 23. The cdf for the derivative of the δ 18 O data during the late Pleistocene. The data are scaled so that the mean is zero and the standard deviation is one. The vertical line indicates the median. cumulativeprobability 1 0.75 0.5 0.25 0 Ͳ3 Ͳ2 Ͳ1 0 1 2 3 slope(standarddeviationsfrommean) Figure 24. The cdf for the derivative of the model output of the ice line during the late Pleistocene. which is just 3 (x̄ − x̃)/σ, where x̄ is the mean, x̃ is the median, and σ is the standard deviation. Since the data in the ﬁgure are scaled so that the mean is zero and the standard deviation is one, this skewness measure is simply −3x̃. The median is shown as the vertical line in Figure 23, occurring at 0.106. Thus Pearson’s second moment of skewness is −0.32, also indicating that the deglaciations are faster than the glaciations. Figure 24 is the corresponding graph for the ice line model. We expect no skewness from the model, since the Milankovitch cycles are symmetric and since the relatively small forcing amplitude means that the system can be expected to respond linearly. Indeed, we see none: the coeﬃcient of skewness is −0.001, while Pearson’s second moment of skewness is +0.03, both eﬀectively zero, meaning no diﬀerence between the rates of glaciation and deglaciation. Thus the model does not reproduce the shape of the δ18 O data for the late Pleistocene. However, the early Pliocene is a diﬀerent story. Figure 25 is the corresponding graph for the δ18 O data for this period. Here the coeﬃcient of skewness is +0.06, while Pearson’s second moment of skewness is +0.007, both eﬀectively zero. Figure 26 is the corresponding graph for the model ice line. Here the coeﬃcient of skewness is +0.001, and Pearson’s second moment of skewness is −0.03, again both eﬀectively zero. Thus, for the early Pliocene, neither the data nor the model indicates an asymmetry between glaciation and deglaciation. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 706 R. McGEHEE AND C. LEHMAN cumulativeprobability 1 0.75 0.5 0.25 0 Ͳ3 Ͳ2 Ͳ1 0 1 2 3 slope(standarddeviationsfrommean) Figure 25. The cdf for the derivative of the δ 18 O data during the early Pliocene. cumulativeprobability 1 0.75 0.5 0.25 0 Ͳ3 Ͳ2 Ͳ1 0 1 2 3 slope(standarddeviationsfrommean) Figure 26. The cdf for the derivative of the model output of the ice line during the early Pliocene. 8. Summary. The Budyko model of ice-albedo feedback accurately reproduces one of the main features of the climate data: the dominant frequency in the climate time series is the frequency associated with the obliquity cycles of the Milankovitch forcing. This feature extends for at least the last ﬁve million years, from the early Pliocene (5.3 Myr to 3.6 Myr ago) to the late Pleistocene (1 Myr to 12 kyr ago). This reproduction, of course, does not imply that ice-albedo feedback is the main cause of the dominance of obliquity in the climate data; it implies only that our analysis does not rule out ice-albedo feedback as a signiﬁcant factor in Earth’s climate. The model fails to reproduce many of the important properties of the climate data, particularly for the late Pleistocene. The amplitude of the δ18 O cycles during the late Pleistocene is much larger than that produced by the model. There is a delay of 7 kyr between the cycles produced by the model and the corresponding cycles present in the data. The data are strongly skewed toward rapid deglaciations and slow glaciations, while the model output is unskewed. On the other hand, the model more accurately reproduces the data from the early Pliocene. The amplitude is approximately correct, the delay is only 2.5 kyr, and neither the data nor the model output is signiﬁcantly skewed. The late Pleistocene is characterized by massive ice sheets on the continents of the northern hemisphere, while the ice of the early Pliocene was mostly in Antarctica . Clearly Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ICE-ALBEDO FEEDBACK 707 something changed between the Pliocene and the Pleistocene that had nothing to do with Milankovitch cycles. Our analysis indicates that ice-albedo feedback may be a factor for the Pleistocene glacial cycles, but it is not the dominant one. However, it might have played a much more signiﬁcant role during the early Pliocene. Acknowledgments. The authors thank Mary Lou Zeeman for reading various drafts of the manuscript and suggesting many important changes. Many other improvements were suggested by the referees. McGehee also thanks Jacques Laskar for enlightening conversations and for his hospitality during a visit to the Observatoire de Paris, and he thanks Ray Pierrehumbert for help with the geological interpretations. REFERENCES  M. I. Budyko, The eﬀect of solar radiation variations on the climate of the earth, Tellus, XXI (1969), pp. 611–619.  J. D. Hays, J. Imbrie, and N. J. Shackleton, Variations in the earth’s orbit: Pacemaker of the ice ages, Science, 194 (1976), pp. 1121–1132.  A. M. Hogg, Glacial cycles and carbon dioxide: A conceptual model, Geophys. Res. Lett., 35 (2008), L01701.  P. Huybers, Early pleistocene glacial cycles and the integrated summer insolation forcing, Science, 313 (2006), pp. 508–511.  P. Huybers, Glacial variability over the last two million years: An extended depth-derived agemodel, continuous obliquity pacing, and the pleistocene progression, Quaternary Sci. Rev., 26 (2007), pp. 37– 55.  P. Huybers and E. Tziperman, Integrated summer insolation forcing and 40,000-year glacial cycles: The perspective from an ice-sheet/energy-balance model, Paleoceanography, 23 (2008), PA1208.  J. Imbrie and K. P. Imbrie, Ice Ages: Solving the Mystery, Harvard University Press, Cambridge, MA, 1979.  J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A. C. M. Correia, and B. Levrard, A longterm numerical solution for the insolation quantities of the earth, Astron. Astrophys., 428 (2004), pp. 261–285.  L. E. Lisiecki and M. E. Raymo, A Pliocene-Pleistocene stack of 57 globally distributed benthic δ 18 O records, Paleoceanography, 20 (2005), PA1003.  D. Paillard, Glacial cycles: Toward a new paradigm, Rev. Geophys., 39 (2001), pp. 325–346.  J. R. Petit, J. Jouzel, D. Raynaud, N. I. Barkov, J.-M. Barnola, I. Basile, J. C. M. Bender, M. Davisk, G. Delaygue, V. M. K. M. Delmotte, M. Legrand, V. Y. Lipenkov, C. Lorius, L. Pepin, C. Ritz, E. Saltzman, and M. Stievenard, Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Nature, 399 (1999), pp. 429–436.  R. T. Pierrehumbert, Principles of Planetary Climate, Cambridge University Press, Cambridge, UK, 2010.  W. D. Sellers, A global climatic model based on the energy balance of the earth-atmosphere system, J. Appl. Meteorol., 8 (1969), pp. 392–400.  K. K. Tung, Topics in Mathematical Modeling, Princeton University Press, Princeton, NJ, 2007.  E. Widiasih, Dynamics of Budyko’s energy balance model, SIAM Appl. Dyn. Syst., submitted.  J. Zachos, M. Pagani, L. Sloan, E. Thomas, and K. Billups, Trends, rhythms, and aberrations in global climate 65 Ma to present, Science, 292 (2001), pp. 686–693. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.