Paleoceanography - Glacial Interglacial Oscillations

1 PALEOCEANOGRAPHY Lecture 6 GLACIAL/INTERGLACIAL 'OSCILLATIONS': WHY? Time-series analysis We now understand something of the course of climate change over the past 700,000 years. Similar basic patterns are seen in most indicators, although we can recognize that geographical variation is significant. We can begin to ask the question: Why? Two approaches: 1. "The physics": derive ice ages from first principles. Good luck! (call me when you get there). 2. "Correlation (not causation)": look for similarities (coincidences?) between certain driving forces known from first principles and the climate record. Examine the odds that such a correspondence is due to chance. If the odds are good, look into the relevant mechanism in more detail. To be successful, this approach needs: (1) reliable climate record not excessively modified by bioturbation; (2) but the record also needs to be reasonably long so that many cycles are available for statistical tests. Must compromise here! Eventually need to come to an accommodation with “the physics”. I. Historically: A. 1840: Agassiz proposed massive continental glaciation; debate ensued which eventually was decided in his favor B. 1860: Croll proposed that changes in the earth's orbital parameters were responsible for glaciation. Theory received a mixed reception. C. 1920 Milankovitch undertook detailed calculations which quantified the expected variations in the earth's orbital parameters. Theory still encountered a mixed reception. D. 1950's: Emiliani found stronger evidence for cyclic oscillations; tried to revive Milankovitch. Problem: time scale. E. 1960's, 70's: Barbados data (and hence correct time scale) revived interest in Milankovitch. Theory was being taken seriously, but it was considered far from proven. F. 1976: Hays, Imbrie, Shackleton paper. Overcame most of the resistance to some degree of orbital influence over climate (at least as a "pacemaker" of ice ages). 3 G. At present, relatively few people doubt that orbital changes play a role in climate change, but there are questions as to their importance relative to other factors. To what extent is climate actually predictable? Are there significant internal resonances that interact with orbital forcing? What is the origin of large amplitude sub-orbital climate variability? Can we parlay our understanding of past climate change into better predictions of future climate change? Some commentary on the last question: Climate models (to be discussed more later) incorporate a lot of first principle physics (e.g. laws of motion and thermodynamics) but because of the enormous complexity of the system and limited computer size and speed, contain less definitive representations of sub-grid scale processes and poorly understood processes (e.g. convective precipitation; sea ice formation and melting; relation between soil, vegetation, evaporation, precipitation and river runoff). In the aggregate, these factors are extremely important, so they cannot be left out. Instead they are handled by empirical parameterizations containing "fudge factors" which are "tuned" to give a reasonable representation of modern climate. This approach is OK for weather forecasting, because many of these factors can't change fast enough to effect tomorrow's weather. But for climate change, these uncertainties are a serious problem. The sensitivity of the models to various forcing factors might be seriously mis-estimated by these parameterizations. One way to test the sensitivity of these models is to see how they behave when driven by well-known forcing functions (e.g. a recent volcanic eruption). One of the best known forcing functions linked to large climate change is the effect of the earth's orbital parameters on incoming radiation. II. Changes in the earth's orbital parameters and their influence on radiation receipt at the top of the atmosphere. Root cause: gravitational interaction of earth with sun and planet. Precession due to the gravitational field of the Sun, and tilt/eccentricity variations due to Sun-Jupiter tug-of-war for the earth. A. The seasons result from the tilt of the earth's axis relative to the plane of the orbit (view from plane of orbit: Edge plane view (circular orbit): no r th ern sum me r Oblique view, with exaggerated eccentricity no r th ern w int er 4 B. Eccentricity = e = (a2 - b2)1/2 / a b a 1. For the earth, e varies between 0 and 0.06 quasi-periodically, with period of ~96,000 years and a longer (and higher-amplitude) period of ~400,000 years. 2. The sun is at one of the foci of the ellipse. The earth moves faster when it is approaching closer to the sun, and slower as it moves away. 3. Looking down at the solar system from outer space (above the plane of the planetary orbits), the ellipse precesses in its orientation: time 2 time 1 eccentricity vs time: 5 0.06 eccentricity 0.05 0.04 0.03 0.02 0.01 0 0 100 200 300 400 500 600 700 800 900 1000 Age, ka bp 2. Change in eccentricity doesn't affect globally averaged net annual solar insolation much, however (~0.1%) - losses in the summer are balanced by gains in the winter. C. Obliquity of the ecliptic (ε) (tilt of the earth axis relative to the plane of the orbit). ε 6 1. Obliquity varies between ~21.8° and 24.4°; at present the obliquity is 23.44°. Affects the latitude of the tropics and the arctic circle. Period is ~41,000 years, and is relatively regular. Think of this as part of a solar system angular momentum oscillator (angular momentum conserved in solar system, but with transfers between objects). obliquity vs. time: obliquity, degrees 24.50 24.00 23.50 23.00 22.50 22.00 0 50 100 150 200 250 300 350 400 450 500 Age, ka bp D. Precession of equinoxes. Physics is similar to a wobbling top; the precession frequency depends on angular momentum of the earth and the torque imposed by the sun on a nonspherical earth. 1. Precession of the earth's axis relative to the fixed stars has a period of 25,700 years. The climatic affect (due to changes in the earth-sun distance as a function of season) depends on an interaction between the precession period of the earth's axis with the precession period of the major axis of the orbit). In ~22,000 years, the earth ends up in a configuration where the radiation distribution cycle is complete: 7 Initial state: Northern Summer Northern Winter ~22,000 years later: Northern Winter Northern Summer Earth’s axis has precessed 22/26 of a cycle, and the eccentricity of the orbit has precessed “backwards” 4/26 of a cycle, bringing us back to a configuration which is equivalent in terms of radiation 8 The net result is an effective ~22,000 year radiation cycle. The effect of precession on radiation depends on the absolute eccentricity of orbit. When the earth has a nearly circular orbit around the sun, precession of the axis produces no net radiation change for either season in either hemisphere. When the earth's orbit is more elliptical, there is a strong difference between radiation received between the hemispheres and between the seasons: southern summer e sin ω is a simple measure of this radiation cycle. 9 Eccentricity "amplitude modulates" the precession radiation cycle. E. These changes in the earth's orbital elements affect radiation received at the top of the atmosphere ("insolation") as a function of hemisphere, latitude, and time of year. 1. Effect of precession/eccentricity: a. Inter-hemispheric radiation receipt: the hemisphere whose summer occurs during perihelion receives more insolation than the other hemisphere). The higher the eccentricity, the stronger the contrast. b. Seasonal contrast within a hemisphere: as a hemisphere alternates between periods of summer perihelion and winter perihelion, the seasonal contrast in insolation is modified. The two hemispheres are out of phase with respect to this effect. 2. Effect of obliquity: high-latitude summer insolation: high tilt leads to more radiation received at high latitudes during summer (that comes at the expense of the mid- and lowlatitudes; high-latitude winters not affected by changes in tilt). F. These orbital characteristics can all be calculated accurately via Newtonian mechanics into the distant past (many millions of years). Accumulation of uncertainties (e.g. in relative masses of planets; shape and density distribution of planets) causes absolute values (i.e., which side of the sun is the earth on at some exact point in time in the future?) to be uncertain beyond several million years, although the frequencies are probably known reasonably well beyond that time. For the last million years, we may regard these calculations as virtually unquestionable. III. How do we make the comparison between these precisely calculable orbital changes and the record of climate change? A. Curve-matching is always popular. But given subjectivity, the possibility of phase lags, and a rubber time scale, this approach has drawbacks. B. Objective methods: many are possible, but for cyclic processes, Fourier analysis is a natural. It is nonetheless possible to be fool oneself with it's mathematical elegance, however; the statistics assume certain characteristics that might not apply to the data at hand. IV. FOURIER POWER SERIES ANALYSIS A. Prism analogy: look at component frequencies of a time-variant process. Fourier showed that any reasonable function of time can be represented as a sum of sines & cosines. B. Theorem: Let G(t) be any well-behaved (continuous; differentiable; continuous derivatives) function. If we sample this function discretely in uniform time intervals (creating a series G(t1), G(t2), ...G(tn), then there is one (and only one) complex function S(f) = A(f) + B(f)i [where s is frequency and i is the square root of -1] such that G(t) = Σ (aj sin fj + bj cos fj) where fj = j/tn for j = 1 ... n These harmonics fj are orthogonal (i.e. one cannot reconstruct fj wave as Σfk≠j and complete (i.e. any "reasonable" signal in t can be fit perfectly & uniquely). 1. Formal definition of fourier transform of the infinite time series G(t) is +∞ S( f ) = ∫ G(t)e −i ωt dt −∞ where ω = 2πf a. Complex short-hand: eiy = cos y + i sin y 0.8sin(x) 1 0.6cos(x) 0.8 0.8 sin(x) + 0.6 cos(x) 0.6 0.4 0.2 0 -0.2 0 2 4 6 8 10 12 14 -0.4 -0.6 -0.8 -1 2. The Fourier transform is reversible: +∞ G(t) = ∫ S( f )eiω f df −∞ (i.e. both the Fourier transform and the time series contain equivalent "information", which also means that we can do frequency-filtering using the inverse fourier transform) b. The Fourier Transform of a single sine wave is a spike at a single frequency. Similarly, the Inverse Fourier Transform of a single spike frequency is a sine wave. C. Computation: 1. Accepting the truth of the theorem, the discrete Fourier Transform is just an application of simultaneous equations. For discrete uniformly-sampled data, we are taking N observations of G(t), and converting to N/2 pairs of sin, cos coefficients. The number of coefficients is equal to the number of observations, so this is an exercise in solving a system of N simultaneous equations. Let A = an N x N square array of the sin and cos time series for each of the harmonics, where each row of the array (for i= 1 to N) is: sin(2πi/N) cos(2πi/N) sin(2.2πi/N) sin(2.2πi/N) … sin(N/2.2πi/N) cos(N/2.2πi/N) e.g.: sin(2π/N) cos(2π/N) sin(4π/N) cos(4π/N) sin(sin(6π/N) cos(6π/N) … sin(4π/N) cos(4π/N) sin(8π/N) cos(8π/N) sin(sin(12π/N) cos(12π/N) … … and let x = an N x 1 vector of the Fourier coefficients (a, b) arranged vertically: a1 b1 a2 b2 … aN/2 bN/2 and let g = an N x 1 vector of the discrete observations of G(t): g(t1) g(t2) … g(tN) Because the Fourier Transform theorem says: (a) g = A x then the Fourier coefficients are given by the simultaneous equations matrix solution: (b) x = A-1 g (You might wonder if the rows must be linearly independent, ensuring that the inverse of X exists. In effect, the Fourier theorem is a statement that they must be, for "reasonable" time series functions) 2. The above method is one way to calculate the Fourier transform. However, the fastest, most efficient route is the Fast Fourier Transform (FFT). This is a simple routine which calculates the solution to the series of simultaneous equations in a highly efficient way. It occupies about 64 lines of a BASIC program; see "The Visible FFT" spreadsheet to understand the algorithm). Typically, you either interpolate your data to give a series with a power of two of evenly spaced (in time) data points (which usually makes sense because paleoclimatic data sets are usually unevenly spaced in time), or add "trailing zeroes" to your data set until you reach a power of two. 3. Older papers use the method of "Blackman and Tukey" which was developed before the FFT. It was commonly used throughout the SPECMAP project. It is based on the "autocorrelation" technique, which involves calculating the correlation coefficient of a time series relative to itself (autocorrelation) after sliding the series by a variable number of "lags". The methods are mathematically equivalent, however, and the results should not depend strongly on the technique used. 14 4. However, expressing the method as in equation IV.C.1 points the way to alternative ways of calculating the spectra in situations involving unevenly spaced data and/or time series with gaps. In those cases, one can convert the time series of IV.C.1. to a time series with sin, cos (X) values corresponding to "times" (t) of the individual measurements, regardless of their spacing. So the N rows in matrix A become sin(2π(t1- t1)/P) sin(2π(t2- t1)/P) cos(2π(t1- t1)/P) cos(2π(t2- t1)/P) sin(2*2π(t1- t1)/P) sin(2*2π(t2- t1)/P) cos(2*2π(t1- t1)/P)…sin(N*2π(t1- t1)/P) cos(2*2π(t2- t1)/P)…sin(N*2π(t2- t1)/P) cos(N*2π(t1- t1)/P) cos(N*2π(t2- t1)/P) cos(2π(tN- t1)/P) sin(2*2π(tN- t1)/P) cos(2*2π(tN- t1)/P)…sin(N*2π(tN- t1)/P) cos(N*2π(tN- t1)/P) …… sin(2π(tN- t1)/P) where ti is the time corresponding to each observation (t1 is the first observation) and P is the time length of the series of observations (tN-t1). This format can result in a situation where there are more Fourier coefficients (which must remain as the harmonics of the time series record length) than observations. This situation can be dealt with by ignoring the highest harmonics (i.e., chopping off the right hand side of IV.C.4 as necessary). In fact, it is possible to ignore most of the higher harmonics entirely [which often makes sense in paleoclimatic time series, where the error on the time scale renders the highest Fourier coefficients meaningless), resulting in a situation where there are more observations than Fourier coefficients. In this case, equation IV.F.1.(a)??? can be solved using a least squares approach: x = (ATA)-1 AT g These methods aren't as computationally efficient as the FFT, but given the power of modern computers, the drain on computing power is trivial when relatively few spectra are being calculated. 5. The finite discretely sampled time series Fourier transform treats the data as though it were periodic, i.e. the measured time series repeats itself over and over again. This leads to some problems: if first and last data points are quite different, the FT must use many highfrequency harmonics of the fundamental frequency in order to create the sharp break. C. Nyquist folding frequency and aliasing: with a discrete sampling frequency, it is not possible to observe frequencies that are higher than half that sampling frequency. Even worse, the higher frequencies come back at you as if they were signals of lower frequency (aliasing): This flaw can be avoided by using a filter to eliminate frequencies that are higher than the Nyquist frequency (such as sampling at a very high rate and integrating the signal for the duration of the sampling interval). Bioturbation is actually helpful here - you won't find daily cycles in oxic marine sediment cores. 15 response D. Harmonics: the spectrum of a sawtooth (repeated linear trend) includes harmonics of the fundamental period P: time 30 Relative Power 25 20 15 10 5 0 P P/2 P/3 P/4 P/5 P/6 P/7 P/8 P/9 P/10 E. Because we stack finite-length time series end-to-end to obtain an effectively "infinite" time series, we must take steps to avoid spurious harmonics due to record-length trends and sudden offsets between the beginning and end of a time series: 16 1. De-trending: remove the long-term trend from the data (by a linear regression (or any other appropriate trend, which need not be linear). time 2. Tapering: in order to ensure that there are no abrupt mismatches between first and last data points (e.g. a glacial connecting to an interglacial with no transition), the time series can be tapered to zero at the beginning and end. For example, a "cosine bell" can be used for this purpose (a cosine maximum to minimum multiplier) for the first and last 10% of the data series. The subject of tapering is complicated, and at present it is generally agreed that the "optimal" method is David Thompson's "multitaper method", which blends the result from a family of tapers. time 3. These methods have some undesirable characteristics. In some cases, it may be preferable to avoid one or the other of these methods. For example, for a low-frequency sine wave, detrending will tend to "pick out" the endpoints of the time series and could produce a "tilted" corrected sine wave. Best to apply thoughtfully while visually examining the effect the procedure has upon on the time series. F. Then: compute power spectrum 1. Electrical analogy: Power = i2 R (where i = electrical current and R = resistance) 17 a. Amplitude of cycle at frequency j is (aj2 + bj2)1/2 (i.e., the length of the vector in (sin,cos) space, or the complex conjugate (a+bi)(a-bi). [where i here is the imaginary number, square root of -1] b. By analogy, then, aj2 + bj2 = power at frequency fj (i.e., the (amplitude squared) of cycle j). c. The measurement units of power are amplitude2/cycle per time unit. d. The logic behind using power spectra for climate processes is debatable. Power spectra arose in the field of electrical engineering, so the choice is obvious in that context. It is not so obvious in the case of sea level; even if you would want to know the change in gravitational potential energy due to changing sea level, it is not proportional to the square of the amplitude. In fact, it might often be more useful to present amplitude spectra, which can simply be related to the wiggles one sees in the data (e.g. meters change of sea level). Typically, this is not done, solely because of the electrical origin of spectral analysis; statistical packages usually derive spectral estimates using programs written for electrical analysis. In any event, it is straightforward to take the square root of the power to get the amplitude. |<-Period=1/frequency->| 1 0.8 0.6 amplitude = sqrt(power) 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 2 4 6 8 10 12 18 2. Periodogram: 95%CI Power Frequency (cycles per time unit) 95%CI log(power) Frequency (cycles per time unit) 3. Problem: Statistical significance of periodogram a. White noise: equal energy (power) per frequency band ∆f (almost equivalent to random number series in time). b. Just as tossing a small number of coins is unlikely to give precisely 50% heads and tails, a finite time series is unlikely to have precisely equal energy per ∆f; in effect, "counting statistics" play a role; the longer the time series, the less the fluctuations within ∆f. c. Hence, a finite time series of a non-periodic function will show peaks and valleys in the periodogram. We would like to define criteria so that we can assess the likelihood that a peak observed in a real time series is due purely to chance. d. There are a number of "rules" around which purport to tell you "how many cycles" you need to prove a periodic process; typically these rules specify 5 or 10 cycles. While there is a grain of truth within these rules, they can be misleading. Is a record with 19 only 9 cycles unconvincing yet one with 11 cycles totally convincing? Of course, the statistics follow a continuum, so there is little statistical difference in the reliability of records with 9 or 11 cycles. Perhaps a more fundamentally misleading characteristic of these "rules" is that they ignore the relationship between the cycle you are trying to see and the characteristics of the rest of the variability. A very weak cycle (which nonetheless exists) can be hidden in the "noise" of the data or of the physical system; it can take many more than 10 cycles to verify the existence of a weak cycle. On the other hand, the characteristics of a very strong cycle can be estimated reliably with just a few cycles (although it is another thing to prove that those cycles couldn't have originated from chance). Basically, these "rules" (and even the more objective statistical tests) are oversimplifications; they do not excuse YOU from understanding what is required to prove the point that you are trying to make. e. For the periodogram, the following statistics apply to assess whether a peak of spectral density Pi is higher than expected from a white noise spectrum: Upper Bound Lower Bound 90% CI 19 Pi 0.33 Pi 95% CI 39 Pi 0.21 Pi Obviously these limits are not very good! If we are to see more subtle characteristics of the spectrum, it is necessary to improve these statistics. 4. Estimating a spectrum using a smoothing window (Daniell Estimator) a. As with most collections of numbers, we can improve error bounds by averaging (at the expense of detail). In a power spectrum, we can average adjacent frequency bands at the expense of frequency resolution. b. Several methods have been proposed to accomplish this averaging; for most geophysical time series, the simplest is adequate: the Daniell Estimator: k= +n Pj = ∑ P j +k k =−n 2n + 1 (Note: in order to keep the spectrum frequency bands the same as the periodogram bands when averaging even numbers of bands, it is common to add the n-1 bands on either side of the central band, and add 1/2 of the bands on either side of these n-1 bands). 20 c. Define "Equivalent Degrees of Freedom" (EDF) r = 2n where n = # of periodogram bands included in average. Chi-square table r 0.025 2 4 6 8 γ 0.05 .051 .484 1.237 1.690 .103 .711 1.635 2.733 _ Bound = (r/γ) Pj 0.95 0.975 5.991 9.488 12.592 15.507 7.378 11.143 14.449 17.535 (1) Example: suppose 2 periodogram bands are averaged (r=4) and we want to know the upper and lower bounds that must be exceeded to define a spectral peak at the 95% confidence interval: _ _ The upper bound is : 4/.484Pj = 8.3 Pj _ _ The lower bound is : 4/11.143 Pj = 0.36 Pj (2) Since the bounds are multiplicative, error bars are fixed on a log-scale diagram: n=2 n=4 _ | | | | | | o | | | _ n=10 _ | | | | o | | _ _ | | o | _ 5. Relationship between variance and power: 2 Σ (ai2 + bi2) ∆f = Σ Xi2 | | freq. bandwidth data 6. The resolution problem: From the definition of the fourier bands, it is obvious that resolution is poor for periods that are only a few multiples of the record length. What is one to do if you want to get a higher resolution spectrum. The safe answer is: Don't try to cheat mathematics: get a longer record! There is a more dangerous alternative, however. Suppose you throw out every other data point in your time series, then add a string of 21 trailing zeros. In effect, you have doubled the length of your time series without having added any data! While this may worry you, it doesn't bother the fourier transform at all, which will still happily compute the spectra. And it doesn't work too badly: in a test case where two sine waves corresponding to adjacent frequency bands in an original time series had the resolution doubled this way, the new FFT in fact was able to resolve the two sine frequencies. However, in doing so, it also generated "side bands" with significant amplitudes. So the extra resolution of these two bands was obtained at the expense of "contaminating" the rest of the spectrum. So in a real time series with a few strong frequency peaks, this method will compromise your ability to look at the rest of the spectrum. Caveat FFTer! 7. Optimal Spectral Estimation. a. Tapering results in some loss of spectral power (most severe at low frequencies) and (in some cases) loss of resolution between frequency bands. But it also minimizes "leakage", the tendency for spectral power to "leak" from one frequency into the estimate for adjacent frequencies. Example: a pure 41 kyr sine wave sampled 256 times during 510 kyr : 1.00 0.50 0.00 -0.50 -1.00 -1.50 0 100 200 300 400 500 22 1.00E+04 Detrended Tapered 1.00E+03 1.00E+02 1.00E+01 1.00E+00 0.000 1.00E-01 0.020 0.040 0.060 0.080 0.100 frequency, cycles per kyr 1.00E-02 1.00E-03 1.00E-04 Note that the detrended (nearly raw) time series shows spectral power at either side of the actual frequency which "isn't really there". Tapering improves this situation, but does not eliminate all of the side-bands. b. This observation leads to a concern over "optimal spectral estimation" i. The simple Fourier transform/periodogram method described below works reasonably well for simple situations where a small number of well-separated dominant frequencies are present, and where you are only interested in dominant power at the dominant frequencies and not at the true weak power in the gaps between the dominant time series. ii. In other situations, where you really care about extracting the best possible estimate of the true spectrum, you must move on to other techniques. Arguments have raged over optimal spectral estimation, but it seems that the winner is David Thompson, who has shown that the best approach is to use a "multi-taper" strategy in which the averaged Fourier transforms (from independent spectra computed using a set of "eigentapers") are used to retain as much information as possible. [Actually what he says is :"...spectral estimates based on the simple periodogram or unwindowed FFT (or equivalently, on sample autocorrelations) including indirect estimates (also known as Blackman-Tukey (1958) estimates) or worse, adding misfit to bias, the autoregressive or 'maximum-entropy' estimates, must be considered both obsolete and extremely dangerous." (Thompson, 1990)]. This approach is beyond the limits of this course, but it is the right thing to do if you must have the best spectra. 23 H. Coherence: How well are two time-series correlated within a frequency band [and what is the phase relation between the two time series] ? Phase Lag time --> highly coherent (with a phase lag) time--> 1. Coherence can be thought of as the "correlation-coefficient in frequency-space"; it is a number ranging from 0 to 1, with 1 reflecting perfect coherence; the statistical significance of a coherence estimate depends on the number of bands averaged. 24 2. Given two 'records' (X,Y) each has a power spectrum a2 + b2 = (a + ib)(a - ib) (i.e., the power spectrum is the complex conjugate of the fourier coefficents) ΦX(S) = Σ X(S) X*(S) ΦY(S) = Σ Y(S) Y*(S) (where * denotes complex conjugate; both ΦX(S) and ΦY(S)are real and >0) 3. The cross-spectrum ΦXY is similarly defined: ΦXY = Σ X(S) Y*(S) (aj+bji)(cj-dji) = ajcj + bjdj - (ajdj-cjbj)i = C + iQ (a complex number) 4. Coherence is the normalized cross-spectrum (the cross-spectral power normalized to the power of the two series) coh = Φxy _________ = C + iQ (a complex number) (ΦxΦy)1/2 0 ≤ coh(x,y) ≤ 1 a. Band-averaging must be done (same definition as Daniell estimator), except that when computing phase, you must band-average the complex numbers before computing phase (think of the beating of nearby frequencies...). This is done as complex addition of vectors in the phase plane coh = _ _ C2 + Q2 _______ _ _ P1 P2 _ _ coherence magnitude = (C2+Q2)1/2 phase _ _ -1 φ = tan (Q/C) 25 Q C Q Φ C 5. Example: X aj + ibj Y cj + idj XY = (a+ib)(c-id) ajcj + bjdj = Cj bjcj - ajdj = Qj Then band-average Cj + iQj (Daniell estimator) 26 Coherence Phase 1. cross-spectrum 2. compute coherence 3. Avg -> coherence 1. cross-spectrum 2. band-average 3. Compute phases 6. Statistical significance of coherence and phase. a. Coherence: critical values for the test of hypothesis q=0 n 2 3 4 5 6 7 10 95% CI 0.975 0.881 0.795 0.726 0.671 0.627 0.533 90% CI 0.948 0.827 0.662 0.475 However, these statistics only apply to time series where the time scales are precisely accurate. Paleoclimatic time series are often "tuned" to match orbital variability - in this case the coherence between orbital parameters and the coherence of the tuned time series is falsely elevated and isn't a particularly useful statistic (ref.). b. Phase error = arcsin [ [(1-coh2)/(coh2(2n-2))]1/2 t2n-2 (α/2)] where t is the t-statistic and α is the desired confidence level, i.e. for 95% CI, α = 1-.05. H. Another useful technique is frequency filtering of data. There are several ways this can be done; the simplest is to Fourier Transform the data, multiply the coefficients by multipliers (0 to 1) (warning: multipliers should change slowly in adjacent bands; i.e. don't go from 0 to 1 in adjacent bands), and then perform the Inverse Fourier Transform. V. Spectral analysis of oceanic climate records. A. Imbrie and others made attempts many years ago, but problems with the time scale and other not-fully-appreciated problems (e.g. bioturbation) precluded useful results. The imperfection of the time scale degraded the statistical significance of peaks that were more significant than they appeared. [This should also caution us not to always accept the criticisms about the appropriate significance levels for geological time series when the inherent limitations the time scale accuracy places on obtainable significance is not appreciated.]. By the early 1970’s, the time scale based on Barbados U/Th coral dates and linear interpolation from the most recent paleomagnetic reversal (Brunhes-Matuyama) in core V28-238 (Shackleton and Opdyke, 1973) provided a time scale that was good enough for the first attack on this problem. [Note: at the time that those papers were written, geophysicists had assigned a date 27 of 730 kyrBP to the B-M reversal based on K-Ar dates. More recent 39Ar/40Ar dating work (Baksi et al. 1992) indicates that this reversal is closer to 780 kyrBP.] B. First widely accepted application was a 1976 paper by Hays, Imbrie, and Shackleton, in which spectral analysis was applied to a 450,000 year record of (radiolarian Imbrie-Kipp) T, δ18O, and the abundance of a particular radiolarian, C. davisiana. δ18 O, ‰PDB RC11-120, E49-18 Patch 1.5 2.0 2.5 3.0 3.5 4.0 0 50 100 150 200 250 300 350 400 450 500 Age, ka bp (SPECMAP) 1. They found that spectral analysis of all of these records produced spectral peaks near 100, 41, and 23-19 kyr (i.e., same as eccentricity, obliquity, and precession. 2. The peaks appeared to be statistically significant (uhm…, sort of: their statistical test, which is differs from that given above, has been criticized and never successfully defended...) 3. This statistical coincidence earned the orbital variations the title of "pacemaker of the ice ages". 4. Because the time scale has relatively large uncertainties in the region between termination I and the B-M magnetic reversal, the phases of the cycles were not coherent with the orbital computations. By making numerous small adjustments to the time scale, well within the uncertainties allowed geologically, it is possible to maintain a consistent phase relationship between the orbital indicators and the climate indicators for the last 450,000 years. a. This finding gave rise to the concept of "tuned" time scales" (e.g., SPECMAP time scale). Basically, in the absence of radiometric age control, the logic is: i. Within the errors of the geologic time scale, and in particular during the intervals where age control is good, orbital parameters seem to maintain a consistent phase relationship with δ18O variations. ii. So, ASSUME that this phase relationship holds for times when the age control is not very good, and use perturbations of the original time scale to bring the precession and obliquity cycles into consistent phase relationships. iii. If we consider this technique only as a chronological tool subject to the eventual proof of the assumptions, this method is a reasonable means for interpolating time between known control points. But SPECMAP then proposed to do something more dangerous: to use data on this new time scale to evaluate the orbital hypothesis! This might seem circular, and it is in part, but it isn't completely, since 28 the aim is to test parts of the orbital hypothesis that weren't tuned: e.g., the eccentricity cycle, and other indicators that weren't specifically phase-locked onto the orbital parameters. Also, to the extent that it is just used to examine the relative timing of events in different regions and of different properties using the δ18O time scale, the tuning isn't really a factor. But there are problems: in a tuned system, the coherence statistic is a poor indicator and was oft-abused in the SPECMAP-era literature. C. Ruddiman and McIntyre (1981) found that (foram Imbrie-Kipp) T near 42°N in the Atlantic showed a very strong precession-driven cycle, with little 41,000 power. One doesn't need an FFT to see this signal! This site and some other regions show strong and obvious orbital signals. D. That the earth's orbital parameters has some influence on the earth's climate is now widely accepted (with a small minority of holdouts). The question that is more controversial is: how significant is the role of orbital parameters on climate change relative to other causes of climate change? It is possible that other processes (internal drifts or resonances) are dominant, but that they simply lock onto the "pacemaker of the ice ages" (e.g., just as the energy for a heartbeat always comes from the body, not the pacemaker). How can we investigate this problem? 29 E. SPECMAP project (spectral mapping): examination of spatial differences in spectral response of climate indicators. 1. The SPECMAP planktonic O18 stack and time scale (Imbrie et al., 1984) 30 3. Introduced concept of phase wheel as a means of summarizing many different climate indicators from different places in a simple diagram. (negative phase) zero phase (in phase with orbital forcing) leading lagging (positive phase) indicators indicators 4. Martinson et al. (1987) benthic stack and time scale 5. Imbrie, Mix, McIntyre (1989): spatial variability of SST spectra in Atlantic 33 6. Imbrie et al. SPECMAP compilation papers (1992, 1993) 36 7. Problem: The earliest observed response is found in the Southern Hemisphere - but the climate system resembles Northern Hemisphere radiation curves! Why? 8. Shackleton et al. (1990) revision of magnetic reversal chronology: An earlier paper by R.G. Johnson had proposed to use the orbital cycles as the “tree rings” of marine sediments, and use them to establish a time scale independent of other constraints. This approach was dismissed by SPECMAP, but later, Shackleton et al. applied spectral analysis tuning techniques towards this technique and proposed that the ages of several paleomagnetic boundaries were significantly older than previously suggested (e.g., B-M at 780 kyrBP). Independent 39Ar/40Ar dating (Spell and McDougall, 1992; Baksi et al., 1992) confirmed these new age estimates suggesting that Ar loss made some (but not all) K/Ar dates too young. 9. “The 100 kyr problem”: Why is dominant cycle ~100 kyr for the past ~million years when there is very little 100K power in insolation spectra? If 100 kyr cycle results from eccentricity, why isn’t there a strong 440 kyr cycle as well? 10. “The Stage 11 Problem”: because stage 11 occurs during a prolonged period of low eccentricity (minimum of both100- and 400-kyr cycles), orbital forcing is weak. But the stage 11 interglaciation is just as strong as those coinciding with much stronger forcing. Why? 11. “The 41K World Problem”: prior to about 1 myr, the dominant cycle in marine δ18O is 41 kyr (Raymo et al, 1989). Why? Why the transition to the 100 kyr mode?