Global vegetation patterns of the past 140,000 years

Insight into global biome responses to climatic and other environmental changes is essential to address key questions about past and future impacts of such changes. By simulating global biome patterns 140 ka to present, we aimed to address important questions about biome changes during this interval.


| INTRODUC TI ON
A recent study of almost 600 palaeovegetation records concluded that climatic changes since the last glacial maximum (LGM) drove moderate-to-large changes in both the composition and structure of vegetation at most locations across the globe (Nolan et al., 2018). Vegetation modelling studies have reached similar conclusions, albeit they also emphasize the role of the increase in atmospheric carbon dioxide concentrations ([CO 2 ] atm ) since the LGM in driving observed vegetation changes (Prentice, Harrison, & Bartlein, 2011). Although such studies contrasting the LGM and recent past provide valuable insights into global vegetation changes across the past 21 kyr, both global climate and [CO 2 ] atm have fluctuated substantially, on time scales of 10 4 -10 5 years, throughout at least the past 800 kyr (EPICA community members, 2004;Jouzel et al., 2007;Luthi et al., 2008), whilst climate has also exhibited strong centennial to millennial fluctuations (Rasmussen et al., 2014). Using the LPJ-GUESS dynamic vegetation model (Smith, Prentice, & Sykes, 2001;Smith et al., 2014),  showed the sensitivity of vegetation in higher latitudes of the Northern Hemisphere to the contrasts between simulated stadial and interstadial climates during the last glacial stage (Singarayer & Valdes, 2010). Other research has shown a strong correlation in southern Africa between high present diversity of endemic bird species and areas where the same biome was simulated to have persisted throughout the past 140 kyr (Huntley et al., 2016).  Smith et al., 2001Smith et al., , 2014 to simulate global vegetation for a series of time slices, 0-140 ka BP, for which climatic conditions had been simulated. The LPJ-GUESS simulations were used to infer the extent and location of a series of biomes for each time slice, allowing global palaeovegetation maps to be drawn, then analysed to address these questions.

| Palaeoclimate experiments
An internally consistent series of 88 palaeoclimate experiments, including seven 'hosing' experiments designed to mimic Heinrich Events H0-H6, in which a 1-Sv fresh water flux was added to the surface layer in the North Atlantic for 200 years, plus a pre-industrial experiment, were performed using the HadCM3 fully coupled atmosphere-ocean general circulation model (AOGCM). The design of the experiments followed Singarayer and Valdes (2010) with respect to the boundary conditions and forcings applied, but with two important exceptions. Firstly, rather than using a static pre-industrial land surface for all experiments, the AOGCM was coupled to Triffid (Cox, 2001), a simple DGVM with five plant functional types (PFTs), enabling the vegetation and climatically important landsurface properties (e.g. forest vs. grassland) to interact dynamically with the atmosphere. The precise model configuration is HadCM3B-M2.1aD, as defined in Valdes et al. (2017). Secondly, rather than using ice sheets of LGM extent but reduced elevation for pre-LGM time slices of the last glacial stage, and modern ice sheets not only for LIG time slices but also for pre-LIG time slices, more realistic representations of pre-LGM ice sheets were used. For pre-LGM time slices of the last glacial stage, ice-sheet configuration used the ICE-5G (Peltier, 2004) ice sheet for the post-LGM time slice with the closest matching eustatic sea-level depression. For LIG and pre-LIG time slices, ice sheets were specified following de Boer, van de Wal, Lourens, Bintanja, and Reerink (2013). More details of the model set-up can be found in Davies-Barnard, Ridgwell, Singarayer, and Valdes (2017).
For input to LPJ-GUESS, anomalies were calculated for each palaeoclimate experiment with respect to the pre-industrial experiment, in each case using the climatology of the last 30 years of the experiment (Singarayer & Valdes, 2010). Anomalies were required for the twelve monthly mean temperatures (°C), monthly precipitation totals (mm) and monthly mean cloud cover (%).

