# Holocene potential natural vegetation in Europe: Evaluating the model spread with three dynamical vegetation models

Isabeau A Bertrix, Hisashi Sato, Nicolas Viovy, Hans Renssen, Didier M Roche

► To cite this version:

Isabeau A Bertrix, Hisashi Sato, Nicolas Viovy, Hans Renssen, Didier M Roche. Holocene potential natural vegetation in Europe: Evaluating the model spread with three dynamical vegetation models. Quaternary Science Reviews, 2026, 376, pp.109813. 10.1016/j.quascirev.2026.109813. hal-05495604

HAL Id: hal-05495604

https://hal.science/hal-05495604v1

Submitted on 5 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Holocene potential natural vegetation in Europe: Evaluating the model spread with three dynamical vegetation models

Isabeau A. Bertrix<sup>a,b,*</sup> <sup>ID</sup>, Hisashi Sato<sup>c</sup> <sup>ID</sup>, Nicolas Viovy<sup>a,b</sup>, Hans Renssen<sup>d</sup> <sup>ID</sup>, Didier M. Roche<sup>a,b,e</sup> <sup>ID</sup>

## ARTICLE INFO

Handling Editor: Qiuzhen Yin

Keywords:

Holocene

DGVM

Paleo-vegetation

Climate modeling

## ABSTRACT

The period of the early Holocene in Europe is marked by climate warming as Earth comes out of the last glacial period and is followed by the emergence of agriculture and animal husbandry in the second half of the period. Increased human influence had profound impacts on the land surface, but the Holocene climate evolution also drove some changes that are intertwined with it. Deciphering the role of each in the vegetation evolution is becoming more difficult as one progresses to the earlier parts of the Holocene here human induced impacts were fainter. Within this general context, we aim at understanding how much Dynamical Vegetation Models (DGVMs) differ in their representation of Potential Natural Vegetation (PNV) in Europe during the mid-to late Holocene (8.5 k.a. BP to 1900 A. D.). We ran three different DGVMs, SEIB-DGVM, ORCHIDEE-DGVM and CARAIB, in Europe, for six time-slices and forced them with identical climatic inputs obtained from the iLOVECLIM Earth system model (downscaled and bias-corrected). Results are then compared with pollen-based reconstructions from the TERRANOVA database. Overall, the three models have a similar performance in representing the pollen-derived vegetation cover at the european scale. However, their results are largely different at regional scales, particularly in mountainous areas and in boreal regions. They also show a very large spread in simulated PFT diversity at the grid cell scale, highlighting the impact of each model's internal dynamics on the results. On a global scale, they all agree on a decreasing match of their results with pollen base reconstructions over time, indicating - as expected - an increase of the human pressure on the landscape.

## 1. Introduction

Climate always had a strong influence on both vegetation and humanity. For example, in Europe, the drop in temperature during the last glacial maximum (LGM) resulted in a drop in human population, as shown by Tallavaara et al. (2015) by using climate envelope modeling tools and modern ethnographic datasets. Climate variation also influenced landscape, which is highly sensitive to temperature, precipitation, CO2 levels, humidity and solar radiation. Another big influencer of landscape are human activities. Indeed, if we look at the Holocene period, human activities such as deforestation, hunting and agriculture, have led to substantial changes in the european vegetation cover

(Nikulina et al., 2024), as well as a rise in CO2 levels more recently (Petit and Raynaud, 2020). Moreover, land cover is an integral component of the climate system that can lead to two types of feedback loops: biogeochemical (albedo, thermal conductivity) and biogeochemical (carbon dioxide absorption or release during photosynthesis or decomposition). Indeed, forests have a lower albedo than crops or pasture, and this difference is even more marked in the presence of snow. Furthermore, when they grow, trees convert atmospheric carbon into carbon that they store in their trunk, branches, leaves, roots and in the soil. This carbon, as well as dust aerosols, is released into the atmosphere during deforestation due to biomass burning. Human societies impact on vegetation thus modifies the climate directly and through the mediation of the

<footnote><sup>a</sup> Laboratoire des Sciences Du Climat et de L'Environnement (LSCE), CEA, CNRS, UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> Institut Pierre-Simon Laplace (IPSL), Université de Versailles Saint-Quentin en Yvelines, Guyancourt, France</footnote>

<footnote><sup>c</sup> Research Institute for Global Change (RIGC), Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokohama, Japan</footnote>

<footnote><sup>d</sup> Department of Natural Sciences and Environmental Health, University of South-Eastern Norway, Bø, Norway</footnote>

<footnote><sup>e</sup> Earth and Climate Cluster, Faculty of Science, Vrije Universiteit Amsterdam, Amsterdam, the Netherlands</footnote>

<footnote>* Corresponding author. Laboratoire des Sciences Du Climat et de L'Environnement (LSCE), CEA, CNRS, UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France. E-mail addresses: isabeau.bertrix@lsce.ipsl.fr (I.A. Bertrix), hsatoscb@gmail.com (H. Sato), nicolas.viovy@lsce.ipsl.fr (N. Viovy), hans.renssen@usn.no (H. Renssen), didier.roche@lsce.ipsl.fr (D.M. Roche).</footnote>
carbon cycle. Additionally, for a given set of climate and human impact conditions the vegetation cover is not ubiquitous but is also a result of the past history. Most of the studies made on the Holocene vegetation have been made on the Green Sahara (Lu et al., 2018).

A good candidate to quantify the human influence on vegetation is the comparison between pollen-based reconstructions and Dynamical Global Vegetation Model (DGVMs) simulations. Indeed, pollen-based reconstructions are a good indicator of the type of plants, even sometimes species, that lived at a given period of time on a land surface. They represent the vegetation effectively affected by the combination of human activities and climate conditions. In the absence of specific land-use forcings (for example, Chini et al., 2014), DGVM simulations represent the Potential Natural Vegetation (PNV) that could have existed in the absence of large scale human intervention. DGVMs use climatic data as input to compute the carbon uptake, vegetation cover and competition and related surface conditions. As such, they are sensitive to climatic variations. Under such natural forcing conditions, the differences between pollen-based reconstructions and DGVM simulations should represent the impact of human activities.

In a previous study using this method Zapolska et al., 2023a has shown that humans had a large impact on vegetation cover even before agriculture started. Indeed, simulations with the CARAIB DGVM model forced by bias-corrected climate model outputs showed that vegetation cover in Europe differed significantly from the state of potential natural vegetation even before 6 k.a. BP (Before Present) and challenges the hypothesis that vegetation during the mid-Holocene was in a relatively natural state. This study shows an increase in the difference between PNV and pollen reconstruction with time. Strandberg et al. (2022) studied the difference between PNV and pollen-based reconstruction in Europe at 6 k.a. BP using LPJ-GUESS as DGVM, and showed that PNV was mostly composed of forests, while pollen-based reconstructions has a large component of open lands. In another study at global scale, Dallmeyer et al. (2022), Dallmeyer et al. (2023) also showed disagreement between numerical simulations with the coupled MPI-ESM1.2 and pollen reconstructions with forest expansion post-deglaciation occurring 4000 years before what the pollen reconstructions indicate in the Northern Hemisphere.

To progress one step further, we need to identify whether these discrepancies are arising from the capability of our current vegetation models to represent fundamental aspects of vegetation evolution during climate change or our ability to apprehend the time-cumulative impacts of early human societies on vegetation, or both. Indeed, vegetation models differ in sensitivities to climate and CO2 levels, which results in uncertainties in simulation results. Comparing results of different DGVM models seems a good way to minimize errors, as it has been used in several studies. For example, on the Holocene period, Li et al. (2019) runs two DGVM models, LPJ-GUESS and VECODE on a global scale, Hopcroft et al. (2017) runs three models, JULES version 4.1 SDGVM, and LPJ on the green Sahara. Our study seeks to further intercompare three different DGVMs over the european area by studying the simulated potential natural vegetation response to climate forcing during the Late Holocene period (8.5 k.a. BP to 1900 A. D.). By comparing the different responses we seek to get a grasp at the common and diverging patterns (in space and time) to assess the common response to climate change; by comparing their responses to an appropriately clustered pollen compilation, we aim at adding a quality assessment for the different regions. We also want to compare the DGVMs in general which has rarely been done.

The first section describes the different methods used during this study, the DGVM models and the adaptations we incorporated, as well as an introduction to the pollen-based reconstruction. The second section provides a detailed description of the numerical experiment setup, the third section presents the results obtained, and the fourth section develops a discussion around the ability of DGVM models to describe PNV.

## 2. Methodology

In order to study the impact of long-term climate changes on the vegetation of the second-half of the Holocene in Europe, we need a climatic dataset. We choose using climate data simulated with the iLOVECLIM Earth system model for its rapidity and flexibility, required for a long term paleo climate study. iLOVECLIM results have to be bias-corrected, before being used as inputs for a DGVM. We also compared the results of the DGVM simulations with pollen-based reconstructions.

### 2.1. Climatic data

Three elements are needed to implement our approach: a climatic dataset for each timeslice we are investigating, the bias correction method and a reference climatic dataset to anchor the bias correction.

#### 2.1.1. iLOVECLIM model description

We used the iLOVECLIM Earth system model (here in version 1.1.5), derived from the original LOVECLIM 1.2 model (Goosse et al., 2010), as revised by Caley and Roche (2013). It belongs to the class of Earth System Models of Intermediate Complexity (EMIC), allowing much faster computation than with the more computationally intensive General Circulation models (GCMs). Efficiency is an integral part of our approach since we want to run the complete Holocene time period in the long run: to keep a consistent approach throughout, we thus keep iLOVECLIM as a base model for our timeslices. iLOVECLIM is run in an atmosphere-ocean-vegetation fully coupled configuration, thus including: the atmospheric model, ECBilt, the sea-ice ocean component, CLIO, and the reduced-form dynamic global vegetation model, VECODE.

