Terrestrial biosphere changes over the last 120 kyr

Abstract

A new global synthesis and biomization of long (> 40 kyr) pollen-data records is presented and used with simulations from the HadCM3 and FAMOUS climate models and the BIOME4 vegetation model to analyse the dynamics of the global terrestrial biosphere and carbon storage over the last glacial–interglacial cycle. Simulated biome distributions using BIOME4 driven by HadCM3 and FAMOUS at the global scale over time generally agree well with those inferred from pollen data. Global average areas of grassland and dry shrubland, desert, and tundra biomes show large-scale increases during the Last Glacial Maximum, between ca. 64 and 74 ka BP and cool substages of Marine Isotope Stage 5, at the expense of the tropical forest, warm-temperate forest, and temperate forest biomes. These changes are reflected in BIOME4 simulations of global net primary productivity, showing good agreement between the two models. Such changes are likely to affect terrestrial carbon storage, which in turn influences the stable carbon isotopic composition of seawater as terrestrial carbon is depleted in${}^{13}$C.

Introduction

Variations in global climate on multi-millennial timescales have caused substantial changes to terrestrial vegetation distribution, productivity, and carbon storage. Periodic variations in the Earth's orbital configuration (axial tilt with a $\sim$ 41 kyr period, precession with $\sim$ 19 and 23 kyr periods, and eccentricity with $\sim$ 100 kyr and longer periods) result in small variations in the seasonal and latitudinal distribution of insolation, amplified by feedback mechanisms (Berger, 1978). For the last $\sim$ 0.8 million years, long glacial periods have been punctuated by short interglacials on roughly a 100 kyr cycle. Glacial periods are associated with low atmospheric CO${}_{\mathrm{2}}$ concentrations, lowered sea level, and extensive continental ice sheets; interglacial periods are associated with high (similar to pre-industrial) CO${}_{\mathrm{2}}$ concentrations, high sea level, and reduced ice sheets (Petit et al., 1999; Peltier et al., 2004; Lüthi et al., 2008).

During glacial–interglacial cycles the productivity and carbon storage of the terrestrial biosphere are influenced by orbitally forced climatic changes and atmospheric CO${}_{\mathrm{2}}$ concentrations. Expansion of ice sheets during glacial periods caused a significant loss of land area available for colonization, but this was largely compensated for by the exposure of continental shelves due to lower sea level. The terrestrial biosphere (vegetation and soil) is estimated to contain around 2000 PgC (Prentice et al., 2001) plus a similar quantity stored in peatlands and permafrost (Ciais et al., 2012). During the last glacial period the terrestrial biosphere was significantly reduced. It has been estimated that the terrestrial biosphere contained 300 to 700 PgC less carbon during the Last Glacial Maximum (LGM, 21 ka BP) compared with pre-industrial times (Bird et al., 1994; Ciais et al., 2012; Crowley et al., 1995; Duplessy et al., 1988; Gosling and Holden, 2011; Köhler and Fischer, 2004; Prentice et al., 2011). As first noted by Shackleton et al. (1977), the oceanic inventory of carbon isotopes (${\mathit{\delta }}^{13}$C) is influenced by terrestrial carbon storage because terrestrial organic carbon has a negative signature, due to isotopic discrimination during photosynthesis. Many of the estimates of the reduction in terrestrial carbon storage at the LGM have therefore been based on the observed LGM lowering of deep-ocean ${\mathit{\delta }}^{13}$C. A reduction in the terrestrial biosphere of this size would have contributed a large amount of CO${}_{\mathrm{2}}$ to the atmosphere, although ocean carbonate compensation would have reduced the expected CO${}_{\mathrm{2}}$ increase to 15 ppm over about 5 to 10 kyr (Sigman and Boyle, 2000).

