600-year global-scale temperature patterns

Global-scale temperature patterns and climate forcing over the past six centuries

Correspondence and requests for materials should be addressed to M.E.M. (e-mail: mann@snow.geo.umass.edu).

Knowing both the spatial and temporal patterns of climate change over the past several centuries remains a key to assessing a possible anthropogenic impact on post-industrial climate¹. In addition to the possibility of warming due to increased concentrations of greenhouse gases during the past century, there is evidence that both solar irradiance and explosive volcanism have played an important part in forcing climate variations over the past several centuries²,³. The unforced 'natural variability' of the climate system may also be quite important on multidecadal and century timescales⁴,⁵. If a faithful empirical description of climate variability could be obtained for the past several centuries, a more confident estimation could be made of the roles of different external forcings and internal sources of variability on past and recent climate. Because widespread instrumental climate data are available for only about one century, we must use proxy climate indicators combined with any very long instrumental records that are available to obtain such an empirical description of large-scale climate variability during past centuries. A variety of studies have sought to use a 'multiproxy' approach to understand long-term climate variations, by analysing a widely distributed set of proxy and instrumental climate indicators¹,^5-8 to yield insights into long-term global climate variations. Building on such past studies, we take a new statistical approach to reconstructing global patterns of annual temperature back to the beginning of the fifteenth century, based on the calibration of multiproxy data networks by the dominant patterns of temperature variability in the instrumental record.

We use a multiproxy network consisting of widely distributed high-quality annual-resolution proxy climate indicators, individually collected and formerly analysed by many palaeoclimate researchers (details and references are available: see Supplementary Information). The network includes (Fig. 1a) the collection of annual-resolution dendroclimatic, ice core, ice melt, and long historical records used by Bradley and Jones⁶ combined with other coral, ice core, dendroclimatic, and long instrumental records. The long instrumental records have been formed into annual mean anomalies relative to the 1902–80 reference period, and gridded onto a 5° 5° grid (yielding 11 temperature grid-point series and 12 precipitation grid-point series dating back to 1820 or earlier) similar to that shown in Fig. 1b. Certain densely sampled regional dendroclimatic data sets have been represented in the network by a smaller number of leading principal components (typically 3–11 depending on the spatial extent and size of the data set). This form of representation ensures a reasonably homogeneous spatial sampling in the multiproxy network (112 indicators back to 1820).

Potential limitations specific to each type of proxy data series must be carefully taken into account in building an appropriate network. Dating errors in a given record (for example, incorrectly assigned annual layers or rings) are particularly detrimental if mutual information is sought to describe climate patterns on a year-by-year basis. Standardization of certain biological proxy records relative to estimated growth trends, and the limits of constituent chronology segment lengths (for example, in dendroclimatic reconstructions), can restrict the maximum timescale of climate variability that is recorded⁹, and only a limited subset of the indicators in the multiproxy network may thus 'anchor in' the longest-term trends (for example, variations on timescales greater than 500 years). However, the dendroclimatic data used were carefully screened for conservative standardization and sizeable segment lengths. Moreover, the mutual information contained in a diverse and widely distributed set of independent climate indicators can more faithfully capture the consistent climate signal that is present, reducing the compromising effects of biases and weaknesses in the individual indicators.

Monthly instrumental land air and sea surface temperature¹⁰ grid-point data (Fig. 1b) from the period 1902–95 are used to calibrate the proxy data set. Although there are notable spatial gaps,this network covers significant enough portions of the globe to form reliable estimates of Northern Hemisphere mean temperature, and certain regional indices of particular importance such as the 'NINO3' eastern tropical Pacific surface temperature index often used to describe the El Niño phenomenon. The NINO3 index is constructed from the eight grid-points available within the conventional NINO3 box (5° S to 5° N, 90–150° W).

Although studies have shown that well chosen regional paleoclimate reconstructions can act as surprisingly representative surrogates for large-scale climate^11-13 multiproxy networks seem to provide the greatest opportunity for large-scale palaeoclimate reconstruction⁶ and climate signal detection¹,⁵. There is a rich tradition of multivariate statistical calibration approaches to palaeoclimate reconstruction, particularly in the field of dendroclimatology where therelative strengths and weaknesses of various approaches to multivariate calibration have been well studied¹⁴,¹⁵. Such approaches have been applied to regional dendroclimatic networks to reconstruct regional patterns of temperature¹⁶,¹⁷ and atmospheric circulation^18-20 or specific climate phenomena such as the Southern Oscillation²¹. Largely because of the inhomogeneity of the information represented by different types of indicators in a true 'multiproxy' network, we found conventional approaches (for example, canonical correlation analysis, CCA, of the proxy and instrumental data sets) to be relatively ineffective. Our approach to climate pattern reconstruction relates closely to statistical approaches which have recently been applied to the problem of filling-in sparse early instrumental climate fields, based on calibration of the sparse sub-networks against the more widespread patterns of variability that can be resolved in shorter data sets²²,²³. We first decompose the twentieth-century instrumental data into its dominant patterns of variability, and subsequently calibrate the individual climate proxy indicators against the time histories of these distinct patterns during their mutual interval of overlap.

