The climate and vegetation of Europe, northern Africa, and the Middle East during the Last Glacial Maximum (21 000 yr BP) based on pollen data

Pollen data represent one of the most widely available and spatially resolved sources of information about the past land cover and climate of the Last Glacial Maximum (LGM; 21 000 yr BP). Previous pollen data compilations for Europe, the Mediterranean, and the Middle East, however, have been limited by small numbers of sites and poor dating control. Here we present a new compilation of pollen data from the region that improves on both the number of sites (63) and the quality of the chronological control. Data were sourced from both public data archives and published (digitized) diagrams. The analysis is presented based on a standardized pollen taxonomy and sum, with maps shown for the major pollen taxa and biomes and the total arboreal pollen (AP), and on quantitative reconstructions of forest cover and of winter, summer, and annual temperatures and precipitation. The reconstructions are based on the modern analogue technique (MAT) adapted using plant functional type (PFT) scores and with a modern pollen dataset taken from the latest Eurasian Modern Pollen Database (EMPD) (∼8000 samples). A site-by-site comparison of the MAT and the inverse modelling method shows little or no significant difference between the methods for the LGM, indicating that the presence of low-CO₂ conditions and no modern analogue during the LGM does not appear to have had a major effect on MAT transfer function performance. Previous pollen-based climate reconstructions using modern pollen datasets show a much colder and drier climate for the LGM than both inverse modelling and climate model simulations do, but our new results suggest much greater agreement. Differences between our latest MAT reconstruction and those in earlier studies can largely be attributed to bias in the small modern dataset previously used and to differences in the method itself (Brewer et al., 2008; Salonen et al., 2019). We also find that quantitative forest cover reconstructions show more forest than previously suggested by biome reconstructions but less forest than suggested by simply the percentage of arboreal pollen, although uncertainties remain large. Overall, we find that LGM climatic cooling and drying were significantly greater in winter than in summer but with large site-to-site variance that emphasizes the importance of topography and other local factors in controlling the climate and vegetation of the LGM.

1 Introduction

During the Last Glacial Maximum (LGM) ∼21 000 yr BP (Mix et al., 2001), the climate, vegetation, and landscape of Europe and its surrounding regions were very different than today. Scandinavia and a large part of the British Isles were covered by a single ice sheet, with separate ice sheets covering the Alps and the Pyrenees, while many smaller and lower mountainous areas were also glaciated (Ehlers et al., 2011). As a result of this global build-up of ice on land, sea levels were around 120 m lower than today, resulting in the retreat of Atlantic and Mediterranean coastlines and the emergence on land of the English Channel and the North Sea basin. Falling sea levels also led to the disconnection of the Black Sea from the Mediterranean and to a subsequent drop in Black Sea water levels as evaporation exceeded inflow (Arslanov et al., 2007). On land, permafrost and periglacial processes occurred immediately to the south of the Scandinavian ice sheet, while the massive discharge of glacial clays and sands provided material to be redeposited by the wind as belts of loess across northern France, Benelux, Germany, and central Europe (Lehmkuhl et al., 2021). Under these cooler and drier climatic conditions, forests are thought to have retreated to the relative shelter of southern Europe and the Mediterranean, while relatively unproductive steppe and tundra dominated the region north of the Alps (Grichuk, 1992).

This traditional view of the LGM has been established for many years, but many details concerning the climate and vegetation of the LGM are still debated. Much of this debate concerns information derived from the pollen record, which represents one of the most widely available and spatially resolved sources of information concerning LGM vegetation and climate, and the primary terrestrial proxy used to evaluate climate models in the Paleoclimate Modelling Intercomparison Project (PMIP) (Bartlein et al., 2011; Braconnot et al., 2012; Harrison et al., 2014).

For example, climate model simulations continue to indicate a climate that is less cold and more humid than pollen-based reconstructions (Jost et al., 2005). These results are similar to reconstructions based on glaciological modelling (Allen et al., 2008b). On the other hand, the pollen-based reconstructions that show the greatest disagreement with climate models have themselves been criticized for not considering the possible effect of low atmospheric CO₂ on the physiological relationship between plants and climate (Ramstein et al., 2007). Methods that use modern pollen samples are based on the assumption that the relationship between vegetation and climate remains the same through time and that this is independent of change in CO₂ concentration. Studies have shown, however, that plant growth processes and plant resilience are sensitive to CO₂ concentration and particularly to water-use efficiency which would make plants more drought-sensitive in low-CO₂ environments (Cowling and Sykes, 1999). Atmospheric CO₂ during the LGM was around 190 ppm, some 100 ppm lower than the pre-industrial period and 200 ppm lower than the levels experienced in the last 50 years. Concerns about the effects of lower CO₂ during the LGM have directly led to the development of pollen–climate reconstruction methods that can take account of CO₂ effects, either through the use of a process-based vegetation model run in inverse mode (Guiot et al., 2000) or through the use of a correction algorithm (Prentice et al., 2017; Cleator et al., 2020). Pollen–climate reconstructions based on inverse modelling that account for these low-CO₂ effects show less cooling and drying and consequently greater agreement with climate models (Ramstein et al., 2007; Wu et al., 2007, 2019; Izumi and Bartlein, 2016).

Further data–model discrepancies have also been highlighted concerning LGM vegetation cover. Earlier pollen synthesis studies, especially those that applied the biomization method (Elenga et al., 2000), give the impression that non-glaciated areas of LGM Europe were dominated by treeless steppe, while vegetation models driven by climate model simulations indicate large areas of forest and woodlands (Binney et al., 2017; Kaplan et al., 2016; Velasquez et al., 2021). The apparent data–model discrepancy associated with steppe has led to the suggestion that early humans, who are not included in vegetation models, could have reduced the forest cover with only a relatively moderate use of fire because of the cold climate and slow speed of vegetation recovery (Kaplan et al., 2016). This debate is important because of studies that have shown the sensitivity of the climate system to vegetation boundary conditions during the LGM (Ludwig et al., 2017; Velasquez et al., 2021). This suggests that accurate knowledge of the vegetation cover during the LGM is a necessary prerequisite to understanding the role of other influences on the climate system at this time.

More recent pollen and macrofossil studies from eastern central Europe have shown that, at least in this region, there existed areas of open boreal forest and woodland with some temperate broadleaf species (Kuneš et al., 2008; Willis and Van Andel, 2004). The evidence of forest, and particularly elements of temperate broadleaf forest, north of the Alps has come to represent a challenge to the traditional view that forest species only survived the LGM in sheltered refugia far to the south of the Fennoscandian ice sheet and close to the moderating influence of the Mediterranean Sea. The presence of micro-refugia north of the Alps is important because it would represent a very different baseline for understanding the later rate and route of plant migrations under the rapid warming that occurred during the Late Glacial to Holocene transition (Douda et al., 2014; Giesecke, 2016; Krebs et al., 2019; Nolan et al., 2018) and for understanding patterns of present-day genetic diversity (Normand et al., 2011; Svenning et al., 2008). Modelling studies have shown difficulty in supporting the very high rates of postglacial expansion that would be necessary for southern refugia (Feurdean et al., 2013; Nogués-Bravo et al., 2018).

Much of this debate has been informed by an increasing number of LGM pollen studies from an ever-broader geographical area and especially by an increasing number of studies north of the Alps. Nevertheless, the synthesis of these studies into a single narrative is made difficult by several factors, for instance, different taxonomic definitions, pollen percentages calculated from non-standardized pollen sums, and quantitative analyses such as climate reconstructions that are based on different training sets and methodologies. This has led to some modelling studies ignoring the pollen record completely, on the basis that data from the LGM are too scarce (Janská et al., 2017). Where standardized methods have been applied to multiple LGM pollen records, poor dating control has resulted in the inclusion of many records that may not actually be from the selected LGM time window. This is particularly important because the 21±2.0 ka time slice commonly used to represent the LGM period in PMIP data–model comparisons and other synthesis studies (MARGO members, 2009; Bartlein et al., 2011) occurs immediately after the glacial maxima in the Alps around 26–23 ka (Heiri et al., 2014; Spötl et al., 2021) and Heinrich Stadial 2 (HS2; 24.3–26.5) whilst also being closely followed by Heinrich Stadial 1 (HS1; 15.6–18.0 ka) (Sanchez-Goñi and Harrison, 2010). These closely associated time periods can therefore be expected to represent both a different vegetation and climate than the LGM itself.

For example, of the 18 European pollen records used in the PMIP benchmarking dataset (Bartlein et al., 2011), 10 fall into the worst class (“poor”) in the COHMAP chronological quality classification scheme if relative dating such as pollen correlation is excluded. More recent synthesis studies have also relied heavily on records from the European Pollen Database (EPD), which currently has 116 records with samples of LGM age (as of June 2022). Many of these records, however, are based on chronologies that are regarded as reliable for the Holocene (Giesecke et al., 2014) but have large uncertainties for the LGM as a result of (1) excessive extrapolation back in time from Holocene-age dates, (2) the use of pollen correlation or other relative dating despite poorly defined regional biostratigraphy, or (3) the inappropriate use of radiocarbon dates contaminated with old carbon. We found that 104 of these 116 EPD records (Williams et al., 2018) fall into the worst class (“poor”) in the COHMAP chronological quality classification.

Here we address these problems using a new synthesis of LGM pollen records from throughout Europe, the Mediterranean, and the Middle East (EurMedMidEst) based on rigorous quality-control criteria. Records were compiled from an extensive review of public databases and archives and of the scientific literature. Pollen records were selected according to the robustness of their chronological control around the PMIP LGM time window (21±2 ka) and combined into a single dataset based on a harmonized taxonomy and standardized pollen sum. The dataset was then analysed so that standardized maps could be produced to show the distribution of the major pollen taxa and biomes and the total arboreal pollen (AP) at the LGM. In addition, quantitative reconstructions of forest cover and of winter, summer, and annual temperatures and precipitation were undertaken using a modified version of the standard modern analogue technique (MAT) (Guiot et al., 1989) utilizing the latest Eurasian Modern Pollen Database (EMPD) v2 dataset. These climate reconstructions are compared and evaluated against previous LGM pollen–climate reconstructions and reconstructions based on other proxies. The dataset and results are fully documented, and the complete data files are provided in the Supplement.

2 Methods

2.1 Pollen data

LGM fossil pollen data from Europe and bordering regions including northern Africa and the Middle East were selected and collated into a single standardized project database. These data were sourced from the EPD/Neotoma database (Williams et al., 2018), the Pangaea data archive, publications in scientific journals, and the original authors. We selected LGM pollen sites and data according to strict quality-control criteria. Where possible, primary raw pollen counts were used when these were available. When the original electronic data were not available, the data were digitized from the published diagram. Overall we have included 63 records in our study, of which 35 were digitized and 28 consisted of the original pollen counts (Table 1).

Figure 1 Site locations and archives (site numbers are as shown in Table 1).

