Unravelling the tree cover dynamics over the last 20 000 years on the Northern Hemisphere

Over the last 20 000 years, Northern Hemisphere vegetation underwent major shifts in response to orbital changes, rising CO₂, and ice sheet retreat. Using the large-scale pollen-based tree cover reconstruction by Schild et al. (2025), we evaluate the performance of the MPI-ESM Earth System Model in simulating tree cover dynamics from the Last Glacial Maximum to the present. Although the model reproduces the broad increase in tree cover during deglaciation and decrease throughout the Holocene, it fails to simulate the mid-Holocene maximum observed in the reconstructions. The model does capture the shift from energy-limited conditions during deglaciation to water-limited conditions in the early to mid-Holocene, and then back to energy-limited conditions in the late Holocene, but regional discrepancies remain substantial. MPI-ESM simulates too much forest in sparsely forested areas and too little forest in densely forested areas, particularly in mid- and high-latitude regions. Statistical analyses indicate that summer temperature dominates simulated high-latitude forests, while precipitation is critical in most other regions, contrasting with reconstructions that highlight cold-season temperature in temperate and boreal forests. Areas of model-data agreement show largely linear responses to climate drivers, whereas regions of disagreement exhibit non-linear dynamics to the temperature of the warmest month and over-sensitivity of the plant-physiological CO₂ response. Employing an emulator with a bias-corrected climate reduces the mismatch in the forest steppe transition zones, but does not lead to an overall improvement of the model-data agreement. In particular, the mismatch in the boreal region remains unresolved, suggesting structural limitations in the model. Improving dynamic vegetation models for simulating climatic transitions in both, past and future contexts, requires integrating realistic soil and permafrost processes, dynamic biome thresholds and disturbance regimes. Trait-based approaches could lead to better representation of the vegetation response to climate changes.

1 Introduction

The climate of the last 20 000 years is characterized by a large transition from the Last Glacial Maximum (LGM, ∼ 21 thousand years ago or ka) through the Holocene onset (∼ 11.7 ka) to the present, accompanied with major shifts in atmospheric boundary conditions. Continuous changes in the Earth's orbital configuration led to a spatial and seasonal redistribution of incoming solar radiation (Berger, 1978), resulting in a minimum in Northern Hemisphere summer insolation around 24 ka and a maximum near 11 ka. These insolation changes were amplified by Earth system feedbacks, including a rise in atmospheric CO₂ concentrations from ∼ 185 to ∼ 280 ppm (Köhler et al., 2017), the complete retreat of the North American and European ice sheets (Batchelor et al., 2019), and large-scale modulations of land, ice, and ocean surfaces (Ivanovic et al., 2016). In response, global temperatures have risen by ∼ 4–7 K since the LGM (Annan et al., 2022; Osman et al., 2021; Shi et al., 2023; Tierney et al., 2020), with most of the warming occurring during the deglaciation. This long-term trend includes spatially heterogeneous responses, with amplified warming at high latitudes due to polar feedbacks and a more gradual increase in the tropics (Shakun et al., 2012; Tierney et al., 2020). Precipitation increased significantly since the LGM (Bartlein et al., 2011; Shi et al., 2023), but regional responses were highly variable, with widespread Northern Hemisphere summer monsoon intensification and shifts in storm tracks (Braconnot et al., 2012; Clark et al., 2012; McGee, 2020). The last 20 000 years represent a key opportunity to study the Earth system's response to large-scale warming, as the natural radiative forcing changes were of similar magnitude to those projected under high emission future anthropogenic scenarios. Therefore, this period serves as a valuable testbed for evaluating the performance of Earth System Models and improving future climate and environmental projections, although the amplitude and patterns of forcing differ between past and future.

Among the most significant consequences of the large-scale climate changes since the LGM was the reorganization of global vegetation cover. However, the dynamics of this transformation – particularly the response of vegetation to radiative forcing and the identification of the dominant climatic drivers – remain only partly understood (Adam et al., 2021; Allen et al., 2020; Dallmeyer et al., 2022; Mottl et al., 2021; Nolan et al., 2018). Pollen-based reconstructions demonstrate a large-scale rearrangement of biomes from the LGM to the present (Li et al., 2025a). During the Last Glacial Maximum, forests were largely restricted to refugial areas, while open vegetation such as tundra and steppe dominated much of the Northern Hemisphere (Binney et al., 2017; Prentice et al., 2000). Tropical forests were often reduced and fragmented (Anhuf et al., 2006; Dupont et al., 2000; Sato et al., 2021; Wurster et al., 2010).

Following the glacial retreat, rising temperatures and retreating ice sheets enabled a gradual northward migration of plant species. This transition was neither uniform nor synchronous across regions; instead, vegetation dynamics were shaped by regional climate, topography, and species-specific dispersal capacities. Boreal forests began to establish in mid-latitude regions, followed by temperate deciduous forests during the early Holocene (Bigelow et al., 2003; MacDonald et al., 2000; Williams et al., 2004). In the tropics, forest recovery was spatially variable, shaped by regional hydroclimate (Wang et al., 2017).

The Holocene epoch, starting around 11 700 years ago, was marked by relatively stable and warmer climatic conditions, leading to the expansion and stabilization of forest ecosystems. According to pollen-based reconstructions, forest extent peaked during the mid-Holocene across large parts of the Northern Hemisphere (Cheng et al., 2021; Dawson et al., 2025; Githumbi et al., 2022; Li et al., 2023a; Sweeney et al., 2025).

During the late Holocene, human activities increasingly influenced vegetation dynamics, leading to widespread deforestation and land-use change – particularly after ∼ 4 ka (Ellis et al., 2013; Kaplan et al., 2011; Stephens et al., 2019). As a result, recent vegetation patterns reflect both natural ecosystem processes and cumulative anthropogenic impacts (Gaillard et al., 2010; Marquer et al., 2017).

So far, quantitative reconstructions of past vegetation were largely restricted to regional scales and often focused on the Holocene epoch (Cao et al., 2019a; Dawson et al., 2025; Githumbi et al., 2022; Li et al., 2023a; Trondman et al., 2015).

However, their comparability across regions has been limited by differences in methodology and spatial coverage. Recently, Schild et al. (2025) have provided the first REVEALS-based vegetation reconstruction for the entire Northern Hemisphere, synthesizing pollen records from ∼ 2800 sites across Eurasia, and North America for the last 14 000 years. This comprehensive dataset opens up new possibilities for model–data comparisons at continental to hemispheric scales, providing a valuable benchmark for evaluating the performance of Earth system models (ESMs) in simulating Holocene vegetation patterns.

In this study, we extend these reconstructions to the last 19 500 years and compare the dataset with a transient climate and vegetation simulation performed with the Max-Planck-Institute Earth-System-Model MPI-ESM1.2. We evaluate the agreement between simulated and reconstructed tree cover patterns and identify the main drivers of simulated vegetation dynamics. By analysing regional climate–vegetation relationships, we assess differences in areas of concordance and divergence between models and data to better understand model-data mismatches. Based on these findings, we outline improvements to dynamic vegetation models to better simulate biogeographical changes, especially under extreme climate conditions or in regions with sparse data.

2 Methods

2.1 MPI-ESM simulation

The model and simulation is explained in detail in Kleinen et al. (2023a) and Dallmeyer et al. (2022). The transient experiment was performed using the Max-Planck-Institute Earth System Model MPI-ESM, version 1.2 (Mauritsen et al., 2019) at T31GR30 spatial resolution, including the JSBACH land surface model and its dynamic vegetation module (Reick et al., 2013, 2021). Vegetation is represented as natural plant functional type (PFT) cover fraction per grid cell. JSBACH distinguishes several tree, shrub, and grass PFTs whose distributions are limited by bioclimatic temperature thresholds that reflect the physiological constraints and bioclimatic tolerance (e.g. chilling or heat requirements and cold resistance) of the plants. Changes in PFT cover are regulated by NPP-based competition, reflecting differences in e.g. moisture availability and physiological requirements. JSBACH also determines the fraction of bare ground, in contrast to the reconstructions in which the vegetation always sums up to 100 %.

The JSBACH version used lacks permafrost, dynamic soil development, and seed-dispersal modules (Dallmeyer et al., 2022). Soils therefore remain fixed, and all PFTs are assumed to have unlimited seed availability. Consequently, JSBACH, like most vegetation modules in Earth System Models, represent potential vegetation in quasi-equilibrium with the climate.

In this model, a transient simulation from 22 to 0 ka has been conducted (Kleinen et al., 2023a), prescribing orbital forcing from Berger (1978), greenhouse gas forcings from Köhler et al. (2017), and ice sheet extent from GLAC-1D (Tarasov et al., 2012). To better match climate reconstructions, meltwater from the Laurentide ice sheet was stored between 15.2 and 12.8 ka to mimic proglacial lakes (e.g., Lake Agassiz). A subsequent 1200-year release into the Mackenzie basin triggered an AMOC collapse and North Atlantic cooling consistent with the occurrence of the Younger Dryas, followed by rapid recovery after 11.8 ka. The simulated temperature trend from the Bølling-Allerød to the Early Holocene aligns well with reconstructions (Dallmeyer et al., 2022).

In line with the reconstructions, we bin the simulation output in 500-year intervals. This means, the time step 19 500 BP contains the simulated vegetation for the years 20 000 to 19 500 BP. For comparing MPI-ESM with reconstructions, we only use grid cells and time periods for which reconstructions are available.

Such temporal aggregation inevitably smooths out short-term variability and may dampen rapid vegetation changes, particularly during periods of abrupt postglacial transitions. As a consequence, the magnitude of submillennial variability, correlations during abrupt centennial-scale climate changes, and the detectability of sub-millennial non-linear responses may be somewhat reduced. Therefore, the comparison focusses on differences in multi-millennial trends and sensitivities, instead of submillennial vegetation changes.

