ABSTRACT

Atmospheric levels of CO₂ are commonly assumed to be a main driver of global climate. Independent empirical evidence suggests that the galactic cosmic ray flux (CRF) is linked to climate variability. Both drivers are presently discussed in the context of daily to millennial variations, although they should also operate over geological time scales. Here we analyze the reconstructed seawater paleotemperature record for the Phanerozoic (past 545 m.y.), and compare it with the variable CRF reaching Earth and with the reconstructed partial pressure of atmospheric CO₂ (pCO₂). We find that at least 66% of the variance in the paleotemperature trend could be attributed to CRF variations likely due to solar system passages through the spiral arms of the galaxy. Assuming that the entire residual variance in temperature is due solely to the CO₂ greenhouse effect, we propose a tentative upper limit to the long-term “equilibrium” warming effect of CO₂, one which is potentially lower than that based on general circulation models.

CLIMATE ON GEOLOGICAL TIME SCALES Return to TOC

The record of climate variations during the Phanerozoic (past 545 m.y.), based on temporal and spatial patterns of climate-sensitive sedimentary indicators, shows intervals of tens of millions of years duration characterized by predominantly colder or predominantly warmer episodes, called icehouses and greenhouses (Frakes et al., 1992 ), respectively (Fig. 1 ). Superimposed on these are higher-order climate oscillations, such as the waning and waxing of ice sheets during the past 1 m.y. The recurring icehouse/greenhouse intervals were postulated to be a consequence of a plethora of causative factors, from celestial to planetary, including geographic distribution of continents, oceanic circulation patterns, atmospheric composition, or any combination of these. Lately, the consensus opinion favors atmospheric CO₂ as a principal climate driver for most time scales, from billions of years (CO₂ supergreenhouse, snowball Earth [Kasting and Ackerman, 1986 ; Hoffman et al., 1998 ]), to decadal and annual (Intergovernmental Panel on Climate Change [IPCC], 2001 ). Past climate variations should therefore correlate positively with coeval atmospheric pCO₂ levels.

For the Phanerozoic, estimates of atmospheric pCO₂ levels (Fig. 1 ) are based on model consideration and proxy data. Presently, three such estimates exist: the GEOCARB III model (Berner and Kothavala, 2001 ) and its precursors, and the reconstructions of Berner and Streif (2001) and Rothman (2002) . These reconstructions rely, to a greater or lesser degree, on the same isotope databases of Veizer et al. (1999) . However, they produce internally inconsistent outcomes and the curves do not show any clear correlation with the paleoclimate record. The poor correlation of the modeled pCO₂ with the observed Phanerozoic climate trends (Frakes et al., 1992 ; Veizer et al., 2000 ; Boucot and Gray, 2001 ) suggest either that the pCO₂ models may be in need of improvement, or, if one of them is validated, that the CO₂ is not likely to be the principal climate driver. In that case, what could be an alternative driving force of climate on geological time scales?

Decompositions of the δ¹⁸O and paleoclimate trends¹ (Veizer et al., 2000 ) display a dominant cyclic component of 135 ± 9 m.y. For δ¹⁸O, this is regardless of the temporal resolution (on m.y. time scales) adopted for deconvolution of the signal. There are no terrestrial phenomena known that recur with this frequency, particularly taking into account the regular near-sinusoidal fashion (Fig. 1 ; Wallmann, 2001 ) of the δ¹⁸O data. This regular pattern implies that we may be looking at a reflection of celestial phenomena in the climate history of Earth.

CELESTIAL CLIMATE DRIVER Return to TOC

