Did the Bronze Age deforestation of Europe affect its climate? A regional climate-model study using pollen-based land-cover reconstructions

. This paper studies the impact of land use and land-cover change (LULCC) on the climate around 2 500 years ago (2.5 ka), a period of rapid transitions across the European landscape. One global climate model was used to force two regional climate models (RCMs). The RCMs used two land cover descriptions. The first was from a dynamical vegetation model representing potential land cover; the second a land cover description reconstructed from pollen data by statistical interpolation. 20 The two different land covers enable us to study the impact of land cover on climate conditions. Since the difference in landscape openness between potential and reconstructed land cover is mostly due to LULCC, this can be taken as a measure of early anthropogenic effects on climate. Since the sensitivity to LULCC is dependent on the choice of climate model, we also use two RCMs. The results show that the simulated 2.5 ka climate was warmer than the simulated pre-industrial (PI) climate. The largest 25 differences are seen in northern Europe, where the 2


Introduction
Increases in the concentration of greenhouse gases in the atmosphere have an impact on the global climate through altering Earth's radiative balance (e.g.Manabe and Wetherhald, 1967).Ruddiman (2003) suggested that humans may have started to affect Earth's climate long before industrialisation and the modern era through human-induced emissions of greenhouse gases, particularly carbon dioxide (CO 2 ) and methane (CH 4 ), since the beginning of agriculture (Neolithic cultures).The "Ruddiman anthropogenic hypothesis" challenged the two "natural hypotheses" (i.e.internal processes of the climate system) proposed by Broecker et al. (1999) and Ridgwell et al. (2003), in an attempt to explain the 20 ppm CO 2 increase during the last 7000 years of the current interglacial (the Holocene) that maintained a warm climate, in contrast to the climate history of earlier interglacial periods.The anthropogenic hypothesis has been tested in a large number of studies over the last 20 years, using scenarios of Holocene land-use and land-cover change (LULCC) (e.g.Kaplan et al., 2011), different population growth data (e.g.Boyle et al., 2011) and/or 14 C dated archaeological data (e.g.Ruddiman et al., 2016).Few tests have used pollen-based reconstructions of plant cover, although Ruddiman et al. (2016) used a pollen-based reconstruction of vegetation "pseudobiomes" (anthropogenic biomes) sensu Fyfe et al. (2015) to provide evidence of sub-continental scale transformation of the vegetation of Europe.Most of the tests have looked primarily on the impact of anthropogenic emissions of CO 2 and CH 4, and none of them have yet been able to reject the anthropogenic hypothesis unequivocally (e.g.see discussions in Stocker et al., 2011 andRuddiman et al., 2011), and Ruddiman et al. (2016) have added a number of convincing arguments in supporting it.
Both anthropogenic and potential land-cover changes can also influence biogeophysical effects by altering the land surface albedo and roughness (e.g.Field et al., 2007).These effects can have large impacts on regional climate (e.g.Jia et al., 2019), which will affect people directly on a regional to local scale.Such impacts were explicitly shown for Europe in a number of studies with high-resolution regional climate models (RCMs) under present-day conditions (e.g.Belušić et al., 2019;Strandberg and Kjellström, 2019;Davin et al., 2020).There is also evidence of anthropogenic impact on climate via land-use changes during the Roman Period (Dermody et al., 2012;Gilgen et al., 2019).For the Mid-Holocene (6000 calibrated years before the present, henceforth abbreviated 6 ka) two studies have investigated the effect of LULCC on the climate of Europe.
The first used two different LULCC scenarios (Strandberg et al., 2014) and the second used pollen-based reconstructions of past quantitative land-cover for the first time (Strandberg et al., 2022).Results from these RCM studies indicate that, in southern Europe where an increase in open land cover (herb and low shrub vegetation) at 6 ka can be in large part assigned to human activities, LULCC can have contributed to a change in the average summer temperature of around 0.5 C. The authors also show that the results are partly model dependent.The response in summer temperature is ultimately decided by the albedo and how the heat fluxes are partitioned between sensible and latent heat.The temperature differences are determined by the model's radiation schemes and land surface physical processes, e.g. the amount of soil moisture available for evapotranspiration.
Although anthropogenic land use by Neolithic cultures impacted the land cover in parts of Europe from ca. 7 ka (e.g.Marquer et al., 2014;2017), and may have altered the regional climate in various ways (Strandberg et al., 2014;2022), major episodes of deforestation caused by human impact did not occur before 3000 years later.The most recent pollen-based reconstruction of land cover in Europe (spatial resolution of 1 ; Githumbi et al., 2022) suggests that the earliest of the two major deforestation episodes before the start of Modern time (1500 CE (0.45 ka) -present) took place between ca. 4 ka and 2.5 ka, i.e. the period during which the Bronze Age culture expanded from South-East (Turkey, Greece) to Central and Western Europe (Mediterranean area included) and Northern Europe (Champion et al., 1994;Coles and Harding, 1979).The second deforestation episode (before the Modern time deforestation) occurred ca.0.9 ka -0.5 ka, during the Middle Ages (ca.500 CE (1.45 ka) -1500 CE in most Europe, started 1050 CE (0.9 ka) in northern Europe)) (Fig. 1a).The difference in open-land cover between 4 ka and 2.5 ka of ca.10% (either in mean or median cover; Fig. 1a) is assumed to represent deforestation of Europe by Bronze Age cultures.This change in the land cover of Europe was also explained by deforestation for agriculture in the study of Marquer et al. (2017).If we consider a mean natural open-land cover of ca.20% for the whole of Europe and an increase in open-land cover of a maximum ca. 5% from 7 ka to 4 ka (Fig. 1a), the Bronze Age deforestation corresponds to an increase of open-land cover by 200% since 4 ka.At a broad regional scale, the mean increase between 4 ka and 2.5 ka is also ca. 10 % in each of Europe three major biomes (boreal, temperate and Mediterranean Europe; Fig. 1b-d).At the finer spatial scale (the individual 1 grid cells), the increase ranges from low (or no change) to sometimes much greater conversion of forest to farmland (Githumbi et al., 2022).The period around 6 ka is used by palaeoecologists and climate modellers as a representative time slice of the Holocene thermal maximum and the start of agriculture in large parts of the world.It is used broadly in model-data comparison to evaluate climate models against palaeoclimate proxy data (e.g.Braconnot et al. 2012), and to understand and explain the 'Holocene temperature conundrum' (e.g.Bader et al., 2020).However, a better understanding of the possible impact of land use on the Holocene climate since the start of agriculture requires that times of major change in land use are studied.The time around 3 ka (the Bronze Age) was also pinpointed as the time when "the planet [was] largely transformed by hunter-gatherers, farmers, and pastoralists", as suggested by an archaeological global assessment of land use from 10 ka to 1850 CE (ArchaeoGLOBE Project, 2019).Nevertheless, this time period has not been studied specifically by the wider palaeoclimate modelling community so far.In this paper, we are testing the following hypothesis: the regional climate of Europe around 2.5 ka was significantly influenced by land use via biogeophysical effects due to changes in the properties of the land surface.From a more general perspective, it is also important to understand the biogeophysical effects of LULCC, given that LULCC is affecting most of the Earth's land surface today (Parmesan et al., 2022).Although landscape management is expected to play a large role in the mitigation of climate warming in the future (Nabuurs et al., 2022), the local impact of such mitigation measures on the regional climate and the net effect of both biogeochemical and biogeophysical effects, in the past as well as today, are still not fully understood/known (e.g.Gaillard et al., 2015).The relative importance of LULCC as a climate forcing factor may be large, particularly in low emissions scenarios (van Vuuren et al., 2011).