iLOVECLIM produces results at approximately 5.625°x 5.625°latitude-longitude spatial resolution grid (T21) over the world. In order to obtain a better spatial resolution over our study area, Europe, we make use of the online interactive downscaling method embedded in iLOVECLIM, first described by Quiquet et al. (2018). It is an online dynamical downscaling of temperature and precipitation included in iLOVECLIM, replicating the processes of surface temperature and precipitation computation on a refined vertically extended grid. It allows for the computation of surface temperature and precipitation at any altitude within a given subgrid. For our European context, we use a 0.25°x 0.25°latitude-longitude spatial resolution grid over mainland Europe.

However, as in all climate simulations, iLOVECLIM presents biases and its outputs need to be corrected using observation records before being used in DGVM models.

#### 2.1.2. Bias-correction with CDF-t

The method we used for bias-correction is the “cumulative distribution function transform” (CDF-t), a statistical method based on Quantile Mapping (QM). But contrary to the QM method, which considers that model and observation distributions keep the same shape with time, CDF-t considers that model and observational distributions can evolve and be different. To bias-correct a climatic dataset using the CDF-t method, an observational reference dataset of the same resolution is needed. Three climatic iLOVECLIM outputs have been corrected using this method: the surface air temperature, the precipitation, and the relative humidity. This method has been tested for the Holocene period and shows good results (Zapolska et al., 2023b).

#### 2.1.3. Reference observation dataset

The EWEMBI climatic dataset was used as the observational reference dataset during the bias-correction (Lange, 2019). It is a re-analysis combination resulting in a daily temporal resolution climatic dataset from 1979 to 2016 at 0.5°horizontal resolution for the entire globe. Before using it, we bilinearly interpolated the data on the 0.25°x 0.25°European grid.
## 2.2. DGVM models

### 2.2.1. SEIB-DGVM

Spatially Explicit Individual Based (SEIB-DGVM, here in version 3.03) is a dynamic vegetation model, which aims to simulate the transient impacts of climate change on the terrestrial ecosystem and land-atmosphere interactions (Sato et al., 2007). It contains mechanical or empirical algorithms describing terrestrial physical processes (hydrology, radiation, air, etc.), plant physiological processes (photosynthesis, respiration, growth, etc.) and plant dynamic processes (establishment, mortality, disturbance). Unlike other existing DGVMs, SEIB-DGVM only simulates the local interactions of individual trees within a spatially explicit virtual forest; several sample plots are placed in each grid box where the growth, competition for light and decay of each individual tree within a group of trees are then calculated. While accurately representing plant dynamics is crucial for capturing the time lag in vegetation distribution responses to climate change, the specific way plant dynamics are incorporated varies significantly among existing DGVMs (Fisher et al., 2018). SEIB-DGVM takes the most direct approach by modeling individual trees competing for space and light, thereby explicitly simulating vegetation dynamics.

SEIB-DGVM simulates 16 PFTs: 6 tropicals, 3 temperates, 5 boreals (with 2 specifically siberian) and 2 herbaceous PFTs. As no tropical or siberian pollen have been recorded in Europe during the Holocene, we chose to disable the 6 tropical PFTs and the 2 siberian PFTs.

As previous studies using SEIB-DGVM were mainly applied either on global scales (Tong et al., 2022) or on african or siberian (Sato et al., 2023) zones, there was no mediterranean PFT (contrary to the pollen records). We decided to create one by modifying the existing "Temperate broad-leaved evergreen" PFT. The parameters associated with this new PFT can be found in Table S8. We also add the drought parameter, which is the state of water satisfactory for photosynthesis (0.0–1.0), from SEIB-DGVM 3.10 version (Sato et al., 2023) to the SEIB-DGVM 3.03 version. Following Sato et al., 2023, we choose to associate coniferous PFTs (temperate or boreal) with an optimal state of water satisfactory of 0.75, a minimum state of water satisfactory of 0.20, and a maximum state of water satisfactory of 0.90; deciduous PFTs (temperate or boreal) with an optimal state of water satisfactory of 0.85, a minimum state of water satisfactory of 0.30, and a maximum state of water satisfactory of 0.99. An other change in SEIB-DGVM is the modification of the net CO2 assimilation rate per needle area for the coniferous. Indeed, it has a strong dependence on leaf age, as measured in Robakowski and Bielinis (2017), which is not taken into consideration in the SEIB-DGVM. Consequently, we decreased the net CO2 assimilation rate per needle area from 9.0 µmol.m⁻²·s⁻¹ to 7.0 µmol.m⁻²·s⁻¹ for both the temperate needle-leaved evergreen and the boreal needle-leaved evergreen PFTs.

Finally, as SEIB-DGVM was never used to study paleoclimate, it had to be adapted to consider astronomical parameters. Indeed, SEIB-DGVM calculates the solar declination and the distance between the Earth and the Sun for a given day of a year, in a way that is only applicable at present day. To make this calculation valid for paleo times, we need to obtain it only using the astronomical parameters (from Berger and Loutre, 1999): the obliquity $e$, the eccentricity $e$, and the climatic precession $\omega$. The details of the calculation are described in supplementary.

### 2.2.2. ORCHIDEE-DGVM

ORCHIDEE is the land surface model of the IPSL (Institut Pierre Simon Laplace) Earth System Model. In this study, we used ORCHIDEE version 2.2 with DGVM activated. ORCHIDEE is based on three different existing models: SVAT SECHIBA, that describes exchanges of energy and water between the atmosphere and the biosphere, and the soil water budget; LPJ-DGVM for the parameterizations of vegetation dynamics (fire, sapling establishment, light competition, tree mortality, and climatic criteria for the introduction or elimination of PFTs); and STO-MATE, for processes such as photosynthesis, carbon allocation, litter decomposition, soil carbon dynamics, maintenance and growth

respiration, and phenology. The DGVM part of ORCHIDEE simulates the growth of 10 PFTs, and have also a bare soil PFT: 2 tropical, 3 temperate, 3 boreal, and 2 herbaceous PFTs. As no tropical pollens as been recorded on Europe during the Holocene, we choosed to disable the two tropical PFTs.

### 2.2.3. CARAIB

Carbon Assimilation In the Biosphere (CARAIB) model is a dynamic vegetation model originally designed to study the role of vegetation in the global carbon cycle, at present and in the past. It is composed of five modules (François et al., 2011) describing respectively (1) the hydrological budget, (2) canopy photosynthesis and stomatal regulation, (3) carbon allocation and plant growth, (4) heterotrophic respiration and litter/soil carbon dynamics, and (5) plant.

Competition and biogeography. In CARAIB, photosynthesis and plant respiration are computed every 2 h, and water and carbon reservoir are updated every day. The model simulates 26 plant functional types (PFTs), which can coexist on the same grid cell: 3 tropical, 3 subtropical, 1 subdesertic, 4 mediterranean, 6 temperate, 6 boreal and 3 herbaceous PFTs. In CARAIB, herbs and shrubs are assumed to grow under the trees, leading to a two layers structure. Each layer has a maximum coverage of 1, meaning that the maximum vegetation fraction of a grid cell is 2. As no tropical pollen has been recorded in Europe during the Holocene, we chose to disable the 3 tropical PFTs and the subdesertic one. Another difference with SEIB-DGVM and ORCHIDEE-DGVM is the presence of shrub PFTs. In order to compare properly the three models, we chose to disable the shrubs PFTs, leaving 16 PFTs (3 herbaceous and 13 tree types). In the post-treatment process, we decided to reclassify CARAIB's remaining PFTs into the 8 SEIB's PFTs. The reclassification is presented below in Table (1).

## 2.3. DGVMs intercomparison methodology

One classical way to intercompare results between different vegetation results is to express them in the form of biomes (Haxeltine and Prentice, 1996; Dallmeyer et al., 2019). Classically, the biomes are computed using the maxvegetfrac of the DGVM model or an equivalent (based on an area fraction of the gridcell). This can't be done in our study as SEIB-DGVM and CARAIB do not compute a fractional PFT cover. This is why we chose to generate biomes according to the methodology from SEIB, which uses as input the dominant PFT and the mean annual maximum Leaf Area Index (LAImax) over the last ten years of simulation at each point. If the LAI max is higher than 2.5 m2.m⁻², the associated biome will be a forest one; if not but still higher than 1.5 m2.m⁻², it will be a woodland one. If the LAI max is lower it will be a grassland/steppe/savanna biome, and if lower than 0.2 m2.m⁻², a desert biome. Polar desert and Arctic/Alpine-tundra biomes are only defined by a threshold in Growing Degree Days with a limit at zero or five Celsius (GDD0 and GDD5) respectively. Biome will be a polar desert if GDD0 <

<table><caption>**Table 1**
CARAIB PFTs classified into SEIB-DGVM PFTs.</caption><thead><tr><td>SEIB-DGVM PFT</td><td>CARAIB PFT</td></tr></thead><tbody><tr><td>Temperate needle-leaved evergreen</td><td>NEg Te cool trees</td></tr><tr><td>Mediterranean</td><td>NEg sMed trees, NEg mMed trees, NEg subTr trees, BEg mMed trees, BEg tMed trees, BEg subTr trees, NSg subTr swamp trees</td></tr><tr><td>Temperate broad-leaved summergreen</td><td>Sg Te cool trees, BSg Te warm trees</td></tr><tr><td>Boreal needle-leaved evergreen</td><td>NEg B/Te cold trees</td></tr><tr><td>Boreal needle-leaved summergreen</td><td>NSg B/Te cold trees</td></tr><tr><td>Boreal broad-leaved summergreen</td><td>BSg B/Te cold trees</td></tr><tr><td>Herbaceous</td><td>C3h, C3d, C4</td></tr></tbody></table>
150, or a Arctic/Alpine-tundra if GDD5 < 370. In order to compare the results of the three DGVM models, we applied the same biome code on ORCHIDEE-DGVM and CARAIB results. As ORCHIDEE-DGVM has a different functioning, the LAlmax needs to be multiplied by the maximum vegetation fraction. The dominant PFT is the PFT having the greatest mean annual Net Primary Production (NPP) over the last ten years of simulation. Again for ORCHIDEE-DGVM, the NPP needs to be multiply by the maximum vegetation fraction.