The distribution of the 63 sites reflects the distribution of suitable archives, with fewer records available from climatically or environmentally challenging regions (Fig. 1). High rates of erosion and a drier and colder climate during the LGM reduced the number of suitable anoxic sediment sinks for pollen preservation, especially in central Europe between the Scandinavian and Alpine ice sheets. Nevertheless, our dataset includes sites from this region and sites in northern Africa and eastern central Europe through to Iran, although most sites are located in an arc across eastern Spain, the Alps, and Italy. Lake sites comprise the most numerous archive and tend to be located in the more sheltered and topographically favourable regions of southern Europe and the Mediterranean. Peat is the next most important archive, followed by alluvial and colluvial sediments and by cave sites, the latter also often being known for their archaeological significance. Sites located at the ice margins that appear to be under the ice reflect uncertainties in the location of the ice margin both in time and space during the LGM and reflect the fact that the selected time window for this study (21±2 ka) is later than the maximum ice advance in some regions (Hughes and Gibbard, 2015). For completeness, we also include seven marine records which have the advantage of more continuous deposition and often better dating over the LGM period but which are prone to taphonomic biases compared to terrestrial records. These biases are discussed later in this section.

LGM pollen records were selected according to a number of quality-control criteria, but primary amongst these was the existence of sufficient independent chronological control points to accurately identify samples that would fall within the 21±2 ka BP time slice of interest. We have used all of the samples within this time frame where the samples have been available in electronic form, or else we have used the sample closest to the target time (21 ka BP). For records taken from the EPD, we have used the latest Bayesian age–depth models where these were available (Giesecke et al., 2014), otherwise we have used the dates and chronology proposed by the original authors. We classified chronologies according to the COHMAP chronological quality scheme for the LGM period (Anderson et al., 1988; Yu and Harrison, 1995), which classifies record quality from 1–6 depending on whether a date falls within 2000 ¹⁴C years (or less) of the time being assessed or whether bracketing dates fall within 6000 and 8000 ¹⁴C years (or less) of the time being assessed (Fig. A1). Chronologies based on dates that fall outside of these limits fall into COHMAP class 7 and are regarded as “poorly dated” with respect to the LGM. Importantly, we have only included radiometric and other absolute dates (such as varves) in this assessment and have excluded dates based on correlation with regional pollen records. These pollen-based stratigraphic dates have been widely used in previous LGM studies but do not include estimates of uncertainty and are generally regarded as unreliable at this time given the sparsity of well-dated pollen sites and samples on which to base any correlation (Giesecke et al., 2014).

All records that were classified as poorly dated (COHMAP class 7) were subsequently excluded from our analysis. This means that many of the pollen records used in previous studies were excluded, including 16 of the 26 LGM records used by PMIP and associated studies in Europe (Bartlein et al., 2011; Elenga et al., 2000; Tarasov et al., 2000; Jost et al., 2005; Peyron et al., 1998; Wu et al., 2007; Cleator et al., 2020). We also excluded 104 of the 116 records in the EPD with samples that fall within our LGM time window. Many of these EPD pollen records have been used in more recent studies, although the exact record (EPD entity number) is often not stated. We estimate that we have excluded 16 of the 17 European sites used by Binney et al. (2017) (this study only included sites above latitude 40° N), 5 of the 6 European sites used by Allen et al. (2010), 28 of the 33 sites used by Cao et al. (2019), and 27 of the 71 sites used by Kaplan et al. (2016).

Table 1 List of selected sites.

Other quality-control criteria were also used in the selection of LGM pollen records. Published pollen diagrams that only included a small part of the terrestrial pollen assemblage or only presented summary taxa were excluded. Records were also excluded where the dating information was incomplete, for instance, where radiocarbon dating uncertainties were not published or where it was not possible to determine if the date shown was in calibrated or uncalibrated radiocarbon years.

The modern pollen data for the climate and tree-cover reconstructions were sourced from the latest version 2 of the Eurasian Modern Pollen Database (Davis et al., 2020), which is managed as part of the EPD. EMPD2 includes 8133 modern pollen samples from across the Palearctic biogeographic region from Europe to the far east of Asia. The taxa from both the fossil and modern pollen data were consolidated into 120 of the most commonly occurring terrestrial taxa types. This taxa list was designed to be compatible with the biomization scheme used in our study (Peyron et al., 1998; Tarasov et al., 2000) and that used in the Holocene mapping study of Brewer et al. (2017). The count of Larix was amplified by a factor of 10 due to its low pollen representation (Edwards et al., 2000; Bigelow et al., 2013; Tarasov et al., 1998, 2000, 2013; Binney et al., 2017).

2.2 Biomization

We converted pollen assemblages to biomes based on the European biomization scheme of Peyron et al. (1998), which in turn is based on Prentice et al. (1996). The method is described in detail in Collins et al. (2012). We expanded the number of taxa included in the biomization procedure proposed by Peyron et al. (1998) to include taxa from the northern Eurasian biomization procedure of Tarasov et al. (1998). The inclusion of additional northern Eurasian taxa reflects recent evidence that modern analogues of LGM vegetation occur in parts of Siberia (Magyari et al., 2014a). The biomization procedure (Prentice et al., 1996) assigns each taxon to a plant functional type (PFT) and calculates a score for each of these PFTs based on the sum of the square root of the percentage of each of the taxa included in that PFT. To reduce the influence of long-distance transport, taxa below 0.5 % are removed at the start of the procedure. Each biome is then assigned one or more PFTs, and a score for each biome is calculated as the sum of the associated PFT scores. The biome with the highest score is then viewed as the dominant biome. Where the highest score is the same for more than one biome, the dominant biome is decided based on a hierarchy of unique PFTs. Peyron et al. (1998) also included a procedure for distinguishing warm and cold steppe biomes based on re-assigning certain steppe PFTs according to the presence or lack of PFTs indicative of cold or warm conditions. Following the BIOME 6000 project (Elenga et al., 2000) and Allen et al. (2010), we did not apply this additional procedure, and we present only the merged steppe biome. In summary, the biomization procedure categorized 39 arboreal pollen taxa and 39 non-arboreal taxa into 22 plant functional types (PFTs), which were then combined into 12 biomes.

2.3 Quantitative climate reconstruction

We reconstructed climate from pollen data based on a modified modern analogue technique (MAT) (Guiot et al., 1989) that used PFT scores to match fossil samples with modern pollen samples (as used by Davis et al., 2003). Other methods using PFT scores and artificial neural network techniques have been developed to reconstruct the climate of Europe during the LGM from pollen data (Peyron et al., 1998; Jost et al., 2005). PFT scores have been used in previous large-scale European pollen-based climate reconstructions for the Holocene (Davis et al., 2003; Mauri et al., 2014, 2015), where performance was found to be better than the conventional approach based on individual taxa (e.g. Marsicek et al., 2018). A particular advantage of the PFT approach for the LGM is that it can help overcome problems associated with vegetation (pollen) assemblages that may have no modern analogue (Davis et al., 2003). This can be a problem during the LGM, when the climate and environment could be expected to be very different from today and when many taxa formed unusual vegetation assemblages as a result of their forced retreat to sheltered refugia locations. The problem of modern analogues is also addressed in our reconstruction by using the latest EMPD2 modern pollen dataset. EMPD2 provides a large number of potential modern analogues for many different LGM vegetation types and climates found today across the Palearctic region. PFT scores were calculated according to the methods already outlined in the Biomization section then normalized so that each sample was proportional to every other sample (Juggins and Birks, 2012).

The MAT method was applied using the rioja package for R (Juggins, 2020). The modern pollen data were taken from the latest version 2 of the EMPD (as detailed earlier). EMPD2 includes 8133 samples, which is considerably larger than the modern datasets used in previous LGM pollen-based reconstructions. For instance, Peyron et al. (1998) used a modern pollen dataset of 683 samples, which was updated by Jost et al. (2005) to include an additional 185 samples. These datasets were also mainly taken from the steppes of Kazakhstan and Mongolia, while EMPD2 covers a much wider area, spanning most of the Eurasian Palearctic region (Davis et al., 2020). The size and distribution of the modern training set in climate and vegetation space is important because, in order for the method to work effectively, it is necessary to have samples representative of the likely vegetation and climate space that could be occupied by the fossil assemblage (Turner et al., 2021; Chevalier et al., 2020; Salonen et al., 2012; Juggins, 2013).

A known problem with the MAT is the role of spatial auto-correlation in providing unrealistically low estimates of uncertainty (Chevalier et al., 2020; Telford and Birks, 2009). This results from the fact that closely analogous modern pollen samples can also be located closely in physical space and therefore in climate space. To reduce this problem, it is possible to exclude closely located samples from the analogue-matching process using a filter based on a set distance (h-block filter) (Telford and Birks, 2009). While this approach can help, there are also three main problems associated with it. The first is error substitution, since removing samples also reduces the number of potential analogues, creating a different source of error that is not easy to categorize. Secondly, multiple samples taken from the same location are actually a strength of pollen training sets, since they are more likely to capture the full range of the assemblage diversity associated with a given climate. Thirdly, current methods that limit spatial range, such as the h-block filter, only do so on the horizontal axis and do not consider the fact that samples can also be found at different elevations. In hilly or mountainous regions, samples can therefore be excluded because they are closely located in horizontal space, but in fact they actually occupy very different climates and vegetation associations, contradicting the logical premise of the h-block filter. It was therefore decided not to apply this filter.

Uncertainties for the pollen–climate reconstructions were calculated using a standard method for the MAT (Juggins, 2020) based on the spread of the climates associated with the best modern pollen analogues used for each fossil sample. The closer the climates of the best modern pollen analogues (six in the case of this study), the smaller the calculated uncertainties assigned to the reconstructed climate of the fossil pollen sample.

Climate reconstructions are presented as anomalies. These have been calculated with respect to modern climate (1970–2000 average) at each core site location using WorldClim 2 (Fick and Hijmans, 2017) (Table A1), which was also used to assign the modern climate for the modern pollen samples in the transfer function (Davis et al., 2020).

2.4 Marine pollen records

We have included marine pollen records in our analysis for reasons explained below, but it is important that these records be viewed with caution, particularly when used for biome and quantitative MAT reconstructions and when compared with terrestrial records from different archives. Biomization methods have been applied to individual marine pollen records (Combourieu Nebout et al., 2009) and to multi-site synthesis studies such as the ACER project (ACER project members et al., 2017). However, marine records were specifically excluded from the BIOME 6000 project (Elenga et al., 2000). Similarly, quantitative climate methods have been applied to individual marine pollen records (Combourieu Nebout et al., 2009; Fletcher et al., 2010) and to multi-site synthesis studies (Sánchez Goñi et al., 2005; Brewer et al., 2008; Salonen et al., 2021). However, marine records have also been specifically excluded from other major pollen–climate studies (Cheddadi et al., 1996; Davis et al., 2003; Marsicek et al., 2018) and from quantitative forest cover reconstructions (Zanon et al., 2018).

