Relative importance of forcings and feedbacks in the Holocene temperature

The pre-industrial Holocene provides the backdrop for the emergence of civilisations and the starting point of anthropogenic climate change. Several reconstructions show an early Holocene warming that was followed by cooling for several thousand years before Industrialisation. In contrast climate simulations show warming throughout the Holocene. Whilst reconstructed trends over ocean can be reconciled with warming through either seasonality or uncertainties, a consistent explanation for cooling trends over some land areas is missing. We present a suite of transient Holocene climate model simulations with a coupled general circulation model and show that a widespread mid-to late-Holocene cooling emerges over some regions of the northern hemisphere with the inclusion of anthropogenic land-use. This is mostly because in regions of prescribed late Holocene deforestation, the simulated early to mid-Holocene is characterised by a lower albedo than the late Holocene. Whilst this cooling through time can quantitatively explain some regional aspects of the reconstructions, the model-data agreement remains imperfect, and diﬀerences between reconstructions also hinders the evaluation. Moreover, model-dependency in the response of several feedbacks, particularly sea-ice, but potentially also clouds, means that it is diﬃcult to uniquely attribute Holocene temperature evolution to speciﬁc factors. Future work should aim to derive a consensus-signal including uncertainties from the available proxy data which could be used to ﬁngerprint simulations covering a range of plausible feedback strengths.


Introduction
The Holocene epoch encompasses 11,700 years of relatively stable climate since the end of the last glacial period.This is the environment in which agriculture and the first civilisations evolved (Roberts, 2014;Scott, 2017).It is the natural background state from which the climate system is shifting away under human influence (Schurer et al., 2017).The Holocene climate evolved under mostly gradual changes in climatic forcings (Roberts, 2014).The resultant changes in surface temperatures can be used to characterise the pre-industrial climate state and hence provide the context for future change (Schurer et al., 2017).The variations are also relevant to aspects of human history (Fagan, 2005), yet the feedbacks and responses are not fully understood.In particular, it remains to be shown whether the temperature peaked early  Renoult et al. (2023) and Hopcroft and Valdes (2015) for HadCM3.
levels (Bereiter et al., 2015).Transient climate model simulations of the Holocene (Liu et al., 2014) largely reflect this in showing a 1 • C warming during the Holocene.This model-data disagreement has been termed the 'Holocene temperature conundrum' (Liu et al., 2014).
A recent compilation of pollen reconstructions from Europe and Eastern North America shows that this early to mid Holocene warmth was confined to the ocean proxy records (Marsicek et al., 2018).However, a synthesis of land and marine records that integrates some of these pollen records shows a much stronger early Holocene thermal maximum (Kaufman et al., 2020).Subsequent analyses of the marine component which dominates the reconstructed warming in the Temperature 12k database (temp12k), suggest that either seasonality (Bova et al., 2021b) or uncertainties and a lack of spatial coverage (Osman et al., 2021) are responsible for differences with model simulations.Elsewhere much of Asia and most of Western North America has warmed throughout most if not all the Holocene (Herzschuh et al., 2023a,b;Zhang et al., 2022).Importantly, both the magnitude and spatial distribution of Holocene temperature anomalies varies between different reconstructions (e.g.Bartlein et al., 2011;Marsicek et al., 2018;Kaufman et al., 2020;Herzschuh et al., 2023b;Zhang et al., 2022), meaning there is not yet a consensus.
A few climate modelling studies have reproduced some aspects of early Holocene warmth.We discuss these examples in the following and provide a summary of key examples in Table 1.
One proposed explanation for the existing model-data discrepancies is that reduced dust from the 'green' Sahara during the early-to mid-Holocene caused a significant warming of the North Atlantic (Liu et al., 2018).However, the model used by Liu et al. (2018) did not include longwave radiative effects of dust (LWRE).More comprehensive models that include LWRE show a range of responses to an idealised greening (Pausata et al., 2016;Hopcroft and Valdes, 2019;Thompson et al., 2019).This range is most likely due to the fact that dust-radiation interactions are sensitive to the dust scattering and absorption coefficients used in models and to the size-distribution of dust particles, since smaller particles mostly cause scattering (hence cooling) and larger particles tend to absorb at longer wavelengths leading to warming (e.g.Miller et al., 2014;Albani and Mahowald, 2019).When prescribed in one Earth System model, a Holocene greening of both the Sahara and Arabia with up to date dust particle properties causes a relatively small net-radiative effect across the North Atlantic (see figure 4 by Hopcroft and Valdes, 2019) which means that the impact on sea-surface temperatures may be small.Accordingly, in Holocene simulations with idealised high and low dust, Thompson et al. (2022) found that the effect of dust on global climate was much smaller than that due to changes in vegeta-tion.However, dust-cloud interactions are expected to have a significant impact in palaeoclimates (Sagoo and Storelvmo, 2017;Thompson et al., 2019), and their role during the Holocene remains to be investigated in detail.
The altered vegetation study by Thompson et al. (2022) shows that a global early to mid-Holocene warm phase could have been caused by changes in the distribution of vegetation as they induce a global signal that is consistent with temp12k by replacing large regions of grasslands and open-land globally with forest.This highlights the potential for vegetation feedbacks on climate at the global scale (Claussen, 2009).However, uncertainties remain because of the idealized nature of the vegetation changes that have been imposed by Thompson et al. (2022).
In transient simulations Bader et al. (2020) showed that Arctic seaice anomalies could enhance early to mid-Holocene warming caused by stronger summer insolation at this time.Summer reductions in sea-ice extent persist into winter in their model, leading to an annual-mean warming signal.This process was also identified in some other models by Park et al. (2019).As the summer insolation subsequently declined, Arctic sea-ice would have covered a larger area, cooling the wider region and providing an explanation for some of the trends in mid-latitude surface temperature reconstructions (Bader et al., 2020).The strength of this Arctic sea-ice feedback is model-dependent (Park et al., 2019;Bader et al., 2020) and the spatial footprint does not correlate particularly well with Northern extra-tropical climate reconstructions (Park et al., 2019), a point we return to below.
An important additional forcing of Holocene climate could have arisen through early land-use (He et al., 2014;Smith et al., 2016).Archaeological evidence shows that agriculture commenced in the early Holocene in many regions of the globe (ArchaeoGLOBE Project Members, 2019;Ellis, 2021), leading by 4 kyr BP, to a globally unprecedented acceleration in the rate of change in natural vegetation (Trondman et al., 2015;Mottl et al., 2021) and in lake sediment accumulation rates (which suggest a widespread deforestation) (Jenny et al., 2019).The spatial intensity of early human land-use is contested (Klein Goldewijk et al., 2011;Kaplan et al., 2011;Klein Goldewijk et al., 2017).Thus, a widespread alteration of the land-surface by humans in recent millennia is plausible (Ellis, 2021) but remains controversial (Klein Goldewijk et al., 2011;Pongratz et al., 2008;Harrison et al., 2020).
The emerging, more complex view of Holocene temperature evolution (Herzschuh et al., 2023a,b;Cartapanis et al., 2022;Zhang et al., 2022) together with additional forcing by volcanoes and land-use and feedbacks from natural vegetation and sea-ice among other processes have yet to be fully evaluated in the context of Holocene temperature conundrum.Here we describe a suite of transient coupled climate model simulations and evaluate the relative importance of forcings and feedbacks over this time interval.