Growing evidence, such as the correlations between paleoclimate records and solar and cosmic ray activity indicators (e.g., ¹⁰Be, ¹⁴C), suggests that extraterrestrial phenomena are responsible for at least some climatic variability on time scales ranging from days to millennia (Friis-Christensen and Lassen, 1991 ; Tinsley and Deen, 1991 ; Soon et al., 1996 ; Svensmark, 1998 ; Beer et al., 2000 ; Egorova et al., 2000 ; Soon et al., 2000 ; Björck et al., 2001 ; Bond et al., 2001 ; Hodell et al., 2001 ; Kromer et al., 2001 ; Labitzke and Weber, 2001 ; Neff et al., 2001 ; Todd and Kniveton, 2001 ; Pang and Yau, 2002 ; Solanki, 2002 ). These correlations mostly surpass those, if any, for the coeval climate and CO₂. Empirical observations indicate that the climate link could be via solar wind modulation of the galactic cosmic ray flux (CRF) (Tinsley and Deen, 1991 ; Svensmark, 1998 ; Marsh and Svensmark, 2000 ; Todd and Kniveton, 2001 ; Shaviv, 2002a , 2002b ) because an increase in solar activity results not only in enhanced thermal energy flux, but also in more intense solar wind that attenuates the CRF reaching Earth. The CRF, in turn, correlates convincingly with the low-altitude cloud cover on time scales from days (Forbush phenomenon) to decades (sun spot cycle). The postulated causation sequence is therefore: brighter sun enhanced thermal flux + solar wind muted CRF less low-level clouds less albedo warmer climate. Diminished solar activity results in an opposite effect. The apparent departure from this pattern in the 1990s (Solanki, 2002 ) may prove to be a satellite calibration problem (Marsh and Svensmark, 2003 ). The CRF–cloud-cover–climate link is also physically feasible because the CRF governs the atmospheric ionization rate (Ney, 1959 ; Svensmark, 1998 ), and because recent theoretical and experimental studies (Dickenson, 1975 ; Harrison and Aplin, 2001 ; Eichkorn et al., 2002 ; Yu, 2002 ; Tinsley and Yu, 2003 ) relate the CRF to the formation of charged aerosols, which could serve as cloud condensation nuclei (CCN), as demonstrated independently by ground-based and airborne experiments (Harrison and Aplin, 2001 ; Eichkorn et al., 2002 ).

Despite all these empirical observations and correlations, the solar-CRF-climate link is still missing a robust physical formulation. It is for this reason that such a link is often understated (IPCC, 2001 ), but this may change when the advocated experimental tests (Kirkby, 2001 ) are carried out. The only solar-climate mechanism that presently has a robust understanding, is change in the integrated solar luminosity, but the centennial increase in solar constant (2–4 W m⁻²: Pang and Yau, 2002 ; Solanki, 2002 ) appears to have been insufficient to account for the observed 0.6 °C temperature increase (IPCC, 2001 ). An amplifier, such as the cloud/CRF link, is therefore required to account for the discrepancy. Note, however, that a similar, albeit not as large, amplifier is implicit also in the CO₂ alternative, because the centennial temperature rise in these models is due mostly to the potential, and to some extent theoretical, positive water vapor feedback (Pierrehumbert, 2002 ) coupled with “parameterized” clouds, not to the CO₂ itself.

In view of the above empirical observations, could it be that the celestial forcing is the primary climate driver on most time scales, including the geological ones? The large stadial-interstadial temperature variations of the latest 420,000 yr, which in the ice cores correlate with 80 ppm variations in atmospheric CO₂ (Petit et al., 1999 ), appear to argue against such an alternative. One should note, however, that it is not clear whether the CO₂ is the driver or is being driven by climate change, particularly since the CO₂ appears to lag by centuries behind the temperature changes (Petit et al., 1999 ; Fischer et al., 1999 ; Mudelsee, 2001 ; Monnin et al., 2001 ; Caillon et al., 2003 ; Clarke, 2003 ), thus potentially acting as an amplifier but not as a driver. Can the geological record shed more light on this conundrum?

Unlike the past century, where solar activity, atmospheric CO₂, and global temperatures were predominantly increasing, and unlike the ice cores with their unresolved cause and effect relationship of CO₂ and climate, the situation over the Phanerozoic is different, with all three variables exhibiting a non-monotonic behavior. This may enable decomposition of the global temperature changes into contributions from CO₂, CRF, and a residual. It may also help to settle the causative sequence because celestial phenomena cannot be driven by terrestrial forcing. Moreover, the inherent time scales required for the global climate system to reach equilibrium can be as large as several millennia, owing to the slow heat exchange between the oceans and the atmosphere, and to the slow ice sheet adjustment time. Thus, by estimating the effects of CO₂ over geological time scales, we may obtain the long-term “equilibrium” response of the global climate system.