Discussions on the advantages and problems associated with marine records can be found elsewhere (Chevalier et al., 2020; Daniau et al., 2019) but are reviewed briefly here where relevant to the methodologies applied in this study. Marine sedimentary records provide continuous and well-dated pollen records for the LGM that are often lacking from many terrestrial regions, especially in arid areas with few alternative anaerobic sediment sinks. Conversely, pollen source areas for marine sites may be many hundreds of kilometres from the coring site and may be liable to change through time in response to changes in distance to the coastline, rates of river discharge, and ocean and atmospheric dynamics. This can theoretically give rise to changes in the vegetation shown in the pollen assemblage recorded at the marine site without any actual change in climate or other environmental pressure. The large and indeterminable source area of marine records also means that it is difficult to apply quantitative MAT reconstruction methods, not least because the mean climate or forest cover of the source area is almost impossible to determine. In addition, the fossil pollen record and the modern pollen dataset to which it is being compared are composed largely of terrestrial lakes and bog sites with much smaller and more homogeneous source areas. This creates a series of problems, the more obvious of which is the calculation of anomalies, since we cannot assume that the modern climate at the (marine) coring site location is representative of the (terrestrial) source area. In this study we have taken the closest point on land as the modern climate for the calculation of anomalies, but we provide the absolute values for all sites so that these can be recalculated if necessary (Table A1). The next problem is that the large source area may capture a combination of different vegetation types that is not going to be represented in a modern pollen dataset based on samples from terrestrial sites with much smaller source areas, for instance, a mixture of coastal and mountain vegetation or even vegetation from different continents (Magri and Parra, 2002). However, in our analysis we did not find any sample from a marine record (or terrestrial record) that did not have a reasonable modern analogue in our training set (chord distance < 0.3) (Huntley, 1990), even though we did not adjust the pollen assemblage for the over-representation of Pinus (and other Pinaceae) in the marine pollen samples.

Typically, the Pinaceae component is excluded from the terrestrial pollen sum when calculating percentages for marine pollen samples, and in some cases it has been excluded entirely from the samples used in marine pollen–climate reconstructions (Combourieu Nebout et al., 2009). The problem with excluding Pinus is twofold: Pinus often represents the main forest-forming tree in the Köppen Csb climate zone on the Atlantic coast where many marine sites are located (García-Amorena et al., 2007), and it also represents the most abundant tree taxa in Europe during the LGM (Fig. A11).

The effect of excluding Pinaceae on the biomization algorithm and MAT climate reconstruction process has not been widely investigated. We therefore decided to evaluate this problem for (1) biomization and (2) pollen–climate reconstruction. In Table A2 we show the biomization results for 8213 modern pollen samples taken from the EMPD2 modern pollen database. Using this as the control, we then artificially varied the number of Pinaceae (Pinus, Abies, and Picea) in the assemblage of each pollen sample and compiled the results (Table A2). This shows quite clearly that removing all of the Pinaceae has a much more profound effect on the biomization process than artificially inflating the number of Pinaceae does (as might be expected in a marine sample where Pinaceae can be over-represented). Even when the number of Pinaceae was artificially inflated by as much as 400 % of the original value, the biomes were changed in only 2348 samples, compared to 5860 samples if all the Pinaceae were removed entirely. In terms of the effects on individual biomes, removing the Pinaceae considerably increased the amount of cold deciduous forest (CLDE), steppe (STEP), and tundra (TUND) while greatly reducing the amount of xerophytic woodlands (XERO), almost eliminating the amount of taiga (TAIG), and completely eliminating the cool conifer forest (COCO) biome. In contrast, the effect of inflating the number of Pinaceae tended to be more evenly distributed between the biomes, with the biggest increase seen in TUND and the biggest decrease seen in STEP. This suggests that, even if the over-representation of Pinaceae was quite extreme in marine pollen samples, the effect on biome classification (and, by definition, the underlying PFT scores) is less than removing Pinaceae completely from the pollen assemblage.

In a second test, we compared the reconstruction of LGM climate from marine pollen samples when Pinaceae were included and excluded. The results are shown in Table A3 and indicate that reconstructed temperatures are generally 1–2 °C cooler and that precipitation is slightly higher when Pinaceae are excluded. The differences between the two methods, however, are small and generally less than half of the uncertainties, suggesting that differences are statistically indistinguishable when considered in the context of the overall uncertainties.

In summary, we find that including Pinaceae in the biomization process is less likely to lead to misassignment of the biome than excluding Pinaceae is, except in extreme cases of over-representation. Percentages of Pinaceae in the LGM marine samples range on average between 23 %–88 %, suggesting that, while the percentage of Pinaceae was high at some sites, it does not appear to completely overwhelm the assemblage as might be expected if over-representation were to be a significant problem. We also find that including Pinaceae in the pollen assemblage of the LGM marine pollen samples gives pollen–climate reconstructions that are statistically indistinguishable from those obtained by excluding Pinaceae from the assemblage. Including Pinaceae in marine samples also provides compatibility with terrestrial samples, particularly when calculating and plotting pollen taxa percentages. For these reasons we have included Pinaceae in the analysis of all marine pollen samples in this study, although it is important to recognize that Pinaceae in such samples can be subject to over-representation and that the results presented here from marine sites should consequently be viewed with caution.

2.5 Quantitative tree-cover reconstruction

It has long been recognized that the proportional representation of individual pollen taxa in a pollen assemblage does not necessarily reflect the proportion of land area covered by those taxa in the pollen source area surrounding the sample site (Davis, 1963; Gaillard et al., 2010; Zanon et al., 2018). These differences can be caused by variations in pollen productivity, differential transport, deposition and preservation of pollen grains, and even the ease or difficulty in identifying pollen grains themselves. This can make the interpretation of pollen taxa percentages difficult, even for relatively simple questions such as the proportion of forest to non-forest in the landscape.

There have been two main methods developed to account for this quantification problem, with one using a physical modelling technique (PMT) based on estimates of pollen production for individual taxa (Gaillard et al., 2010) and the other using a MAT very similar to that used in pollen–climate reconstructions (Williams and Jackson, 2003). Both approaches have been widely applied during the Holocene in Europe (Zanon et al., 2018), but we know of no previous study that has applied either of these approaches to the LGM. The LGM presents a number of challenges, not least the problem of potential missing vegetation analogues and of low atmospheric CO₂, which has been shown to influence pollen productivity (Leroy and Arpe, 2007).