Methods
We performed a suite of transient coupled general circulation model simulations covering the last 10,000 years.We use the coupled general circulation model HadCM3BB-M2.1aD(Valdes et al., 2017) which has been widely used to study past (e.g.Singarayer and Valdes, 2010;Hunter et al., 2019;Roberts and Hopcroft, 2020), present (Stott et al., 2000) and future (Murphy et al., 2004) climate and is among the best models of its generation (Reichler and Kim, 2008).HadCM3's atmospheric general circulation model has a horizontal resolution of 3.75×2.5• (longitude-latitude) with 19 unequally spaced vertical levels and the ocean has a horizontal resolution of 1.25×1.25 • with 20 vertical levels.The land surface is a tiled patchwork of land-cover types including fully dynamic vegetation (Cox, 2001).HadCM3 (Valdes et al., 2017) is used here in a new configuration (-BB: with updated model parameters) that interactively reproduces the greening of the Sahara during the Holocene through changes to the convection scheme (Hopcroft et al., 2021) and the parameterisation of dynamic vegetation (Hopcroft and Valdes, 2021).Although the model has been superseded by more recent Hadley Centre configurations, these models are orders of magnitude slower than HadCM3.HadCM3's computational efficiency and relatively skillful climatology therefore make it a suitable choice for multi-millennial simulation lengths needed to study the timedependence of Holocene climate.
All HadCM3 simulations here have been driven with calculated changes in Earth's orbit (Berger, 1978), variations in the major greenhouse gases reconstructed from ice-core records (Bereiter et al., 2015;Loulergue et al., 2008;Schilt et al., 2010;Köhler et al., 2017) and a reconstruction of ice-sheets and sea-level (Argus et al., 2014;Peltier et al., 2015).Ocean salinity follows the 'melt-uniform' scenario of Ivanovic et al. (2016).At each timestep a global freshwater flux is calculated such that any difference between the simulated globally-averaged salinity and that implied by the ice-sheet reconstruction is redistributed over the volume of the ocean model.Initially the solar constant is fixed at 1361 Wm −2 and volcanic eruptions are not considered.Land-use is set to zero globally and variations in tropospheric aerosols are not considered.This simulation is initialised from the 10 kyr BP state of a transient deglaciation simulation with HadCM3B-M2.1d (Valdes et al., 2017) which was run from 23 kyr BP to 10 kyrBP following the PMIP deglaciation protocol described by Ivanovic et al. (2016).A 10,000 year simulation from 10 kyr BP with orbital, greenhouse gas and icesheet forcings is labelled ORB+GHG+ICE (Hopcroft and Valdes, 2021).From this simulation four further simulations were branched, incrementally adding forcings due to variations in solar irradiance (Vieira et al., 2012) (+SOLAR), volcanic eruptions (Sigl et al., 2022) (+SO-LAR+VOLC) and anthropogenic land-use (ALU).We evaluated a low (+SOLAR+VOLC+ALU(HYDE) Klein Goldewijk et al., 2017) and a highintensity land-use reconstruction (KK10) (+SOLAR+VOLC+ALU(KK10): Kaplan et al., 2011).These simulations are analysed with respect to the hydroclimate of northern Africa by Hopcroft and Valdes (2022).
More details about the volcanic aerosols are provided in the Supporting Information.Areas of anthropogenic crop and pasture are imposed in the model as a fractional area of disturbed land in which only grasses are allowed to develop.This disturbance mask is updated every 500 model years up until 500 yr BP after which it was updated at 200 and 100 years BP.The KK10 land-use is left constant at CE 1850 values from CE 1850-1950 as this reconstruction does not cover the later period.The different forcings are compared in Figs. 1 and S1.
Currently it remains an open question on the true intensity of palaeoland-use (Harrison et al., 2020).This is because the conversion of population densities to land-use area through time are highly uncertain (Klein Goldewijk et al., 2011, 2017), among other uncertainties.KK10 predicts a more gradual increase in land-use over time, so that the spa-  (Houghton et al., 1983;Ramankutty and Foley, 1999;Pongratz et al., 2008;Klein Goldewijk et al., 2017).
Pollen-inferred land-cover (Trondman et al., 2015) is one of the few independent ways to validate these estimates.Over Europe this comparison shows that the upper-bound KK10 reconstructed from populationbased method (Kaplan et al., 2017) underestimates open-land compared with the pollen-based estimate.Consistent with this, a more qualitative comparison between a synthesis of archaeological data and the lower-bound HYDE3.2 reconstructed from population-based method (Klein Goldewijk et al., 2017) suggests that HYDE3.2 may underestimate land-use globally in the millennia preceding Industrialisation (ArchaeoGLOBE Project Members, 2019).
The potential human-induced fluxes of CO 2 , CH 4 and N 2 O are implicitly included in the ice-core greenhouse gas records used in all simulations.Therefore, only the biophysical land-use effects are discussed in the following.The five simulations are compared in Table 2 and further details are provided in Supporting Information Text S1.We compare our modelling approach with previous studies in Table 1.