2.2 REVEALS-based reconstructions

We use the REVEALS-based quantitative tree cover reconstruction by Schild et al. (2025) that contains 2752 records, distributed over the Northern Hemisphere. This dataset is based on the LegacyPollen 2.0 synthesis (Li et al., 2025b), that was compiled from individual publications and the Neotoma Paleoecology Database, which includes data from the European Pollen Database and the North American Pollen Database (Fyfe et al., 2009; Giesecke et al., 2014; Williams et al., 2018) and underwent standardized age modelling and taxonomic harmonization. Vegetation cover was reconstructed using the REVEALS model (Sugita, 2007), which corrects pollen counts for taxon-specific productivity and dispersal characteristics. The implementation followed the REVEALSinR routine (Theuerkauf et al., 2016), which incorporates repeated simulations with randomized input parameters to estimate uncertainty. The resulting dataset provides time-series of vegetation cover per taxon, including associated standard deviations and percentile-based confidence intervals. In general, REVEALS-based forest cover is characterized by a slight overestimation. This arises partly from persistent taxon-specific biases and partly from the presence of unvegetated areas. REVEALS-based vegetation cover reflects the composition of vegetated areas but does not capture the overall extent of vegetation. As a result, forest cover can be overestimated when large unvegetated areas occur within the pollen source area. Further details are described in Schild et al. (2025).

In this study, we extend the dataset by Schild et al. (2025) to the last 19 500 years and map them to the T31 grid used in the MPI-ESM simulation. Our analysis is based on all available reconstructions (see Appendix A). The distribution of the reconstructions in space and time is seen in Fig. 1.

Figure 1 (a) Temporal coverage of the REVEALS reconstructions [% of 41 time-steps between 19.5 and 0 ka] and (b) number of available grid cells per continent for each time-step.

2.3 Evaluation of the modern tree cover with remote sensing data

The present-day pattern of simulated and reconstructed tree cover fractions is compared to remote sensing data (MODIS, year 2000–2019; DiMiceli et al., 2022) (Fig. 2). Tree cover in Europe and China is relatively well-represented by REVEALS, whereas REVEALS strongly overestimates tree cover in Alaska, Northeastern North America, and Eastern Siberia, often exceeding 20 %. These discrepancies between MODIS and REVEALS tree cover may provide a first-order indication of the magnitude by which REVEALS could overestimate tree cover due to the reconstruction of vegetation composition only. However, this potential bias is unlikely to be temporally constant, as it depends on changing environmental conditions, vegetation dynamics, and the extent of non-vegetated surfaces.

Figure 2 Comparison of the (a) simulated and (b) reconstructed tree cover fraction of the modern time-slice with MODIS tree cover (DiMiceli et al., 2022). For MPI-ESM only grid cells for which reconstructions exist are shown.

The MPI-ESM model shows too little tree cover in Alaska and Western North America. However, it significantly overestimates tree cover along the East Asian Summer Monsoon (EASM) margin and in Central and Eastern North America, as well as in Europe. It is important to note that land-use is not included in this version of MPI-ESM, and many of these regions are heavily influenced by historical and contemporary land-use practices, which contributes to overestimated tree cover. Furthermore, validation studies indicate that the MODIS tree cover exhibits substantial uncertainties in boreal and transition regions, particularly at low tree cover fractions (< 20 %) (Montesano et al., 2009), further challenging the comparison with the simulation and reconstructions.

2.4 Comparison metrics

We use three statistical measures for evaluating the similarity of the REVEALS data and the MPI-ESM simulation, the Pearson correlation coefficient, the mean absolute error and the variance difference. The Pearson correlation coefficient (r) is used to assess the similarity in the temporal evolution of the tree cover in each grid cell. It is invariant to differences in the mean and variance. As such, it is well-suited for assessing if the REVEALS data and the MPI-ESM tree cover simulation follow a similar trajectory over time. We consider all grid-cells where the correlation between reconstructions and model data was significant at p < 0.1 in order to maintain spatial representativeness, particularly in regions with sparse data availability and low tree cover where correlations tend to be less statistically significant due to lower signal-to-noise ratios.

In contrast, the mean absolute error (MAE) quantifies the average absolute deviation between the two datasets. As a scale-dependent metric, it is highly sensitive to differences in the temporal mean and amplitude of changes. Consequently, MAE provides a direct measure of absolute agreement and reflects how closely the tree cover values simulated by MPI-ESM aligns in magnitude with the REVEALS-based estimates.

Third, the temporal variance captures the degree of fluctuation or variability within a signal. It complements the other measures by highlighting whether MPI-ESM captures not only the correct magnitude and trend but also the extent of variability inherent in the REVEALS data. We here consider the difference between the variance in the simulated and reconstructed tree cover.

Table 1 Different metric groups used in this study. For mean absolute error and correlation, the groups are based on the median (x, here x = 0.45), for the variance difference, the groups are based on a treshold (y) defined ad hoc as 10 % of the overall range (maximum minus minimum) of the variance difference.

For clarity of presentation, we assign the grid cells for each metric into different groups. For MAE and the correlation coefficient, we categorize the group by thresholds that are multiplier of the median (x) of all values. The variance difference values are categorized into three groups based on a threshold (y) defined as 10 % of the overall range (maximum minus minimum) of the variable. Grid cells with variance difference below this threshold indicate an underestimation of the variance in MPI-ESM compared to REVEALS, those exceeding the positive threshold are indicating an overestimation of the variance compared to REVEALS, and values within the interval ± threshold are defined to be of similar variance. The group definitions are provided in Table 1. Please note that grid cells with non-significant r-values are assigned to the group “no-correlation”.

2.5 Emulation with MPI-ESM simulated climate

To disentangle the effects of precipitation, soil moisture, temperature, and CO₂ on the simulated tree cover, we additionally train a statistical emulator to the simulation output of MPI-ESM, following the procedure of Adam et al. (2021). The emulator is based on generalized additive models (GAMs; Wood, 2017), which are non-parametric, non-linear statistical models, in which tree cover is the response variable and mean temperature of the coldest month (Tc), mean temperature of the warmest month (Tw), annual mean precipitation (P), annual soil moisture (Sw), and atmospheric CO₂-concentration are used as predictors.

The selection of these variables is motivated by their known importance for vegetation productivity and the implementation of bioclimatic thresholds for Tw and Tc in JSBACH. We train one GAM, which we call emulator in the following, using the data from all grid points, time steps, and climate predictors. Thereby, we estimate the modeled climate-CO₂-vegetation relationship across the full climate, CO₂, and vegetation gradients, independent of any specific location. This approach is justified by JSBACH modeling the vegetation independently for each grid box, without considering horizontal vegetation interactions such as seed dispersal. We use the function gam in the R package mgcv to fit the emulator (Wood, 2017). Following (Adam et al., 2021), we use a binomial response function due to tree cover being restricted to the interval 0 % to 100 %, and a full tensor product smooth with a cubic regression spline basis for modeling the vegetation-CO₂-climate relationship (Wood, 2006). Thereby, we estimate not just the non-linear effects of the individual predictors but also interactions between them. The number of smooth terms was determined with the default optimization routine in mgcv (Wood, 2017). To keep the fitting procedure computationally feasible, we use only every second time step for training.

We simulate sets of surrogate tree cover time series using transient climate forcing by all five variables, and by each individual variable alone, to disentangle the main climatic driver. For the latter, we employ the full emulator but set all predictors to time mean value in the respective grid box. The emulator explain the major vegetation changes over the deglaciation as indicated by predominantly high grid box correlations between the simulated and emulated tree cover (Fig. 3). Lower correlations are mainly concentrated in subtropical South America, the Sahara, arid regions in Asia and the Mediterranean region, and thus include mostly regions with marginal (Sahara) or low changes in tree cover. Marginal tree cover change may lead to spurious correlations, while low correlations in regions with low tree cover changes points to unresolved climate-vegetation dynamics in these regions.

Figure 3 Pearson correlation of the emulated tree cover performed with all forcings and the simulated tree cover.

In a second set of emulations, we employ the emulator trained on the MPI-ESM climate and force it with a bias-corrected climate, assuming that grid cell-specific model biases (e.g., related to orography) remain constant over time (Maraun and Widmann, 2018). To construct the bias-corrected climate forcing, we apply the delta method following Beyer et al. (2020), using additive corrections for temperature and multiplicative corrections for precipitation, both downscaled to a 0.5°-grid to align with the CRU-TS 4.08 baseline (University of East Anglia Climatic Research Unit et al., 2024; Harris et al., 2020). Soil moisture was only interpolated to the 0.5°-grid due to the lack of an equivalent observational dataset in CRU-TS 4.08. More details on the preparation of the climate forcing data, time-series and maps for some time-slices are given in the Appendix B.

2.6 Assessing the climate drivers

Complementing the emulation approach, which mimic the response of trees in the MPI-ESM model to climate changes, we assess the non-linearity of climate–vegetation relationships by distinguishing linear from non-linear behaviour. Linear relationships are quantified using Pearson correlation coefficients between simulated or reconstructed vegetation and the individual climate variables. Non-linear relationships are explored using generalized additive models (GAMs; Wood, 2017), with simulated tree cover as the response and the full set of climate variables as predictors. In contrast to the emulations, GAMs are fitted separately for each grid cell here, using 500-year binned datasets smoothed with a Loess filter. Each predictor is represented by a smooth, flexible function capable of capturing curves and other complex patterns. The effective degrees of freedom (edf) associated with each smooth term indicate the degree of non-linearity, with edf ≈ 1 reflecting near-linear relationships and higher values indicating stronger non-linear effects. F-values quantify the importance of each term in explaining variation in the response. For our analysis, we consider only grid cells with edf > 1.2 and p < 0.05, so that we can use the corresponding F-values as indicators of the non-linear climate–vegetation relationships. The most influential variable is identified as the predictor with the highest F-value.