Model chain
This study builds upon a chain of model simulations (see detailed model descriptions below).First, 2.5 ka and latest preindustrial (1850 CE, hereafter PI) climate conditions are simulated by the ocean-atmosphere coupled global model EC-Earth using prescribed potential PI land cover (Fig. 2 and Table 1).In this context, 'potential' refers to land cover without human intervention.The climate conditions for 2.5 ka are then used to force the DVM LPJ-GUESS over the European domain to simulate the potential land cover consistent with the climate, labelled L (2.5k-L in Table 1).The climate conditions from EC-Earth and the land cover from LPJ-GUESS are used to force the RCMs RCA4 and HCLIM.The RCMs simulate the same periods as the GCM, but at higher horizontal resolution and with their physical parameterisations.For each RCM, the model output includes a high-resolution climate simulation for 2.5 ka and PI respectively (2.5k-L and PI in Table 1).Land cover (i.e.plant cover) at 2.5 ka was reconstructed at a 1 spatial scale using >1000 pollen records and the REVEALS model (Sugita, 2007;Githumbi et al., 2022).Because pollen records are not available in all grid cells, a spatial statistics method was used to interpolate the pollen-based land cover over the entire grid covering Europe using the approach described in e.g.Pirzamenbein et al. (2018).The method uses additional co-variates, i.e. simulated potential land cover from LPJ-GUESS (driven by the climate from the EC-Earth simulation, Fig. 1) and the KK10 LULCC scenario for 2.5 ka (Kaplan et al., 2009).This spatially continuous land-cover reconstruction is labelled R (2.5k-R in Table 1) and represents the climate-and humaninduced land cover.
A similar approach was used in Strandberg et al. (2022).Here we use a more straightforward approach where the LPJ-GUESS simulations use the climate from EC-Earth only, instead of also doing an iteration where the RCMs force LPJ-GUESS.
Although some studies indicate that more than one iteration between a climate model and the vegetation model are needed to reach equilibrium (e.g.Velasquez et al., 2021), and simulated land cover is sensitive to the representation of climate used and vice versa (Strandberg et al., 2022), we chose to use only one iteration.This is because the difference between reconstructed (based on pollen data) and simulated (vegetation-model based) land cover is larger than the differences between alternative vegetation-model simulations of land cover (Strandberg et al., 2022).Thus, differences in potential vegetation (L) between RCM run iterations would probably not have a large impact on the response to "actual" land cover (R).That response is to a large degree determined by the scale of land-cover changes and the climate model used (Davin et al., 2020;Strandberg et al., 2022).Contrasting the simulated potential land cover with the pollen-based reconstruction allows us to explore the effect of LULCC on regional climate.The only difference between the regional climate simulations is the representation of land cover.This would not be possible to achieve with coupled DVM-RCM simulations.Therefore, we catch some of the uncertainty by using the two different land-cover representations L and R (Fig. 2).Such sensitivity tests give us the possibility to study the effect of land cover on climate, and how two different RCMs respond to land-cover change.The approach also allows to produce and evaluate a multi-model estimate of 2.5 ka climate.We perform the PI and 2.5 ka simulations following the PMIP4 protocol (Otto-Bliesner et al., 2017).The model setup is described in Zhang et al. (2021).For the 2.5 ka simulation, the climate forcing factors that differ from PI are orbital parameters, and concentrations of carbon dioxide (CO 2 ) and methane (CH 4 ).The orbital forcing for PI and 2.5 ka is calculated following Berger (1978) are the same in PI and 2.5 ka.The land cover used in the PI and 2.5 ka simulations is prescribed using simulated potential land cover for 1850 CE from LPJ-GUESS implemented offline.The quasi-equilibrium representing the 2.5 ka climate conditions is reached after 200 years of simulation.Note that this LPJ-GUESS simulation is not part of the model runs performed specifically for this study (for more details, see Zhang et al., 2021).Because the choice of vegetation plays a role in the resulting simulated PI (1850 CE) and 2.5ka climates, it will in turn have an impact on the comparison between PI and 2.5 ka climates.
However, the expected climate response of using 'actual' rather than potential PI vegetation in the EC-Earth simulations is much smaller than the difference in climate between PI and 2.5 ka.Strandberg et al. (2014) show a response to anthropogenic land cover (instead of potential land cover) at 0.2 ka of maximum +0.5 C. Therefore, we assumed that the choice of PI vegetation (either potential or 'actual') in the EC-Earth simulations should not affect the simulated 2.5 ka vegetation or climate to such a degree that it would significantly influence the difference in climate between 2.5 ka and PI.Moreover, given that the major focus of our paper is on the difference in the 2.5 ka simulated regional climate when either potential vegetation or 'actual pollen-based vegetation' is used in the RCMs runs, we considered that the use of the same potential land-cover in the EC-Earth climate runs for both PI and 2.5 ka would have little impact on the boundary conditions used for the RCMs runs at 2.5 ka BP.