This methodology is summarized in Fig. 1.

## 2.4. Pollen-based reconstruction: the TERRANOVA database

### 2.4.1. The *reveals* method

In order to compare the results of our European DGVM Holocene simulations, to land-cover reconstructions, we used the TERRANOVA pollen-based reconstruction dataset. It comes from the aggregation of several raw pollen databases (1607 records, mainly from lakes and peatlands) which are interpreted quantitatively using the REVEALS (Regional Estimation of VEgetation Abundance from Large Sites) model (Serge et al., 2023a). The taxa used are only anemophilous taxa, their dispersion is assumed to be isotropic (Serge et al., 2023b). To limit representation biases in the model, data is weighted by the Relative Pollen Productivity (corresponding to the number of pollen grains produced per species), by the average pollen fall speed, by the diameter of the sampling zone and by climatic conditions, following the protocol of Githumbi et al. (2022).

The REVEALS model output, TERRANOVA, is a dataset comprising a total of 378 different coordinates for 31 taxa, on a 1°x1° grid across 30°–71°N, 20°W–47°E (north-western, central Europe, Mediterranean area, and part of the East until 47°E) for 25 contiguous time slices of 100–500 years covering all of the Holocene.

In the TERRANOVA dataset, total cover of plant taxa is always 100 %, as it is challenging to estimate the proportion of bare soil in pollen reconstruction.

### 2.4.2. Classification of pollen taxons into PFTs

As DGVMs simulate the growth of PFTs and not species, we classified TERRANOVA taxa into PFTs in order to compare models and data efficiently. The TERRANOVA taxa are given as Taxons in Tables 2 and 3. We classified the PFTs of CARAIB into those of SEIB-DGVM in Table 1. But SEIB-DGVM and ORCHIDEE-DGVM do not exactly have the same PFTs. Instead of a mediterranean PFT, which can contain both broad-leaved and needle-leaved species, ORCHIDEE-DGVM has a temperate broad-leaved evergreen PFT. This is why Pinus is classified as mediterranean in SEIB-DGVM, and as Temperate needle-leaved evergreen in ORCHIDEE-DGVM. All the other PFTs classified as mediterranean in SEIB-DGVM are classified as Temperate broad-leaved evergreen in ORCHIDEE-DGVM. We reclassified each taxa into one or two of the SEIB-DGVM PFTs, as shown in Table 2, and into one or two of the ORCHIDEE-DGVM PFTs in Table 3. For Juniperus, Ericaceae, Pinus, and Salix, we chose to classify them into two different PFTs, as their ecology

<table><caption>**Table 2**

TERRANOVA taxons classified into SEIB-DGVM PFTs.</caption><thead><tr><th>SEIB-DGVM and CARAIB PFT</th><th>Taxons</th></tr></thead><tbody><tr><td>Temperate needle-leaved evergreen</td><td>Abies alba, Juniperus</td></tr><tr><td>Mediterranean</td><td>Buxus sempervirens, Carpinus orientalis, Ericaceae, Phillyrea, Pinus, Pistacia, Quercus evergreen</td></tr><tr><td>Temperate broad-leaved summergreen</td><td>Alnus glutinosa, Carpinus betulus, Castanea, Corylus avellana, Fagus sylvatica, Fraxinus, Quercus deciduous, Salix, Tilia, Ulmus</td></tr><tr><td>Boreal needle-leaved evergreen</td><td>Ericaceae, Juniperus, Picea, Pinus</td></tr><tr><td>Boreal needle-leaved summergreen</td><td>Betula, Salix</td></tr><tr><td>Boreal broad-leaved summergreen</td><td>Amanranthaceae chenopodiaceae, Artemisia, Calluna vulgaris, Cerealia t, Cyperaceae, Filipendula, Plantago lanceolata type, Poaceae, Rumex acetosa t, Secale</td></tr><tr><td>Herbaceous</td><td></td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Biome conversion of the model results.</figcaption></figure>
<table><caption>**Table 3**

Terranova taxons classified into ORCHIDEE-DGVM PFTs.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>ORCHIDEE-DGVM PFT</td><td>Taxons</td></tr><tr><td>Temperate needle-leaved evergreen</td><td>Abies alba, Juniperus, Pinus</td></tr><tr><td>Temperate broad-leaved evergreen</td><td>Buxus sempervirens, Carpinus orientalis, Ericaceae, Phillyrea, Pistacia, Quercus evergreen</td></tr><tr><td>Temperate broad-leaved summergreen</td><td>Alnus glutinosa, Carpinus betulus, Castanea, Corylus avellana, Fagus sylvatica, Fraxinus, Quercus deciduous, Salix, Tilia, Ulmus</td></tr><tr><td>Boreal needle-leaved evergreen</td><td>Ericaceae, Juniperus, Picea, Pinus</td></tr><tr><td>Boreal needle-leaved summergreen</td><td>Betula, Salix</td></tr><tr><td>Boreal broad-leaved summergreen</td><td>Amanranthaceae chenopodiaceae, Artemisia, Calluna vulgaris, Cerealia t, Cyperaceae, Filipendula, Plantago lanceolata type, Poaceae, Rumex acetosa t, Secale</td></tr><tr><td>Herbaceous</td><td></td></tr></tbody></table>

can be associated with both. For each grid cell, those taxa are classified in one of the two PFTs, according to the fraction of these two PFTs already present on the grid cell considered. At the end, we obtained the percentage of each PFT for a given box. That way, every box can be composed of the sum of the 6 PFTs, for a total of 100 %. No taxa corresponded to the boreal needle-leaved summergreen PFT. Also, the herbaceous PFT contain not only C3 herbaceous taxa, but also heather (Calluna Vulgaris) and C4 herbaceous.

### 2.4.3. Distance between model and data

In order to compare model and data, we also create a set of distance values between PFTs in Table 4 for SEIB-DGVM and CARAIB PFTs compared to TERRANOVA PFTs, and 5 for ORCHIDEE-DGVM PFTs compared to TERRANOVA PFTs. This distance will be used later to estimate how much each model differs with the pollen data. Basically, a distance 0 is when PFTs are exactly the same, 1 is for same climate (temperate or boreal) but different PFTs (for example, a Boreal Broad-Leaved Summergreen instead of a Boreal Needle-Leaved Evergreen), 2 is for different climate but still tree PFT (for example a Temperate Broad-Leaved Summergreen instead of a Boreal Broad-Leaved Summergreen), and.

3 is for herbaceous instead of tree or vice-verça. Note that the combination TeNEg/BNEg has a distance of 1 even if one is boreal and the other one temperate, because in the pollen, the taxa classified as TeNEg or BNEg can all belong to both of those PFTs.

## 3. Numerical experiment setup

During this study, we ran six different time-slice DGVM simulations within the Holocene: 1900 A.D., 1 k.a. 3 k.a., 4

k.a., 6 k.a. and 8.5 k.a. BP (Before Present). To force the different DGVMs, we create climatic inputs consistently

Generated by iLOVECLIM (downscaled and bias-corrected). Each model uses a set of closely related climate fields, as

Shown in Table 7. We also provided appropriate CO2 concentration relative to the Holocene to each model, as well as.

<table><caption>**Table 4**

Distance PFTs SEIB-CARAIB/TERRANOVA.</caption><thead><tr><th></th><th>TeNEg</th><th>Med</th><th>TeBSg</th><th>BNEg</th><th>BNSg</th><th>BBSg</th><th>Herbaceous</th></tr></thead><tbody><tr><td>TeNEg</td><td>0</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>3</td></tr><tr><td>Med</td><td>1</td><td>0</td><td>2</td><td>3</td><td>3</td><td>3</td><td>3</td></tr><tr><td>TeBSg</td><td>1</td><td>2</td><td>0</td><td>2</td><td>2</td><td>2</td><td>3</td></tr><tr><td>BNEg</td><td>1</td><td>3</td><td>2</td><td>0</td><td>1</td><td>1</td><td>3</td></tr><tr><td>BNSg</td><td>2</td><td>3</td><td>2</td><td>1</td><td>0</td><td>1</td><td>3</td></tr><tr><td>BBSg</td><td>2</td><td>3</td><td>2</td><td>1</td><td>1</td><td>0</td><td>3</td></tr><tr><td>C3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>0</td></tr><tr><td>C4</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>0</td></tr></tbody></table>

Astronomical parameters (obliquity, eccentricity and climatic precession). All is summarized in Table 6.

It should be noted that iLOVECLIM runs already with a simplified vegetation model included (VECODE), to provide the first order response of climate in a coupled mode. Hence, our climate forcing is dependent upon the vegetation results of VECODE. This type of setup has been already used (Li et al., 2019) and has shown that as long as we concentrate on the first order response (which is the case when looking at the climate impact in iLOVECLIM), it was comparable in vegetation models of different complexities. In addition, since the process here involves a further step with bias correction, the obtained climate forcing is not simply a direct response to that vegetation coupling.

Here, the daily surface air temperature amplitude and the wind speed come from the EWEMBI dataset and were kept constant for each period. Indeed, we couldn't bias-corrected the wind as it is not a spatial scalar field but a spatial vector field, hence we chose to use the wind from observation dataset. Concerning the daily temperature amplitude, this is a fundamental aspect of the iLOVECLIM model which does not compute day/night variations (nor any of the processes attached to this time scale) but a permanent gray insolation, resulting in the absence of a daily temperature range. We thus need to fix the daily temperature amplitude and the simplest is to keep it constant to observational values. Each simulation was run until reaching equilibrium for the simulated vegetation (300 years for SEIB-DGVM, 250 years for ORCHIDEE-DGVM, 380 years for CARAIB), using the iLOVECLIM climatic dataset loop-repeated.