3 Results and Discussion

We structure our analysis around a sequence of interconnected questions that progressively link model evaluation with process understanding. We begin by comparing simulated and reconstructed tree cover dynamics from the Last Glacial Maximum to the present (Sect. 3.1), establishing the overall level of agreement between MPI-ESM and REVEALS. Building on this, we quantify and characterize the spatio-temporal differences between both datasets (Sect. 3.2), thereby identifying where and when discrepancies occur.

In a next step, we investigate the underlying mechanisms by analysing the main drivers of tree cover change (Sect. 3.3). This process-based perspective is further refined by explicitly contrasting regions with high and low model–data agreement to assess which factors constrain tree cover dynamics in the model (Sect. 3.4). Finally, we address the question whether a bias-corrected climate can reduce the identified model–data differences (Sect. 3.5), thereby assessing the influence of systematic climate biases on simulated tree cover dynamics.

3.1 Tree cover dynamics in REVEALS and MPI-ESM from the LGM to Present

For the Last Glacial Maximum (LGM), only ∼ 50 vegetation records are available (Fig. 4). REVEALS indicates a predominantly open landscape, consistent with previous reconstructions (Cao et al., 2019a; Davis et al., 2024; Kern et al., 2025; Li et al., 2025a; Prentice et al., 2000). Dense forests occurred mainly in eastern North America (up to 75 % tree cover) and eastern Asia (up to 100 % in Japan), with moderate cover in western North America (up to 60 %). Isolated cells with moderate or high tree cover appear in western central Europe (up to 35 %), the Mediterranean (up to 20 %), Alaska (up to 50 %), and eastern Siberia (up to 20 %). MPI-ESM broadly agrees but shows higher tree cover in Europe and temperate North America, and lower cover in most other regions, including Siberia and the rim of the Tibetan Plateau. The isolated high-cover cells in REVEALS are not reproduced by MPI-ESM, likely reflecting glacial refugia not captured at the model's spatial scale (Leroy et al., 2020).

Figure 4 Tree cover fraction in MPI-ESM (upper row), REVEALS (second row) and tree cover difference between MPI-ESM and REVEALS (last row) for selected time-slices. A video showing all time slices can be found here (Dallmeyer, 2026).

By 10.5 ka, regions of reconstructed high tree cover expanded markedly. Dense forests developed across eastern North America (70 %–80 %) and western Canada (Pacific coast), and most Siberian records show 60 %–80 % tree cover. Southern China changed little from the LGM, while the northern part of China has experienced a moderate increase in tree cover (30 %–50 %) and up to 70 % in the north eastern part along the Pacific coast. In Europe, REVEALS indicates only modest increases to 30 %–60 %. MPI-ESM captures the general trends but deviates regionally, underestimating forest cover along the North American Pacific coast, in eastern North America's temperate zone, the Mediterranean, and parts of Siberia, while overestimating it in central North America, western and central Europe, northern Siberia, and the East Asian monsoon region.

Towards 7 ka, reconstructions indicate a further increase in tree cover fractions in the Northern Hemisphere north of approximately 50° N, particularly across boreal latitudes and Europe, similar to previous studies (Sweeney et al., 2025). In contrast, MPI-ESM simulates a small but widespread decline, since the maximum Northern Hemispheric tree cover is already reached during early Holocene (Dallmeyer et al., 2022). Exceptions to this trend include periglacial areas in North America and Scandinavia, where tree coverage remains relatively stable or increases slightly. In boreal regions of North America and parts of Siberia, MPI-ESM simulates much less tree cover (∼ 30 % to 80 %) than suggested by the reconstructions. The pattern of disagreement is similar as in 10.5 ka, but much more pronounced in 7 ka. Only in some parts of Europe and eastern North America, tree cover is consistent with the reconstructions, showing dense forests.

At 0 ka, the absence of land-use effects in MPI-ESM likely leads to overestimated tree cover, especially in Europe. Overall, the spatial mismatch between REVEALS and MPI-ESM remains similar to 7 ka: the model simulates too little forest in boreal Asia and North America, but too much in East Asia and central North America.

Overall, MPI-ESM captures the continental-scale patterns and temporal evolution of tree cover reasonably well, but its performance decreases substantially at regional scales.

3.2 Spatio-temporal differences in the tree cover dynamics between MPI-ESM and REVEALS

To identify regional differences between the REVEALS-based reconstructions and the MPI-ESM simulation, we calculate three different metrics (Fig. 5): (a) the Pearson correlation coefficient measuring the synchronicity of the time-series, (b) the MAE measuring the absolute difference in the tree cover estimates and (c) the tree cover variance difference as estimate of differences in the amplitude of tree cover changes. For easier accessibility, we group these measures into statistical categories, aligned on the median value over all grid cells (MAE and r) or the range of the variance difference values (see methods, Table 1).

Figure 5 Results of different metrics for determining the differences between the simulated and reconstructed tree cover time-series per grid-cells. These are (a) Pearson correlation coefficient, (b) mean absolute error (MAE), and (c) differences in the tree cover variance (Model − Reconstructions). For easier comparability, the results have been grouped in statistical categories by the median value (see Table 1). In panel (b) the regions with highest MAE are marked with rectangles Higher correlation and lower MAE indicate better agreement between simulation and reconstruction, while larger variance differences highlight regions where the model fails to capture temporal variability (blue = model overestimates variance, red = model underestimates variance).

Overall, most of Europe demonstrates a highly positive correlation between REVEALS and MPI-ESM, indicating that the model generally captures the European climate and vegetation dynamics. The Mediterranean region stands out as an exception, experiencing a low or even highly negative correlation. The MAE is moderate for most of Europe, but the model has a higher variance than the reconstructions. These differences could be enhanced by the ignorance of land-use change in Europe which leads to systematic differences in the late Holocene tree cover. However, previous studies report no substantial large-scale human impact on vegetation composition before 2 ka (Roberts et al., 2018; Trondman et al., 2015), thus land use cannot be the only driver of model-data differences. Estimates of potential tree cover for the present time slice indicate more than 80 % tree cover in most parts of Europe except for Spain (Roebroek et al., 2025). This is in line with MPI-ESM.

In North America, mid-latitude regions (40–60°) display a strongly positive correlation, particularly in the eastern part, where also the MAE is low despite a mostly lower variance in the model. However, Central North America exhibits notable weak or even negative correlations and exceptionally high MAE. As the variance difference is low in most grid cells, the high MAE points to a systematic bias and/or diverging trends of the tree coverage in this region. Most grid cells in Western Canada show high MAEs with differing correlations. Coastal regions, also in Alaska, experience negative or weak correlations, while in inland Western Canada correlations are positive in some grid cells only.

In Asia, the correlation indicates no clear pattern, reflecting no robust alignment between REVEALS and MPI-ESM. Siberia displays rather weak or negative correlations, while positive correlations can be observed in non-forested regions, in Western China and along the Pacific coast. Positive correlations are also observed in certain grid cells in Western Siberia. The MAE is high in most grid cells in East Asia, particularly in Central and Eastern Siberia and along the northern rim of the East Asian monsoon region. The variance difference is low in most grid cells, which eliminates amplitude differences as reasons for high MAEs.

These regional discrepancies between the model and the reconstructions point to several key challenges in comparing the REVEALS-based reconstructions with vegetation models. Possible reasons include methodological caveats such as shortcomings in the process of reconstructions, including the REVEALS model (Githumbi et al., 2022; Kaplan et al., 2017; Li et al., 2023b; Schild et al., 2025), biases in the simulated climate, missing or idealized model components (Dallmeyer et al., 2022, 2023; Li et al., 2025a; Zhang et al., 2025) and differences in the vegetation-climate relationship. These factors may lead to local ecological processes, migration lags, and disturbance regimes affecting the reconstructions, while being inadequately represented in the model. This may be particularly problematic in ecotones, where small climatic changes or disturbances can cause large shifts in vegetation composition because multiple vegetation types coexist near their ecological limits.

One important source of divergence lies in the fundamental assumptions underlying each approach: REVEALS estimates plant cover fractions that always sum to 100 % vegetation cover, as it is based on pollen assemblages that represent only vegetated areas without accounting for the proportion of bare ground. Consequently, REVEALS reflects the relative composition of vegetation types within the vegetated area, rather than the absolute area (f_i) covered by each vegetation type (i). In contrast, MPI-ESM explicitly calculates both, a bare soil fraction and the composition of plant functional types (c_i) on the vegetated area (V), with $c_i=f_i/V$. While in principle this allows a comparison to REVEALS cover fractions, we have not applied it in this study, because this approach is only valid for regions and time periods where overall tree and vegetation cover is high, such as during the Holocene in Europe (Dallmeyer et al., 2023). In rather open landscapes, this comparison suffers from the artefact of high relative PFT fractions (c_i), because for numerical reasons, PFT cover fractions always have to be non-zero in the model. We therefore use the absolute area covered by the PFTs (f_i). This fundamental difference (vegetated area vs. vegetation ratios) influences the MAE.

We therefore analyse the differences in tree cover estimates (MPI − REVEALS) grouped by REVEALS tree cover fraction classes. The results demonstrate a clear, non-linear trend across the vegetation gradient (Fig. 6): In areas/time-slices for which REVEALS indicates low tree cover (0 %–20 %), MPI tends to overestimate vegetation cover, with a mean positive bias of approximately 14 %, turning into a negative bias for high REVEALS tree cover classes (−32 % for the 80 %–100 % class). Part of this pattern, however, is methodological, because the difference metric (MPI − REVEALS) is analysed as a function of REVEALS tree cover. This mathematical dependence tends to produce increasingly negative values at high REVEALS fractions even in the absence of ecological effects. The figure should therefore primarily be interpreted as illustrating the magnitude of model–data deviations along the vegetation gradient rather than a purely ecological non-linear response.