Implicit in our approach are at least three fundamental assumptions. (1) The indicators in our multiproxy trainee network are linearly related to one or more of the instrumental training patterns. In the relatively unlikely event that a proxy indicator represents a truly local climate phenomenon which is uncorrelated with larger-scale climate variations, or represents a highly nonlinear response toclimate variations, this assumption will not be satisfied. (2) A relatively sparse but widely distributed sampling of long proxy and instrumental records may nonetheless sample most of the relatively small number of degrees of freedom in climate patterns at interannual and longer timescales. Regions not directly represented in the trainee network may nonetheless be indirectly represented through teleconnections with regions that are. The El Niño/Southern Oscillation (ENSO), for example, shows global-scale patterns of climatic influence²⁴, and is an example of a prominent pattern of variability which, if captured, can potentially describe variability in regions not directly sampled by the trainee data. (3) Patterns of variability captured by the multiproxy network have analogues in the patterns we resolve in the shorter instrumental data. This last assumption represents a fairly weak 'stationarity' requirement — we do not require that the climate itself be stationary. In fact, we expect that some sizeable trends in the climate may be resolved by our reconstructions. We do, however, assume that the fundamental spatial patterns of variation which the climate has shown during the past century are similar to those by which it has varied during past recent centuries. Studies of instrumental surface-temperature patterns suggest that such a form of stationarity holds up at least on multidecadal timescales, during the past century²³. The statistical cross-validation exercises we describe later provide the best evidence that these key underlying assumptions hold.

We isolate the dominant patterns of the instrumental surface-temperature data through principal component analysis²⁵ (PCA). PCA provides a natural smoothing of the temperature field in terms of a small number of dominant patterns of variability or 'empirical eigenvectors'. Each of these eigenvectors is associated with a characteristic spatial pattern or 'empirical orthogonal function' (EOF) and its characteristic evolution in time or 'principal component' (PC). The ranking of the eigenvectors orders the fraction of variance they describe in the (standardized) multivariate data during the calibration period. The first five of these eigenvectors describe a fraction = 0.93 (that is, 93%) of the global-mean (GLB) temperature variations, 85% of the Northern Hemisphere-mean (NH) variations, 67% of the NINO3 index, and 76% of the non-trend-related (DETR) NH variance (see Methods for a description of the statistic used here as a measure of resolved variance). A sizeable fraction of the total multivariate spatiotemporal variance (MULT) in the raw (instrumental) data (27%) is described by these five eigenvectors, or about 30% of the standardized variance (no. 1 = 12% , no. 2 = 6.5%, no. 3 = 5%, no. 4 = 4%, no. 5 = 3.5%). Figure 2 shows the EOFs of the first five eigenvectors. The associated PCs and their reconstructed counterparts (RPCs) are discussed in the next section. The first eigenvector, associated with the significant global warming trend of the past century, describes much of the variability in the global (GLB = 88%) and hemispheric (NH = 73%) means. Subsequent eigenvectors, in contrast, describe much of the spatial variability relative to the large-scale means (that is, much of the remaining MULT). The second eigenvector is the dominant ENSO-related component, describing 41% of the variance in the NINO3 index. This eigenvector shows a modest negative trend which, in the eastern tropical Pacific, describes a 'La Niña'-like cooling trend²⁶, which opposes warming in the same region associated with the global warming pattern of the first eigenvector. The third eigenvector is associated largely with interannual-to-decadal scale variability in the Atlantic basin and carries the well-known temperature signature of the North Atlantic Oscillation (NAO)²⁷ and decadal tropical Atlantic dipole²⁸. The fourth eigenvector describes a primarily multidecadal timescale variation with ENSO-scale and tropical/subtropical Atlantic features, while the fifth eigenvector is dominated by multidecadal variability in the entire Atlantic basin and neighbouring regions that has been widely noted elsewhere^29-34.