Calculation of radiative forcing and feedbacks
We calculate the radiative forcing over the Holocene based on the changes in Earth's orbit, mixing ratios for CO 2 , CH 4 and N 2 O, icesheet extent and orography, reconstructed variations in solar irradiance (Vieira et al., 2012), stratospheric aerosol optical depth due to volcanic eruptions (Sigl et al., 2022) and anthropogenic land-use.The radiative forcing due to orbital variations (f  ) is given by: where   is the planetary albedo for the present day (t 0 ) and S↓   (t) is the incoming solar radiation at the top of the atmosphere at time t before present and Δ implies the change from the present-day values.
We use simulated monthly-mean S↓   and   .The radiative forcing due to solar irradiance variations (f  ) is calculated in a similar manner as: where S ⊙ (t) is the solar constant at time t before present and Δ implies the change from the present-day value.
Greenhouse gas radiative forcing is calculated using validated polynomial approximations of the results from line-by-line radiative transfer calculations (Etminan et al., 2017).The radiative forcing due to icesheets follows an approach used before (Hopcroft and Valdes, 2015).The short-wave (SW) component is estimated using the Approximate Partial Radiative Perturbation Method (APRP) (Taylor et al., 2007) applied to monthly model outputs.The APRP is also used to calculate radiative effects arising from land-use.APRP constructs a simplified representation of the land-surface and atmospheric column in each gridbox based on monthly-mean cloud cover and radiation diagnostics for the surface and top of the atmosphere for both all-sky and clear-sky conditions.This is used to derive the radiative effects due to changes in cloud-cover, surface properties and non-cloud atmospheric constituents and is shown to reproduce more detailed radiative kernel methods to within around 15% (Taylor et al., 2007).The long-wave emission due to orographic change over ice-sheets is calculated from monthly model outputs as: where  is emissivity and is assigned a value of 0.85 for ice,  is the Stefan-Boltzmann constant and T are the monthly-mean simulated surface air temperatures in the pre-industrial simulation.T  is the temperature field corrected for elevation changes (Δℎ) due to palaeo ice-sheet elevation change.T  is calculated using a relatively conservative lapse rate =6.5 Kkm −1 .The radiative forcing due to land-use induced changes in surface albedo is also estimated using monthly model outputs and the APRP method and is based on the anomalies between simulations with and without prescribed land-use.This was calculated at 100-year intervals using 50-year mean climatologies.The radiative forcing due to volcanic eruptions is calculated by multiplying the global-mean stratospheric aerosol optical depth by -19 Wm −2 (Gregory et al., 2016).The shortwave feedback effects for the land-surface, clouds and clear-sky are also calculated using the APRP method described above.