# 4. Results

## 4.1. Model intercomparison

### 4.1.1. Biome comparison

As seen in Fig. 2, the biome results have some common patterns in all DGVM models: boreal biomes are located at the high latitudes and temperate biomes at the mid-latitudes as expected. The spatial coverage of the different biomes is however different, highlighting potential different processes behind the distribution: SEIB-DGVM generates a biome coverage that is heterogeneous, with more variability in the biome determination within a single region. For example, where Scandinavia is covered by forests, it is a mixture of boreal forests, both deciduous and evergreen. Conversely, in ORCHIDEE-DGVM and CARAIB, zones are much more homogeneous, with a single biome over wide areas. Another obvious difference is that SEIB-DGVM has zones with woodland at the mid latitudes whereas ORCHIDEE-DGVM has none and CARAIB have few in the high latitudes.

In the Alps, all models differ: CARAIB grows a boreal deciduous forest, ORCHIDEE-DGVM a steppe, SEIB-DGVM a boreal conifer forest.

In boreal regions, CARAIB grows very few conifer forests but clearly separated steppes and deciduous forest, SEIB-DGVM grows a mix of conifer and deciduous forest (dominated by conifer forest), and ORCHIDEE-DGVM grows separate conifer and deciduous forest, as well as steppes on the Atlantic coast of Norway.

Also, the limit between boreal and temperate biomes is shaped differently: SEIB-DGVM seems to have a limit determined by latitude, whereas ORCHIDEE-DGVM and CARAIB seems to have a northwest-southeast limit, resulting in regions like southern Sweden and Belarus, Estonia, Lithuania and Latvia with temperate biomes in ORCHIDEE-DGVM and CARAIB, but boreal in SEIB-DGVM. One probable reason for this is that boreal and temperate PFTs geographical repartition in SEIB-DGVM is less sensitive to temperature than in ORCHIDEE-DGVM and CARAIB and more to incoming solar radiation. Indeed, the shape of this limit in ORCHIDEE-DGVM and CARAIB seems similar to the average temperature gradient for a typical year (taken here as an example from the 6 k.a. dataset in Fig. 3) and might be the dominant factor as indicated by perturbation experiments (not shown) while the incoming radiation is the one forcing that has a clear latitudinal
<table><caption>**Table 5**
Distance PFTs ORCHIDEE/TERRANOVA.</caption><thead><tr><th></th><th>TeNEg</th><th>TeBEg</th><th>TeBSg</th><th>BNEg</th><th>BNSg</th><th>BBSg</th><th>Herbaceous</th></tr></thead><tbody><tr><th>TeNEg</th><td>0</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>3</td></tr><tr><th>TeBEg</th><td>1</td><td>0</td><td>2</td><td>3</td><td>3</td><td>3</td><td>3</td></tr><tr><th>TeBSg</th><td>1</td><td>2</td><td>0</td><td>2</td><td>2</td><td>2</td><td>3</td></tr><tr><th>BNEg</th><td>1</td><td>3</td><td>2</td><td>0</td><td>1</td><td>1</td><td>3</td></tr><tr><th>BNSg</th><td>2</td><td>3</td><td>2</td><td>1</td><td>0</td><td>1</td><td>3</td></tr><tr><th>BBSg</th><td>2</td><td>3</td><td>2</td><td>1</td><td>1</td><td>0</td><td>3</td></tr><tr><th>TeC3</th><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>0</td></tr><tr><th>C4</th><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>0</td></tr><tr><th>TrC3</th><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>0</td></tr><tr><th>BC3</th><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>0</td></tr></tbody></table>

<table><caption>**Table 6**
Parameters of each simulation run.</caption><thead><tr><th>Simulation name</th><th>Astronomical parameters</th><th>co2 concentration [pm]</th></tr></thead><tbody><tr><td>1900 A.D.</td><td>obliquity = 23.45, eccentricity = 0.01724, precession = 101.35</td><td>280</td></tr><tr><td>1 k.a. BP</td><td>obliquity = 23.45, eccentricity = 0.01724, precession = 101.35</td><td>280</td></tr><tr><td>3 k.a. BP</td><td>obliquity = 23.8205, eccentricity = 0.01838, precession = 50.30</td><td>275.14</td></tr><tr><td>4 k.a. BP</td><td>obliquity = 23.9294, eccentricity = 0.01870, precession = 33.45</td><td>273.32</td></tr><tr><td>6 k.a. BP</td><td>obliquity = 24.1054, eccentricity = 0.019755, precession = 359.99</td><td>263.98</td></tr><tr><td>8.5 k.a. BP</td><td>obliquity = 24.2221, eccentricity = 0.1976, precession = 318.57</td><td>261.78</td></tr></tbody></table>

<table><caption>**Table 7**
Climate fields for each DGVMs.</caption><thead><tr><th>Variable</th><th>Unit</th><th>SEIB</th><th>CARAIB</th><th>ORCHIDEE</th></tr></thead><tbody><tr><td>Surface Air Temperature</td><td>°C</td><td>X</td><td>X</td><td>X</td></tr><tr><td>Surface Air Temperature daily amplitude Minimum surface Air Temperature</td><td>°C</td><td>X</td><td>X</td><td>X</td></tr><tr><td>Maximum surface Air Temperature</td><td>°C</td><td></td><td></td><td>X</td></tr><tr><td>Minimum surface Air Temperature</td><td>°C</td><td></td><td></td><td>X</td></tr><tr><td>Precipitation</td><td>mm. day-1</td><td>X</td><td>X</td><td>X</td></tr><tr><td>1 - cloud cover</td><td>%</td><td></td><td>X</td><td></td></tr><tr><td>Downward shortwave radiation at midday</td><td>W.m-2</td><td>X</td><td></td><td></td></tr><tr><td>Downward longwave radiation</td><td>W.m-2</td><td>X</td><td></td><td></td></tr><tr><td>Wind speed</td><td>m.s-1</td><td>X</td><td>X</td><td>X</td></tr><tr><td>Relative humidity</td><td>%</td><td>X</td><td>X</td><td>X</td></tr></tbody></table>

component.

Now having a look at the differences across the different time windows in Fig. 2, we can clearly see that all models have very similar simulated biomes maps between 8.5 k.a. and 6 k.a. (with a clearer difference between those and 1900 A.D.), and that all models see Europe as being temperate deciduous forest dominated. Similarly, all models grow more steppes at 1900 A.D. than at 6 k.a. or 8.5 k.a. In particular, ORCHIDEE-DGVM grow more herbaceous in Scandinavia with time, with Norway and Sweden mainly covered by steppes at 1900 A.D., SEIB-DGVM grows more steppes and woodlands in central Europe and Iberian peninsula with time and CARAIB grow steppes in the north at all time scale and more in central Europe with time. SEIB-DGVM and CARAIB grow less temperate broad-leaved evergreen forest at 1900 A.D. than at 6 k.a. or 8.5 k.a., and ORCHIDEE-DGVM just a little bit less.

One big difference is the boreal deciduous forest, which is present in all models at all time scales, but with very different evolution with time: SEIB-DGVM grows the most boreal deciduous forest at 6 k.a. and more at 1900 A.D. than at 8.5 k.a., ORCHIDEE-DGVM and CARAIB less with

time.

### 4.1.2. Analysis of the differences

As biomes are computed from LAImax and dominant PFT, we hereafter investigate how these variables are simulated in the different models. In Fig. S18, the LAImax per number of point for each model, at each timescale is represented on histograms. Those LAImax values are very different: LAImax is much higher in ORCHIDEE-DGVM and in CARAIB than in SEIB-DGVM, which can explain why there is no woodland biome in ORCHIDEE-DGVM, and a few in.

CARAIB (located in areas surrounded by steppes: nordic regions and south-east of Turkey). Indeed, trees never have a low enough LAImax for that. In the case of CARAIB, approximately half of the points have a very high LAImax, higher than $5m^2/m^2$. The high variability in the dominant PFT in certain geographical regions in SEIB-DGVM questions the concept of dominance. In the following, we investigate how dominant a PFT is with respect to the others.

In Fig. 4, the number of PFTs which represent more than 10 % of the parameter used to determine the dominant PFT (the mean annual NPP for SEIB-DGVM and CARAIB or the mean annual maximum vegetation fraction for ORCHIDEE-DGVM) is represented on maps for each model and each time window, representing the number of PFTs sharing the same gridcell.

The three models give quite different results. CARAIB is the model which has the largest number of PFTs per grid cell, hence the highest PFT diversity, followed by ORCHIDEE-DGVM and then SEIB-DGVM. CARAIB has more diversity than the other two except in boreal regions, where between 1 and 3 PFTs are represented.

SEIB-DGVM has again a substantial variability of PFT distribution between adjacent cells in boreal regions, which is not the case in ORCHIDEE-DGVM and CARAIB. This high spatial heterogeneity can be attributable to its being an individual-based model. In opposite, SEIB-DGVM depict a very homogeneous distribution in temperate regions, with only 1 to 3 PFTs present in the same gridcell.

Whereas SEIB-DGVM and ORCHIDEE-DGVM have a the same range of PFT diversity in temperate regions and in boreal regions, CARAIB has, on the contrary a large diversity with between 3 and 5 PFTs growing in the same gridcell in temperate regions, but seems to have a very clear latitudinal border at around 60°N: beyond that limit to the north, there are only 1 or 2 PFTs growing. We can notice that, except for the 1900 A. D., SEIB-DGVM simulates a little bit more PFT diversity in boreal regions than temperate regions, which is not a common feature shared with CARAIB, which has a very poor diversity in Scandinavia, Finland and Russia.

