The biogeophysical climatic impacts of anthropogenic land use change during the Holocene

. The ﬁrst agricultural societies were established around 10 ka BP and had spread across much of Europe and southern Asia by 5.5 ka BP with resultant anthropogenic deforestation for crop and pasture land. Various studies (e.g. Joos et al., 2004; Kaplan et al., 2011; Mitchell et al., 2013) have attempted to assess the biogeochemical implications for Holocene climate in terms of increased carbon dioxide and methane emissions. However, less work has been done to examine the biogeophysical impacts of this early land use change. In this study, global climate model simulations with Hadley Centre Coupled Model version 3 (HadCM3) were used to examine the biogeophysical effects of Holocene land cover change on climate, both globally and regionally, from the early Holocene (8 ka BP) to the early industrial era (1850 CE). Two experiments were performed with alternative descrip-tions of past vegetation: (i) one in which potential natural vegetation was simulated by Top-down Representation of Interactive Foliage and Flora Including Dynamics (TRIFFID)

The most important anthropogenic alteration of the natural environment was the clearing of forests to establish cropland and pasture and the exploitation of forests for fuel and construction materials (Darby, 1956). This long history of anthropogenic land cover change (ALCC) has implications for regional hydrology and climate and possibly for global climate. Deforestation results in both biogeochemical and biogeophysical changes. The biogeochemical changes tend to increase temperature by the emission of greenhouse gases such as CO 2 and CH 4 (CH 4 emissions are influenced not just directly by deforestation but by irrigation in rice agriculture and by emissions from livestock and humans). The impacts of biogeophysical changes are many and varied, being dependent on the local climate, soil, and the natural vegetation that is being replaced, e.g. if natural savannah or grassland is replaced by crops the impact will not be as great as if woodland is replaced.
There are several mechanisms by which biogeophysical changes due to deforestation can affect regional climate. A combination of a reduction in aerodynamic roughness, in the root extraction of moisture, and in the capture of precipitation on the canopy leads to reduced evaporation and thus decreases the fluxes of moisture and latent heat from the surface to the atmosphere. These changes work to increase the local surface temperature (Lean and Rowntree, 1993). Conversely, the increase in surface albedo due to deforestation acts to decrease surface temperature by increasing the reflection of shortwave radiation. This is particularly true at high latitudes where snow on the ground is a factor for some of the year, and the snow-covered ground is no longer masked by the canopy of the forest. Generally, in mid-to high latitudes the albedo increase is considered to be the dominant effect, leading to a net cooling of the regional surface temperature, whereas in the moist tropics, the evaporation is more important and, therefore, a localised overall warming may result (Betts et al., 2007).
During the Holocene the climate has been influenced by natural forcings. Orbital variations have caused a decline in summer solar insolation in the Northern Hemisphere over the last 6000 years. During the same period concentrations of greenhouse gases such as CO 2 and CH 4 have been increasing. On decadal to centennial timescales, fluctuations in solar and volcanic activity have also had a climatic impact (Wanner et al., 2008;Schmidt et al., 2004). The impact of ALCC is superimposed on these natural forcings. The extent and timing of these early anthropogenic land surface changes is the subject of much debate, as is their role in changing Holocene climate. Ruddiman (2003) proposed the idea that anthropogenic impacts on greenhouse gases, and consequently climate change, began thousands of years ago as a consequence of early agriculture and have been increasing in amplitude ever since, which he termed "the early anthropogenic hypothesis". The idea has been hotly debated in the literature (see, e.g., Broecker and Stocker, 2006;Joos et al., 2004;Singarayer et al., 2011;Mitchell et al., 2013;Ka-plan et al., 2011). Whilst the early anthropogenic hypothesis may likely not account entirely for the pre-industrial rises in CO 2 and CH 4 , there is no doubt that land use changes do have a climatic impact on both regional and global scales. The real debate is the scale of these effects of early agriculture.
Whilst paleoecological and archaeological evidence of anthropogenic land use changes exists, there are not enough sites to comprehensively determine continental-scale impacts of deforestation (Kaplan et al., 2009). Therefore, in order to better estimate impacts of anthropogenic land use, several databases of land use change have been developed. Examples of these include the HYDE 3.1 (History Database of the Global Environment; Goldewijk et al., 2011), Kaplan and Krumhardt 2010 (KK10;Kaplan et al., 2009Kaplan et al., , 2011, and Pongratz et al. (2008) models. Although the methodologies differ in the details, the basic premise of these models is that from an estimated database of historical population trends, anthropogenic deforestation is calculated based on population density and the suitability of land for crops or pasture.
To quantify the impact of ALCC on climate, data sets of past ALCC can be used in conjunction with climate models. Several studies have estimated the influence of pre-industrial ALCC on a global climate (He et al., 2014;Kutzbach et al., 2011;Pongratz et al., 2010). Globally, the biogeophysical effects of anthropogenic land use change have been estimated to cause a slight cooling that is offset by the biogeochemical warming, giving a net global warming (He et al., 2014;Pongratz et al., 2010). On the local to regional scale, in the most intensively altered landscapes of Europe, Asia, and North America, the biogeophysical effects can be comparable with the biogeochemical (He et al., 2014;Pongratz et al., 2010). In addition, Strandberg et al. (2014) used a regional climate model to evaluate the climatic effect of anthropogenic deforestation in Europe at 6 and 0.2 ka BP with both the HYDE 3.1 and KK10 ALCC scenarios. For the KK10 scenario at 6 ka BP, small but significant temperature differences were found in summer and, at 0.2 ka BP, changes up to ±1 • C were found over widespread areas in both summer and winter. Other authors (e.g. Oglesby et al., 2010;Cook et al., 2012) have modelled a decrease in precipitation in response to deforestation in Mesoamerica.
These existing studies are, however, limited in either temporal or spatial extent and do not address the question of when anthropogenically induced climate change first occurs. In this study global climate model simulations are used to provide a comprehensive evaluation of the influence that the biogeophysical effects of regional human-induced land cover change have had on the climate both globally and regionally throughout much of the Holocene. As described in detail in Sect. 2, the period under consideration is from 8 ka BP to the pre-industrial period (1850) HadCM3, is a coupled atmospheric, ocean, and sea ice model. The atmospheric component has a horizontal resolution of 2.5 • latitude and 3.75 • longitude with 19 unequally spaced levels in the vertical and a 30 min time step. It has an Eulerian advection scheme and includes effects of CO 2 , N 2 O, CH 4 , CFC11, and CFC12. The spatial resolution over the ocean is 1.25 • by 1.25 • with 20 unequally spaced layers extending to a depth of 5200 m. It also includes a 1.25 • by 1.25 • resolution model for the formation of sea ice with simple dynamics whereby the sea ice drifts on the ocean currents (Cattle and Crossley, 1995).
TRIFFID is coupled to the GCM (general circulation model) via MOSES every 10 days of the model run. Within TRIFFID, nine surface types are specified: five plant functional types (PFTs) and four non-vegetation types.
HadCM3 was widely used in both the third and fourth assessment reports of the Intergovernmental Panel on Climate Change (IPCC, 2001(IPCC, , 2007 and still performs well in a number of tests relative to other global GCMs (Covey et al., 2003;IPCC, 2007). For the fifth IPCC assessment, it has been superseded by Hadley Centre Global Environment Model version 2 (HadGEM2; Collins et al., 2011), but being relatively computationally efficient, HadCM3 can be the better choice for some palaeoclimate modelling applications as it allows more and/or longer runs to be conducted than would be possible with a higher-resolution model.

Project-specific model configuration
The version of HadCM3 used does not include interactive ice, carbon cycle, or methane and thus must be forced with prescribed changes in orbit, greenhouse gases, and ice-sheet evolution. Orbital parameters are taken from Berger and Loutre (1991), atmospheric concentrations of gases are determined from ice cores (CO 2 from Vostok (Petit et al., 1999;Loulergue et al., 2008) and CH 4 , and N 2 O from EPICA; Spahni et al., 2005), and the ice-sheet evolution is estimated using the ICE5G model of Peltier (2004). For further details of these natural forcings of the climate model, readers are referred to Singarayer et al. (2011).
To prescribe Holocene ALCC the KK10 data set of Kaplan et al. (2009Kaplan et al. ( , 2011 was used. The original data set has a 5' spatial resolution and has modelled crop and pasture land use for every year from 8 ka BP to present. For this study the data at 1000-year intervals were taken (8, 7 ka BP, etc.) for both crop and pasture combined and upscaled to the spatial resolution of HadCM3 to formulate a time series of cropland masks (Fig. 1). Within TRIFFID the global crop area is designated by a cropland mask, which can only be occupied by agricultural-type vegetation (i.e. C 3 and C 4 grasses) or bare soil (Betts et al., 2007). Hence, the actual cropland is equivalent to the mask area, less inland water, urban, and ice tiles and less the area covered by non-grassland vegetation or bare soil. The crop mask area is not dynamically updated by climate data. The cropland area incorporates the natural C 3 and C 4 grass fractional areas before converting tree fractions.
For each simulation, the boundary condition forcings (orbit, greenhouse gases, and ice sheets) were specified, and in all simulations the initial conditions were the same, based on a spun-up early industrial simulation. Simulations were run for 1000 years. By the final 500 years of the simulation, the climate system had adjusted to a new surface equilibrium, and thus these final 500 years were averaged to result in the mean altered climatic conditions. The relatively long averaging period increases the signal-to-noise ratio between the modified and control climates and thus distinguishes between differences that are statistically significant but which can be hidden by decadal or multidecadal variability in shorter averaging periods. This is especially important for assessing the impact of agriculture in the earlier time slices of the Holocene when the land use change is small and localised.

Local impacts of land use
In the regions where ALCC was significant, surface air temperature changes can be seen in all the time slice simulations (Figs. 2, 3, and 4), and for all time slices except 8 ka BP (not shown) the temperature anomalies in most regions are outside the normal range of variability, which is considered to be within 2 standard deviations of the mean. The anomalies are more pronounced in the JJA season ( Fig. 2) than in DJF (Fig. 3). This is due to a combination of the land imbalance between the Northern and Southern Hemisphere, the lack of land surface changes in the extratropical Southern Hemisphere, and the enhanced effect of land surface changes during the season of greatest solar insolation and plant growth (which is JJA in the Northern Hemisphere).

Initial conditions
Based on a spun-up early industrial simulation. Simulations KK10 Dynamic vegetation from TRIFFID with crop mask based on KK10 land use data set (Kaplan et al., 2009 Natural dynamic vegetation from TRIFFID Simulation dates 8, 7, 6, 5, 4, 3, 2, 1 ka BP; 1850 CE Length of simulation 1000 years Averaging period 500 years as it depends on local climate, soil, and the natural vegetation. In the extratropics, where the albedo effect is generally dominant, there is a trend towards increasing (negative) anomalies, with an increase in disturbance fraction ( Fig. 5a for PI). At 7 ka BP the range of extratropical temperature anomalies within the areas of land disturbance is +0.1 to −1.2 • C (JJA) and +0.6 to −0.5 • C (DJF), by 4 ka BP it has increased to +0.4 to −2 • C (JJA) and +0.9 to −1 • C (DJF), and by the pre-industrial period (PI; 1850) it had reached +0.7 to −4 • C (JJA) and +0.3 to −2 • C (DJF). Although there are some positive temperature anomalies, the vast majority of grid points shows a negative temperature trend. Regions with the highest ALCC intensity show the largest negative temperature anomalies, in particular Europe and E Asia and China, where agricultural land use occurs earliest and has the highest concentration of land conversion (Figs. 2, 3, 4, and 5a). The tropical response shows less of a trend because the impact of reduced evaporation is more significant, and thus there are conflicting signals between the cooling effect of increased albedo and the warming effect of reduced evaporation. In some tropical areas ALCC leads to net cooling, while in other areas net warming is simulated (Figs. 2, 3, 4, and 5a), partly dependent on the availability of moisture at the surface and partly on cloud cover changes (not shown). The main areas that show a warm anomaly response are southern Africa and India in JJA from 5 ka BP and the area bordering the Bay of Bengal in DJF where E India is warmer from 6 ka BP with warming extending to the east coast of the Bay of Bengal by 2 ka BP. The Indian JJA warming is enhanced by cloud feedbacks; a decrease in monsoon circulation leads to decreased cloudiness, thus increasing the shortwave radiation reaching the surface and warming the lower atmosphere. In contrast, tropical South America generally shows a net cooling in response to ALCC. Around the mid-to late Holocene the tropical temperature anomaly range within the areas of land disturbance is ±0.5 • C, and by the early industrial era (1850 CE) it has reached ±1 • C (Fig. 5a).
Analysis of the standard deviation of both the KK10 and control simulations indicated no significant changes in the amplitude of interannual variability of surface temperature or precipitation (using the f test statistic).

Remote impacts of land use
In addition to the local temperature changes described above, cooling can also be observed in regions remote from the areas of major ALCC, particularly in the Northern Hemisphere. The most intense cooling is always in the regions of ALCC but even as early as 7 ka BP in the JJA season our model simulations show a band of cooling that stretches across much of the extratropical Northern Hemisphere and the North Atlantic (Fig. 2). This cooling starts influencing the northern Pacific regions by 5 ka BP, and by the early industrial era the surface air temperature over most of the world's land masses and much of the ocean is cooler due to the effects of ALCC in remote areas.
In the DJF season in the Northern Hemisphere, ALCC leads to cooling both locally and regionally starting at  (Kaplan et al., 2009(Kaplan et al., , 2011 7 ka BP (Fig. 3). The model simulations also show cooling in the Arctic and warming in Siberia. Cooling remote from the areas of major ALCC becomes more extensive by 3 ka BP, and most landmasses of the Northern Hemisphere are cooler than the control simulation by 2 ka BP. The Siberian warm anomaly has ceased by 3 ka BP, but it remains less affected by the cooling than other regions. In the Southern Hemisphere, cooling remains more localised until 2 ka BP, by which time the majority of the land surface is cooler than the control.
There is an increased temperature anomaly response for the same level of disturbance fraction in the later time slices (Fig. 5b), implying that the responses to the land use changes are not just due to the local effects. Some possible mechanisms for these remote impacts are large-scale circulation changes (such as stationary waves in the upper troposphere at mid-to high latitudes and monsoonal circulation changes), near-surface advection, and the amplifying factors of snow cover. These mechanisms will be discussed in more detail in Sect. 3.2 for atmospheric dynamics and Sect. 3.3.2 for snow cover changes. Changes to the natural vegetation cover (outside the regions of land use) due to the climatic impacts of land use were also investigated as a potential mechanism of further feedbacks but were not found to be significant.

Upper tropospheric dynamics
There are several factors that can impact the dynamics of the upper troposphere in response to ALCC. Cooler surface air temperature means that the density of the air is greater and, therefore, the geopotential height in cooler regions will be lower (See Fig. 7c and d). The changes to mean sea level pressure (MSLP) described in Sect. 3.2.2 can change the areas of convergence and divergence leading to changes in regions of vertical motions. The unequal latitudinal distribution of the temperature anomalies, with the regions of greatest cooling in the midlatitudes of the Northern Hemisphere, affects the meridional temperature gradient leading to a change in baroclinicity, which has been shown to impact storm tracks (Yin, 2005).
In the JJA season from 7 ka BP there is a reduction in the 500 hPa geopotential height over the extratropical Northern Hemisphere in a pattern similar to the temperature pattern but more extensive, completely encircling the globe. From 3 ka BP onwards the height reduction expands southwards, so the geopotential height is lowered almost everywhere by pre-industrial times. The most intense reduction is always in a zonal belt across Europe and E Asia. The Southern Hemi-  sphere response from 5 ka BP appears to show a standing wave pattern affecting the subtropical highs.
In the DJF season a stationary wave in the anomaly field in both the Northern and Southern Hemisphere is apparent in all time slices but most pronounced at 4 ka BP (Fig. 7d). This is a recognised response to surface temperature anomalies as described in Hoskins and Karoly (1981) although, in this case, there are multiple thermal anomalies caused by ALCC. By 2 ka BP there is a reduction in 500 hPa geopotential height over most of the globe. Note that there are several areas that are not statistically significant in Fig. 7c and d, implying a large amount of variability. These geopoten-tial height changes contribute to the simulated remote temperature changes by altering the regions of vorticity, which in turn influence the regions of ascent and descent and thus the MSLP, and surface climatic conditions. The geopotential height anomalies can also alter the pattern of the upper level winds, thus influencing surface storm tracks.
In particular, the positioning of the geopotential height anomalies in the earlier time slices (up to 4 ka BP, Fig.7d) indicates an increased tendency towards a positive tropical-Northern Hemisphere (TNH) pattern (Barnston et al., 1991) with above average heights over the Bering Strait and Gulf of Alaska and northeastward of the Gulf of Mexico and be-  low average heights over eastern Canada. This would be expected to cause cooler temperatures over the continental United States by increasing the transport of cold polar air into the United States. In several time slices it can be seen that DJF temperature anomalies over the Bering Strait and Alaska (e.g. Fig. 3) show a warming pattern where there are positive geopotential height anomalies in Fig. 7c and d. The DJF temperature anomalies (warming) over Siberia up to 4 ka BP (Fig. 3)

Mean sea level pressure
Changes to MSLP can also have an effect on the climate system. The colder surface air temperature in the region of disturbance means reduced ascent in those regions and thus higher MSLP, and as discussed in Sect. 3.2.2 the changes in the upper troposphere can influence the surface conditions. The MSLP changes can be seen quite clearly in China from 6 ka BP and in Europe by 4 ka BP (Fig. 6  ably due to more remote influences such as the North Atlantic storm track. These MSLP changes could play a part in steering weather systems and thus influencing the climate in remote regions. For example, the JJA MSLP anomaly pattern ( Fig. 7a and b) over the North and central Atlantic is indicative of a negative phase of the North Atlantic Oscillation (NAO), with above-normal pressure over the North Atlantic and below-normal pressure over the central Atlantic. This pattern is apparent from 5 ka BP. The NAO alters the intensity and location of the North Atlantic jet stream and storm track and thus the patterns of heat and moisture transport (Hurrell, 1995). A negative NAO would contribute to a tendency to wetter summers in all but the southernmost regions of Europe (Folland et al., 2009), and this was seen in the results described in Sect. 3.3.1. The DJF NAO shows a trend towards a positive NAO which could result in drier winters over the Mediterranean region , although it is difficult to ascertain whether this is the case as the pattern is not as consistent as in JJA.

Surface advection
Some of the cooling in regions adjacent to the areas of land disturbance is due to advection by low-level winds. In regions with a prevailing wind direction an advection pattern can be clearly seen. Figure 6 shows advection of cold air from the areas of land use change in east Asia to the east across the west Pacific and from the region of land disturbance in Mexico to the west across the east Pacific; this also affects the sea surface temperatures (SSTs) in these regions of the Pacific (not shown). In other areas where the surface wind direction is more variable the effect is more difficult to detect but probably does contribute to the spread of the cold anomaly outward from the region of land disturbance.

Hydroclimate
Our model simulations show precipitation changes in response to ALCC, but it should be noted that there are larger uncertainties in climate-model-simulated precipitation and other variables related to model dynamics than in temperature, which is primarily controlled by thermodynamics (Shepherd, 2014). Different climate models show a wide range of responses in their dynamics to palaeo-and future climate change scenarios, and we acknowledge that aspects of the precipitation anomaly patterns in this study may be less robust than those of other climate variables.
The precipitation responses to ALCC (Fig. 8 for JJA and Fig. 9 for DJF) tend to be caused by a response to large-scale circulation changes rather than being directly attributable to local land use. The European precipitation response in the DJF season is not entirely consistent throughout the time slices, but the general response is a slight decrease in precipitation around the western and Mediterranean coasts, with this dryness extending further into the continent by 1 ka BP. The simulations show that, in the JJA season, Europe has an increase in precipitation from 7 ka BP onwards compared to the control; the anomaly gradually grows in magnitude and extent, possibly influenced by the increased tendency to a negative North Atlantic Oscillation (NAO). Positive anomalies begin in the warm pool of the Gulf of Mexico and extend across the North Atlantic, following the track of positive anomalies to the 850hPa wind field (Fig. 10). In addition, as the cooler temperature anomalies extend quite high in the troposphere over Europe, this increases relative humidity throughout the low to mid-troposphere (not shown) and thus the likelihood of large-scale precipitation.
In India there is a decrease in monsoon precipitation from 5 ka BP, which then gradually increases in intensity. This is partly driven by the slightly cooler Indian subcontinent temperatures (Fig. 11) in the critical months for monsoon development and also by cooler temperatures in Europe and east Asia and increased snow cover on the Tibetan plateau. This leads to decreased monsoonal circulation and decreased cloudiness. There are also changes to the east Asian monsoon, with wetter conditions to the north and south of the region and drier conditions in the centre. This pattern is seen reasonably consistently from 4 ka BP onwards.

Intertropical Convergence Zone (ITCZ)
Analysis of the precipitation fields (Figs. 9 and 10) shows an overall southward migration of the ITCZ. These changes are most obvious in the Atlantic and Pacific oceans and over the continent of Africa. In the DJF season there are changes to the ITCZ in the western Pacific, but a consistent pattern is not seen until 3 ka BP when there is southward shift in www.clim-past.net/12/923/2016/ Clim. Past, 12, 923-941, 2016 P I Figure 6. DJF surface temperature anomalies for KK10−control and KK10 surface winds for the late pre-industrial (1850 CE) and 2 ka BP. strengthening of the northward cross-equatorial energy transport (e.g. Kang et al., 2008). This shift south of the ITCZ to transport heat to the cooler Northern Hemisphere is seen in both the DJF and JJA seasons.

Snow cover
Lower surface air temperatures in the ALCC scenario relative to the control lead to an increase in winter snow accumulation (Fig. 12). This increase is seen by 5 ka BP, mostly in northern and mountainous regions. The areas affected gradually increase so that by 3 ka BP more temperate and lowerlying areas see increases in snow depth. The effects are most pronounced in North America and Europe. In regions outside the areas of permanent snow cover, increases in snow depth will delay the melting of the snow pack and thus result in a longer period of snow cover. The increased snow cover due to the cold temperature anomalies will cause additional cooling due to the increased albedo. This will be greatest in regions of deforestation where the snow-covered ground is no longer masked by the canopy of the forest. This increased snow cover would also lead to decreases in precipitation due to lower rates of moisture recycling over land. It should be noted that the modelling of snowfall is subject to the same uncertainties as described in Sect. 3.3.

Temporal evolution of Holocene climate
In the control experiment the changes in orbital configuration, greenhouse gases (GHG), and ice sheets and sea level lead to monotonically increasing global temperatures through the Holocene (Fig. 13a). Analysis of previous experiments to assess the sensitivity to different natural forcings (data from Singarayer and Valdes, 2010;Singarayer et al., 2011) suggests that while the changes to orbital configuration effect a cooling in global temperature over the Holocene, this is outweighed by increases in greenhouse gases (∼ 17 ppm CO 2 from 8 ka BP to late pre-industrial time), which result in overall warming. Cooling through the Holocene occurs in Northern Hemisphere summer, when forced with orbit and GHG variations, but is not as pronounced as when only forced by orbital variations. In winter, when HadCM3 is forced with orbit-only variation, there is little change in temperature, but when GHG increases are included, this becomes a warming over the Holocene, which then outweighs the reduced summer cooling. Whilst this contrasts with recent data compilations that suggest a general decline in global temperatures since the mid-Holocene (Marcott et al., 2013), it is within the range of other climate model responses when compared with the Paleoclimate Model Intercomparison Project 3 (PMIP3) mid-Holocene (MH) minus late pre-industrial (PI) temperature anomalies. Although palaeodata syntheses

Wi n d S t r e n g t h A n o ma l y ( ms -1 )
Wind strength anomaly ( may suggest a cooling of Northern Hemisphere temperatures, there are regional and seasonal variations in the data such as that from the Bartlein et al. (2011;Fig. 14c) and Mauri et al. (2015) data compilations. In both these compilations and the combined proxy reconstructions of Wanner et al. (2008), the cooling is most evident over the higherlatitude Northern Hemisphere. For our specific configuration of HadCM3 the inclusion of land use changes through the Holocene has a significant impact on the progression of modelled global average temperatures, such as to alter the direction of the multi-millennial trend described above, thus reducing the mismatch between the model simulations and the proxy reconstructions which show an extratropical cooling trend throughout the Holocene. The increasing magnitude and spread of ALCC through the Holocene reconstructed in KK10 counteracts the influence of increasing greenhouse gases, so that temperatures are effectively steady from 3 ka BP in HadCM3 (Fig. 14a).
These global trends are composed of considerable heterogeneity on the regional scale. When the trends are broken down into zonal regions, it can be clearly seen that the difference in trends with/without land use is greatest in the northern extratropics (Fig. 13b), where the Holocene trend is modified from increasing temperatures to decreasing temperatures in the late Holocene by the addition of ALCC. There are impacts on mean temperatures in the tropics (Fig. 13c) and southern extratropics (Fig. 13d) but not sufficient to influence the direction of the Holocene trend.
MH minus PI anomalies of annual mean surface air temperature (Fig. 14a) show a near-global distribution of cooling, except over high-latitude sea ice regions, which are particularly influenced by changes in obliquity (higher in the MH than in the PI). The PMIP3 suite of models shows a similar pattern of surface air temperature anomalies. The cooling is most dramatic over the tropics and monsoonal regions, where changes in the seasonality of insolation (due to orbital precession variation) intensify monsoon circulation in the early and mid-Holocene and the resulting addi-

T e mp e r a t u r e A n o ma l y (°C )
Temperature anomaly ( C) 0 Figure 11. Indian temperature seasonality; KK10−control surface temperature anomalies for the late pre-industrial (1850 CE) simulation; the vertical bars indicate the normal range of variability which is considered to be within 2 standard deviations of the mean. tional cloud cover reduces incoming shortwave radiation as well as increased surface water, altering the balance of sensible to latent heat fluxes. The inclusion of land use change in the KK10 experiment reduces the magnitude of MH cooling, especially in the midlatitudes. Over Europe and eastern North America the anomaly is reversed to a warming (i.e. over these regions the cooling from deforestation shown in Fig. 2 increases and outstrips the warming from greenhouse gases). These are also regions where there is the highest concentration of pollen reconstructions ( Fig. 14c; Bartlein et al., 2011;Mauri et al., 2015). The influence of land use change improves the data-model comparison with the reconstruction by Bartlein et al. (2011) over these key areas (Fig. 14a-c). Likewise, when simulated top-level ocean temperatures are compared with SST data used within the Marcott et al. (2013) compilation, the inclusion of land use improves the data-model comparison (Fig. 14d and e). However, the largest MH warming in the model is in the summer months, whereas recent seasonal temperature reconstructions (also using pollen; Mauri et al., 2015) suggest that the largest and most widespread MH warming may have occurred in winter.
In contrast, using the KK10 ALCC scenario as a boundary condition to the climate model does not improve the agreement in annual mean MH−PI precipitation anomalies when compared to the palaeoclimate reconstruction of Bartlein Figure 12. DJF snow depth anomalies (cm) for KK10 minus control for the late pre-industrial period (1850 CE). The stippling indicates grid boxes where the anomalies are significant at the 95 % level using Wilcoxon rank sum statistical analysis.
al. (2011; Fig. 15). Model and data are in reasonable agreement in most regions except for Europe and the temperate regions of Asia where the data implies a wetter MH than PI, which is not seen in the model runs. The difference over Europe is exacerbated by the inclusion of land use, which results in a drier MH.

Discussion
Anthropogenic land cover change leads to climate change well beyond the core regions of land use early in the Holocene. These results suggest that regional ALCC has an effect on the atmospheric circulation, e.g. the ITCZ shift is a remote response on global scale. The implications of this finding are that regional models or atmospheric-only models would not simulate these atmospheric circulation changes as well as a global coupled model. In this study we observed multiple thermal anomalies (from intense regions of cooling directly over anthropogenic land use change), but the standing wave response of the geopotential height field would likely also be seen even for a single thermal source from just one region (Hoskins and Karoly, 1981). Regional models have the advantages of higher spatial resolution and more detailed orography, but they may not include these potential remote atmospheric changes, possibly resulting in different impacts from the same land cover forcing in regional and global model simulations. The positioning of the major temperature anomalies in the midlatitudes and at similar latitudes may be particularly significant in producing the stationary wave pattern. Whilst these are the results from only one model, there are many similarities in the distribution of the temperature anomalies with those found by He et al. (2014) for 1850 CE and Pongratz et al. (2010) for the 20th century although the temperature changes found in this study were greater, e.g., a pre-industrial global annual mean temperature anomaly of −0.23 • C as opposed to the −0.17 • C estimated by He et al. (2014). Running similar simulations with a greater number of models would improve the robustness of the results particularly with respect to hydroclimate due to the high uncertainties involved. The variability in the results from different models can be greater than the variability of the property that is being assessed (Pitman et al., 2009;de Noblet-Ducoudre et al., 2012;Brovkin et al., 2013). These inconsistencies have been attributed to disagreements in how land use change is implemented, the parameterisation of albedo, the representation of crop phenology and evapotranspiration, and the partitioning of available energy between latent and sensible heat fluxes (Pitman et al., 2009;de Noblet-Ducoudre et al., 2012;Boisier et al., 2012). The albedo and turbulent heat fluxes from our model simulations for the North America-Eurasia region (Fig. 16) are within the range of other climate model responses when compared with those from the Land-Use and Climate, IDentification of robust impacts (LUCID) set of experiments . The negative turbulent and latent heat fluxes would offset some of the cooling due to the increased albedo. Although the largest albedo changes are in the DJF season the impact of this will be lessened due to the lower levels of incoming solar radiation in this season. Results could be further improved by the running of transient simulations that could capture events such as the Maunder minimum. A transient simulation response could result in different biogeophysical impacts to the ones achieved when using a long equilibrium simulation where the ocean-land-atmosphere system can reach more of a steady state and the climate sensitivity is different.  From late preindustrial era simulations, one using observed atmospheric greenhouse gas concentrations and the other using greenhouse gas concentrations in a world with no anthropogenic emissions (based on linear projection from earlier Holocene trends from Kutzbach et al., 2011), He et al. (2014 estimated a net global warming of 0.9 • C due to the biogeochemical effects of ALCC, with between 0.5 and 1.5 • C warming in the areas of most intense land use changes. If we incorporate this degree of warming into our early industrial era (1850 CE) simulations, there would still be a net cooling in Europe, E Asia, and NE America with, e.g., a net cooling of up to 2 • C in parts of Europe. To put this in perspective, the IPCC (IPCC, 2014a) consider a temperature rise of more than 2 • C to be undesirable and that changes of 1 • C could have an impact on vulnerable ecosystems. However, the temperature changes in this study took place over a much longer period than the time frame considered in the IPCC and ecosystems and human societies would have had more time to adapt. The consequences of these changes for agricultural societies would vary depending on the pre-existing conditions. For example, in drier regions, where crops are more likely to be water limited, cooler, wetter summer conditions may have been beneficial to the agricultural output although the risk of erosion would have been increased. The generally lower temperatures might also make societies more vulnerable to further transient cooling effects such as volcanic activity.
There are discrepancies between simulations of the mid-Holocene climate and the independent data-based reconstructions. Both the simulations from this study and virtually all the PMIP3 (Palaeoclimate Model Intercomparison 3; https://pmip3.lsce.ipsl.fr) models show a temperature increase from the mid-Holocene to the PI, whereas the Marcott et al. (2013) and Mann et al. (2008) (on a shorter timescale); reconstructions show a decrease. It is interesting to note that ALCC reduces this mismatch for HadCM3, especially in key areas such as Europe. Other factors that could lead to this discrepancy are uncertainties in the proxy reconstructions and deficiencies in climate models (e.g. Lohmann et al., 2013). These climate model deficiencies include low resolution and sensitivity and, importantly, their dependence on soil moisture, whereby energy is utilised for evaporation rather than for temperature increase.
This study shows a significant increase in precipitation over Europe with increasing land use, which means that the PI becomes wetter than the mid-Holocene, which leads to increases in soil moisture and changes in the sensible to latent heat flux balance, and, in combination with increased albedo caused by deforestation, this results in cooler temperatures for PI than MH. If the land-atmosphere coupling strength was different and soil moisture was strongly reduced with deforestation, it is likely that the cooling effect would be smaller (cf. Strandberg, 2014).
Further uncertainties arise from the robustness of the land use reconstructions, which is difficult to evaluate due to the lack of global-scale evidence for human impact on the Earth's land surface. Much of the uncertainty comes from the lack of knowledge about the magnitude and distribution of the global human population and the rate of technological evolution and intensification through time. As part of our initial investigations simulations were also run using an alternative land use scenario (the HYDE 3.1 data set; Goldewijk et al., 2011). The HYDE 3.1 reconstruction has substantially lower levels of land use early in the Holocene (as compared with KK10), which resulted in a later development of consistent temperature anomalies at 4 ka BP (not shown) in comparison with the KK10 land use scenario. The decision was taken to proceed with the KK10 data due to its assumptions of a larger per capita land use earlier in the Holocene when agricultural methods were less efficient. Several ongoing international initiatives that aim to synthesise palaeoecological and archaeological data promise to lead to more robust reconstructions of Holocene ALCC in the future (e.g. the Past Global Changes (PAGES) LandCover6k project; http: //www.pages-igbp.org/workinggroups/landcover6k/intro).
By the early industrial period, simulated biogeophysical temperature changes in the regions of land disturbance are of the same order of magnitude (e.g. 0.83 • C annual anomaly in the main agricultural areas of Europe) as the changes seen due to CO 2 increases during the industrial period (0.85 • C; IPCC, 2014b). Part of Ruddiman's original hypothesis (Ruddiman, 2003) is that pre-industrial global warming caused by anthropogenic CO 2 or CH 4 emissions should have been ∼2 • C at higher latitudes, but there was no evidence for this warming. Ruddiman (2003) attributed this to a natural cooling trend caused by decreasing summer insolation. This study suggests that biogeophysical effects of the land use changes may also have played a part in counteracting the warming due to anthropogenic greenhouse gas emissions, as acknowledged in Ruddiman (2013). The precipitation changes might also have an impact on the availability of water for rice irrigation and on natural wetlands, thus affecting the production of methane.

Conclusions
In our global model simulations that use a Holocene ALCC scenario as a boundary condition, a surface temperature response to the biogeophysical effects of ALCC is seen in regions of early land use such as Europe and SE Asia as early as 7 ka BP in the JJA season and throughout the entire annual cycle by 2-3 ka BP. Areas outside the major regions of ALCC are also affected, with virtually the whole globe experiencing significant temperature changes with a net global cooling of 0.22 • C by the pre-industrial period. Although the temperature changes are predominantly cooling, some regions such as India, southern Africa, and Siberia show warming as a response to ALCC. The greatest changes are generally seen in the JJA season, with a mean regional cooling of 1.4 • C experienced in Europe and 1 • C in E Asia in the early industrial period (1850 CE). Much of the precipitation response to the land use tends to be due to large-scale circulation changes such as a decrease in the intensity of the Indian monsoon, the southward movement of the ITCZ, and changes to the North Clim. Past, 12, 923-941, 2016 www.clim-past.net/12/923/2016/ Atlantic storm track. In Europe there is a slight decrease in precipitation in the DJF season and a more substantial increase in the JJA season. Some causal factors for the teleconnections are advection by surface winds, MSLP anomalies, and tropospheric stationary wave train disturbances in the mid-to high latitudes. The potential for an early global impact of ALCC on climate strongly implied by this study suggests that due consideration of this should be taken in simulations covering the Holocene. The inclusion of ALCC in the model improves the model comparison for surface air temperature with the data-driven palaeoclimate reconstructions, especially in key areas such as Europe. The remote teleconnections seen in this study have implications for the regional modelling of land use change due to circulation changes that occur outside the domain of the regional model. Overall, our model simulations indicate an increase in global surface air temperatures through the Holocene. Globally, the inclusion of ALCC data reduces the magnitude of this warming especially in the late Holocene when the temperatures remain relatively constant. Regionally, in the northern extratropics, this warming is reversed in the late Holocene. It should be noted that in this study it is not possible to distinguish the anthropogenic component of the biogeochemical changes as the same atmospheric CO 2 and CH 4 concentrations (from ice core measurements) are prescribed for both the KK10 and control simulations. However, the level of early industrial warming due to the biogeochemical impacts of ALCC predicted by He et al. (2014) would negate the early industrial biogeophysical cooling seen in this study in all regions except for the most intensively altered landscapes of Europe, E Asia, and NE America.
Other caveats are the large uncertainties in the land use data and, therefore, in our understanding of the Holocene evolution of land surface-climate interactions as well as our ability to evaluate climate models. To reduce these uncertainties there is an urgent need to extend land cover reconstructions and the prehistory of land use globally (cf. PAGES LandCover6k initiative). It must also be stressed that these are the results from only one model and the use of equilibrium simulations may produce different results from a transient simulation. However, these uncertainties are more likely to impact the timing, intensity, and regional details of the changes rather than the overall trends and patterns.

Data availability
Data are available from the Bristol Research Initiative (BRIDGE, 2016) for the Dynamic Global Environment website: http://www.bridge.bris.ac.uk/resources/simulations. The relevant experimental runs are tdna (control) and tdpo (KK10).