Existing palaeo-temperature reconstructions
The marine compilation by Marcott et al. (2013) (M13) shows peak temperatures occurring at the beginning of the Holocene.This signal is mostly driven by records from the North Atlantic (Marsicek et al., 2018).Over land, Marsicek et al. (2018) (M18) reconstructed a warming up until around 4000 years before present (yr BP) followed by around 0.4 K of cooling based on pollen records.A global synthesis (temp12k) incorporates some of this pollen database with other terrestrial, ice-core and the M13 records and also shows a two-phase evolution (Kaufman et al., 2020), but with a much stronger warming early in the Holocene referred to as the Holocene thermal maximum (HTM), see Fig. 2.An expanded pollen dataset (Herzschuh et al., 2023a,b, H22) confirms the results from M18 over Europe.However, the temperature trends are opposite in Western North America and in many regions of Asia where there is essentially no early Holocene warmth.This implies that the conundrum is not global.
As temp12k, M18 and H22 share much of the same source data over the extra-tropical continents it is important to understand the reasons for these differences.In Fig. 2 we show that temp12k is very similar to a blend of the M13 (ocean-dominated) and M18 (land-only) reconstructions.From the marine records, M13 imparts a strong cooling signal which is only partially offset by the warming over land in M18.A simple average of the two reconstructions is very nearly equal to temp12k for the northern hemisphere mid-latitudes.Thus, despite the large amount of overlap in data between M18 and temp12k, they differ substantially.H22 show that blending the Legacy1.0 (from Herzschuh et al., 2023a) with the ocean component of temp12k results in a much weaker mid-Holocene warming than in temp12k, but this is also masking sometimes opposing temperature trends in different continental regions.
The potential for a seasonal bias in Holocene temperature reconstructions has long been recognised (e.g.Liu et al., 2014, and references therein).A potential seasonal bias in the marine component of M13 and temp12k has been modelled empirically by Bova et al. (2021b) (B21) for tropical sites although their approach has generated much debate (Zhang and Chen, 2021;Bova et al., 2021a;Laepple et al., 2022;Bova et al., 2022).B21's correction transforms the 0.5 K cooling signal in temp12k to a continuous warming trend that is in agreement with climate simulations with CCSM3 (Bova et al., 2021b).A statistically-based palaeoclimate assimilation of a global marine dataset (using forward models rather than calibrated palaeo-temperatures) has also shown a warming throughout the Holocene (Osman et al., 2021).The main reason identified are uncertainties in key records that show the strongest early Holocene warmth and spatial gaps in the ocean reconstructions.Another reanalysis using calibrated temperatures shows the opposite but with only a modest early Holocene peak (Erb et al., 2022), potentially highlighting the importance of the model-based prior information in the data-assimilation methodology (e.g.Annan et al., 2022).
If marine reconstructions are biased by either seasonality (B21) or by uncertainties (Osman et al., 2021), then it is very likely that much of the Earth System has been gradually warming over the last 10,000 years, in agreement with climate model simulations.The contrasting trends over different land areas potentially point to regional forcing or circulation trends, which are unlikely to be explained by seasonal biases in reconstructions because unlike marine and aquatic records that mostly track peak seasonal warmth (e.g.Lorenz et al., 2006;Samartin et al., 2017), terrestrial records are thought to register integrated growing season conditions (Marsicek et al., 2018;Herzschuh et al., 2023a,b).