Many palaeoclimate data and modelling studies have focused on the contrasts between the LGM, the mid-Holocene (6 ka BP), and the pre-industrial period. The BIOME 6000 project (http://www.bridge.bris.ac.uk/resources/Databases/BIOMES_data) synthesized palaeovegetation records from many sites to provide global data sets for the LGM and mid-Holocene. Data syntheses are valuable in allowing researchers to see the global picture from scattered, individual records, and to enable model–data comparisons. The data can be interpreted in the context of a global, physically based model that allows the point-wise data to be seen in a coherent way. There are continuous, multi-millennial pollen records that stretch much further back in time than the LGM but they have not previously been brought together in a global synthesis to study changes of the last glacial–interglacial cycle. These records can provide a global picture of transient change in the biosphere and the climate system. Here we have synthesized and biomized (Prentice et al., 1996) a number of these records (for locations see Fig. 1), providing a new data set of land biosphere change that covers the last glacial–interglacial cycle. In Sect. 2.1 we outline the biomization procedures applied to reconstruct land biosphere changes.

Locations and altitudes of pollen records superimposed on pre-industrial HadCM3 orography (m).

[Figure omitted. See PDF]

To improve understanding of land biosphere interactions with the ocean–atmospheric reservoir, we have modelled the terrestrial biosphere for the last 120 kyr, from the previous (Eemian) interglacial to the pre-industrial period. Details of the atmosphere–ocean general circulation model (AOGCM) climate and vegetation model simulations are provided in Sect. 2.2. In Sect. 3 we evaluate biome reconstructions based on our model outputs using the BIOME 6000 project (www.bridge.bris.ac.uk/resources/Databases/BIOMES_data), and our new biomized synthesis of terrestrial pollen-data records, focusing on the pre-industrial period and 6 (mid-Holocene), 21 (LGM), 54 (a relatively warm interval in the last glacial period), 64 (a relatively cool interval in the glacial period), 84 (the early part of the glacial cycle), and 120 ka BP (the Eemian interglacial). The effects of rapid millennial-scale climate fluctuations were not simulated. Finally, in Sect. 4 we use our biome simulations to estimate net primary productivity.

Methods

Biomization

Biomization assigns pollen taxa to one or more plant functional types (PFTs). The PFTs are assigned to their respective biomes and affinity scores are calculated for each biome (sum of the square roots of pollen percentages contributed by the PFTs in each biome). This method was first developed for Europe (Prentice et al., 1996) and versions of it have been applied to most regions of the world (Jolly et al., 1998; Takahara et al., 1999; Tarasov et al., 2000; Thompson and Anderson, 2000; Williams et al., 2000; Elenga et al., 2004; Pickett et al., 2004; Marchant et al., 2009). We apply these regional PFT schemes (Table 1) to pollen records that generally extend > 40 kyr, assigning the pollen data to megabiomes (tropical forest, warm-temperate forest, temperate forest, boreal forest, savanna/dry woodland, grassland/dry shrubland, desert, and tundra) as defined by Harrison and Prentice (2003) in order to harmonize regional variations in PFT to biome assignments and to allow globally consistent model–data comparisons.

Details of the various biomization schemes applied for the different regions.

Africa	Jolly et al. (1998)
Southeast Asia, Australia	Pickett et al. (2004)
Japan	Takahara et al. (1999)
Southern Europe	Elenga et al. (2004)
Northeastern Europe	Tarasov et al. (2000)
North America: western north	Thompson and Anderson (2000)
North America: east and northeast	Williams et al. (2000)
Latin America	Marchant et al. (2009)

Table 2 lists the pollen records used. Biomization matrices and megabiome score data can be found in the Supplement. For taxa with no PFT listing, the family PFT was used if part of the regional biomization scheme. Plant taxonomy was checked using itis.gov, tropicos.org, and the African Pollen Database. Pollen taxa can be assigned to more than one PFT either because they include several species in the genus or family, with different ecologies, or because they comprise species that can adopt different habitats in different environments.

Age models provided with the individual records were used. However, in cases where radiocarbon ages were only provided for specific depths (e.g. Mfabeni, CUX), linear interpolations between dates were used to estimate ages for the remaining depths. Some age models may be less certain, especially at sites which experience variable sedimentation rates and/or erosion. Sometimes more than one age model accompanies the data, illustrating the range of ages and also that there can be large uncertainties. To aid comparison, for several southern European sites (e.g. Italy and Greece) it has been assumed that vegetation changes occurred synchronously within the age uncertainties of their respective chronologies, for which there is evidence (e.g. Tzedakis et al., 2004a).

Model simulations

Global simulations of vegetation changes over the last glacial cycle were produced using a vegetation model (BIOME4) forced offline using previously published climate simulations from two AOGCMs (HadCM3 and FAMOUS). By using two models we test the robustness of the reconstructions to different climate forcings.

HadCM3

HadCM3 is a general circulation model, consisting of coupled atmospheric model, ocean, and sea ice models (Gordon et al., 2000; Pope et al., 2000). The resolution of the atmospheric model is 2.5${}^{\circ }$ in latitude by 3.75${}^{\circ }$ in longitude by 19 unequally spaced levels in the vertical. The resolution of the ocean is 1.25 by 1.25${}^{\circ }$ with 20 unequally spaced layers in the ocean extending to a depth of 5200 m. The model contains a range of parameterizations, including a detailed radiation scheme that can represent the effects of minor trace gases (Edwards and Slingo, 1996). The land surface scheme used is the Met Office Surface Exchange Scheme 1 (MOSES1; Cox et al., 1999). In this version of the model, interactive vegetation is not included. The ocean model uses the Gent–McWilliams mixing scheme (Gent and McWilliams, 1990), and sea ice is a thermodynamic scheme with parameterization of ice drift and leads (Cattle and Crossley, 1995).

Multiple “snapshot” simulations covering the last 120 kyr have been performed with HadCM3. The boundary conditions and setup of the original set of simulations have been previously documented in detail in Singarayer and Valdes (2010). The snapshots were done at intervals of every 1 ka between the pre-industrial (PI) and LGM (21 ka BP), every 2 ka between the LGM and 80 ka BP, and every 4 ka between 80 and 120 ka BP. Boundary conditions are variable between snapshots but constant for each simulation. Orbital parameters are taken from Berger and Loutre (1991). Atmospheric concentrations of CO${}_{\mathrm{2}}$ were taken from a stacked ice core record of Vostok (Petit et al., 1999) prior to 62 kyr, incorporating Taylor Dome (Indermühle et al., 2000) to 22 kyr and EDC96 (Monnin et al., 2001) up to 0 kyr. CH${}_{\mathrm{4}}$ and N${}_{\mathrm{2}}$O were taken from EPICA (Spahni et al., 2005; Loulergue et al., 2008), and all greenhouse gas concentrations were on the EDC3 timescale (Parrenin et al., 2007). The prescription of ice sheets was achieved with ICE-5G (Peltier, 2004) for 0–21 ka BP, and extrapolated to the pre-LGM period from the ICE-5G reconstruction using the method described in Eriksson et al. (2012). The simulations were each spun up from the end of previous runs described in Singarayer and Valdes (2010) to adjust to the modified ice-sheet boundary conditions for 470 years. The monthly climatologies described hereafter are of model years 470–499. The model performs reasonably well in terms of glacial–interglacial global temperature anomaly (HadCM3 is in the middle of the distribution of global climate models and palaeoclimate reconstructions) and high-latitude temperature trends (although as with all models, the magnitude of the temperature anomalies in the glacial is underestimated), as well as at lower latitudes (Singarayer and Valdes, 2010; Singarayer and Burrough, 2015).

FAMOUS

FAMOUS (Smith, 2012) is an Earth System model, derived from HadCM3. It is run at approximately half the spatial resolution of HadCM3 to reduce the computational expense associated with AOGCM simulations without fundamentally sacrificing the range of climate system feedbacks of which it is capable. Pre-industrial control simulations of FAMOUS have both an equilibrium climate and global climate sensitivity similar to that of HadCM3. A suite of transient FAMOUS simulations of the last glacial cycle, conducted with specified atmospheric CO${}_{\mathrm{2}}$, ice sheets, and changes in solar insolation resulting from variation in the Earth's orbit, compare well with the NGRIP, EPICA, and MARGO proxy reconstructions of glacial surface temperatures (Smith and Gregory, 2012). For the present study, we use the most realistically forced simulation of the Smith and Gregory (2012) suite (experiment ALL-ZH), forced with Northern Hemisphere ice sheets taken from the physical ice-sheet modelling work of Zweck and Huybrechts (2005), atmospheric CO${}_{\mathrm{2}}$, CH${}_{\mathrm{4}}$, and N${}_{\mathrm{2}}$O concentrations from the EPICA project (Lüthi et al., 2008 and Spahni et al., 2005, mapped onto the EDC3 Parrenin et al., 2007, age scale) and orbital forcing from Berger (1978). The composite CO${}_{\mathrm{2}}$ record contained in Lüthi et al. (2008) uses data from the Vostok core (Petit et al., 1999) between 22 and 393 kyr. The Vostok record is now believed (Bereiter et al., 2012) to be erroneously low during the early part of Marine Isotope Stage 3. For this reason, the FAMOUS results during this period are likely biased too cold. Although of a lower spatial resolution than HadCM3, these FAMOUS simulations have the benefit of being transient and representing low-frequency variability within the climate system, as well as using more physically plausible ice-sheet extents before the LGM than were used in the HadCM3 simulations. To allow the transient experiments to be conducted in a tractable amount of time, these forcings were all “accelerated” by a factor of 10, so that the 120 kyr of climate are simulated in 12model kyr – this method has been shown to have little effect on the surface climate (Timm and Timmerman, 2007; Ganopolski et al., 2010), although it does distort the response of the deep ocean. In addition, we did not include changes in sea level, Antarctic ice volume, or meltwater from ice sheets to enable the smooth operation of the transient simulations. The impact of ignoring the continental shelves exposed by lower sea levels will be discussed later; the latter two approximations are unlikely to have an impact over the timescales considered here. Although within the published capabilities of the model, interactive vegetation was not used during this simulation, with (ice sheets aside) the land surface characteristics of the model being specified as for a pre-industrial simulation.

BIOME4

BIOME4 (Kaplan et al., 2003) is a biogeochemistry–biogeography model that predicts the global vegetation distribution based on monthly mean temperature, precipitation, and sunshine fraction, as well as information on soil texture, depth, and atmospheric CO${}_{\mathrm{2}}$. It derives a seasonal maximum leaf area index that maximizes NPP for a given PFT by simulating canopy conductance, photosynthesis, respiration, and phenological state. Model grid boxes are then assigned biome types based on a set of rules that use dominant and sub-dominant PFTs, as well as environmental limits.

Two reconstructions of the evolution of the climate over the last glacial cycle were obtained by calculating monthly climate anomalies with respect to the simulated pre-industrial for the HadCM3 and FAMOUS glacial climate simulations, respectively, then adding these anomalies, on the native FAMOUS and HadCM3 grids, to an area averaged interpolation of the Leemans and Cramer (1991) observed climatology provided with the BIOME4 distribution. These climate reconstructions were then used to force two BIOME4 simulations. The climate anomaly method allows us to correct for known systematic errors in the climates of HadCM3 and FAMOUS and produce more accurate results from BIOME4, although the method assumes that the pre-industrial errors in each model are systematically present and unchanged over ice-free regions throughout the whole glacial cycle. We chose to use the actual climate model grids for the BIOME4 simulations, rather than interpolating onto the higher-resolution observational climatology grid, to avoid concealing the significant impact that the climate model resolution has on the vegetation simulation, and to highlight the differences between the physical representation of the climate between the two different models. Because of its lower resolution, FAMOUS cannot represent geographic variation at the same scale as HadCM3, which affects not only the areal extent of individual biomes but also how altitude is represented in the model, which can have a significant effect on the local climate and resulting biome affinity. The frequency of data available from the FAMOUS run also limits the accuracy of the minimum surface air temperature it can force BIOME4 with, as only monthly average temperatures were available. This results in some aspects of the FAMOUS-forced BIOME4 simulation seeing a less extreme climate than it should, and may artificially favour more temperate vegetation in some locations.

Soil properties on exposed shelves were extrapolated from the nearest pre-industrial land points. There is no special correction for the input climate anomalies over this exposed land, which results in a slightly subdued seasonal cycle at these points (due to smaller inter-seasonal variation of ocean temperatures). The version of the observational climatology distributed with BIOME4 includes climate values for these areas. The BIOME4 runs used the time-varying CO${}_{\mathrm{2}}$ records that were used to force the corresponding climate models, as described in Sects. 2.2.1 and 2.2.2. As well as affecting productivity, the lower CO${}_{\mathrm{2}}$ concentrations found during the last glacial favour the growth of plants that use the C${}_{\mathrm{4}}$ photosynthetic pathway (Ehleringer et al., 1997), which can affect the distribution of biomes as well. All other BIOME4 parameters as well as soil characteristics were held constant at pre-industrial values.

The results of the HadCM3-forced BIOME4 simulation will be referred to in this paper as B4H, and those from the FAMOUS-forced BIOME4 simulation as B4F.

Results

In this section, the results of both the pollen-based biomization for individual regions and the biome reconstructions based on the GCM climate simulations will be outlined. The biomized records and biomization matrix can be found in the Supplement. Biome changes relating to millennial-scale climate oscillations are discussed elsewhere (e.g. Harrison and Sanchez Goñi, 2010, and references therein).

Biomization

This method translates fossil pollen assemblages into a form that allows direct data–model comparison and allows the reconstruction of past vegetation conditions. Biome affinity scores for each location are shown in the Supplement.

North America

Two regional PFT schemes were used for sites from North America: the scheme of Williams et al. (2000) for northern and eastern North America and the scheme of Thompson and Anderson (2000) for the western USA. For their study of biome response to millennial climate oscillations between 10 and 80 ka BP, Jiménez-Moreno et al. (2010) applied one scheme for the whole of North America, with a subdivision for southeastern pine forest. All biomization matrices and scores for individual sites used in our study, generally at 1 kyr resolution, as well as explanatory files can be found in the Supplement. The Arctic Baffin Island sites (Amarok and Brother of Fog) have highest affinity scores for tundra during the ice-free Holocene and last interglacial.

At Lake Tulane (Florida) the grassland and dry shrubland biome has the highest affinity scores for the last 52 kyr, apart from two short intervals ($\sim$ 14.5 to 15.5 ka BP and $\sim$ 36.5 to 37.5 ka BP) where warm-temperate forest and temperate forest have highest scores. According to Williams et al. (2000), present day, 6 ka BP, and LGM records of most of Florida and the southeast of the USA should be characterized by highest affinity scores for the warm-temperate forest biome (Williams et al., 2000). The discrepancy in our biomization results with those of the regional biomization results of Williams et al. (2000) is due to high percentages of Quercus, Pinus undiff. (both are in the grassland and dry shrubland and warm-temperate forest biomes), and Cyperaceae and Poaceae that contribute to highest affinity scores of the grassland and shrubland biome. Interestingly, the temperate forest biome has highest affinity scores in a short interval ($\sim$ 15 ka BP) during the deglaciation. In Jiménez-Morene et al. (2010) Pinus does not feature in the grassland and dry shrubland biome, but comprises a major component of the southeastern pine forest; hence their biomized Lake Tulane record fluctuates between the “grassland and dry shrubland” biome and “southeastern pine forest biome”.

In western North America pollen data from San Felipe (16 to 47 ka BP), Potato Lake (last 35 kyr), and Bear Lake (last 150 kyr) all show highest scores for the grassland and dry shrubland biome. Potato Lake is currently situated within a forest (Anderson, 1993). In our biomizations Pinus pollen equally contribute to scores of boreal forest, temperate forest, warm-temperate forest, and the grassland and dry shrubland biomes. In addition, high contributions of Poaceae occur, so that the grassland and dry shrubland biome has highest affinity scores throughout the last 35 kyr. Again, in the Jiménez-Morene et al. (2010) biomizations Pinus does not feature in the grassland and dry shrubland biome and hence the forest biomes have highest affinity scores in their biomizations. At Carp Lake the Holocene is characterized by alternating highest affinity scores between the temperate forest and grassland and dry shrubland biomes, whereas during the glacial only the grassland and dry shrubland biome attains highest affinity scores. The age model of Carp Lake suggests this record goes back to the Eemian, and if so, the last interglacial climate was lacking the alternation between the temperate forest and grassland and dry shrubland biomes as found during the late Holocene. Modern and LGM biomizations at Carp Lake and Bear Lake are similar to those of Thompson and Anderson (2000). Biomizations for Carp Lake between 10 and 80 ka BP by Jiménez-Morene et al. (2010) generally look similar to ours, apart from 36, 57–70, and 72–80 ka BP, where the temperate forest biome shows highest affinity scores because Pinus undiff. is treated as insignificant in their biomization. Biomizations of Bear Lake between 10 and 80 ka BP are similar to Jiménez-Morene et al. (2010).

Details of the locations of pollen-data records synthesized in this study.

	Core	Latitude	Longitude	Altitude (m a.s.l.)	Age $\sim$/(ka BP)	Reference	Biomization reference
North America
Canada (short)	Brother of Fog	67.18	$-$63.25	380	Last interglacial	Frechette et al. (2006)	Williams et al. (2000)
Canada (short)	Amarok	66.27	$-$65.75	848	Holocene and lastinterglacial	Frechette et al. (2006)	Williams et al. (2000)
USA	Carp Lake	45.92	$-$120.88	714	0 to ca. 130	Whitlock and Bartlein (1997)	Thompson and Anderson (2000)
USA	Bear Lake	41.95	$-$111.31	1805	0 to 150	Jiménez-Moreno et al. (2007)	Thompson and Anderson (2000)
USA	Potato lake	34.4	$-$111.3	2222	2 to ca. 35	Anderson et al. (1993)	Thompson and Anderson (2000)
USA	San Felipe	31	$-$115.25	400	16 to 42	Lozano-Garcia et al. (2002)	Thompson and Anderson (2000)
USA	Lake Tulane	27.59	$-$81.50	36	0 to 52	Grimm et al. (2006)	Williams et al. (2000)
Latin America
Mexico	Lake Patzcuaro	19.58	$-$101.58	2044	3 to 44	Watts and Bradbury (1982)	Marchant et al. (2009)
Guatemala	Lake Petén-Itzá	16.92	$-$89.83	110	0 to 86	Correa-Metrio et al. (2012)	Marchant et al. (2009)
Colombia	Ciudad Universitaria X	$-$4.75	$-$74.18	2560	0 to 35	van der Hammenand González (1960)	Marchant et al. (2009)
Peru	Laguna Junin	$-$11.00	$-$76.18	4100	0 to 36 (LAPD1?)	Hansen et al. (1984)	Marchant et al. (2009)
Peru/Bolivia	Lake Titicaca	$-$15.9	$-$69.10	3810	3 to 370(shown until 140)	Gosling et al. (2008),Hanselman et al. (2011),Fritz et al. (2007)	Marchant et al. (2009)
Guatemala	Lago Quexil	16.92	$-$89.88	110	9 to 36	Leyden (1984),Leyden et al. (1993, 1994)	Marchant et al. (2009)
Brazil	Salitre	$-$19.00	$-$46.77	970	2 to 50 (LAPD1)	Ledru (1992, 1993),Ledru et al. (1994, 1996)	Marchant et al. (2009)
Brazil	Colonia	$-$23.87	$-$46.71	900	0 to 120	Ledru et al. (2009)	Marchant et al. (2009)
Brazil	Cambara	$-$29.05	$-$50.10	1040	0 to 38	Behling et al. (2004)	Marchant et al. (2009)
Peru/Bolivia	Lake Titicaca	$\sim$ $-$16 to $-$17.5	$\sim$ $-$68.5 to $-$70	3810	3 to 138	Hanselman et al. (2011),Fritz et al. (2007)	Marchant et al. (2009)
Bolivia	Uyuni	$-$20.00	$-$68.00	653	17 to 108	Chepstow Lusty et al. (2005)	Marchant et al. (2009)
Europe
Russia	Butovka	55.17	36.42	198	Holocene, earlyglacial and Eemian	Borisova (2005)	Tarasov et al. (2000)
Russia	Ilinskoye	53	37	167	early glacial & Eemian	Grichuk et al. (1983),Velichko et al. (2005)	Tarasov et al. (2000)
Poland	Horoszki Duze	52.27	23		$\sim$ 75 to Eemian	Granoszewski (2003)	Tarasov et al. (2000)
Germany	Klinge	51.75	14.51	80	early glacial,Eemian & Saalian(penultimate glacial)	Novenko et al. (2008)	Tarasov et al. (2000)
Germany	Füramoos	47.59	9.53	662	0 to 120	Müller et al. (2003)	Prentice et al. (1992)
Germany	Jammertal	48.10	9.73	578	Eemian	Müller (2000)	Prentice et al. (1992)
Germany	Samerberg	47.75	12.2	595	Eemian and earlyWürmian	Grüger (1979a, b)	Prentice et al. (1992)
Germany	Wurzach	47.93	9.89	650	Eemian and earlyWürmian	Grüger and Schreiner (1993)	Prentice et al. (1992)
Italy	Lagaccione	42.57	11.85	355	0 to 100	Magri (1999)	Elenga et al. (2004)
Italy	Lago di Vico	42.32	12.17	510	0 to 90	Magri and Sadori (1999)	Elenga et al. (2004)
Italy	Valle di Castiglione	41.89	12.75	44	0 to 120	Magri and Tzedakis (2000)	Elenga et al. (2004)
Italy	Monticchio	40.94	15.60	656	0 to 120	Allen et al. (1999)	Elenga et al. (2004)
Greece	Ioannina	39.76	20.73	470	0 to 120	Tzedakis et al. (2002, 2004b)	Elenga et al. (2004)
Greece	Tenaghi Philippon	41.17	24.30	40	0 to 120	Wijmstra (1969),Wijmstra and Smith (1976),Tzedakis et al. (2006)	Elenga et al. (2004)
Africa
Uganda	ALBERT-F	1.52	30.57	619	0 to 30	Beuning et al. (1997)	Jolly et al. (1998)
Uganda	Mubwindi Swamp 3	$-$1.08	29.46	2150	0 to 40	Marchant et al. (1997)	Jolly et al. (1998)
Rwanda	Kamiranzovy Swamp 1	$-$2.47	29.12	1950	13 to 40	Bonnefille and Chalie (2000)	Jolly et al. (1998)
Burundi	Rusaka	$-$3.43	29.61	2070	0 to 47	Bonnefille and Chalie (2000)	Jolly et al. (1998)
Burundi	Kashiru Swamp A1	$-$3.45	29.53	2240	0 to 40	Bonnefille and Chalie (2000)	Jolly et al. (1998)
Burundi	Kashiru Swamp A3	$-$3.45	29.53	2240	0 to 40	Bonnefille and Chalie (2000)	Jolly et al. (1998)
Tanzania	Uluguru	$-$7.08	37.62	2600	0 to > 45	Finch et al. (2009)	Jolly et al. (1998)
Madagascar	Lake Tritrivakely	$-$19.78	46.92	1778	0 to 40	Gasse and Van Campo (1998)	Jolly et al. (1998)
South Africa	Tswaing (Saltpan) Crater	$-$25.57	28.07	1100	0 to 120 (althoughafter 35 probablyless secure based)	Scott (1988),Partridge et al. (1993),Scott (1999a, 1999b)	Jolly et al. (1998)
South Africa	Mfabeni Swamp	$-$28.13	32.52	11	0 to 43	Finch and Hill (2008)	Jolly et al. (1998)
Asia/Australasia
Russia	Lake Baikal	53.95	108.9		114 to 130
Japan	Lake Biwa	35	135	85.6	0 to 120	Nakagawa et al. (2008)	Takahara et al. (1999)
Japan	Lake Suigetsu	35.58	135.88	$\sim$ 0	0 to 120	Nakagawa et al. (2008)	Takahara et al. (1999)
Thailand	Khorat Plateau	17	103	$\sim$ 180	0 to 40	Penny (2001)	Pickett et al. (2004)
Australia	Lynch's Crater	$-$17.37	145.7	760	0 to 120	Kershaw (1986)	Pickett et al. (2004)
New Caledonia	Xero Wapo	$-$22.28	166.97	220	0 to 120	Stevenson and Hope (2005)	Pickett et al. (2004)
Australia	Caldeonia Fen	$-$37.33	146.73	1280	0 to 120	Kershaw et al. (2007)	Pickett et al. (2004)
New Zealand	Okarito	$-$43.24	170.22	70	0 to 120	Vandergoes et al. (2005)	Pickett et al. (2004)

Latin America

The regional biomization scheme of Marchant et al. (2009) was used for Latin American locations. Hessler et al. (2010) discuss the effects of millennial climate variability on the vegetation of tropical Latin America and Africa between 23${}^{\circ }$ N and 23${}^{\circ }$ S, using similar biomization schemes. In our study eleven sites from Central and South America are considered covering a latitudinal gradient of 49${}^{\circ }$ (from 20${}^{\circ }$ to $-$29${}^{\circ }$) and an elevation range of 3900 m (from 110 to 4010 m a.s.l.) (Table 2). Five of the sites are from relatively low elevations (< 1500 m a.s.l.); from north to south these are Lago Quexil and Petén-Itzá in Guatemala and Salitre, Colonia, and Cambara in southeastern Brazil. The high-elevation records (> 1500 m a.s.l.), with the exception of the most northerly site in Mexico (Lake Patzcuaro), are distributed along the Andean chain: Ciudad Universitaria X (Colombia), Laguna Junin (Peru), Lake Titicaca (Bolivia/Peru), and Salar de Uyuni (Bolivia).

The five lowland sites indicate the persistence of forest biomes for much of the last 130 kyr. In Central America, the Lago Quexil record stretches back to 36 ka BP and has highest affinity scores for the warm-temperate forest biome during the early Holocene. During glacial times the temperate forest biome dominates, intercalated with mainly the grassland and dry shrubland and desert biomes during the LGM and last deglaciation. At Lago Petén-Itzá (also Guatemala) highest affinity scores for the warm-temperate forest biome are recorded for the last 86 kyr. The Salitre and Colonia records are the only Latin American sites that fall within the tropical forest biome today. The majority of the Salitre record shows high affinities for tropical forest from $\sim$ 64 ka BP to present day, apart from an interval coinciding with the Younger Dryas which displays highest affinity scores for the warm-temperate forest biome. The southern-most Brazilian record, at Colonia, has highest affinity scores for tropical forest for the last 40 kyr, except between 28 and 21 ka BP ($\sim$ coincident with the LGM), when scores were highest for the warm-temperate forest biome. Between 120 and 40 ka BP, highest affinity scores alternate between the tropical forest and warm-temperate forest biome at Colonia. The biomized Colonia record of Hessler et al. (2010) generally shows the same features, apart from an increase in affinity scores of the drier biomes between 10 and 18 ka BP. To the south, at Cambara (Brazil), highest affinity scores are found for warm-temperate forest during the Holocene and between 38 and 29 ka BP, whilst during the interval in between they alternate between warm-temperate forest and grassland and dry shrubland.

Apart from Laguna Junin, higher-elevation sites (> 1500 m: Lake Patzcuaro, Titicaca, Uyuni, and CUX) do not show a strong glacial–interglacial cycling in their affinity scores; Mexican site Lake Patzcuaro (2240 m) and Colombian site CUX (2560 m) have highest affinity scores mainly for warm-temperate forest over the last 35 kyr, although they alternate between warm-temperate forest and temperate forest during the Holocene and at CUX also during the LGM. Lake Patzcuaro and CUX biomization results for the Holocene, 6 ka BP, and LGM compare well with those derived by Marchant et al. (2009). At Uyuni (3643 m) highest affinity scores are for temperate forest and grassland and dry shrubland biome between 108 and 18 ka BP. At Titicaca (3810 m) high affinity scores are found for temperate forest over the last 130 kyr, apart from during the previous interglacial (Eemian), when highest affinity scores for the desert biome occur. Finally, at Lago Junin highest affinity scores alternate between warm-temperate forest and temperate forest during the Holocene and temperate forest and grassland and dry shrubland during the glacial period.

Africa

For the biomization of African pollen records the scheme of Elenga et al. (2004) was applied. What is specifically different from southern European biomizations is that Cyperaceae are not included as this taxon generally occurs in high abundances in association with wetland environments, where they represent a local signal (Elenga et al., 2004). It is noted that most African sites are from highland or mountain settings, with the exception of Mfabeni (11 m a.s.l.).

At the mountain site Kashiru Swamp in Burundi the Holocene is characterized by an alternation of highest affinity scores for tropical forest, warm-temperate forest, and the grassland and dry shrubland biomes. During most of the glacial, scores are highest for the grassland and dry shrubland biome, preceded by an interval where warm-temperate forest showed highest scores. Our results are similar to those obtained by Hessler et al. (2010). Highest affinity scores for tropical forest and warm-temperate forest are found during the Holocene at the Rusaka Burundi mountain site, whereas those of the last glacial again have highest scores for the grassland and dry shrubland biome. At the Rwandan Kamiranzovy site the grassland and dry shrubland biome display highest scores during the last glacial (from $\sim$ 30 ka BP) and deglaciation, occasionally alternating with the warm-temperate forest biome. In Uganda at the low mountain site Albert F (619 m), the Holocene and potentially Bølling–Allerød is dominated by highest affinity scores for tropical forest, whereas the Younger Dryas and last glacial show highest affinity scores for the grassland and dry shrubland biome. In the higher-elevation Ugandan mountain site Mubwindi Swamp (2150 m), the Holocene pollen record shows alternating highest affinity scores between tropical forest and the grassland and dry shrubland biome, whereas the glacial situation is similar to the Albert F site (e.g. dominated by highest scores for the grassland and dry shrubland biome). In South Africa, the Mfabeni Swamp record shows highest affinity scores for the grassland and dry shrubland biome for the last 46 kyr, occasionally alternated with the savanna and dry woodland and tropical forest biome during the late Holocene. At the Deva Deva Swamp in the Uluguru Mountains highest affinity scores are for the grassland and dry shrubland biome for the last $\sim$ 48 kyr. At Tswaing Crater the grassland and dry shrubland biome dominates throughout the succession, including the Holocene and glacial. At Lake Tritrivakely (Madagascar) the grassland and dry shrubland biome dominates, apart from between 3 and 0.6 ka BP, when the tropical forest biome shows highest affinity scores. Our results compare well with those of Elenga et al. (2004), who show a LGM reduction in tropical rainforest and lowering of mountain vegetation zones in major parts of Africa.

Europe

For European pollen records three biomization methods were used that are region-specific. For southern Europe the biomization scheme of Elenga et al. (2004) was used, where Cyperaceae are included in the biomization as they can occur as an “upland” species characteristic of tundra. For sites from the Alps the biomization scheme of Prentice et al. (1992) was used, and for northern European records the biomization scheme of Tarasov et al. (2000). Fletcher et al. (2010) use one uniform biomization scheme to discuss millennial climate in European vegetation records between 10 and 80 ka BP.

In southern Europe at the four Italian sites (Monticchio, Lago di Vico, Lagaccione, and Valle di Castiglione) the Holocene and last interglacial show highest affinity scores for warm-temperate forest and temperate forest biomes. During most of the glacial and also cold interglacial substages the grassland and dry shrubland biome has highest affinity scores, whereas during warmer interstadial intervals of the last glacial the temperate forest biome had highest affinity scores. At Tenaghi Phillipon and Ioannina a similar biome sequence may be observed, with highest affinity scores for temperate forest and warm-temperate forest biomes during interglacials. During the last glacial and cool substages of the previous interglacial the grassland and dry shrubland biome showed highest affinity scores at Tenaghi Philippon. At Ioannina the LGM and last glacial cool stadial intervals have highest affinity scores for grassland and dry shrubland, whereas affinity scores of glacial interstadial periods are highest for temperate forest. Our biomization results for southern European sites agree well with those of Elenga et al. (2004), who also found a shift to drier grassland and dry shrubland biomes during glacial times. Instead of a desert and tundra biome, Fletcher et al. (2010) define a xyrophytic steppe and eurythermic conifer biome in their biomizations for Europe, giving subtle differences in the biomization records, with the Fletcher et al. (2010) biomized records showing an important contribution of affinity scores to the xerophytic steppe biome. Characteristic species for the xerophytica steppe biome include Artemisia, Chenopodiaceae and Ephedra, which in the southern Europe biomization scheme of Elenga et al. (2000) feature in the desert biome and grassland and dry shrubland biome (only Ephedra).

All four alpine sites are from altitudes between 570 and 670 m and for all four sites the last interglacial period was characterized by having highest scores for the temperate forest biome. At Füramoos the last glacial showed highest affinity scores for the tundra biome, whilst during the Holocene the temperate forest biome shows highest affinity scores. In the Fletcher scheme, characteristic pollen for the eurythermic conifer biome includes Pinus and Juniperus. In our biomization Pinus and Juniperus contribute to all biomes except for the desert and tundra biomes.

Most northern European sites are mainly represented for the last interglacial period, apart from Horoszki Duze in Poland. At most sites the temperate forest biome and boreal forest biome show highest affinity scores during the last interglacial (Eemian), whereas cool substages and early glacial (Butovka, Horoszki Duze) show high affinity scores for the grass and dry shrubland biome These results compare well with Prentice et al. (2000), who suggest a southward displacement of the Northern Hemisphere forest biomes and more extensive tundra- and steppe-like vegetation during the LGM.

Asia

For the higher latitude site Lake Baikal the biomization scheme of Tarasov et al. (2000) was used. For the two Japanese pollen sites we used the biomization scheme of Takahara et al. (1999). At Lake Baikal, during the Eemian the highest affinity scores are for the boreal and temperate forest biomes; the penultimate deglaciation and cool substage show highest affinity scores for the grassland and dry shrubland biome, similar to northern European sites. Pollen taxa such as Carpinus, Pterocarya, Tilia cordata, and Quercus have probably been redeposited or transported over a large distance; however, they all make up less than 1 % of the pollen spectrum and therefore did not influence the biomization much.

At Lake Suigetsu in Japan the warm-temperate forest biome shows highest affinity scores over the last 120 kyr; those of other biomes (including tundra) show increasing affinity scores during glacial times but never exceed those of the warm-temperate forest biome. At Lake Biwa the warm-temperate forest biome shows highest affinity scores during interglacial times, whilst in between they alternate between the warm-temperate forest biome and the temperate forest biome. These results agree well with those of Takahara et al. (1999) and Takahara et al. (2010).

East Asia/Australasia

For East Asian and Australasian sites the scheme of Pickett et al. (2004) was used. In Thailand the Khorat Plateau site shows highest affinity scores for the tropical forest biome over the last $\sim$ 40 kyr. At New Caledonia's Xero Wapa, the warm-temperate forest and tropical forest biomes show highest affinity scores over the last 127 kyr. In Australia's Caledonia Fen interglacial times (Holocene and previous interglacial) show highest affinity scores for the savanna and dry woodland biome. During the glacial the grassland and dry shrubland biome generally shows highest affinity scores, occasionally alternated with highest scores for the savanna and dry woodland biome during the early part of Marine Isotope Stage (MIS) 3 and what would be MIS 5a (ca. 80–85 ka BP). Over most of the last glacial–interglacial cycle highest affinity scores at Lynch's Crater are for the tropical forest and warm-temperate forest biomes. The savanna and dry forest biome becomes important during MIS 4 to 2 and generally shows highest affinity scores between 40 and 7 ka BP, probably as a result of increased biomass burning (human activities) causing the replacement of dry rainforest by savanna. In addition, the significance of what is considered to be tundra from MIS 4 is due to an increase in Cyperaceae with the expansion of swamp vegetation over what was previously a lake. At Okarito (New Zealand), the temperate forest biome has highest affinity scores throughout (occasionally alternated with warm-temperate forest), apart from during the LGM and deglaciation ($\sim$ 25–14 ka BP), where those of savanna and dry woodland as well as grassland and dry shrubland show highest affinity scores. Biomization results for the Australian mainland and Thailand agree well with those obtained by Pickett et al. (2004) for the Holocene and LGM.

HadCM3–FAMOUS model comparison

Although the source codes of HadCM3 and FAMOUS are very similar, differences in the resolution of the models and the setup of their simulations result in a number of differences in both the climates they produce and the vegetation patterns seen in B4H and B4F over the last glacial cycle. Specific regions and times where they disagree on the dominant biome type will be discussed later, but there are a number of features that apply throughout the simulations.

Both B4H and B4F keep the underlying soil types constant as for the pre-industrial throughout the glacial cycle. The HadCM3 snapshot simulations allowed for the exposure of coastal shelves as sea level changed through the glacial cycle, with reconstructions based on Peltier and Fairbanks (2006), who used the SPECMAP ${\mathit{\delta }}^{18}$O record (Martinson et al., 1987) to constrain ice volume/sea level change from the last interglacial to the LGM. FAMOUS, on the other hand, kept global mean sea level as for the present day throughout the whole transient simulation. As a consequence the area of land available to vegetation expands and contracts with falling and rising sea level in B4H but remains unchanged in B4F. Inclusion of changing land exposure with sea level therefore allows for significant additional vegetation changes as will be discussed further later.

Full details of the climates produced by FAMOUS and HadCM3 in these simulations can be found in Smith and Gregory (2012) and Singarayer and Valdes (2010). In general, land surface temperature anomalies in the HadCM3 simulations are a degree or so colder than in FAMOUS. This difference in temperature, present to some degree throughout most of the simulation, is attributed mainly to differences in surface height and ice-sheet extent, although differences in the CO${}_{\mathrm{2}}$ forcing play a role in MIS 3. FAMOUS model results are also, on average, slightly drier compared with those of HadCM3. This is additionally related to the model resolution, with HadCM3 showing much more regional variation (some areas become wetter and some drier), whilst FAMOUS produces a more spatially uniform drying as the climate cools. A notable exception to this general difference is in northwestern Europe, where FAMOUS more closely reproduces the temperatures reconstructed from Greenland ice cores (Masson-Delmotte et al., 2005), compared to the HadCM3 simulations used here which have a significant warm bias at the LGM. Millennial-scale cooling events and effects of ice rafting are not features of our model runs, which present a relatively temporally smoothed simulation of the last glacial cycle.

Data–model comparison

We present here an overview of the vegetation reconstructions for the last glacial–interglacial cycle simulated in B4H and B4F. We compare the simulated biomes in B4H and B4F with each other and with the dominant megabiome derived from the pollen-based biomizations, restricting our description of the results to major areas of agreement and disagreement. Maps of the dominant megabiomes produced by B4H and B4F with superimposed reconstructed dominant megabiomes for these periods are shown in Fig. 2.

We focus on a few specific periods, detailed below, since reviewing every detail present in this comparison is unfeasible. The pre-industrial period serves as a test bed to identify biases inherent in our model setup, before climate anomalies have been added. The 6 ka BP mid-Holocene period represents an orbital and ice-sheet configuration favouring generally warm Northern Hemisphere climate (Berger and Loutre, 1991). The LGM simulation at 21 ka BP is at the height of the last glacial cycle, when ice sheets were at their fullest extent, orbital insolation seasonality was similar to present and CO${}_{\mathrm{2}}$ was at its lowest concentration ($\sim$ 185 ppm), and the resulting climate was cold and dry in most regions. These three time periods form the basis of the standard PMIP2 simulations and were used in the BIOME 6000 project. We thus additionally compare our simulations with the BIOME 6000 results for these time periods. The 54 ka BP interval is representative of peak warm conditions during Marine Isotope Stage 3 (MIS 3), where both the model climates and some proxy evidence suggest relatively warm conditions, at least for Europe (Voelker et al., 2002), associated with temporarily higher levels of greenhouse gases, an orbital configuration that favours warmer Northern Hemisphere summers, and Northern Hemisphere ice sheet volume roughly half that of the LGM. The time slice 64 ka BP represents MIS 4, both greenhouse gases and Northern Hemisphere insolation were lower, and Northern Hemisphere ice volume was two-thirds higher than at 54 ka BP, resulting in significantly cooler global climate. The time slice 84 ka BP is representative of stadial conditions of the early part of the glacial (at the end of MIS 5), after both global temperatures and atmospheric concentrations of CO${}_{\mathrm{2}}$ have fallen significantly and the Laurentide ice sheet has expanded to a significant size but before the Fennoscandian ice sheet can have a major influence on climate. The 84 ka BP period can be compared with the Eemian (120 ka BP, the earliest climate simulation used here), which represents the end of the last interglacial warmth (MIS 5e), before glacial inception. The Eemian period (120 ka BP) differs from the pre-industrial mainly in insolation. The earlier parts of the Eemian (e.g. 125 ka BP) are often studied due to their higher temperature and sea level compared to the Holocene (Dutton and Lambeck, 2012), but 120 ka BP is the oldest point for which both FAMOUS and HadCM3 climates were available.

 

[Figure omitted. See PDF]

Reconstructed biomes (defined through highest affinity score) superimposed on simulated biomes using FAMOUS (B4F, left) and HadCM3 (B4H, right) climates for selected marine isotope stages (denoted in ka BP).

[Figure omitted. See PDF]

Pre-industrial

Our BIOME4 simulations were forced using anomalies from the pre-industrial climates produced by HadCM3 and FAMOUS. Differences between B4H and B4F for this period thus only arise from the way the pre-industrial climate forcing has been interpolated onto the two different model grids we used. Differences between B4H and B4F and the pollen-based reconstructions for this period highlight biases that are not directly derived from climates of HadCM3 and FAMOUS but are inherent to BIOME4, the pollen-based reconstruction method, or simply the limitations of the models' geographical resolution.

Although few of the long pollen records synthesized in this study extend to the modern period and their geographical coverage is sparse, a comparison with previous high-resolution biomizations of BIOME6000 (see Table 1 for details; these studies include the sites synthesized here amongst many others) and Marchant et al. (2009) show that they are generally representative of the regionally dominant biome. The biomized records of Carp Lake and Lake Tulane in North America are exceptions, showing dry grassland conditions rather than the forests (conifer and warm mixed, respectively) that are more typical of their regions (Williams et al., 2000).

There is generally very good agreement between B4H and B4F for this period and the high-resolution BIOME6000 and Marchant et al. (2009) studies. A notable exception, common to both B4H and B4F, can be seen in the southwest USA being misclassified compared to the regional biomization of Thompson and Anderson (2000). The open-conifer woodland biome they assign to sites in this region appears to be sparsely distributed (their Fig. 2) amongst larger areas likely to favour grassland and desert, and thus may be unrepresentative of areas on the scale of the climate model grid boxes. The limitations of HadCM3 and FAMOUS's spatial resolution appear most evident in South America, where the topographically influenced mix of forest and grassland biomes found by Marchant et al. (2009) cannot be correctly reproduced, with disagreement at the grid-box scale between B4F and B4H. Eurasia is generally well reproduced, although the Asian boreal forest biome does not extend far enough north, and overruns what should be a broad band of steppe around 50${}^{\circ }$ N on its southern boundary. Australia, with a strong gradient in climate from the coasts to the continental areas also shows the influence of the coarse model resolutions, with B4F more accurately reproducing the southern woodlands but neither simulation reproducing the full extent of the desert interior. Both Australian records are from the eastern coastal ranges; there are no long continuous records in the interior because of the very dry conditions. Overall, our comparison with the full BIOME6000 data set gives reasonable support to our working hypothesis that BIOME4, operating on the relatively coarse climate model grids we use here, is capable of producing a realistic reconstruction of global biomes, although local differences may occur.

6 ka BP mid-Holocene

As for the pre-industrial period, in both the mid-Holocene and LGM periods the high-resolution biomizations of the BIOME6000 project (see Table 1) provide a better base for comparison of our model results than the relatively sparse, long time-series pollen records synthesized in this study. A common thread in the BIOME 6000 studies is the global similarity between the reconstructions for 6 ka BP and the pre-industrial period, and this is, by and large, also the result seen in B4H and B4F. An increase in vegetation on the northern boundary of the central Africa vegetation band is the most notable difference compared to the pre-industrial in the regional biomizations (Jolly et al., 1998), which is also suggested by the long central African pollen records synthesized here. Both climate model-based reconstructions show grassland on the borders of pre-industrial desert areas in North Africa, although the additional amount of rainfall in both models is too low, and the model resolution insufficient to represent any significant “greening” of the desert. B4F shows a smaller change in tropical forest area in central Africa than B4H does, agreeing better with the regional biome reconstructions. Both HadCM3 and FAMOUS predict similar patterns and changes in precipitation for this period, but the magnitude of the rainfall anomaly in FAMOUS is slightly lower. The reduction in forest biomes at the tip of South Africa in B4F has some support from Jolly et al. (1998), although B4F initially overestimates forest in this area.

B4H and B4F show limited changes elsewhere too. In North America, FAMOUS's increase in rainfall anomalies produces more woodland in the west in B4F compared to the pre-industrial period, which is not seen in B4H. This is not a widespread difference shown in the regional biomization, although individual sites do change. Marchant et al. (2009) suggest drier biomes than the pre-industrial for some northern sites in Latin America, agreeing with B4F but not B4H. For Eurasia and into China, Prentice (1996), Tarasov et al. (2000), and Yu et al. (2000), all suggest greater areas of warmer forest biomes to the north and west across the whole continent, with less tundra in the north. However, neither BIOME4 simulation shows these differences, with some additional grassland at the expense of forest on the southern boundary in B4H, and B4F predicting more tundra in the north. Although both FAMOUS and HadCM3 produce warmer summers for this period, in line with the increased seasonal insolation from the obliquity of the Earth's orbit at this time, the colder winters they also predict for Eurasia skew annual average temperatures to a mild cooling which appears to prevent the additional forest growth to the north and west seen in the pollen-based reconstructions.

21 ka BP (Last Glacial Maximum)

For the LGM, both the BIOME4 simulations and pollen-data-based reconstructions predict a global increase in grasslands at the expense of forest, with more tundra in northern Eurasia and desert area in the tropics than during the Holocene. Along with the cooler, drier climate, lower levels of atmospheric CO${}_{\mathrm{2}}$ also favour larger areas of these biomes. Our long pollen records do not have sufficient spatial coverage to fully describe these differences, showing only smaller areas of forest biomes in southern Europe, central Africa, and Australia, but there is again good general agreement between our two BIOME4 simulations and the regional biomizations of the BIOME6000 project.

The FAMOUS and HadCM3 grids do not seem to have sufficient resolution to reproduce much of the band of tundra directly around the Laurentide ice sheet in either B4H or B4F, but the forest biomes the simulations show for North America are largely supported by Williams et al. (2000). However, Thompson and Anderson (2000) suggest larger areas of the open-conifer biome in the southwestern USA than in the Holocene that the BIOME4 simulations again do not show. Both B4H and B4F predict a smaller Amazon rainforest area. Marchant et al. (2009) suggest that the Holocene rainforest was preceded by cooler forest biomes, whereas both HadCM3 and FAMOUS simulate climates that favours grasslands. Marchant et al. (2009) also provide evidence for cool, dry grasslands in the south of the continent; FAMOUS follows this climatic trend but B4F suggests desert or tundra conditions, whilst B4H shows a smaller area of the desert biome. For Africa, Elenga et al. (2004) show widespread grassland areas where the Holocene has forest, with which the simulations agree, and dry woodland in the southeast, which neither B4H or B4F show; HadCM3 and FAMOUS appear to be too cold for BIOME4 to retain this biome. Elenga et al. (2000) also shows increased grassland area in southern Europe, which is not strongly indicated by either B4H or B4F, which have some degree of forest cover here.

The large areas of tundra shown by Tarasov et al. (2000) in northern Eurasia to the east of the Fennoscandian ice sheet are well reproduced by the BIOME4 simulations, although HadCM3's slightly wetter conditions produce more of the boreal forest in the centre of the continent in B4H. The generally smaller amounts of forest cover in Europe in B4F agree with the distribution of tree populations in Europe at the LGM proposed by Tzedakis et al. (2013) better than those from B4H, possibly due to HadCM3's warm bias at the glacial maximum. Both B4H and B4F agree with the smaller areas of tropical forest in China and southeast Asia reconstructed by Yu et al. (2000) and Pickett et al. (2004) compared to the Holocene but have too much forest area in China compared to the biomization of Yu et al. (2000). Neither BIOME4 simulation reproduces the reconstructed areas of xerophytic biomes in southern Australia, or the tropical forest in the north (Pickett et al., 2004).

54 ka BP (early Marine Isotope Stage 3)

There are fewer published biomization results for periods before the LGM, so our model–data comparison is restricted to the pollen-based biomization results at sites synthesized in this paper. Of these sites, only two sites show a different megabiome affiliation when compared to the LGM: in South America, Uyuni shows highest affinity scores for the forest biome, and in Australia, Caledonia Fen shows highest affinity scores for the dry woodland biome (both sites show highest affinity score for grassland during the LGM). Overall, the few sites where data are available show few differences compared with the LGM. This is perhaps a surprise given the evidence that this was a relatively warm interval within the last glacial, at least in Europe (Voelker et al., 2002). These mostly unchanged biome assignments derived from our pollen-data records are supported by our BIOME4 simulations in that, although both FAMOUS and HadCM3 do produce relatively warm anomalies compared to the LGM, both B4H and B4F simulations at 54 ka BP are similar to the LGM close to the pollen sites in the Americas, most of southern Europe (apart from Ioannina, where the data show highest affinity scores for temperate forest) and east Africa.

In other parts of the world, the biomes simulated at 54 ka BP in B4H and B4F do differ significantly from those of the LGM. Both BIOME4 simulations show increased vegetation in Europe and central Eurasia due to the climate influenced by the smaller Fennoscandian ice sheet, as well as reduced desert areas in North Africa and Australia, generally reflecting a warmer and wetter climate under higher CO${}_{\mathrm{2}}$ availability than at the LGM. However, our simulations disagree on both the climate anomalies and the likely impact on the vegetation in several areas in this period. These include differences, both local and far-field, related to prescribed ice sheets, particularly in North America, where the ice-sheet configuration in FAMOUS shows largely separate Cordilleran and Laurentide ice sheets compared to the more uniform ice coverage of the continent in HadCM3. Further afield, B4H has significantly more tropical rainforest, especially in Latin America, and predicts widespread boreal forest cover right across Eurasia. B4F, however, reproduces a more limited forest extent, with more grassland in central Eurasia. The differences in the tropics appear to be linked to larger rainfall anomalies in HadCM3 than FAMOUS, whilst the west and interior of northern Eurasia is cooler in FAMOUS than HadCM3. This may be due to the erroneously variable and low CO${}_{\mathrm{2}}$ applied to FAMOUS from the Vostok record around this period, or it may indicate a stronger response to precessional forcing in FAMOUS, with a greater influence from the Fennoscandian ice sheet.

64 ka BP (Marine Isotope Stage 4)

There are only a few differences between biomized records at the LGM and 54 and 64 ka BP. Apart from one southern European site (Ioannina), which has a highest affiliation with grassland (compared with temperate forest during the LGM), the pollen biome affiliations are much the same as at the LGM for the sites presented here. The two sites in northern Australasia show a highest affiliation with the warm-temperate forest biome during this period, compared with tropical forest at 54 ka BP; however, affinity scores between the two types are close, so this is unlikely to be related to different climates. The BIOME4 simulations support this as they also do not show major differences at the pollen sites.

Both B4H and B4F are, in general, similar for 64 and 54 ka BP. The 64ka BP climate in HadCM3 is cooler and drier than for 54 ka BP, with B4H producing larger areas of tundra in north and east Eurasia and patchy tropical forests. There is less difference between 64 and 54 ka BP in the FAMOUS reconstructions, which simulates a cooler climate at 54 ka BP compared to HadCM3, so B4F and B4H agree better in this earlier period than at 54 ka BP. North American vegetation distributions primarily differ between B4H and B4F in this period due to the different configurations of the Laurentide ice sheet imposed on the climate models.

84 ka BP (Marine Isotope Stage 5b)

The pollen-based biomization for 84 ka BP clearly reflects the warmer and wetter conditions with more CO${}_{\mathrm{2}}$ available than at 64 ka BP, especially in Europe, with the majority of sites showing highest affinity scores for the temperate forest biomes. Sites in other parts of the world show similar affinity scores to those at the 64 ka BP time slice, although there are not many sites and it is less clear whether they reflect widespread climatic conditions.

The BIOME4 simulations reflect the warmer European climate resulting from the smaller Fennoscandian ice sheet at 84 ka BP than 64 ka BP, with B4F showing some European forest cover, and B4H extending Eurasian vegetation up to the Arctic coast. However, B4H shows more of this vegetation to be grassland rather than forest, probably a result of a slightly cooler climate in HadCM3. Around the southern European pollen sites themselves, however, B4H shows little difference from the dry grassland biomes present at 64 ka BP and B4F predicts dry woodlands, perhaps a result of the models' representation of the Mediterranean storm tracks that would bring moisture inland which are often poorly reproduced in lower-resolution models (Brayshaw et al., 2010).

Although there are differences in the configuration of the Laurentide ice sheet between HadCM3 and FAMOUS, both B4H and B4F reproduce dry vegetation types in the US Midwest and significant boreal forest further north at 84 ka BP. Both BIOME4 simulations show significantly smaller desert areas in North Africa and larger areas of forest in the tropical belt than at 64 ka BP, reflecting significant precipitation and higher CO${}_{\mathrm{2}}$ levels here, although both also show a dry anomaly over Latin America. Because of increased rainfall in Australia, B4H shows a smaller desert compared with 54 ka BP.

120 ka BP (last interglacial period, Marine Isotope Stage 5e)

This time slice represents the previous interglacial, and would be expected to have the smallest anomalies from the pre-industrial control climate of the climate models. The pollen-based biomization shows widespread forest cover for Eurasia, with the only other difference from both the 84 ka BP period and the pre-industrial control being Lake Titicaca, which has the highest affinity toward desert for this period. The affinity scores for temperate forest are almost as high for this site, and neither HadCM3 nor FAMOUS has the resolution to reproduce the local climate for this altitude well (Bush et al., 2010), although both do reflect dry conditions near the coast here.

The models do indeed produce relatively small climate anomalies and vegetation similar to the pre-industrial control and each other. Both models produce widespread forest cover north of 40${}^{\circ }$ N, much as for the pre-industrial climate, although FAMOUS is slightly too wet over North America for B4F to produce Midwest grasslands as seen in B4H. Both B4H and B4F increase the extent of their tropical forests, although FAMOUS has a relative dry anomaly over central Africa, and B4F has less tropical forest than for the pre-industrial or B4H, which once again appears to have a stronger response to precessional forcing.

Global terrestrial vegetation changes

The BIOME4 simulations compare reasonably with biomizations of BIOME 6000 (pre-industrial period, 6 ka BP, and LGM) and those presented in this paper. Below we calculate changes in the biome areas and net primary productivity, keeping in mind that these calculations carry some uncertainties relating to several mismatches. As is discussed in Sect. 3.1 there are several occasions where the modern biomized pollen data do not agree with actual biome presence, for example Potato Lake and Lake Tulane in North America. In both cases high contributions of Pinus and some other taxa skewed the affinity scores towards drier biomes (grassland and dry woodland). For the past, not knowing whether a pollen distribution is representative for an area puts restrictions on the biomization method. It is, however, noted that in most cases the biomized pre-industrial pollen agree with pre-industrial biomes. The climate models produce some differences in climate forcing of the vegetation due to (1) difference in resolution, affecting the biome areal extent and altitude and (2) ice-sheet extent, affecting temperature (Sect. 3.2). We can use the pre-industrial period as a test bed to compare model outputs and pollen reconstructions (using the BIOME 6000 database), showing that there are some biases that can be attributed to biases in BIOME4, some to the biomization method, and some to the models' limiting geographical resolution.

Global area coverage of megabiome types in the model reconstructions. $S$ indicates the inclusion of potentially vegetated continental shelves after sea level lowering, and NS indicates no vegetated continental shelves following sea level lowering. FAMOUS megabiome areas are represented by a dotted line between 30 and 60 ka BP in the period where the Vostok CO${}_{\mathrm{2}}$ data used to force the simulation is thought to be erroneously low.

[Figure omitted. See PDF]

Biome areas

Whilst there is general agreement between B4H and B4F, there are also areas and periods with significant regional differences. A clearer picture of the effect on the global biosphere can be seen by using the global total areas of each megabiome for the two simulations (Fig. 3). Cooler temperatures, reduced moisture, and lower levels of CO${}_{\mathrm{2}}$ through the glacial result in a general reduction of forest biomes and increases in grassland, desert, and tundra. Lower levels of atmospheric CO${}_{\mathrm{2}}$ also preferentially favour plants using the C4 photosynthetic pathway (Ehleringer et al., 1997), contributing to the expansion of the grassland and desert biomes during the glacial. The changes in atmospheric CO${}_{\mathrm{2}}$ levels through the glacial cycle are largely common to all the BIOME4 simulations, so CO${}_{\mathrm{2}}$ fertilization effects and C3/C4 competition are generally not responsible for differences in vegetation response between B4F and B4H. The exception to this rule comes between 40 and 60 kyr BP, where the FAMOUS runs sees erroneously strong CO${}_{\mathrm{2}}$ variations in this time interval from the Vostok record which may affect both the climate used to force B4F and the fertilization effects. B4F predicts consistently lower areas of warm-temperate and boreal forest than B4H, and higher amounts of grassland and desert. FAMOUS also neglects the additional area of land that HadCM3 sees as continental shelves are exposed, reducing the area of land available to the biosphere, although some of this additional land is occupied by the Northern Hemisphere ice sheets in HadCM3. The global total areas of biomes highlights a significant oscillation in the areas of the different megabiomes of $\sim$ 20 kyr in length – this is particularly notable between 60 and 120 ka BP in the grassland megabiome and results from the 23 kyr cycle in the precession of the Earth's orbit. The precession cycle exerts a significant influence on the seasonality of the climate, as noted in tropical precipitation records (e.g. the East Asian monsoon; Wang et al., 2008; Carolin et al., 2013). Such variations are not explicitly evident in the dominant megabiome types at any of the pollen sites, but the precession oscillation does appear in the individual biome affinity scores of several sites (Supplement), lending support to this feature of the model reconstructions.

Net primary production throughout the last glacial cycle derived from the model-based biome reconstructions. B4H includes the additional influence of land exposed by sea-level changes, while B4H_NS and B4F do not. B4F NPP is represented by a dotted line between 30 and 60 ka BP in the period where the Vostok CO${}_{\mathrm{2}}$ data used to force the FAMOUS simulation is thought to be erroneously low.

[Figure omitted. See PDF]

Net primary productivity

Net primary productivity (NPP) is the net flux of carbon into green plants (in this case terrestrial plants) due to photosynthesis, after accounting for plant respiration. Global NPP derived from our BIOME4 simulations for the PI is 74 for B4H and 78 PgC yr${}^{-\mathrm{1}}$ for B4F (Fig. 4). These values are somewhat higher than previously estimated present-day ranges of 46 to 62 PgC yr${}^{-\mathrm{1}}$ (Tinker and Ineson, 1990; Nemani et al., 2003). Recent estimates using eddy covariance flux data estimate global NPP as $\sim$ 62 PgC yr${}^{-\mathrm{1}}$ (assuming 50 % carbon use efficiency to convert from GPP to NPP; Beer et al., 2010).

Some other model estimates for the PI are also lower (e.g. Prentice et al., 2011: 59.2 PgC yr${}^{-\mathrm{1}})$. As mentioned in Sect. 3.3.1, BIOME4 is driven solely by an observational climate data set for the pre-industrial due to the anomaly approach used to reduce the impact of climate model biases (see Methods Sect. 2.2.3). Therefore, any overestimate in NPP is not a result of the climate model forcing but possibly due to biases in the vegetation model, and/or biases in the observational climatology used to drive the model, and the spatial resolution used. For example, the lower-resolution topography does not represent mountainous regions such as the Andes well nor its topographically induced variation in vegetation (see Sect. 3.3.1), which may positively skew NPP values. The model may also overestimate NPP compared to observationally based techniques for the modern or pre-industrial period, partly because it does not contain any representation of non-climatically induced changes, e.g. cultivation or land degradation.

The LGM BIOME4 simulations show a global NPP decline to $\sim$ 42 in B4F and 48 PgC yr${}^{-\mathrm{1}}$ in B4H (Fig. 4). While these are also higher than some other model-based estimates of 28–40 PgC yr${}^{-\mathrm{1}}$ (e.g. François et al., 1999, 2002), the relative decrease in the LGM in our simulations to approximately two-thirds of PI is consistent with several previous studies. A calculation based primarily on isotopic evidence has produced an even lower estimate of LGM NPP of 20 $\pm$ 10 PgC yr${}^{-\mathrm{1}}$ (Ciais et al., 2012), with LGM primary productivity approximately 50 % lower than their PI estimate.

The PI–LGM difference is greater in B4F than in B4H (Fig. 4), primarily due to the fact that HadCM3's glacial land area increases as sea level lowers, enabling additional NPP on continental shelf regions, whereas FAMOUS land area remains the same. This is demonstrated by recalculating global NPP for B4H neglecting exposed shelf regions (B4H_NS), which then matches the values from B4F (Fig. 5a, green line). The effect of vegetating continental shelves on global NPP is small in comparison to the overall decrease during the glacial period; NPP reduction at the LGM is 40 % for B4H_NS and 35 % for B4H compared to the PI. The impact of large continental ice sheets reducing the land surface area available for primary production has a negligible effect on NPP compared to reduced CO${}_{\mathrm{2}}$ and glacial climate change. These high-latitude areas only contribute a small fraction of global NPP in any case, and if the area covered in ice at the LGM is excluded from NPP calculations of the PI, global NPP only decreases by a maximum of $\sim$ 5 PgC yr${}^{-\mathrm{1}}$. In addition, sensitivity tests with B4H, with and without CO${}_{\mathrm{2}}$, variation suggests that CO${}_{\mathrm{2}}$ fertilization, rather than climate, is the primary driver of lower glacial NPP in the model (accounting for around 85 % of the reduction in global NPP at the LGM).

Some differences in the timing of some multi-millennial peaks/troughs in NPP between B4H and B4F are apparent, especially in the earlier half of the simulation. These differences, all of the order of a few thousand years, can largely be ascribed to the different CO${}_{\mathrm{2}}$ forcings used for B4H and B4F as well as the multiple snapshot setup of the HadCM3 run, which only produces simulations at 2 or 4 ka intervals, compared to the 1 ka resolution of B4F. Differences in the forcing provided by the ice-sheet reconstructions used in the models, as well as in the strength of their responses to orbital forcing in the early part of the glacial (e.g the oscillations in area coverage of various biomes in Fig. 3), may also play a role.

Both BIOME4 simulations predict slightly lower NPP during the Eemian (3–5 PgC yr${}^{-\mathrm{1}}$ lower), compared with pre-industrial times. The first large-scale decrease in NPP occurs during the initial glaciation following the Eemian, between 120 and 110 ka BP (in both simulations). There is then a second large drop of $-$10 (HadCM3_S) to $-$20 PgC yr${}^{-\mathrm{1}}$ (B4H_NS, B4F) between 75 and 60 ka BP, associated with MIS 4. NPP then increases during MIS 3, followed by the final reduction ($-$10 PgC yr${}^{-\mathrm{1}}$) to lowest values during the LGM (Fig. 4). We note here that the details of the magnitude and timing of the NPP variations will be highly dependent on the prescribed CO${}_{\mathrm{2}}$ curve given that CO${}_{\mathrm{2}}$ fertilization is the predominant factor driving the changes. A recent composite CO${}_{\mathrm{2}}$ curve derived from several ice core records (Bereiter et al., 2013) has CO${}_{\mathrm{2}}$ that is 5–20 ppm higher during MIS 3 and MIS 4 than Vostok. Further sensitivity tests with B4F forced with higher CO${}_{\mathrm{2}}$ levels suggest that NPP could be up to 8 PgC yr${}^{-\mathrm{1}}$ higher at certain time slices (see Supplement Fig. S1).

Changes in NPP will likely affect terrestrial carbon storage, which in turn influences the stable carbon isotope composition (${\mathit{\delta }}^{13}$C) of seawater because terrestrial organic carbon is depleted in ${}^{13}$C. Various early modelling studies and data-based reconstructions produced a range of 270–1100 PgC decrease in terrestrial carbon storage during the LGM compared with pre-industrial time (see summary Table 1 in Köhler and Fischer, 2004). These early estimates were unreliable, however, because (a) they do not account for variation in carbon storage within biomes and (b) they neglect the substantial influence of atmospheric CO${}_{\mathrm{2}}$ concentration on carbon storage (see Prentice and Harrison, 2009, for a fuller discussion). More recent studies have narrowed the range of LGM terrestrial carbon storage decreases to 300–700 PgC. Using isotopic and modelling methods Ciais et al. (2012) suggested that only 330 PgC less carbon was stored in the terrestrial biosphere at the LGM than PI Holocene. They included a large inert carbon pool to represent permafrost and peatland carbon storage in their modelling, and their optimization procedure suggested that this inert carbon pool was larger by 700 PgC at the LGM than PI, meaning the reduction in their active terrestrial biosphere was therefore larger than most other studies have suggested, at approximately 1000 PgC.

Globally decreased LGM deep ocean stable carbon isotope ratios (${\mathit{\delta }}^{13}$C), as recorded by benthic foraminifera at $-$0.3 to $-$0.4 ‰, have also been used as an alternative method to calculate the decrease in global LGM terrestrial carbon storage compared with the PI (e.g. Broecker and Peng, 1993; Duplessy et al., 1988; Bird et al., 1996; Kaplan et al., 2002; Beerling et al., 1999). A recent estimate derived from a compilation of 133 ocean cores is $-$0.34 $\pm$ 0.13 ‰ (Ciais et al., 2012), and an ensemble of ocean circulation model simulations suggests a similar decrease of $-$0.31 $\pm$ 0.2 ‰ (Tagliabue et al., 2009). Robust reconstructions of terrestrial carbon storage could be used in a similar but inverted approach to estimate global ocean ${\mathit{\delta }}^{13}$C changes over the same time period.

From our simulations of changes in NPP over the glacial cycle we would expect lower terrestrial carbon storage shortly following the last interglacial period, with lowest values during the LGM. We would also expect, given the compensation in terms of NPP, that the vegetation on the exposed continental shelves would be an important consideration for changes in total terrestrial carbon storage. However, the large uncertainties associated with both the climate and biome models and their forcings, as well as those involved in deriving full estimates of carbon storage and ocean ${\mathit{\delta }}^{13}$C from the variables that are explicitly produced in the models, currently prohibit the robust quantitative reconstruction of these quantities from our results.

Summary

We have created a new global synthesis and biomization of long pollen records, and used it in conjunction with model simulations to analyse the sensitivity of the global terrestrial biosphere to climate change over the last glacial–interglacial cycle. Model output and biomized pollen data generally agree, showing a reduction in the global average areas of tropical forest, warm-temperate forest, and temperate forest biomes during the LGM, MIS 4, and cool substages of MIS 5, whilst showing an increase in the global average areas of the grassland and dry shrubland, desert, and tundra biomes. BIOME4 simulations of global net primary productivity also indicate significant reductions at those intervals, driven by changes in vegetated land area and CO${}_{\mathrm{2}}$ fertilization.

Existing data coverage is still low, and so there are still large areas of uncertainty in our knowledge of the palaeo-Earth system. Better spatial and temporal coverage for all parts of the globe, especially lowland areas, is required, and for this we need data from new sites incorporated into global data sets that are easily accessible by the scientific community.

The synthesized biomized data set presented in this paper can be downloaded as a supplement to this paper, or may be obtained by contacting the authors. Outputs from the climate and biome model simulations are also available from the authors.

The Supplement related to this article is available online at doi:10.5194/cp-12-51-2016-supplement.

Acknowledgements

W. A. Watts is thanked for his contribution to counting the Monticchio pollen record. The authors would like to acknowledge funding from the NERC QUEST programme (grant NE/D001803/1) and NCAS-Climate. J. S. Singarayer and P. J. Valdes acknowledge the British Broadcasting Corporation (BBC) for funding initial climate simulations. Edited by: H. Fischer