1 Introduction

Climate models provide both the means to explore potential future climatic changes, and to test hypotheses concerning the causes of past climate. Although the former application has received the most publicity, the latter has the potential to provide challenging tests of modelling capability by providing a range of cases in which climate boundary conditions have varied well outside the limits of contemporary (or historically recent) observations (for further discussion on model benchmarking using palaeoclimate simulations see e.g. Braconnot et al. 2007a, b; Edwards et al. 2007; Schmidt 2010).

The Last Glacial Maximum (LGM, ca. 21,000 years ago, 21 ka) and middle Holocene (MH, ca. 6,000 years ago, 6 ka) have been the main target periods for the Palaeoclimate Modelling Intercomparison Project (PMIP http://pmip2.lsce.ipsl.fr/) and continue to be important standard benchmarks for model evaluation in the current phase of this project (Otto-Bliesner et al. 2009). Both time periods are included in the suite of simulations being carried out as part of the fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change, and thus will contribute to the benchmarking of climate models used to predict potential future climate changes. At the LGM, large ice sheets covered substantial areas of the northern continents; the major greenhouse gases (carbon dioxide, methane and nitrous oxide) were at historic minimum atmospheric concentrations; and terrestrial ecosystems were less productive, of different composition, and generally shifted toward the equator. The MH was much more similar to the present than the LGM, and greenhouse gases were at levels close to those before the industrial revolution. However, the Earth’s orbital configuration differed from that of today, with perihelion in boreal summer/autumn (implying greater seasonality of insolation in the northern hemisphere) and greater obliquity, implying higher summer (and annual) insolation in high latitudes. Although neither period is a strict analogue or “anti-analogue” for future conditions, the LGM is of considerable interest as a period with strongly negative radiative forcing (including low carbon dioxide concentration) while the MH offers the possibility of studying a period with warmer conditions and longer growing seasons in many regions, comparable in some respects to climates expected as a result of increasing radiative forcing in the future.

The development of well-documented, quantitative global palaeodata sets portraying the spatial climatic patterns of the LGM and MH time periods is a central objective of the PMIP3 research programme (Otto-Bliesner et al. 2009). Terrestrial environments are diverse and many sources of palaeodata provide information on specific areas and ecosystems. Palaeovegetation (fossil pollen and plant macrofossil) data provide a unique combination of near-global coverage, multivariate responses (i.e. vegetation composition carries information about several distinct aspects of climate), and well-documented methods of providing quantitative palaeoclimatic estimates. This paper presents a global synthesis of reconstructed terrestrial climatic estimates for the LGM and MH, based on the application of common methods.

The palaeoceanographic community has recently completed a global reconstruction of LGM sea surface temperatures (MARGO Project Members 2009), which supersedes the widely used CLIMAP data set developed in the 1980s. Terrestrial environments are at least as diverse as marine environments, and as in the marine realm, many methods and data sources have been developed in order to generate palaeoclimate reconstructions that can be compared with modelling results. However, the most recent synthesis of quantitative climate reconstructions based on terrestrial palaeodata (COHMAP Members 1988; Wright et al. 1993) fell short of generating a unified global data set. Here we present such a data set, for the first time.

2 Methods

The MH and LGM reconstructions presented here were produced by combining existing quantitative reconstructions based on pollen and plant macrofossils that have been made for different regions of the world and using a variety of different methods. Although a few climate reconstructions have been made using other types of palaeoenvironmental records, such as chironomid or coleopteran assemblages (e.g. Atkinson et al. 1987; Heiri et al. 2003), only pollen and plant macrofossil data are extensive enough to provide quantitative palaeoclimate estimates at continental-to-regional scales. The procedure for selection, screening and combining the existing data sets is discussed in Sect. 2.6 after discussion of the general principles and specifics of the different quantitative reconstruction techniques.

2.1 Principles

Geographic distributions of plant species can be approximated by “envelopes” representing conterminous regions in climate space. The climate space required to describe species distributions is not a 1D space (i.e. species distributions cannot be represented accurately along a single axis such as mean annual temperature) but rather, at least a 2D (Iversen 1944; Hintikka 1963), 3D (Thompson et al. 1999) or 4D (Bartlein et al. 1998; Williams et al. 2001) space, with dimensions related to temperature seasonality and water availability. This situation is understood to arise because of the existence of bioclimatic limits on the distribution of taxa (Woodward 1987; Harrison et al. 2009). If these limits acted separately, they would define a volume in climate space, where the axes are provided by variables that relate to the relevant mechanisms (e.g. growing-degree days, GDD, rather than mean July temperature as an index of accumulated warmth during the growing-season, etc.) This assumption underpins the way in which the limits of population viability for taxa (Sykes et al. 1996; Smith et al. 2001) or plant functional types (PFTs, e.g. Sitch et al. 2003; Gerten et al. 2004; Prentice et al. 2007) are specified in regional and global modelling of vegetation dynamics.

Fossil pollen records provide evidence of changing plant distributions through time (e.g. Prentice et al. 2000; Bigelow et al. 2003; Pickett et al. 2004; Marchant et al. 2009). Pollen are often identifiable to genus and family level, but only rarely to species level. These higher taxa tend to have distributions as coherent as lower taxa or species in climate space, even when species from different continents are included: this was first noted by Huntley et al. (1989) for Fagus, and has been confirmed by more recent work which has established the general tendency for lineages to be evolutionarily conservative in terms of climate preferences (Ackerly 2003; Moles et al. 2005; Wake et al. 2009). Bioclimatic limits determine the fundamental niche of a taxon, but taxa (whether at species or higher level) commonly have overlapping fundamental niches. They also vary in abundance, with maximum values either near physiological optima, or displaced by competition. We therefore find more-or-less continuous variation in taxon abundance along environmental gradients in space, and we can construct relationships between composition and location in climate space (Jackson and Williams 2004).

In time, continuous changes in composition are observed. These can be explained by continuous changes in climate. The observed response in pollen assemblages is due in part to quantitative compositional shifts, and in part to range boundary shifts. When climate change is rapid, so is compositional change. The composition of pollen assemblages tracks the climate with lags of the order of 100 years or less (Gajewski 1987; Prentice et al. 1991; Birks and Ammann 2000; Davis and Shaw 2001; Tinner and Lotter 2001; Williams et al. 2002; Wehrli et al. 2007; Pross et al. 2009). Thus, climate changes (on multi-centennial to glacial-interglacial time scales) can be reconstructed by applying the modern spatial relationships between climate and plant assemblages to the “downcore” temporal variations in pollen, ignoring plant succession and persistence effects.

This substitution is usually done by using large surface sample data sets as “training sets” to establish empirical relationships between taxon abundances and climate variables, which can then be applied to fossil samples. It is also possible to start from contemporary species distribution directly (Kühl et al. 2002), which allows surface samples to be reserved for testing. Composition can also be predicted (most successfully in terms of plant functional types, PFTs) using process-based modelling (see e.g. Lischke et al. 2002; Koca et al. 2006; Miller et al. 2008), which is independent of either surface samples or species distributions.

2.2 Main categories of methods

2.2.1 Statistical methods

In inverse regression (Huntley and Prentice 1988; Bartlein et al. 1984), each climate variable is treated as a dependent variable in a multiple regression with the taxon abundances as predictors. (The regression is said to be “inverse”, in the sense that in nature climate is the control and taxon abundance is the response.) Improved fit can be obtained by Box-Cox transformation of the taxon abundances (Bartlein et al. 1984). The regression equations are developed on the surface sample data set, and then applied to the fossil samples directly.

Weighted averaging methods (ter Braak and Juggins 1993; Birks 1995) are a more recent type of regression, better suited to data (like pollen data, when large regions of climate space are considered) where the taxa display optima around which they decline in abundance. Optima are estimated from modern climate-abundance relationships, and estimates for fossil assemblages are made as a weighted average of the optima of all taxa found.

Artificial neural networks are widely used to calibrate the non-linear relationships between two or more parameters, and have been adapted to pollen-based climate reconstruction (e.g. by Peyron et al. 1998). A network is established as a series of coefficients or nodes that are adjusted during calibration to provide the closest possible fit, and then climate variables may be estimated from any pollen assemblage.

The modern analogue technique (Overpeck et al. 1985; Guiot 1990; Jackson and Williams 2004) searches through the surface sample data set to find the samples that give the lowest values of a dissimilarity measure, such as squared chord distance, between the fossil pollen sample and the modern ones. The method then attributes the average climate associated with these “modern analogues” to the fossil pollen sample.

The response surface technique (Bartlein et al. 1986) is conceptually related to the modern analogue technique. A response surface is a smooth surface fitted to each taxon’s abundance in climate space. The first applications fitted quadratic surfaces (paraboloids); later applications have used a more flexible surface-fitting procedure that allows more complicated shapes, and both procedures allow for mild extrapolation of the fitted surfaces. The technique seeks the point in climate space where the fossil-pollen sample shows lowest dissimilarity to the assemblage constructed from the fitted values of all of the taxa. The climate at this point is attributed to the fossil pollen sample. Waelbroeck et al. (1998) have extended this by combining response surfaces with a modern analogue approach to avoid excessive smoothing of the taxon abundances, and Gonzales et al. (2009) describe an approach for extrapolating the response surfaces by fitting symmetrical unimodal distributions to the taxon abundances in order to deal with the situation in which fossil-pollen samples have no analogues in the modern data set.

Inversion of species envelope models provides an alternative approach conceptually related to the response surface technique. Instead of using surface sample data as the training set, this approach starts with species range maps or biogeographic atlas data, usually recording species presence. Thompson et al. (2008) have used the modern analogue technique to select points from a continent-wide grid of species presence/absence distributions and associated climate that are the most similar to a fossil assemblage, assigning the associated climate as the reconstruction. The assignment may be weighted by an index of similarity such as the Jaccard coefficient. In the mutual climatic range approach, the climatic ranges for each species found in a fossil assemblage are superimposed, and the estimate is given as the climate interval common to all species (Atkinson et al. 1987). The reconstructed interval may be refined by using the 5–95 percentile range of the climate envelope to preferentially weight the estimated climate (Thompson et al. 2008) or, more rigorously, a probability density function (pdf) may be fitted to the species presence data (Kühl et al. 2002). The first applications of the pdf method fitted Gaussian surfaces; later applications have used a mixture model to obtain a more generalized form (e.g. Gavin and Hu 2006). As with response surfaces, the pdfs are fitted separately for each taxon, and fossil pollen samples are assigned the “most probable climate” to have produced the sample, based on multiplying together the pdfs of the constituent taxa.

2.2.2 Process-based methods

State-of-the-art biogeography models (Haxeltine and Prentice 1996; Kaplan et al. 2003) predict the geographic distribution and quantitative properties such as net primary production (NPP) and leaf area index (LAI) of plant functional types (PFTs), as a function of climate, soil properties (which affect water retention) and atmospheric CO2 content (which affects photosynthesis). When such models are “inverted” for climate reconstruction (Guiot et al. 2000), pollen taxa are assigned to PFTs using standard tables developed for the BIOME 6000 project and subsequent work in vegetation reconstruction (Prentice et al. 2000; Pickett et al. 2004; Marchant et al. 2009), or using independently generated bioclimatic affinity groups (Laurent et al. 2004). A climate space is constructed by systematically perturbing the input variables of the model. The method searches for the point in this climate space providing the best match between the PFT profile of the fossil pollen sample (in terms of abundances) and the PFT profiles generated by the model (in terms of NPP or LAI). A Monte Carlo Markov Chain algorithm has been used to accelerate the search and thus make the approach computationally tractable.

2.3 Advantages and problems of different methods

Inverse regression has difficulty in accommodating unimodal (as opposed to monotonic) taxon responses. As a result, it is necessary in practice to develop different regression equations applying to different regions of climate space (Bartlein and Webb 1985; Huntley and Prentice 1988). This leads to a degree of arbitrariness in the reconstruction because the best-fitting equation for a given point may not be the most suitable when applied to the same point in the past. It is encouraging that major features of early reconstructions using inverse regression have been upheld by newer techniques. But the method proved too cumbersome, and was abandoned.

Weighted averaging methods are designed to deal with unimodal responses and thus they deal with the most problematic features of conventional inverse regression. However, they have difficulty in handling more complex (e.g. bimodal) species responses. And in common with all regression methods they will always give an answer. This could be the wrong answer in no-analogue situations, i.e. when the pollen samples considered fall outside the compositional range represented by the training set.

Artificial neural networks adapt well to non-linear problems, such as vegetation climate relationships, but require a number of parameters to be chosen, including the number of nodes and the number of training iterations. Further, this remains a black-box approach in which the coefficients obtained are difficult to interpret.

Modern analogue and response surface methods can deal with arbitrarily complex niche shapes. Although conceptually simple, the application of the modern analogue technique requires several decisions that can affect the palaeoclimate estimates. These include decisions about the number of pollen taxa used to compute the squared chord distance, the threshold value of the squared chord distance that distinguishes a “good analogue” from a “no-analogue” fossil pollen sample, and the number of analogues averaged to compute the fossil palaeoclimate estimate (see e.g. Jackson and Williams 2004; Sawada et al. 2004; Viau et al. 2006, 2008; Williams and Shuman 2008). Although methods have been developed to identify “no-analogue” samples, e.g. comparison of between and within-vegetation zone dissimilarity values, the solutions depend on the vegetation classification scheme and taxonomic diversity (e.g. Overpeck et al. 1985; Anderson et al. 1989; Waelbroeck et al. 1998; Sawada et al. 2004; Williams and Shuman 2008). Assemblages dominated by widespread taxa may generate probably unrealistic “jumps” as analogues may be selected from geographically and climatically disparate regions (Brewer et al. 2008). These jumps can occur when a pollen taxon includes several species with highly distinct geographic and climatic distributions (Williams and Shuman 2008), and are also a problem when the modern calibration datasets are sparse. This feature may have led to misleading artifacts when the methods have been applied to continuous time series. No-analogue situations can be identified (by high dissimilarity values), but in these cases no reconstruction should be made; or at least the reconstruction has to have a caveat attached to it (Cheddadi et al. 1997). Furthermore, “wrong-analogue” situations can occur, and go undetected. These are situations where similar pollen assemblages occur across a wide range of a climate variable. The methods can then be biased by the available surface sample data, because it is possible that a close analogue will be found within the range of the surface samples even if the true palaeoclimate lies outside the range. Alternatively, wrong analogues can occur if the relative abundances of species represented by a particular pollen morphotype change substantially over time (Jackson and Williams 2004).

Inversion of species envelope models represents a partial solution to the no-analogue problem in that it uses species distributions only, rather than abundances, to produce the pdfs. No-analogue situations (in terms of species presence) will still produce no estimate, and wrong-analogue situations could still go undetected.

Inversion of process-based biogeographic models should, in principle, solve both the no-analogue and wrong-analogue problems by (a) using process-based logic to predict no-analogue assemblages, and (b) providing realistic error bars based on the full range of possible climates. A further strength of this approach is that it provides a way to include effects of CO2 concentration on competition between PFTs and therefore also the relative abundances of taxa (Harrison and Prentice 2003; Prentice and Harrison 2009). All other techniques are likely to give biased results for CO2 concentrations unlike those of the Late Holocene. An inevitable weakness however is that the method is model-dependent, and in practice, the forward-simulation requires the generation of potentially non-unique values of a large number of input variables (i.e., more than one reconstructed climate could correctly simulate a particular assemblage).

2.4 Including multiple data sources

Modern analogue methods can be easily adapted to apply constraints derived from another data source, e.g. lake level changes (Cheddadi et al. 1997). Inversion of biogeographic models can also be adapted in a natural way because process-based models predict ancillary variables such as runoff and plant δ13C (Hatté and Guiot 2005), the latter being a powerful additional indicator of C3 plant water status and C3/C4 plant competition (Harrison and Prentice 2003; Prentice and Harrison 2009). Hatté et al. (2009) showed that the use of δ13C as an additional constraint considerably reduces the uncertainties of climatic reconstructions obtained from inverse modelling compared to using pollen alone. Guiot et al. (2009) have indicated that inverse modelling could provide an elegant solution of combining multiple data types, including lake-levels, δ13C and pollen.

2.5 Comparative studies

Although these techniques have been employed in palaeoclimatology for a number of years, there are only a few studies in which different methods are used at the same site (see e.g. Bartlein and Whitlock 1993). In a recent study of the Eemian in Europe, Brewer et al. (2008) reconstructed climate for a set of sites using seven methods, include many of those described above. The methods were largely coherent during the interglacial period, characterised by quite diverse and information-rich assemblages. In glacial periods, however, there is increasing divergence between the values estimated by the different methods, and analogue-based methods show a marked increase in the variability of reconstructed values. These periods are characterised by information-poor assemblages that are often dominated by Pinus and Artemisia, both widespread taxa found in a range of environments. This situation poses problems for many of the methods: for example it becomes harder to obtain either true analogues or precise estimate from regression techniques. The level of agreement or divergence between different techniques does however provide a measure of the robustness of the reconstruction.

2.6 Methods used for the synthesis

We provide gridded reconstructions of climate variables for the MH (defined here as the period around 6,000 calendar years BP, 6 ka) and the LGM (defined here as the period around 21,000 calendar years BP, 21 ka) based on combining all existing quantitative reconstructions, subject to availability of the primary data (i.e. the reconstructions) and a transparent screening procedure. Quantitative estimates of past terrestrial climates were provided by the authors, from online public-access datasets, or from the literature.

Regional- to continental-scale palaeoclimate reconstructions have been made for Europe, Eurasia, North America and Africa (Table 1) using extensive archive pollen databases (the Global Pollen Database http://www.ncdc.noaa.gov/paleo/gpd.html, the North American Pollen Database http://www.ncdc.noaa.gov/paleo/napd.html, NEOTOMA http://www.neotomadb.org/, the European Pollen Database http://www.europeanpollendatabase.net/, the African Pollen Database http://medias3.mediasfrance.org/apd/) and modern calibration datasets (e.g. Gajewski et al. 2002; Whitmore et al. 2005). No large-scale regional climate reconstructions have been made for South America and Australasia, despite the existence of archive pollen databases for both regions (see http://www.ncdc.noaa.gov/paleo/lapd.html, http://palaeoworks.anu.edu.au/databases.html).

Table 1 Characteristics of existing regional pollen-based reconstruction data sets
Full size table

Several different reconstruction techniques have been applied to data sets extracted from the archival databases (Table 1; Fig. 1), ranging from simple regression through analogue techniques to inverse modelling methods. We made no attempt to select a single data set for a given region––in part because it would be difficult to evaluate regional patterns given that the maps would be based on different methods for different regions, and in part because (as will be demonstrated later) different techniques appear to yield remarkably similar results at least to first order.

Fig. 1
figure 1

Distribution of sites included in the benchmark data sets for 6 ka (upper panels) and 21 ka (middle panels) in geographic space by reconstruction method. Reconstructions made with different methods are available for some sites. The sites (shown as black dots) are also plotted in climate and vegetation space (lower panel), where climate space is represented by accumulated warmth during the growing season (GDD0) and the ratio of actual to potential evapotranspiration (AE/PE) and vegetation space by the distribution of major biomes within this climate space

Full size image

The data coverage for Europe and North America is considerably denser than that for Russia, Asia and Africa. For regions with no large-scale reconstructions, we searched the literature for individual studies where researchers made time-series reconstructions for one or more individual cores based on an explicit statistical calibration method and a regional calibration dataset. These were retained only in regions from where we had no other reconstructions (Table 2).

Table 2 Characteristics of existing individual pollen-based reconstructions
Full size table

We have used climate reconstructions based on pollen samples assigned by the original authors to 6 or 21 ka. In those cases where the original authors provide reconstructions for all pollen samples in a sequence, we have used the climate reconstructions based on the pollen sample closest to 6 ka, provided it fell within the range 5.5–6.5 ka. Some early datasets used the same convention in radiocarbon years, equivalent to 6.3–7.45 ka. As it was difficult to determine the age of the sample in order to calibrate them, we instead retained the original non-calibrated sample. In the worse case, this means the samples used could range in age from 5.5 to 7.45 ka. However the bulk of the samples (79%) definitely fall within the range of 5.5–6.5 ka, and only 3.5% of the samples definitely fall outside this range. A similar situation pertains to 21 ka, where we have taken the closest sample to 21 cal ka (±1 ka), but earlier studies would have been based on samples around 18,000 14C year BP, i.e. ranging from 20.25 to 22.9 cal year BP. In the worse case, this means the samples could range from 20 to 22.9 ka. Again, the majority of the samples (63%) definitely fall in the range 21 ± 1 cal ka.

Regional scale reconstructions have been mapped as values at individual sites (e.g. Jackson et al. 2000; Wu et al. 2007) or in some cases gridded (e.g. Cheddadi et al. 1997; Webb et al. 1998). Although in some cases only gridded fields were published, we obtained the original point data for this study. However, in one case (Webb et al. 1998), the reconstructions were done after the pollen data were gridded and thus point data reconstructions were not available. Gridding data involves interpolation and loss of the original data values, and in some cases extrapolation beyond the limits of the data (see discussion in Sawada et al. 2004). We have not used the gridded datasets in our analyses because we wanted a dataset that could subsequently be gridded at the appropriate scale for analysis and comparison with climate model output.

The datasets we have used in this study include reconstructions of many different climatic variables. However, in any one study, only a subset of these were typically reconstructed. Many climate variables are highly correlated; there are insufficient degrees of freedom to reconstruct all of the potential climate variables of interest. We retained those that are considered functionally important to plant growth (Harrison et al. 2009), generally called bioclimatic variables. Some climatic variables are functionally equivalent to these bioclimatic variables (e.g. July temperature in the northern hemisphere is functionally equivalent to mean temperature of the warmest month, MTWA: Harrison et al. 2009). So, in the northern hemisphere we have amalgamated July temperature (December temperature in the southern hemisphere) with reconstructions of MTWA. Similarly, we have amalgamated December temperature in the northern hemisphere (July temperature in the southern hemisphere) with its functional equivalent, mean temperature of the coldest month (MTCO). Six variables are included in the final analysis: mean temperature of the warmest month (MTWA), mean temperature of the coldest month (MTCO), accumulated temperature during the growing season (an index of growing season warmth) measure by growing degree days above a base level of 5°C (GDD5), mean annual temperature (MAT), mean annual precipitation (MAP) and a moisture index (alpha, α) defined as the ratio of actual to potential evaporation (Table 3).

Table 3 Reconstruction variables used in this study and those considered functionally equivalent
Full size table

We estimated climate anomalies (i.e. the difference between the target period and the present day, e.g. 6–0 ka) at each site. We do this because palaeoclimate reconstructions can differ from the present due to changes in the climate but also due to differences in the calibration set used as a baseline. The use of anomalies avoids this source of error. Wherever possible, we used the original climate reconstructions for 6 and 0 ka to estimate the climate anomaly. In the cases where modern values were not included in the database, we extracted the modern climate values from the Climate Research Unit (CRU) historical climatology data set (CRU CL 2.0, New et al. 2002).

Different studies have expressed reconstruction uncertainties in different ways. In this study, uncertainties are re-expressed as standard error of estimates where there was sufficient information to calculate these. Where the authors did not report uncertainty, or where it was impossible to compute it from available data, we have used the average global uncertainty (i.e. the RMSE of regression) for each bioclimatic variable derived from the regressions of observed versus reconstructed bioclimate variables shown in Fig. 2.

Fig. 2
figure 2

The number of sites in each grid cell contributing to the benchmark reconstructions of individual climatic variables, after the screening procedures, at 6 and 21 ka

Full size image

There are a few instances where the reconstructions of modern climate variables are very different from the observed climate (i.e. CRU). These were most often sites in topographically diverse regions. We used the Mahalanobis Distance to identify points which differed significantly from the observations in an objective way (Bartlein et al. 1984). These sites were not included in the subsequent gridding procedure.

Finally, we produced gridded data sets of the anomalies of each bioclimatic variable. We used a regular latitude/longitude grid of 2° by 2°. The grid size was selected to be comparable to the grid size typical of the models used in state-of-the-art palaeoclimate simulations. This size of grid also avoids over-smoothing the resulting patterns of regional climate change and limits the number of grid cell values based on a single reconstruction. The grid-cell value of the anomaly was obtained by simple averaging, and that for the uncertainty as a pooled estimate of the standard error obtained in the usual way. The data (Electronic Supplementary Material) were plotted as symbols, where the colour reflects the magnitude of the anomaly and the symbol size reflects the significance of the grid-cell average of that anomaly: a large symbol is plotted when the absolute value of a t-statistic calculated using the anomaly values and the pooled standard error exceeds 2.0.

The product of our synthesis is a set of data sets each consisting of grid-point average anomalies, their uncertainties and associated t-statistics for a particular variable. This product is analogous to a “composite anomaly” map in synoptic climatology, and uses the same general concepts to establish robustness (i.e. invariance with respect to the inclusion of individual cases) and interpretability. At a particular location or grid-point on a composite-anomaly map, if the individual cases that are being composited produce anomalies that are similar in sign and magnitude, then the composite (usually the mean, but other statistics can be used) will take on that sign and magnitude, and change little if a single (or even several) cases are dropped. In contrast, if the individual anomalies differ from case to case in sign or magnitude, then they will cancel each other out in the composite. If the individual cases share the same map pattern (as well as anomaly values at individual grid cells) then that pattern will emerge in the composite-anomaly map, and again, will change little if a single (or several) cases are dropped; if the individual map patterns differ, than no resolvable (or interpretable) pattern emerges. In our analysis, if the individual reconstructions (from different sites and/or studies) are similar within a grid cell, then a significant anomaly will emerge, and if the reconstructions are similar among nearby grid cells, then a coherent map pattern will emerge. Conversely, if the reconstructions differ from one another within and among grid cells, then no discernable large-scale pattern will emerge. Note that there is no requirement that the individual grid-point anomalies that make up broad (regional-to-continental) scale anomaly patterns must individually be significant in order for the pattern as a whole to be interpretable; it would be extremely unlikely that broad-scale patterns would emerge by chance.

3 Results

3.1 Data in climate and vegetation space

The original data set contains 4,806 individual observations for 6 ka, of which 4,706 were used in the final data set after applying age-selection and statistical outlier criteria (Fig. 1; Table 1). For 21 ka, the original data set contains 357 individual observations, of which 349 were used in the final data set (Fig. 1; Tables 1, 2). Gridding results in the creation of climate estimates for 715 cells (out of a possible 3,687 non-ice-covered land cells) at 6 ka and for 153 cells (out of a possible 3,618 non-ice-covered land cells) at 21 ka. The number of individual reconstructions contributing to the estimate for each grid cell is shown in Fig. 2. These numbers reflect the number of sites, samples at each site, and individual studies that have produced reconstructions.

The spatial coverage is not even: there is very good coverage for Europe (north of 35N and west of 45E) (197 out of 354 non-ice-covered land cells at 6 ka) and North America (220 out of 735 non-ice-covered land cells at 6 ka) but coverage is sparse for Eurasia and Africa, and there is information for only one or at best a few cells for other tropical regions. There are no data from extratropical South America or Australia, and only seven points from New Zealand. Despite the uneven geographic coverage, the data set provides a reasonably comprehensive representation of modern climate and vegetation space (from tropical to tundra, and from maritime to continental climates, see Fig. 1, bottom panel). There is also uneven coverage of individual climate variables, with counts ranging from 1,408 for GDD5 to 4,232 for MTWA at 6 ka (yielding 516 and 628 gridded values, respectively) and from 163 for alpha to 333 for MTCO at 21 ka (yielding 88 and 145 gridded values, respectively). Thus, although the data set may fall short of providing a detailed picture of climate for every region of the globe, the data set encompasses a wide range of different climates at present, and so should provide an adequate sampling of the kinds of climate changes that we expect climate models to be able to reproduce.

3.2 Assessment of performance

A basic test of any reconstruction method is to determine how well it predicts observed climates. Here, we compare the reconstructed modern climate variables at individual sites (using the results from all of the different reconstruction methods) with observations derived from the CRU historical climatology data set (CRU CL 2.0, New et al. 2002). The CRU data set provided mean monthly temperature and precipitation data for the period 1961–1990 on a 10′ grid. GDD5 was estimated by direct accounting and alpha was calculated using the approach described by Prentice et al. (1992a, b). The 0 ka reconstructed and CRU values are highly correlated (p ≪ 0.001) for all bioclimatic variables (Fig. 3). For those variables for which reconstructions are available from both North America and Europe/Africa (MTCO, MTWA, MAT and GDD5), the correlation between the 0 ka reconstructed and CRU values do not differ significantly. There are a small number of obvious outliers in the MTCO and MTWA plots; these are mostly from sites at higher elevations or at higher latitudes, where the climate data have been extrapolated from comparatively few climate stations and positional errors are greater.

Fig. 3
figure 3

Correlation between observed and reconstructed modern climate variables. The observed values are derived from the Climate Research Unit (CRU) historical climatology data set (CRU CL 2.0, New et al. 2002). The reconstructed values are derived from the various data sets included in the current synthesis which provided point-based estimates of modern climate. North America includes all sites between 0° and 75°N, and 50°–175° E), Europe and northern Africa includes sites between 0° and 75°N, and 2°E to 60°W, and all the remaining sites are included in the category “other regions”

Full size image

The spatial variability of correlation for the data rich areas of North America and Europe/North Africa can be demonstrated by calculating Pearson’s correlation coefficient between the CRU data (Fig. 4, top panels) and the reconstruction using a moving spatial-window of 20° × 20° (Fig. 4, bottom panels). In North America, MTCO, MAT, MAP and GDD5 show high correlations with little spatial variability. On the other hand, the correlations for MTWA in the west are lower than those in the east, and for alpha are higher on the west, particularly at the mid latitudes (20°–60°N), than those on the east and in the arctic regions. These results indicate that MTWA reconstruction is less reliable in topographically complex regions (such as the American West) and also in regions where summer temperature does not limit vegetation, and that alpha is less reliable in regions where moisture is not a major limiting factor for vegetation. In Europe and North Africa, the correlations are generally high for MTCO, MTWA and MAT with little spatial variability. GDD5 in North Africa shows high correlations in the west and lower correlations to the east.

Fig. 4
figure 4

Maps showing regional climates (upper panels) based on data from the Climate Research Unit (CRU) historical climatology data set (CRU CL 2.0, New et al. 2002) and spatial variations in the strength of the correlations between observed and reconstructed climate variables (lower panels). The spatial variability of correlation is estimated by calculating Pearson’s correlation coefficient between the reconstruction and the CRU data using a moving spatial-window of 20° × 20°. The spatial window shifts 5° at each step in both longitudinal and latitudinal directions. The correlation coefficient is only calculated when the number of available 0 ka reconstruction estimates in each spatial window is equal to or larger than fifteen. The size of the squares are proportional to the correlation values; the center of each square is the location of the center of the corresponding 20° × 20° window

Full size image

Standard goodness-of-fit statistics such as R 2 may overestimate the predictive power of climate reconstructions (especially those made with the modern analogue technique) due to unaccounted-for spatial autocorrelation in the response variables (e.g. Telford and Birks 2005). The extent of this effect in published pollen-based climate reconstructions cannot easily be quantified. However, it should be noted that the spatial autocorrelation of vegetation composition at a regional scale derives almost entirely from its causal relation to climate, provided that attention is confined to variables that influence the growth, establishment and regeneration of plants (Harrison et al. 2009). Spatial pattern in pollen data thus constitutes valuable information for the reconstruction, to be retained rather than rejected (Legendre 1993; Legendre and Legendre 1998). In any case, spatial autocorrelation in pollen data becomes non-significant at length scales of 200–300 km, and is slight at any scale when full taxon lists are used (Sawada et al. 2004). Guiot and de Vernal (2007) showed that the goodness-of-fit (as measured by the R 2 statistic) is an appropriate measure when spatial autocorrelation in the pollen data arises from the underlying climate and not from processes internal to the vegetation system.

3.3 Robustness and uncertainties

As a simple test of the robustness of the final reconstruction, we analysed an early version of the data set to determine whether the choice of sampling points and averaging period could have a substantial impact on the resulting reconstruction maps. We used palaeoclimatic anomalies from all sites in the database having at least three samples within each of the following time windows: 5.5–6.5, 6.5–7.5 and 5.5–7.45 cal kyr BP. Although there are 4,806 6 ka sites in the database, most of these sites only have a single reconstruction for 6 ka. The requirement that a site should have at least three samples within the 6 ka time windows therefore resulted in the use of a more limited number of sites for this analysis. In the first test, reconstructions from the single 6.0 ka sample used for mapping were compared to the average of all samples from that site within the broadest possible temporal envelope defined for ‘6.0 ka’ samples: 5.5–7.45 ka (total 1,097 sites). The second test compared the same 6.0 ka point with the average of all samples between 5.5 and 6.5 ka (1,046 sites). The third test was to compare the differences between average values from 5.5 to 6.5 and 6.5 to 7.5 ka (978 sites).

All three tests indicate a strong linear relationship between each of the sets of values compared (Table 4). Mean values for each of the climatic parameters indicate that sampling choice has a small overall impact on the results (Table 4). Large standard deviations are expected since the dataset includes sites from the tropics to the polar regions and 6 ka anomalies cannot be assumed to be similar in all areas of the globe. Generally the correlation between point and average data was better for temperature parameters than for mean annual precipitation. This is simply a data size issue, since temperature reconstructions exist for almost all sites, while precipitation has been reconstructed for about 15% of sites used in testing. These analyses demonstrate that uncertainties about the exact temporal window-width used in selecting samples have little impact on the final reconstructions.

Table 4 Comparison of means and standard deviations for single point and averaged palaeoclimatic anomalies from the global dataset for 6 ka
Full size table

Figure 5 exploits the overlap of reconstruction data sets to illustrate the robustness of reconstructed spatial patterns with respect to the choice of reconstruction technique. Cheddadi et al. (1997) used a modern analogue approach, while Wu et al. (2007) used model inversion to produce reconstructions of MH climate across Europe. The two data sets share ca. 85% of MH sites for the European window, although Wu et al. (2007) contains many more sites from Russia. The broad scale spatial patterns in GDD5 and alpha obtained in these two studies are surprisingly similar, and the correlations between the two sets of reconstructions at points where they overlap are significant (p < 0.001). Both show warmer than present growing-season conditions over most of Europe giving way to cooler than present conditions around the Mediterranean. Both show wetter than present conditions across Europe with the exception of a region centred on the Baltic Sea. However, there are consistent differences in the magnitude of reconstructed anomalies; for example, Wu et al. (2007) show a generally larger warming, and a greater increase in moisture availability, over northern Europe. There are also differences at smaller spatial scales, such as along the northwestern seaboard (for GDD) and southern Finland (for alpha). This comparison suggests that the robustness of the reconstructions is greater for broader spatial scales than for individual sites, and that analysis and model-data comparison should therefore focus on these scales.

Fig. 5
figure 5

Reconstructed climates (GDD5, alpha) for Europe at 6 ka using different data sets and two different methods: Cheddadi et al. (1997) used a constrained modern analogue approach, while Wu et al. (2007) use inverse modelling. Large symbols are used to indicate grid points with significant anomalies (i.e. those that exceed twice the standard error of the reconstructions) while small symbols are used to indicate anomalies that are not significant by this measure

Full size image

3.4 Climate patterns at 6 ka

There are coherent large-scale patterns in the gridded reconstructions for 6 and 21 ka (Figs. 6, 7). The robustness of these patterns can be assessed in two ways: from the statistical significance of the individual grid cell values (shown by the size of the circle) and by the coherency of the sign and magnitude of the anomalies (shown by circle colour). Thus, features of these patterns which are non-significant statistically (because the change is small or there is no change, or because there are only a limited number of observations) may still be physically plausible and should be considered as robust features of the reconstructed changes in regional climate for the purposes of model evaluation.

Fig. 6
figure 6

Reconstructed a GDD5, b MTWA, c MTCO, d MAT, e MAP, and f alpha at 6 ka. Large symbols are used to indicate grid points with significant anomalies (i.e. those that exceed twice the pooled standard error of the reconstructions) while small symbols are used to indicate anomalies that are not significant by this measure. Note that the anomaly patterns are spatially coherent and include patches of both significant and not-significant values

Full size image
Fig. 7
figure 7

Reconstructed a GDD5, b MTWA, c MTCO, d MAT, e MAP, and f alpha at 21 ka. Large symbols are used to indicate grid points with significant anomalies (i.e. those that exceed twice the pooled standard error of the reconstructions) while small symbols are used to indicate anomalies that are not significant by this measure. Note that the anomaly patterns are spatially coherent and include patches of both significant and not-significant values

Full size image

The 6 ka GDD5 reconstructions (Fig. 6a) show large and spatially coherent differences from present. Northern Europe is characterised by an increase in growing season warmth compared to present, with anomalies in GDD5 of >500 kdays. In contrast, growing season warmth was reduced by >1,000 K day compared to present over most of southern Europe. A tri-partite latitudinal pattern is seen in Eurasia, with moderately increased growing season warmth compared to present along the Arctic coast, a reduction in GDD5 in the zone between ca. 45°–65°N, and increased growing season warmth in the continental interior. The northern coast of North America is also characterised by positive anomalies in GDD5, as is the Pacific Northwest away from the coast and most of the northeastern part of the continent. The American southwest and the southern part of the continent is characterised by a reduction in GDD5 compared to present, as is the coastal strip of the Pacific Northwest and the interior of northwestern Canada. Equatorial Africa and some sites in the western Sahara are characterised by reduced GDD5, whereas the surrounding zone including the eastern Sahara, eastern Africa and southern Africa was characterised by >1,000 K day increase in GDD5 compared to present.

The accumulated temperature sum during the growing season, as expressed by GDD5, is partially determined by changes in the length of the growing season and partially by changes in summer temperatures. It is therefore unsurprising that the 6 ka MTWA reconstructions (Fig. 6b) show similar, though somewhat less coherent, patterns of change relative to present to those evident in the GDD5 reconstructions. The similarity between the GDD and MTWA is clearly evident in Europe, Eurasia and Africa, and reflects the spatial coherence of the anomalies of both variables at a regional scale. The reconstructions of MTWA over much of North America have large uncertainties and this reduces the apparent coherence with the GDD5 reconstructions over this continent. Nevertheless, the contrast between warmer-than-present conditions in the northeast and colder-than-present conditions to the south and west is apparent in both data sets.

The reconstructed changes in 6 ka MTCO (Fig. 6c) are in general smaller and less significant than the changes in GDD5. Nevertheless, the spatial patterns are regionally coherent. The maritime region of western Europe is characterised by cooler conditions than present, by ca. 1°C in the north and up to 4°C in the south. The Mediterranean region as a whole is also significantly cooler (>2°C) than present. Central Europe and western Eurasia, however, are characterised by warmer conditions than today, with temperatures increased by as much as 2–4°C at some sites. The reconstructed MTCO changes in eastern Siberia are not significant but show a slight cooling compared to present. Reconstructed temperature changes in North America show a dipole structure, with positive MTCO anomalies in the northeastern part of the continent and negative MTCO anomalies in western and southern regions. In contrast to this regional pattern, however, some sites around Hudson Bay show a reduction in MTCO compared to present by ca. 1°C. Equatorial Africa is characterised by reduced MTCO, whereas the surrounding zone including the eastern Sahara, eastern Africa and southern Africa was characterised by an increase in MTCO.

In many regions of the northern extratropics, the differences between the MTWA and MTCO reconstructions show the expected impact of precession-driven insolation changes (i.e. increased seasonality). However, in some regions the reconstructions show no apparent change in seasonality or even changes in the opposite direction from expected (e.g. around Hudson Bay in eastern North America, the peri-Baltic region of Europe, and central Eurasia). As will be discussed later, these are predictable, second-order responses to atmospheric circulation changes.

The change in growing season conditions, as shown by GDD5, is the strongest influence on reconstructed changes in MAT (Fig. 6d). The largest and most significant changes in MAT occur where the changes in summer and winter temperatures are in the same direction: thus, the significant reduction in MAT in southern Europe compared to present reflects cooler conditions in both seasons, while the significant increase in MAT in central Europe reflects warmer conditions in both seasons. In the maritime regions of northeastern Europe and northern Scandinavia, the reconstructed change in MAT is smaller, because the seasonal changes in these regions are of opposite sign. A similar situation is apparent in North America, with large changes in MAT in regions where both seasons are characterised by warming or cooling compared to today, and smaller (or non-significant changes) in regions where the temperature changes result in increased seasonality. The changes in MAT across Africa reflect the changes in GDD5, with cooling (>2°C) compared to present in the western Sahara and the equatorial zone and warming elsewhere. Cooler conditions compared to present are shown by the MAT reconstructions over most of China; the exceptions mostly occur in the extreme north and northwest. In the absence of MTCO or MTWA reconstructions, it is unclear whether this is due to changes in both seasons but the coincidence of cooler conditions with increased precipitation (see below) suggests that the change in MAT reflects both orbitally induced winter cooling and summer cooling consequent on increased cloudiness in regions affected by the enhanced Asian monsoon during the MH.

The reconstructions show increased precipitation compared to present across Europe and Eurasia (MAP: Fig. 6e); the changes are small (<100 mm) and largely non-significant across most of Europe, but spatially very coherent. There are significant increases of between 100 and 500 mm compared than today in the Mediterranean region. Significantly wetter conditions (ca. 100 mm) than today are also registered across Eurasia. Increased precipitation is also shown across much of western and southern North America. The increase in the Great Basin is not significant although it is spatially coherent. Significant increases in precipitation compared to present, generally of >200 mm and in some localities of >500 mm, are registered across northern Africa. The equatorial zone is characterised by reduced precipitation compared to present, as is the northern part of the East African Rift Valley and Madagascar. Most of southern Africa is characterised by increased precipitation compared to present, except for sites near the south coast. An increase in precipitation compared to present is shown across most of China.

These reconstructed changes in precipitation are closely related to reconstructed changes in plant-available moisture (alpha: Fig. 6f). However, alpha is also affected by changes in growing-season temperature and evapotranspiration. Thus, the reconstructed patterns of change in alpha are in general less coherent and less significant than those shown by precipitation. Nevertheless, in regions characterised by very large changes in precipitation, there are congruent, large and significant changes in alpha. The Mediterranean region, for example, is characterised by a significant increase in precipitation compared to present and this translates into an increase in alpha of >0.2. Similarly, eastern North America is characterised by a decrease in precipitation compared to present, which translates into a decrease in alpha of ca. 0.5. The interplay between precipitation and temperature-controlled changes in evaporation, however, can result in significant changes in alpha. In southwestern North America, for example, the reconstructed increases in precipitation compared to today are small and non-significant. However, as a result of cooler temperature throughout the year, this negligible change in precipitation results in a significant increase in alpha.

One indication of the reliability and robustness of these climate reconstructions is the fact that the changes in individual variables are mutually consistent with one another, spatially coherent and statistically significant. Furthermore, the reconstructed changes are physically plausible and can be explained in terms of current understanding of the hierarchy of controls on regional climates (see e.g. Bartlein et al. 1998). The ultimate cause of changes in regional climate during the mid-Holocene is differences in the seasonal and latitudinal distribution of incoming solar radiation (insolation). Broadly speaking, as a result of changes in precession, insolation was higher during summer and lower during winter in the northern hemisphere. This resulted in increased summer temperatures in the mid- to high-latitudes of the northern hemisphere (see e.g. Braconnot et al. 2007a, b), a major feature of the reconstructions for northern Europe, and the Arctic coasts of both northern Eurasia and North America (Fig. 6a, b). Heating over the northern continents increased land–sea temperature contrast and resulted in the enhancement of the northern hemisphere monsoons (see e.g. Kutzbach et al. 2001: Braconnot et al. 2007a, b), again a major feature of the reconstructions from northern Africa and Asia, and apparent to a lesser extent in the increased precipitation found in southern and southwestern North America (Fig. 6e). Monsoon expansion in Africa results in a latitudinal shift in the rainbelts, and thus the reconstruction of reduced precipitation in the equatorial zone is consistent with orbitally induced changes in monsoon circulation. In North America, upward motion associated with the expanded monsoon is compensated by subsidence in the regions surrounding the monsoon core: thus the reconstruction of decreased precipitation (Fig. 6e) and moisture availability (Fig. 6f) in the Pacific Northwest and in the mid-continent of North America is consistent with orbitally induced circulation changes (Harrison et al. 2003).

Although the insolation-induced expansion of the northern hemisphere monsoons is apparent in reconstructed precipitation (Fig. 6e) and alpha (Fig. 6f), changes in summer temperature (measured either by GDD5 or MTWA) in the monsoon regions are less coherent. Increased summer temperature is registered at some sites in northern Africa but is not a noticeable feature of the American southwest. Furthermore, we have interpreted the reconstruction of cooler temperature year-round in China as indicating that the winter cooling associated with the stronger winter monsoon was not noticeably offset by a marked summer warming. The absence of a strong warming signal in the monsoon regions might seem surprising. However, although enhanced land-sea temperature contrast is a necessary prerequisite for monsoon enhancement (see e.g. Kutzbach et al. 2001), maintenance of monsoon precipitation is strongly dependent on latent heat flux and the energetics of the system, and their influence on atmospheric dynamics (Ramage 1971; Trenberth et al. 2000; Webster et al. 2002; Anderson et al. 2009; Taylor et al. 2010). Modern observations show that monsoon rains are accompanied by a strong regional temperature depression, and a similar phenomenon is presumably the cause of the weak and spatially less coherent structure of reconstructed summer temperature changes during the mid-Holocene.

One direct consequence of northern-hemisphere insolation changes during the mid-Holocene is increased seasonality in temperatures, but this is not a prominent feature of the reconstructions. Although the combination of colder winters and warmer summers is found in e.g. northwestern Europe, most regions show either an increase in temperature in both seasons (e.g. central Eurasia, northeastern North America) or a decrease in temperature in both seasons (e.g. southern Europe, southern and southwestern North America). In these regions, the direct effects of insolation changes are overwhelmed by changes in atmospheric circulation regimes. Thus, cooler temperatures year-round found in southern and southwestern North America reflect the impact of increased onshore flow of moist air and increased cloudiness associated with the strengthened North American monsoon (Harrison et al. 2003). Similarly, the increased winter temperatures characteristic of central Eurasia are a result of increased warm-air advection from the North Atlantic (Harrison et al. 1992; Bonfils et al. 2004).

At a more local scale, changes in land-sea geography also have an impact on regional climates that can be seen in the reconstructions. Both Hudson Bay and the Baltic Sea were more extensive during the mid-Holocene than they are today. In both regions, this resulted in local and year-round cooling (Fig. 6b, c) and, at least in the case of the peri-Baltic regions, it is likely that the enhancement of anticyclonic conditions over the region led to increased aridity (Fig. 6f).

Insolation changes in the northern and southern hemisphere between the mid-Holocene and today are not perfectly symmetrical. Although insolation is less in southern-hemisphere summer and more in southern-hemisphere winter, the largest increase in insolation occurs during the spring months (September through November). This asymmetry is reflected in the climate reconstructions for southern Africa, the only area of the southern hemisphere extratropics for which we have quantitative reconstructions. The increase in insolation during winter is responsible for the increase in MTCO relative to present (Fig. 6d). In general, the onset of the growing season is earlier when winter temperatures are higher, and this effect combined with a significant increase in insolation during the spring results in an increase in both GDD5 (Fig. 6a) and MTWA (Fig. 6b). Increased temperatures year-round offset, to some extent, the reconstructed increase in precipitation (Fig. 6e) such that the area which shows increased alpha compared to today (Fig. 6f) is less extensive than the area with increased precipitation. The reconstructed changes in southern Africa regional climates are, to first order, consistent with orbital forcing.

3.5 Climate patterns at 21 ka

The basic picture emerging from the 21 ka reconstructions is of colder conditions year-round across all regions (Fig. 7a–d). The change in MAT was largest (>8°C, Fig. 7d) close to the northern-hemisphere ice sheets and comparatively muted (<3°C, Fig. 7d) in the tropics. In general, the changes in winter temperatures as represented by MTCO (Fig. 7c) are larger than the changes in summer temperatures as represented by either GDD5 (Fig. 7a) or MTWA (Fig. 7b).

There are elements of the reconstructions that are not congruent with the general picture of year-round cooling: the reconstruction of MTWA warmer than today in Africa (Fig. 7b), and the reconstruction of warmer conditions than today as shown by one or more of the temperature-related variables at a few sites in the northern hemisphere extra-tropics, most noticeably in Alaska.

The MTWA reconstructions for Africa show temperatures between 1 and 4°C warmer than today (Fig. 7b, see also Wu et al. 2007). This is inconsistent with the reconstruction of reduced GDD5 (Fig. 7a) and MAT (Fig. 7d), unless the climate was characterised by a shift towards a short extremely hot summer in both hemispheres. The insolation patterns at 21 ka are inconsistent with such a change in seasonality in both hemispheres. The increased warm-season temperatures are unlikely to be a response to changes in advection, given the reconstruction of reduced precipitation compared to today across most of the continent (Fig. 7e). They are similarly unlikely to be a response to changes in the local partitioning of latent and sensible heat, due to reduction in cloudiness, because of their widespread expression. No physically plausible mechanism has been advanced to explain the warming.

The reconstruction of conditions similar or slightly warmer than today in Alaska during the LGM is registered chiefly in MTCO (Fig. 7c) and MAT (Fig. 7d) although single points show warmer values for other seasonal temperature variables. In contrast to the situation in Africa, there is a plausible mechanism for why this region might not experience the major cooling shown over the rest of the continent (Bartlein et al. 1998). The Laurentide ice sheet was sufficiently large to cause a major reorganisation in the atmospheric circulation pattern. This reorganisation could have resulted in more southerly onshore flow into Alaska, which would have produced advective warming of the region year-round.

The reconstructions (Fig. 7e) show a significant decrease in precipitation compared to present across Europe (>200 mm), Eurasia (>100 mm), the Pacific Northwest (>500 mm), eastern North America (>500 mm) and Africa (>200 mm) at 21 ka. Reduced vigour of the hydrological cycle and less precipitation are consistent with the colder glacial conditions at 21 ka. Despite the generally colder conditions, the reduction in precipitation was translated into a reduction in plant-available moisture (Fig. 7f). The only regions that appear to have registered an increase in both precipitation and plant-available moisture compared to today are the Great basin of North America and the southernmost part of southern Africa. Shifts in atmospheric circulation, specifically the displacement of the track of the westerlies by the Laurentide ice sheet, have been invoked to explain increased precipitation in the Great Basin during the LGM (Benson and Thompson 1987; Hostetler and Benson 1990). Changes in the southern westerlies have been invoked to explain increased precipitation during the LGM in southern Africa (e.g. Stuut and Lamyb 2004; Justino et al. 2005).

Low CO2 levels during the glacial affect vegetation by reducing the growth of C3 plants, and allowing C4 species to compete more successfully (Harrison and Prentice 2003: Prentice and Harrison 2009). This can lead to a shift towards more open vegetation, independent of climate change (Jolly and Haxeltine 1997; Cowling and Sykes 1999) and could result in reconstructions that are drier and/or colder than the true glacial climate (Prentice and Harrison 2009). The impact of CO2 has been taken into account in the inverse-modelling reconstructions for Europe, Eurasia and Africa (Wu et al. 2007) and China (Guiot et al. 2008). The impact of CO2 could mean that the reconstructions for North America, which are primarily based on analogue or response surface approaches, overestimate the degree of drying and cooling. Nevertheless, the magnitude of changes e.g. in eastern North America and Europe are similar. This suggests that the impact of CO2 does not dominate the reconstructed signals in North America, a conclusion consistent with recent palaeovegetation reconstructions for eastern North America based on pollen data and BIOME4 simulations (Gonzales et al. 2008; Williams et al. 2008).

4 Discussion and future perspectives

The data presented here are a resource for model-data comparison and are available at http://pmip3.lsce.ipsl.fr/. Comparisons should focus on the large-scale spatial patterns that are significant, robust and climatologically interpretable. Given that there is some redundancy among the reconstructed climate variables, we recommend focusing comparisons on GDD5, MTCO and MAP. GDD5 and MTCO represent independent controls on the distribution and abundance of plant taxa (Harrison et al. 2009). MTWA and MAT provide little additional information beyond what is already contained in GDD and MTCO. MAP reflects a third independent control (plant-available moisture) less directly, but we show here that it can be estimated with greater accuracy (Fig. 3) and that it yields more strongly coherent patterns than the alpha index (Figs. 6, 7), which might be preferred on theoretical grounds.

We have of necessity adopted rather broad time windows to represent the MH (potentially 5.5–7.45 ka) and LGM (potentially 20–22.9 ka) in order to include information from as many studies as possible. Some of the breadth is due to the inclusion of older data sets in which the time window had been defined around uncalibrated 14C ages, as well as data sets in which the time window had been defined around calendar ages. It would be desirable for future data compilations to be based on time windows defined around calendar ages only. It would also be desirable for a narrower window to be adopted for the LGM, to avoid potential inclusion of samples from associated Heinrich events (HS2 at 26.5–24.3 ka and HS1 at15.6–18 ka, Sanchez-Goñi and Harrison, in press).

A possible source of ambiguity in climate reconstruction is the choice of modern climate data for calibration. This has been recognized as an issue for the reconstruction of sea surface temperatures (Schäfer-Neth and Paul 2003) and it is equally an issue for climates over land. There is no strong basis for recommending one global climatology over others, but many analyses have used the CRU climatological data set for 1961–1990 as we have here. There is considerable merit in adopting a standard, even if it is imperfect. The CRU data set has the merits of applying to a consistent observation period and being widely used.

We have demonstrated that reconstructions made using different methods yield similar results for large-scale spatial patterns (Figs. 6, 7), even if the reconstructions for individual sites can differ. We suggest that the way forward is to use inverse modelling in preference to analogue or response-surface techniques, because inverse modelling (a) provides sound logic to cope with the no-analogue problem, (b) makes the wrong-analogue problem explicit, (c) allows the consequences of changing atmospheric CO2 concentration to be taken into account (Harrison and Prentice 2003; Prentice and Harrison 2009) and (d) allows the consequences of changing seasonality of climate forcing or climate responses to be taken into account (Guiot et al. 2000, 2009). However there are many detailed methodological choices to be made in this approach, such as the assignment of taxa to PFTs, choice of model, and the choice of model variable to compare with PFT abundance and more work is required before this technique will entirely supersede other methods.

An important limitation of the data presented here is the absence of quantitative reconstructions from South America and the southern-hemisphere extratropics generally (except for southern Africa). The primary pollen data, on which reconstructions are based, is also generally sparse in the tropics. To some extent, these limitations reflect the history of pollen analysis and the compilation of pollen data sets. However, the BIOME 6000 project (Prentice et al. 2000) has stimulated the compilation of new regional pollen data sets, e.g. from Australasia and southeast Asia (Pickett et al. 2004), South America (Marchant et al. 2009) and the Indian subcontinent (Sutra et al., in prep.). It is now possible and timely to exploit these compilations for the reconstruction of regional climate changes. There are additional ongoing efforts to expand pollen data coverage in some regions where existing data are sparse, such as central Asia (e.g. Cordova et al. 2007).

Although the greatest volume of quantitative climate reconstructions exists for the two time periods analysed here, some of the underlying data sets include reconstructions at more closely spaced time intervals; synthesis of these for, say, 1,000-year intervals would be worthwhile for comparison with a new generation of transient climate model simulations, as proposed in PMIP3. Further analyses of the data presented here are also possible, for example examination of elevational gradients in climate anomalies which can give further insight into underlying climate processes (cf. Farrera et al. 1999; Kageyama et al. 2005).

We have focused on quantitative reconstructions based on pollen and/or plant macrofossil data, but site-based climate reconstructions have been made using chironomid or coleopteran assemblages (e.g. Atkinson et al. 1987; Heiri et al. 2003), and using geochemical techniques (e.g. Stute et al. 1992; Jones et al. 2004). Synthesis of such data could provide a way of checking the reconstructions presented here, but more importantly of expanding the regional coverage of quantitative reconstructions.

The changes in regional climates inferred from pollen and/or plant macrofossil data could also be compared with records that provide more qualitative indications of past climate changes, for example lake-level records (see e.g. Street-Perrott and Harrison 1984; Harrison and Digerfeldt 1993; Harrison et al. 1996). However, some caution needs to be exercised here, to ensure that the sensors are responding to the same aspect of seasonal climate. Vegetation, for example, responds to plant-available moisture whereas lake-levels respond to runoff (i.e. precipitation in excess of that taken up by plants) (Cheddadi et al. 1997). In cases where there is a large enough change in precipitation, and/or plant growth is limited by seasonal temperatures, both plant-available moisture and runoff may increase and hence both sensors will appear to register wetter conditions (see e.g. Thompson et al. 1993; Yu and Harrison 1996; Jolly et al. 1998; Shuman et al. 2002). However, in some situations, lake levels may be high at the same time as the vegetation is registering drier conditions. The classic example is the co-existence of high lake levels and steppic vegetation in southern Europe during the Last Glacial Maximum (see Prentice et al. 1992a, b), where the vegetation is responding to a reduction of precipitation during the growing season while the lakes reflect increased runoff during the winter.

Finally, recent advances in data assimilation open up the possibility that palaeodata could be used to constrain the parameters of climate models and thereby improve reliability and realism, and assign quantitative uncertainties, in future climate prediction. Systematic approaches to data assimilation have been developed for present day climate and Earth system model components with reduced spatio-temporal resolution, simplified process representations, or reduced sets of parameters (e.g. Marotzke et al. 1999; Rayner et al. 2005). A system to apply these approaches to climate modelling has been demonstrated (Kaminski et al. 2007), and a further natural development of this system would allow fully coupled, dynamic climate system models to be constrained with palaeoclimate reconstructions using this variational approach. A rigorous method to combine models and palaeodata should help to provide a better understanding of the interactions between different components of the climate system, and establish a tighter link between palaeoclimate observation and Earth system modeling.