We calibrate each of the indicators in the multiproxy data network against these empirical eigenvectors at annual-mean resolution during the 1902–80 training interval. Although the seasonality of variability is potentially important — many extratropical proxy indicators, for example, reflect primarily warm-season variability⁶,⁷ — we seek in the present study to resolve only annual-mean conditions, exploiting the seasonal climate persistence, and the fact that the mutual information from data reflecting various seasonal windows should provide complementary information regarding annual mean climate conditions¹⁰. Following this calibration, we apply an overdetermined optimization procedure to determine the best combination of eigenvectors represented by the multiproxy network back in time on a year-by-year basis, with a spatial coverage dictated only by the spatial extent of the instrumental training data. From the RPCs, spatial patterns and all relevant averages or indices can be readily determined. The details of the entire statistical approach are described in the Methods section.

(1) In the reconstructions from 1820 onwards based on the full multiproxy network of 112 indicators, 11 eigenvectors are skilfully resolved (nos 1–5, 7, 9, 11, 14–16) describing 70–80% of the variance in NH and GLB mean series in both calibration and verification. (Verification is based here on the independent 1854–1901 data set which was withheld; see Methods.) Figure 3 shows the spatial patterns of calibration and verification and the squared correlation statistic r², demonstrating highly significant reconstructive skill over widespread regions of the reconstructed spatial domain. 30% of the full spatiotemporal variance in the gridded data set is captured in calibration, and 22% of the variance is verified in cross-validation. Some of the degradation in the verification score relative to the calibration score may reflect the decrease in instrumental data quality in many regions before the twentieth century rather than a true decrease in resolved variance. These scores thus compare favourably to the 40% total spatiotemporal variance that is described by simply filtering the raw 1902–80 instrumental data with 11 eigenvectors used in calibration, suggesting that the multiproxy calibrations are describing a level of variance in the data reasonably close to the optimal 'target' value. Although a verification NINO3 index is not available from 1854 to 1901, correlation ofthe reconstructed NINO3 index with the available Southern Oscillation index (SOI) data from 1865 to 1901 of r = -0.38 (r² = 0.14) compares reasonably with its target value given by the correlation between the actual instrumental NINO3 and SOI index from 1902 to 1980 (r = -0.72 ). Furthermore, the correspondence between the reconstructed NINO3 index warm events and historical³⁵ El Niño chronology back to 1820 (see Methods) is significant at the 98% level.

(3) The network available back to 1700 of 74 indicators (including only two instrumental or historical indicators) skilfully resolves five eigenvectors (nos 1, 2, 5, 11, 15) and shows some significant signs of decrease in reconstructive skill. In this case, 60–70% of NH variance is resolved in calibration and verification, 14–18% of MULT in calibration, and 10–12% of MULT in verification. The verification r of NINO3 with the SOI is in the range of r -0.25 to -0.35, which is statistically significant (as is the correspondence with the historical³⁵ chronology back to 1700) but notably inferior to the later calibrations. In short, both spatial patterns and large-scale means are skilfully resolved, but with significantly less resolved variance than in later calibrations.

The reconstructions discussed here are derived using all indicators available, and using the optimal eigenvector subsets determined in the calibration experiments described above (11 from 1780–1980, 9 from 1760–1779, 8 from 1750–1759, 5 from 1700–1749, 4 from 1600–1699, 2 from 1450–1599, 1 from 1400–1449). To better illustrate the workings and effectiveness of the proxy pattern reconstruction procedure, we show as an example (Fig. 4) the actual, the EOF-filtered, and the reconstructed temperature patterns for a year (1941) during the calibration interval. This year was a known ENSO year, associated with a warm eastern tropical Pacific and a cold central North Pacific. Pronounced cold anomalies were also found over large parts of Eurasia. The proxy-reconstructed pattern captures these features, although in a relatively smoothed sense (describing 30% of the full variance in that pattern), and is remarkably similar to the raw data once it has been filtered by retaining only the 11 eigenvectors (nos 1–5, 7, 9, 11, 14–16) used in pattern reconstruction. It is thus visually apparent that the multiproxy network is quite capable of resolving much of the structure resolved by the eigenvectors retained in the calibration process.

We consider the temporal variations in the first five RPCs ( Fig. 5a). The positive trend in RPC no. 1 during the twentieth century is clearly exceptional in the context of the long-term variability in the associated eigenvector, and indeed describes much of the unprecedented warming trend evident in the NH reconstruction. The negative trend in RPC no. 2 during the past century is also anomalous in the context of the longer-term evolution of the associated eigenvector. The recent negative trend is associated with a pattern of cooling in the eastern tropical Pacific (superimposed on warming associated with the pattern of eigenvector no.1) which may be a modulating negative feedback on global warming²⁶. RPC no. 5 shows notable multidecadal variability throughout both the modern and pre-calibration interval, associated with the wavelike trend of warming and subsequent cooling of the North Atlantic this century discussed earlier^29-34 and the longer-term multidecadal oscillations in this region detected in a previous analysis of proxy climate networks⁵. This variability may be associated with ocean–atmosphere processes related to the North Atlantic thermohaline circulation⁴,³⁴.

