...

A Paleoclimate Model of Ice-Albedo Feedback Forced by Variations in... Orbit Richard McGehee Clarence Lehman

by user

on
Category: Documents
5

views

Report

Comments

Transcript

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 confined 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 affect Earth’s climate [7]. In a seminal paper, Hays, Imbrie, and Shackleton [2] firmly
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 [2]. 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 [2]. (See Figure 1.) They ascribed this reversal of
the contributions to possible nonlinear effects in climate dynamics [2].
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. [16] reached a somewhat different conclusion. They found the strongest contribution
to climate response in frequencies associated with obliquity, followed by eccentricity, with
negligible contribution from precession [16]; 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 [2] 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 effects [3, 5, 10].
Other studies have replaced Q65 with the total insolation during the summer [4, 6].
One feedback mechanism is the ice-albedo effect. The reflectivity of Earth’s surface, called
the albedo, is different for ice than for water or land. If more of Earth is covered with ice, then
more of the solar energy is reflected back into space, leaving less to heat Earth. On the other
hand, if less ice covers Earth, more solar energy is absorbed. This effect constitutes positive
feedback. As the ice advances, more solar energy is reflected, 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 amplifies the
effect of Milankovitch cycles on climate originated in the mid–nineteenth century as well [7].
The ice-albedo feedback can be modeled using ideas dating back to Budyko [1] and Sellers [13]. 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 [5] and that they are driven by the summer average insolation [6].
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 effect 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 fixed 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 effect due to temperature is smaller
than the effect due to concentration during the glacial cycles [12], 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 fixed 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 [9].
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 [9]).
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 influence of the other planets. As
will be discussed below, the main component of Earth’s orbit affecting 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. [8].
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 effect 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. [8].
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 differences 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 fixed 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 affects 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 effect 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 affect Earth’s weather and
climate. It is sometimes useful to measure their effect 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 superficial inspection of the graph leads to the speculation that this insolation closely
tracks the precession index. This observation is confirmed 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 [2], 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.” [2]
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. [16] 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 effect 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 effect
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 effect of ice-albedo feedback on Earth’s ice
cover date back to Budyko [1] and Sellers [13]. Many authors have explored similar models.
Here we use a version due to Tung [14].
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 specific 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 − reflection − 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 reflected back into space is given by the term
reflection = 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 reflected 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 effect
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 [14].
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 simplification to assume that all parts of the surface are attempting to reach the
GMT and that the annual heat transfer is proportional to the difference between the annual
GMT and the annual mean temperature at a particular latitude.
4.5. Equilibrium solution. Following Tung [14], 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 simplifies 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 profile,
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 profile that is discontinuous across the ice boundary,
despite the assumption of heat transport between latitudes. The heat transport is not modeled
as a diffusion process, which would produce a continuous temperature profile, but is instead
modeled as a relaxation toward the GMT.
Continuing to follow Tung [14], 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 profile
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 different for different 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 specific 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 find 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. [8], 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 [14],
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 [14] or Widiasih [15].
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. [8] 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 differential 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 first glance, these graphs seem
unrelated. The differences will be discussed below, but first 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, reflecting the effect of
the eccentricity. A comparison of this figure with Figure 8 shows an absence of frequencies
associated with precession—not a surprise, since the effect 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 different, 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 affected 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 reflected 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 effect on climate than the obliquity, its effect on the GMT is
stronger than its effect 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 sufficient 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 [9]. 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 [9]. 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 [15], 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 figure, 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 quantified this observation
by computing the statistical quantity of “skewness” for the derivative of the data [5]. 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 definition of coefficient of skewness as the third moment
scaled by the standard deviation, we find 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 figure 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 coefficient of skewness is −0.001, while Pearson’s second moment of skewness is +0.03,
both effectively zero, meaning no difference 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 different story. Figure 25 is the corresponding graph for
the δ18 O data for this period. Here the coefficient of skewness is +0.06, while Pearson’s second
moment of skewness is +0.007, both effectively zero. Figure 26 is the corresponding graph for
the model ice line. Here the coefficient of skewness is +0.001, and Pearson’s second moment
of skewness is −0.03, again both effectively 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 five 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 significant
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 significantly 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 [16]. 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 significant 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
[1] M. I. Budyko, The effect of solar radiation variations on the climate of the earth, Tellus, XXI (1969),
pp. 611–619.
[2] 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.
[3] A. M. Hogg, Glacial cycles and carbon dioxide: A conceptual model, Geophys. Res. Lett., 35 (2008),
L01701.
[4] P. Huybers, Early pleistocene glacial cycles and the integrated summer insolation forcing, Science, 313
(2006), pp. 508–511.
[5] 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.
[6] 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.
[7] J. Imbrie and K. P. Imbrie, Ice Ages: Solving the Mystery, Harvard University Press, Cambridge, MA,
1979.
[8] 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.
[9] L. E. Lisiecki and M. E. Raymo, A Pliocene-Pleistocene stack of 57 globally distributed benthic δ 18 O
records, Paleoceanography, 20 (2005), PA1003.
[10] D. Paillard, Glacial cycles: Toward a new paradigm, Rev. Geophys., 39 (2001), pp. 325–346.
[11] 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.
[12] R. T. Pierrehumbert, Principles of Planetary Climate, Cambridge University Press, Cambridge, UK,
2010.
[13] 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.
[14] K. K. Tung, Topics in Mathematical Modeling, Princeton University Press, Princeton, NJ, 2007.
[15] E. Widiasih, Dynamics of Budyko’s energy balance model, SIAM Appl. Dyn. Syst., submitted.
[16] 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.
Fly UP