The effect of calibration data set selection on quantitative palaeoclimatic reconstructions

Introduction

Quantitative reconstructions based on fossil biological proxies are one of the major sources of long-term palaeoclimatic information. Such reconstructions are particularly important for validating climate-model results (e.g. Braconnot et al., 2012). Assessing the robustness and reliability of these reconstructions is thus essential (Juggins, 2013; Telford and Birks, 2011). In such reconstructions, palaeoclimatic estimates are derived from the observed variability in fossil assemblages using numerical calibration methods (Birks et al., 2010; Brewer et al., 2008). Typically, a modern calibration data set (CD) is used to describe the variability of the studied biological group in modern climate space. In the case of Quaternary microfossils (e.g. pollen, chironomids), CDs generally consist of modern assemblages analysed from surface-sediment samples from lakes, bogs or soils, with modern climate data associated with each sample. These modern assemblages and climate values are used to model climate–taxon relationships in regression and calibration methods such as multivariate transfer functions (TFs; Birks et al., 2010) or in the modern-analogue technique (MAT; Birks et al., 2010; Jackson and Williams, 2004; Overpeck et al., 1985; Simpson, 2012).

Highly diverse CDs are used in palaeoclimatic reconstructions differing, for example, in their geographic location, spatial extent, type and number of modern samples and level of taxonomic detail. The CD samples, commonly numbering 30–200 and up to thousands in some pollen data sets, are collected along major climatic gradients, in a geographical region containing the fossil sites. However, CD sampling occurs under many constraints. The most fundamental constraint is set by the distribution of land and sea in modern climate space. Only a subset of all possible climates is thus available for sampling in any given region (Jackson and Overpeck, 2000). Also, in some areas, sampling sites may not be available due to the lack of suitable lakes or bogs. Thus, the selection of CDs is often constrained, and samples are selected from where they can be readily obtained. As preparation of CDs is time-consuming and expensive, extensive use is often made of available modern data and CDs are used with fossil sequences from sites in different parts of the modern calibration region. In some cases, due to lack of calibration data, reconstructions are made with ‘extra-regional’ CDs, located hundreds or thousands of kilometres from the fossil site (cf. Brooks and Birks, 2000; Engels et al., 2008; Self et al., 2011). There is thus great variability in the geographical and climatic extents of the CDs used, and the spatial displacements between the CDs and fossil sites. This variability is, to a major extent, dictated by the availability of data, and can thus be due to chance rather than to design.

It is possible that the geographical region and climate space represented in the CD may have a major effect on climate reconstructions, and CD selection is thus a potential cause of between-reconstruction variability. For example, expanding a pollen–climate CD into the arctic tundra may significantly lower the estimates of climatic optima of many non-arboreal pollen types, subsequently lowering palaeotemperatures reconstructed from fossil assemblages containing abundant non-arboreal pollen. Also, changing the position of a CD along a secondary gradient (e.g. moving a summer temperature CD to a different seasonality, winter temperature or precipitation regime) may affect the climatic values indicated by particular taxa, as the variability of the proxy can also be affected by secondary climatic determinants (Jackson et al., 2009; Juggins, 2013; Salonen et al., 2012b; Self et al., 2011). It is thus probable that, to some degree, taxa may indicate different climates as the spatial position and extent of the CD vary. Consequentially, the choice of CD may be a major source of the variability seen between different reconstructions from fossil biological proxies. Despite this, the effect of CD selection on reconstructions has received relatively little attention (but see Bjune et al., 2010; Salonen et al., 2012a; Velle et al., 2011; Williams and Shuman, 2008).

Here, we explore the role of CD selection as a source of bias, by deriving July mean temperature (T_jul) reconstructions from fossil pollen sequences, using four randomly selected CDs from different parts of northern Europe. We thus use these randomly selected, spatially shifted CDs to simulate the effect of spatially constrained CD selection on climate reconstructions. We find significant variation in the resulting reconstructions and suggest new approaches in the construction and use of CDs and in the interpretation of reconstructions.

Materials and methods