Recently, Fields and Ellis (1999) and Shaviv (2002a , 2002b) proposed that the CRF reaching the planet has not only an extrinsic variability due to its attenuation by solar wind, but also an intrinsic one arising from a variable interstellar environment. For example, a nearby supernova could bathe the solar system with a higher CRF for many millennia, leave a detectable ⁶⁰Fe imprint in ocean-floor deposits, and perhaps even give rise to a “cosmic ray winter” (Fields and Ellis, 1999 ) due to increased cloudiness and planetary albedo. Shaviv (2002a , 2002b) proposed that a particularly large CRF variability should arise from passages of the solar system through the Milky Way's spiral arms that harbor most of the star formation activity. Such passages recur at 143 ± 10 m.y. intervals, similar to the 135 ± 9 m.y. recurrence of the paleoclimate data (Veizer et al., 2000 ). Unlike the extrinsic solar-induced CRF modulations, which change the ionization rate at the bottom of the troposphere by typically <10%, the galactic flux variations are much larger and are expected to be about an order of magnitude more effective. It is these intrinsic CRF variations that may be responsible for the long-term climate changes over the past 1 Ga. Specifically, the “icehouses” and the oxygen isotope cold intervals (Fig. 1 ) appear to coincide with times of high CRF (Fig. 2 ), as deconvolved from galactic diffusion models and exposure ages in iron meteorites (Shaviv, 2002a , 2002b ). The shorter-term annual to multi-millennial climatic effects, superimposed on this long-term baseline, would then reflect the extrinsic modulations of the CRF due to variable solar activity. Changes in orbital parameters and in solar and terrestrial magnetic fields may also potentially modulate this superimposed CRF-solar impact.

CORRELATION OF THE CRF AND PALEOTEMPERATURE DATA Return to TOC

In order to estimate the intrinsic CRF reaching Earth, we used a diffusion model that takes into account the geometry and dynamics of the spiral arms, and considers that cosmic rays are generated preferentially in these arms. We chose three sets of diffusion model parameters (Fig. 2 )², which span the entire range of CRF histories that are consistent with observational constraints, the latter limiting the period of CRF oscillations to P₀ = 143 ± 10 m.y. (Shaviv, 2002a , 2002b ). Because the statistical record of exposure ages for iron meteorites has Poisson noise, the CRF histories we used are not directly extracted from it but they are the smoothed output of the galactic diffusion models constrained to fit the meteoritic record (see Shaviv 2002b for further caveats).

We model the temperature anomaly using the generalized form of:

The constant A normalizes for the average ΔT while the term B × t describes a linear temporal trend in ΔT. A term of this form is expected due to the increasing solar luminosity during the Phanerozoic, but may also arise from a possible secular variation in the CRF reaching the solar system; for example, from a changing star formation rate. A contribution to this term may also arise from systematic errors in the detrending procedure of the δ¹⁸O data. The third term considers the possibility that CO₂ variations affect ΔT , but at this stage we assume that the term is zero and defer its discussion to subsequent text. The fourth term arises from the variable CRF , where g() describes the functional dependence between ΔT and , and D is the actual normalization.³

All data (temperature, CRF, and the CO₂ discussed later) are binned into 10 m.y. intervals and averaged using a 50 m.y. window running average. This is because the temporal resolution of the isotope databases and the derivative pCO₂ models are in the 10⁶ yr range, while that of the CRF is in the 10⁷ yr range. Although Shaviv (2002a , 2002b) discussed the secular variations in CRF for the entire planetary history, the complementary δ¹⁸O record is available only for the Phanerozoic. We therefore truncate our comparison at 520 m.y. B.P. (560 m.y. for the Berner and Streif reconstruction). This gives us N_meas = 53 (57) correlated ΔT_i and their corresponding predicted ΔT_model(t_i). Utilizing the three limiting models of CRF variations (Fig. 2 ), we tested our models by minimizing the residual variance between the model ΔT_model(t_i) and the observed ΔT_i. We find that models that include solely the terms A and B result in a large σ_min² of 95(°C)². However, once the term D, the CRF normalization, is included, the σ_min² reduces to 32–36(°C)², in accord with the remarkable inverse correlation of CRF with the paleotemperature (Fig. 2 ). The CRF alone can explain 66% of the total variance in the temperature data. Can we further constrain the uncertainties in these models?

