Dependence of global temperatures on atmospheric CO_{2} and solar irradiance
See allHide authors and affiliations
Abstract
Changes in global average temperatures and of the seasonal cycle are strongly coupled to the concentration of atmospheric CO_{2}. I estimate transfer functions from changes in atmospheric CO_{2} and from changes in solar irradiance to hemispheric temperatures that have been corrected for the effects of precession. They show that changes from CO_{2} over the last century are about three times larger than those from changes in solar irradiance. The increase in global average temperature during the last century is at least 20 times the SD of the residual temperature series left when the effects of CO_{2} and changes in solar irradiance are subtracted.
Although it is generally conceded that the average surface temperature of the Earth has increased by about 0.6°C during the last century, there is little agreement on the cause of this warming. The primary cause of this disagreement is uncertainty about the relative contribution to this warming of atmospheric CO_{2} and changes in solar irradiance. The purpose of this paper is to describe some data analysis that may help to discriminate between solar and CO_{2} effects, and to give estimates of the relative magnitudes of these two effects. The difference between analysis such as those described in ref. 1 and those here is that this data analysis is based on deseasonalized temperature time series where the effects of precession were included.* The detection of precession in instrumental temperature series and the necessity of including it when removing the annual cycle from temperature data was demonstrated in ref. 2. I also describe some of the statistical peculiarities and limitations of these data series and suggest where better data are needed.
The paper begins with a discussion of the data being analyzed and, to delineate the issues, presents some ordinary leastsquares fits of the temperature data with atmospheric CO_{2} concentration and changes in solar irradiance. I next discuss the mathematical methods used and describe some statistical properties of the various data series. This is followed by some simple estimates of the transfer functions between fossil fuel consumption and atmospheric CO_{2} levels and from CO_{2} levels and changes in solar irradiance to temperature. The penultimate section summarizes recent findings on destabilization of the annual cycle, followed by conclusions.
In these analyses I do not directly take into account the effects of stratospheric aerosols nor various internal feedback mechanisms such as cloud cover. Stratospheric aerosols are generally believed (3, 4) to result in cooling, so their omission makes the estimates for sensitivity conservative. Similarly, while a detailed understanding of internal feedback mechanisms, such as water vapor, is necessary to predict temperature changes from first principles, one may use measurements to assess the general climate response to forcing without having to consider the internal feedbacks explicitly, much as one can design a filter using operational amplifiers without detailed consideration of the quantum mechanics, or even current flow, in the individual transistors in the amplifiers.
The estimates given here depend neither on general circulation models nor on the assumptions that underlie such models. The transfer functions are estimated directly from observations of temperature and CO_{2} and, for solar irradiance, a physically based proxy data series.
Data Sources and Preparation
For measurements of surface air temperature I use the lowpass filtered Jones–Wigley Land plus Marine data (5) shown in figures 9 and 10 of ref. 2. The bandwidth of the lowpass filter was 0.5 cycle/year so the Nyquist rate is one sample per year. These series differ from the ones usually seen in two important aspects: First, I replaced the standard “deseasonalizing” procedure used to produce temperature anomaly series with a projection filter separation into lowpass, annual, and highfrequency components so, implicitly, the usual “boxcar” runningmean smoother has been replaced with a lowpass filter. Second, instead of assuming a constant amplitude climatology with a period of 1 calendar year, I allowed the phase of the annual components to track the observed phase. Thus, the significant changes in the annual cycle caused by the changing balance between direct insolation, periodic at one cycle per tropical year, and transported heat, periodic at one cycle per anomalistic year, has been removed from the data, eliminating the spurious monthly trends associated with temperature anomaly series (2, 6). Note that, although the timeresolution of these series is one year, the series is as smooth as that given by the usual “boxcar” procedure at decadescale resolution. The “Global” temperature series used here is the arithmetic average of the Northern and Southern Hemisphere series. The Northern Hemisphere, Southern Hemisphere, and Global series are denoted by T_{n}(t), T_{s}(t), and T_{g}(t) respectively, with t the Gregorian calendar date.
In this paper I use the average temperature over the 65 years from 1854 to 1918 as a base reference. There are several reasons to prefer this period to the usual 1951–1980 reference period. First, based on the sunspot record, solar activity in the 1854–1918 period appears to be representative of the 245year available record and 65 years covers most of the 88year Gleisberg cycle (7). Second, median fossil fuel consumption in this period was only about 6% of the current rate and the concentration of atmospheric CO_{2} increased by 4.7% between 1854 and 1918 (8), only onequarter of the 1951–1980 rate. Third, chloroflourocarbons and similar ozonedepleting chemicals were not yet in use, and as there is considerable evidence (9, 10) for stratospheric control of climate, it is desirable to use a reference period before changes in stratospheric chemistry by chloroflourocarbons began. Fourth, the major erratic changes in the timing of the seasons described in ref. 2 appear to have begun about 1920, so the early reference period is largely free of their effects. Opposed to these considerations are poorer spatial coverage of the temperature series and the presence of several major volcanic events (11) during the earlier period. There are no solar irradiance measurements without the confounding influence of the atmosphere in the early period and only a few near the end of the 1951–1980 period. Keeling’s CO_{2} measurements began in 1958, so the 1951–1980 period has better CO_{2} data than the earlier period. The average and median values of the 1854–1918 data are 171.9 mK and 180.0 mK (Northern Hemisphere) and 150.7 mK and 159.0 mK (Southern Hemisphere) below the 1951–1980 reference. Because the SD (σ) of the raw T_{g}(t) series during the 65year reference period is 57.7 mK, the 1990 temperature of 635 mK above the base temperature is, at a minimum, an 11 σ increase.
The CO_{2} data are as listed in table A.6 of ref. 8 up to 1955, and the Mauna Loa averages through 1994 from the Oak Ridge National Laboratory data set ndp001r5 (12) since 1955. The early data has been interpolated from irregular and inhomogeneous observations and, statistically, is too smooth. I denote this series by C(t), and log_{2}CO_{2}(t) by C_{L}(t). [Radiation theory predicts that temperature is a logarithmic function of atmospheric CO_{2} concentration (13). Use of the base 2 logarithm of the CO_{2} data gives coefficients that directly describe the effect of doubling CO_{2}.]
I have also used Marland’s fossilfuel production series (14), denoted F(t), as an adjunct to the CO_{2} measurements. This series starts in 1860, later than the temperature data, but the early data appear to be better than the corresponding CO_{2} data. The combination of these data with the CO_{2} data is described later.
For solar irradiance I used the Foucal–Lean (15) reconstruction. This series, L(t), is independent of the temperature data, matches direct solar irradiance measurements since they have been available, and is a reconstruction from other solar measurements before then. I emphasize, however, that this series is a proxy, not direct measurements, so that inferences drawn from it may not have the reliability of inferences from direct observations.
Changes in the period of the sunspot cycle were suggested (16) as a solar irradiance proxy. Although there is a high apparent correlation between this proxy and a heavily smoothed version of the Hansen temperature series (17), this correlation is not reliable. The jackknife variance (see below) of the lowfrequency coherence between the temperature and the sunspot period is large, more than four times that expected under Gaussian theory. Because of this, lack of a physical basis for the proxy, and the failure of other statistical tests described in ref. 2, this proxy is not used here.
LeastSquares Fits
Both CO_{2} and changes in solar irradiance have been invoked to explain the observed increase in global temperature. Fig. 1 shows the filtered Jones–Wigley global temperature series, T_{g}(t) together with a leastsquares fit to it using C(t) plus a constant, and a second leastsquares fit to T_{g}(t) using L(t) plus a constant. Each of these fits explains more than 75% of the variance over the full 1854–1990 interval (the residual SDs are 74.2 mK and 83.8 mK, respectively). Including both C(t) and L(t) simultaneously as explanatory variables further reduces the residual SD to 62.5 mK. Nearly 87% of the variance is explained, and both partial F statistics (18) are highly significant, so neither variable can be dropped. Examining these residuals further, one finds that their autocorrelation at a 1year lag is 0.914 so that the conditions required for the Gauss–Markov theorem (the basis of leastsquares) to be valid (19) are not satisfied. This is not simply a technical mathematical quibble, but indicates serious fitting problems whose existence may be verified by repeating the fitting process over different time intervals and observing the change in the estimated coefficients. For example, a leastsquares fit to T_{g}(t) with just the proxy solar irradiance L(t) and a constant the interval 1854–1918 gives a negative temperature response to increasing proxy solar irradiance. Consequently, one must use statistical timeseries methods that are reliable when serially correlated residuals are present.
Mathematical Preliminaries
Many of the problems in the analysis of climate data require new methods for timeseries analysis. Most of the commonly used timeseries methods were derived under the assumption of stationarity, that is, their derivations assume that the statistics of the observed process are independent of the choice of time origin. The climate seems to be nonstationary; over a few years the annual cycle of the seasons makes cyclostationary (or periodically correlated) processes a better model than stationary processes, implying that Fourier transforms of frequencies offset by multiples of one cycle per year will be correlated. On longer time scales, evolution of the Earth’s orbit obviously results in vast shifts in the climate, and I recently have shown (2) that these changes in the orbit must be considered in the analysis of instrumental data series as well. Here these effects are accounted for in the way the data were filtered. Because both solar and CO_{2} effects alter the seasonal cycle as well as low frequencies, I have not removed common modulation terms at 0 and 1 cycle/year. Finally, anthropogenic changes are altering the composition of the atmosphere and the climate system. Thus, analysis methods that contain an implicit assumption of stationarity should be used with caution.
Similarly, confidence intervals for parameters derived from time series are often based on an assumed Gaussian distribution. It appears that this assumption is invalid for climate data: there are outliers, predominantly from volcanic effects, and, surprisingly, some quantities appear to be more stable than one would expect from a Gaussian distribution.
Even discounting the intrinsic complications of the climate and solar data, the available data is far from homogeneous. There are, for example, no direct measurements of solar irradiance, and only a few direct CO_{2} measurements before 1958. In addition, the spatial coverage of the temperature averages has changed markedly. Fluctuations in the fossilfuel production series (14) appear to be roughly proportional to the level of production, so the absolute errors in the 19th century were lower than they are at present.
The methods used in this paper are a mixture of timedomain and frequencydomain techniques. The multiplewindow (20–22) method of spectrum estimation is used as an indirect basis for most estimates. This basic method is formally extended in several ways: Repeating the analysis with each window omitted in turn, or jackknifing over windows (23), gives estimates of the variances of spectra and coherences that do not depend on a Gaussian assumption. Quadraticinverse theory (24) tests for correlations between closely spaced frequencies, and a variant (25) allows checking for nonstationary effects. The singular value decomposition of a multiplewindow spectrogram (21) provides insight into some of the stationarity issues. In addition, a variety of prewhitening and robust estimation methods typified by those in (26, 27) have been used throughout. The Jones–Wigley data, however, has been so carefully prepared that robust procedures detect little more than volcanic activity and do not significantly change the results. Space limitations preclude including many of these checks and comparisons, and I have omitted those that give expected results.
The stochastic structure of both climate and solar data is amazingly complicated when compared to that commonly encountered in textbook examples or even most engineering problems. Because analysis methods for nonstationary time series do not appear to be well enough developed to cope with these series, nor the stochastic structure of the data well enough understood, I use methods that allow simple models for the trends and forcings while leaving nearly stationary residuals.
Stochastic Properties of the Data Series
Because one cannot reliably estimate the relationships between random series without understanding their internal structure, I began this study with a brief examination of the properties of the individual data series. Power spectra of the various series have the usual “red” characteristics of much geophysical data and are not shown. Most of the highenergy, lowfrequency part of the spectrum is a result of the trends under consideration, and the stochastic structure is more obvious if the data are prewhitened by filtering. Because the spectra appear similar, I computed the geometric mean of the spectra of T_{g}(t), L(t), and C_{L}(t), Fourier transformed this average spectrum to give an autocorrelation sequence, then computed a fourthorder autoregressive prewhitening filter from the autocorrelations. I then use this AR4 filter on all the data series and recompute the spectra. The spectra have been computed with a multiplewindow (20) estimate typically using a timebandwidth product of 6.0 and the 10 lowestorder windows so the spectrum estimates have a nominal χsquare distribution with 20 degrees of freedom. Their 5% and 95% confidence intervals have been determined by jackknifing over windows (23). The jackknife method is known analytically to give good performance over a range of distributions, to be reasonably distribution free, and, empirically, to be sensitive both to unresolved structure in the spectrum and to nonstationarity. The jackknife variances of the spectrum of the prewhitened temperature series are, on average, about 87% of that expected under Gaussian theory (28).
In contrast, the jackknife variance of the spectrum of L(t), Fig. 2, is nearly 10 times that expected near the 22year Hale cycle and considerably lower than expected around the 11year sunspot cycle. In a stationary Gaussian series such excursions would be extremely rare.
Examining the cause for the low jackknife variance, the quadraticinverse test for unresolved spectral details (24), Fig. 3, for the Northern Hemisphere is lower than expected except at low frequencies. (The geometric mean, across frequencies, is 3.1, compared to an expected value of 4, and the estimate is below the expected value over much of the frequency range.) This implies that the sampling variations within the analysis bandwidth in the estimated spectrum are smaller than one would expect in a Gaussian random process with a constant spectrum. From this, I infer that much of the apparently random structure in the temperature data is not random but, more probably, a consequence of either deterministic forcing or internal oscillations. Because of the similarity between the statistics of the solar and temperature data, and the stationary temperature residuals (described below), solar forcing is a more likely explanation than internal oscillations. The unresolved structure test in both hemispheres have local maxima near the 22year Hale cycle and significant maxima at low frequencies.
Fig. 4 is a plot of an estimate of the frequency derivative d/df ln{S(f)} [strictly, the relative quadraticinverse coefficients, b_{1}(f)/S̄(f) in the notation of ref. 24] for the two hemispheric temperature series and L(t). The higherorder coefficients have similar characteristics, but lower amplitudes. The similarity of the lowfrequency features in the temperature and irradiance data lend credence to many of the suggested solarclimate relations and also imply that extracting the statistical details of these relationships will be complicated and difficult.
Turning to the lowfrequency maximum in Fig. 3, narrower bandwidth spectrum estimates hint at unusual activity at low frequencies in both T_{n} and T_{s}. The harmonic Ftest shows moderate evidence for a periodic component with an estimated period of 180 years. The estimated period is longer than the length of the data series, and while the nominal significance level is 99.7%, this is not much higher than 1 − 1/N = 0.993, the level where one expects at least one false detection of a line in noise. However, because the ±1 − σ confidence interval on the period of 153 to 220 years includes the 208year Suess period, the signal is detected in both hemispheres, and the magnitudesquared coherence between them is about 0.6, the possibility that there is a complicated natural signal in the data near a period of 200 years should not be dismissed. The similarity of the filtered amplitude estimates, 4.4 mK and 4.0 mK, and phases, −135.0° and −138.5° (referenced to 1900.0), in the Northern and Southern Hemispheres, respectively, further support this conclusion. Analysis of long series of ^{14}C data (21) shows significant components at periods of 230.6, 208.4, and 199.4 years, so because the frequency difference between the latter two periods corresponds to 4,600 years, one cannot expect these components to be resolved in the 137year series, nor the single line Ftest to work reliably.
A possible alternative explanation to deterministic forcing as the cause of the low variances of the temperature spectra is nonstationarity, and this was tested for using the methods in ref. 25. This test shows less variability in time than expected at periods longer than 8 years. Thus, near the 11year solar cycle, both the tests for unresolved structure and stationarity show more stability than one would expect in a random series. Examination of the nonstationary quadraticinverse coefficients, a_{m}(f), shows that the a_{2}(f)s, approximately, the second timederivatives of the spectra† are reasonably significant and similar in both hemispheres.
Looking at details of the spectrum of the solar irradiance series, L(t), the predominant feature is the character of the jackknife variance (Fig. 2). There is significant excess variability near 22year periods and a paucity around 11year periods. The unresolved structure test is uniformly higher than expected out to a frequency of about 0.13 cycle/year, that is, until the highfrequency edge of the 11year band, followed by a low region. The nonstationary quadraticinverse coefficient a_{1}(f), roughly the timederivative of the spectrum, or Ṡ(f) in both hemispheres shows decreasing power near 22year periods and increasing power in the 7 to 11year period band.
The singular value decomposition of the logspectrogram described in ref. 21 suggests mild nonstationarity in T_{n}(t), T_{s}(t), and L(t) concentrated about the 104year Suess period, and possibly related frequencies. The phases of the annual cycle of the temperature data shows similar characteristics. Note that this 104year period is not a simple additive term, but a periodic modulation of the stochastic structure of the process.
I emphasize that the statistical peculiarities just described are those of the individual data series; they contrast markedly with the statistics of the temperature residuals to be described. To the extent that the climate has a linear response to forcing, nonstationarity in the forcing should appear as a similar nonstationarity in the temperature series so, internal oscillations aside, removing the effects of such forcing should leave stationary residuals.
The estimate of the spectrum of the log_{2} CO_{2}(t) data, shown in Fig. 5. shows two important features: first, the range of the spectrum, seven decades, is larger than that of either the temperature or irradiance data, and second, the jackknife variance estimate of ln Ŝ(f) has a maximum of 12.06 at periods near 62 years, compared to an expected value of about 0.13, or 90 times expected. Thus the 5% to 95% confidence region for the CO_{2} spectrum estimated by the jackknife method covers a range of about 10^{4}, instead of the factor of 2 that would be expected from Gaussian theory. This uncertainty in the spectrum of CO_{2} variations carries through to the coherences and transfer function estimates and, consequently, standard frequencydomain estimates of transfer functions are not presented here.
Temperature Response Models and Transfer Function Estimates
I attempted to assess the simultaneous effects of solar variability and greenhouse gases by estimating transfer functions. Because both effects have been advanced as plausible explanations for the increasing temperature, one should expect that changes in solar irradiance and greenhouse gas concentrations should appear coherent as well. It has been shown (29) that, considered as a pair, changes in temperature and CO_{2} were coherent, with changes in temperature leading those in CO_{2}, instead of viceversa, as popularly supposed. Marland’s fossilfuel production series, in contrast, leads both the temperature and CO_{2} data. Changes in solar irradiance also lead those in both CO_{2} and temperature.
In addition to these complicated relationships, there are other difficulties to consider; the physically distributed and multicomponent nature of the problem, combined with spatially and temporally discrete data, makes the appropriate form of the transfer functions obscure. The roughly periodic character of the solar cycle makes leadlag relations ambiguous. In an attempt to untangle these relationships, I tried several different forms of estimates; the most successful of these is described below. Here I define “successful” as low prediction error, few free parameters, and physically realistic implications. For example, I consider a negative lowfrequency response to either increases in solar irradiance or CO_{2} unacceptable. Transfer functions implying noncausal response to changes in solar irradiance imply either that the transfer function is unacceptable, or that there is a problem with the assumed irradiance. I also require that the estimates be reasonably insensitive to linear filtering operations applied to the raw data, which eliminates many of the possible models. The rationale for this requirement is twofold: first, if there is a linear relationship between the input and output of a system, estimates of the transfer function should not change significantly when the same linear filter is applied to both input and output, and second, residuals in the series are autocorrelated, so if the filtering reduces the residual autocorrelation, the requirements of the Gauss–Markov theorem are more nearly satisfied. The estimated transfer functions are clearly not unique, but other acceptable estimates with similar prediction errors have similar implications.
A simple model for the atmospheric CO_{2} concentration at time t is 1 where the coefficient ζ on the previous concentration allows for exchange and sequestration by the oceans and biosphere, β is the part of the carbon present in fossil fuels that contributes to atmospheric CO_{2}, and γ allows for the temperature forcing described in ref. 29. Doing a robust fit with the more accurate Keeling data weighted a factor of 6 higher than the pre1958 data gives ζ̂ = 0.999974, implying that the atmosphere acts, as expected, as an integrator. The estimated coefficients are β̂ = 0.524, somewhat lower than the estimate in ref. 8 or direct estimates between CO_{2} and the integrated fossil fuel record. This is offset by the temperature feedback, γ̂ = 1.196 GT/K. Replacing the T_{g}(t) term with L(t) gives a slightly poorer fit with a coefficient ≈ 0.38 GT/(W/m^{2}).
In the following paragraphs, I describe some simple timedomain estimates of the transfer functions. If one takes the simplest conceptual model 2 and substitutes Eq. 1 for C(t) one obtains, approximately, the equation 3 with suitable definitions of the constants, and F(t) incorporated into C_{L}(t). This can be generalized to the class of models with suitable choices of A, P, and Q. When A = 0 the autoregressive feedback terms are absent, and the model that of a direct finite impulse response filter from L(t) and C(t) to temperature. In particular the models A, P, Q = 0, 1, 0 are the direct leastsquares fits, mentioned earlier, to L(t) only; 0, 0, 1 to C(t) only; and 0, 1, 1 to both simultaneously.
I choose the interval 1854–1965 to estimate the coefficients for the particular model and use the 25year interval, 1966–1990, for validation. Reconsider the direct leastsquares estimate model 0, 1, 1. The coefficients obtained by fitting from 1854–1965 are ŝ_{0} = 0.1115, ĉ_{0} = 1.984. The meansquare error in the fitting period is 6,695 (mK)^{2}, or an rms error of 81.8 mK. This fit is not particularly impressive but, nonetheless, it suggests that changes in solar irradiance lead those in temperature by a few years. Using the data for the 25 years 1966–1990, with the coefficients found for the earlier period gives a prediction error of 39,820 (mK)^{2}, for an rms error of 199.5 mK. Because the autocorrelations of the residuals in all these cases are large, the assumptions of the Gauss–Markov theorem are violated, and ordinary leastsquares does not give reliable results. To see how badly, replace Eq. 2 with 4 where Δ_{r}T(t) = T(t) − rT(t − 1) and similarly with L and C. Using the 1year autocorrelation of the temperature data, r = 0.881 for all three series gives s̃_{0} = −0.0018, c̃ = 3.121 with corresponding fitting and prediction errors of 1,369 (mK)^{2} and 2,436 (mK)^{2} respectively. Because the innovations variance of the temperature data where B(v) ≈ 1.07 is a bias correction depending on the degreesoffreedom (30) of the estimate and Ŝ(f) the estimated spectrum at frequency f, is about 1,303 (mK)^{2} and the variance of Δ_{r}T is 1,835 (mK)^{2} the errors in Eq. 4 are far less correlated than are those of Eq. 2, so the conditions of the Gauss–Markov theorem are closer to being satisfied. Note that using this simplest of filters on the data increases the estimated sensitivity to C_{L} by 57% and changes the sign of ŝ_{0}. Because a decrease in average temperature as a response to an increase in solar irradiance is implausible, I reject such simple models.
Direct frequencydomain estimates of the transfer functions appear to vary too rapidly to be reliably estimated from the available data as their Fourier transforms are significantly noncausal. For example, the impulse response estimated between T_{n}(t) and L(t) is noncausal with a width about the length of a solar cycle. This is an artifact of the near periodicity of the solar cycle as, without other knowledge or assumptions, one cannot assign the response to a periodic process to a given time without ambiguities of N periods. Similar problems are common in transfer function estimation problems ranging from magnetotellurics (31) to the design of telephone equalizers.
The estimated impulse response between C_{L}(t) and T(t) is causal and oscillatory with about a 20year duration. If this characteristic persists with improved early CO_{2} data, it implies that large transient responses may be possible.
In addition to thermal inertia, the observed feedback from temperature to CO_{2} also generates persistence. To model this together with the effects of CO_{2} and solar irradiance changes on the temperature record, I fit the filtered Jones–Wigley series as a 1, 1, 1 model, explicitly Eq. 3. Here the feedback term αT(t − 1) results in a rational (as opposed to polynomial) approximation of the transfer functions allowing more rapid lowfrequency changes in them than possible with simple models, it includes some persistence, and generally permits simple feedback effects. Direct leastsquares estimates of the coefficients and some performance measures are given in Table 1.
To assess the model’s predictive power I start at 1966 and iterate forward to t = 1990 using the observed solar irradiance and CO_{2} data. The predictions were run in two ways: in the first, the temperature data from the previous year was used; in the second, the model was started with the 1965 temperatures and the prediction of the previous year’s temperature used from there on. The 1966– errors in Table 1 refer to the first case. The second is not as accurate, but is not ridiculously inaccurate either; for example, the 25year prediction of the 1990 Northern Hemisphere temperature is 622 mK, not too dissimilar to the 648 mK observed.
The data, fit, and extrapolations obtained with these models are shown in Fig. 6 (Northern Hemisphere), and Fig. 7 (Southern Hemisphere). The general 25year prediction agrees reasonably well with the observed temperature in level, if not in exact details. It might be noted that the autoregressive nature of the feedback term in these models causes them to be biased slightly toward zero.
Examination of the residuals from this model shows them to be mundane and the exotic characteristics of the individual series missing.
First, the residuals show no significant departures from normality; their range is ±2.35 σ for the global series, there is no detectable skewness, and the standardized fourth moment is 2.65 for the 1854–1965 data, and similar in the validation interval. The Southern Hemisphere residuals are slightly longtailed. The extreme residual, −107 mK, is a result of the Krakatau eruption in 1883. Thus the errors of 30.4 mK and 26.9 mK obtained for the global temperature with a model having only four parameters imply that the 1990 temperature of 636 mK is about 21 SD above the 1854–1918 base level measured on the scale of the residual SD. The feedback term, αT(t − 1) in Eq. 3, reduces this to about a 9.3σ increase above what persistence would normally generate. The probability of a Gaussian random variable being 9σ from the mean is about 10^{−19}.
Second, the range of the spectrum of the residuals from this model is only 18, just slightly larger than expected for a whitenoise process. This spectrum hints at a weak 4.6year echolike term.
Third, the residuals were tested for stationarity using the quadraticinverse methods of ref. 25 and, generally, see Fig. 8, the test statistics are slightly less than their expected value, suggesting that the residuals due not contain any serious nonstationary terms. There are, nonetheless, some features of interest. The “time derivative” of the spectrum a_{1}(f)/a_{0}(f) is mostly negative, implying that the residual spectrum is decreasing as a function of time. This probably reflects nothing more than the improvements in spatial coverage and data quality that have occurred over the course of the record. The tests for both hemispheres show significant nonstationarity at periods of 4 years, perhaps an artifact of leap years in the original monthly averaging process. In addition, the Southern Hemisphere test has a significant peak in the ElNino band, and the Northern Hemisphere data are quite nonstationary at the quasibiennial oscillation frequency. At low frequencies, the test statistics from both hemispheres are well below their expected value. Thus the tests for stationarity reveal no evidence for longterm unexplained variability.
The transfer functions for this model are for solar irradiance, and similarly for C_{L}. The longterm response is given by the lowfrequency response; evaluated at f = 0 with the coefficients estimated above one obtains 0.109K/(W/m^{2}) for solar irradiance changes and 2.075 K/(2 × CO_{2}). To simplify the comparison of the different units one may examine the lowfrequency variance explained for the two components and H_{C}(0)^{2}σ_{c}^{2} = 22,536 (mK)^{2} and H_{S}(0)^{2}σ_{I}^{2} = 6,067 (mK)^{2} in the Northern Hemisphere, with corresponding figures of 28,320 (mK)^{2} and 317 (mK)^{2} in the Southern Hemisphere. Thus CO_{2} explains over 3 times as much variance in changes in solar irradiance in the Northern Hemisphere, over 100 times as much in the Southern. Using simply the change in CO_{2} concentration from 287.7 in 1854 to 352.7 parts per million in 1990 the lowfrequency transfer function predicts an increase of 0.61 K for CO_{2} while an irradiance change of 2.4 W/m^{2} (from 1365.7 to 1368.1 W/m^{2} results in an estimated temperature change of 0.26 K.
Given the great improvement that the addition of the feedback term makes on the direct leastsquares estimates, it seems reasonable to investigate related models. A suite of such models, incorporating a range of time delays and prewhitening filters was tested. None of these simple models significantly improve on Eq. 3.
Changes in the Timing of the Annual Cycle
I presented evidence (2) that the frequency of the annual temperature cycle at many locations is closer to the anomalistic year than it is to the tropical year, and, in addition, the character of the phase of the annual cycle began to change rapidly near the middle of this century with the average change in phase coherent with the changes in atmospheric concentration of CO_{2}. I also suggested that the cause of the rapid phase shift was the changing balance between direct insolation at a given latitude and transported energy. Recall that, at temperate latitudes, the periodicity of insolation is the tropical year and peaks at the summer solstice. The total radiation received by the earth, on the other hand, has the periodicity of the anomalistic year and is maximum at perihelion, currently in early January. In the Northern Hemisphere the phases of the two components are thus nearly opposed and, consequently, the phase of the resultant is easily perturbed. The Sable Island example shown in ref. 6 demonstrates that the proposed mechanism can explain the rapidly changing phase seen at individual stations. The geographic distribution of these phase changes is available in (32, 33), and confirmatory changes are directly observable in the CO_{2} record (34).
The observed phase changes must be taken seriously as an indicator of climate change caused by increasing concentrations of atmospheric CO_{2}. First, the increasing variability of the phases and the change in distribution (Fig. 4 of ref. 2) is so large that the probability of such an occurrence in a stationary process is nearly zero. Second, because the effects of observational errors and similar noiselike effects on the data will be uncorrelated between the low frequencies where the increase in temperature is observed and the annual cycle, one must consider them as independent evidence for a changing climate.
Discussion and Conclusions
More effort must be made to obtain improved CO_{2} data for the 19th and early 20th centuries if possible. Similarly, extensive effort to develop better, or independent, solar irradiance proxies is desirable.
Analysis of the lowfrequency and annual parts of the temperature records yield at least three largely independent indicators of climate change; the change in distribution of individual station phase trends about the mean, the change in average phase, and the increase in average temperature. All are unprecedented in the instrumental record.
The probability of the observed changes occurring through natural, as opposed to anthropogenic, causes appears to be exceedingly small. First, although a major ice age causes a larger temperature change than has happened so far in response to CO_{2}, the temperature increase that occurred between 1920 and 1990 would have taken more than 2,000 years even at the historically rapid rate of the last deglaciation. During deglaciation, transient warming in the North Atlantic after a Heinrich event was faster (35), but there is no evidence for a Heinrich event during the last few centuries. Second, internal climate oscillations, shifts in storm tracks, and the like obviously can change local and regional climate by much larger amounts than observed in the hemispheric averages. These regional fluctuations appear to be superimposed on the general global trend. Additionally, such internal oscillations produce warm and cool regions that interchange over decade to century time scales (32, 36), but whose effects largely cancel in hemispheric averages. Third, while there is reasonable evidence for greater climate variability during the Holocene than has been observed during the period where instrumental data are available (37, 38), there is no evidence in the statistics that a major unidentified source of natural variation is present during the instrumental record. Such a source would have to mimic, perversely, either solar irradiance changes or the changes in atmospheric CO_{2} to cause the observed temperature changes and to be mistaken for them. Similarly, while mindful of the many caveats on data quality, spatial coverage, etc. given in ref. 1, the appearance of possible leapyear artifacts at a level below 10 mK in the residuals suggests that the data cannot be as untrustworthy as is occasionally implied. The residual temperature variation remaining once the known effects of precession, solar irradiance changes, and atmospheric CO_{2} concentration are removed bound unknown effects to about 200 mK peakto peak in the hemispheric average series during the last century.
Consider the null hypothesis that the observed temperature fluctuations and atmospheric CO_{2} levels are independent: The probability that the hemispheric temperatures would fluctuate purely by chance in such a way to produce the observed coherences with CO_{2} is exceedingly low. Given that the records encompass more than a century, the probability is so low that one would not expect to see such an event by chance during the age of the earth. The probability of the observed coherence between atmospheric CO_{2} and changes in the timing of the seasons shown in figure 13 of ref. 2 without a causal connection is similarly low. Consequently one must strongly reject the hypothesis of independence between atmospheric CO_{2} and temperature. The alternative hypothesis, that increasing levels of atmospheric CO_{2} plus a slight change in solar irradiance are causally responsible for the observed changes in temperature, in contrast, results in test statistics that are ordinary in every way. Because major changes in climate as a response to human use of fossil fuels have been predicted for more than a century (39, 40), their detection can hardly be considered surprising.
From examining the data records I conclude: Changes in solar irradiance explain perhaps onequarter of the increase in temperature during the last century. The changes in atmospheric CO_{2} concentration resulting from human consumption of fossil fuels cause most of both the temperature increase and the changes in the seasonal cycle.
Footnotes

↵* Briefly, the annual temperature cycle at a given latitude consists of direct insolation plus transport effects. The direct insolation components vary with the tropical year, the time from equinox to equinox, 365.2422 days, while the net radiation received by earth and hence the mean transport vary as the anomalistic year, the time from perihelion to perihelion, or 365.2596 days. Precession, loosely defined, is the change in the longitude of perihelion measured from the vernal equinox.

↵† In a nonstationary process the spectrum S(f_{1}, f_{2}) is a function of two frequencies and several representations are useful. (Stationary processes are much simpler because the spectrum is concentrated on the line f_{1} = f_{2} and, consequently, a function of a single frequency.) A onedimensional Fourier transform of S(f_{1}, f_{2}) taking the variable f_{1} − f_{2} into t_{o} gives a dynamic spectrum D(f, t_{o}) where the frequency f = (f_{1} + f_{2})/2. The “time derivatives” of the spectrum referred to are expansion coefficients of lnD(f, t) and, approximately, of the form ∂/∂t lnD(f, t).
 Copyright © 1997, The National Academy of Sciences of the USA
References
 ↵
 Houghton J T,
 Meira Fihlo L G,
 Callender B A,
 Harris N,
 Kattenberg A,
 Maskell K
 ↵
 Thomson D J
 ↵
 Santer B D,
 Taylor K E,
 Wigley T M L,
 Penner J E,
 Jones P D,
 Cubasch U
 ↵
 ↵
 ↵
 ↵
 Sonett C P,
 Giampapa M S,
 Matthews M S
 ↵
 Peterson D H
 Keeling C D,
 Bacastow R B,
 Carter A F,
 Piper S C,
 Whorf T P,
 Heimann M,
 Mook W G,
 Roeloffzen H
 ↵
 ↵
 ↵
 Bradley R S,
 Jones P D
 ↵
Boden, T. A., Kaiser, D. P., Sepanski, R. J. & Stoss, F. W. (1994) Trends ’93: A Compendium of Data on Global Change Pub. ORNL/CDIAC65, (Oak Ridge National Laboratory, Oak Ridge, TN).
 ↵
 ↵
Marland, G., Anders, R. J. & Boden, T. A. (1994) Trends ’93: A Compendium of Data on Global Change Pub. ORNL/CDIAC65, (Oak Ridge National Laboratory, Oak Ridge, TN), pp. 505–584.
 ↵
 Foukal P,
 Lean J
 ↵
 FriisChristensen E,
 Lassen K
 ↵
 Hansen J,
 Lebedeff S
 ↵
 Draper N R,
 Smith H
 ↵
 Seber G A F
 ↵
 ↵
 Thomson D J
 ↵
 Percival D B,
 Walden A T
 ↵
 Haykin S
 Thomson D J,
 Chave A D
 ↵
 Thomson D J
 ↵
 Thomson D J
 ↵
 Thomson D J
 ↵
 Kleiner B,
 Martin R D,
 Thomson D J
 ↵
Thomson, D. J. (1994) in Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing, Adelaide, Australia 6, 73–76.
 ↵
 ↵
 Jones R H
 ↵
 Larsen J C
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 O’Brien S R,
 Mayewski P A,
 Meeker L D,
 Meese D A,
 Twickler M S,
 Whillow S I
 ↵
 Arrhenius S
 ↵