RCA4 and HCLIM
Two regional climate models (RCM) are used to downscale the simulations from EC-Earth, thereby adding geographical details and improving the simulation of climatic processes as the horizontal grid spacing is denser in the applied RCMs (~50 km) than in EC-Earth (~100 km).In this study, the Rossby Centre Atmosphere model (RCA4, Strandberg et al., 2015;Kjellström et al., 2016) and the HCLIM38-ALADIN (HCLIM, Belušić et al., 2020) are run on a horizontal grid spacing of 0.44 (corresponding to approximately 50 km) across Europe (the CORDEX EUR-44 domain (Jacob et al., 2014)).They have different model physics and can not be expected to respond to LULCC in the same way.Both RCMs have been extensively used for simulations of past and future climates.RCA4 and HCLIM were used in the same set-up to simulate 6 ka climate (Strandberg et al. 2022).Details on model parameterisations and how the RCMs read data from EC-Earth can be found in that study.
PI and 2.5 ka are represented by 30-year periods in the RCMs.Differences in temperatures and heat fluxes between the L and R simulations are tested using a student's t-test with Bonferoni (1936) correction for multiple testing.Accounting for multiple tests is important to reduce the risk of incorrect conclusions when testing for effects across all grid cells (Wilks, 2016).The resulting procedure has a 5 % family-wise error rate (FWER), i.e. the probability of one or more false positives among all grid cells is 5 % instead of the 5 % false positive rate for each grid cell obtained when no correction is applied.An alternative would be to use a procedure that corrects for the fraction of false positives, the false discovery rate (FDR) (e.g.Benjamini & Hochberg, 1995).However, FWER tests are more conservative compared to the FDR alternatives.Only significant differences are shown in Figures 7, 8 and 10.

LPJ-GUESS
The dynamic vegetation model (DVM) LPJ-GUESS (Lund-Potsdam-Jena General Ecosystem Simulator) used in this study is an individual-based ecosystem model optimized for regional studies (Smith et al., 2001;Sitch et al., 2003;Smith et al., 2014).
LPJ-GUESS simulates the potential land cover for Europe at 2.5 ka using climate forcing (temperature, precipitation, and radiation) from EC-Earth, and the same CO 2 level as applied in EC-Earth.The land cover input to the RCM is generated by summing the plant functional types into three land-cover types: summer-green trees (ST), evergreen trees (ET) and open land (OL) (Table A1).The fraction of the non-vegetated area (or bare land, BL) was estimated by subtracting the summed land cover from one.

2.5
Pollen-based REVEALS estimates of plant cover and their spatial interpolation REVEALS (Regional Estimates of Vegetation Abundance from Large Sites) is a model that estimates regional plant cover using pollen records (Sugita, 2007).To achieve reconstructions of past land cover that can be used in climate modelling, pollen-based REVEALS estimates of plant cover are produced at a 1˚grid scale using all suitable pollen records available in each grid cell (Gaillard et al., 2010;Gaillard et al., 2018).The first REVEALS reconstruction of this kind was performed for Europe (Trondman et al., 2015) and was used to evaluate the LULCC scenarios KK10 and HYDE (Kaplan et al., 2017) and discuss the effect of LULCC on simulated regional climate in Strandberg et al. (2014).For this study and the earlier one on the climate at 6 ka (Strandberg et al., 2022), we increased the coverage of the REVEALS reconstruction southwards to the Mediterranean area and eastwards to western Russia and the Middle East (Githumbi et al., 2022).Because the REVEALS estimates of plant cover are not obtained for all grid cells, it is necessary to perform an interpolation to produce a complete pollen-based land cover at a 1 resolution across Europe.This interpolation applies a spatial statistical model described in Pirzamanbein et al. (2018).
The spatial interpolation uses large-scale features consisting of elevation and the LPJ-GUESS land cover fractions adjusted by LULCC at 2.5 ka as derived from the KK10 scenarios (Kaplan and Krumhardt, 2011).The correlation between these largescale features and the REVEALS estimates of plant cover is computed and used for interpolation.It is important to note that by estimating correlations, only the patterns in the features are used for interpolation.The effect on the interpolation of using large-scale features (i.e.LPJ-GUESS land cover) derived from different climate models was small in a sensitivity study (Pirzamanbein et al., 2020), which suggests that a single iteration between the climate model and the vegetation model is sufficient for the spatial interpolation described here.Small-scale spatial patterns in the REVEALS estimates that are not captured by the large-scale features are interpolated by a Gaussian Markov Random Field (Lindgren et al., 2011).The model uses Dirichlet distributions for compositional data (Aitchison, 1986) to ensure that the interpolated land cover fractions are between 0 and 1 and that the total land cover in each grid cell sums to 1.
The REVEALS reconstructions and the spatial interpolation provide land cover for vegetated areas, with total land cover summing to 1.To account for bare land (BL), the statistical interpolation is adjusted using the fraction of bare ground in the simulated potential land cover obtained from LPJ-GUESS.The adjustment increases the land-cover type open land (OL) (see LPJ-GUESS modelling above) in sparsely vegetated areas, such as mountainous or dry regions.
To compare the potential and reconstructed landcover used in the simulations we computed differences between each landcover class for each pixel as well as joint measures of distance between the complete landcover compositions for each pixel.Let p k and q k be the proportions of land cover in the 4 different classes, indexed by k, for the potential and reconstructed landcover in single pixels.The distances between the two compositions were computed using standard Euclidian distances as, and using the compositional distance suggested in Aitchison (1992) and Aitchison et al. (2000) as where K in the root is the total number of classes (K=4 here).
Compared to the Euclidian distance the Aitchison compositional distance (ACD) puts greater emphasis on relative differences between compositions.As an example, if the landcover fractions change from 40% evergreen, 20% summergreen, 20% open land, 20% bare land to 49%, 17%, 17%, 17% respectively the relative abundance of the last three classes stays the same (they are each 1/3 of the remaining landcover) and the ACD will correctly account for this.As a consequence of this focus on relative differences an increase from 1% to 2% will be seen as very similar (but not identical) to an increase from 10% to 20% (both represent a doubling of that composition).Thus the ACD shows larger differences for areas where one, or more, compositions are very rare.