The only error on which we have a good handle is the statistical variance arising from the experimental δ¹⁸O data of Veizer et al. (1999) . From the internal variance of the δ¹⁸O data within the bins, we can calculate σ_min² expected from this source of error.⁴ This would be the minimum residual statistically attainable if we had perfect knowledge of all sources of climatic factors, exact CRF history, and no other error. This minimum variance, σ¹⁸_min², is found to be about 12(°C)² for models including the CRF. Thus, once we introduce CRF as a driver and remove the intrinsic δ¹⁸O measurement variance, we can explain 75% of the paleotemperature variability.

In addition to the δ¹⁸O measurement errors, additional errors may arise, for example, from translation of the δ¹⁸O data into ΔTs that required assumptions on the ice sheet volumes (Veizer et al., 2000 ), from an inaccurate CRF (e.g., inaccurate knowledge of spiral arm width, amplitude, and exact phase), or from additional factors that may affect the climate (e.g., CO₂, continental geography, oceanic circulation). The magnitude of such a “compound error” and its statistics can be estimated by the bootstrap method.⁵ Using this method, we can rule out a fluke correlation between the CRF and temperature at the 99.5% level. That is, we can rule out with a high confidence level models that do not include the effects of a variable CRF. This conclusion rests on the reasonable assumption that at least one of the two “celestial” data sets with the apparent 150 m.y. periodicity, the galactic spiral arm analyses or the iron meteorites exposure ages, is valid. While the above correlations are unlikely to be statistical flukes, we do emphasize that the data sets come with some caveats (see Shaviv, 2002b ). For example, although the variable meteoritic CRF signal is statistically significant, it could still be generated in 1.2% of random realizations. In another example, it appears that actually two spiral arm pattern speeds emerge from various astronomical analyses. While the number that fits the geological and meteoritic data is supported by a strong theoretical argument (Shaviv, 2002b ), the meaning of the second number is not yet resolved. Both numbers may be real, however, their meaning hinges on astrophysical considerations that are beyond the scope of this paper.⁶

Armed with the above statistics, we can then place quantitative limits on the CRF-climate connection. We tested 11 models (see Table 1 ), varying each variable, to find the range of values that gives reasonable fits at the 68% confidence level. The normalization parameter D for all these models varied between 3 and 12 °C. Almost all the error in D arises because we have no good limit on the amplitude of the variation of the CRF itself, except for the lower limit of 2.5 for its maximum/minimum ratio (Shaviv, 2002b ). We also find an average spiral arm passage period of P₀ = 137 ± 4 m.y., or 137 ± 7 m.y. if we consider the “jitter” from the epicyclic motion of the solar system (i.e., the noncircular motion around the Milky Way). This is consistent with the meteoritic data showing a periodicity of 143 ± 10 m.y. and the paleoclimate and paleotemperature data with a recurrence at 135 ± 9 m.y. To further check for consistency, we artificially add a lag to the predicted CRF. We find that the best lag is −3 ± 18 m.y. This implies that the results are consistent with our CRF diffusion model and astronomical data on the spiral arm location. They are only marginally consistent with other possible galactic models, which predict (Shaviv, 2002b ) that the actual spiral arm crossing took place 30 m.y. before the midpoint of the high CRF-climate episode. If we include an independent analysis of the lag in the correlation between the spiral arm passages and apexes of icehouses (Shaviv, 2002b ), we can exclude these alternative models at the 98% confidence level.

THE MODEL IMPACT OF CO₂ Return to TOC

Realizing that the pCO₂ reconstructions are internally inconsistent, the conservative point of view is to assume at the outset that the entire residual variance that is not explained by the measurement error is due to pCO₂ variations. From the model fit, we find that the temperature variance⁷ σT(CO₂)² attributable to such a “pCO₂” is at most (0.62°C)². To further quantify the effect of “pCO₂,” we need to know its variance. Considering that we are not aware of any mechanisms that would stabilize the Phanerozoic pCO₂ at today's values, particularly in view of the large sources and sinks, we assume that these variations span the entire range of the existing pCO₂ models (Fig. 1 ). With variations of this magnitude, the doubling of CO₂ can account for about (σT(CO₂)²/σln(pCO₂)²)^1/2ln(2) 0.5°C. A higher impact could be possible only if it is assumed that the Phanerozoic pCO₂ oscillations were limited to values close to the present-day levels.