Those results suggest that the concept of PFT dominance is largely model dependent. In SEIB-DGVM, in general one PFT dominate all the others whereas in CARAIB several PFTs can be present with close PFT fractions.

Another way to look at this relative dominance could be via the magnitude of NPP fraction that is represented by the first dominant PFT as shown in Fig. 5.

Here again, the three models yield very different results. CARAIB has
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Biomes maps for SEIB-DGVM (left) and ORCHIDEE-DGVM (middle) and CARAIB (right) for 8.5 k.a. (top), 6 k.a. (middle) and 1900 A.D. (bottom) climatic dataset.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Average temperature map for a typical year in the 6 k.a. climatic dataset.</figcaption></figure>

low overall values, consistent with a high number of dominant PFTs and most of the map has a NPP proportion of the first PFT of only 20–40 %. SEIB-DGVM has very few regions with less than 50 % of NPP dominance. If we consider that a true dominance is more than half the NPP or maximum vegetation fraction, this means that the dominant PFT is actually a dominant PFT for SEIB, most of the map in ORCHIDEE and only a few zones in CARAIB.

To have a better view on how mixed the type of plants are in a same grid cell, we further simplified the PFTs in three classes: broad-leaved forest, needle-leaved forest, and herbaceous. Fig. 6 shows the proportion of each of those three groups of PFTs for each pixel for SEIB-DGVM, ORCHIDEE-DGVM and CARAIB and for the three time periods.

Again, we can see that there is not much difference between 8.5 k.a. and 6 k.a. for all models. From 8.5 k.a. to 1900 A.D. period, ORCHIDEE-DGVM progressively gains needle-leaved trees (especially in Russia), while SEIB-DGVM has a quite constant distribution, again showing a less mixed forest per point than the other models. There are more

herbaceous in nordic countries at 1900 A.D. than at 6 k.a. and 8.5 k.a. in all model. CARAIB has herbaceous everywhere, unlike ORCHIDEE-DGVM and SEIB-DGVM which have them mainly in Scandinavia and the Alps, especially at 1900 A.D. Overall, needle-leaved are more abundant in SEIB-DGVM.

## 4.2. Model/data comparison

To evaluate how model simulations compare with data, we used the TERRANOVA pollen database. As a first step, we expressed the TERRANOVA dataset in PFTs fraction for each gridcells. Figs. 7–9 show respectively the results and analysis for 8.5 k.a., 6k.a., and 1900 A.D. On each figures, SEIB-DGVM is on the left, ORCHIDEE-DGVM on the middle, CARAIB on the right, pollen map reconstruction on the top, dominant PFTs map for each model on the middle, and distance map between pollen data and models on the bottom (see Tables 4 and 5).

The dominant PFT distribution for each of our three models is hereafter computed following the same methodology as before, using averages over the last ten years. To quantify the model data disagreement, we use the distance matrix introduced in 2.4.3, Tables 4 and 5. The first main global result from such a comparison is that the distance between model and data became higher with time for all models. This is coherent with the increase of land use activities (Zapolska et al., 2023a). Indeed, the pollen reconstruction maps contain more herbaceous PFTs at 1900 A.D. than at 6 k.a. and 8.5 k.a.

The 8.5 k.a. simulations are overall in good agreement with the data, except for Iberian Peninsula and south-east of Europe, dominated by herbaceous in the pollen data. ORCHIDEE-DGVM has a better score than SEIB-DGVM and CARAIB, and SEIB-DGVM doesn't represent the temperate PFT in southern Sweden, again probably because.

Of its distinction between boreal and temperate PFTs which may not be temperature dependent enough (unlike ORCHIDEE-DGVM and CARAIB).

At 6 k.a., the same result is found regarding the limit between boreal and temperate trees, again with ORCHIDEE-DGVM and CARAIB in better agreement than SEIB-DGVM. Also, there are a lot of scores with a
<figure><img src="image_5.png" /><figcaption>Fig. 4. Number of PFTs with a NPP greater than 10 % of the mean annual sum for all PFTs for SEIB-DGVM (left), ORCHIDEE-DGVM (middle) and CARAIB (right), at 8.5 k.a. (top), 6 k.a. (middle) and 1900 A.D. (bottom).</figcaption></figure>

value of one in the distance map (orange color), which correspond to needle-leaved versus broad-leaved disagreement. As seen before, this is due to the mixed forest in ORCHIDEE-DGVM. We consider that difference to be a rather small one. Globally, 6 k.a. results and distance maps are very similar to 8.5 k.a.

Concerning the Alps, all models are in disagreement with the pollen reconstructions: TERRANOVA indicates a dominance of temperate needle-leaved trees, whereas SEIB-DGVM gives boreal needle-leaved trees and the output of ORCHIDEE-DGVM and CARAIB yield a dominant herbaceous cover. The classification between temperate and boreal needle-leaved is somewhat arbitrary; taking that aspect into consideration, SEIB-DGVM is probably the model closest to the reconstructions in the alpine realm.

As expected, the 1900 A.D. map is not data compliant at all for all three models, because of the large proportion of herbaceous dominance in numerous parts of Europe (caused by crop-land and urbanisation), except for the boreal regions where pollen-based reconstruction show a broad forest dominance.

In order to have a better perspective at the time evolution of the data to model agreement, we introduce a generic matching score between the two as follow:

$$Ms = 100 - \frac{\text{distance}}{\text{maxdistance} \times 100}$$

where distance is represented by the values arising from the bottom panels in Figs. 7–9 and maxdistance is three.

In addition to 8.5 & 6 ka B.P. and 1900 A.D., we also performed similar computations for the 1 k.a., 3 k.a. and the 4 k.a. time windows, and reported all distance values as a matching score in Fig. 10. As already mentioned, the first order signal is a decrease in the matching score between model and data with time. At 8.5 k.a. and 6 k.a., the results are not significantly different in each of the three models. From 6 k.

a. BP onwards, all models agree on an accelerating decrease of the matching score towards the 1900 A.D. The difference between models is about 3–8 %. For the 1900 A.D. period, the results show a matching score of ≈25 – 28 %.

Since the three models are yielding quite similar and comparatively high matching results at 8.5 k.a. BP and 6 k.a. BP, we will as following step analyse where the models give a coherent spatial response for those two time windows. Fig. 11 shows the number of models having a distance to pollen data below or equal to one, which corresponds to a perfect agreement with the dominant PFT, or a small error, at 8.5 k.a. and 6 k.a.

This map defines robust areas for which all models agree and are consistent with pollen data. To a first order, the results are spatially coherent in temperate areas (western and central Europe) and most coastal areas but highlight strong divergences in the Iberian Peninsula, the Alps, Scandinavia and Eastern Europe.

Since the largest human impact within this time period is the expansion of agriculture, we expect the difference between trees and herbaceous PFTs to be more telling than the difference between different tree PFTs. It is therefore useful to further simplify our comparison using a simplified weighting matrix, only counting a distance of 1 between herbaceous and trees and 0 in other cases. It is interesting to see (Fig. 12) that for 6 k.a., the difference between the highest and the lowest score is still around 10 % as in Fig. 10 and that the three models have a reduced spread at 8.5 k.a. BP and 1900 A.D.

## 5. Discussion

In this study, we assumed that the DGVM results are relevant to represent potential natural vegetation. More precisely, the vegetation functioning part of DGVMs, which is the most mechanistic part of the models, has been well tested in previous studies and is probably more
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Percentage of total NPP of the dominant PFT for SEIB-DGVM (left), ORCHIDEE-DGVM (middle) and CARAIB (right) at 8.5 k.a. (top), 6 k.a. (middle) and 1900 A.D. (bottom).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Color map (red = herbaceous [%], green = broad-leaved trees [%], blue = needle-leaved trees, minimaps with only red, green or blue colors are showing the individual components and the larger map the combined values.) for SEIB-DGVM (left), ORCHIDEE-DGVM (middle) and CARAIB (right), for 8.5 k.a. (top), 6 k.a. (middle) and 1900 A.D. (bottom) climatic dataset.</figcaption></figure>

reliable than parts related directly to regeneration and survival. However, as those models are based on empirical relationships between vegetation and climate (Levavasseur et al., 2013), they do not take into account the role of non climatic factors such as megafauna on vegetation. It is important to note that megafauna may have had an impact on the opening up of environments, as suggested by Feurdean et al. (2018)

who studied actual big mammals and their role in land cover in central Europe and Zhu et al. (2018), who shows how taking herbivores into account modifies productivity in ORCHIDEE. That said, we assumed here that climate and human activities are the main driving factors of vegetation cover variability during the Holocene in Europe.

Another caveat is arising from the climate forcing used. The
<figure><img src="image_8.png" /><figcaption>Fig. 7. Dominant PFT map for pollen-based reconstruction (top), dominant PFT map for simulation (middle), and distance between model and data (bottom), for SEIB-DGVM (left), ORCHIDEE-DGVM (middle) and CARAIB (right), for 8.5 k.a.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Dominant PFT map for pollen-based reconstruction (top), dominant PFT map for simulation (middle), and distance between model and data (bottom), for SEIB-DGVM (left), ORCHIDEE-DGVM (middle) and CARAIB (right), for 6 k.a.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Dominant PFT map for pollen-based reconstruction (top), dominant PFT map for simulation (middle), and distance between model and data (bottom), for SEIB-DGVM (left), ORCHIDEE-DGVM (middle) and CARAIB (right) for 1900 A.D.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Data/model matching score as a function of time for SEIB (blue), ORCHIDEE (green), and CARAIB (red).</figcaption></figure>