Here we use the MAT to provide quantitative estimates of forest cover, following the approach of Zanon et al. (2018), who applied this method to the Holocene pollen record of Europe. We apply the MAT in exactly the same way as for the climate reconstructions described earlier, including the use of PFT scores to match fossil and modern pollen samples. Instead of modern climate values, we assigned an estimate of modern forest cover to each of our modern pollen sites. To do this, we used a high-resolution (∼100 m) remote-sensing dataset derived from satellite observations (Hansen et al., 2013). Zanon et al. (2018) showed that the MAT calibrated in this way gives comparable results to the PMT approach in Europe, at least for the Holocene. One of the main differences, however, is that the PMT is designed to provide estimates of the proportions of different taxa, whereas the MAT (as applied here) is designed to provide estimates of the proportion of forest cover. Where the PMT can only reconstruct the proportion of forest-forming trees, irrespective of their size, the MAT (following Zanon et al., 2018) is calibrated specifically to reconstruct forest composed of trees over 5 m tall. This follows the FAO definition of forest as “land spanning more than 0.5 ha with trees higher than 5 m and a canopy cover of more than 10 %, or trees able to reach these thresholds in situ” (FAO Terms and definitions 2020, http://www.fao.org/3/I8661EN/i8661en.pdf, last access: 24 July 2024).

2.6 Maps

We present our results in the form of maps that include the main physiographic features of the LGM in the study area. The maps are based on the WGS84 projection. Coastlines reflect LGM sea level at 120 m below present, while ice sheets are based on Ehlers et al. (2011). Modern national country boundaries are also included for reference.

3 Results

3.1 Vegetation and biomes

Results of the biomization analysis show that steppe (STEP) was the most common biome at the LGM across the study area, occurring at 36 out of 63 sites, indicating that the landscape was largely dominated by cool temperate grasslands across much of western central Europe, the central and eastern Mediterranean, and northern Africa and the Middle East (Fig. 2). However, at the same time, we also find that there were a significant number of sites where woody and forest biomes occurred, more particularly in southern and eastern Iberia, northern Italy, and central eastern Europe. The most dominant of these forest and woody biomes are taiga (TAIG) in the north and cool mixed forest (COMX) and xerophytic woodlands (XERO) in the south.

Figure 2 Pollen biomes.

As would be expected, the dominance of STEP biomes is generally reflected in low arboreal pollen percentages across the same areas/sites (Figs. 3 and 4). Exceptions to this rule can be found at marine sites such as MD99-2331 site no. 3 and MD01-2430 site no. 58, where STEP is reconstructed despite arboreal pollen percentages of 71 % and 80 % respectively. This apparent contradiction illustrates some of the idiosyncrasies of the biomization method, especially when applying the method to marine pollen samples. In this case it is important to remember that the AP percentage is calculated from the sum of the percentages of each relevant taxon, but the score for each biome is the sum of the square root of the percentages of each of its constituent taxa. This results in biomes with taxa with large percentage values scoring proportionally lower and biomes with taxa with small percentage values scoring proportionally higher. For example, a single taxon at 50 % has a square root of 7.07, but the sum of the square roots of 10 taxa each at 5 % is 22.36 even though the sum of the percentages is the same 50 %. This effect can be particularly pronounced in marine pollen samples because they are usually dominated by a single taxon (Pinus) that forms a high percentage of the total assemblage. Since there are often more non-arboreal taxa than arboreal taxa in a pollen assemblage, the non-arboreal taxa can dominate in the biomization process even if collectively their percentage of the assemblage is a lot less than the arboreal taxa, resulting in a non-arboreal biome such as STEP having the highest biome score.

Figure 3 Percentage of arboreal pollen (AP) forest cover.

Figure 4 Percentage of modern analogue technique (MAT) forest cover.

Of the main arboreal biomes, taiga (TAIG) is the dominant biome at eight sites at the eastern end of the Alpine ice sheet and at a site just to the north in northern Germany and a site in Slovakia, while cool conifer forest (COCO) is found at one site close to the Scandinavian ice sheet in Lithuania. Cool mixed forest (COMX) is found much more widely at eight sites south of the Alps from southwestern Iberia to Romania, with xerophytic woodlands (XERO) occurring at eight sites with a similar distribution but not as far east or west. Cold mixed forest (CLMX) occurs at just two sites in Georgia and in the Alboran Sea at the far east and west of the study area, while warm mixed forest (WAMX) is the dominant biome at just one site in southern Spain. We do not record temperate deciduous forest (TEDE), tundra (TUND), or desert (DESE) as the dominant biome at any site at the LGM, although they do occur as sub-dominant biomes.

An alternative picture of LGM tree cover is provided by the MAT reconstructions (Fig. 4). MAT performance statistics for tree cover are shown in Table 2, based on an evaluation using the modern training set. This shows a relatively large root-mean-square error (RMSE) of 21.03 and an R² of 0.52 that is not as good as for the MAT climate analysis, but overall the results are comparable with previous MAT tree-cover studies (Zanon et al., 2018). In general, the MAT values (site average 34 %) show forest cover around 16 % less than that suggested from the AP percentage (site average 50 %; Fig. A2), although sites with very low AP percentage also show higher values based on the MAT. These differences are consistent with comparisons between the MAT and AP percentages in Zanon et al. (2018), although it should be noted that uncertainties related to the MAT reconstructions are large (±23 %). Zanon et al. (2018) found that the differences between the MAT and AP percentages were greatest over northern Europe and in Arctic and sub-Arctic climate regions that are likely to be comparable to many areas of Europe during the LGM. Today, these regions are associated with tree-forming taxa such as birch that fail to grow to a height of 5 m or more, developing only as shrubs or krummholz forms.

Table 2 MAT performance statistics based on the modern pollen sample training set. This includes mean annual temperature and precipitation (TANN and PANN), mean winter temperature and precipitation (TDJF and PDJF), and mean summer temperature and precipitation (TJJA and PJJA).

Pollen taxa percentages are shown in Fig. A3, and distribution maps of the 33 most common taxa are shown in Figs. A4–A9. Of the 21 arboreal taxa, Pinus generally has the highest values and is the most widespread, being present at all 63 sites. Other acicular arboreal taxa include Juniperus, which also has a wide distribution across EurMedMidEst, although at lower values. The rest of the acicular arboreal taxa have more regional distributions. Picea is found mainly to the north of the study region, away from the Mediterranean, whilst Abies is generally found more to the south. Larix occurs only in the central European area, including the northern edge of the Po plain just south of the Alps, whilst Cedrus is found mainly across southern and western Europe in locations much further north than its Holocene and modern distributions, which are confined mainly to Morocco and Lebanon (Collins et al., 2012). Temperate broadleaf arboreal taxa, which also include cold-tolerant species such as Betula and Salix, are relatively widely spread across the EurMedMidEst during the LGM, while taxa that are less drought-tolerant, such as Alnus, Carpinus, and Corylus, are found more to the southwest through to the northeast. Other temperate broadleaf arboreal taxa such as Quercus (deciduous) and Ulmus have a much more southern distribution, with Fraxinus, Olea, and Quercus (evergreen) being more prevalent in the southwest. In contrast, Fagus occurs more in the centre and in the east, while Tilia is found even in more northern locations of central Europe. The remaining arboreal taxa are more shrubby and drought-adapted, with Ephedra and particularly Ephedra fragilis having a southern distribution, whilst the more cold-adapted Hippophae are found even in the north of central Europe (similar to Tilia).

The main non-arboreal taxa generally indicate cool, dry, and environmentally disturbed conditions across much of the EurMedMidEst. The most widely distributed taxon is Poaceae, which, like Pinus, is found in all records. Other non-arboreal taxa with a widespread distribution include Rubiaceae, Apiaceae, and Asteraceae (Asteroideae), while Plantago, Caryophyllaceae, Brassicaceae, and Asteraceae (Cichorioideae) have a more southern and western distribution. Thalictrum can be found mostly at sites in the centre of the EurMedMidEst along with Helianthemum, which also extends to sites in the southwest. Other taxa such as Chenopodiaceae and Artemisia have a more southern distribution, reflecting their preference for climates that are drier and less cold.

3.2 Climate reconstruction evaluation

Evaluation of transfer function performance based on the modern training set is presented in Table 2. This shows that root-mean-square error of prediction (RMSEP) values were smallest for summer temperatures (2.21 °C) and largest for winter temperatures (3.35 °C), with mean annual temperatures in between (2.28 °C). The weaker performance for winter temperatures largely reflects the much greater range of winter temperatures in the training set. In turn, this contributes to a better R² performance for winter temperatures (0.91) than annual temperatures (0.9) and summer temperatures (0.81). Overall R² performance for precipitation is weaker than for temperature, which is typical because of the higher spatial variability in precipitation compared to temperature. Summer precipitation has the strongest R² performance (0.75) compared to winter and annual precipitation (both 0.69) and has smaller RMSE values (52 mm) than winter does (78 mm).

Given the widespread occurrence of steppe during the LGM, we also undertook a separate evaluation of transfer function performance in this type of environment. For this we used a subset of 1588 pollen samples from EMPD2 that are classified with the steppe pollen biome (Davis et al., 2020). The results indicate (Table A4) little difference in performance compared to the full dataset, with a small decrease in performance in annual and summer seasons in both precipitation and temperature and a slight increase in performance in winter.

The results indicate good transfer function performance overall, especially for temperature, and are comparable with those found in other continental-scale pollen–climate studies (Bartlein et al., 2011). It is important to remember, though, that comparisons between studies can only be made with caution because results are often heavily dependent on the nature of the modern pollen dataset used as the training set (Juggins, 2013) and on the method used (Salonen et al., 2019; Brewer et al., 2008; Peyron et al., 2013).

3.3 Climate reconstruction

Reconstructed LGM temperatures indicate an overall mean annual cooling of $-\mathrm{7.2}\pm \mathrm{3.3}$ °C, with a greater cooling of around $-\mathrm{9.3}\pm \mathrm{4.5}$ °C in winter and $-\mathrm{5.0}\pm \mathrm{3.2}$ °C in summer (Fig. 5). All sites apart from Lake Van (site no. 62) in eastern Turkey show cooler temperatures at the LGM compared to in modern times (Fig. 6), and, even at this site, cooler conditions fall within the uncertainties. With greater cooling in winter compared to summer, the difference in temperature between winter and summer also increased (shown by positive anomalies) at most (but not all) sites (Fig. 6). This increase in continentality was around +4.2 °C on average across all sites (Fig. 5).

Figure 5 Pollen-based MAT reconstructions for LGM annual, winter, and summer temperature anomalies (uncertainties represent 1 standard deviation). Continentality represents the difference in temperature between summer and winter, with positive anomalies indicating an increase in the temperature difference between summer and winter. All values are expressed as anomalies compared with the present day. The green line indicates the mean for all the sites.

Figure 6 Maps of pollen-based MAT reconstructions for LGM annual, winter, and summer temperature anomalies (as shown in Fig. 9). Continentality represents the difference in temperature between summer and winter, with positive anomalies indicating an increase in the temperature difference between summer and winter. All values are expressed as anomalies compared with the present day.

We reconstruct an overall decline in mean annual precipitation of around $-\mathrm{91}\pm \mathrm{270}$ mm (−13 %) at the LGM. Most of this decline is in winter ($-\mathrm{38}\pm \mathrm{90}$ mm) (−21 %), while in summer a small increase is shown (10±57 mm) (6 %), although uncertainties are large (Fig. 7). Compared to temperature, there is significant seasonal and spatial variability in positive and negative precipitation anomalies (Fig. 8). Positive anomalies appear more predominant in eastern and southern Spain and in central eastern Europe in both summer and winter, while positive anomalies are found more generally in summer across sites in southern Europe and the Mediterranean. These more positive summer anomalies also reflect a relative shift from winter to summer in the seasonality of precipitation in this region.

Figure 7 Pollen-based MAT reconstructions for LGM annual, winter, and summer precipitation anomalies (uncertainties represent 1 standard deviation). Seasonality represents the difference in precipitation between summer and winter, with positive anomalies indicating an increase in summer precipitation compared to winter. All values are expressed as anomalies compared with the present day. The green line indicates the mean for all the sites.

Figure 8 Maps of pollen-based MAT reconstructions for LGM annual, winter, and summer precipitation anomalies (as shown in Fig. 11). Seasonality represents the difference in precipitation between summer and winter, with positive anomalies indicating an increase in summer precipitation compared to winter. All values are expressed as anomalies compared with the present day.

4 Discussion

Before we consider the results of our analysis, it is important to provide some context in terms of European LGM geography and environment, which was very different from today (Fig. 1). Major ice sheets covered Scandinavia and much of the UK, the Alps, and the Pyrenees. The sea level was 120 m lower, resulting in much of the North Sea and English Channel becoming dry land and in the European coastline extending over 100 km out into the Atlantic and Mediterranean, especially around the Bay of Biscay and the Adriatic Sea. The Black Sea was no longer connected to the Mediterranean, and it was smaller, with a water level around 100 m lower than today (Genov, 2016). These changes in sea- or water level had two main consequences, the first being that the marine sites were closer to land and therefore closer to (low-lying) terrestrial vegetation and (pollen-carrying) river discharge points than they are today. The second consequence of lower sea levels is that terrestrial pollen sites were located further from the moderating effect of the ocean than they are today, resulting in a localized modification of the climate experienced by the site irrespective of regional or global changes (Geiger, 1960).

The maps used in our analysis show the maximum ice sheet at 21k±2k (Ehlers et al., 2011). The precise geographical location of the ice sheet is difficult to resolve at a fine spatial scale, however, which explains why some sites close to the ice margin appear to be actually located under the ice (for example, Kersdorf-Briesen site no. 46 and Mickūnai site no. 54). The resolution of the map also shows the occurrence of permanent ice not only to the north and over the Alps but also on many subsidiary areas of high ground across central and southern Europe, including areas such as the Pyrenees, Massif Central, and the Vosges and Carpathian mountains. While global ice volume may have peaked ∼21 ka, individual ice sheets in Europe and other areas are known to have reached their maximum extent at different times (Hughes et al., 2016). The larger ice sheets are likely to have had a significant influence on regional climate and environmental conditions across Europe, but the smaller ice sheets had similar if more localized impacts as well. Surrounding each ice sheet would have been an unglaciated area of active periglacial processes and newly created and unstable ground. This would include outwash plains, impounded lakes and recently drained lake beds, seasonally and sporadically flooded areas, moraines, kettle holes, and other glaciological and periglacial features. Soils in these areas would have been non-existent or skeletal, and vegetation would have found it difficult to obtain nutrients and water for survival, irrespective of the prevailing climatic conditions. Outside of these areas, permafrost is also likely to have been present, particularly north of the Alps (Vandenberghe et al., 2014), which would also have acted as an impediment to vegetation growth.

In terms of regional climate, the major ice sheets would have provided significant barriers to westerly atmospheric circulation or even north–south circulation in the case of the Alps and the Pyrenees. As well as representing a physical obstruction, the thermodynamic response of the atmosphere to these high, cold obstructions would have been to encourage the formation of areas of semi-permanent high pressure, similar to those found today, for instance, over the Greenland ice sheet. In addition, the Laurentide ice sheet located over North America would have generated downstream effects over Europe (COHMAP, 1988). These physical and thermodynamic effects would have affected the direction of storm tracks and more local climatic effects commonly associated with ice sheets such as strong katabatic winds (Kageyama et al., 2021; Velasquez et al., 2021; Luetscher et al., 2015; Lefort et al., 2019)

4.1 Vegetation cover

The nature and extent of forest cover during the LGM remains a matter of considerable debate. Vegetation models driven by LGM climate model simulations generally indicate extensive areas of boreal forest north of the Alps and a mix of temperate and warm-temperate woodland to the south across southern Europe and much of the Mediterranean. Treeless areas such as steppe are mainly confined to those areas where it is also found today, namely inland Iberia, Ukraine, southern Russia, and Turkey, while tundra is found to the north close to the Scandinavian Ice Sheet (Allen et al., 2010; Cao et al., 2019; Prentice et al., 2011; Velasquez et al., 2021).

The evaluation of these vegetation model simulations against data has been largely based on comparisons with compilations of pollen biome reconstructions (Prentice et al., 2011; Allen et al., 2010; Cao et al., 2019; Velasquez et al., 2021). Early studies were based on only a limited number of sites in southern Europe and showed steppe at all sites in contradiction with model simulations (Elenga et al., 2000). More recent pollen compilations have included more sites, especially to the north, that have revealed a more mixed picture of vegetation cover, with forest biomes at some sites both south and north of the Alps that appear more consistent with model simulations (Binney et al., 2017; Cao et al., 2019). However, many of the pollen sites used in these studies were assigned an LGM age based on poor or incorrect dating control and likely date to MIS3, the Late Glacial, or even the Holocene. Nevertheless, based on our compilation of more securely dated LGM pollen sites, we also show a wider distribution of forest biomes particularly in Iberia, northern Italy, and central Europe, although with greater areas of steppe over the remaining regions than suggested by the models.

However, the interpretation of biome reconstructions requires care, since the forest cover and vegetation composition may not be as clear as the dominant biome suggests. For instance, we find that steppe is still reconstructed as the dominant biome at some sites despite arboreal pollen forming 70 %–80 % of the pollen assemblage. In addition, it is important to remember that pollen biomes are based only on the proportion of taxa that can form forest and woodland, while these taxa may in fact exist only as shrubs or stunted krummholz forms in the challenging climate and environment of the LGM. Alternatively, similar conditions may favour low-lying non-arboreal taxa forms with poor pollen dispersion or even insect-pollinated taxa forms that may be poorly represented in the pollen assemblage, giving greater prominence to arboreal taxa whose pollen may be the result of long-distance transport, particularly Pinus. However, there also appear to be plenty of samples with low or even very low (<20 %) arboreal percentages, so not all sites in open areas may be affected by long-distance transport of Pinus in the same way.

Quantitative MAT-based reconstructions of forest cover can overcome some of these problems, where they can be detected, based on the composition of the pollen assemblage when compared with the modern land cover. Chord distance measurements of the match between fossil and modern pollen assemblages indicate that good LGM analogues exist in our large Eurasian modern pollen dataset. The results of the MAT forest cover reconstruction indicate that forest cover was low but not entirely devoid of woodland in most areas, similar to the modern boreal forests of Siberia and consistent with a steppe–tundra–woodland mosaic proposed by many authors (e.g. Birks and Willis, 2008; Willis and Van Andel, 2004). This is confirmed in an analysis of the most commonly found modern analogue ecoregions for LGM pollen samples at each site (Table A5). Uncertainties are large, but for comparison the MAT site average of 33 % forest cover is slightly less than the average today over the Boreal region of Europe (43 %) and slightly more than the average today over the Mediterranean region (27 %) (Zanon et al., 2018).

By calculating the percentage of each of the taxa in each LGM pollen sample using a standardized pollen sum, we are able to make direct comparisons between different LGM pollen records and their taxa percentages (Figs. A3 and A4). The results show a preponderance of boreal forest taxa to the north of the Alps, consistent with biome results mentioned earlier. Pinus is the most common forest-forming taxa in this boreal zone, together with Picea, and including Larix to the east and Abies to the west. The occurrence of Betula and Juniperus also suggests shrubby elements consistent with arctic shrub tundra, although high Poaceae and other herbaceous taxa such as Artemisia and Chenopodiaceae indicate more steppe than tundra. Other deciduous taxa found north of the Alps include cold-tolerant generalists, such as Corylus and Alnus, and low percentages of relatively thermophilous taxa in the east, such as Carpinus and Tilia.

These results are consistent with charcoal (Magyari et al., 2014a; Willis and Van Andel, 2004), malacology (Juřičková et al., 2014), biomarkers (Zech et al., 2010), and genetic evidence (Stivrins et al., 2016; Willis and Van Andel, 2004) that the main forest region north of the Alps was in the eastern region of central Europe around the Carpathian Basin. This was also an area where cold- and moisture-sensitive deciduous taxa were also able to survive (Magyari et al., 2014b), although evidence of temperate taxa found in the pollen record has yet to be supported by charcoal and macrofossil records (Feurdean et al., 2014). Our pollen evidence indicates an open taiga or cool mixed forest that extended in central and eastern Europe to areas close to the Scandinavian and Alpine ice caps, as proposed by Willis and Van Andel (2004) and Huntley and Allen (2003), although whether this represents isolated pockets of forest or an extended open steppe forest is difficult to determine (Kuneš et al., 2008). Even steppe or tundra areas in western Europe show a low but significant presence of the pollen of tree taxa at sites close to the ice sheets that are unlikely to be solely the result of long-distance transport or reworking (Kelly et al., 2010). The presence of woodland in these areas is also supported by mammalian remains, for instance, at Kents Cavern in SW England (Stewart and Lister, 2001).

Overall, however, our results clearly show a much greater predominance of thermophilous and moisture-sensitive deciduous taxa south of the Alps, particularly in Iberia and northern Italy, where temperate broadleaf forests survived in sheltered refugia (Kaltenrieder et al., 2009). Most of these appear to be in hilly areas with the ability to generate orographic rainfall (Monegato et al., 2015) on south-facing slopes to make the most of the sun's radiant energy and are located above the valley floor to escape frost and flooding. We might also expect these areas to be sheltered from cold northerly winds and benefit from relatively mild and moisture-laden winds coming from the Mediterranean Sea. For instance, the presence of woodland and low glacier altitudes along the southern slopes of the Alps around the Po Valley and the Trentino region is consistent with strong orographic rains generated by southerly and easterly winds that today can be generated by low pressure located south of the Alps in the Gulf of Genoa and is consistent with a southerly storm track around the Alps (Kehrwald et al., 2010; Luetscher et al., 2015). Generally, as might be expected, areas of forest reconstruct similar or increased precipitation compared to today, and areas of steppe indicate decreased precipitation (see next section).

Independent evidence of LGM vegetation is provided by archaeozoological data. These data support the palynological evidence for the existence of forest and woodland refugia across the ice-free areas of Europe at latitudes north of the Alps. For instance, large vertebrates in these areas show patterns of extirpation and extinction in response to shifts in climate and vegetation cover that is different for different species, indicating a variety of environments and niches (Lister and Stuart, 2008; Stewart and Lister, 2001). As with the pollen record, the presence of temperate-adapted large-vertebrate taxa within the glacial landscape of western Europe also suggests the existence of temperate “micro-refugia” (Stewart and Lister, 2001), consistent with suggestions that temperate arboreal taxa were not entirely extirpated from the region during the LGM (Magri, 2010). Further east, mammal assemblages indicate generalized loss of forest components in the East European Plain (Demay et al., 2021; Puzachenko et al., 2021), which is consistent with our data indicating low forest cover in this region. In other areas, evidence of the prevailing land cover at the LGM comes from studies of small-vertebrate communities, which have a closer affinity to the prevailing environment than large vertebrates (López-García and Blain, 2020) that have the propensity to migrate large distances, often on a seasonal basis. These studies of small-vertebrate assemblages also support the existence of temperate “micro-refugia” in France (Royer et al., 2016) and the existence of woodland components in many regions across southern Europe, including parts of Iberia (Bañuls-Cardona et al., 2014) Italy (Berto et al., 2019), and the Balkan Peninsula (Mauch Lenardić et al., 2018).

Other palaeobotanical evidence also supports our land-cover reconstruction. Schäfer et al. (2016) suggest that leaf wax patterns from palaeosols in Spain may indicate the presence of drought-intolerant deciduous trees and more humid conditions during the LGM. Significantly, none of the pollen sites indicate that temperate broadleaf forests were dominant, and broadleaf temperate taxa always appear as part of a mixed woodland together with cold- or aridity-adapted evergreen and needleleaf taxa, including typical Mediterranean taxa. This type of mixed vegetation probably extended to the Balkans, where the hilly terrain and proximity to the Mediterranean would appear to have provided favourable climatic conditions, although we still lack LGM sites from this region. At sites in central and southern Italy and east through Greece and Turkey to the Middle East (and including northern Africa), the vegetation appears drier with a greater prevalence of steppe. Only a site in Georgia at the edge of the Caucasus Mountains indicates the presence of significant amounts of forest (mainly Pinus), a result that was also found by Tarasov et al. (2000) and is probably linked to favourable orographic precipitation and proximity to the Black Sea.

Comparisons with LGM land cover from vegetation modelling studies driven by climate model simulations indicate a much wider presence of forest than that shown by the pollen data (Kaplan et al., 2016). Data–model agreement appears to be closest over eastern central Europe, where pollen indicates the presence of open boreal forest, and over southwestern Europe with the presence of cool mixed temperate forest, including broadleaf deciduous and thermophilous elements (Prentice et al., 2011; Allen et al., 2010; Cao et al., 2019; Velasquez et al., 2021). Nevertheless, agreement still appears to be weak over western central Europe and southern and eastern Europe through to the Middle East, where pollen data continue to indicate widespread steppe. One proposed explanation for this data–model discrepancy is the role of fire (including man-made fire) in maintaining forest openness, a factor influencing forest cover that is not included in most vegetation models (Kaplan et al., 2016). In the Carpathian Basin, Magyari et al. (2014a) noted that charcoal increased as forest cover declined, suggesting that wildfires played a role in decreasing forest cover during the LGM. Other studies have noted low levels of charcoal and therefore fires during the LGM, although these tend to be from steppe areas with low biomass and fuel availability (Connor et al., 2013; Kaltenrieder et al., 2009). Recent LGM vegetation simulations that include fire indicate much lower values of forest cover than those without fire over western central Europe, while forest remains in central eastern Europe (Velasquez et al., 2021). This appears closer to the data, but the values are perhaps too low compared with our MAT reconstructions here (Fig. 4).

4.2 Climate

4.2.1 Comparison with previous pollen-based reconstructions

The climate of the LGM is generally regarded as having been cooler and drier than today, but data–model comparisons continue to highlight important discrepancies, not only in the degree of cooling and drying but also in their seasonal and spatial distribution. Data–model comparisons over Europe have mainly used pollen-based climate reconstructions, especially the Paleoclimate Modelling Intercomparison Project (PMIP/CMIP) (Kageyama et al., 2006, 2021; Bartlein et al., 2011; Harrison et al., 2015; Braconnot et al., 2012; Cleator et al., 2020). The most commonly used reconstructions are based on two main methods, the artificial neural network (ANN) methodology of Peyron et al. (1998) and Jost et al. (2005) and an inverse modelling (INV) approach applied by Wu et al. (2007). The ANN method uses modern pollen samples and does not include any correction for CO₂ effects, being similar in these respects to the MAT method used in this study. In contrast, the INV method does not use modern pollen samples but instead uses a process-based vegetation model run in inverse mode (Guiot et al., 2000). Ordinarily, a vegetation model will use climate as an input to generate a type of vegetation as an output, but in inverse mode the model is reconfigured so that the input climate (and CO₂) can be varied iteratively until the closest match is found between the vegetation simulated by the model (represented by PFT scores) and the fossil pollen assemblage (also represented by PFT scores). One of the advantages of the INV method is that CO₂ can also be varied as an input, and therefore the effect of changes in CO₂ on the vegetation, and therefore reconstructed climate, can be investigated. Comparisons of these ANN and INV reconstructions have shown important differences, with the INV reconstruction generally not as cold as and somewhat drier than ANN (Wu et al., 2007). These differences between pollen–climate methods have often been attributed to CO₂ effects (Wu et al., 2007) but this is not clear, since there may be other factors, such as the size and location of the training set used in the ANN reconstruction.

Figure 9 A site-by-site comparison between LGM pollen–climate reconstructions based on the modern analogue technique (MAT; this study), artificial neural networks (ANNs; Peyron et al., 1998), and inverse modelling (INV; Wu et al., 2007). The results show that MAT and INV give similar climate reconstructions but that ANN is significantly cooler and drier.

We make a comparison with these earlier reconstructions based on 10 sites/records in our dataset which we identified as also being included in these earlier studies (Fig. 9). While we were able to identify the site and data source and the time window, we were unable to establish if the the data represented a single sample or the mean of multiple samples within a time window, the exact depth of those samples, or the actual sediment core in the case of multiple cores from the same site. While these aspects are unknown, it seems likely that the pollen data we used in our analysis were very similar, if not identical, in most cases and that reconstructed biomes for these sites from our pollen dataset are identical to the biomes reconstructed using the earlier pollen dataset (Elenga et al., 2000).

We compare our MAT with the ANN and INV reconstructions in Fig. 9. On average across all 10 records, the MAT and INV methods give almost identical results for anomalies of both mean annual temperature (MAT −6.6 °C, INV −7.2 °C) and precipitation (MAT 158 mm, INV 165 mm). Uncertainties are also similar for both methods. In contrast, the ANN method gives much cooler mean annual temperature anomalies (ANN −13.9 °C) and drier precipitation anomalies (ANN −474 mm). On a site-by-site basis, the MAT and INV methods show closer agreement for temperatures than precipitation, although precipitation has proportionally larger uncertainties. The reconstructions based on these two methods are close enough that the uncertainties overlap at all sites for both temperature and precipitation, except the precipitation reconstruction at Lac de Bouchet (site no. 25). The reason for this is not clear, but there could easily be minor differences with the pollen data analysed by Wu et al. (2007) in their INV reconstruction, since the pollen record (Reille and de Beaulieu, 1988) includes multiple cores, each with many different samples covering the LGM period.

This comparison shows that our MAT reconstructions are very similar to the INV method but are not as cold or dry as the ANN method. This has two main implications. The first is that our reconstructions indicate greater agreement with the results of climate model simulations, since climate models indicate temperatures closer to the INV reconstructions (Latombe et al., 2018) than to the ANN reconstructions (Jost et al., 2005; Kageyama et al., 2006). The difference between our MAT and earlier ANN reconstructions is likely the result of the modern pollen datasets used, since the ANN reconstruction was based on a considerably smaller number of samples taken mainly from the cold dry steppes of Kazakhstan and Mongolia.

The second implication is that the MAT method may not be significantly impacted by the effects of lower CO₂ (Cowling and Sykes, 1999; Prentice and Harrison, 2009; Williams et al., 2000) or indeed insolation changes during the LGM, since the MAT results are similar to those based on the INV method, which specifically takes account of these non-climatic factors (Wu et al., 2007). This would suggest that the MAT could also work well for pollen-based climate reconstructions on longer glacial–interglacial timescales where insolation and CO₂ vary significantly from their modern values. This is consistent with the findings of Pini et al. (2022), who applied a correction algorithm developed by Prentice et al. (2017) and Cleator et al. (2020) to a MAT reconstruction of mean annual precipitation at Lake Fimon in northern Italy. This shows a very small correction of 0 to 30 mm for samples across the LGM time window, which indicates that CO₂ is not a very significant factor in influencing this type of reconstruction, at least compared to the overall uncertainties (±200 mm) of the reconstruction itself. The uncertainties associated with the correction algorithm are not discussed, but, given that inputs include estimates of both LGM temperature and cloud cover, it seems likely that these could be significant. Importantly, both Pini et al. (2022) and Cleator et al. (2020) specifically exclude the necessity of applying a correction algorithm to temperature reconstructions, since they only regard hydrological variables as being affected by changes in atmospheric CO₂.

4.2.2 Comparison with climate reconstructions based on other proxies

Temperature

Proxies that are not based on plants should remain unaffected by the CO₂ problem during the LGM and provide an alternative basis for evaluating pollen-based reconstructions. Samartin et al. (2016) reconstructed LGM summer temperatures based on chironomid remains from Lago della Costa (site no. 34) in northern Italy. They also undertook pollen analysis on the same samples down the core, allowing us to make a sample-by-sample comparison between the chironomid temperature record and our MAT reconstruction (Fig. 10). Our pollen–climate reconstruction is for JJA mean temperature, while the chironomid reconstruction is for July mean temperature, with the anomalies based on the modern equivalent JJA and July mean temperatures respectively. The average anomaly values for all eight samples reconstructed by the pollen–climate MAT are $-\mathrm{10.2}\pm \mathrm{3.5}$ °C and are $-\mathrm{9.5}\pm \mathrm{3.0}$ °C for the chironomids. This indicates that pollen and chironomid average summer temperature reconstructions are very similar on average, taking into account the overlapping uncertainties, while also showing a strong similarity on a sample-by-sample basis throughout the time series.

Figure 10 Comparison between the LGM pollen–climate MAT and chironomid summer temperature reconstructions at Lago della Costa, Italy (chironomid reconstruction and pollen data from Samartin et al., 2016). Dashed lines show uncertainties.

Other reconstructions based on other proxies provide a basis for more general regional comparisons (Figs. A6 and A10). We reconstruct both summer and winter temperatures and show that cooling in winter was greater than in summer at most sites, associated with an increase in continentality (increased temperature difference between summer and winter). A similar seasonal pattern of temperature change has also been shown in other studies that reconstruct both summer and winter LGM temperatures, including Prud'homme et al. (2016) using δ¹⁸O analysis of earthworm calcite granules at Nußloch near the French–German border; Bañuls-Cardona et al. (2014) using faunal remains of small mammals at four locations in western Spain; and Ferguson et al. (2011), who examined seasonal temperature change using δ¹⁸O and Mg–Ca analysis of limpet shells at Gibraltar in southern Spain. The increase in continentality at Nußloch (Prud'homme et al., 2016) was reconstructed at between 11.6 to 15.6 °C, comparable at the lower end with nearby pollen sites: La Grotte Walou, site no. 28, with an increase of 10.4±5.8 °C and Bergsee, site no. 29, with an increase of 7.9±5.7 °C. The faunal sites in western Spain studied by Bañuls-Cardona et al. (2014) gave much reduced increases in continentality that were nevertheless similar to nearby pollen sites. They include the following: Valdavara, 5.1 °C (similar to MD99-2331, site no. 3, 5.2±3.1 °C); El Mirón, 1.2 °C (similar to Tourbière de l'Estarrès, site no. 19, 5.1±6.2 °C); El Portalón, 0.9 °C (similar to Torrecilla de Valmadrid, site no. 16, 2.8 ±1.8 °C); and Cueva de Maltravieso, 6.1 °C (similar to SU81-18, site no. 2, 4.8±3.4 °C). Further south at Gibraltar, the limpet-based study of Ferguson et al. (2011) also shows a relatively small increase of 2 °C. The nearest pollen site (Gorham's Cave, site no. 5), however, shows a larger increase of 4.7±2.3 °C, although differences could be expected given the different temporal resolution of annual laminae on mollusc shells compared to pollen assemblages that reflect much slower changes in trees and other long-lived flora.