Potential and reconstructed land cover at 2.5 ka
Below, potential land cover (L) is the climate-induced land cover at 2.5 ka simulated by LPJ-GUESS using the EC-Earth derived climate.Reconstructed land (R) is the interpolated pollen-based REVEALS-estimated plant cover, i.e. both climateand human-induced land cover (Fig. 2 and Table 1).
Potential land cover (L) is dominated by evergreen trees (ET) in northern and eastern Europe (Fig. 3).Summer green trees Both standard Euclidian distances (ED) and Aitchison compositional distances (ACD) (Fig. 3) show that the largest differences between R and L land cover are found in the westernmost parts of Europe in westernmost Norway, Ireland, Britain, western France, and northern Spain.ACDs are larger than EDs for areas where one, or several, land-cover types have small % cover (section 2.5), which is the case in Ireland and northernmost Britain.In contrast, EDs are larger than ACDs where none of the land-cover types have very low % cover (most of northern, central, and eastern Europe).Difference maps (R-L) for the three land-cover types (ET, ST, and OL), BL being the same in both R and L (Fig. A1), provide information on which land-cover type(s) most influence the differences between R and L. ET cover is smaller in R than in L land cover (-10% to -40%) in most of Boreal and northernmost Europe, in particular along the Norwegian coast (up to -40% to -80%), while it is larger in R than in L (+10% to +40%) in the southern part of temperate and Mediterranean Europe, and more so in western Spain, SW France, the French riviera, Italy, and an area in eastern Europe from the Czech Republic down to Greece and western Turkey ( + 20% to + 40%).ST cover is smaller in R than in L land cover in temperate and Mediterranean Europe, especially in Spain, France, and the British Isles (-60% to -80%), while it is larger in R than in L mainly along the coasts of Boreal and northernmost Europe (+ 10% to + 20%).OL cover is larger in R than in L in most of Europe (mainly + 10% to + 40%; ≥ 40% in W Europe) except in an area from Hungry down to northern Greece, Bulgaria, Romania, and Central Turkey where R-L is around zero or -10% to -40%.The largest R-L differences in OL are found in NW Spain, NW France, the British Islands, and the northwestern coast of Norway (+ 40% to + 60 (80)%) and in NE Greece, central Romania, southern Russia (Volgograd area), and western Kazakhstan (-20% to -40%).

Simulated climate
Figures 4 and 5 show the differences between the simulated 2.5 ka and PI temperature (at 2 m above surface) and precipitation, respectively, in winter and summer.Here we describe the results for the 2.5k-R runs because they use land-cover data inferred from empirical proxies of past land cover, i.e. pollen-based REVEALS estimates of plant cover, so far the best possible reconstruction of actual land cover in the past (e.g.Trondman et al., 2015;Kaplan et al., 2017;Roberts et al., 2018;Githumbi et al., 2022).We describe the impact of different land-cover data on climate in section 3.3.Henceforth, the difference of temperatures or precipitations between 2.5 ka (with land cover R in both RCMs) and PI are abbreviated (T 2.5k -PI) and (P 2.5k -PI), respectively.Note that the EC-Earth simulations do not use the land-cover data R for 2.5 ka, but prescribed potential land cover for 1850 CE (see Methods for more details).As expected, the general similarities between the GCM and the RCMs simulated climates are larger over the ocean, where the model resolution is less important, and smaller over land, where the land surface is treated differently in the GCM and the RCMs.
In winter, all simulations indicate that 2.5 ka was warmer than PI and show a similar pattern, with the smallest temperature differences between 2.5 ka and PI (T 2.5k -PI) on the Iberian Peninsula (0-1 C) and the largest ones in north-eastern Europe (4->6 C) (Fig. 4).The RCMs exhibit similar patterns as EC-Earth, however with smaller (T 2.5k -PI) with ca.0.5 C on the Iberian Peninsula.Moreover, HCLIM shows smaller (T 2.5k -PI) than EC-Earth and RCA4 in most of central Europe.In summer, the smallest (T 2.5k -PI) are also found on the Iberian Peninsula and in western Europe, while the largest (T 2.5k -PI) occur in the northern parts of the domain, especially in areas close to the sea-ice margin in the far north (Fig. 4).The difference in (T 2.5k -PI) between RCA4 and HCLIM is larger in summer, although patterns are similar.While (T 2.5k -PI) are 1-2 C in the HCLIM simulations, they are close to zero in central and eastern Europe in the RCA4 ones.Both RCMs simulate a (T 2.5k -PI) of 0.5 to 1.5 C in the eastern and northern parts of the domain, which is less than in EC-Earth.Overall, and in both winter and summer, the difference in precipitation between the two time periods (P 2.5k -PI) is small (within ±10 %), except for areas close to the northern rim of the domain in winter, where 2.5k-R is wetter than PI with (P 2.5k -PI) ≥50 % in all models (Fig. 5).Most regions in central Europe and the Mediterranean are characterised by low positive or 325 negative (P 2.5k -PI) (up to ±10 %).All three models show similar precipitation patterns, which suggests that precipitation is mainly governed by the large-scale circulation as simulated by the driving GCM.However, there are notable differences in (P 2.5k -PI) between the RCMs and EC-Earth, which suggests that changes in the precipitation amount results also from local processes including small-scale convection.Figure 6 shows sea level pressure (SLP) at PI and the difference between 2.5 ka and PI (SLP 2.5k -PI).In winter, the RCMs show similar patterns as EC-Earth, however with in general slightly lower pressure over the continent.Comparing the two 335 time periods, it is clear that SLP is lower at 2.5 ka over Scandinavia, the British Isles and the North Atlantic, which implies enhanced cyclone activity.Lower SLP values at 2.5 ka than PI occur farther south, and are more common over Scandinavia, in EC-Earth than in the RCM simulations.Such differences between the RCMs and the driving GCM partly explain differences in (P 2.5k -PI) between the models and between time periods.For instance, the negative SLP difference in Scandinavia and the North Atlantic may be related to that region's larger positive precipitation difference (Fig. 5).