It is entirely possible that none of the reconstructed Phanerozoic pCO₂ curves (Fig. 1 ) is a true representation of reality. Nonetheless, we also tested eight scenarios that assume that one of the reconstructed Phanerozoic pCO₂ trends (Fig. 1 ) is validated. To do this, we reintroduced the third term into equation 1 and considered the impact on temperature by a combined CRF and CO₂ forcing.⁸ We find that, depending on the model, the introduction of CO₂ as a driver reduces the σ_min² [32.1(°C)²] by only 0.5–1.5(°C)², compared to 0.2–5(°C)² for models that do not include CRF as a driver. That is, there is no statistically significant correlation between pCO₂ and reconstructed temperature, and we cannot therefore estimate the actual driving impact of CO₂. We can, however, estimate the upper bounds of model uncertainty in terms of temperature that, potentially, could be attributable to CO₂ forcing. This we can do by looking at the errors on the parameter C. Such formal 90% confidence limits are 0.91, 0.92, and 1.14 °C for the Berner and Streif, GEOCARB III, and Rothman reconstructions, respectively. At the 99% confidence limit they are 1.67, 1.46, and 1.93 °C (Table 2 ). In summary, we find that with none of the CO₂ reconstructions can the doubling effect of CO₂ on low-latitude sea temperatures be larger than 1.9 °C, with the expected value being closer to 0.5 °C. These results differ somewhat from the predictions of the general circulation models (GCMs) (IPCC, 2001 ), which typically imply a CO₂ doubling effect of 1.5–5.5 °C global warming, but they are consistent with alternative lower estimates of 0.6–1.6 °C (Lindzen, 1997 ).

As a qualifier, one should note that global temperature changes should exceed the tropical ones because the largest temperature variations are in the high-latitude regions for which we do not have any isotope record. A review of GCMs (IPCC, 2001 ) shows that the globally averaged warming from CO₂ is expected to be typically 1.5 times larger than that of the tropical temperatures, and our model uncertainty limits should therefore be modified accordingly. Note also that the bootstrapping “compound error” includes, among others, any error associated with the ice volume correction. Taking an unrealistic ultimate scenario that assumes no ice volume correction at all, the amplitude of temperature oscillations in Figure 2 could be almost doubled. While these and similar considerations may help in expanding somewhat the above calculated temperature limits potentially attributable to CO₂, they will not alter the relative importance of the celestial-CO₂ forcings. The model impact of CRF will increase in tandem with that of CO₂ for any change in the amplitude of ΔT. As a final qualification, we emphasize that our conclusion about the dominance of the CRF over climate variability is valid only on multimillion year time scales. At shorter time scales, other climatic factors may play an important role, but note that many authors (see previous references) suggest a decisive role for the celestial driver also on multi-millennial to less than annual time scales.

POTENTIAL IMPLICATIONS Return to TOC

Our approach, based on entirely independent studies from astrophysics and geosciences, yields a surprisingly consistent picture of climate evolution on geological time scales. At a minimum, the results demonstrate that the approach is potentially viable, as is the proposition that celestial phenomena may be important for understanding the vagaries of the planetary climate. Pending further confirmation, one interpretation of the above result could be that the global climate possesses a stabilizing negative feedback. A likely candidate for such a feedback is cloud cover (Lindzen, 1997 ; Ou, 2001 ). If so, it would imply that the water cycle is the thermostat of climate dynamics, acting both as a positive (water vapor) and negative (clouds) feedback, with the carbon cycle “piggybacking” on, and being modified by, the water cycle (Nemani et al., 2002 ; Lovett, 2002 ; Lee and Veizer, 2003 ). It is our hope that this study may contribute to our understanding of the complexities of climate dynamics and ultimately to quantification of its response to potential anthropogenic impact.

ACKNOWLEDGMENTS

This research was supported by the F.I.R.S.T. (Bikura) program of the Israel Science Foundation (grant no. 4048/03), the Deutsche Forschungsgemeinschaft, the Natural Sciences and Engineering Council of Canada, and the Canadian Institute for Advanced Research.

REFERENCES CITED Return to TOC