Transient climate model simulated temperature changes
In Fig. 3 the simulated surface air temperatures in the five simulations are compared with reconstructions of Holocene temperatures.The results are plotted as deviations from the pre-industrial.The modelled trends are sampled from the gridcells that correspond with the site locations in each of the reconstructions, as shown in lower panels of Fig. 3.
All three simulations without land-use show a long-term warming by approximately 1.2 • C which is slightly larger in amplitude than with previous simulations with other models (Smith and Gregory, 2012;Liu et al., 2014;Bader et al., 2020) or data-model assimilation (Osman et al., 2021).HadCM3 has a higher climate sensitivity (around 3 K) than many models used before in transient Holocene simulations (e.g.CCSM3 is around 2.5 K) and it includes dynamic vegetation which can act to amplify warming over land (see below).The simulated temperatures stabilise in the last two millennia as polar regions begin to cool under reducing annual mean (and Boreal summer) insolation.The inclusion of variations in solar irradiance and volcanic eruptions do not appreciably alter the long-term temperature trend.
The forcing by anthropogenic land-use mitigates warming in the final millennia of the simulations, especially over Europe to some extent in Eastern North America and Asia.With the more intense KK10 land-use the simulated temperatures stabilise at around 5000 years BP over northern mid-latitudes (Fig. 3a) and then declines from around 3000 years BP.The HYDE3.2 land-use simulation begins to cool very slightly after around 3000 years BP over NH mid-latitudes but with a much smaller magnitude than indicated by the M18 pollen-based estimate (Marsicek et al., 2018).The no-land-use and HYDE3.2 simulations begin to exceed the uncertainty range of the M18 reconstruction  2021b) and d)-g) regional reconstructions by Herzschuh et al. (2023bHerzschuh et al. ( , 2022)).The model outputs have been averaged only at the gridcells where there are sites from each reconstruction.The spatial masks applied to the simulations are shown below each timeseries.
(Fig. 1a) by around 4000 years BP and show too much warming in all regions except western North America up to the present day.The HYDE3.2 land-use slightly improves the comparison, but the warming is still too great in comparison with M18 (Fig. 3a)).With the more intense land-use reconstructed in KK10 the simulation remains mostly within the reconstruction 95% limits for almost the entire time period in Fig. 1a, and shows a cooling over the last three thousand years in agreement with this reconstruction.
The temp12k reconstruction (Kaufman et al., 2020) shows a gradual Holocene cooling from around 6.5 kyr BP which is not matched by any of the simulations (Fig. 3b).The model-data agreement is very slightly improved with the inclusion of the more intense KK10 land-use, but only after 3 kyr BP after which the rate of cooling in the reconstruction and the model (KK10) are approximately equal.The simulations without land-use or with the HYDE3.2land-use show a long-term warming during the early-and mid-Holocene which is not in agreement with the reconstruction.
The regional components of H22 are compared in the lower panels of Fig. 3.Over Europe as described above, the KK10 simulation produces a long-term cooling towards the present day after a period of relatively stable temperatures that is in agreement with the record, although the early Holocene warmth in the reconstruction is not fully captured by the model.Over western North America the modelled and reconstructed trends are reasonably similar with a modest warming of 1 C mostly in the early Holocene when ice-sheets were retreating in east of the continent.Over eastern North America the more extreme KK10 landuse brings the simulation closer to the H22 reconstruction but, as with Asia, the simulations warm too much during the early to mid-Holocene.
The simulated sea-surface temperature (SST) trends show a warming by around 0.7 K for the sites used by B21.The B21 'corrected' annual mean reconstruction shows around half this magnitude of warming (Fig. 3c).The original temp12k reconstruction at these sites shows substantial cooling by around 0.5 K, which is not seen in any of the present model simulations or earlier studies with a different model (Liu et al., 2014;Bova et al., 2021b).It is only in the last 1000 years that the simulations without land-use or with the HYDE3.2land-use begin to exceed the 95% uncertainty limit of the 'corrected' reconstruction.Although as described above KK10 cools over mid-latitudes, here at the tropical marine locations, this simulation continues to warm gradually.However, this warming is muted relative to the other simulations.In this way the KK10 temperatures remain within the reconstructed uncertainty limits (Fig. 3c) throughout the reconstruction time interval.
Recently, Zhang et al. (2022) showed that northern hemisphere pollen records indicate a stable seasonality throughout the Holocene.In Fig. 4 we show the change in seasonality of the simulated temperature.In contrast to the pollen-inferred stable seasonality, all of the simulations analysed here show a significant reduction in temperature seasonality in the northern hemisphere extra-tropics which follows the trend in insolation seasonality.This is consistent with the results from the TRACE-21k simulation analysed by Zhang et al. (2022).This modeldata discrepancy highlights the differential trends in seasons which may be at the root of the Holocene conundrum, since the model simulations resemble evolving insolation, whereas the reconstructions indicate a more complex response of the system.This warrants further investigation.
One additional feature that is evident in the comparisons in Figs. 3  and 4 is that no single model simulation is able to reconcile all (or even most) of the different reconstructions.This is also seen in a comparison of simulated precipitation and reconstructions.The modelled annualmean precipitation and precipitation minus evaporation are higher during the mid-Holocene than the pre-industrial over northern midlatitudes.This is the opposite of a recent compilation by Routson et al. (2019) who argued for dryer mid-latitudes in the early to mid-Holocene.In contrast, pollen-based reconstructions show wetter conditions across much of Eurasia but the opposite over North America (Bartlein et al., 2011(Bartlein et al., , 2017;;Herzschuh et al., 2023b).Thus, just as for temperature, disagreements between reconstruction and between some of these and simulations are evident.Further evaluation of the hydrological cycle is beyond the scope of this work, but would very likely benefit from the inclusion of water isotopes in the model which is the subject of ongoing development.Taken together these mismatches show that either the model is missing some critical combination of forcings and feedbacks, or that the reconstructions are not consistent with each other.If the latter is true then inter-reconstruction divergences need to be better understood before the Holocene temperature conundrum can be resolved.

Model-data comparison for the mid-Holocene
Spatial patterns of temperature change for the mid-Holocene (6 kyr BP) are shown in Fig. 5 and Fig. S2.This comparison identifies significant differences between the three reconstructions and the datamodel assimilation of last glacial to Holocene climate (LGMR: Osman et al., 2021) where the latter shows opposite anomalies over Eurasia and North America.The only feature common to all three reconstructions is a warming centred over central and north eastern Europe and eastern North America.Elsewhere Bartlein et al. (2011) and H22 both show some evidence for cooler conditions in much of Asia and western North America.The H22 reconstruction is more spatially coherent than the compilation of Bartlein et al. (2011).In many locations H22 and temp12k show differing sign of change.Potential reasons for these differences are discussed by H22.Over land, the LGMR data assimilation shows a cooling over North America and warming over all of Eurasia throughout the Holocene (Fig. S2).This pattern differs substantially from all of the other reconstructions and model simulations.This warmcool North America-Eurasia pattern affects all Holocene time-slices in LGMR but the cause of this signal is unclear.
Without land-use none of the simulations show any warming over Europe or eastern North America (Fig. 5d).In the HYDE3.2simulation (Fig. 5e) there is a weak signal that is only around half the magnitude of the reconstructed anomaly.The KK10 simulation is able to reproduce the magnitude of warmer conditions over Europe and thereby reconcile the differential spatial anomalies evidence in the reconstructions.The absence of the land-use effect (Smith et al., 2016) is probably why existing model simulations show a poor correlation with reconstructions regardless of whether they simulated warmer or cooler Arctic during the mid-Holocene (Park et al., 2019).Elsewhere, the simulations are generally overly cool although the disagreement among the records makes firm conclusions difficult.The inclusion of more intense land-use mitigates this cooling over much of the Northern Hemisphere (Smith et al., 2016) but does not fully reconcile the model with the data.
The simulated warming over Europe in the HYDE3.2 and KK10 simulations is predominantly a response to deforestation (shown in Figs.S1  and S3) which causes the surface albedo to increase by around 0.05 as more open landscapes are generally higher in albedo than forested regions (Fig. 6).In the HYDE3.2simulation this albedo increase and the associated cooling signal are much weaker.

Holocene radiative forcing and response
To better understand the simulated responses, the radiative forcing was calculated (see Methods).The individual timeseries are compared with the simulated global mean temperature trends in Fig. 7 for the period from 10 kyr BP to present.The gradual change in orbit results in a global-mean increase in radiative forcing of around 0.1 Wm −2 over the past 10 kyr (Fig. 7a).Regionally there is a redistribution of insolation from the high latitudes to the tropics (Park et al., 2019;Bader et al., 2020).Over the poles annual-mean insolation reduced by 3.5 Wm −2 between 10 kyr BP and the present day.Over the mid-latitudes (30-60 • N) insolation remained approximately constant and it increased by around 1.0 Wm −2 over tropics.
Between 10-6 kyr BP greenhouse gas radiative forcing does not vary much, but rising sea-level and waning ice-sheets induce a radiative forcing increase of approximately 0.8 Wm −2 (Fig. 7b).After this, CO 2 and CH 4 begin to rise and dominate the global mean radiative forcing, while the contribution from N 2 O remains small.The forcing due to reconstructed variations in solar irradiance are modest and are mostly smaller than ±0.2 Wm −2 (Fig. 7c).Volcanic eruptions show an early Holocene Fig. 5. Reconstructed and simulated annual-mean temperature anomalies during the mid-Holocene (6±0.5 kyr BP with respect to the pre-industrial (0-1000 yr BP).a) Bartlein et al. (2011) reconstruction, b) temp12k (Kaufman et al., 2020) reconstruction, c) Legacy1.0 (Herzschuh et al., 2023a,b) reconstruction, d) the +SOL+VOLC simulation, e) the +SOL+VOLC+ALU(HYDE) simulation and f) the +SOL+VOLC+ALU(KK10) simulation.For the simulated anomalies, the stippling indicates where the difference between the two time periods is significant at the 95% level according to a two-sided Student's t test.period of enhanced activity (Fig. 7d) but otherwise show a highly variable forcing through time.
The radiative forcing due to anthropogenic land-use is the only forcing other than insolation over the Arctic with a negative trend over the Holocene (Fig. 7e).Averaged globally, the KK10 reconstruction causes a 0.5 Wm −2 reduction in radiative forcing which mostly occurs during the last 3000 years (Fig. 7d).The HYDE3.2 reconstruction shows a much smaller impact (Figs.7e).The figure shows differences relative to the late Holocene pre-industrial and so differences between the two scenarios cause the land-use forcings to diverge at earlier times.In all simulations except KK10 and HYDE3.2, the total radiative forcing (Fig. 7f) increases until the end of the simulation.In the HYDE3.2simulation the negative land-use forcing causes a departure from the simulations without land-use from around 1500 years BP.From around 3 kyr BP in the KK10 simulation there is a gradual forcing reduction that continues until the onset of Industrialisation.
The simulated global mean temperature response (Fig. 7g) shows a close correspondence with the global RF and implies an Earth System sensitivity (ESS) of around 4 K for a doubling of CO 2 but it should be remembered that the land use forcing is very regionally variable.This ESS value is caused by feedbacks due to clouds, other atmospheric constituents such as water vapour, the cryosphere and vegetation.In all of the simulations except the KK10, land-surface changes due to dynamic vegetation and changes in snow-cover act to amplify simulated warming in the extra-tropics (Fig. 8).The simulations show a more open vegetation distribution in the early-to mid-Holocene in mid-latitudes (40-60 • N) of central Asia and North America (Fig. S3), but little northwards expansion of the tree-line, in agreement with palaeo-botanical data from North America (Williams, 2002;Williams et al., 2009) and the Arctic (Bigelow et al., 2003) but not elsewhere.There is also no significant contraction of the grasslands in mid-latitudes over Eurasia (Binney et al., 2017).The more open canopy in mid-latitudes allows snow to produce a more comprehensive coverage of the canopy during the winter months.The resultant cooling is supported by the B11 reconstruction shown in Fig. 5 (Bartlein et al., 2011) but agreement with H22 or temp12k is less clear.
Over the tropics, dynamic vegetation induces the opposite feedback mostly in response to the simulated desertification of the western Sahara at around 6 kyr BP (Hopcroft andValdes, 2021, 2022).The simulated early to mid-Holocene green Sahara causes a regional cool-   (Berger, 1978), b) greenhouse gas mixing ratios (black) and ice-sheets (blue), c) solar irradiance (Vieira et al., 2012), d) volcanic eruptions (Sigl et al., 2022) smoothed with a 10-year running-average, e) anthropogenic land-use (Kaplan et al., 2011;Klein Goldewijk et al., 2017), and f) the sum of contributions due to orbit, GHGs and ice (blue) and incrementally adding solar irradiance (red), volcanic eruptions (cyan) and land-use relative (green or brown), all relative to the mean over the last 500 years.g) Simulated time-series of global mean surface air temperature ( • C) for the same combinations of forcings as shown in f), both smoothed with a 10-year running-mean.c) Global-mean short-wave (SW) radiative effects in HadCM3 calculated with the Approximate Partial Radiative perturbation method (APRP).APRP was applied every 100 years using repeated 50-year climatologies.The clear-sky timeseries for the three simulations that include volcanic eruptions show substantial anomalies due to variations in the 50-year mean values of the volcanic stratospheric aerosol optical depth.d)-f) anomalies of the three SW components in ORB+GHG+ICE at 6 kyr BP relative to 0 kyr.ing in HadCM3 (Fig. 5).This is because, despite the darker land surface arising from the expansion of vegetation cover, the expansion of the rainfall systems into this region produces more cloud cover and increases latent heat fluxes which both act to cool the surface (as shown in Fig. 8d) and e)).This monsoon-cooling relationship is a robust result across GCM simulations of the Holocene (Braconnot et al., 2007;Brierley et al., 2020).The end of the green Sahara could therefore cause a modest warming signal that is possibly abrupt (de Menocal et al., 2000).A small step-increase in temperature does in fact occur around 5.5 kyr BP in the M18 reconstruction, but this does not contribute to a reversal of the overall warming trend in that reconstruction.Such a step-increase is not evident in the globally averaged simulation results.HadCM3 shows the opposite impact of Sahara greening to the simu-lations described Thompson et al. (2022) in which a prescribed 100% vegetated Sahara and Arabia induces a global warming that halves the magnitude of the discrepancy with the temp12k reconstruction as compared with the Holocene simulations employing present-day vegetation.
Averaged globally (Fig. 8), the strongest feedback in HadCM3 is due to short-wave effects from clouds.The cloud feedback acts to amplify warming with a strength of 0.7 Wm −2 over the course of the Holocene.The land-surface and clear-sky feedbacks are approximately half this magnitude and both also act to amplify warming.In the simulation with KK10 land-use the natural vegetation feedback is more than compensated for by the land-use forcing which leads to a net negative contribution to the radiation budget (Fig. 8a).Changes over the oceans are much smaller in magnitude, mostly because HadCM3 does not pro- Coloured shading shows 6kyr minus pre-industrial cover, the black contours indicate 6kyr BP 10-50% ice margins, and markers show the locations of marine cores showing perennial (black triangles) and seasonally-free (red triangles) early to mid-Holocene conditions (from de Vernal et al., 2021).duce a significant change in sea-ice (Fig. S4).Recent evidence suggests that the central Siberian sector of the Arctic was seasonally ice-free during the early to mid-Holocene (de Vernal et al., 2021).HadCM3 simulates very little change in sea-ice cover here and so most probably underestimates the strength of sea-ice feedbacks (Wadhams, 2016).
In Fig. 9 we compare the JJA mid-Holocene Arctic sea-ice cover anomalies in three of the HadCM3 simulations with existing model simulations from TRACE21 Liu et al. (2014) and MPI-ESM1.2Dallmeyer et al. (2021).JJA-averages are chosen to enable the identification of potential seasonally ice-free conditions.All models show a decline in JJA sea-ice cover which is consistent with northern hemisphere summer orbital forcing.Only MPI-ESM shows a mid-Holocene sea-ice cover that is consistent with the study of de Vernal et al. (2021) who found that sites in the Eastern Arctic were seasonally-ice free but sites closer to Greenland remained perennially ice-covered.Their site locations are shown by black versus red markers in Fig. 9.The contours lines of 6kyr BP sea-ice coverage (shown in black) show that only MPI-ESM correctly captures the transition between the perennial and seasonal sea-ice cover zones.HadCM3 also simulates a thinning of the sea-ice, particularly in the KK10 simulation and this is in the correct location in the eastern sector of the Arctic, but the sea-ice is still too thick with perennial cover that reaches very close to the Eurasian coastline much like in the TRACE simulation.These two models do not therefore agree well with the records of de Vernal et al. (2021).It is unclear whether this disagreement between models is due to differences in the sea-ice schemes, for example because in HadCM3 meltponds are parameterised rather than directly represented or because of other differences in the overall strength of seasonal polar amplification.In future work, simulations forced with different distributions of sea-ice could help to clarify the potential importance of sea-ice in Holocene climate.This is the subject of ongoing work.