Figure 6 Difference between simulated tree cover (MPI-ESM) and reconstructed tree cover (REVEALS) across REVEALS tree cover classes. The boxplots show the distribution of tree cover differences (MP − REVEALS) at any time in relation to the reconstructed REVEALS tree cover class (0 %–20 %, 20 %–40 %, …, 80 %–100 %). The horizontal line in each box indicates the median difference, while the black dot represents the mean. Positive values indicate that the simulation overestimates tree cover compared to REVEALS, while negative values indicate underestimation.

This pattern is nevertheless consistent with the increase in MAE from the LGM (open landscapes) to the Holocene (forested period) (Fig. A3 in the Appendix A) and underlines that the comparison of tree cover fraction of grid cells in the model against tree cover fractions of total vegetation in the reconstructions increase the MAE, particularly in highly forested landscapes.

The positive bias in sparsely wooded regions may result from the model's omission of environmental constraints such as wetlands, permafrost, or water-logged soils, while the negative bias in regions of dense forest may arise from limited number of PFTs, the global, relatively simple definitions of PFTs and their interaction, as well as tuning parameters calibrated for modern rather than past conditions (Dallmeyer et al., 2023; Horvath et al., 2021; Zhao et al., 2025). In addition, the coarse grid resolution of the MPI-ESM leads to a smoothing of regional vegetation patterns. Local maxima, such as tree refugia or areas of particularly dense forest due to local climate conditions, are not adequately captured, resulting in a more homogenized vegetation distribution that overlooks the fine-scale structural heterogeneity evident in the reconstructions.

The differences in absolute tree cover levels expressed by the MAE show a pattern that shares features with the modern deviations to the MODIS observations (Fig. 2, method part). This could point to systematic errors in the REVEALS data and simulation results. REVEALS rather seems to overestimate tree cover in boreal latitudes (Schild et al., 2025) and thus may be too sensitive in cold climates such as glacial periods with open landscapes. Some taxa can occur as trees or shrubs, and their differentiation is not possible based on pollen, especially as the pollen taxonomy was harmonized at the genus level (Herzschuh et al., 2022). More generally, each taxon is assigned to a single plant functional type (PFT), which inherently departs from ecological reality, as many taxa comprise species belonging to different PFTs. Although the reconstruction method by Schild et al. (2025) applies a static threshold that assigns shrub status to common tree taxa (e.g. Betula, Alnus, Corylus) north of 60° N, this adjustment is unlikely to eliminate all biases in tree and shrub PFT representation at northern latitudes. Consequently, the distinction between shrubs and trees in REVEALS may remain imprecise, potentially leading to overestimated tree cover in boreal regions.

On the other hand, the MPI-ESM simulation might suffer from the fixed bioclimatic limits and a tendency towards too cold conditions in the boreal latitudes. For Central North America and the East Asian monsoon margin, the spatial pattern of the MAE diverges from the pattern of disagreement with the MODIS data. In terms of MAE, these regions, which are both forest–steppe transition zones, thus represent the most challenging areas in the comparison.

Many methodological limitations in both the model and the REVEALS reconstructions cannot be quantified directly, and their individual contributions to the model–data tree cover discrepancies therefore remain unresolved. However, we can assess the extent to which systematic climatic biases or differences in the climate–vegetation relationships contribute to the observed mismatches in tree cover patterns. To this end, we analyze the regional climatic drivers of the simulated vegetation dynamics and evaluate which factors distinguish grid cells with close model–data agreement from those exhibiting substantial deviations. Additionally, we employ the trained emulator with bias-corrected climate inputs to further isolate the influence of climatic biases.

3.3 What are the main driver of the tree cover change?

To evaluate the climate–vegetation relationships, we apply a statistical emulator (see Methods) that captures both linear and non-linear responses of vegetation to climatic forcing. To maximize the training data we use the complete global model domain, although reconstructions are limited to the Northern Hemisphere. Pearson correlation coefficients and generalized additive models (GAMs) are additionally used to disentangle the linearity and non-linearity of the climate–vegetation relationship (see Methods). The tested predictors include temperature of the warmest and coldest month (Tw, Tc), annual precipitation (P), annual mean soil moisture (Sw), and atmospheric CO₂ (Fig. B1) and we determine the locally most influential predictor as the emulated tree cover curve with the highest correlation to the simulated/reconstructed tree cover. We emphasize that the emulator diagnoses sensitivities intrinsic to the model itself, and thus characterizes model behavior rather than real-world sensitivities.

The emulator results show coherent regional patterns. In high-latitude and high elevation regions (Siberia, Alaska, northern Canada, Tibetan Plateau), the emulation forced with Tw (orange) shows the strongest correlation with simulated tree cover (Fig. 7), consistent with the well-established control of growing-season warmth on vegetation in cold environments (Adam et al., 2021; Dallmeyer et al., 2021; Li et al., 2022a; Seddon et al., 2016; Tucker et al., 2001). In the tropics and subtropics (Africa, South America, Southeast Asia), tree cover is predominantly governed by precipitation (P, light blue), reflecting the central role of water availability (Adam et al., 2021; Gosling et al., 2022; Zhao et al., 2018).

Figure 7 Spatial distribution of the most influential climatic variables determining vegetation dynamics across different regions of the world. Each grid cell is color-coded according to the dominant driver of vegetation change in the emulations, which is here based on the emulation with single variable forcing that show the highest Pearson correlation to the simulated tree cover. The variables considered are atmospheric CO₂ concentration (yellow), annual precipitation (P, light blue), temperature of the coldest month (Tc, dark blue), temperature of the warmest month (Tw, orange) and annual mean soil moisture (Sw, purple). Grid cells marked with black dots indicate regions where the difference in correlation coefficients between the leading predictor and the second most influential predictor exceeds 0.1 (the median difference), highlighting areas with a comparatively strong dominance of a single driver.

Temperate regions display mixed controls. In Central North America, Tw dominates, whereas in Europe both P and Tc are influential. In large parts of northwestern Asia, precipitation emerges as the primary driver. CO₂ (yellow) dominance is spatially scattered, with occurrences in South America, East Africa, Australia, and southern Asia, particularly in the East Asian monsoon region. This strong impact of CO₂ has also been found in time-slice-based biome change analysis for the deglaciation period (Tian et al., 2018). However, this relationship should be interpreted as a model-derived sensitivity rather than direct evidence of physiological CO₂ fertilisation. The apparent CO₂ control may partly reflect a proxy for hemispheric deglacial trends that are not fully captured by the other climate predictors. Soil moisture shows weak influence overall, except in arid regions (e.g. North East Africa) and in tropical monsoon margins.

The spatial patterns inferred from the emulations, correlations, and GAMs are largely consistent and remain robust when the analysis is repeated with a data subset matched to the spatial and temporal structure of the reconstructions. In addition to the dependency of tree cover on Tw in the boreal latitudes, GAMs indicate a non-linear influence of Tc in temperate regions (Central North America and Europe) and around the Tibetan Plateau (Fig. 8c), consistent with physiological evidence that winter cold constitutes a major limitation for tree survival and recruitment. Extremely low temperatures can cause cellular freezing, xylem embolism, and bud mortality, while repeated freeze–thaw cycles and deep soil frost disproportionately damage roots and seedlings (Charrier et al., 2017; Sakai and Larcher, 1987). Juvenile trees, with shallow roots and limited carbohydrate reserves, are more vulnerable than grasses or shrubs, favoring cold-tolerant vegetation. Many temperate species require a chilling period to break dormancy, but excessively low temperatures can exceed frost tolerance and prevent establishment, creating a narrow climatic window for regeneration. In JSBACH, these constraints are represented by PFT-specific minimum-temperature thresholds for the coldest month, below which tree establishment is suppressed, contributing to the non-linear response of trees to Tc.

Figure 8 Spatial distribution of the most influential climatic variables determining vegetation dynamics according to different methods (Emulation, linear Pearson correlation and GAM) for all considered grid cells in MPI-ESM and REVEALS. The key variables include atmospheric CO₂ concentration (yellow), annual precipitation (P, light blue), temperature of the coldest month (Tc, dark blue), and temperature of the warmest month (Tw, orange) and annual mean soil moisture (Sw, purple). Please note that “most influential” means here the variable with the highest R (Pearson correlation) or F value (GAM) or the highest correlation between the emulated and simulated/reconstructed tree cover.

In contrast to the spatially coherent patterns in MPI-ESM, the relationships between reconstructed REVEALS tree cover and emulated tree cover show a much more fragmented structure of climate drivers (Fig. 8d). Moreover, the dominant climatic drivers inferred from REVEALS and the simulated climate differ markedly from those in the model. Most notably, the strong influence of Tw in boreal regions is absent in the reconstructions. Instead, REVEALS tree cover correlates more strongly with P and Tc, even at high latitudes. This finding contradicts the commonly assumed dominance of growing-season temperature in controlling boreal vegetation dynamics and indicate that MPI-ESM does not sufficiently capture the past real-world climate change or that the strong restriction of the extratropical trees by warm season temperature limitations are oversimplified and regionally inappropriate.

A substantial number of grid cells shows the highest correlation of REVEALS tree cover with CO₂. Although time-slice experiments have demonstrated a strong impact of the ecophysiological CO₂ effect on vegetation differences between the LGM and the present (Claussen et al., 2013; Woillez et al., 2011), the strong correlation between REVEALS-derived tree cover and the CO₂ signal in this study rather suggests that the local REVEALS reconstructions follow the overall deglacial warming trend – reflected in the CO₂ record – more closely than they follow the simulated local climate or climate-vegetation relationship. This statistical agreement may therefore primarily reflect the shared large-scale deglacial trend in both variables, rather than a direct plant-physiological CO₂ response in the real-world.