Climate response to differences in land cover-the effect of human-induced land cover
In this section, the RCM 2.5k-R runs described in section 3.2 are used as the reference and compared with the 2.5k-L runs, i.e.
Figure 7 shows the difference in climate at 2.5 ka resulting from changing land-cover from the potential (climate-induced) land cover (L) to the pollen-based reconstructed land cover (both climate-and human-induced) (R) (Figure 2).The difference is expressed as (2.5k-R) -(2.5k-L), henceforth abbreviated R-L (Table 1).Here we present the surface temperature instead of the diagnostic 2 m-temperature, the latter being defined in different ways depending on the model (Breil et al., 2020), while surface temperature has a standard definition (Fig 7).In most of the area temperature differences are not significant except for some areas, most notably for HCLIM.In winter the R-L temperature differences are negative (-1 to -0.5 C) over the Alps in both RCMs (Fig. 7a, e).HCLIM simulates larger R-L temperature differences than RCA4 in large parts of Scandinavia and western Europe, suggesting a stronger sensitivity of HCLIM than RCA4 to the difference in land-cover data.Also in spring, differences are small.A common feature for both RCMs is that R-L temperature differences are negative (-3 to -1 C) in large areas in Scandinavia including the mountain range (Fig. 7b, f).
In summer, RCA4 responds strongly to the inferred land-cover differences.The R-L temperature difference is 0.5-2°C in southwestern Europe (Fig. 7c, g).HCLIM, on the other hand, shows very small R-L temperature differences.The only significant differences are seen over limited regions in western Europe, and over scattered points in the Iberian Peninsula and along coastal areas in the southeast.

Differences in land-cover descriptions
Simulated potential (L) and pollen-based reconstructed (R) land cover exhibit compositional differences (Figure 3, description in Results).Major differences between LPJGUESS-simulated Holocene vegetation and pollen-based REVEALS reconstructions of land cover were discussed in a study using 19 pollen records along N-S and W-E transects through Europe (Marquer et al., 2017) based on pollen-based REVEALS reconstructions of open-land cover.Other earlier studies have compared the results of DVM simulations in Europe with pollen-accumulation rates and pollen percentages (e.g., Miller et al., 2008), or pollen-based REVEALS estimates of land cover (e.g., Pirzamanbein et al., 2014).
It must be kept in mind that the prescribed simulated climate entirely determines the DVM simulated potential land cover in this study, i.e., the simulations do not account for the effects of LULCC.The reconstructed land cover is a pollen-based reconstruction of the climate-and human-induced land cover, i. e. a product of complex interactions between several natural and anthropogenic factors.Disentangling climate-from human-induced pollen-based vegetation is a challenge.Marquer et al.
(2017) used a novel approach including the use of variation partitioning to quantify the percentage of variation in vegetation composition explained by climate and land-use variables.The results suggest that climate is the major driver of vegetation when the Holocene is considered as a whole or at a sub-continental scale, but land use appears to be important regionally.Land use becomes a major control of vegetation from ca. 4.5 ka although climate is still the principal driver; however, its influence decreases gradually.Land use is the cause of the more or less continuous increase in OL and decrease in ST from ca. 7 ka (southern Europe) and 6 ka (temperate and most of boreal Europe) until 1850 CE (e.g., Marquer et al., 2014Marquer et al., , 2017;;Githumbi et al. 2022).However, changes in pollen-based OL do not necessarily represent LULCC but may be due to climate-induced land-cover changes such as spread or decrease of open wetland areas as a consequence of water-level changes (e.g.Pinke et al., 2017).Pollen from the plant families Poaceae (grasses) and Cyperaceae (sedge family) are either difficult or impossible to identify to genera or species, which hampers categorization into dry-fresh and wet meadows, i.e., all types of open land, natural or human-induced, include large amounts of Poaceae and/or Cyperaceae.Consequently, pollen-based OL cover may overestimate human-induced land-cover change.In contrast, cropland (or agricultural land, AL in e.g., Trondman et al., 2015), also included in OL, is always underestimated in pollen-based AL cover as it is represented by cereals only.This is because pollen productivity estimates are not available so far for other crops, and/or the pollen grains of many crops are not distinguishable from the pollen of wild species belonging to the same plant genera or families.Further details on the interpretation of OL cover in terms of human-induced vegetation openness can be found in Trondman et al. (2015).Despite these issues, it is assumed that most of the increase in pollen-based OL cover since the start of the Neolithic period is humaninduced, and that OL cover before the introduction of agriculture was in general < 20% (ca.10-20%; Fig. 1).
Differences between R and L land cover could, besides LULCC, be caused at least partly by differences between the "actual" climate (responsible for the climate-induced part of pollen-based vegetation) and the model-simulated climate.Causes of differences between R and L other than discrepancies in climate and LULCC include processes involved in tree migration (e.g. the Holocene migration lag in Picea and Fagus: Giesecke and Bennett, 2004;Giesecke et al., 2007;Giesecke et al. 2017) and soil development (e.g.Miller et al., 2008).Despite continuous progress in effectively modelling such processes in DVMs/ecosystem models, it is not yet implemented as a routine (Lehsten D. et al., 2014, Lehsten, V. et al. 2019).Lower ET values in R than in L land cover in most of northern Europe at 2.5 k BP are partly due to overly large abundances of Picea in the LPJ-GUESS simulation.The latter is known from earlier comparisons between LPJ-GUESS simulated vegetation and records of pollen accumulation rates (Miller et al., 2008) or pollen-based REVEALS reconstructions of plant cover (Marquer et al., 2017) over the Holocene.However, the latest comparison between LPJ-GUESS and REVEALS total tree cover in boreal Europe shows that LPJ-GUESS mean total tree cover (MTTC) is in good agreement with REVEALS MTTC in boreal Europe at 2 k BP when standard deviation is considered (Dallmeyer et al., 2023).At a subcontinental spatial scale and for PI, MTTC is larger in LPJ-GUESS simulations than in REVEALS estimates in the British Isles only, and it is lower in LPJ-GUESS than REVEALS in northernmost Europe and the alpine region.In continental and boreal Europe, the difference in MTTC is insignificant when the standard deviations are considered.The largest disagreement between LPJGUESS and REVEALS estimates in large parts of northern Europe is found instead in the relationship between evergreen and deciduous tree cover where LPJGUESS simulates larger cover of evergreen trees than their REVEALS estimated cover (Fig. A1; Marquer et al., 2017).