Discussion
Reconstructions of Holocene climate show either a limited warming or a cooling until Industrialisation.Since much of the simulated warming is driven by the increase in atmospheric carbon dioxide, this raises the question of whether the reconstructed temperature change could be used to constrain or learn about the climate sensitivity of the Earth System.A very high climate sensitivity (ECS) value e.g.> 5 K has been found in some recent Earth System Models Meehl et al. (2020).We can speculate about how such a model would perform were it integrated for the Holocene as in this study.For an ECS of around 5 K (i.e. S  2 = 5/3.7=1.35 K/Wm −2 ), the peak warming in the KK10 scenario (at 4 kyr BP) could be as much as 1.3 K instead of around 0.9 K in HadCM3.Regional variations would also be amplified and it is very likely that model-data comparisons shown in Fig. 3 would be significantly degraded with a higher ECS.Conversely, a lower ECS value may improve the model-data agreement especially considering the lessintense land-use scenario.The ECS value in HadCM3 is around 3 K which is lower than the ESS (Earth System sensitivity) value of 4 K calculated herein.The HadCM3 ECS value is comparable to the central estimates from the uncertainty distribution on ECS derived from multiple lines of evidence (Sherwood et al., 2020).
The Holocene Arctic cooling mode identified in MPI-ESM by Bader et al. ( 2020) is likely to be important for understanding Holocene temperature changes.We did not find a strong signature of this process in HadCM3.This could mean that the process itself is model dependent, possibly because of varying complexity in the schemes used to represent melt-ponds or the accumulation of snow on sea-ice.The spatial footprint of the Arctic warming mode is also an important constraint, but it is not clear that models that simulate a warmer Arctic during the mid-Holocene are in better agreement with spatially distributed temperature reconstructions (Park et al., 2019) and it may be that land-use forcing discussed above is also required.More recent Hadley Centre models (HadGEM2-ES and HadGEM3-GC3.1)both show warmer conditions in the Arctic in response to increased summer insolation (Williams et al., 2020).The CMIP5 model HadGEM2-ES also simulates warmer conditions over Europe and eastern North America in good agreement with reconstructions (e.g.Bartlein et al., 2011;Herzschuh et al., 2023b).However, this response is not driven by sea-ice but arises because of a simulated expansion of forest cover (see Fig. S5).
The above analysis of HadGEM2 lends some additional support to the role of land-use forcing in the spatial climate anomalies during this time period.However, in both of the land-use scenarios imposed, the HadCM3 simulations presented here still have discrepancies with the available reconstructions.This calls for further improvements in the modelling of Holocene climate.An obvious candidate is the cryosphere where further work is required to understand the sea-ice feedback.Emerging reconstructions (e.g. de Vernal et al., 2021;Axford et al., 2021;Detlef et al., 2023) show that sea-ice in the Arctic in particular was strongly impacted by climatic conditions in the early to mid-Holocene and HadCM3 fails to replicate this.Other factors to be considered include the wider sensitivity of models to orbital forcing (beyond the sea-ice response), and the initial conditions particularly in the ocean which have not been evaluated in detail.
Although not directly considered here, another potential consequence of anthropogenic land-use is the impact on atmospheric greenhouse gas levels, particularly carbon dioxide (Ruddiman et al., 2016).The early anthropocene hypothesis suggests that more intense early land-use globally altered the natural course of atmospheric constituents and thereby prolonged the current interglacial (Ruddiman and Thomson, 2001;Ruddiman, 2003;Kaplan et al., 2011;Ruddiman, 2013;Ruddiman et al., 2020).This has proven controversial and has not been supported by modelling of the methane (Singarayer et al., 2011) or by carbon isotope modelling (e.g.Elsig et al., 2009;Stocker et al., 2017) and is hindered by uncertainties around the magnitude of anthropogenic, peatland and oceanic carbon fluxes (e.g.Brovkin et al., 2019).Further refining our understanding of Holocene climate requires a better understanding of the drivers of atmospheric CO 2 and CH 4 but this beyond the scope of the present study.
If there were a robust target for Holocene temperature evolution, such differences between models could be evaluated so that competing representations of the climate system are benchmarked against longterm climate trajectory.This would allow the strength of key feedbacks (i.e. the Earth System Sensitivity, ESS) to be effectively inferred from palaeoclimate information.Whether the sum of these effects, i.e. the Earth System sensitivity (ESS) could be meaningfully constrained by the Holocene temperature record will depend on the uncertainties and the spatial coverage of the records.For example, sea-ice feedbacks would likely require a well-resolved record of high-latitude climate through time, whereas uncertainties in cloud feedbacks are probably more important over low-latitude ocean regions (Sherwood et al., 2020).
However, while some of the forcings are relatively well known, such as arising from Earth's orbital variations or greenhouse gases, others are subject to much wider uncertainties.Bayesian statistical inference provides a natural way to integrate the uncertainties in both the forcings and the palaeoclimate reconstructions (Rougier, 2007).In Bayesian methods each is represented as time-dependent probability distributions that account for the varying knowledge through time and for the relative uncertainty levels.The result is a conditioned probability distribution on the feedback strengths (or Earth System Sensitivity).Incorporating water isotopes and other physical variables like precipitation or even direct simulation of proxies (Dee et al., 2015;Osman et al., 2021) could provide valuable additional constraints in such a model.