We select four random, 60-sample CDs from a 526-sample northern European CD (Figure 1a). The random selection is done with climatic constraints to shift the CD locations along two climatic gradients: the reconstructed variable T_jul, and the Gorczynski (1920, 1922) continentality index (K_G), a major secondary climatic gradient in the region (Figure 1b). The base set 60S (Figure 1c) is located in a low-continentality regime (K_G < 25) with T_jul of 11–17°C. 60S_high-KG and 60S_wide-KG have the same T_jul constraints as 60S, but 60S_high-KG samples a higher-continentality regime (K_G > 25), while 60S_wide-KG samples the entire K_G gradient. 60S_low-Tjul has the same K_G constraint as 60S but is shifted into a cooler T_jul regime (9–15°C). We use a stratified randomization to produce CDs which are highly comparable apart from their positions along climatic gradients. All CDs have identical sample counts (n = 60) and near-identical T_jul gradient lengths and distributions of samples along the T_jul gradient (Figure 1d). For each CD, five samples meeting the K_G criterion were picked from each 0.5°C T_jul segment (in two cases, one sample was picked from an adjacent segment due to lack of samples), in order to establish an even sampling along the climatic gradient. 60S_wide-KG was additionally stratified by K_G to sample the K_G gradient evenly, with one sample selected from each of the five equal-length K_G segments, within each T_jul stratum.

OPEN IN VIEWER

Figure 1.

(a) Locations of all modern pollen samples (Salonen et al., 2012b). The locations of the fossil data sites used are also indicated. (b) Climate maps for modern mean July temperature (T_jul) and Gorczynski (1920, 1922) continentality index (K_G), calculated based on the WorldClim (Hijmans et al., 2005) modern climate grids. (c) Four 60-sample subsets are selected from all 526 samples, from different locations along the T_jul and K_G gradients. Subsets are selected randomly from all samples meeting the T_jul and K_G selection criteria for each set. All sets are stratified by T_jul, with five samples selected from each 0.5°C T_jul interval. Subset 60S_wide-KG is also stratified by K_G to sample evenly the wider K_G gradient. (d) Distribution of calibration data set samples along modern T_jul and K_G gradients.

T_jul reconstructions are presented based on fossil pollen sequences from four lakes (see Figure 1a for locations). Two of the fossil data sets have been published earlier, Laihalampi in Heikkilä and Seppä (2003) and Arapisto in Sarmaja-Korjonen and Seppä (2007). The Svartvatnet and Reiarsdalvatnet fossil data sets are previously unpublished. Svartvatnet and Reiarsdalvatnet pollen diagrams are shown in Supplementary Figures 1 and 2 and radiocarbon ages in Supplementary Tables 1 and 2.

Reconstructions are prepared from each fossil site with all four 60-sample CDs, using weighted-averaging (WA) regression and calibration with inverse deshrinking (Birks et al., 1990) (see Supplementary Figure 3 for cross-validation results). We use WA because it is a commonly used and robust reconstruction method and also allows the examination of the taxon optima and tolerances which underpin the reconstructions (Birks et al., 2010). We calculate the WA taxon optima and tolerances in each CD for all taxa which occur in at least 12 samples in all CDs. Our stratified and evenly sampled CD selection procedure permits the robust comparison of WA taxon responses. The unimodal taxon-response models of WA are sensitive to sample placement along the gradient. However, with these evenly sampled CDs, the WA taxon-response models approximate those calculated with Gaussian logit regression, a generally more robust method of estimating taxon optima (ter Braak and Looman, 1986). Although we here use pollen data and WA, the issues and results discussed here are relevant to reconstructions based on other fossil proxies and other reconstruction methods that involve modelling climate–taxon relationships with modern CDs.

To validate the different WA-based reconstructions, we compare them with MAT reconstructions based on the weighted mean of the 10 closest analogues from all 526 samples, with squared-chord distance (Overpeck et al., 1985) as the dissimilarity measure (see Supplementary Figure 4 for cross-validation results). All reconstructions are also compared with modern T_jul values for the fossil sites, extracted from the WorldClim (Hijmans et al., 2005) mean T_jul grid (lapse-rate corrected as in Salonen et al. (2012b)).