Response in the RCMs to land-cover changes
The temperature response to different land-cover data in winter and spring results from the snow albedo feedback and the insulating impact of the snow cover.Open land generally has a higher albedo than forests.Open land is also more readily covered by snow which makes the differences in albedo between open land and forest larger during the snow season.Thus, the temperature response is expected to be the strongest in regions of high altitude or high latitude where there is snow in winter and where the fraction of open land is relatively large.
At high latitudes the amount of incoming solar radiation is small in winter making the snow albedo response weak.The only model consistent R-L temperatures in winter are seen over the Alps.In the Scandinavian mountains, the snow season is long, stretching through spring.Here, the temperature response is stronger in spring when the amount of incoming radiation is larger, even though the difference in albedo between 2.5k-R and 2.5k-L is smaller in spring than in winter (Fig. 7a, b, e, f; see also Fig. 9).The response in temperature to land-cover differences is coherent with the differences in albedo in winter (Fig. 8a ,b).
The differences in albedo follow closely the differences in land cover, since land cover to a large extent determines albedo.
This agrees with earlier studies of idealised land-cover changes in RCMs showing that albedo explains a large part of the temperature signal in winter in Scandinavia (Strandberg & Kjellström, 2019;Davin et al., 2020).Note that differences in albedo can be small, but still significant, because the variability is small.Some of the differences in the range of ±0.1 in Figure 8a,b are < 0.05 and should be considered as not different from zero (i.e.no difference).The areas with the largest R-L temperature differences in summer (Fig. 7c, d, g, h) correspond to the areas with the largest significant differences in latent heat flux as simulated by RCA4 (Fig. 8c).Trees have deeper roots and more leaves than grass, which implies more evapotranspiration (EVT).A more open landscape, such as 2.5k-R, is expected to be warmer in summer due to decreased EVT.In contrast, HCLIM shows few significant temperature differences, as a result of very small differences in latent heat flux (Fig. 8d).Davin et al. (2020) demonstrated that RCMs could disagree in their temperature response to idealised changes in land cover in summer.There is not a direct relationship between differences in land cover and simulated changes in EVT.EVT is determined in the models by the magnitude of land-cover differences, as well as the absolute amount and type of vegetation.Furthermore, local conditions, such as available soil moisture, determine how much water is available for EVT.How these factors are translated into temperature responses depends on the individual models' physics and parameterisations.
The roughness length (Z0) is also changed by LULCC, therefore differences in Z0 could potentially explain some of the climatic responses.As Z0 differences are constant over the year, the same scale of vegetation changes would give the same Z0 difference across the domain.Consequently, roughness changes can not explain why the response to LULCC is different depending on the season or the region.We consider albedo and latent heat flux as the most likely explanations because they best correlate to temperature differences.
To better understand the relationships between surface albedo, latent heat flux and temperature, we present the annual cycles in absolute values for three different small regions (3×3 grid points) representing Scandinavia (SCA), west central Europe (WCE) and the Iberian Peninsula (IBP) (regions shown in Fig. 8).Overall, all four simulations (RCA4 and HCLIM with R or L land cover) agree in terms of shape and amplitude of the annual cycles (Fig. 9).The largest discrepancies between RCA4 and HCLIM are found in the annual cycles of EVT in WCE and IBP.InWCE, RCA4 simulates a larger amplitude and slightly higher values of EVT in the second half of the year.In IBP, it is HCLIM that simulates a larger amplitude and a longer period of high EVT.Trees have a longer period of EVT than herbs, but not necessarily a higher annual maximum of EVT (Kelliher et al., 1993;Strandberg & Kjellström, 2019).Resulting differences are seen for WCE and IBP in RCA4, and to some extent HCLIM simulations.A possible explanation for the higher latent heat fluxes in the 2.5k-R runs is that grasslands do not retain water as efficiently as woodlands, which would mean that the latent heat flux increases faster and reaches higher levels, to then drop faster.Such a pattern is seen in both RCA4 and HCLIM simulations, but it is clear that RCA4 tends to dry out more rapidly than HCLIM (Fig 9f ).

RCM simulated climate compared to proxies
To evaluate the climate-model simulations, we compare simulated temperatures to proxy data from the Temperature 12k Database (Kaufman et al., 2020).We let all proxy data represent the overall patterns of 2.5 ka temperature and compare that with the overall patterns from the models.Figure 10 shows winter and summer temperatures at 2.5 ka as given by proxy data and the climate models with R land cover in this study.
In winter, there is a general agreement among the models that temperatures range, along a northeast to southwest gradient, from -10 C and below in northeastern Europe to around 5 C in the Iberian Peninsula.There is also agreement that eastern Europe is colder than western Europe at the same latitudes, and that the Scandes mountains, the Alps and the Carpathians are colder than the surrounding areas.Compared to proxies, the models are colder on Iceland, except perhaps HCLIM, and warmer along the coasts of western Europe and parts of central and eastern Europe.In summer, there is a general agreement that temperatures range from around 10 C in the north to approximately 24 C in the south.Compared to proxies, the models are warmer in Iceland and central eastern Europe, especially in EC-Earth.
Even though there is general agreement between models and proxies there are also differences.Some proxies show considerable differences to the model data, and even if proxies and data have the same colour in Fig. 12, the difference may be up to 2 C. In addition, the proxy data is not internally consistent, which makes it difficult to evaluate their utility in datamodel comparisons.The proxy data come from lake sediment, marine sediment, peat, glacier ice, and other natural archives, made with different methods and varying uncertainties (Kaufman et al., 2020).There is also a problem with the spatial scale  Strandberg et al., 2011).For this study, it is sufficient to conclude that the models give a picture of 2.5 ka temperatures that largely agrees with that from proxy data.To decide whether the inclusion of reconstructed land cover improves the comparison with proxy data, we compare each proxy 500 data point with the value from the corresponding grid cell in all four RCM simulations.Figure 11 shows how the proxy data correlate to the simulated temperature in HCLIM in winter and summer.The general agreement between models and proxies seen in Figure 10 is seen also here as fairly strong correlations both in winter and summer.The differences between HCLIM-L and HCLIM-R are small, both the root mean square error (RMSE) and the correlation coefficient are more or less the same for both simulations (Table 3).RCA4 and EC-Earth also give similar values.It is a challenge to assess the effect of adding human-induced land cover on simulated climates in comparison to having climate-induced land-cover only, since the uncertainties in climate proxy data are large and the land-cover induced differences in climate are small compared to the absolute temperatures and temperature differences due to other forcing (e.g.Strandberg et al., 2011).Moreover, the model-data comparison (Fig. 11; Table 3) is hampered by the irregular spatial distribution of climate proxy data.They do not represent the entire European domain and are often missing in the regions where the largest differences between L and R simulations are suggested by our results.Several sites are located along the coast or in the ocean.
Even if we consider only the points in regions where we suspect large differences, such as Scandinavia in winter and western Europe in summer, the conclusion is the same.The climate proxy dataset at hand has not the necessary spatial distribution to assess whether the inclusion of human-induced land cover improves the simulated climate, and whether the RCMs improve simulated regional climate compared to the GCM EC-Earth (Table 3).