| LPJ-GUESS simulations
The LPJ-GUESS DGVM simulates the dynamics of stands of vegetation comprised of a series of PFTs defined by the user. PFTs that can coexist in a given climate compete for resources, their functional traits determining their competitive balance and, hence, their relative biomass and fractional cover. The fire module used here (Smith  Thonicke, Venevsky, Sitch, & Cramer, 2001) simulates fire disturbance and its effects based on soil moisture (as a proxy for fuel moisture), fuel load and a PFT-specific fire resistance. In addition, a stand-replacing stochastic disturbance with pre-defined average return interval represents other disturbances, such as wind storms and pest attacks. As nitrogen deposition for past periods is uncertain, the model was run without the nitrogen cycle and nitrogen limitation (i.e. the C-only version of Smith et al., 2014). The vegetation of each cell of a 0.5 × 0.5° longitude × latitude grid, comprising all grid cells that had non-zero land area for at least one time slice, was simulated. Following Huntley, Allen, Barnard, Collingham, and Holliday (2013), the potential 1901-30 monthly mean temperature, monthly precipitation and monthly mean cloudiness of grid cells on shelf areas exposed by lowered glacial sea levels were estimated using thin-plate spline surfaces fitted to data for land grid cells in the CRU TS 3.0 dataset (Harris, Jones, Osborn, & Lister, 2014). Palaeoclimatic conditions for each grid cell, including those on shelf areas, were obtained following the approach of Miller et al. (2008). Anomalies derived from the palaeoclimate experiments were applied to the CRU TS 3.0 annual time series of climatic data for 1901-30, monthly mean temperature time series first being detrended to remove the warming trend seen over those three decades. This time series was then cycled through repeatedly as required. Atmospheric carbon dioxide concentration for each time slice was determined according to Vostok ice-core data (Loulergue et al., 2008;Petit et al., 1999) using the EDC3 chronology (Parrenin et al., 2007). Obliquity was specified following Laskar et al. (2004). Twenty PFTs were used (Table 1).
Eleven tree and two grass PFTs were specified following Forrest