Summer temperatures were warm enough during the LGM over the Alpine areas that Swiss lakes were largely ice-free in summer, while glacier equilibrium-line altitude (ELAs). around the time of the LGM suggest summers were −6.5 to −7.7 °C cooler compared to the Little Ice Age (LIA) (Heiri et al., 2014). This cooling was similar to that found at Nußloch some 200 km north of the Swiss border by Prud'homme et al. (2016), who reconstructed anomalies of −6 to $-\mathrm{8}\pm \mathrm{4}$ °C from δ¹⁸O analysis of earthworm calcite granules, representing warm-season temperatures (May–September). Slightly less cooling was found close by at the nearby site of Achenheim, where the analysis of mollusc assemblages gave summer (August) cooling estimates of −3.5 to −6.5 °C based on the MAT (Rousseau, 1991) and −5.5 to −9.5 °C based on the mutual climatic range method (Moine et al., 2002). These reconstructions appear somewhat cooler than nearby pollen sites (La Grotte Walou, site no. 28, at $-\mathrm{1.4}\pm \mathrm{3.6}$ °C and Bergsee, site no. 29, at $-\mathrm{2.7}\pm \mathrm{5.1}$ °C) although comparable with the pollen site Pilsensee, site no. 32 ($-\mathrm{7.3}\pm \mathrm{5.0}$ °C), 200 km further east. Similar differences also occur at the site of Les Échets on the western edge of the Alps, where a diatom-based reconstruction of summer (July) temperatures (Ampel et al., 2010) indicated a greater cooling (−10.5 to −11.5 °C) than our pollen reconstruction at Les Échets G, site no. 27 ($-\mathrm{4}\pm \mathrm{2.7}$ °C). However, the authors caution that the results were based on poor analogues and rare taxa and on a small training set of only 90 lakes in Switzerland.