Robustness of the results
The natural variability in Europe is large (e.g.Deser et al., 2020), which influences the comparison between 2.5 ka and PI climates.Three adjacent 30-year periods in EC-Earth, considered here to represent the studied 2.5ka climate, differ by up to 1 C.A wider description of the 2.5 ka climate, including natural variability, would require ensemble simulations including more time slices.The difference in climate between 2.5 ka and PI is a result of orbital forcing and greenhouse gas concentrations (Singh et al., 2023).Since our model experiment does not include sensitivity runs changing one forcing at a time, we do not know what forcing plays the largest role.
The large-scale climate conditions as represented by the driving GCM are important for the RCM simulated regional climate (e.g.Rummukainen, 2010).The dependence of GCM-simulated changes in large-scale circulation on precipitation in the RCMs is a well-known feature seen in projections of future climate change (e.g.Kjellström et al., 2018;Christensen and Kjellström, 2020).Also, different GCMs and/or different simulations with the same GCM reflecting natural variability tend to yield large regional differences (e.g.Deser et al., 2020).For instance, natural variability alone may explain a large part of differences between models and why one model study differs from another, or from reconstructions using proxy data.
Even with identical forcing, the response to regional differences in land cover can differ significantly between model simulations.The response in several RCMs (of which RCA4 was one) to idealized land-cover changes was studied for present climate conditions (Breil et al., 2020;Davin et al., 2020).They showed large agreement across models in the temperature response to albedo in winter.In summer the analysed RCMs showed different signals in the response to heat flux and temperature.Russo et al., (2022) suggest that this could be explained by how the land-surface scheme works in the individual model and how sensitive they are to e.g.perturbations in soil moisture.The small response in latent heat flux and temperature in HCLIM in summer relative to that in RCA4 shows that sensitivity to land-cover changes differs across models, regions and seasons; however, it is not necessarily a sign of conflict or disagreement between models.Larger land-cover changes may trigger a similar response in both RCA4 and HCLIM, although testing this would require a series of simulations with incrementally larger land cover changes, and these tests have not been performed.
Consequently, using other models would have given different results.We acknowledge that we only represent a small part of the total uncertainty, and that a larger variety of climate models and/or initial conditions representing natural variability would enable a better assessment of these differences.Our way of handling model uncertainty is to use two models with different model physics.As we get a similar response as that obtained in studies of idealised climate-land-cover interactions using idealised land-cover forcing studies of under present climate conditions (Strandberg & Kjellström, 2019;Breil et al., 2020;Davin et al., 2020) and comparable studies of past climate (Russo et al., 2022), we ascribe most of the uncertainty associated with the response to differences in land-cover in our study to differences in the two models' physics.To improve the estimate of uncertainty of our results will require larger multi-model ensembles.Such ensembles are not available for 2.5 ka to date.5

Conclusions
This study describes the climate and land cover 2500 years ago (at 2.5 ka) as simulated by one GCM, two RCMs, one DVM and using one reconstruction of 2.5 ka land cover based on pollen data, statistical interpolation methods and climate model results.The results provide some insights on how simulated climate is influenced by LULCC, and how sensitive RCMs are to differences in land cover.
The results presented here suggest that LULCC at 2.5 ka impacted the climate in parts of Europe.Simulations including LULCC give up to 1 C higher seasonal mean temperature in parts of northern Europe in winter and up to 1.5 C warmer in southern Europe in summer than simulations with potential land cover.This relatively strong response suggests that anthropogenic land-cover changes may have affected European climate during the Bronze Age.This also implies that LULCC is important in future scenarios, especially in low-emissions scenarios where the greenhouse gas forcing is relatively small.
The models simulate a 2.5 ka climate that was warmer than the PI climate, as a result of different orbital forcing and feedbacks.
The largest differences are seen in northern Europe, where the 2.5 ka climate is 2-4 C warmer than PI.In summer, the difference between the simulated 2.5 ka and PI climates is smaller (0-3 C), with the smallest differences in southern Europe.
Differences in precipitation are mostly within ±10 %.In parts of northern Europe, the climate in 2.5 ka is up to 30 % wetter in winter than that of PI.In summer, there is a tendency for the 2.5 ka climate to be drier than PI in the Mediterranean region.The choice of land cover has a significant impact on the simulated temperature.Winter and spring temperatures are closely related to albedo, which is largely the same in both RCMs, and is strongly affected by land cover.In summer, the RCMs used in this study respond somewhat differently to land-cover differences, showing that the choice of land cover and climate model is important for the resulting simulated climate.Summer temperatures are strongly related to differences in heat fluxes between the atmosphere and the ground.Since the response in heat fluxes to differences in land cover depends on model physics, it is more likely that models respond differently in summer than in winter.For summer, RCA4 responds more strongly to the imposed differences in land cover than HCLIM.This explains some differences between the climate conditions simulated by RCA4 and HCLIM.It is difficult to assess which model has the most realistic response.Model performance is dependent on many other factors including, but not limited to, large-scale circulation, parameterisations and resolution.The best way to understand this model uncertainty is to use multi-model ensembles to capture the range of possible climates.
The importance of model and land-cover choice calls for caution when designing palaeoclimate experiments.Here we show that it is essential to have a good, well-motivated description of land cover to simulate the same climate with different models in a model ensemble.