Conclusions
The temperature history of the Holocene is determined by a range of forcings and feedbacks.Many of these are relatively small in amplitude.Although the previously fundamental disagreement between model simulations and reconstructions of Holocene temperature has been reduced (Marsicek et al., 2018;Bova et al., 2021b;Osman et al., 2021), there is no widely accepted explanation for the divergent trends between land and ocean and different reconstructions disagree on the regional temperature amplitudes.Our simulations show that anthropogenic land-use could provide an important additional regional forcing with a negative trend on millennial timescales.The regional expressions could help to explain divergent temperature trends in some parts of the globe whilst allowing for a warming over in other locations (Bova et al., 2021b;Osman et al., 2021).However, further work is required to understand model-data discrepancies in some regions and to quantify whether differences between models could be benchmarked against reconstructions.

CRediT authorship contribution statement
POH designed the study, performed simulations and analysed results; PJV set up the transient simulation boundary conditions; MT & MS produced the volcanic aerosol optical depth reconstruction; BS analysed proxy reconstructions and produced Fig. 2. All authors contributed to the manuscript and commented on the text.

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.
Model output generated in this study is available for further analysis from http://www .bridge.bristol.ac.uk/resources /simulations.

Fig. 1 .
Fig. 1.Timeseries of relevant variables used in forcing the Holocene climate simulations.a) Annual mean insolation over the northern tropics, mixing ratios of b) CO 2 , c) CH 4 and d) N 2 O, e) sea-level contributions from ice-sheets and f) area of anthropogenic land-use.