South of the Alps, other proxies show the opposite relationship with the pollen reconstructions. For instance, at Lago della Costa in the Po Valley, a summer (July) temperature chironomid reconstruction by Samartin et al. (2016) is around 1–2 °C less cool than the pollen reconstruction (JJA) for the same site, Lago della Costa, site no. 34 ($-\mathrm{11.4}\pm \mathrm{2.7}$ °C), although both reconstructions fall within their respective uncertainty ranges (Fig. 8). In the Pindus mountain range in Greece, Hughes et al. (2006) estimated LGM summer temperature anomalies of −7 °C based on glacier modelling, which is comparable with those reconstructed at the nearest pollen site, Ioannina, site no. 51 ($-\mathrm{7.7}\pm \mathrm{2.8}$ °C). In Spain, the analysis of small-mammal remains by Bañuls-Cardona et al. (2014) shows similarly less cooling in summer or even warmer conditions than present positive anomalies compared to the nearest pollen sites, such as Valdavara, 1.4 °C (similar to MD99-2331, site no. 3, $-\mathrm{2.3}\pm \mathrm{2.8}$ °C); El Mirón, −2.3 °C (similar to Tourbière de l'Estarrès, site no. 19, $-\mathrm{5.7}\pm \mathrm{5.4}$ °C); El Portalón, 0.8 °C (similar to Torrecilla de Valmadrid, site no. 16, $-\mathrm{2.6}\pm \mathrm{1.1}$ °C); and Cueva de Maltravieso, −1.1 °C (similar to SU81-18, site no. 2, $-\mathrm{10.4}\pm \mathrm{2.8}$ °C). Further south at Gibraltar, the limpet-based study of Ferguson et al. (2011) suggests an anomaly of around −7 °C, which is a greater cooling than the pollen reconstruction from this location (Gorham's Cave, site no. 5, $-\mathrm{1.3}\pm \mathrm{2.2}$ °C), although it is comparable with other pollen sites slightly further east.

Winter temperature reconstructions from non-pollen proxies show a similar pattern to summer temperatures in relation to pollen reconstructions. North of the Alps at Achenheim, Prud'homme et al. (2016) use δ¹⁸O in earthworm remains to reconstruct particularly cold winter anomalies of −17.6 to −23.6 °C compared to nearby pollen sites (La Grotte Walou, site no. 28, with anomalies of $-\mathrm{11.8}\pm \mathrm{8.0}$ °C and Bergsee, site no. 29, with anomalies of $-\mathrm{10.6}\pm \mathrm{6.3}$ °C). South of the Alps in Spain, the analysis by Bañuls-Cardona et al. (2014) based on the remains of small mammals shows less cooling in winter compared to the nearest pollen sites, in particular, Valdavara, −3.7 °C (similar to MD99-2331, site no. 3, $-\mathrm{7.5}\pm \mathrm{3.4}$ °C); El Mirón, −3.5 °C (similar to Tourbière de l'Estarrès, site no. 19, $-\mathrm{10.8}\pm \mathrm{7.0}$ °C); El Portalón, −0.1 °C (similar to Torrecilla de Valmadrid, site no. 16, $-\mathrm{5.4}\pm \mathrm{2.5}$ °C); and Cueva de Maltravieso, −7.2 °C (similar to SU81-18, site no. 2, $-\mathrm{15.2}\pm \mathrm{4.0}$ °C). And, again, in southern Spain at Gibraltar, the analysis of limpet shells by Ferguson et al. (2011) suggests winter cooling of around −9 °C, while the pollen reconstruction suggests $-\mathrm{6.0}\pm \mathrm{2.5}$ °C for Gorham's Cave, site no. 5, although sites further east indicate cooler conditions.

