Climate Physics and Chemistry - Wiggle Matching and the Milankovitch Hypotheses

Class Notes: Wiggle Matching and the Milankovitch Hypotheses 1 Introduction The study of past climates encounters a number of severe problems not arising, or not as severe, in the study of the present climate state. The “big three” of these issues are: (1) the physical interpretation of sometimes extremely indirect climate proxies (e.g., are changes in [δ18O] really an accurate local thermometer?). (2) The accuracy of complex numerical models, that always contain a variety of errors when run over indefinitely long periods. (3) The problem of assigning dates to measurements that are commonly functions of depth and not time (the “age-model” problem). Here we look at the age-model problem and in particular, the inference that astronomical forcing controls much of climate change in general, and the Pleistocene ice ages in particular. Perhaps unexpectedly, this subject leads one into problems of human psychology and evolution. The human eye clearly evolved so as to be a powerful instrument for detecting weak patterns in noisy backgrounds. The rationale derives from a survival scenario: a human hunter in the jungle can in turn be stalked by a tiger. There are four main outcomes: (1) no tiger is present, and no tiger is seen; (2) no tiger is present, but a tiger is falsely perceived–leading to unnecessary termination of that day’s hunt. (3) A tiger is present and is correctly perceived, leading to successful escape; (4) A tiger is present but not detected, and the hunter is himself a meal. The reproductive cost of a false alarm (2–a false positive) is much less than a false negative (4– failure to detect). Thus people see patterns–everywhere. The most famous such patterns are the constellations to which the most ancient people assigned shapes and stories. A more modern example is in Fig. 1–showing on the “canali” (“lines” in Italian, not “canals”) as drawn by a sophisticated astronomer in the late 19th century. Also shown is a modern image from orbiting spacecraft. Many other such examples, where the eye picks out patterns that may or may not be real, exist in science. (Astronomers in particular have worried about this phenomeon; see the ref-erences in Wunsch, 2006). Statisticians have therefore developed objective tests to determine 1 Figure 1: On the left is a view of Mars drawn in 1894 by Giovanni Virginio Schiaparelli and on the right is a recent Mars satellite compositve (from NASA website). whether some pattern (a “signal” is likely real or an accident of data or the eye’s ability to see patterns everywhere). A closely related problem is the great difficulty most people have in making correct inferences from statistical data (Kahneman, et al., 1982). Even experts are fooled.1 2 Wiggle Matching In the paleoclimate literature, the most common manifestation of the wish to see patterns lies in the practice known to its critics as “wiggle-matching”–the hypothesis that if two records show similar fluctuations, that they are identical, and one is permitted to infer that the fluctuations can be aligned within any absolute dating uncertainty. Some examples are shown in Figs. 2 5. In some of these cases, it is known that the record excursions are pure accidents, and there is no relationship between them. In the remaining cases, it is an hypothesis that the records should be aligned. Often the hypothesis cannot be proved or disproved, and one must prevent the assumption from becoming a published fact. It is actually easy to show that under certain circumstances two unrelated records will inevitably appear to have related oscillations. That circumstance arises when two records have similar frequency content. (See Appendix 1 for a brief discussion of frequency content and what is known as spectral estimation.) If two records have similar frequency content, they must have, 1 The class is asked to answer some simple questions. (A) A game (the game of Peter and Paul) is played by flipping a true coin. If the coin comes up heads, Peter pays Paul $1. If tails, Paul pays Peter $1. Sketch a “typical” curve of Peter’s winnings through time. (B) You are told that a couple has two children, of which one is a girl. What is the probability that the second child is a boy? For (A) see Feller (1957); for (B) Gauch (2003). 2 Figure 2: Left panel is taken from Wunsch (1999) and shows two statistically independent records that show strong visual resemblance over finite time spans because they have similar frequency content. An objective test of similarity would reject the hypothesis that these records are related. Middle panel shows a plot of sunspot number and the Dow Jones Industrial average through time (author unknown). They do seem to track impressively well. Right hand panel, top, is from Brooks (1923) and shows the apparent correlation between Central African lake levels and sunspot numbers. Lower panel shows that the oscillations on the 11-year sunspot cycle had disappeared by the 1920s, and the entire character was different toward the end of the record. Figure 3: (Re-drawn from Hendy et al., 2002) showing the apparent correspondence between the δ 18 O record in Santa Barbara Basin and that in the GISP2 record. That an equivalent degree of high frequency variability exists in both records is evident; whether the oscillations actually correspond as the dashed lines indicate, is much less obvious. 3 Figure 4: Identification of supposedly corresponding events in the Hulu Cave record and in Greenland ice core (re-drawn from Alley, 2005). Notice, e.g., that the large excursions in the Hulu cave record near -45KY and -30KY have no counterpart in the Greenland record. Figure 5: Maximum monthly temperature records for Oxford, England taken from completely different time spans. But as the physics of temperature change there has probably remained fixed, the two records have similar frequency content. Note how easy it would be to align the various maxima and minima if there should be any uncertainty in the record timing (there isn’t any here). 4 on average, the same number of positive and negative excursions about the mean in any interval (see for example, Vanmarcke, 1983). Thus it is almost inevitable that one can draw lines between corresponding maxima and minima in any finite record interval. But demonstrating that the maxima and minima represent the same event requires much stronger evidence. (A bit more discussion is given in Wunsch, 2006). 3 The Milankovitch Hypotheses Another part of human psychology that has been much studied is the understandable wish that the world should be predictable. That is, when some event occurs, one should be able to rationalize why it occurred, and preferably be able to explain the chain of reasoning. Some events are, however, very difficult to discuss in causal terms: e.g., why a hurricane occurred in a particlar location on a particular day; why a field of ocean surface waves led to a breaker at a particular time and place. In principle, one might explain the latter as having been caused by a superposition of waves coming from several directions that just added up to become unstable at that moment and place. In practice, such events are much more easily understood in statistical terms, because the chain of events could never have been predicted in detail. The wish for causality pervades the discussion of paleoclimate, and it may well be true that much of what is inferred about the past depends upon a readily rationalizable chain of events. Again, the main issue is less that one can demonstrate that past climate is primarily stochastic in nature, and more that one should not rush to claim causality just because of its psychological appeal. That leads us to the hypothesis that major climate shifts arise from small shifts in the Earth’s orientation relative to the sun–as predictable a process as exists in nature. There are several Milankovitch hypotheses and not all authors are careful in defining which they are referring to: "It is widely accepted that climate variability on time scales of 103 to 105 years is driven primarily by orbital, or so-called Milankovitch, forcing." (McDermott et al., Science, 2001). "...it is now quite clear that orbital forcing played a key role in pacing glaciations during the Quaternary...." (Bradley, R. S., Paleoclimatology, Academic Press, 1999, p. 281) "The orbital theory of climate is the prevailing theory of glacial-interglacial climate change over tens of thousands to hundreds of thousands of years." (Cronin, T. M., Principles of Paleoclimatology, Columbia Un. Press, 1999, p. 131) 5 "...we confirm that moisture source temperature signal recorded in Vostok deuterium excess over the last 150ka fully reflects the obliquity time-varying relative contribu tion of low and high latitudes to Vostok precipitation." (F. Vimeux et al., Earth and Plan. Sci. Letts., 203, 2002, p. 829) "...a strong case has been made that on the time scale of tens of thousands of years, the Earth’s climate is being paced by the so-called Milankovitch cycles..." (W. Broecker, Earth Sci. Revs., 51, 137-154, 2000). Note that there have been some strongly dissenting views, including Winograd et al. (1992); Roe and Allen (1999); Wunsch (2004), but they have been largely ignored or rationalized away. 3.1 A Bit About Earth Orbit To understand these inferences, we need to briefly review the orbital motion of the Earth about the sun. Because insolation depends only upon the separation of an observer on the surface of the earth and the solar position, the only relationships needed are geometric. For this purpose, it is adequate to think of the sun as in orbit about the Earth; see Fig. 6) and I will tend to use this Ptolemaic point of view, although switching back and forth where it is convenient. The dynamics producing the orbital position of the sun are done separately. This approach is widely used in the closely related problems of ocean tides and geodesy. By Kepler’s laws, to a very high approximation, the sun appears to be in an elliptical orbit about the sun with the Earth at one of the foci and with the ellipse today having a small eccentricity e = 0.0167 (Smart, 1962). The plane of the ecliptic in which the sun appears to move is tilted relative to the Earth’s equator by about 23◦ (23◦ 260 21” in the year 2000, e.g. Siedelmann, 2006, P. 114). Suppose the solar orbit were perfectly circular, e = 0. Then as seen in Fig. 7, the northern hemisphere is tilted toward the sun with the northern hemisphere receiving radiation through a shorter atmospheric path than in the southern. During the remaining six months of the year, the situation is reversed. This change accounts for the seasons as we know them. With finite e, the Earth is closer to the sun during some times of the year than at others. At present, the time of closest approach corresponds to northern hemisphere (“boreal”) winter (now about 12 days after the winter solstice on 21 December) and is furthest from the sun during northern hemisphere summer. Thus the northern hemisphere gets slightly more radiation each day during winter than it would in a circular orbit, and slightly less in summer (reversed in the southern hemisphere). There are two effects: (1) the northern hemisphere winters receive a bit more radiation than they would be in a circular orbit, and the summers might therefore 6 Image courtesy of US Government. Figure 6: Geometry of the solar orbit taken about a fixed Earth. The vernal equinox, Υ, is the point where the sun crosses into the northern hemisphere each year, but moves slightly westward taking 26,000 years to complete 360◦ . Υ is used as a coordinate because it is important in climate, among other reasons. The point of summer solstice is 90◦ further in the orbit than Υ. At the present time, that near to the time of closest approach to the Earth, which moves slowly eastward toward the westward moving Υ. (From Bowditch, 1962). 7 Figure 7: (From ?). Basic elements of the orbit, showing the (greatly exaggerated) effects of eccentricity (and the true time scales are closer to 95,000 and 400,000y of tilt (obliquity) which controls the intensity of high and low latitude insolation, and of precession, which with non-zero eccentricity, produces a slowly changing amplitude of the annual cycle. Figure 8: Kepler’s second law says that equal areas are swept out in equal times, and so the Earth is moving most rapidly at the time of closest approach. (Or, as used here, the sun appears to be moving most rapidly.) Eccentricity is greatly exaggerated. 8 Image courtesy of US Government. Figure 9: The constellations in a band about the ecliptic plane. They can be used as a fixed reference from for measuring directions from the Earth (or sun). Thus Υ once lay in the constellation Aries, but has subseequently moved westward into Pisces. (Bowditch, 1962). be a bit cooler than they would be in a circular orbit, whereas in the southern hemisphere, the winters may be slightly colder and the summers slightly warmer than in the circular orbit. (2) Kepler’s law of equal-areas in equal times shows that when the earth is closest to the sun, its orbital motion is faster than when it is farther away (Fig. 8). Thus northern hemisphere winter is a few days shorter than northern hemisphere summer, and again reversed in the southern hemisphere. The degree of compensation of the higher insolation in northern hemisphere winter by the shorter season has to be determined in a detailed calculation (Huybers, 2006). Looked at it from the point of view of an observer on the Earth, the sun crosses a line overhead of the equator twice/year, once headed northward (the spring or vernal equinox) and once headed southward (fall or autumnal equinox). The line in space defined by the intersection of the orbital plane (Fig. 6) and the Earth’s equator is where the sun crosses the equator. Astronomers long ago chose the direction of the northern hemisphere vernal (spring) equinox as the origin of a coordinate system, and gave it the symbol Υ (upsilon), and the name “first point of Aries,” because in ancient times it lay in the direction of the constellation Aries (the Ram; see Fig. 9). If this were the whole story, the only forcing by the sun of the climate system would be strictly periodic with two components: (1) the day, (2) the year. The day would be (and is) defined as slightly longer than the Earth’s rotation period, by the amount required, on average over the year, for the sun to be overhead at noon (the sun advances in its orbit while the Earth is rotating and so there is a catch-up interval). The year would be defined as the interval in days between passage of the sun through Υ. Because of the complications that now enter, it is important to remember that daily and seasonal forcing are by far the largest periodic drivers of 9 the system–by many orders of magnitude, and that is easily forgotten. The system is already a bit confusing, when one recalls e.g., that above the polar circles, daylight lasts 24 hours during summer and vanishes altogether during winter, with increasing amounts of winter daylight as one moves southward toward the equator. Now consider that while the earth orbits in a near-perfect ellipse about the sun, it is not quite perfect. The major deviation from perfection is that the ellipse does not close exactly–the point of closest approach moving counter-clockwise (that is in the direction of the earth’s orbital movement) retreating from the advancing sun in the geocentric picture. Thus it takes slightly more than an orbital period for the sun to catch up to the retreating point of closest approach (direction in space of perihelion, p1 ). If one defines the orbital period as the “sidereal year” (today, 365 days, 6 hours, 9 minutes, 10s= 365.2563d; Seidelmann, 2006), then there is another, slightly longer, year defined by the interval for the sun to pass through the retreating perigee (the Earth is at “perihelion”, but the hypothetical orbiting sun is at perigee) of 365d, 6hours, 13min, 53s = 365.2596d, and called the “anomalistic year”. Perihelion moves through 360◦ in about 114,000 sidereal years so the deviation from a perfect ellipse in any given year is very small. It should be apparent that perigee, which now occurs in boreal winter with northern hemisphere winter getting slightly more insolation, will occur later and later in the year, eventually being in the spring, and then in northern hemisphere summer, reversing the present configuration. These periods slowly change over time. For example (Bowditch, 1962, Appendix D), the sidereal year can be written approximately as 365.24219879-0.6.14×10−8 (τ − 1900)d where τ is the year and the first term represents the value in 1900 (different authors use different reference dates). For the very long intervals of time of interest in climate, however, further corrections are necessary. This story is still not, however, quite complete. We have told it as though the direction of the Earth’s spin axis (strictly its principle axis of angular momentum) defined as a fixed direction in space. But in fact it moves: the Earth has an equatorial bulge, and as a result of the action of the moon and sun on that bulge, the spin axis describes a slow nearly circular motion in space called the “precession of the equinoxes” to distinguish it from the orbital precession just described. The period of that circular motion is about 26,000 years (often calculated in freshman mechanics courses). It has the effect of shifting the direction of Υ slowly clockwise in space, toward the oncoming orbiting sun with a period of about 26,000 sidereal years2 . If we measure the year in terms of the time between earth or sun passages through the moving Υ, it is then slightly shorter than the sidereal year at 365d, 5hours, 48min, 46s =365.2422d. and is called the “tropical” year. Usually the tropical year is the one people care about because it denotes 2 The motion is not completely regular. In particular there is a superimposed oscillation in direction called “nutation” with about an 18.6y period (see Siegemann, 2006). 10 the onset of northern hemisphere spring and so is the year most directly related to climate, even though it has no special dynamical significance. When simply a “year” is designated, it is usually the tropical one. Here, the day, hour, etc. are all defined as in normal civil time. Bowditch (1962) and Seidemann (2006), have good summary discussions of time-keeping. The fundamental unit if time is the second, defined in terms of an atomic standard. The “common year” has 365 days and leap year has 366 days, and whose mean might be taken as approximately 365.25 days (the Julian day) with leap-year rules used to keep the calendar nearly synchronous with the tropical year so that the seasons appear at the expected times. Note that leap years are skipped every century with an exception every 400 years, so the true mean civil year is not exactly 365.25d). The civil day is defined in terms of an average through the year; astronomers also distinguish days of various length. We now have a picture of the sun orbiting the earth, catching up to the retreating perigee point, and meeting the oncoming Υ. Thus p1 is advancing toward Υ. For climate, the issue is when perihelion occurs relative to springtime. With the two positions advancing towards each other, they meet approximately every 21,000 tropical years , and thus the intervals when northern hemisphere winters have larger insolation than the southern (as now) occur roughly every 21,000 years, not the much longer interval that one would get from orbital precession if there were no equinoctial precession, nor the 26,000 years one would get if there were equinoctial precession and no orbital precession. An important point is that there is no forcing periodicity near 21,000 years period, rather the amplitude of the seasonal cycle in the northern and southern hemisphere varies on this period and the amplitude variations are 180 ◦ out of phase across the equator. Because of the orbital motions, there are two dominant frequency bands near 23,000 and 19,000 years. As shorthand, these will be referred to together as the 21,000 year-band. There are a couple more complications. The seasons exist because of the tilt (obliquity) of the earth’s axis relative to the orbital plane. This tilt varies slightly (about ±1◦ ) with periods clustering near 41,000 years. (Keep in mind that nothing is actually strictly periodic–orbital motions are too complicated for that. See Appendix 1 for a discussion of narrow-band processes of different types.) When the tilt increases, there is an increase in radiation at both poles on the annual average, and a consequent slight diminution at low latitudes clustered with periods around 41,000y (crossover latitude about 44◦ latitude). There is one more effect to account for: the eccentricity e varies between about 0 and 0.06 on time scales of about 400,000 and 100,000 years. If the eccentricity vanishes, recall that one loses the 21,000 year-band modulation of the annual cycle. The annual average radiation effects of changing eccentricity are essentially negligible at 100 and 400ky–despite a large literature claiming the opposite; we will return to that later. 11 Figure 10: Geometry of the subsolar point Q, relative to that of an observer at P. Only the angle α relative to the center of the earth, and the distance ρS are required to calculate the gravitational forcing or insolation at P. (Adapted from Munk and Cartwright, 1966). ________________ (Some textbooks and papers mistakenly refer to spin precession as “wobble”. In astronomy and geophysics, “wobble” refers to the movement of the spin axis within the earth about the north pole. That is, the earth does not normally spin about the geometric north pole, but at a point nearby which moves within the earth around the nominal north pole with periods around 12 and 14 months, slightly changing latitudes of positions fixed on the Earth and introducing a small “pole tide” via the varying centrifugal force; see Lambeck, 1980.) It is useful to understand a bit about how astronomers specify the position of the sun. First consider Fig. 10 defining the sub-solar co-latitude, θS , and the sub-solar longitude, λS , and the co-latitude, θ, and longitude, λ, of an observer at P (co-latitude is 90◦ −latitude). Fig. 11 shows the solar longitude, lS , measured along the ecliptic plane from Υ. From the point of view of determining insolation at P as a function of time, what matters is the distance ρS and angle α to the observer. The distance to the observer can be expressed as, ³ ´1/2 ρS = RS 1 − 2 (a/RS ) cos α + (a/RS )2 from the law of cosines, in terms of the distance RS between the center of earth and the center of the sun, and the radius of the earth a. Following Munk and Cartwright (1966), define an “insolation function”, that is, the incoming 12 Figure 11: Latitude and longitude of the subsolar point in Greenwich coordinates (λs , θs ). Angles ls , ψs are measured from the vernal equinox Υ. (Adapted from Munk and Cartwright, 1966). ________________ radiation, as R=S0 µ R̄S ρS ¶ cos α R̄S = S0 ¡ ¢1/2 , cos α > 0 RS2 − 2 (aRS ) cos α + a2 = 0, cos α < 0 (1) (2) where α is the angle between the sub-solar point and the observer, S0 is the solar constant and R̄S is the time mean solar distance. Eq. (1) is non-linear in RS and cos α, so that one expects all of the astronomical frequencies governing the relative position of observer and sub-solar point to interact. The abrupt cutoff where there is no illumination is another form of non-linearity.3 One then expands the denominator of (1) as, ! à µ ¶ ¯S R a 2 a R=S0 P1 (cos α) + P2 (cos α) + ... , cos α > 0 1+ RS RS RS = 0, cos α < 0. (3) (4) ¡ ¢ Here the Pn are Legendre polynomials, P1 (cos α) = cos α, P2 (cos α) = 1/2 3 cos2 α − 1 (e.g., Jackson, 1975). 3 This expression is approximate in treating the sun as a point source–ignoring the finite diameter which leads to a slow and not abrupt insolation rise and fall at sunrise and sunset. Similarly, it fails to account for the effects of the atmosphere in generating an extended twilight (see Bowditch, 1962 for a discussion of the definitions of twilight), and which lasts for extended periods near and above the polar circles. 13 From elementary spherical trigonometry in Fig. 11, cos α = cos θS cos θ + sin θS sin θ cos (λS − λ) . (5) which can be substituted into Eqs. 3 and then all one needs is the sub-solar coordinates θS , λS . Smart (1962) or Green (1985) have clear accounts of how these are produced in terms of the angle, M relative to Υ, of a uniformly moving fictitious sun. λS will clearly increase by about 360◦ every day as the observer spins under the sun. 3.2 Non-linear functions of periodic terms Eq. (1) is nonlinear in the angle α, and in the solar distance RS (equivalently, in the observer-sun distance ρS .) α, RS contain near-periodic terms in the cosine (or sine) of a cycle per day, a cycle per year (any year or day definition will do), a cycle per 21,000 years and a cycle per 41,000 years (all periods approximate), plus even longer period terms. Sums of periodic components in non-linear functions have a quite general behavior–they generate products of the periodic terms which in turn lead to “line splitting.” (A “line” is a pure frequency.) Consider the very simple behavior of a rule for insolation depending upon an inverse distance: y (t) = 1/r, r = 1 + δ1 cos (2πs1 t) + δ2 cos (2πs2 t) . Then if δi are both small, ¡ ¢ y (t) = 1 − (δ1 cos (2πs1 t) + δ2 cos (2πs2 t)) + (δ1 cos (2πs1 t) + δ2 cos (2πs2 t))2 + O δi3 (6) ≈ 1 − (δ1 cos (2πs1 t) + δ2 cos (2πs2 t)) + δ12 cos (2πs1 t)2 + δ22 cos (2πs2 t)2 + δ1 δ2 cos (2πs1 t) cos (2πs2 t) = 1 − (δ1 cos (2πs1 t) + δ2 cos (2πs2 t)) + δ12 + δ12 cos (4πs1 t) + δ22 + δ22 cos (4πs1 t) + δ1 δ2 cos (2πs1 t) cos (2πs2 t) Notice that the terms in δi2 generate one term with no time-dependence, and a second term which has twice the underlying frequency. The last term in Eq. (6) can be written as, δ1 δ2 cos (2πs1 t) cos (2πs2 t) = δ1 δ2 δ1 δ2 cos (2π (s1 − s2 ) t) + cos (2π (s1 + s2 ) t) , 2 2 (7) that is, the sum and difference frequencies s1 ±s2 . Variations in amplitude of one nearly periodic process caused by another is called “amplitude modulation” Let s2 be a frequency of one cy cle/year, and let the modulation frequency s1 be one cycle in 21,000 years. Then the amplitude modulated seasonal cycle is split into two neighboring frequencies separated from s2 by being one cycle in 21,000y higher and lower on either side. Unless one has a very long record, a Fourier 14 Figure 12: Example of a sine wave “split” into two nearby frequencies, which alternately constructively and destructively interfere. A fourier analysis of such a record shows no energy on the timescale of the interference change (modulation time scale). But any process sensitive to the amplitude of the high frequency will show such energy. The insolation curves include this kind of effect where there are approximately 20,000 annual cycles during one “beat”–far too many to plot here, and with an irregularity because of the splitting by more than one quasi-periodic term in the precession band. analysis as in Appendix 1, will not be able to distinguish the two different frequencies from a single one at one cycle/year. An analogue of this behavior is shown in Fig. 12. In general almost any nonlinear function, including the insolation function, can be expected to display all combinations of sum and difference frequencies (including self-interactions at zero and twice the frequency) and one must carry out an analysis to find the relative magnitudes. In computing the insolation, one inserts Eq. (5) and uses expressions for λs , θS , as quasi-periodic functions of time and works out the result.4 (See for example, Rubincam, 2004; exactly the same process is used to derive tide disturbing potentials, and the discussion in Doodson, 1921, is helpful.) The function in Eq. (7) has no contribution at or near frequency s1 because s2 >> s1 . Because the effects of precession are proportional to e, they vanish at times when the eccentricity vanishes (the amplitude modulation is of the form 1 + e cos (2πs1 t) with e << 1, so that when e vanishes, the seasonal cycle does not–only its 21,000y band variation). 4 I use the expression “quasi-periodic” for elements of a process exhibiting a clustering of nearby frequencies in a narrow-band, with the understanding that they are deterministic, but irrational multiples of each other so that the process is never strictly periodic one. It may come arbitrarily close to being periodic on a very long time scale. 15 Image removed due to copyright restrictions. Citation: Huybers, P., and W. Curry. "Links between Annual, Milankovitch and Continuum Temperature Variability." Nature 441 (2006): 329-332.�� Figure 13: From Huybers and Curry (2005), showing (lower panel) the logarithm of the squared Fourier series coefficients of insolation at 65◦ N samped every month from one day to about 1 million years periods (plus a small amount of white noise for numerical reasons). Vertical dashed lines are shown at periods of 100,000, 41,000 and 21,000 years. There is energy at all periods, with excess energy in the band around 41,000 years, a slight enhancement at 100,000 years, and not at all near 21,000 years. Diurnal and seasonal peaks are by far the largest. Upper panel shows a synthesis of various observations leading to a rough estimate of the energy as a function of frequency in the climate system. There is energy at all frequencies, a great excess relative to neighboring frequencies at diurnal (not shown) and seasonal and semi-annual time scales, a bit at 41,000 years and in a wide-band around 100,000 years, the latter not much like a narrow peak. The monthly sampling is important to avoid aliasing the annual cycle. The diurnal variation is filtered out by explicitly suppressing terms dependent upon the longitude, M, of the mean sun in the ecliptic plane. 16 The insolation at any fixed point on the Earth will then be a non-linear function involving near-periodic terms at one day, one year (any choice of definition of day or year will suffice), near 21,000 years, near 41,000 years and very weak terms near 100,000 and 400,000 years. One thus expects to see all harmonics of all of these frequencies, and all interaction terms. The coefficients latitudes and longitudes of the observer appear in the insolation equation so that one gets products of trigonometric terms in latitude and longitude as well. Rubincam (2004) has given an explicit expression for the leading terms in the insolation (his Eq. 14) omitting only the terms that have a daily variation. (Averaging out the daily variation is perhaps acceptable for climate purposes, but see the discussion below of rectification.) 3.3 Insolation Variations One needs to calculate the amount of insolation falling on any given point on the Earth as the orientation of that point relative to the sun varies day-by-day, over the course of a year, and then with superimposed longer period variations. Thousands of years of effort have been devoted by astronomers to finding the orientation of observers on the spinning Earth relative to the sun. A sketch of the basic approach is given in Appendix 2. As already said, during the course of the year, θS varies as the sun moves in its tilted orbit (relative to the equator), and observer P rotates under the sun every day. Appendix 2 describes the method by which α and ρs can be determined from purely geometric considerations and Eq. (18) in Appendix 2 reduces to, µ ¸ ¶µ ¶¾ µ ¶2 ½ ∙ R̄s e2 1 e2 5 3 2 1 1+ + P2 (cos θ) 1 + sin εS − , R|non−seasonal = S0 as 4 4 16 2 4 2 (8) so that the only non-seasonal effects enter from the variations of eccentricity e, and obliq uity, εS . aS is the almost-constant semi-major axis of the solar ellipse and S0 is the solar constant. Because e is very small, and appears to remain so over the history of the solar system, the dominant non-seasonal forcing arises from εS at about 41,000 years. If the ef fects of precession are to be retained, one must include the terms containing the annual vari ation. 4 Milankovitch Hypotheses If the climate system were linear, the response to variations in insolation through time would occur only at the quasi-periodicities appearing in R, and Eqs. (18–in Appendix 2, 8) show that these would be daily (although omitted in the expressions above), annually, near 41,000 years and at the frequencies present in e, near 100,000 and 400,000y. One would expect the annual 17 Image removed due to copyright restrictions. Citation: See Figures 4 and 6. Hays, J. D., J. Imbrie, and N. J. Shackleton. "Variations in the Earth's Orbit, Pacemake of the Ice Ages." Science 194 (1976): 1121-1132. Figure 14: From Hays, Imbrie, and Shackleton (1976) showing spectra of proxies in deep-sea cores. Plots on the left are linear in the power, those on the right are logarithmic. cycle to be split (as above) approximately into two frequency bands, of 1 cycle/tropical year ±1 cycle/21,000 years. Fig. 13 shows, however, that all frequencies are present in the climate system, not just those in the insolation forcing. One can distinguish (see Wunsch, 2004) several different Milankovitch hypotheses: 1. Obliquity and precessional band energy is discernible in spectra of climate proxies. 2. Obliquity and precessional band energy dominate climate variability, as measured by their variances, between about 19,000Y and 42,000Y periods 3. Obliquity and/or precessional band energy, irrespective of (2) control (“pace”) the 100,000Y variability characteristic of the glacial-interglacial shifts in the last about 800,000 years. What is the evidence? Imbrie and Imbrie (1986) provide some historical background. The intense modern interest in the hypothesis that orbital changes are important in the climate system can be traced back to an influential paper by Hays et al. (1976). By doing a Fourier (spectral) analysis of a number of deep-sea cores, they produced the results seen in Fig. 14. The spectral analysis is shown in two forms, as both linear power and a logarithmic power plots. Their inference was that energy above the “background” (defined to be those frequencies distintinguishable from the ones dominating the insolation forcing) was very important. Notice that the linear power plot emphasises the peaks, whereas the logarithimic plot suggests the peaks are only marginally above the background. Hays et al (1976) concluded that obliquity and precessional periods forcing controlled the larger energy levels near 100,000 years. (Displacements of the peaks from the expected periods are ascribed to dating errors in the depth versus age conversion. How differing types of dating errors affect different signals has to be worked out in detail; see Wunsch, 18 2000.) In the list of versions of hypotheses above, the Hays et al. results and many that followed strongly support the inference that there is some excess energy in the climate records at the Milankovitch frequencies, keeping in mind the power explained issue in terms of the Parseval relationship in Appendix 1. Whether hypothesis (2) is true is much less clear, as the amounts of energy appearing above the background variability is small. The main issues that have emerged relate to: (A) the presence of the precession band energy and above all, (B) the inference that the 100,000 year-band energy is controlled by the 41,000 and 21,000 year-bands. The hypothesis, however, has led to the assumption that it is true, and then core agemodels determined by forcing the obliquity and precession-band energies in the core to follow fluctuations in the insolation curves–usually at some fixed latitude and time such as 65◦ N on 21 June. A summer day or range of days is chosen on the grounds that summer-time melting will determine whether a continental glacier can grow through time. That is, the insolation curve is “wiggle-matched” as described above, to the proxy record, shifting the maxima and minima generally toward some aspects of the insolation curve with fixed phases between the records being matched. Once records are so adjusted, it is difficult to further explore the relationship between insolation and climate change, as a relationship has then been imposed. (A) Why is there energy in the precession band, when the effect of precession is to split the annual cycle into several nearby frequencies of one cycle/year? A general framework for answering this question lies with radio engineering, where the annual cycle would be the carrier frequency and the precessional variations the amplitude modulation. An AM radio works by detecting the evelope (“demodulating” the carrier, typically by using a “rectifier”). This problem was discussed by Huybers and Wunsch (2003) and we will only summarize their results. First, any component of the climate system responding to forcing in one season only will serve as a demodulator. For example, if summer insolation alone during some particular month controls the growth of an ice sheet, one can rewrite that forcing as, F = R (t) , t ∈ τj , where τj are e.g., July values alone. The “standard” modern Milankovitch hypothesis, repeated in hundreds of papers, is that insolation at 65◦ N typically on specific dates like June 21 (summer solstice), is the controlling forcing. Fourier analysis of such a “gapped” time series will contain energy at the envelope frequency. That is, the gaps act as a rectifier of the envelope frequency. Another possibility is that parts of the climate system, e.g., the Asian monsoon strength, only depend upon summertime forcing, with winter values having no particular influence. Hypotheses such as these two are examples of nonlinearities in the climate system and in general, it is easy to 19 show that almost any nonlinearity (e.g., another example is the assertion that |R (t) − constant| controls climate) will generate energy in the precession band. One can infer that the existence of precession band energy implies that the climate system is non-linear – not a surprise. There is, however, one very serious problem with this straightforward interpretation (Huybers and Wunsch, 2003). Many of the proxies used to infer the climate state (e.g., δ 18 O; dust deposits, etc.) also have a seasonal dependence. So for example, suppose that in summertime when temperatures exceed some threshold, dust is carried by the wind from the Gobi desert into the open Pacific Ocean, but that in wintertime the prevailing winds are in the wrong direction. One would thus see increasing dust deposits in a Pacific deep-sea core as the annual cycle increased, and may be led to infer that the world was warming. But because one would not see a reduced amount of dust from the colder winter (one cannot go below zero dust), no inference about the annual mean can be made. One is entitled to infer that summer temperatures were warmer, but not anything about the annual average. A similar effect arises from the biological processes. If foraminifera grow only in springtime, they do not form an annual average. One can discuss springtime temperatures only. To go beyond that season, one must therefore understand the annual cycle of the climate system–an important complication for someone interested in myriadic (10,000years) time scales and longer fluctuations–as well as the annual cycle of the proxies including, particularly, the biological ones. It is important to keep in mind too, that the precessional effect has the reversed sign across the equator, and so any dependence of the system on the global average forcing vanishes in the precession band. Assuming the driving is restricted to the northern hemisphere is another way of introducing a rectifier. Little attention has been paid to rectification of diurnal insolation variability, presumably because its frequency is so high. But as we have seen, self-interactions give rise to zero-frequency effects, and given the line splitting present of the diurnal frequency, there seems little doubt that rectification of the day-night cycle must occur to some extent (coastal sea-breezes are an example). Above the polar circles, the day/night cycle combines with the seasonal one. Whether such effects are of sufficient amplitude to be of climate concern is unknown. Turn now to hypothesis (B), that obliquity variations and modulation of the annual cycle in the precession band control the 100,000 year glacial cycles of the past 800,000y. Linear climate systems are much easier to analyze and understand, and a number of authors have claimed that the 100,000 year-band was a linear response to the minute change in eccentricity (which controls the global average forcing). There are several reasons why this explanation is very unattractive. If the system is linear, one would anticipate seeing a 400,000 year response where the variation in e has even more energy, and the nominally 100,000 year response period energy 20 Image removed due to copyright restrictions. Citation: Raymo, M. E., and K. Nisancioglu. "The 41 Kyr World: Milankovitch's other Unsolved Mystery." Paleoceanog 18, no. 1 (2003), doi:10.1029/2002PA000791. Figure 15: From Raymo and Nisancioglu (2003) showing a North Atlantic Ocean deep sea core extending back almost 3.5 million years. There is a transition from a variability with excess energy around 40,000 years to one having excess energy at around 100,000yrs. The transition occurs approximately 800,000 years ago. Because the insolation time scales did not change over the whole period, some physical change in the climate system itself is likely responsible. See Wunsch (2004), Tziperman et al. (2006). lies in a frequency band somewhat lower than the dominant change in e. But most important, the hypothesis requires a linear resonance mechanism in the global climate system, a possibility that flies in the face of almost everything known about the behavior of complicated turbulentfluid/cryospheric systems. (See Wunsch, 2000 for a discussion of linear resonances.) How does one obtain a nominally 100,000 year period from forcing at 41,000 and 21,000 years (the latter assuming rectification of the annual cycle)? If one simply assumes that obliquity dominates, then 3 times 41,000 years is about 120,000 years, and 2 times 41,000 years is about 80,000 years. So one hypothesis is that the 100,000 year period is actually an average of glacial cycles occuring at either 121,000 and 80,000 years (this hypothesis was tested by Huybers and Wunsch, 2005, who concluded it was likely correct. Alternatively, 5 times 21,000 years is about 100,000 years; or, instead, an average of 6 times 21,000 and 4 times 21,000 also produces a nominal 100,000 year period (Ridgewell et al., 1999). Both these suppositions, or both in concert, imply that the system is highly non-linear in some way. (Huybers and Wunsch, 2005 concluded that there was significant evidence for the obliquity hypothesis, but that the timing uncertainties in the data were inadequate to say anything about the precession one.) To proceed, one turns to making specific models claiming to describe the climate system in terms of simple physics. This approach is discussed by Tziperman et al. (2006). What they show is that almost any non-linear model can entrain the Milankovitch forcing so that 21 it shows some form of synchrony with insolation. (See Strogatz, 2003 for a general discussion of the phenomenon.) Because the number of non-linear models is essentially infinite, fitting Milankovitch cycles to the last million years of record cannot distinguish between them. At the present time, one must infer that the late Pleistocene ice ages appear to be linked at least in part with insolation variations, but that it is not possible to distinguish the controlling mechanisms (the nonlinear physics) that is acting. Before about -800,000y, the records (see Fig. 15) show little energy in the 100,000y band, but an excess in the obliquity band. Something had to have happened to the climate system about 800,000 years ago, to produce an excess at the lower frequency where there is little direct forcing. No consensus exists on what that change was. Note that the 40,000 year-band energy prior to -800,000 years is only a fraction of the total variability present then, and that much of it persists right through the recent period. (The eye tends to see the dominant pattern wiggles and not much else.) A Stochastic Mechanism Before about -800,000 years, the predominant pattern of variation in Fig. 15 and many other such records, is an approximately 41,000 year time-scale oscillation followed by the recent period (late Pleistocene) in which the variability is dominated instead by an approximately 100,000 year time scale. The 41,000 year time scale is plausibly connected with a direct response to obliquity; that there is no particular energy at 21,000 years would suggest little or no direct dependence upon the seasonal cycle modulation. A large literature has grown up attempting to rationalize the 100,000 year time scale including the nonlinear synchronization process already described. But the system is clearly quite noisy, with variability at all frequencies (see e.g., Fig. 13), and one needs to understand that purely random driving can also produce a predominant time scale. The most basic version (Wunsch, 2004) is simply a random walk in one-dimension. Consider a drunk at one end of a room. He starts at time zero, at position ξ (t = 0) = 0, with his back to one wall and randomly takes a step of length ∆x forward or back at time t = n∆t, where n is an integer and ∆t is a dimensional time-step (but he cannot penetrate the wall at his back; if ξ(n∆t) = 0, a backward step results in no movement). There is another wall at x = L; when the drunk hits x (m∆t) = L, he is flung back to ξ ((m + 1) ∆t) = 0. It is an easy matter to calculate the average time interval between encounters with the wall at x = L. (Feller, 1957 is a classical account of random walks; the subject is closely related to the so-called gambler’s ruin problem.) A time series of the drunk’s position looks visually much like the time series in ice cores (Fig. 16). To make this idea physically plausible, the position becomes ice volume, V (t = n∆t) . When ice volume reaches a threshold, it is postulated that it becomes unstable (there is speculation 22 Image removed due to copyright restrictions. Citation: EPICA Community Members. "Eight Glacial Cycles from an Antarctic Ice Core." Nature 429 (2004): 623-628. Figure 16: (From EPICA Community Members, 2004). Second and third panels show climate proxies extending back approximately 800,000 years (note backwards running time). The records are far from periodic. Note that the insolation curves in the top panel show very different behavior depending upon which elements of insolation are regarded as most important. There is much freedom here to pick and choose best fits to observed records. 23 Image removed due to copyright restrictions. Citation: Wunsch, C. "Quantitative Estimate of the Milankovitch-forced Contribution to Observed Quaternary Climate Change." Quaternary Sci Rev 23 (2004): 9-10, 1001-1012. Figure 17: Synthetic time series generated from a simple random walk and a threshold collapse requirement. All units are arbitrary here. (Wunsch, 2004). about such behavior). The non-negative value of ξ corresponds to V = 0, no negative ice volumes. Fig. 17 shows a simulation of such a simple random process (Wunsch, 2004). At the present time, it is very difficult to determine if this process could be dominant in the climate system. A Summary There is convincing evidence that excess (relative to the background) 41,000 year-band energy exists in the spectrum of climate variability. In many, but not all, cores there is some evidence of 21,000 year-band energy, and whose significance is ambiguous because of the seasonality in many proxies. The amount of energy in these Milankovitch bands is a small fraction of the record variance. Nonetheless, one can construct simple nonlinear, conceptual, climate models that lead to oscillations with an average time scale near 100,000 years determined by the insolation variations. One can also construct purely random models that do so as well. Presumably, both deterministic and stochastic driving are present, and it seems likely that both mechanisms are at work. The observational record is insufficient to distinguish between the multiplicity of nonlinear mechanisms that could be operative. Some evidence (Huybers and Wunsch, 2005) exists for control by obliquity; control by precession cannot be ruled out, but there is no positive evidence in its favor. A major remaining enigma concerns the transition from 40,000 year time-scale to 100,000 year time-scale dominance, with rationalizations including a major change in CO2 content of the atmosphere and the expected behavior of an intrinsically chaotic system. (See e.g., Huybers and Tziperman, 2006). 24 5 Appendix 1. Fourier Representations Most scientists encounter Fourier analysis as undergraduates studying periodic functions. So, for example, let x (t) be periodic with period, T. Then one learns that x (t) can be written exactly as a Fourier series, µ µ ¶ X ¶ ∞ ∞ 2πnt 2πnt a0 X an cos + bn sin , 0≤t≤∞ x (t) = + 2 T T n=1 n=1 µ ¶ µ ¶ Z Z 1 T 2πnt 1 T 2πnt an = x (t) cos dt, bn = x (t) sin dt. T 0 T T 0 T (9) But, and this is an essential conceptual point, one can represent any record exactly as in Eq. (9) merely by restricting 0 ≤ t ≤ T, that is defining T as the period of observation, making no statement about how the record would behave outside the interval of measurement. There is then no implication of periodicity, but the representation remains exact within the restricted interval, and the values of an , bn can still convey a great deal of information. (Of course, there could be truly periodic elements present, but their existence as such would need to be separately demonstrated (by a variety of means)). It is an easy matter to prove the Parseval-Rayleigh theorem: 1 T Z T x (t)2 dt = 0 ∞ X ¡ 2 ¢ an + b2n (10) n=1 (I have here assumed the mean, a0 /2 = 0, which can always be arranged; it saves a bit of writing and removal of the sample mean should always be done prior to any practical Fourier analysis.) The left-hand-side of Eq. (10) is commonly used as a measure of the “energy” in a record (the mean-square value). The right-hand-side shows that the energy is the sum of the squares of the Fourier series coefficients. Thus the explanatory power, PE , of any group of Fourier coefficients of a record can be measured as the fraction ¢ ¢ Pp2 ¡ 2 Pp2 ¡ 2 2 2 p=p1 ap + bp p=p1 ap + bp = P∞ PE = R 2 2 2 1 T x (t) dt n=1 (an + bn ) T 0 In many records, a2n + b2n clusters around specific values of n, showing that the power is concentrated in particular frequency bands (see Fig. 18). Within that band several different physical behaviors can be present. (1) There is an isolated value a2p + b2p >> a2n + b2n where p represents a true pure frequency (a “line spectrum”). Such physical processes are extremely rare in nature. The closest examples are tidal records, and other astronomically-driven phenomena. But even then, isolated pure frequencies are the exception rather than the rule. Instead one has, 25 anton I. Sea Level Power Density Estimate 17−Nov−2006 12:29:22 CW−lptp 6 10 12 hour band 4−days mm /CPH 4−days 5 5 10 10 4 10 0 3 10 0 10 10 −3 10 −2 10 CYCLES/HOUR CYCLES/HOUR lunar fortnightly tide mm /CPH 12 hour band 5 5 10 10 0 0 10 −1 10 10 0.07 0.072 CYCLES/HOUR 0.074 0.076 CYCLES/HOUR 0.078 Figure 18: Upper left panel shows a power density spectral estimate from a sea level record at Canton Island (almost on the equator in the central Pacific Ocean). The Fourier coefficients used in this estimate have a variety of characters. Note that about 95% of the variance, as measured by the Parseval relationship, lies in the semi-diurnal tide (labelled as 12 hour band). Upper right panel shows an expanded view of the power density estimate in the vicinity of the peak near 4 days (discussed at length by Wunsch and Gill, 1976). Here the Fourier coefficients are essentially random. Lower left panel shows the Fourier coefficients in the vicinity of 12 hours where there are a series of nearby almost pure tidal frequencies (pure sinusoids) arising from the solar and lunar tides. Lower right shows the spectral density near a two-week period (note the linear frequency scale) where there are a series of comparatively weak lunar tides (Mf–the lunar fortnightly tides) also a cluster of periodic terms. In all frequency bands there is a background of fourier coefficients best understood as arising from a random processes on which these other processes are superimposed. For a Gaussian random process the squared fourier coefficients are a χ2 random process. (For display purposes, the raw a2n + b2n have been averaged over immediate neighbors to render the result somewhat smoother in frequency than would be the case otherwise.) The red line interval denotes an approximate 95% confidence interval for those Fourier estimates that are random variables. 26 0.08 (2) A process represented by Eq. (7) giving rise to a “fine structure” or, generalizing, µ ¶ δ1 δ2 δ1 δ2 cos (2π (s1 − s2 ) t) + cos (2π (s1 + s2 ) t) cos (2πs3 t) δ1 δ2 cos (2πs1 t) cos (2πs3 t) cos (2πs2 t) = 2 2 (11) = δ1 δ2 cos (2π (s1 − s2 + s3 )) + .... 4 (12) representing a further splitting into a “hyper-fine” structure, assuming s1 , s3 << s2 , etc. This situation describes the tidal frequencies (see Munk and Cartwright, 1966) as well as the Mi lankovitch forcing. The ability to determine whether an apparent line is actually split depends directly upon the record length, (frequencies are separable in the Fourier analysis by 1/T ). In addition to this amplitude modulation, one can modulate the frequency, with similar results. Such components of records are deterministic, but not necessarily periodic if the underlying real frequencies are generally not commensurable (not integer multiples of each other), the result is nearly, but not actually periodic. Consider by way of example, the Fourier series over a finite interval of a function whose exact form is, µ ¶ µ ¶ 2πt 2πt + α2 cos , x (t) = α1 cos T1 T2 where T1 ≈ T2 , but T1 /T2 cannot be written as the ratio of two integers (it is irrational).) Such almost periodic terms appear in the demodulated insolation forcing at periods around 41,000 and 21,000 years. (3) x (t) has random components. That implies the an , bn are also random at least in part. Interpretation of their particular values leads to what is known as spectral analysis (Percival and Walden, 1993 is a good textbook; there are many others). They may well cluster about particular frequencies p/T , yet there is no implication that the physics are periodic. Such phenomena are very common in nature (see Fig. 18). One way of distinguishing such “narrow band random processes” from those that are truly or nearly periodic as in Eqs. (7, 11) involves understanding the relationships of the phases of nearby Fourier coefficients, in the form tan−1 (bp /ap ) among others. With modern computers, t is almost always taken to be discrete, and all integrals are replaced by sums. The zero-order message is that almost all non-periodic records can be written as Fourier series over finite time intervals; there is no implication of periodicity, and the values of the Fourier coefficients, which almost never vanish, contain a great deal of useful information. For the orbital motions described above, the an , bn are not random (until records exceed many tens of millions of years and motions become chaotic). Aliasing. 27 Aliasing refers to the effect of not sampling a record rapidly enough, so that high frequencies appear fictitiously as low frequencies. Many textbooks discuss this insidious phenomenon. It is widely ignored, but is extremely important and is the source of some published blunders. As an example consider a climate record having an annual cycle at precisely the tropical year of 365.2422 civil days and which is sampled nominally on 1 July every year, that is at the interval of 1 Julian year of 365.25 days. The true record is, ¶ µ 2πt , y (t) = cos 365.2422 and the sampled record is, yn = y (n∆t) = cos µ 2π (365.25)n) 365.2422 ¶ n = 0, 1, ..., N. (Try plotting it for N = 100000.) It is not hard to show (e.g., Percival and Walden, 1993), that the apparent periodicity is near 47,000 years and the apparent period is extremely sensitive to the numerical definition of the sampling year. It might appear that sampling at 365 day intervals rather than 365.25 days would be a very small difference, but the apparent period when the sampling is at ∆t = 365 days is very different. The expression for the aliasing frequency, sa , is that one must determine n such that ¯ 1 n ¯¯ ¯ . sa = ¯s − ¯≤ ∆t 2∆t Frequency Content. ­ ® Two records whose average Fourier coefficients, a2n + b2n as a function of frequency, n/T, are similar must tend, on average, to have the same number of zero crossings, or of maxima and minima in any finite interval. Even if they are in truth completely unrelated, the human eye tends to align the peaks and troughs, permitting “tuning” if there is any age-model uncertainty to render records even more alike (e.g., Fig. 5). Because of this tendency for similarity to be perceived, statisticians have developed objective techniques, such as coherence analysis–e.g. Priestley, 1982– to measure the true similarity of records. 6 Appendix 2 In Fig. 11, the solar longitude, lS , measured along the ecliptic from Υ, is not an angle which changes uniformly with time, because the angular movement of a body in an elliptic orbit is not uniform. Consider Fig. 19 (adapted from Smart, 1962, p. 111), The position of the sun in the ellipse is given by its radius, RS , and the angle v, called the “true anomaly”. Kepler’s Second Law says that in an elliptical orbit, “equal areas are swept out in equal time”, but this implies that 28 Figure 19: Geometry of a body at S moving in an exaggerated ellipse about the earth at a focus. ν is the true anomaly of the body (the sun), E measured relative to the center of the ellipse is the eccentric anomaly. ω defines the angle between the vernal equinox Υ and perihelion. aS , b are the semi-major and minor axes of the ellipse, and aS is also the radius of the bounding circle. ν, E, M (the solar anomalies) are all measured from perigee. M is the angle of the position of the fictitious mean sun moving uniformly along the circle. (Adopted from Smart, 1962.) the rate of change of v varies with time, and is inconvenient to compute. Astronomers prefer to use instead the angle E, measured from the ellipse center, called the “eccentric anomaly” which can then in turn, be related to another angle, called the “mean anomaly” which has a uniform rate of change. Then with some elementary trigonometry and ellipse geometry, Smart (1962) shows that Furthermore, the radius vector, (1 + e)1/2 v tan tan = 2 (1 − e)1/2 µ E 2 ¶ . RS = aS (1 − e cos E) , (13) (14) where aS is the semi-major axis of the ellipse. Define a fictitious mean sun moving at the average angular rate of the real sun at orbital longitude M , the mean anomaly, M = ω (t − τ ) (15) where τ is a time origin and ω/2π is a frequency defined as one cycle/year. Some more elementary geometry and Kepler’s Law shows that the relationship between E and M is given by the “Kepler Equation”, E − e sin E = M = ω (t − τ ) . 29 (16) For small e (typically about .02 for the sun), the Kepler Equation can be solved iteratively for E in terms of t, or M ; then v is known from (13) , and RS from (14) . Ignoring all of the details, it is apparent that v and RS will carry periodicities at one year, 6 months, and all of the higher harmonics of the year Astronomical references give the fundamental parameters as functions of time, typically in Taylor series in terms of Julian Centuries, T , of 36525 days. So for example, Anonymous (1961, p. 98), gives the mean anomaly, M = 358◦ 280 33.0000 + 12959657900 .10T − 0.5400 T 2 + ... (17) relative to the vernal equinox of 1 January 1950. Note that the mean anomaly is not exactly a uniformly increasing angle with time. The Milankovitch theory commonly uses the position of Υ in 1850, but the geometry can be developed for moving positions and Eq. (17) is not adequate for paleoclimate work. Starting with an expression identical to (3), Rubincam (1994) omits all terms dependent upon λ, thus suppressing all of the daily and higher frequency variations, leaving only the long-period terms. He then finds, to lowest order (following a correction, Rubincam, 2004), ⎧ h i 1 e2 5 2 ⎪ P (cos θ) 1 + + 2e cos M + e cos 2M + ⎪ 0 4 2 ⎪ ⎪ ⎤ ⎡ 4³ ´ ⎪ 2 ⎪ e ⎪ ⎪ sin (ω + M ) + 2e sin ε sin (ω + 2M ) + sin ε 1 − S S ⎪ 2 1 ⎪ ⎦+ ⎣ ⎪ 2 P1 (cos θ) ⎪ e2 27 2 ⎪ sin ε sin (ω − M ) + e sin ε sin (ω + 3M ) µ ¶2 ⎪ S S ⎨ 8 ´ 8 ⎡ ³ R̄S R ¡3 2 ¢ 3¡ ¢ 2 e2 1 7 2 = − 1 − sin εS cos (2 (ω + M )) sin ε − e 1 + S 2 4 2 4 2 ⎪ S as ⎪ ⎢ ⎪ ¢ ¡ ⎪ ⎢ ⎪ +2e 34 sin2 εS − 21 cos M + 34 e sin2 εS cos (2ω + M ) ⎪ ⎢ 5 ⎪ P (cos θ) + ⎪ ⎢ 2 ¡ ¢ 16 ⎪ 9 5 2 3 1 2 2 ⎪ ⎢ ⎪ e sin ε cos (2ω + 3M ) + e sin ε − − S S ⎪ ⎣ 4 2 4 2 cos 2M ⎪ ⎪ ⎩ 39 2 2 − e sin ε cos (2ω + 4M ) + ... 8 S ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (18) If averaged over a year (all terms in M , which has a one year periodicity) disappear, one obtains Eq. (8). Here, M is again the true anomaly, ω the angle of perigee, as is the semi-major axis of the solar orbit. 30 References Alley, R. B., 2005: Abrupt climate changes, oceans, ice, and us. Oceanography, 17, 194-206. Anonymous, 1961: Explanatory Supplement to The Astronomical Ephemeris and The Amer ican Ephemeris and Nautical Almanac. Her Majesty’s Stationery Office, London, 533pp. Bowditch, N., 1962: American Practical Navigator. United States Government Printing Office, Washington, 1524 pp. Bradley, R. S., 1999: Paleoclimatology, 2nd Ed. Academic, San Diego, 610 pp. Broecker, W.„ 2000: Abrupt climate change: causal constraints provided by the paleoclimate record. Earth Sci. Revs., 51, 137-154. Brooks, C. E. P., 1923: Variations in the levels of the central African lakes Victoria and Albert. Geophys. Mem. London, 2, 337-344. Cronin, T. M., 1999: Principles of Paleoclimatology. Columbia Un. Press, New York, 560 pp. Doodson, A. T., 1921: The harmonic development of the tide-generating potential. Proc. Roy. Soc., A, 100, 305-329. EPICA Community Members, 2004: Eight glacial cycles from an Antarctic ice core. Nature, 429, 623-628. Feller, W., 1957: An Introduction to Probability Theory and Its Applications, Second Ed. Wiley, New York, 461pp. Gauch, H. G., 2003: Scientific Method in Practice. Cambridge Un. Press, 435 pp. Green, R. M., 1985: Spherical Astronomy. Cambridge Un. Press, Cambridge, 520pp. Hays, J. D., J. Imbrie, and N. J. Shackleton, 1976: Variations in the Earth’s orbit, pacemake of the ice ages. Science, 194, 1121-1132. Hendy, I. L., J. P. Kennett, E. B. Roark, and B. L. Ingram, 2002: Apparent synchorneity of submillennial scale climate events between Greenland and Santa Barbara Basin, California from 20-10ka. Quat. Sci. Rev., 21, 1167-1184. Huybers, P. and W. Curry, 2006: Links between annual, Milankovitch and continuum tem perature variability. Nature, 441, 329-332. Huybers, P. and E. Tziperman, 2006: Integrated summer insolation forcing and 40,000 year glacial cycles from an icesheet/energy-balance model perspective. unpublished document, . Huybers, P. and C. Wunsch, 2005: Obliquity pacing of the glacial cycles. Nature, 414, 491-494. Huybers, P. and C. Wunsch, 2003: Rectification and prcession signals in the climate system. Geophys. Res. Letts., 30(19), 2011, doi:10.1029/2003GL017875, . 31 Huybers, P., 2006: Early Pleistocene glacial cycles and the integrated summer insolation forcing. Science, 313, 508-511. Imbrie, J. and K, P. Imbrie, 1986: Ice Ages. Solving the Mystery. Harvard Un. Press, Cambridge, 224 pp. Jackson, D. D., 1975: Classical Electrodynamics, 2nd ed. John Wiley, New York, 848 pp. Kahneman, D., P. Slovic, A. Tversky, 1982: Judgment Under Uncertainty: Heuristics and Biases. Cambridge Un. Press, Cambridge, 555pp. Lambeck, K., 1980: The Earth’s Variable Rotation: Geophysical Causes and Consequences. Cambridge University Press, Cambridge, 449pp. McDermott, F., D. P. Mattey, C. Hawkesworth, 2001: Centennial-scale Holocene climate variability revealed by a high-resolution speleotherm δ 18 O record from SW Ireland. Science, 294, 1328-1331. Munk, W. H. and D. E. Cartwright, 1966: Tidal spectroscopy and prediction. Phil. Trans. Roy. Soc., A, 259, 533-581. Percival, D. B. and A. T. Walden, 1993: Spectral Analysis for Physical Applications. Multitaper and Conventional Univariate Techniques. Cambridge Un. Press, Cambridge, 583 pp. Priestley, M. B., 1982: Spectral Analysis and Time Series. Volume 1: Univariate Series. Volume 2: Multivariate Series, Prediction and Control, 890pp. plus appendices (combined edition), Academic, London. Raymo, M. E. and K. Nisancioglu, 2003: The 41 Kyr world: Milankovitch’s other unsolved mystery. Paleoceanog., 18(1), doi:10.1029/2002PA000791. Ridgwell, A. J. A. J. Watson and M. E. Raymo, 1999: Is the spectral signature of the 100 kyr glacial cycle consistent with a Milankovitch origin? Paleoceanog., 14, 437-440. Roe, G. H. and M. R. Allen, 1999: A comparison of competing explanations for the 100,000 yr ice age cycle. Geophys. Res. Letts., 26, 2259-2262. Rubincam, D. P., 1994: Insolation in terms of Earth’s orbital parameters. Theor. Appl. Climatol., 48, 195-202. Rubincam, D. P., 2004: Black body temperature, orbital elements, the Milankovitch preces sion index, and the Seversmith psychroterms. Theor. Appl. Clim., 79, 111-131. Siedelmann, P. K., Ed., 2006: Explanatory Supplement to the Astronomical Almanac. Un. Science Books, Mill Valley, CA, 752pp. Smart, W. M., 1962: Textbook on Spherical Astronomy, Fifth Ed. Cambridge Un. Press, Cambridge, 430 pp. Strogatz, S. H., 2003: Sync: The Emerging Science of Spontaneous Order. Theia, New York, 338pp. 32 Vanmarcke, E., 1983: Random Fields: Analysis and Synthesis. The MIT Press, Cambridge, 382 pp. Vimeux, F., K. M. Cuffey, and J. Jouzel, 2002: New insights into Souther Hemisphere temperature changes from Vostok ice cores using deuterium excess correction. Earth Plan. Sci. Letts., 203, 829-843. Winograd, I. J., T. B. Coplen, J. M. Landwehr, A. C. Riggs, K. R. Ludwig, B. J. Szabo, P. T. Kolesar, K. M. Revesz, 1992: Continuous 500,000-year climate record from vein calcite in Devils Hole, Nevada. Science, 258, 255-260. Wunsch, C. and A. E. Gill, 1976: Observations of equatorially trapped waves in Pacific sea level variations. Deep-Sea Res., 23, 371-390. Wunsch, C., 1999: The interpretation of short climate records, with comments on the North Atlantic and Southern Oscillations. Bull. Am. Met. Soc., 80, 245-255. Wunsch, C., 2000: On sharp spectral lines in the climate record and the millennial peak. Paleoceanog., 15, 417-424. Wunsch, C. , 2004: Quantitative estimate of the Milankovitch-forced contribution to observed quaternary climate change. Quaternary Sci. Rev., 23/9-10, 1001-1012. Wunsch, C., 2006: Abrupt climate change: an alternative view. Quat. Res., 65, 191-203. 33