Fig. 4 .
Fig. 4. Simulated change in the seasonality of surface air temperature (JJA minus DJF) averaged over mid-latitude land of the northern hemisphere (a), southern hemisphere (b), and ocean mid-latitudes in the northern hemisphere (c) and southern hemisphere (d).The gridcells averaged are shown by the masks under each panel.

Fig. 7 .
Fig.7.Estimated difference in radiative forcings over the Holocene with respect to the late Holocene pre-industrial.These are arising from variations in a) Earth's orbit around the Sun(Berger, 1978), b) greenhouse gas mixing ratios (black) and ice-sheets (blue), c) solar irradiance(Vieira et al., 2012), d) volcanic eruptions(Sigl et al., 2022) smoothed with a 10-year running-average, e) anthropogenic land-use(Kaplan et al., 2011;Klein Goldewijk et al., 2017), and f) the sum of contributions due to orbit, GHGs and ice (blue) and incrementally adding solar irradiance (red), volcanic eruptions (cyan) and land-use relative (green or brown), all relative to the mean over the last 500 years.g) Simulated time-series of global mean surface air temperature ( • C) for the same combinations of forcings as shown in f), both smoothed with a 10-year running-mean.

Fig. 8 .
Fig.8.a)-c) Global-mean short-wave (SW) radiative effects in HadCM3 calculated with the Approximate Partial Radiative perturbation method (APRP).APRP was applied every 100 years using repeated 50-year climatologies.The clear-sky timeseries for the three simulations that include volcanic eruptions show substantial anomalies due to variations in the 50-year mean values of the volcanic stratospheric aerosol optical depth.d)-f) anomalies of the three SW components in ORB+GHG+ICE at 6 kyr BP relative to 0 kyr.

Table 1
Selected data-model comparisons of Holocene temperature change.
* Equilibrium climate sensitivity values (ECS) from