| Biome inference
The biome for each 0.5° grid cell and during each time slice was inferred from the simulated values of LAI and carbon mass for the 20 PFTs. A total of 21 biomes was recognized (Table 2). Inference was performed using a series of rules applied sequentially (expressed in the form of a dichotomous key in Appendix 1) and implemented in a Fortran program written for the purpose (see Appendix S3). In addition to reporting the inferred biome for each grid cell, this program also reported for each time slice: the extent of each biome, both as a total area (km 2 ) and as percentage of ice-free land area; the number and percentage of grid cells with ice-free land assigned to each biome; and the total carbon mass (kg) for each biome. The global total of ice-free land area, number of grid cells with ice-free land and global total carbon mass were also calculated and reported for each time slice. For each biome, the total carbon mass of each PFT, and the percentage of the biome's total carbon mass contributed by each PFT, were also calculated and reported for each time slice.
Using maps of the simulated carbon mass and LAI for the present for the 20 PFTs and 11 PFT aggregates, as well as a map of total simulated carbon mass, and a proposed set of 19 biomes, the initial rule set was developed, guided principally by expectations as to which PFTs would characterize each biome, and applied. Whereas most biomes could be distinguished using LAI criteria for associated PFTs, a carbon mass threshold was needed to distinguish desert and tundra biomes from the remainder, the initial threshold value being assigned by inspection of the map of total simulated carbon mass.
The resulting global distribution of simulated biomes was mapped and compared visually with published maps of simulated potential pre-industrial biomes (e.g. Prentice et al., 1992Prentice et al., , 2011 and of present biomes inferred from remaining natural and semi-natural vegetation (e.g. Olson, Watts, & Allison, 1983;Olson et al., 2001). Discrepancies

| Simulated present biomes
Tropical Evergreen Forest (TrEF-see Table 2 for all biome acronyms used below) is the most spatially extensive biome at present (

| Impacts of last glacial stage millennial climatic fluctuations on global biomes
The last glacial stage is characterized by repeated alternations between relatively warm interstadials and severely cold stadials, as demonstrated by a range of Northern Hemisphere palaeoclimatic records spanning the last glacial stage and that have sufficient temporal resolution to record these centennial to millennial fluctuations (e.g. Allen et al., 1999;Allen, Watts, & Huntley, 2000). The magnitude and rate of climatic warming from stadials to succeeding interstadials are particularly clearly established for Greenland from ice-core data (Wolff, Chappellaz, Blunier, Rasmussen, & Svensson, 2010), with the most rapid and large magnitude warming episodes being associated with Heinrich Events (Andrews, 1998). Figure 3 (see also

| Holocene and last interglacial global biome patterns
Comparisons Although there are obvious differences among all six maps presented in Figure 4, three are more similar, and also more similar to that for the present (Figure 1), namely those for 7 ka (Figure 4b  The other three biome distributions have substantial differences from the present (Figure 1), including the following:

| Last and penultimate glacial maximum global biome patterns
Biome patterns for 21 ka, the time slice most frequently used to rep-

| Global biome extents and carbon mass 140 ka to present
Over  Events H0-H6 is more than sufficient to account for the increased [CO 2 ] atm recorded in ice-core data if this carbon was returned to the atmosphere, and such a source is also consistent with observed shifts in δ 13 C of atmospheric CO 2 (Bauska et al., 2018).
As stated in Section 1, our study aimed to address several ques-  was simulated for 60 ka and of the latter for H2 (24 ka); maximum extents of both were simulated for 118 ka. Thus, our results indicate that forests everywhere have experienced moderate-to-substantial changes in global extent over the last 140 kyr. In terms of their biodiversity legacy, changing extents of tropical forests may be of the greatest importance; it is noteworthy that TrRF experienced greater changes in extent (maximum-to-minimum extent ratio 3.53) than did the more diverse TrEF (maximum-to-minimum extent ratio 1.65).
Similarly, the three categories of temperate forest show markedly different maximum-to-minimum ratios, the highest being for TeMxF (6.10). The ratio of 3.95 for TeSF may have a significant biodiversity legacy, whilst that of 1.54 for TeBEF, the most globally extensive temperate forest biome, may at least in part account for its generally higher species richness than that of TeSF. The most extreme ratio of maximum-to-minimum global extent, a value of 206, is seen for the generally least extensive of the Boreal forest biomes, BSNF, its minimum extent being reached at 120 ka and its maximum at 86 ka.  (Axmanová et al., 2020), furthermore, indicates that they consumed woody as well as herbaceous material, leading to an inference that they occupied landscapes with mosaics of grasslands, shrublands and woodlands. This suggests that ShT which, despite being dominated by the BESh PFT, has a C3G carbon mass per unit area >50% higher than Tun, and a total carbon mass per unit area two orders of magnitude greater than Tun and an order of magnitude greater than St, may have provided suitable habitat for these creatures. The extent to which LIG vegetation patterns paralleled those of the Holocene was the subject of our fourth question. As discussed above and illustrated in Figure 4, whilst there are, as expected, broad similarities, there are also numerous differences both regionally and in terms of the temporal evolution of the vegetation patterns during the two interglacials. Furthermore, these differences are neither restricted to nor most striking in higher latitudes of the Northern Hemisphere, but are obvious on all vegetated continents. Once again, it is the climatic consequences of the quantitatively different forcing that underlie the differences in simulated spatial patterns and temporal evolution of the vegetation. Both the greater magnitude of the Northern Hemisphere insolation peak, and the higher [CO 2 ] atm , at the onset of the LIG than at the onset of the Holocene (Figure 7) resulted not just in globally warmer conditions (Fischer et al., 2018), but in regional differences in climatic conditions and climatic evolution that are reflected by the simulated vegetation differences.
Finally, we asked how similar were global vegetation patterns at the PGM and LGM. Figure 5 illustrates our biome simulations for the LGM (21 ka) and PGM (140 ka), showing marked differences in many regions; these differences were described above. Once again the quantitatively different forcing, principally the more extreme Northern Hemisphere summer insolation minimum at 140 ka than at 21 ka, leads to differences in simulated climatic conditions that in turn result in the differences in biome patterns. Thus, as with the Holocene and LIG, climate, and the resulting extent and distribution of biomes, differed between the LGM and PGM. is simulated to occur around the middle of the LIG (see Figure 4d, 118 ka), the greatest reduction in the extent of Des in the Sahara, Arabia and the Indian sub-continent is simulated early in the LIG (see Figure 4b, 126 ka). Testing these model-generated hypotheses awaits advances in dating and new records from these regions.

| CON CLUS IONS
Secondly, simulated biome patterns also differ between the last and penultimate glacial maxima, emphasizing the more general conclusion that, although orbital forcing has a series of characteristic Events had a larger impact upon the simulated global extent of the biome than did the overall orbitally forced glacial-interglacial cycle of climatic conditions ( Figure 6).
Fourthly, most of Earth's land surface has experienced at least one, and in many regions multiple, biome changes since 140 ka.
Primary exceptions are the core areas of TrEF, although the extent of TrEF in each of its principal centres, South America, sub-Saharan Africa and south-east Asia, is nonetheless simulated to have varied substantially, with marked reductions in extent at some times in all three regions.
Despite being in general accord with palaeovegetation evidence, one striking feature that our simulations do not show is a vegetated Sahara during either the early Holocene or LIG, thus conflicting with inferences of a 'green' Sahara at these times on the basis of geomorphic and fossil evidence (Larrasoaña, Roberts, & Rohling, 2013). Although all modelling studies show an enhanced West African monsoon during the early Holocene, they struggle to simulate the magnitude of change; this is a general challenge in palaeoclimatic modelling (Valdes, 2011). Early modelling studies suggested that without incorporation of the strong land surface to atmosphere feedbacks as the vegetation zones expand in response to enhanced monsoon precipitation, climate models are unable to generate sufficient precipitation in the Sahara to support the early Holocene 'greening' (Claussen & Gayler, 1997;de Noblet-Ducoudre, Claussen, & Prentice, 2000). However, this result was not replicated by fully coupled atmosphere-ocean-vegetation models participating in PMIP2 and PMIP3, including HadCM3 with TRIFFID (Braconnot et al., 2012;Harrison et al., 2014). Recent climate model simulations suggest that changes in dust could also be important (Pausata, Messori, & Zhang, 2016), although this result may have been an artefact of that particular model (Hopcroft & Valdes, 2019;Thompson, Skinner, Poulsen, & Zhu, 2019). A lack of sensitivity of modelled vegetation to rainfall in sub-Saharan Africa may also contribute to the problem (Hopcroft, Valdes, Harper, & Beerling, 2017). Collingham provided key support with respect to generating the climatic anomalies required to run LPJ-GUESS from the HadCM3 outputs, and also generated the estimates of 'present' climate for shelf grid cells exposed at the LGM. HadCM3 simulations were performed on the facilities of the Advanced Computing Research Centre at the University of Bristol (http://www.bris.ac.uk/acrc/). We are extremely grateful for very helpful comments from Jack Williams, and for those of an anonymous referee, on our initial submission that enabled us greatly to improve the paper.

DATA AVA I L A B I L I T Y S TAT E M E N T
The

Rule-based biome inference
The biome occupying each grid cell was inferred using the following rules, applied in the order given below where they are expressed in the form of a dichotomous key. The majority of the rules applied relate to the simulated leaf area index (LAI) of one or more plant functional types (PFTs -see Table 1). The initial rule, however, considers the simulated total carbon mass (C-Mass) for the grid cell, relative to a threshold value of 1.0 × 10 9 kg for a half-degree grid cell at the equator, area 3,087.57 km 2 , all of which is ice-free land. The threshold applied to a given grid cell is adjusted to account for its ice-free land area, thus taking into account the reduced areas of grid cells at higher latitudes, of coastal grid cells and of grid cells partially covered by ice sheets.