Because the emulator diagnoses model-internal sensitivities within the MPI-ESM framework, discrepancies between simulated and reconstructed patterns do not necessarily indicate model failure; they may also arise from characteristics of the reconstructions, including their temporal resolution, chronological uncertainties, and spatial aggregation, which can enhance alignment with long-term trends while smoothing local variability.

The dominant climatic drivers in MPI-ESM according to the emulation vary through time. While energy limitations (i.e. Tw or CO₂ as most influential variables) are a key driver during the deglaciation, moisture availability (i.e. precipitation and soil moisture dominant driver) becomes increasingly important toward the Holocene (Fig. 9). The moisture availability influence peaks at 7.5 ka, indicating a widespread shift from energy- to water-limited environments in the early Holocene before shifting back to higher energy-limited influences in the late Holocene. The reconstructions reflect the transitions between energy- and water-limited regimes (in terms of their correlation with the emulations), with a maximum area of water-limited regimes reached at 8.5 ka. The influence of Tc stays relatively constant through time in the simulations, whereas it increases by ∼ 10 % in the Holocene in the reconstructions.

Figure 9 Summary of the climate drivers (averaged over all grid cells) for different periods and their evolution over time (moving average over 15 time-steps = 7500 years) for (a) MPI-ESM and (b) the REVEALS-based reconstructions (vs. simulated climate). For easier interpretation we merged the climate drivers into the groups representing energy-limited (Tw and CO₂), moisture-limited (Sw and P) and winter-cold-limited (Tc) regimes.

3.4 What constrains tree cover dynamic in the model – High- vs. low-agreement grid cells?

Based on the correlation, mean absolute error (MAE), and variance differences between simulated and reconstructed tree cover, we classify all grid cells into statistical categories (cf. Fig. 5 and Table 1) and analyse differences in tree cover responses to climate change among these groups. Additionally, we assess whether the observed group differences arise from insufficient data (e.g., incomplete or short time series) or are related to geographic position (latitude or longitude), and whether they affect the representability of the tree cover dynamic by the emulation. To this end, we pre-select variables that show significant differences between the groups based on pairwise comparisons using the Wilcoxon rank-sum test (Appendix C, Table C1 in the Appendix C). Then, we use this preselection of predictors for running a PCA-analysis to determine the most influential variables for group-separation.

It should be noted that high MAE values and low correlations do not necessarily indicate deficiencies in the climate-vegetation relationship in the model. They may also reflect the fact that the reconstructions capture local ecological processes, migration lags, or microclimatic conditions that are not represented in the coarse-resolution simulation. Since these factors cannot be filtered out, all grid cells are included in the analysis to examine the dependence on climatic conditions.

Grid cells exhibiting a synchronous dynamic in simulated and reconstructed tree cover (positive correlation-group) differ markedly in their tree cover response to climate change from those with negative correlations (Fig. C1 in the Appendix C). High agreement occurs where both simulated and reconstructed tree cover follow the simulated climate linearly as indicated by strong Pearson correlations with CO₂, precipitation (P), soil moisture (Sw), and temperature of the coldest month (Tc). In these cells, Tc shows high correlations with precipitation (r ≈ 0.85), which itself closely follows CO₂ changes (r ≈ 0.7). The temperature of the warmest month, in contrast, is only moderately correlated with CO₂ and precipitation (r ≈ 0.45 and r ≈ 0.3). Good agreement can also arise in cells with more non-linear responses to precipitation, which the GAM analysis identifies as the main driver in significantly more cells of the positive-correlation-group than in the negative-corelation-group. In contrast, disagreement (negative correlation) is linked stronger to non-linear responses to Tw, which emerges significantly more often as the main driver in grid cells of the negative-correlation–group according to the GAM analysis. In grid cells of the negative-correlation-group, Tw is furthermore rather negatively correlated with CO₂. Neither the variance difference nor the MAE are deviating significantly between grid cells in the negative- and the positive-correlation-group. The agreement is furthermore independent of the time-series length, including the timing of the first data point and the number of missing values.

Figure 10 Principal Component Analysis (PCA) biplot of climate–tree cover relationships for key-climate variables in grid cells with (a) positive (blue) and negative (red) correlations, (b) low (blue) and high (red) MAE and (c) overestimated (blue), underestimated (red) and similar (grey) variance between simulated and reconstructed tree cover. Arrows indicate the direction and strength of climate-related variables contributing to the ordination. Hereby, “F_” corresponds to the F-coefficients in the GAM-analysis (i.e. non-linear relationship) and “R_ ” to the Pearson correlation coefficients (i.e. linear relationship) in the “model”, the reconstructions (“_rec”) or in the emulation (“_emu”). The considered variables are precipitation (P), soil moisture (Sw), temperature of the coldest (Tc) and warmest (Tw) month and CO₂. For instance, R_Tw_rec corresponds to the Pearson correlation coefficient for reconstructed tree cover and simulated Tw, while R_CO₂ corresponds to the Pearson correlation coefficient for the simulated tree cover and CO₂. Significance is tested via Wilcoxon test (a, b) or Kruskal-Wallis test (c).