A number of additional proxies have also been used to reconstruct LGM mean annual temperature. Heyman et al. (2013) applied glacier mass balance modelling at sites located in the smaller mountain regions north of the Alps. These are generally slightly cooler than our pollen-based reconstructions at sites close to the Vosges Mountains, $-\mathrm{12.7}\pm \mathrm{2.0}$ °C, and the Black Forest, $-\mathrm{11.4}\pm \mathrm{2.3}$ °C (similar to Bergsee, site no. 29, $-\mathrm{8.2}\pm \mathrm{3.3}$ °C); the Bavarian Forest, $-\mathrm{10.7}\pm \mathrm{2.2}$ °C (similar to Pilsensee, site no. 32, $-\mathrm{9.2}\pm \mathrm{1.2}$ °C); and the Giant Mountains, $-\mathrm{8.5}\pm \mathrm{1.8}$ °C (similar to Kersdorf-Briesen, site no. 46, $-\mathrm{7.3}\pm \mathrm{0.3}$ °C). These values obtained by Heyman et al. (2013) are warmer than those obtained by Prud'homme et al. (2016), who estimated annual mean temperature anomalies of −15.1 to −19.1 °C based on δ¹⁸O of earthworm calcite at the Nußloch site just north of the Vosges mountains and the Black Forest. The annual temperatures reconstructed by Heyman et al. (2013) are also around 2 °C warmer than those reconstructed by Allen et al. (2008a), who applied a similar, although simpler, method to over 29 different mountainous regions across Europe that were glaciated during the LGM. Since glacier mass balance is a function of both snowfall and temperature, these estimated temperatures vary according to estimated changes in precipitation. For instance, mean annual temperature estimates by Allen et al. (2008a) are much cooler than those reconstructed by pollen, with an average anomaly of −13.2 °C for the 29 sites assuming a 40 % reduction in precipitation, but this is reduced to −11.8 °C assuming the same precipitation as during the modern day. These average anomalies across all sites calculated by Allen et al (2008a) compare with an average temperature anomaly of −7.2 °C across all 63 of our pollen sites. The glacier mass balance modelling by Allen et al. (2008a) assumes a seasonal distribution of precipitation that is similar to the present day, and it does not consider increases in winter precipitation or mean annual precipitation above present-day levels. Both of these are suggested by the pollen data in some regions, and both could explain glacier extent found during the LGM based on less-extreme temperature anomalies more comparable with the pollen data. Mean annual temperatures have also been reconstructed from the Paris Basin area in eastern France by Bekaert et al. (2023) using the noble gas proxy. The authors suggest an LGM temperature anomaly of $-\mathrm{9.1}\pm \mathrm{0.9}$ °C, although this is actually dated to 25.6±0.5k, which is earlier than our 21k±2.0k time window that we adopt here. The sample closest to 21 k is at 21.9±0.5k and suggests slightly warmer temperatures at −7.77 °C, which compares well with our pollen reconstructions nearby at Bergsee, site no. 29 ($-\mathrm{8.2}\pm \mathrm{3.3}$ °C), and La Grotte Walou, site no. 28 ($-\mathrm{6.6}\pm \mathrm{3.1}$ °C).

To the east of the Alps in the Pannonian Basin, mean annual temperature anomaly estimates have been made from noble gas measurements on groundwater ranging from −2 to −4 °C (Stute and Deak, 1990) up to −9 °C (Varsányi et al., 2011). These are similar to estimates ranging from −2 to −9 °C from oxygen isotope ratios from mammoth tooth enamel (Kovács et al., 2012) and are comparable with nearby pollen sites: Lake Fehér, site no. 50 ($-\mathrm{8.2}\pm \mathrm{3.3}$ °C), and Kokad, site no. 52 ($-\mathrm{4.5}\pm \mathrm{2.3}$ °C). On a broader scale, Sanchi et al. (2014) estimated LGM cooling in the Danube and Dnieper basins based on lipid biomarkers in a core from the Black Sea and came up with similar mean annual temperature anomalies between −6 to −10 °C, which again are comparable with pollen sites from the region that range from Nagymohos, site no. 48 ($-\mathrm{10.5}\pm \mathrm{4.1}$ °C), to Straldzha, site no. 57 ($-\mathrm{4.3}\pm \mathrm{5.8}$ °C).

Further south and west, García-Amorena et al. (2007) reported mean annual temperature anomalies of −2.0 to −11.3 °C at LGM sites along the Portuguese coast, based on an indicator species method using plant macrofossils. This is similar to the closest marine pollen sites off the coast, which recorded values of $-\mathrm{10.5}\pm \mathrm{4.6}$ °C at MD95-2039, site no. 1, and $-\mathrm{5.3}\pm \mathrm{2.9}$ °C at MD99-2331, site no. 3. Meanwhile, in the far east of the study area, Zaarur et al. (2016) estimated a mean annual temperature anomaly of around −3 °C based on a clumped isotope analysis of Melanopsis shells from LGM sediments in the Sea of Galilee. This limited cooling appears similar to the nearest pollen site (Lake Zeribar, site no. 63), where we reconstruct a cooling of $-\mathrm{2.2}\pm \mathrm{4.6}$ °C.

Reconstructions of LGM sea surface temperatures (SSTs) provide yet another source of comparison with our terrestrial pollen-based reconstructions, although many of the physical processes controlling sea surface temperatures, such as upwelling, surface mixing, surface currents, stratification, and thermal inertia through the seasonal cycle, represent quite different processes to those controlling surface temperatures over land, particularly at the sub-regional scale. Nevertheless, the Atlantic coastal waters of Iberia and the waters throughout the Mediterranean Sea include many SST sites that lie in relative proximity to our terrestrial pollen sites, allowing us to make a comparison at the largest scale. Within this area, the MARGO database (MARGO Project Members, 2009) includes 13 alkenone-, 2 Mg–Ca-, and 41 foraminifera-based SST records of mean annual temperature, with the foraminifera records also providing an additional 41 winter (JFM) and summer (JAS) SST estimates. We compare the SST records with the 36 closest terrestrial pollen records which fall within a box of −11 to 35° longitude and 32 to 43° latitude containing all of the SST records. A simple site average indicates a mean annual SST anomaly of $-\mathrm{5.5}\pm \mathrm{1.0}$ °C, which is relatively close to the value of $-\mathrm{7.2}\pm \mathrm{3.4}$ °C obtained from the terrestrial pollen sites (site nos. 1–4, 5, 7–24, 25, 26, 30, 35–38, 41, 47, 51, 53, and 56–59). Interestingly the inter-site variance (standard deviation of the reconstructed temperatures across all sites) is almost identical for the two datasets, 2.57 °C for the SST sites and 2.63 °C for the pollen sites, despite representing very different environments, proxies, and uncertainties. However, when we look at the seasonal temperature anomalies, we find very different results. Site-averaged winter SST anomalies are $-\mathrm{3.7}\pm \mathrm{1.1}$ °C compared to $-\mathrm{9.3}\pm \mathrm{4.2}$ °C for winter temperatures from terrestrial pollen sites, while in summer the values are reversed ($-\mathrm{7.0}\pm \mathrm{0.8}$ °C compared to $-\mathrm{5.38}\pm \mathrm{3.3}$ °C respectively). This suggests that SSTs experienced greater cooling in summer compared to winter, which is the opposite to what is generally found in terrestrial seasonal temperature reconstructions throughout the region, although this is consistent with model simulations (Mikolajewicz, 2011).

Precipitation

Few proxies apart from pollen provide quantitative reconstructions of precipitation during the LGM. Glacier mass balance modelling includes assumptions about precipitation in order to derive temperatures (Allen et al., 2008a), but neither is independent of the other. Hughes et al. (2006) estimate from glacier modelling that the mean annual precipitation during the LGM at sites in the Pindus mountains in Greece was around 2300±200 mm, which they regard as similar to the present day (>2000 mm). A small change in precipitation compared to modern values is also indicated by the nearest pollen site, which is around 47 km to the south (Ioannina, site no. 51), and indicates a mean annual precipitation anomaly of $-\mathrm{152}\pm \mathrm{294}$ mm, representing just 15 % of the modern value. A larger reduction in mean annual precipitation of −45 % (maximum) is reconstructed by García-Amorena et al. (2007) based on plant macrofossil remains from sites on the Portuguese coast. In comparison, the closest pollen sites record values which are a little lower, ranging from −22 % at MD95-2039, site no. 1, to −34 % at MD99-2331 site no. 3. Further north, in southwestern Germany, Prud'homme et al. (2018) reconstructed mean annual precipitation from the δ¹³C of earthworm calcite granules at Nußloch. They estimated a field site average of 333 (159–574) mm yr⁻¹ at the LGM, which represents an anomaly of −503 mm yr⁻¹ (−60 %) relative to the modern precipitation of 836 mm yr⁻¹. This is comparable with the closest pollen site (Bergsee, site no. 29) with an anomaly of −540 mm yr⁻¹.

As with glaciers, lake levels reflect changes in moisture balance that include the effects of both temperature (via evapotranspiration) and precipitation rather than just precipitation. They also represent semi-quantitative data at best, with changes often described relative to the modern or other baseline. There are few lake-level records available north of the Alps, but, to the south, many records indicate high lake levels in areas such as Spain (Lacey et al., 2016; Moreno et al., 2012; Vegas et al., 2010), Italy (Belis et al., 1999; Giraudi, 2017), Greece and Turkey (Harrison et al., 1996; Reimer et al., 2009), and the Middle East (Kolodny et al., 2005; Lev et al., 2019). These lake records are also supported by the evidence of higher river levels in Morocco (El Amrani et al., 2008). The cause of the higher lake levels has been the subject of some debate, since many pollen records (and especially early biome reconstructions) show steppe vegetation that would suggest aridity that appears incompatible with higher lake levels. Prentice et al. (1992) proposed that the co-existence of steppe vegetation and high lake levels could be possible if precipitation increased outside of the summer growing season, while summers themselves were drier and cooler with decreased evaporation. However, the results of our analysis tend to indicate the opposite in regions with higher lake levels, with increased summer rainfall and decreased winter rainfall. In addition, the increase in summer precipitation was enough to compensate for the decrease in winter rainfall, leading to an overall increase in mean annual precipitation at many pollen sites in Spain and Greece, for instance. This, together with depressed temperatures and consequently decreased evaporation, could explain the higher lake levels whilst also limiting the growth of trees as a result of cooler temperatures and prolonged aridity outside of the summer season. Davis and Stevenson (2007) also note a differential hydrological response between summer and winter rainfall in the Mediterranean during the Holocene that may also provide an explanation. In this case, sporadic summer storms may result in high rates of runoff that may fill runoff-fed lakes, but they may also result in low rates of soil moisture recharge that fail to benefit vegetation in the same way winter rainfall does.