¹The oxygen isotope record is based on 4500 measurements of shells (brachiopods, belemnites, foraminifera) composed of low-Mg calcite, the carbonate phase most resistant to post-depositional overprint of the signal. The data show a secular ¹⁸O depletion trend with age, with superimposed higher order oscillations and it is these that are in phase with reconstructions of the Phanerozoic climate. The major features of these oscillations represent a robust signal in and of themselves, which would be reproduced even if only a fraction of the samples, those near the upper envelope of the secular trend, were taken into consideration. It is this detrended oscillating pattern, with correction scaled for ice volume effect, that yields the Phanerozoic ΔT [°C] variations for contemporaneous low-latitude shallow sea water in

²The observational constraints (for P₀ = 143 ± 10 m.y.) include the cosmic ray ¹⁰Be age, and limits on CRF variations derived from iron meteorites. The three models that we utilize (

³g() is defined such that g() = 0, 1 for (t) = 0,(today) respectively. (t,P₀) itself is one of the three CRF histories used (

⁴

σ_min²Σ_i(ΔT_i−ΔT_model(t_i))².

⁵If we had a perfect model and knowledge of the errors, then the χ² of the fit should, on average, be the number of actual degrees of freedom. We therefore add errors quadratically to increase the error in the data until the modified χ² per degree of freedom is 1. If we further assume that the unknown measurement and model errors have a Gaussian distribution, we can estimate the errors in the fit parameters. To check the assumption of Gaussianity, we look at the distribution of the residual differences between the model and the best fit (model 6 in

⁶Although the astronomical data points to two possible spiral arm pattern speeds, a theoretical argument based on the observed outer extent of the galactic spiral arm and the spiral arm density wave theory can be made (

⁷

⁸ In order to calculate the impact of combined CRF and CO₂ forcing we used the functional form f = 1.236 [ln(c + 0.0005c²) – ln(c₀ + 0.0005c₀²)] to which the radiative driving is expected to be proportional (

Click on thumbnail for full-sized image.

Figure 1. Phanerozoic climatic indicators and reconstructed pCO₂ levels. The bottom set of curves are the detrended running means of δ¹⁸O values of calcitic shells over the Phanerozoic (Veizer et al., 2000 ). 3/6, 5/10, 10/20 and 10/50 indicate running means at various temporal resolutions (e.g., 3/6 means step 3 m.y., window 6 m.y. averaging). The paleolatitudinal distribution of ice rafted debris (PIRD) is on the right-hand vertical axis. The available, Paleozoic, frequency histograms of other glacial deposits (OGD)—such as tillites and glacial marine strata— are dimensionless. The blue bars at the top represent cool climate modes (icehouses) and the white bars are the warm modes (greenhouses), as established from sedimentological criteria (Frakes and Francis, 1998 ; Frakes et al., 1992 ). The lighter blue shading for the Jurassic- Cretaceous icehouse reflects the fact that true polar ice caps have not been documented for this time interval. The upper set of curves describes the reconstructed histories of the past pCO₂ variations (GEOCARB III) by Berner and Kothavala (2001) , Berner and Streif (2001) and Rothman (2002) . The pCO₂(0) is the present-day atmospheric CO₂ concentration. All data are smoothed using a running average of 50 m.y. with 10 m.y. bins. The hatched regions depict the uncertainties quoted in the Rothman and the GEOCARB reconstructions.

Click on thumbnail for full-sized image.

Figure 2. The cosmic ray flux () and tropical temperature anomaly (ΔT) variations over the Phanerozoic. The upper curves describe the reconstructed CRF using iron meteorite exposure age data (Shaviv, 2002b ). The blue line depicts the nominal CRF, while the yellow shading delineates the allowed error range. The two dashed curves are additional CRF reconstructions that fit within the acceptable range (together with the blue line, these three curves denote the three CRF reconstructions used in the model simulations). The red curve describes the nominal CRF reconstruction after its period was fine tuned to best fit the low-latitude temperature anomaly (i.e., it is the “blue” reconstruction, after the exact CRF periodicity was fine tuned, within the CRF reconstruction error). The bottom black curve depicts the 10/50 m.y. (see Fig. 1 ) smoothed temperature anomaly (ΔT) from Veizer et al. (2000) . The red line is the predicted ΔT_model for the red curve above, taking into account also the secular long-term linear contribution (term B × t in equation 1). The green line is the residual. The largest residual is at 250 m.y. B.P., where only a few measurements of δ¹⁸O exist due to the dearth of fossils subsequent to the largest extinction event in Earth history. The top blue bars are as in Figure 1 .