The PCA on the preselected key variables (F-statistics and Pearson correlations) confirms the significant (Wilcoxon test, p = 6.75 × 10⁻¹¹) separation between the correlation-groups, with PC1 and PC2 jointly explaining 49.5 % of the variance (Fig. 10a). High agreement aligns with strong positive loadings of the Pearson correlation coefficients between tree cover and CO₂ and Tc, as well as the strong positive precipitation-CO₂- and precipitation-Tc-correlations, which most strongly differentiate the positive-correlation-group. In contrast, the non-linear response to Tw is clearly oriented towards the centroid of the negative-correlation-group (negative loadings on PC1), substantially contributing to its separation from the positive-correlation-group. Other variables (including the linear and non-linear relationship to precipitation (R_P, F_P), and correlation to soil moisture (R_Sw) load mainly on PC2 and play only a minor role in distinguishing the groups.

MAE groups do not match with any correlation-group classification. However, grid cells with low MAE differ in climate response from high-MAE cells typically located at mid-latitudes (∼ 50° N) (Fig. C2 in the Appendix C). In the mean, grid cells with early starting reconstructions have a significantly lower MAE than grid cells with later starting reconstructions. High-MAE cells exhibit significantly higher linear correlations of the simulated tree cover with Tc and CO₂. No significant differences in dominant climate drivers are found between the MAE groups, indicating that the prevalence of a particular variable as the main driver does not systematically vary with model error. The PCA of MAE groups (PC1 66.8 %, PC2 26.4 %) confirms the statistically significant (Wilcoxon test, p = 2.04 × 10⁻³) separation, mainly by the linear correlation to Tc, and CO₂.

For the tree cover variance difference, we compare all three groups. Grid cells with overestimated variance in the model are characterized by earlier, more complete time series, stronger linear correlation with CO₂, and Tc, and also higher correlations between precipitation and CO₂ and Tc compared to grid cells with similar variance (Fig. C3 in the Appendix C). Moreover, overestimated variance co-occurs with much higher linear correlation between Tw and CO₂ (r ≈ 0.6), but can be less well represented by the full-forcing-emulator.

Grid cells with underestimated variance in the model mainly cluster in Eastern North America and show higher linear correlations with Tw, Tc, P, and CO₂ and lower non-linear relationship to Sw, Tw and CO₂ compared to grid cells with similar variance. Consistently, correlation between P and CO₂, Tw and Tc are significantly higher compared to grid cells with similar variance. Interestingly, grid cells within both groups (overestimated and underestimated variance) exhibit higher correlations between simulated and reconstructed tree cover than grid cells with similar variance, despite larger MAEs. This indicates that the temporal structure can be preserved even when variability and amplitude is misrepresented. The lower correlations and low MAE in the group with similar variance points to poor temporal coherence despite a reasonable magnitude fit.

The PCA of the preselected variables confirms the significant separation between the variance groups (Kruskal-Wallis, p = 1.66 × 10⁻⁵) in their multivariate environmental characteristics, with PC1 and PC2 explaining 43.6 % of the variance. Along PC1, the most influential contributors to group separation are the contrasting non-linear vs. linear responses of the tree cover to CO₂. The non-linear relationship to CO₂ (F_CO2) loads positively and aligns with the overestimated-variance–group, while the linear correlation to CO₂ exhibits negative loadings associated with the underestimated-variance–group. Additional separation along the horizontal axis arises from variation in linear correlations with Tc, Tw, and P, which load negatively towards the underestimated-variance-group. Separation along PC2 is driven primarily by differences in the Tw–CO₂ relationship and by the length of the time-series.

Overall, these analyses show that the agreement between simulated and reconstructed tree cover is strongly influenced by the choice of evaluation metric and the specific ways in which vegetation responds to climate drivers. The temporal evolution is best captured when tree cover responds linearly to the hemispheric warming and moistening trends, reflecting the coherent signal shared by precipitation, cold-season temperature, and the global CO₂ trajectory. In these grid cells, reconstructed tree cover appears to follow the large-scale trend closely, probably indicating regions where vegetation dynamics are primarily controlled by broad-scale climatic forcing and are therefore well captured by the model.

In contrast, grid cells in which a non-linear response to summer temperature fosters poor temporal correlations could coincide with regions in which the reconstructions indicate strong regional or local vegetation trajectories that are not governed primarily by hemispheric climate forcing. This pattern suggests that the mismatches to REVEALS do not necessarily arise from shortcomings in the model's representation of large-scale vegetation–climate interaction, but rather highlights areas where regional ecological or climatic processes play a comparatively larger role, influencing the reconstructions.

MAE reflects differences in simulated tree cover sensitivity to cold-season temperature and CO₂, with greater disagreement when tree cover is more linearly aligned with Tc and CO₂. Deviations in variance do not increase model–data asynchronicity but tend to coincide with higher MAE, indicating that the timing of changes can be accurately reproduced even when amplitude is misrepresented. Strong linear correlations with CO₂ lead to overestimated variance, whereas non-linear relationships tend to produce underestimated variance. Overall, the response of the climate–vegetation system to CO₂ forcing emerges as one key controlling factor. Linear alignment of climate with the CO₂ signal co-occurs with synchronous tree-cover changes in reconstructions and MPI-ESM, but simultaneously increases MAE and exaggerates amplitudes in the model, pointing to a too strong sensitivity of forests to CO₂ variations in MPI-ESM. This high sensitivity to CO₂ in MPI-ESM has also been observed in other studies (Arora et al., 2020; Kleinen et al., 2023a; Rogers et al., 2017).

3.5 Can we reduce the biases by prescribing a bias-corrected climate?

The trained emulator is forced with a bias-corrected climate to test the impact of potential, systematic climate offsets on the tree cover change. The bias-corrected climate and examples of the resulting tree cover emulation in comparison with the original simulation is provided in Figs. B1–B3 in the Appendix B. To evaluate where the bias-corrected emulation enhances model-data agreement (i.e. it shows higher r or lower MAE values or variance differences), we statistically compare the distribution of grid cells with relative improvement (> 25 % percentile) or worsenings (< −25 % percentile) for the correlation-, MAE-, and variance-groups, using the same classification as in Fig. 5 and Table 1. Note that a relative improvement in model-data agreement does not necessarily imply a good absolute agreement, as r or MAE values may still indicate substantial discrepancies.

The bias correction enhances the model–data agreement at the regional scale, but not in the mean across all grid cells. Correlation coefficients and variance differences improve in only about one third of the grid cells, and the MAE improves in roughly 40 %, while around 25 % of the grid cells in all groups remain unaffected, respectively. The remaining cells show deteriorated performance. MAE improvements (blue areas, Fig. 11) concentrate in Western/Central North America Scandinavia, Southern Europe and Eastern Asia and mostly coincide with an increase in the tree cover correlation, except for Asia. Variance differences increase across much of North America, Europe, and China, but decrease north of 40° N in Asia and south of 40° N in North America, pointing to a high degree of regional sensitivity in the model's ability to reproduce reconstructed variability patterns. The bias correction is most effective in grid cells with initially poor performance, such as those exhibiting negative correlations, high MAE, or large variance mismatches. For correlation and variance, the improvements show no systematic association with the dominant climate driver, suggesting that errors in these metrics are broadly distributed and more difficult to correct with this kind of bias correction. In contrast, MAE improvements occur disproportionately more in grid cells where temperature (Tw) is the primary driver, while deteriorations are more frequent in precipitation-dominated (P) cells.

Figure 11 Impact of the bias correction on the three performance metrics (correlation, MAE, and variance). The left panels display spatial distributions of significant changes in these metrics after bias correction. Blue indicates relative improvement compared to MPI-ESM (> 25 % percentile), orange relative worsening (< −25 % percentile), and grey no change. The right panels summarize these effects by showing the proportion of grid cells with each type of change across the different model-performance categories (mae-, correlation- and variance-groups). Note that a relative improvement (blue) does not necessarily imply good absolute agreement; it only indicates better agreement in the emulation compared to the original MPI-ESM simulation.

The results indicate that MPI-ESM performs reasonably well at the hemispheric scale, so that regional bias corrections often come at the expense of other areas, producing localized improvements but limited overall gain. The corrections primarily reduce large systematic errors, leading to improvements in MAE and correlation, but they do not consistently preserve the amplitude of variability, and variance is degraded across large areas. Consequently, the bias-correction methods tend to redistribute errors rather than uniformly improving the model performance across the domain. Furthermore, tree-cover dynamics driven by temperature are more amenable to correction, whereas precipitation-driven dynamics remain more sensitive and are more likely to deteriorate under bias correction, reflecting the differing spatial and temporal characteristics of these climate drivers.

To better understand the regional behaviour, we compare the time series of tree cover fractions (Fig. 12) from reconstructions (red), the MPI-ESM simulation (black), and statistical emulations (blue) with bias correction for the six regions with the largest MAE. Overall, the (bias-corrected) emulated tree cover aligns better with the reconstructions than the simulation in most high-MAE regions. Particularly, in the temperate forest–steppe transition zones in Central North America and Northern China (EASM margin), the emulation produces substantially lower tree cover and aligns more with the reconstructions than the simulation, indicating that the large simulation–reconstruction discrepancies in these regions mainly stem from systematic climate biases.

Figure 12 Comparison of reconstructed (red), simulated (black), and emulated (blue) tree cover fraction over the past 19 500 years in six selected regions, i.e. the regions indicating the highest MAEs between the simulated and reconstructed tree cover. Shaded areas represent uncertainty, assessed via bootstrapping. The central panel shows the spatial distribution of the MAE between simulated and reconstructed tree cover (see Fig.  5), classified into four categories: very low, low, moderate, and high. Regions outlined with black boxes correspond to the locations of the time series plots. Please note that for Northeastern North America, no time-series are shown since the time-series are relatively short there.

In Western Canada, bias-corrected warmer summers advance postglacial tree-cover growth and reduce offsets with reconstructions by 30 %, but neither simulation nor emulation capture the mid-Holocene maximum. The model appears to miss the prolonged early- to mid-Holocene warming, pointing to deficiencies in the simulated climate trajectory

In Central Siberia, the bias correction increases emulated tree cover too strongly during deglaciation and leads to a too early peak, indicating an overestimated sensitivity to warm-month temperature and explaining the worsening of the correlation despite an improved MAE in this region. The simulation with much lower temperatures agrees with reconstructions until 9.5 ka, after which their trends diverge, while the emulation and reconstructions only converge in the late Holocene. Neither simulation nor emulation reproduce the reconstructed Holocene increase, pointing to a wrong driver or structural vegetation deficiencies in the model. This mismatch may partly be caused by the missing permafrost processes in this version of MPI-ESM. Permafrost is a major regulator of tree cover across northern high-latitude regions, including Siberia, western Canada, and Alaska. It controls soil thermal and hydrological conditions that limit tree establishment and growth (Cao et al., 2019b; Helbig et al., 2016; Jin et al., 2022, 2021; Kruse et al., 2022; Stuenzi et al., 2021), influences fire regimes (Kim et al., 2024) and shapes forest–tundra mosaics (Viereck, 1973) that may decrease the absolute forest cover in the REVEALS reconstructions. Thus, permafrost can maintain unsuitable ground conditions and can cause significant temporal lags in postglacial forest expansion (Herzschuh et al., 2016; Li et al., 2022b; Tian et al., 2018).

In Eastern Europe, bias correction slightly improves pre-Holocene agreement, but the emulation substantially underestimates Holocene tree cover and fails to capture the reconstructed mid-Holocene rise and 4 ka maximum. This indicates that the model-data mismatch in this region is not related to systematic climate biases. Tree cover in the emulation closely follows the precipitation that is much lower in the bias-correction than in the simulation. This strong relationship may be incorrect, but also none of the other considered climate variables indicate a mid- to late-Holocene peak that could explain the reconstructed late tree cover maximum.

On the Tibetan Plateau, warmer conditions in the emulation would raise tree cover, but reduced precipitation limits its growth. Nonetheless, emulated trends better match reconstructions, highlighting the combined influence of multiple climate variables on tree-cover dynamics and that JSBACH in principle captures the vegetation-climate relationships in this region well.

Overall, bias correction is most effective in regions with strong model–data mismatches, where systematic climate biases are the main cause of discrepancies, such as in Western Canada, Central North America or Northern China. In these regions, the MAE can be substantially reduced, but whether the correlation to the reconstructed tree cover is improved, depends on the relationship between the climate variables and there joint effect on the tree cover. The interactions among climate variables (or with CO₂) are often non-linear, so that changes in one variable can alter the importance or sensitivity of the tree cover response to another. Significant mismatches persist in regions where non-climatic factors, such as missing processes in the model, cause the mismatch (Eastern Europe, Central Siberia). The effectiveness of bias correction is further limited in regions with inappropriate climate-temperature relationships, either due to incorrect dependencies or an overestimated tree-cover sensitivity to temperature changes. In addition, the applied bias correction can not reduce errors in the temporal temperature evolution. Therefore, the mismatch in the timing of the Holocene tree cover maximum persists across large parts of the Northern Hemisphere in the emulation with bias-corrected climate, particularly in boreal latitudes. This can partly be expected because the applied bias correction removes only stationary offsets.

At the same time, the results indicate that the absence of a mid-Holocene vegetation maximum is unlikely to be explained solely by biases in the simulated climate. Instead, they point toward structural limitations in the representation of boreal vegetation dynamics and associated land–surface feedbacks in the model, which the emulator does not correct for.

This interpretation is consistent with findings from transient Holocene simulations. Several studies highlight the importance of vegetation feedbacks for Holocene climate evolution, in particular mid-Holocene warmth in the Northern Hemisphere (Thompson et al., 2022). However, models with dynamic vegetation frequently underestimate the magnitude and spatial extent of vegetation changes inferred from proxy records (Dallmeyer et al., 2019; Harrison et al., 2015; Tian and Jiang, 2013), potentially limiting the strength of these feedbacks.

Transient Holocene simulations also show strong regional differences in the timing and magnitude of the Holocene Thermal Maximum, particularly in boreal regions. While models agree reasonably well in regions strongly influenced by residual ice sheets, substantial model–proxy discrepancies persist in areas such as Alaska and Siberia (Zhang et al., 2017, 2018). In these regions, simulated temperature trends during the Early Holocene often differ from pollen-based reconstructions (i.e. warming in reconstructions, cooling in models), highlighting persistent uncertainties in both climatic forcing of vegetation and regional feedbacks (Hao et al., 2025).

In addition, the varying impact of bias correction across regions and time periods may at least partly result from the fact that different climatic factors constrain vegetation at different times. While temperature corrections may improve model–data agreement particularly during the deglaciation, Holocene mismatches may relate to insufficient representation of precipitation dynamics.

4 Summary and Conclusion

Over the last 20 000 years, global warming, driven by orbital changes, rising CO₂, and ice sheet retreat, caused major shifts in climate and a reorganisation of vegetation. While pollen-based reconstructions show significant biome shifts since the Last Glacial Maximum (LGM), the vegetation dynamics and the impact of the climate drivers on the vegetation dynamics remain unclear. The new large-scale pollen-based reconstruction for the Northern Hemisphere by Schild et al. (2025) provides a valuable basis for evaluating Earth system models as it provides quantitative, gridded tree cover reconstructions, improving the comparability to the model output and therewith our understanding of vegetation responses, model-data mismatches, and model performance. We have compared this dataset with a transient simulation by the Max-Planck-Institute Earth System Model MPI-ESM (Kleinen et al., 2023a, b) to answer the following questions:

• a.

  How well does the model reproduce reconstructed tree cover patterns?

  Pollen-based REVEALS reconstructions show regionally variable tree cover during the Last Glacial Maximum (LGM), with high values in parts of eastern North America and Asia. Tree cover expanded markedly during deglaciation, especially north of 50° N. While MPI-ESM broadly captures these trends, it simulates consistently higher tree cover in Europe and East Asia and lower tree cover in boreal North America and Siberia. MPI-ESM tends to simulate too high tree cover fractions in sparse forest and too little tree cover fractions in dense forests. Discrepancies are most pronounced from mid-Holocene to present, partly due to missing land-use in MPI-ESM. Reconstructed and simulated tree cover show overall strong correlations in Europe, except in the Mediterranean. North America and Asia indicate regional mismatches, with high MAE in Siberia, the boreal latitudes in America and the forest-steppe transition zones in Central North America and East Asia.

• b.

  What are the drivers of the tree cover dynamic?

  Statistical emulator and correlation analyses demonstrate clear regional patterns in the climatic controls of simulated vegetation. In high latitudes, summer temperature is the dominant driver, while in tropical and subtropical regions, precipitation is key. In temperate zones, both temperature and precipitation are important, with regional differences. MPI-ESM shows consistent spatial patterns, but the comparison with REVEALS reconstructions indicates notable discrepancies. For REVEALS, simulated precipitation and cold season temperature and the CO₂ change better explain tree cover changes – even in boreal regions – contrasting with the model's emphasis on warm season temperature. This suggests limitations in MPI-ESM's climate simulation or bioclimatic assumptions. Both, reconstructions and simulations show a shift from energy-limited regimes during the deglaciation to moisture-limited regimes during the early and mid-Holocene, followed by regrowth of energy-limited areas in the late Holocene.

• c.

  Can we infer differences in the climate-vegetation relationships between areas with model-data agreement and model-data-disagreement?

  We note that low model–data agreement does not necessarily indicate model deficiencies, but may also reflect local processes mirrored in the reconstructions. Moreover, our analysis indicates that grid cells with high model-data agreement in the tree cover dynamic differ significantly in their climate sensitivity from those with poor agreement. High correlations are associated with a largely linear vegetation response to multiple climate drivers in the model, including CO₂, precipitation, soil moisture, and cold-season temperature. Vegetation dynamics in these grid cells are dominated by the overall (hemispheric) deglacial warming trend. The precipitation and cold season temperature themselves are strongly correlated with the CO₂ signal. In contrast, poor agreement is often linked to more complex, non-linear responses dominated by summer temperature, and rather indicate regional tree cover trajectories.

  Grid cells with high MAE cluster in the mid-latitudes and the time-series start later. Interestingly, high-MAE cells exhibit stronger correlations with cold-season temperature and CO₂, suggesting a too strong sensitivity of the tree cover to CO₂ in the model. The dominant climatic drivers, however, do not differ systematically between the high- and low-MAE- groups.

  Grid cells where the model overestimates variance show strong linear correlations between tree cover and CO₂ or Tc, as well as stronger P–CO₂ and Tw–CO₂ linkages than grid cells with similar variance. Likewise, cells where variance is underestimated exhibit stronger correlations with Tw, Tc, P, and CO₂ than cells with similar variance. Variance biases mainly reflect the tree cover response to CO₂ changes: Linear responses lead to overestimation, non-linear responses to underestimation, suggesting a too strong sensitivity of the tree cover to CO₂ changes in MPI-ESM. Notably, grid cells with similar model-data variance have lower tree cover correlations and lower MAEs than the other groups, indicating low temporal coherence and demonstrating that correct amplitude does not guarantee realistic dynamics.

• d.

  What can we learn from this study?

  Global vegetation modelling with Earth System Models remains challenging, as these models prioritize broad climatic functioning over detailed regional patterns. While the model captures the major large-scale shifts from open LGM landscapes to Holocene forests, the regional dynamics are often misrepresented, showing low correlation with reconstructions or mismatches in variability and absolute cover. Reconstructions typically use regionally differentiated definitions and vegetation characteristics to aggregate thousands of taxa, whereas global models rely on broad plant functional types constrained by universal, static bioclimatic limits derived from modern human-affected vegetation. These fixed bioclimatic limits in the models seems to lead to a too strong dependence of the simulated vegetation on the warm season temperature and to a misrepresentation of the non-linear tree cover response to climate. Furthermore, the sensitivity of the tree cover response to CO₂-changes seem to be too strong in this model version, expressed by higher MAEs and overestimated variance under strong linear coupling of the tree cover to CO₂. Our study also shows that the climate–vegetation relationship can regionally be well captured but that systematic climate biases can still produce large mismatches between simulated and reconstructed tree cover. Two mechanisms may contribute to these discrepancies. First, deficiencies in simulated climate trajectories can lead to incorrect climate forcing of the vegetation model. In the temperate forest–steppe transition zones of Central North America and East Asia, the model is systematically too wet, enabling tree growth much earlier than observed and causing major offsets. Bias correction can reduce the mismatch to the reconstructions here. At high northern latitudes, the model is too cold and warm season temperatures peak during the early Holocene, which distorts the simulated climate forcing of tree cover and contributes to unrealistic vegetation dynamics, e.g. in Western Canada.

  Second, limitations in vegetation sensitivity to climate may further amplify these mismatches. In combination with the global and fixed bioclimatic thresholds used in the vegetation model, the cold bias likely alters the tree cover sensitivities to temperature in JSBACH (and thus in the emulator trained on it). As a result, bias correction does not improve the agreement with the reconstructions in all boreal regions, but instead increases the model-data mismatch, e.g. in Siberia, indicating a misrepresented climate forcing and an unrealistic representation of the boreal forest dynamic. Missing processes in the model, such as the impact of permafrost on the tree cover dynamic, further increase the model-data differences. Consequently, the model is not able to capture the mid-Holocene forest maximum observed across large parts of the Northern Hemisphere. Instead, tree cover peaks during the Early Holocene in the model due to the strong sensitivity of the simulated tree cover to warm season temperature, which closely follows the insolation forcing. However, the underlying mechanisms and consequences require further analysis and will be addressed in a follow-up study.

  In addition to the sources of error in the model, the reconstruction method has shortcomings that are not explicitly discussed here. However, the question arises to which extent REVEALS is suitable for reconstructing forests in open landscapes. Especially in regions with sparse forests, such as the boreal latitudes, REVEALS greatly overestimates the proportion of forest. The fact that bare ground cannot be estimated is problematic and contributes to the model-data mismatch.

• e.

  How can we improve Dynamic Vegetation Models to better simulate past and future tree cover changes

  Current Dynamic Vegetation Models (DVMs) capture global vegetation responses and their role in the climate system, but they are less reliable for detailed regional or local vegetation patterns, particularly with the low spatial resolutions used in ESMs. Most DVMs employ overly simplified soil schemes and lack explicit representation of permafrost processes. In cold and arid regions, vegetation dynamics are often constrained by soil processes and not just by climate. Without realistic soil and permafrost modeling, DVMs tend to overpredict tree establishment at treelines or in permafrost peatlands, leading to inaccurate biogeographical predictions.

  Vegetation dynamics are rarely linear and are often governed by unknown or context-dependent thresholds. Transition zones between biomes are particularly difficult to model due to complex interactions of multiple limiting factors. Dynamic thresholds, rather than static ones, are essential to capture these processes more accurately.

  Plant Functional Types (PFTs) are often too coarse to capture regional vegetation responses, especially under strongly changing climatic conditions or in transition zones between biomes. By grouping species into broad categories, PFTs can overlook important differences in traits that determine survival, growth, and reproduction under stress, such as drought tolerance, frost resistance, rooting depth, or fire-adaptive characteristics. Trait-based approaches, in contrast, allow models to explicitly represent these functional differences and the variability within species or populations, but are still in an early phase of development. They would enable a more mechanistic simulation of how vegetation responds to environmental changes.

  Improving DVMs for palaeoclimate studies, thus, requires (i) integrating realistic soil and permafrost dynamics, (ii) accounting for non-linear processes and dynamic thresholds, (iii) enhancing disturbance modules, (iv) moving from static PFTs to trait-based or individual-based approaches, and (v) increasing regional resolution and flexibility.

Appendix A

The dataset of REVEALS-based tree-cover reconstructions includes the reconstructed fraction of tree cover as well as a mask identifying valid grid cells. A grid cell is considered valid if the reconstruction meets all requirements of the REVEALS model. The REVEALS approach is subject to several constraints and assumptions; for example, pollen records must originate from lakes that are sufficiently large for the model to produce reliable estimates. Grid cells that do not fulfil all criteria are marked as invalid. After proving and comparing the dataset containing all reconstructions and the dataset containing valid grid cells only with the MPI-ESM simulation, we decided to use the entire dataset, i.e. valid and invalid grid cells. We therefore consider the inclusion of these additional grid cells to be justified, as it increases the spatial coverage of the dataset by retaining approximately 23 % more grid cells (Fig. A1c) that would otherwise be excluded, with the strongest impact in Asia, the continent with the fewest available reconstructions. Excluding these grid cells would disproportionately reduce data coverage in this region and further limit its representation in large-scale analyses. Importantly, the deviations from the model remain small (Fig. A1a). In grid cells that fully comply with all REVEALS requirements, the mean absolute error is even smaller when taking all grid cells into account than only in the valid grid cells (Fig. A1b). This indicates that incorporating a larger number of records improves the representation of hemispheric tree cover patterns and diversity, leading to more robust large-scale reconstructions.

We furthermore decided to use the dataset starting at 19 500 years before present, which is an extension of the dataset presented in Schild et al. (2025).

This allows us to compare tree cover dynamics during the deglaciation and the Holocene. Although the applicability of the REVEALS model to periods prior to 14 ka remains debated – primarily due to uncertainties regarding the temporal stability of taxon-specific relative pollen productivity (RPP) estimates – changes in taxa productivity and dispersal are unlikely to be large enough to substantially alter the reconstruction results. Consequently, other studies have extended REVEALS-based reconstructions beyond the Holocene into the last glacial period (Kern et al., 2025). Thus, we do not consider the pre-14 ka reconstructions to be inherently less reliable than those for the Holocene. In particular, the deviations from the model, quantified as mean absolute error (MAE), increase towards more recent times (Fig. A2). This pattern argues against systematically larger distortions in the earlier interval and supports the robustness of the reconstructions for this period.

Figure A1 (a) Difference in MAE between the simulated and reconstructed tree cover between the comparison with the REVEALS-datasets with only valid grid cells and with all grid cells, (b) mean difference over continents. (c) differences in grid cell number [%] per continent.

Figure A2 MAE between MPI-ESM and REVEALS over time (red) and number of available grid cells per time-steps (black).

Appendix B: Preparation of the bias corrected climate prescribed to the trained emulator

In a second set of emulations, we use the emulator trained with the MPI-ESM climate and force it with an idealised, bias-corrected climate, assuming that grid cell specific model biases (e.g. due to differences in the orography) do not change over time (Maraun and Widmann, 2018). We use the delta method following (Beyer et al., 2020) to downscale and bias-correct the simulated climate. For temperature variables, we calculated anomalies between past time-slices and the modern time-slice in the model, bilinearly interpolated these anomalies to a 0.5°-grid and added these to the CRU-TS 4.08 (University of East Anglia Climatic Research Unit et al., 2024; Harris et al., 2020) derived temperature for the period 1931 to 1970. To bias-correct the precipitation, we applied the delta method in a multiplicative way, which is seen to be more appropriate for moisture changes as precipitation covers different orders of magnitudes, depending on the region (Beyer et al., 2020; Maraun and Widmann, 2018). This method typically involves calculating relative anomalies by dividing simulated precipitation at each timestep by a reference period, and then scaling these factors by a modern observational baseline. However, in arid regions, this approach can produce unrealistically high precipitation values, when the present-day simulated precipitation is extremely low, resulting in large delta factors when dividing by these small reference values. When these exaggerated anomalies are then applied to modern observational baselines (which may already be low), the resulting scaled precipitation becomes unreasonably inflated. To cope with this, we adopted a modified delta approach, using the simulated climate at 4 ka as the reference period in regions with dry climates. This adjustment provides a more stable denominator and reduces the risk of over-amplifying anomalies in dry regions. To harmonize the results with present-day conditions, and ensure consistency across the dataset, we applied a correction term based on the relative difference between the simulated precipitation at 0 and 4 ka This effectively normalizes the delta factors to present-day conditions while still preserving the benefits of using the 4 ka baseline in dry regions. The final anomalies were then applied to the modern observational dataset to produce bias-adjusted precipitation fields. As the CRUTS4.0.8 data set does not contain a soil moisture field, we just interpolated MPI-ESM soil moisture to the 0.5°-grid, without correcting biases. Time series of Northern Hemisphere mean precipitation, warm- and cold month temperature for both the original and bias-corrected climate, together with the corresponding difference in mean tree cover between the bias-corrected emulation and the original simulation are shown in Fig. B1 in the Appendix B. In addition, spatial maps of the climate variables at the time of maximum change (∼ 12 ka BP), as well as maps of emulated tree cover, simulated tree cover, and their difference for three representative time slices, are provided in Figs. B2 and B3 in the Appendix B.

Figure B1 Time series of Northern Hemisphere mean precipitation (precip), temperature of the warmest (Twarm) and coldest month (Tcold) for the original (MPI) and bias-corrected (EMU) climate, together with the corresponding difference in mean tree cover between the bias-corrected emulation and the original simulation. Shown is also the global atmospheric CO₂ concentration that is used as predictor in the emulation. The dotted lines mark the three time-slices for which maps of tree cover are provided in Fig. B2.

Figure B2 Difference in annual mean precipitation (Pann) [mm/year] and in the temperature of the warmest (Tw) and coldest (Tc) month [K] at 12 000 years BP between the bias-corrected and the original simulation. The difference was calculated after remapping the simulation onto a 0.5° grid.

Figure B3 Tree cover fraction at different time steps according to the emulation with bias-corrected climate and the MPI-ESM simulation. The bottom panel also shows the difference between the two; for this, the emulation has been remapped to a T31 grid.

Appendix C: Testing significant difference between the metric-groups

To analyse whether the distributions of climate correlations and climate–tree cover dependencies differ among the mae-, r- and var-groups, we used a Kruskal–Wallis test. There results are shown in Figs. C1, C2 and C3. To select important variables for the PCA-analysis, we additionally perform Wilcoxon rank-sum tests to determine variables with significant differences between r-group 1 and 3, mae-group 1 and 4. The results are shown in Table C1.

Figure C1 Distribution of the coefficients/variables in the different r-groups (1 = positive correlation, 2 = no correlation, 3 = negative correlation) with significance (Kruskal-Wallis, ${}^{\ast \ast \ast }$ 0.001, ${}^{\ast \ast }$ 0.01, ^∗ 0.05, • 0.1), Hereby, “F_” corresponds to the F-coefficients in the GAM-analysis and “R_” to the Pearson correlation coefficients in the model, or in the reconstructions (“_rec”). The climate variables are Tw: Temperature of the warmest month, Tc: temperature of the coldest month, P: precipitation, Sw: Soil moisture. Valid data describes the number of available time-steps and First_valid_year the first available time-step in the reconstructions. R_Trees is the correlation between the reconstructed and simulated tree cover. “BiasCorr_” denotes the Correlation, MAE and variance differences between the bias-corrected emulation and the reconstructions.

Figure C2 Distribution of the coefficients/variables in the different mae-groups (1 = positive correlation, 2 = no correlation, 3 = negative correlation) with significance (Kruskal-Wallis, ${}^{\ast \ast \ast }$ 0.001, ${}^{\ast \ast }$ 0.01, ^∗ 0.05, • 0.1), Hereby, “F_” corresponds to the F-coefficients in the GAM-analysis and “R_” to the Pearson correlation coefficients in the model, or in the reconstructions (“_rec”). The climate variables are Tw: Temperature of the warmest month, Tc: temperature of the coldest month, P: precipitation, Sw: Soil moisture. Valid data describes the number of available time-steps and First_valid_year the first available time-step in the reconstructions. R_Trees is the correlation between the reconstructed and simulated tree cover. “BiasCorr_” denotes the Correlation, MAE and variance differences between the bias-corrected emulation and the reconstructions.

Figure C3 Distribution of the coefficients/variables in the different variance-groups (1 = overestimated in MPI-ESM, 2 = similar variance, 3 = underestimated) with significance (Kruskal-Wallis, ${}^{\ast \ast \ast }$ 0.001, ${}^{\ast \ast }$ 0.01, ^∗ 0.05, • 0.1), Hereby, “F_” corresponds to the F-coefficients in the GAM-analysis and “R_” to the Pearson correlation coefficients in the model, or in the reconstructions (“_rec”). The climate variables are Tw: Temperature of the warmest month, Tc: temperature of the coldest month, P: precipitation, Sw: Soil moisture. Valid data describes the number of available time-steps and First_valid_year the first available time-step in the reconstructions. R_Trees is the correlation between the reconstructed and simulated tree cover. “BiasCorr_” denotes the Correlation, MAE and variance differences between the bias-corrected emulation and the reconstructions.

Table C1 Significant difference in the pairwise comparison between the correlation, mae and variance-groups using the Wilcoxon rank-sum test and Kruskal-Wallis test (${}^{\ast \ast \ast }$ 0.001, ${}^{\ast \ast }$ 0.01, ^∗ 0.05, • 0.1, – not significant). Hereby, “F_” corresponds to the F-coefficients in the GAM-analysis representing the non-linear relationship and “R_” to the Pearson correlation coefficients in the model (i.e. linear relationship), or in the reconstructions (“_rec”). The climate variables are Tw: Temperature of the warmest month, Tc: temperature of the coldest month, P: precipitation, and Sw: Soil moisture. Valid data describes the number of available time-steps and First_valid_year the first available time-step in the reconstructions. Lon and lat define the geographical position. R_Trees is the correlation between the reconstructed and simulated tree cover. R_emu is the correlation between the simulated and emulated tree cover.