Overall, we reconstruct only a small reduction in precipitation during the LGM of around 91 mm (13 %) averaged over all sites, which is less than the ∼200 mm reduction based on the sites in the pollen–climate compilation used by PMIP (Bartlein et al., 2011). Since our precipitation reconstruction on average matches that of the INV reconstruction by Wu et al. (2007), we can attribute much of the difference to the greater aridity shown in the ANN reconstruction by Peyron et al. (1998) and Jost et al. (2005) (see Fig. 9). As with temperature, this is probably a reflection of the modern training set used in the ANN reconstruction, which is much smaller than our training set and is largely taken from the arid steppes of Kazakhstan and Mongolia. However, it is also important to recognize the significant spatial variability in precipitation, which means that a simple average of different sets of sites from different regions may not accurately reflect the change in LGM precipitation at the European scale. Nevertheless, one of the most consistent signals in our dataset is for an increase in summer precipitation over many areas of southern Europe and the Mediterranean. This is also found in climate models, where it has been attributed to an increase in convection-driven precipitation, although the amount of precipitation generated by this mechanism varies significantly between models (Beghin et al., 2016). It may seem counter-intuitive to see an increase in reconstructed precipitation in the same regions where we also find a preponderance of steppe or xerophytic biomes and taxa, including Artemisia and Chenopodiaceae. This is attributable to the fact that climate can change quite markedly, necessarily invoking a major change in vegetation and especially in the pollen biome. For instance, a semi-arid climate ranges from 250–500 mm of rainfall a year, so we could expect semi-arid vegetation to be dominant even if the rainfall increases by 250 mm (100 %).

A more consistent response in models is for an increase in winter precipitation across southern Europe and the Mediterranean related to a stronger and more southerly displaced jet stream, with winter precipitation also accounting for much of the change in mean annual precipitation (Beghin et al., 2016). Our reconstruction of winter precipitation, however, shows less support for this scenario, with a more general decrease in winter precipitation apart from in southern and eastern Iberia and with summer precipitation generally more important in those sites that show an overall increase in mean annual precipitation. This may not necessarily contradict the models in terms of the strength and position of the winter jet stream, but it may instead indicate that models overestimate the amount of moisture being carried westward from the cold North Atlantic along the storm track, especially across the far northern Mediterranean. The increase in winter precipitation across southern and eastern Iberia is, however, entirely consistent with a strengthened and more southerly jet stream, which also brings increased winter precipitation to the region today as a result of blocking over northern Europe and the Atlantic and a negative North Atlantic Oscillation (NAO) (Vicente-Serrano et al., 2011).

Other areas that show an increase in winter precipitation include pollen sites around the eastern end of the Alps. This is consistent with a recent study by Spötl et al. (2021), who argued, on the basis of cryogenic carbonates preserved in a cave in Austria, that heavy winter (and autumn) precipitation was a significant factor in driving LGM glaciation in the region. The seasonally specific nature of this precipitation is also supported by the same pollen sites, which do not show any increase in summer precipitation at this time.

5 Conclusions

We have reconstructed the climate and vegetation cover across Europe, northern Africa, and the Middle East at the time of the LGM based on 63 pollen records. These records were selected using strict quality-control criteria, with particular attention paid to dating control, which led to the exclusion of many records that have been used in previous studies. This fully documented dataset represents the most chronologically precise and spatially resolved view of LGM climate and vegetation during the PMIP benchmarking time window at 21±2 ka. Nevertheless, it is important to recognize that there are still significant spatial gaps in pollen sites, especially north of the Alps, the Balkans, Turkey, and the Middle East, and we continue to have only a partial understanding of the LGM over these areas.

One of the key questions concerning the vegetation landscape of the LGM in Europe has been the extent to which forest rather than steppe covered the continent and to what extent temperate elements could be found north of the classical refugia areas of southern Europe and the Mediterranean. Our results show that, although steppe and tundra were extensive at the time of the LGM, areas of open forest also occurred in many regions, particularly (but not exclusively) in Iberia, northern Italy, and central Europe. These forest or woodland stands are likely to have been located in environmentally favourable areas, with good soils, elevated rainfall, and shelter from cold, desiccating winds. In those areas where woodland existed, boreal taxa generally dominated north and east of the Alps, while temperate and thermophilous (mainly drought-adapted) taxa were generally confined to areas south of the Alps and around the Mediterranean. The temperate deciduous forests that compose the climax community in many areas of Europe today were displaced to the south and reduced to a partnership role with boreal elements. Overall, our new reconstruction indicates greater agreement with model land-cover simulations, but models still appear to overestimate the amount of forest and woodland over areas such as France, the Benelux, Greece, Turkey, and the Far East.

Another key question about the LGM concerns the ability of climate models to simulate the climate of this period and whether pollen-based climate reconstructions which show disagreement with models have been biased by the effects of low CO₂ on plant physiology. We find that our new pollen–climate reconstruction shows much closer agreement with climate models than previous reconstructions that did not take account of low-CO₂ effects. We also find close agreement with previous reconstructions that did take account of CO₂ effects. Since our MAT method itself does not specifically take account of low-CO₂ effects, this would suggest that this problem is not a significant hindrance to MAT performance at the time of the LGM, at least not compared to other uncertainties. Instead, we suggest that the main factor in the performance of pollen–climate transfer functions that use modern analogue methods is the provision of a large enough modern pollen dataset with suitable LGM analogues.

This conclusion is supported by comparison with climate reconstructions based on other proxies. We found little difference between our MAT reconstruction and a chironomid-based summer temperature record based on a downcore sample-by-sample comparison and on comparisons with records from a variety of other proxies at a regional scale. However, it is notable that some studies using glacier mass balance modelling methods indicate LGM temperatures that are much cooler than our pollen-based reconstruction. The reasons behind this are not clear, but our pollen-based results indicate higher-than-present precipitation in some areas that could potentially explain low-elevation glacier ELAs without the need for such cold temperatures.

We also find that, although our pollen-based reconstruction and reconstructions of SSTs generally agree in terms of mean annual temperatures, SSTs indicate greater cooling in summer compared to winter, while terrestrial records indicate greater cooling in winter compared to summer. These seasonal differences are also reproduced in climate models and probably reflect the different processes driving seasonal temperature change in the terrestrial and marine domains.

Our reconstructions of precipitation show large spatial and seasonal variability but generally indicate less overall aridity than previously suggested from smaller-scale studies which sampled less of the spatial domain. We find that, in some regions of southern Europe, precipitation may actually have been greater than at present, especially in summer but also in winter in southern and eastern Iberia and around the southern slopes of the Alps. This may have important implications in understanding the development of LGM glaciation, which may be less a function of temperature than previously supposed. This could also help better explain the observed asynchronous nature of glaciation, even within relatively small regions such as Europe, as a result of more localized controls on ice sheet development such as precipitation.

We hope that this new continental-scale dataset of climate and vegetation reconstructions will provide an improved baseline for data–model comparisons and other studies that will allow us to better understand the complex LGM environment.

Appendix A

Table A1 Modern climate values for each site used in the calculation of anomalies (taken from WorldClim 2, Fick and Hijmans, 2017).

Table A2 This shows the results of the experiment to test the sensitivity of pollen biomes to changes in the number of Pinaceae in the pollen assemblage using 8213 modern pollen samples from EMPD2. Pinaceae can be over-represented in marine samples, and it has been proposed that removing all Pinaceae from these samples is better than leaving the Pinaceae in the pollen assemblage. The “Control” column on the left shows the number of samples that were classified for each biome without changing the number of Pinaceae (i.e. using the original pollen assemblage). The other eight columns to the right show the number of samples where the biome changed relative to the number shown in the control column as a result of either the removal of all Pinaceae (“0 Pinaceae”) or by artificially increasing the number of Pinaceae from 5 % to 400 % of the original count (“+5 % Pinaceae” to “+400 % Pinaceae”). For instance, for the CLDE (cold deciduous) biome, 25 pollen samples were classified as CLDE without any change in Pinaceae (“Control”), but 454 more samples were classified as CLDE when all Pinaceae were removed (“0 Pinaceae”) compared to 4 fewer samples that were classified as CLDE when the number of Pinaceae was increased by as much as 400 % (“+400 % Pinaceae”). The totals along the bottom show that, out of the 8213 pollen samples included in the experiment, 5860 biomes changed when all Pinaceae were removed, compared to up to 2348 when the number of Pinaceae was artificially increased by up to 400 %.

Table A3 A comparison of the LGM reconstructed climate for marine sites showing the effect of excluding Pinaceae (shaded) from the pollen assemblage compared to the results of including Pinaceae (unshaded; also presented in Figs. 6–8). It has been proposed that, because of the potential for over-representation of Pinaceae in marine pollen samples, it is better to exclude Pinaceae completely from marine pollen samples. Comparing the two approaches, temperatures are generally 1–2 °C cooler and precipitation is slightly higher when Pinaceae are excluded. The differences for both temperature and precipitation are significantly less than the standard deviation of their uncertainties.

Table A4 A comparison of MAT performance statistics based on the modern pollen sample training set using all surface samples from EMPD2 used in the LGM reconstruction (as shown in Table 3) and a subset of 1588 samples from EMPD2 that were classified as steppe. The results show little difference between the two different types of sample. The table includes mean annual temperature and precipitation (TANN and PANN), mean winter temperature and precipitation (TDJF and PDJF), and mean summer temperature and precipitation (TJJA and PJJA).

Table A5 The biome and ecoregion of the modern surface samples used as analogues in the pollen–climate reconstructions.

Notes: modern analogue biomes and ecoregions were calculated as the most commonly occurring amongst all six best modern analogue pollen samples in all LGM samples for each pollen site/record. These are taken from EMPD2 (Davis et al., 2020) using the classification of Olson et al. (2001).

Figure A1 Chronological control.

Figure A2 Pollen biomes (see Fig. 2 for key), percentage of arboreal pollen (AP) forest cover, percentage of MAT forest cover, and residuals (percentage of AP compared to percentage of MAT forest).

Figure A3 Pollen taxa percentages for all LGM sites/records.

Figure A4 Percentage maps of Abies, Alnus, Betula, Carpinus, Cedrus, and Corylus.

Figure A5 Percentage maps of Ephedra, Ephedra fragilis, Fagus, Fraxinus, Hippophae, and Juniperus.

Figure A6 Percentage maps of Larix (×10 exaggeration), Olea, Picea, Pinus, Quercus (deciduous), and Quercus (evergreen).

Figure A7 Percentage maps of Salix, Tilia, Ulmus, Apiaceae, Artemisia, and Asteraceae (Asteroideae).

Figure A8 Percentage maps of Asteraceae (Cichorioideae), Brassicaceae, Caryophyllaceae, Chenopodiaceae, Helianthemum, and Plantago.

Figure A9 Percentage maps of Poaceae, Rubiaceae, and Thalictrum.

Figure A10 Maps of pollen-based MAT reconstructions for LGM annual, winter, and summer temperature anomalies (as shown in Fig. 5), shown with the results of other published studies. Continentality represents the difference in temperature between summer and winter, with positive anomalies indicating an increase in the temperature difference between summer and winter. All values are expressed as anomalies compared with the present day unless otherwise indicated.

Figure A11 Map of pollen-based MAT reconstructions for LGM annual precipitation anomalies (as shown in Fig. 7), shown with the results of other published studies. All values are expressed as anomalies compared with the present day.