All T_jul reconstructions were calculated in C2 software (Juggins, 2007), with pollen data square-root transformed to reduce noise (Prentice, 1980). Most reconstructions are statistically significant (p < 0.05) according to the Telford and Birks (2011) randomization test (see Supplementary Table 3 for details). All terrestrial pollen and spore taxa were used.

Results and discussion

All CDs and MAT produce palaeotemperature reconstructions with similar shapes for T_jul (Figure 2). However, the absolute temperatures are markedly affected by the location of the CD. For example, although 60S and 60S_high-KG sample identical T_jul gradients, the 60S_high-KG-based reconstructions have higher values throughout, with a difference of c. 0.5–2.5°C compared with the 60S-based reconstructions. Comparing the reconstructed T_jul from the uppermost samples with modern observed T_jul, the match is closest with local CDs (60S for Svartvatnet and Reiarsdalvatnet and 60S_high-KG for Arapisto and Laihalampi). Reconstructions with locally representative CDs also tend to show closer agreement with MAT-based reconstructions. This suggests that WA reconstructions based on CDs from the modern continentality regime of the fossil site are more robust. When comparing WA taxon optima in 60S and 60S_high-KG (Figure 3), most taxa indicate higher temperatures in 60S_high-KG. These higher-taxon optima contribute to the higher reconstructed temperatures with 60S_high-KG. Reconstructions with the full-continentality 60S_wide-KG are similar to 60S-based reconstructions (Figure 2), and 60S_wide-KG taxon optima are also closer to 60S than with 60S_high-KG (Figure 3). Comparing 60S_low-Tjul with 60S, the 60S_low-Tjul produces lower palaeotemperatures (Figure 2) due to the lower WA taxon optima (Figure 3).

OPEN IN VIEWER

Figure 2.

Reconstructions from the four fossil sites using all four calibration data sets. The present-day continentality index (K_G) value at each fossil site is indicated. Error bars show sample-specific standard errors estimated by a 1000-cycle bootstrapping procedure (Birks et al., 1990). A LOESS smoother (span 0.15, one robustness iteration) is fitted through each reconstruction. The modern observed T_jul values (dashed lines) and reconstructions based on the MAT using all 526 calibration samples are shown for comparison.

OPEN IN VIEWER

Figure 3.

Comparison of WA taxon optima for T_jul in the calibration data sets. A total of 28 taxa, which occur in at least 12 samples in all four calibration sets, are considered. Box-plots show the WA optima of the taxa based on calibration sets 60S_high-KG, 60S_wide-KG and 60S_low-Tjul, expressed as differences to the optima of the same taxa based on 60S. Also shown are Spearman rank correlation coefficients (ρ) between the taxon optima based on 60S and the optima based on 60S_high-KG, 60S_wide-KG and 60S_low-TjulThe WA optimum and tolerance estimates of the individual taxa are listed in Supplementary Table 4.

The differences between 60S_low-Tjul with 60S are not surprising considering the general sensitivity of WA to the placement of CD samples along the reconstructed gradient (ter Braak and Looman, 1986). However, the effect of CD location along the secondary K_G gradient (60S vs 60S_high-KG) is remarkably large. The higher WA optima and higher reconstructed T_jul in the more continental east are probably due, in part, to the effect of other seasonal temperatures or composite climatic factors apart from summer temperatures. This is ecologically credible, as it is well established that in northern Europe, plants may require higher summer temperatures in the more continental east to compensate for the shorter growing season there. Conversely, plants may tolerate lower summer temperatures in the west where winters are mild (Giesecke et al., 2008). A similar pattern has been discovered in chironomid–T_jul CDs from Norway and European Russia, with taxon optima systematically higher in Russia compared with the more oceanic Norway (Self et al., 2011).

The same plant and insect taxa can, in effect, indicate different values for the reconstructed climatic variable given different distributions along the secondary gradients. This has some major implications for reconstructions. First, when using reconstruction methods based on fitting taxon-response models to modern data, it is vital that the CD represents a climatic regime similar to that of the palaeoclimatic period being studied. Second, reconstructions using extra-regional CDs may produce significantly biased absolute palaeoclimatic values if the CD represents a different climatic regime (e.g. different continentality) compared with the fossil site. Third, reconstructions over long timescales and covering climate regimes markedly different from those of today are likely to present considerable challenges, as the joint distribution of primary and secondary variables increasingly deviates from the present-day pattern. For example, the decoupled changes in summer and winter insolation are likely to have caused significant temperature seasonality changes during the last Milankovitch precessional cycle (Jackson and Overpeck, 2000). In reconstructions reaching into early Holocene and glacial climates, these independent seasonal shifts present significant challenges to finding representative calibration data for different time windows (Birks et al., 2010; Jackson et al., 2009; Salonen et al., 2013).