Holocene period has an overall weak climate signal and is notoriously
difficult to consistently capture in climate models as has been shown
previously for example in Hargreaves et al. (2013). The iLOVECLIM
model results used here provide one climate forcing that is specific to
that model. While our approach correct the large scale model biases
through the bias correction, the climate anomalies between present-day
and the difference past periods are conserved and are characteristic of
the iLOVECLIM model simulations. These have been analyzed elsewhere
(Arthur et al., 2023). To further analyse the impact of model specific
simulations of climate change, the use of a range of climate models re-
sults could be used in a future study. A period such as the mid-Holocene
(Brierley et al., 2020) for which these simulations are available could be

a nice target to do so.

In terms of comparison with data, the DGVM simulations are statis-
tically similar for all models since they have overall similar percentage
of correlation to the pollen data. Matching score decrease with time,
which is expected as human pressure increase, consistently with what as
been shown in Zapolska et al. (2023a) on similar time period. The
matching score indicates that 72–75 % of the studied land cover differs
from the PNV at 1900 A.D. according to Fig. 10. Again, those results are
consistent with the results found in Zapolska et al. (2023a). Looking at
Fig. 10, even at 8.5 k.a, the maximum matching score is about 54 %
(SEIB-DGVM), which either means that the human pressure was already
important, or the models are not reliable at reconstructing the
<figure><img src="image_12.png" /><figcaption>Fig. 11. Number of models with a distance lower or equal to one with the pollen data for 8.5 k.a. (left), and 6 k.a. (right).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Data/model matching score as a function of time for SEIB (blue), ORCHIDEE (green), and CARAIB (red) for simplified pollen vs model distance.</figcaption></figure>

vegetation cover. Looking back at the spatial results, the distance between pollen reconstructions and SEIB-DGVM in Fig. 7 shows that the strongest divergences (score of 3) are located in Turkey and Spain – locations where agriculture was likely already present Gronenborn et al., 2023–, part of Ireland, Scandinavia and a few sporadic zones in temperate Europe. Those results could suggest that the presence of herbaceous taxons in Turkey and Spain have an anthropic origin, but ORCHIDEE-DGVM grows mainly herbaceous in Turkey as PNV. Looking at the TERRANOVA database, the anthropic origin of the herbaceous PFT is confirmed: herbaceous points in France, Ireland or Scandinavia are not dominated by human related herbaceous species (such as Castanea, Cerealia, Plantago lanceolata, Secale), and could be considered as natural, whereas Spain and Turkey are yet dominated by human related species. Please note that here, we did not count the taxon Poaceae, in either of the herbaceous category (natural or human-related), as Poaceae are a very big group containing both human-related species such as cerealia and natural species. Concerning Scandinavia, the climate-vegetation models do not account for the long-lasting impact of the ice cover, including soil development and succession of vegetation, since the models calculate the potential vegetation that is in equilibrium with the climate at that moment. In reality the vegetation could have lagged the climate response by up to a few thousand years and this lag effect can be expected to be present in the TERRANOVA pollen database but not in the models. So part of the difference between the models and pollen data could maybe be explained by this lag-effect (Moen and Lilthun, 1999). We can highlight that none of the DGVM models is able to reproduce the herbaceous dominance in the center of France (taken as a local example) and in Spain (Fig. 11). As discussed above, the presence of herbaceous dominance in Spain at 8.5 k.a. BP could be explained by

human land use. For the few points in the center of France, it could be due to the nature of the soil not being taken into account in the models - or to specific eco-environnements.

A way to understand if the decrease in data/model matching score can be attributed to human activities is to look at the land-use database (Hurtt et al., 2020). We used the C3 annual crops, the C3 nitrogen-fixing crops, the C3 perennial crops, the C4 annual crops, the C4 perennial crops, the urban areas and the managed pasture proportion in 1900 A.D. to determine a score of human land-use from 0 to 3 (0: 0–20 % of the land is used, 1: 20–40 % of the land is used, 2: 40–60 % of the land is used, 3: 60–100 % of the land is used) (see panel a) in Fig. 13). We then subtracted the score computed from the distance between models and pollen data in 1900 for the three models to it (panels b,c,d in Fig. 13). The difference in score obtained ranges from -3 to +3 and can be interpreted as follows: a value of 3 indicates that the.

Simulated PNV agrees with the pollen reconstructions (score of 0) and the LUH2-derived human-pressure is maximal (score of 3), thus that the human-modified vegetation is in the same class as the natural one (i. e. type of herbaceous cover); a value of 0 indicates the same score for both PNV-to-pollens reconstructions and human-pressure, thus that the difference between simulated PNV and pollen reconstruction is coherent with the human-pressure reconstructed in LUH-2; a value of -3 indicates that the score of PNV to pollen reconstructions is maximal (score of 3) while there is little to no human-pressure (score of 0) indicating that the large difference between simulated vegetation cover and reconstructed is probably not to be attributed to anthropogenic action but to a misfit between the model simulations and the natural reconstructed vegetation cover.

Fig. 13b shows this difference in score for SEIB-DGVM. Most of
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** (a) score of human land use from 0 to 3 calculated from the LUH2 dataset; (b) subtraction of the score of human land use and the distance between TERRANOVA pollen data and the SEIB-DGVM results; (c) same for ORCHIDEE-DGVM; (d) same for CARAIB.</figcaption></figure>

Europe has a near zero score, meaning that the difference between model and data can be attributed to human pressure. A big part of the UK and Ireland has a score of 3, because even if the land-use was important, it did not change the vegetation class which was herbaceous. Concerning northern Europe (Norway, Sweden and Finland) the value suggests a bad representation of the natural vegetation, as there was little land use. On a global scale, the same can be said for ORCHIDEE-DGVM on Fig. 13c and CARAIB on Fig. 13d.

However, if we look closely at northern Europe (Fig. 14), the LUH2 database shows that there were three types of land cover: non forested primary land, forested primary land, and forested secondary land. As

Fig. 14b shows, Norway and north-west Russia were mainly covered by primary forest, meaning that indeed, the ORCHIDEE-DGVM and.

CARAIB (Fig. 14d and e) did not succeed in simulating the PNV. For Sweden it is more complicated, because in that case the forest was man-made and we cannot be sure the PNV would have also been forest. Even so, this is a clear example of the limitation of our methodology, because this kind of land-use cannot be seen with the comparison score. SEIB-DGVM is the closest when it comes to comparing PNV to LUH2 reconstruction.

Nevertheless, Figs. 14a and b have some differences. According to LUH2, Norway is covered by forest, which is not the case in

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** (a) TERRANOVA pft dominant for northern Europe in 1900 A. D.; (b) color map (red = non forested primary land [%], green = forested primary land [%], blue = forested secondary land [%], minimaps with only red, green or blue colors are showing the individual components and the larger map the combined values.) according to the LUH2 dataset in 1900 A. D.; (c) biomes map obtained with SEIB-DGVM at 1900 A.D.; (d) biomes map obtained with ORCHIDEE-DGVM at 1900 A.D; (e) biomes map obtained with CARAIB at 1900 A.D.</figcaption></figure>
TERRANOVA, where some data suggest the presence of steppes. In the same way, LUH2 indicates that Finland is covered by non forested primary land, in other words, by natural steppes, whereas TERRANOVA indicates forests. This could be explained by the fact that LUH2 doesn't rely on pollen fossils but on scenarios of anthropogenic land-cover change (ALCCs) named HYDE. It is based on geographic information system (GIS) models and use historical data to estimate the size of the human population at a given time, as well as environmental information, such as altitude and proximity to sources of water (aridity index), in order to determine whether a location is suitable for human settlement or resource exploitation. Kaplan et al. (2017) found that HYDE underestimate land use when compared to REVEALS with increasing magnitude with time in the past, and that REVEALS estimates of open land fractions are realistic in magnitude and spatial distribution for the present-day time windows (0–100 BP), which is the case of our 1900 A.D study. They compared it to land cover observed with satellite remote sensing, combined with national land-use statistics on crop and pasture. Thus, where LUH2 indicates a strong human impact (and is therefore reliable), TERRANOVA indicates a dominance of herbaceous vegetation, and the model shows forest growth, we can conclude that there is human impact. However, when LUH2 indicates forest and TERRANOVA indicates herbaceous vegetation, it is not possible to draw any conclusions.

Indeed, in our study, the difference between TERRANOVA and LUH2 is even more noticeable at 1 k.a. BP.: according to LUH2, France, Benelux, Germany and Northern Italy are the only places where the human impact is significant (landuse superior to 20 % of the surface, see Fig. 15a), and Europe is mainly covered by natural steppes (in red in Fig. 15b). The pollen data for this time period (Fig. 15c) agree with the presence of herbaceous dominance in France, Spain, and Turkey, but indicate that the east of Europe is dominated by forest. On the opposite, LUH2 presents Ireland and.

Scotland as primary forested, whereas TERRANOVA presents them as herbaceous dominated. Those differences highlight the necessity to evaluate and improve ALCC scenarios in order to have a better representation of past landscape changes (Gaillard et al., 2018; Harrison et al., 2020). However, the comparison with our DGVM results at 1900 A.D. is still valuable, as this time period seems more accurate and that we are looking at the first order signal. LUH2 and pollen data agree on European landscape largely modified by human activities except for northern countries where natural vegetation was still significantly present. That said, ORCHIDEE-DGVM and CARAIB didn't grow forest but mainly steppes at those latitudes, meaning SEIB-DGVM is better at representing nordic PNV.

The three models outputs are very different. SEIB-DGVM have a very low diversity of PFTs, whereas CARAIB and ORCHIDEE-DGVM have a higher diversity. The high diversity of PFTs of CARAIB and the low importance of its dominant PFT could be explained by its two layers structure, leading to an absence of competition between the herbs and the trees. There will be only herbs or trees if the other cannot grow.

In Figs. 4 and 5, the three models differ a lot concerning both the