Table A1
Groups of land-cover types used in this study.Ericaceae*(MTSE): the pollen productivity used for Ericaceae pollen in the REVEALS reconstruction represents the mean pollen productivity of several species of which Arbutus unedo, Erica arborea, E. cinerea and E. multiflora are dominant.The genus Calluna vulgaris (heather, LSE) also belongs to the Ericaceae family but its pollen productivity has been estimated separately (Githumbi et al., 2021).Cerealia t.: all cereals except Secale cereale (rye) that is easily that is easily separated on the basis of pollen morphology and for which pollen productivity was estimated separately.Abbreviation: t = type.of land-cover types (LCTs) and Plant Functionnal Types (PFTs) used in this study.**The most recent plant taxonomy has merged the family Chenopodiaceae into the family Amaranthaceae, i.e. "new" Amaranthaceae = "former" Amaranthaceae + Chenopodiaceae.Pollen analysts have mostly used the name Chenopodiaceae for this pollen-morphological type, but it includes all species from the two former families, therefore the name Amaranthaceae/Chenopodiaceae.

Figure 1 :
Figure 1: density plots of REVEALS-based open land (% cover) for (A) all European grid cells; (B) the boreal zone; (C) temperate Europe and (D) Mediterranean Europe.Full black lines represent the mean value of open land cover; dashed green lines represent the median value.Panel (E) is illustrative of how the density plots are derived for one time window, showing data values for all grid cells (opaque grey circles), the derived density function as a curve and expressed on the colour gradient, and the mean and median values.Key time intervals (6ka, 2.5ka and PI) are highlighted on panels A-D by the grey bars.

Figure 2
Figure 2 Description of the model chain for a) PI and b) 2.5 ka.EC-Earth is run for both periods with prescribed potential land cover for PI(i.e. 1850 CE).All RCM simulations read boundary conditions from EC-Earth.For PI, the RCMs use modern land use.For 2.5 ka, the EC-Earth climate is used in LPJ-GUESS to provide the potential land cover (L) subsequently used in the RCMs.A

Figure 3
Figure 3 Upper panels: Composite maps of (L) LPJ-GUESS simulated potential land cover using climate input from EC-Earth and (R) interpolated pollen-based REVEALS estimated land cover in Europe at 2.5 ka.Legend: ET -evergreen trees; ST -summer-295 green trees; OL -open land; BL -bare land.Lower panels: Two measures of the difference between the composite maps L and R,

Figure
Figure 7 Difference (2.5k R-L) in surface temperature (°C) in RCA4 (a -d) and HCLIM (e -h) for winter (DJF; a, e), spring (MAM; b, f), summer (JJA; c, g) and autumn (SON; d, h).Only grid cells that show a significant difference on a 0.05 level are coloured.

Figure 8
Figure 8 Difference between 2.5k-R and 2.5k-L in albedo in winter (DJF) for a) RCA4 and b) HCLIM and in latent heat flux in summer (JJA) for c) RCA4 and d) HCLIM.Black squares show the areas representing Scandinavia (SCA), West Central Europe (WCE) and the Iberian Peninsula (IBP).Only grid cells that show a significant difference (p < 0.05) are coloured.

Figure
Figure 9 Annual cycles of surface albedo (alb; a-c), latent heat flux (hfls (W m -2 ); d-f) and surface temperature (ts (K); g-i) for Scandinavia (SCA; a, d, g), West Central Europe (WCE; b, e, h) and the Iberian Peninsula (IBP; c, f, i).Data from RCA4 are 470 represented by black lines, data from HCLIM by red lines.Simulations using potential land cover from LPJ-GUESS are labelled with 'L' and shown as full lines, and those using pollen-based REVEALS reconstructed land cover are labelled with 'R' and shown as dashed lines.
of the reconstructed versus simulated climate.Proxy data may represent local climate conditions that do not necessarily match the spatial scale of the model simulated climate.Making an in-depth model-proxy comparison would require an assessment of each data point (what they represent, uncertainty ranges etc.) which is outside the scope of this paper.Future studies could improve the model-proxy comparison by e.g.compensating for differences in altitude between model and proxy data (e.g.

Figure 11
Figure 11 Simulated (y axes) and proxy (x axes) temperatures at 2.5 ka in winter (a) and summer (b).Temperatures from HCLIM-L are represented by black triangles, and from HCLIM-R by red circles.The locations of the data points are seen in Fig 12.

Fig
Fig A1 Upper panels: Cover (in fraction of grid cell) of open land (OL), evergreen trees (ET), and summer-green trees (ST) at 2.5 935 ka in Europe as, first row: simulated by LPJ-GUESS using climate input from EC-Earth (L), and second row: estimated using pollen data and the REVEALS model (R).Lower panels: Difference between L and R cover (R-L).See main text under methods and results for details.

Table 1 .
Combination of models and land cover (LC) used in each simulation.EC-Earth is used as the driving GCM in all simulations.The DVM is driven by climate conditions from EC-Earth and simulates potential LC labelled 'L'.
between the boreal and temperate biomes) starting in southern Norway and becoming larger eastwards over southern Sweden, the Baltic countries and Russia, and an area in the lowlands north of the Alps that meets the hemiboreal zone in Poland.Most mountainous areas of south-central Europe (Alps, Carpathians, Caucasus) have higher ET proportions than the surrounding plains.At high elevations and latitudes, the Scandinavian mountains are characterised by a higher proportion of open land cover (OL) and bare land (BL) than at lower elevations and latitudes, at the expense of trees.Besides those areas, OL proportions are negligible in most of northern and central Europe.OL, often with high BL proportion, has its largest proportions on the western fringe of Europe and in the south-eastern part of the continent.The reconstructed land cover (R) shows similar spatial patterns of tree cover to those in the potential land cover, with ET dominant in northern Europe and the mountainous areas of northern and central Europe, and ST largest proportions at low latitudes in central and western Europe.However, the proportions of mixed ET and ST woodland and ET/ST mixed with OL in the reconstructed land cover (R) are much larger than in the potential land cover (L).
(ST) prevail in western Europe, and mixed evergreen and summer green woodland (ET + ST, orange colours) occupy the hemiboreal zone (

Table 3 Correlation (Corr.) and root mean square error (RMSE, °C) between simulated 2.5 ka temperatures and proxy temperatures.
Two different forest types dominate the simulated potential land cover at 2.5 ka: evergreen coniferous forests in northern and eastern Europe, and deciduous broadleaved forests in western Europe.The portion of open land cover is negligible in all forested areas and dominates only in the steppe region of southeast Europe.Reconstructed land cover, however, shows mixed forests in most of Europe.Furthermore, there is a high portion of open land in forested regions, and, in addition to the steppe region, open land cover dominates in most of western and southern Europe.Altogether the reconstructed land cover is considerably more open than the simulated potential land cover in most of Europe.