Thus, the methodological challenge is that the estimated climatic optima of taxa in relation to different climatic parameters – the cornerstone of many reconstruction methods – are not constant but instead can vary in both space and time (Veloz et al., 2012). The spatial variability of climatic optima causes reconstruction methods based on taxon-response models to be sensitive to CD spatial patterns (Juggins, 2013). Temporal variability causes taxa in the same region to peak at different climate parameter values in different time periods, making it difficult to rely on any single CD over long timescales (Jackson et al., 2009; Salonen et al., 2013).

One common reconstruction method – MAT – does not fit taxon-response models, but rather, the palaeoclimatic estimates are based on modern samples with assemblages most similar to the fossil sample (Birks et al., 2010; Overpeck et al., 1985; Simpson, 2012). MAT reconstructions typically use large, continental-scale CDs as a pool of possible modern analogues for different palaeoclimates. MAT, combined with large CDs, thus circumvents the problem of finding applicable CDs for different palaeoclimatic regimes, which is an advantage of MAT, although the method also has its own weaknesses (Birks et al., 2010; Simpson, 2012) and challenges in CD selection (Williams and Shuman, 2008), most particularly the frequent absence of modern analogues for particular fossil assemblages (Jackson and Williams, 2004; Simpson, 2012).

We find reconstructed temperatures to be affected significantly by CD location and the climatic setting represented in the CD used. This needs to be considered when interpreting reconstructions and in climate-model validation with palaeodata. By comparison, while the absolute values are highly sensitive to CD location, the shapes of the relative palaeotemperature curves seem comparatively robust, as the curve shapes mostly remain similar as the CD is spatially shifted. However, significant variation remains in the relative temperature curves. For example, the magnitude of late Holocene cooling since 6000 cal. yr BP varies considerably, by up to 1–1.5°C, depending on the CD used (Figure 2). This significant between-reconstruction variability stresses the need for ensemble reconstructions combining different proxies, reconstruction methods and fossil sites (e.g. Brewer et al., 2008; Kaufman et al., 2009; Mann et al., 2009), to control for reconstruction-specific biases.

The relative robustness of the curve shape is likely due to the taxa occupying similar relative positions along the sampled T_jul gradient, regardless of CD location, as shown by the high Spearman rank correlation coefficients between taxon optima of the CDs (Figure 3). Thus, cold indicators remain cold indicators and warm indicators remain warm indicators, even if the CD is spatially shifted. As the relative palaeoclimatic curves can thus be comparatively robust to CD selection, this opens up the possibility of a new, more robust approach for estimating absolute values. The reconstructed values can first be expressed as deviations from the reconstructed value for the most recent sample, or from the mean of a few most recent samples (Birks et al., 2010; Birks and Seppä, 2004). These deviations from the modern reconstructed value can then be added to the modern observed value to estimate pseudo-absolute palaeoclimatic values. This method thus shifts the entire curve, fixing the top-most reconstructions with the observed value, thus adjusting for possible reconstruction-wide bias introduced by CD location. A limitation of this correction is that the bias is assumed to remain constant through time, which may not be true especially for reconstructions covering long timescales (cf. Jackson et al., 2009; Salonen et al., 2013). Our proposed correction procedure is analogous to an approach that has been used for many years in climate modelling. In the climate modelling application, one first calculates the difference between the modelled past or future climate and the modelled modern climate, and these anomalies are then added to the modern observed climate, thus debiasing the past or future climate projection (cf. Tabor and Williams, 2010). While this debiasing method has not, to our knowledge, been used when expressing reconstructions from fossil proxies, our results suggest that it may be valuable in increasing the robustness of quantitative palaeoclimatic reconstructions.