number of PFTs per grid cells and the relative importance of the dominant PFT. The number of PFTs is not necessarily an advantage or disadvantage; it is by.

Comparing them to pollen data that we can conclude whether we are observing, for example, a mixed forest or one representative of a single PFT. Like for Figs. 4 and 5, in Fig. 16 we show the number of PFTs and percentage of the dominant one but estimated from the pollen dataset. Indeed, all the taxons are given as a percentage, so as the reclassified PFTs we created. Those results suggest that the pollen-based reconstruction favors a dominance of between 1 and 4 PFTs per grid cell, mainly 2 or 3. For this particular aspect, SEIB-DGVM is in better agreement with the data than CARAIB and above all ORCHIDEE-DGVM, which grows many more. However, ORCHIDEE-DGVM is better at reproducing the relative importance (in terms of NPP) of a dominant PFT. Given those results, we can notice that the best model on an aspect is not necessarily the best for all the studied parameters; using two or more models can thus add confidence into the simulated fields.

On a regional basis, pollen data suggest a higher number of PFTs in boreal zones, a tendency not reproduced by SEIB-DGVM and ORCHIDEE-DGVM; CARAIB shows exactly the opposite. CARAIB is also the only model which does not correspond to either the number of PFTs, the geographical distribution, and the importance of the dominant PFT. From this perspective, CARAIB is more of an outlier.

One classical framework to look at vegetation changes over time, especially in relation to the human occupation (Nikulina et al., 2024) is to express the vegetation in terms of openness.

Fig. 17 represents the percentage of herbaceous PFT in the TERRANOVA dataset, taken here as an approximation for vegetation openness. The mean value in Europe is 23.7 % at 8.5 k.a. (which is the period with the lowest openness), and 40.2 % for the 1900 A.D. The only model approaching the 8.5 k.a. openness is ORCHIDEE-DGVM. The vegetation fraction of CARAIB cannot be used, as it is half tree half herbaceous covered everywhere, the openness is always equal to 50 %, as a consequence of its two layers structure.

8.5 k.a. and 6 k.a. DGVM simulations are very similar for all models, which should be excepted, as the climate was very similar in those period in the northern hemisphere (Borgatti and Soldati, 2013). Temperate zones seems to be well represented by all models but, models show big differences in boreal zones. In particular, CARAIB shows almost no trees northward from a latitude of about 60°N.

Another aspect to consider is the differences in representation of the PFTs which greatly condition the differences between models. Indeed, in ORCHIDEE-DGVM, the PFTs are represented by a surface covered by each and spatially separated (at least partially because in reality, in the DGVM there is a competition for light which therefore supposes an overlap) and by internally representing the functioning of an average plant. This is why the NPP of ORCHIDEE-DGVM needs to be multiplied by the maximum vegetation fraction to be compared with the NPP of SEIB-DGVM and CARAIB. On the other hand, SEIB-DGVM has a totally different way of representing vegetation since it does not have explicit

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** (a) score of human land use from 0 to 3 calculated from the LUH2 dataset; (b) color map (red = non forested primary land [%], green = forested primary land [%], blue = human related landscape, minimaps with only red, green or blue colors are showing the individual components and the larger map the combined values.); (c) dominant PFT of the TERRANOVA pollen data for 1 k.a. BP.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Number of PFTs with a fraction of the cells greater than 10 % in the TERRANOVA dataset (top), and percentage occupied by the dominant PFT in the TERRANOVA dataset (bottom) at 8.5 k.a. (left), 6 k.a. (middle) and 1900 A.D. (right).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Openness of pollen data at 8.5 k.a. (left), 6 k.a. (middle), and 1900 A.D. (right).</figcaption></figure>

an VEGETMAX the PFTs being intermingled. For CARAIB, there is no VEGETMAX either, and there is also as we saw a two-layer structure: one with herbs and shrubs, and one with trees.

## 6. Summary and perspectives

In our study, we ran three different DGVMs forced by simulated climates for six different time-slices, from 8.5 k.a. B.P to 1900 A.D., and compared their results to reclassified pollen-based reconstructions. The three models have statistically similar results compared to the pollen-based reconstructions, but they also differ a lot on certain regions. The Alps and Scandinavia are poorly described by the models, suggesting a problem in the description of cold and mountainous regions. However, SEIB-DGVM seems to be the closest to pollen data in northern Europe, which can be considered as PNV according to LUH2.

All those results suggest that SEIB-DGVM and ORCHIDEE-DGVM are complementary and that there is a benefit of them being used together in order to have a better description of the natural vegetation cover of the Holocene. Overall, our study suggests that CARAIB performs less well for simulating the PNV over the Holocene, as, by and large, it reproduces neither the PFT diversity nor the PFTs dominance geographical distribution. An interesting perspective would be to make the same study under the first part of the Holocene (from 12 k.a. BP to 9 k.a. BP), to see if the pathway is the same for the matching score between pollen-based reconstructions and DGVM models.

To improve on the latter, an objective methodology to build a continuous spatial reconstruction from the TERRANOVA database will be the subject of a future study.

With such differences between DGVM model results, we can ask the

impact such differences could have on climate models. Indeed, results such as snow layer, leading to a change in albedo, and canopy height, impacting the wind, could be different depending on the DGVM used in climate model. On further interesting aspect could be to work on adapting SEIB-DGVM and/or ORCHIDEE-DGVM to run online within the same climatic model to assess the importance of climate vegetation coupling.

The already mentioned impact of the climate forcing could be also worth investigating by targeting one particular time window such as the mid-Holocene (similar to PMIP) and perform a large ensemble of multi-climate/multi-vegetation model combinations to assess in full the uncertainty arising from the modeling components on the simulated vegetation cover.

## Credit author statement

All authors contributed to the study design and to the writing of the manuscript. Hisashi Sato additionally provided support regarding the SEIB-DGVM model, Nicolas Viovy ran the ORCHIDEE-DGVM simulations and Didier M. Roche ran the iLOVECLIM simulations. All analysis and further technical setup was done by Isabeau Bertrix. All authors contributed to the writing and approved the final version of the manuscript.

## Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors in addition to base support of all co-authors from their respective institutions.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment

The authors would like to thank Louis M. François for providing the CARAIB global dynamic vegetation model, Maria Fernanda Sanchez Goñi for her help in the classification of the REVEALS pollen taxa into PFTs, Didier Paillard for his support on the insolation questions, Yoann Robin for his help in running the CDF-t bias correction with the SBCK package and Jean-Yves Peterschmitt for general coding support.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.quascirev.2026.109813.

Data/model availability

The REVEALS dataset is available from Fyfe et al. (2021). The software for the CDF-t bias correction method is available at https://github.com/yrobink/SBCK. The SEIB-DGVM model (version 3.03) was downloaded from: https://seib-dgvm.com/en/. The CARAIB model and the version of the ORCHIDEE model used here are not available publicly.

References

Arthur, F., Roche, D.M., Fyfe, R., Quiquet, A., Renssen, H., 2023. Simulations of the Holocene climate in Europe using an interactive downscaling within the iLOVECLIM model (version 1.1). Clim. Past 19 (1), 87-106.