The long-term trends in the reconstructed annual mean NH series (Fig. 5b) are quite similar to those of decadal Northern Hemisphere summer temperature reconstructions⁶, showing pronounced cold periods during the mid-seventeenth and nineteenth centuries, and somewhat warmer intervals during the mid-sixteenth and late eighteenth centuries, with almost all years before the twentieth century well below the twentieth-century climatological mean. Taking into account the uncertainties in our NH reconstruction (see Methods), it appears that the years 1990, 1995 and now 1997 (this value recently calculated and not shown) each show anomalies that are greater than any other year back to 1400 at 3 standard errors, or roughly a 99.7% level of certainty. We note that hemispheric mean values are not associated with globally or hemispherically uniform trends. An example of the global pattern for anhistorically documented³⁵ "very strong" El Niño year (1791) is shown in Fig. 6 top panel, demonstrating the classic warm eastern tropical Pacific and cold central North Pacific sea surface temperature patterns. Analysis of ENSO variability in these reconstructions is discussed in more detail elsewhere³⁶. We also show the reconstructed pattern for 1816 (Fig. 6 bottom panel). Quite anomalous cold is evident throughout much of the Northern Hemisphere (even relative to this generally cold decade) but with a quadrupole pattern of warmth near Newfoundland and the Near East, and enhanced cold in the eastern United States and Europe consistent with the anomalous atmospheric circulation associated with the NAO pattern. Such a pattern is indeed observed in empirical³⁷ and model-based studies³⁸ of the atmospheric response to volcanic forcing. We infer in the 1816 temperature pattern a climatic response to the explosive Tambora eruption of April 1815 based on both the anomalous hemispheric coolness and the superimposed NAO-like pattern. Reconstructed time series RPCs nos 1–5, the NH series, the NINO3 index and reconstructions for specific grid-points can be obtained through the NOAA palaeoclimatology Web site (http://www.ngdc.noaa.gov/paleo/paleo.html ).

We take an empirical approach to detecting the possible effects of external forcings on the climate. The reconstructed NH series is taken as a diagnostic of the global climate, and we examine its relationship with three candidate external forcings during the period 1610–1995 including (1) CO ₂ measurements³⁹ as a proxy for total greenhouse-gas changes, (2) reconstructed solar irradiance variations² and (3) the weighted historical 'dust veil index' (DVI) of explosive volcanism (see Fig. 31.1 in ref. 40) updated with recent data⁴¹. While we warn that historical series for these forcing agents are imperfectly known or measured, they do nonetheless represent our best estimates of the time-histories of the corresponding forcings. More detailed discussions of the estimation of, and potential sources of uncertainty or bias in, these series are available²,³⁹,⁴⁰. Industrial-aerosol forcing of the climate has also been suggested as an important forcing of recent climate⁴²,⁴³, but its physical basis is still controversial⁴⁴, and difficult to estimate observationally. Noting that in any case, this forcing is not believed to be important before about 1940, its omission should be inconsequential in our long-term detection approach. Our empirical signal detection is complementary to that of model-based "fingerprint" signal detection studies⁴²,⁴³,⁴⁵; while our empirical approach relies on the faithfulness of the reconstructed forcing series and on the assumption of a linear and contemporaneous response to forcings, it does not suffer the potential weaknesses of incomplete representations of internal feedback processes²⁶, poorly constrained parametrizations of climatic responses⁴⁴, and underestimated natural variability¹ in model-based studies. To the extent that the response to forcing is not contemporaneous, but rather is delayed owing to the inertia of the slow-response components of the climate system (for example, the ocean and cryosphere), our detection approach will tend to underestimate the response to forcings, making the approach a relatively conservative one.

We estimate the response of the climate to the three forcings based on an evolving multivariate regression method (Fig. 7). This time-dependent correlation approach generalizes onprevious studies of (fixed) correlations between long-term Northern Hemisphere temperature records and possible forcing agents²,³. Normalized regression (that is, correlation) coefficients r are simultaneously estimated between each of the three forcing series and the NH series from 1610 to 1995 in a 200-year moving window. The first calculated value centred at 1710 is based on data from 1610 to 1809, and the last value, centred at 1895, is based on data from 1796 to 1995 — that is, the most recent 200 years. A window width of 200 yr was chosen to ensure that any given window contains enough samples to provide good signal-to-noise ratios in correlation estimates. Nonetheless, all of the important conclusions drawn below are robust to choosing other reasonable (for example, 100-year) window widths.

We test the significance of the correlation coefficients (r) relative to a null hypothesis of random correlation arising from natural climate variability, taking into account the reduced degrees of freedom in the correlations owing to substantial trends and low-frequency variability in the NH series. The reduced degrees of freedom are modelled in terms of first-order markovian 'red noise' correlation structure of the data series, described by the lag-one autocorrelation coefficient during a 200-year window. This parameter ranges from 0.48 in the first window (1610–1809) to 0.77 in the final window (1796–1995) of the moving correlation, the considerably larger recent value associated with the substantial global warming trend of the past century. This latter trend has been shown to be inconsistent with red noise⁴⁶ and could thus itself be argued as indicative of externally forced variability. An argument could in this sense, be made for using the smaller pre-industrial value = 0.48 of the NH series in estimating the statistical degrees of freedom appropriate for the null hypothesis of natural variability. Nonetheless, we make the conservative choice of adopting the largest value = 0.77 as representative of the natural serial correlation in the series. We use Monte Carlo simulations to estimate the likelihood of chance spurious correlations of such serially correlated noise with each of the three actual forcing series. The associated confidence limits are approximately constant between sliding 200-year windows. For (positive) correlations with both CO₂ and solar irradiance, the confidence levels are both approximately 0.24 (90%), 0.31 (95%), 0.41 (99%), while for the 'whiter', relatively trendless, DVI index, the confidence levels for (negative) correlations are somewhat lower (-0.16, -0.20, -0.27 respectively). A one-sided significance test is used in each case because the physical nature of the forcing dictates a unique expected sign to the correlations (positive for CO₂ and solar irradiance variations, negative for the DVI fluctuations).

The correlation statistics indicate highly significant detection of solar irradiance forcing in the NH series during the 'Maunder Minimum' of solar activity from the mid-seventeenth to early eighteenth century which corresponds to an especially cold period. In turn, the steady increase in solar irradiance from the early nineteenth century through to the mid-twentieth century coincides with the general warming over the period, showing peak correlation during the mid-nineteenth century. The regression against solar irradiance indicates a sensitivity to changes in the 'solar constant' of 0.1 K W ^-1 m^-2, which is consistent with recent model-based studies⁴². Greenhouse forcing, on the other hand, shows no sign of significance until a large positive correlation sharply emerges as the moving window slides into the twentieth century. The partial correlation with CO₂ indeed dominates over that of solar irradiance for the most recent 200-year interval, as increases in temperature and CO ₂ simultaneously accelerate through to the end of 1995, while solar irradiance levels off after the mid-twentieth century. It is reasonable to infer that greenhouse-gas forcing is now the dominant external forcing of the climate system. Explosive volcanism exhibits the expected marginally significant negative correlation with temperature during much of 1610–1995 period, most pronounced in the 200-year window centred near 1830 which includes the most explosive volcanic events.

A variety of general circulation⁴²,⁴⁷ and energy-balance model experiments⁴³,⁴⁵,⁴⁸ as well as statistical comparisons of twentieth-century global temperatures with forcing series⁴⁹ suggest that, although both solar and greenhouse-gas forcings play some role in explaining twentieth-century climate trends, greenhouse gases appear to play an increasingly dominant role during this century. Such a proposition is consistent with the results of this study.

As larger numbers of high-quality proxy reconstructions become available in diverse regions of the globe, it may be possible to assimilate a more globally representative multiproxy data network. Given the high level of skill possible in large-scale reconstruction back to 1400 with the present network, it is reasonable to hope that it may soon be possible to faithfully reconstruct mean global temperatures back over the entire millennium, resolving for example the enigmatic⁷ medieval period. Geothermal measurements from boreholes⁵⁰ recover long-term temperature trends without many of the complications of traditional proxy indicators and, in combination with traditional multiproxy networks, may prove helpful in better resolving trends over many centuries. With a better knowledge of how the climate has varied before the twentieth century, we will be able to place even better constraints on the importance of natural and anthropogenic factors governing the climate of the past few centuries, factors which will no doubt continue to affect climate variability in the future, in addition to any anthropogenic effects.

Supplementary information is available on Nature's World-Wide Web site (http://www.nature.com) or as paper copy from the London editorial office of Nature.

Received 9 May 1997; accepted 27 February 1998