Berger, A., Loutre, M.-F., 1999. Parameters of the earths orbit for the last 5 million years in 1 kyr resolution [supplement to: Berger, A; Loutre, M-F (1991): insolation values for the climate of the last 10 million of years. Quat. Sci. Rev. 10 (4), 297-317.
https://doi.org/10.1016/0277-3791(91)90033-Q.

Borgatti, L., Soldati, M., 2013. Hillslope processes and climate change. Treatise on Geomorphology 7, 306-319.

Brierley, C.M., Zhao, A., Harrison, S.P., Braconnot, P., Williams, C.J.R., Thornalley, D.J. R., Shi, X., Peterschmitt, J.-Y., Ohgaito, R., Kaufman, D.S., Kageyama, M., Hargreaves, J.C., Erb, M.P., Emile-Geay, J., D'Agostino, R., Chandan, D., Carré, M., Bartlein, P.J., Zheng, W., et al., 2020. Large-scale features and evaluation of the PMIP4-CMIP6 midHolocene simulations. Clim. Past 16 (5), 1847-1872.

Caley, T., Roche, D.M., 2013. δ18O water isotope in the iLOVECLIM model (version 1.0) – part 3: a palaeo- perspective based on present-day data-model comparison for oxygen stable isotopes in carbonates. Geosci. Model Dev. (GMD) 6 (5), 1505–1516.

Chini, L., Hurtt, G., Froliking, S., 2014. LUH1: Harmonized Global Land Use for Years 1500-2100, p. VI.

Dallmeyer, A., Claussen, M., Brovkin, V., 2019. Harmonising plant functional type distributions for evaluating Earth system models. Clim. Past 15 (1), 335-366.

Dallmeyer, A., Kleinen, T., Weitzel, N., Cornou, M.E., 2022. The Deglacial Forest Conundrum.

Dallmeyer, A., Poska, A., Marquer, L., Seim, A., Gaillard, M.-J., 2023. The challenge of comparing pollen-based quantitative vegetation reconstructions with outputs from vegetation models – a European perspective. Clim. Past 19 (7), 1531-1557.

Feurdean, A., Ruprecht, E., Molnár, Z., Hutchinson, S.M., Hickler, T., 2018. Biodiversity-rich European grasslands: ancient, forgotten ecosystems. Biol. Conserv. 228, 224-232.

Fisher, R.A., Koven, C.D., Anderegg, W.R.L., Christoffersen, B.O., Dietze, M.C., Farrier, C. E., Holm, J.A., Hurtt, G.C., Knox, R.G., Lawrence, P.J., Lichstein, J.W., Longo, M., Matheny, A.M., Medvigy, D., Muller-Landau, H.C., Powell, T.L., Serbin, S.P., Sato, H., Shuman, J.K., et al., 2018. Vegetation demographics in Earth system models: a review of progress and priorities. Glob. Change Biol. 24 (1), 35-54.

François, L., Utescher, T., Favre, E., Henrot, A.-J., Warnant, P., Micheals, A., Erdei, B., Suc, J.-P., Cheddadi, R., Mosbrugger, V., 2011. Modelling late Miocene vegetation in Europe: results of the CARAIB model and comparison with palaeovegetation data. Palaeogeogr. Palaeoclimatol. Palaeoecol. 304 (3), 359-378.

Fyfe, R.M., Githumbi, E., Trondmann, A.-K., Mazier, F., Nielsen, A.B., Poska, A., Sugita, S., Woodbridge, J., contributors, LandClimII, Gaillard, M.-J., 2021. A Full Holocene Record of Transient Gridded Vegetation Cover in Europe.

Gaillard, M.-J., Morrison, K.D., Madella, M., Whitehouse, N., 2018. Past land-use and land-cover Change: The Challenge of Quantification at the Subcontinental to Global Scales.

Githumbi, E., Fyfe, R., Gaillard, M.-J., Trondman, A.-K., Mazier, F., Nielsen, A.-B., Poska, A., Sugita, S., Woodbridge, J., Azuara, J., Feurdean, A., Grindean, R., Lebreton, V., Marquer, L., Nebout-Combourieu, N., Stačikaitė, M., Tačtāu, I.,

Tonkov, S., Shumilovskikh, L., data contributors, L., 2022. European pollen-based REVEALS land-cover reconstructions for the Holocene: methodology, mapping and potentials. Earth Syst. Sci. Data 14 (4), 1581-1619.

Goosse, H., Brovkin, V., Fichefet, T., Haarsma, R., Huybrechts, P., Jongma, J., Mouchet, A., Selten, F., Barriat, P.-Y., Campin, J.-M., Deleersnijder, E., Driesschaert, E., Goelzer, H., Janssens, I., Loutre, M.-F., Morales Maqueda, M.-A., Opsteegh, T., Mathieu, P.-P., Munhoven, G., et al., 2010. Description of the Earth system model of intermediate complexity LOVECLIM version 1.2. Geosci. Model Dev. (GMD) 3 (2), 603–633.

Gronenborn, D., Horejs, B., Börner, M., Ober, M., 2023. Expansion of Farming in Western Eurasia, 9600 - 4000 Cal BC (Update Vers. 2023.1).

Hargreaves, J.C., Annan, J.D., Ohgaito, R., Paul, A., Abe-Ouchi, A., 2013. Skill and reliability of climate model ensembles at the last glacial maximum and mid-Holocene. Clim. Past 9 (2), 811–823.

Harrison, S.P., Gaillard, M.-J., Stocker, B.D., Vander Linden, M., Klein Goldewijk, K., Boles, O., Braconnot, P., Dawson, A., Fluet-Chouinard, E., Kaplan, J.O., Kastner, T., Pausata, F.S.R., Robinson, E., Whitehouse, N.J., Madella, M., Morrison, K.D., 2020. Development and testing scenarios for implementing land use and land cover changes during the Holocene in Earth system model experiments. Geosci. Model Dev. (GMD) 13 (2), 805–824.

Haxeltine, A., Prentice, I.C., 1996. BIOME3: an equilibrium terrestrial biosphere model based on ecophysiological constraints, resource availability, and competition among plant functional types. Glob. Biogeochem. Cycles 10 (4), 693-709.

Hopcroft, P.O., Valdes, P.J., Harper, A.B., Beerling, D.J., 2017. Multi vegetation model evaluation of the green sahara climate regime. Geophys. Res. Lett. 44 (13), 6804-6813.

Hurtt, G.C., Chini, L., Sahajpal, R., Frolking, S., Bodirsky, B.L., Calvin, K., Doelman, J.C., Fisk, J., Fujimori, S., Klein Goldewijk, K., Hasegawa, T., Havlik, P., Heinimann, A., Humpenöder, F., Jungclaus, J., Kaplan, J.O., Kennedy, J., Krisztin, T., Lawrence, D., et al., 2020. Harmonization of global land use change and management for the period 850-2100 (LUH2) for CMIP6. Geosci. Model Dev. (GMD) 13 (11), 5425-5464.

Kaplan, J.O., Krumhardt, K.M., Gaillard, M.-J., Sugita, S., Trondman, A.-K., Fyfe, R., Marquer, L., Mazier, F., Nielsen, A.B., 2017. Constraining the deforestation history of Europe: evaluation of historical land use scenarios with pollen-based land cover reconstructions. Land 6.

Lange, S., EartH2Observe, WFDEI and ERA-Interim data Merged and Bias-corrected for ISIMIP (EWEMBI). V. 1.1. GFZ Data Services. https://dataservices.gfz-potsdam.de/pik/showshort.php?id=escidoc:1809891.

Levavasseur, G., Vrac, M., Roche, D.M., Paillard, D., Guiot, J., 2013. An objective methodology for potential vegetation reconstruction constrained by climate. Global Planet. Change 104, 7-22.

Li, H., Renssen, H., Roche, D.M., 2019. Global vegetation distribution driving factors in two dynamic global Vegetation Models of contrasting complexities. Global Planet. Change 180, 51-65.

Lu, Z., Miller, P.A., Zhang, Q., Zhang, Q., Wärlind, D., Nieradzik, L., Sjolte, J., Smith, B., 2018. Dynamic vegetation simulations of the Mid-Holocene green sahara. Geophysical Research Letters 45 (16), 8294-8303.

Moen, A., Lillethun, A., 1999. National Atlas of Norway : Vegetation. Norwegian Mapping Authority.

Nikulina, A., MacDonald, K., Zapsolska, A., Serge, M.A., Roche, D.M., Mazier, F., Davoli, M., Svenning, J.-C., van Wees, D., Pearce, E.A., Fyfe, R., Roebroeks, W., Scherjon, F., 2024. Hunter-gatherer impact on European interglacial vegetation: a modelling approach. Quat. Sci. Rev. 324, 108439.

Petit, J.-R., Raynaud, D., 2020. Forty years of ice-core records of CO2. Nature 579 (7800), 505-506.

Quiquet, A., Roche, D.M., Dumas, C., Paillard, D., 2018. Online dynamical downscaling of temperature and precipitation within the iLOVECLIM model (version 1.1). Geosci. Model Dev. (GMD) 11 (1), 453-466.

Robakowski, P., Bielinis, E., 2017. Needle age dependence of photosynthesis along a light gradient within an Abies alba crown. Acta Physiol. Plant. 39 (3).

Sato, H., Itoh, A., Kohyama, T., 2007. SEIB-DGVM: a new Dynamic Global vegetation model using a spatially explicit individual-based approach. Ecol. Model. 200 (3), 279-307.

Sato, H., Shibuya, M., Hiura, T., 2023. Reconstructing spatiotemporal dynamics of mixed conifer and broad- leaved forests with a spatially explicit individual-based dynamic vegetation model. Ecol 38 (3), 465-478.

Serge, M., Fyfe, R., Gaillard, M.-J., Klein, T., Lagnoux, A., Galop, D., Githumbi, E., Mindrescu, M., Nielsen, A.-B., Trondman, A.-K., Poska, A., Sugita, S., Woodbridge, J., 2023a. Spatially extensive and temporally continuous three REVEALS pollen-based vegetation reconstructions. In: Europe over the Holocene.

Serge, M.-A., Mazier, F., Fyfe, R., Gaillard, M.-J., Klein, T., Lagnoux, A., Galop, D., Githumbi, E., Mindrescu, M., Nielsen, A., Trondman, A.-K., Poska, A., Sugita, S., Woodbridge, J., Abel-Schaad, D., Åkesson, C., Alenius, T., Ammann, B., Andersen, S., et al., 2023b. Testing the effect of relative pollen productivity on the REVEALS model: a validated reconstruction of Europe-Wide Holocene vegetation. Land 12 (5).

Strandberg, G., Lindström, J., Poska, A., Zhang, Q., Fyfe, R., Githumbi, E., Kjellström, E., Mazier, F., Nielsen, A.B., Sugita, S., Trondman, A.-K., Woodbridge, J., Gaillard, M.-J., 2022. Mid-Holocene European climate revisited: new high-resolution regional climate model simulations using pollen-based land-cover. Quat. Sci. Rev. 281, 107431.

Tallavaara, M., Luoto, M., Korhonen, N., Järvinen, H., Seppä, H., 2015. Human population dynamics in Europe over the Last glacial maximum. Proc. Natl. Acad. Sci. U. S. A 112 (27), 8232-8237.

Tong, S., Wang, W., Chen, X., Jie, Chong-Yu, Sato, H., Wang, G., 2022. Impact of changes in climate and CO2 on the carbon storage potential of vegetation under limited water
availability using SEIB-DGVM version 3.02. Geosci. Model Dev. (GMD) 15 (18), 7075–7098.

Zapolska, A., Vrac, M., Quiquet, A., Extier, T., Arthur, F., Renssen, H., Roche, D.M., 2023a. Improving biome and climate modelling for a set of past climate conditions: evaluating bias correction using the CDF-t approach. Environ. Res.: Climate 2 (2), 025004.

Zapolska, A., Serge, M.A., Mazier, F., Quiquet, A., Renssen, H., Vrac, M., Fyfe, R., Roche, D.M., 2023b. More than agriculture: analysing time-cumulative human impact on European land-cover of second half of. Holocene 314, 108227.

Quaternary Science Reviews, 314, 108227, Zhu, D., Ciais, P., Chang, J., Krinner, G., Peng, S., Viovy, N., Peñuelas, J., Zimov, S., 2018. The large mean body size of mammalian herbivores explains the productivity paradox during the last Glacial maximum. Nat. Ecol. Evol. 640–649.