Interactive comment on “The impact of Sahara desertification on Arctic cooling during the Holocene” by F. J. Davies et al

Abstract. Since the start of the Holocene, temperatures in the Arctic have steadily declined. This has been accredited to the orbitally forced decrease in summer insolation reconstructed over the same period. However, here we present climate modelling results from an Earth model of intermediate complexity (EMIC) that indicate that 17–40% of the cooling in the Arctic, over the period 9–0 ka, was a direct result of the desertification that occurred in the Sahara after the termination of the African Humid Period. We have performed a suite of sensitivity experiments to analyse the impact of different combinations of forcings, including various vegetation covers in the Sahara. Our simulations suggest that over the course of the Holocene, a strong increase in surface albedo in the Sahara as a result of desertification led to a regional increase in surface pressure, a weakening of the trade winds, the westerlies and the polar easterlies, which in turn reduced the meridional heat transported by the atmosphere to the Arctic. We conclude that during interglacials, the climate of the Northern Hemisphere is sensitive to changes in Sahara vegetation type.


Introduction
The Holocene is characterized by an early thermal maximum ( ∼ 11-6ka) in the Northern Hemisphere, followed by gradually declining global temperatures that persisted up until the recent period of anthropogenically induced warming.This cooling was most prevalent in the high northern latitudes with July temperatures, north of 60 • N, decreasing by 3-4 • C from 9 to 0 ka (Renssen et al., 2005), which is attributable to the orbitally forced reduction in summer insolation (Renssen et al., 2005(Renssen et al., , 2009) which, at 65 • N, decreased by 42 Wm −2 over the same period (Berger, 1978).The early Holocene positive summer insolation anomaly (declining by 36 Wm −2 at 30 • N from 9 to 0 ka) also had a strong impact on the vegetation in northern Africa through a strengthening of the summer monsoons, leading to enhanced precipitation and grassland vegetation in the Sahara region (Kutzbach and Street-Perrott, 1985).This phase is often referred to as the African Humid Period (AHP), which lasted until the mid-Holocene, although the exact timing of its demise and the subsequent rate at which it occurred is still a contentious issue (deMenocal et al., 2000;Kröpelin et al., 2008).Following the AHP, the Sahara, under the influence of the long-term decline in summer insolation, evolved into a desert environment.
The termination of the AHP and the associated vegetationclimate interactions have been studied in numerous climate model studies (Claussen et al., 1999;Renssen et al., 2003;Liu et al., 2006;Notaro et al., 2008;Timm et al., 2010).However, these simulations were mostly focused on the regional climate effects in northern Africa, implying that the impact of the Saharan desertification on the global climate has not yet been addressed in model studies.It is important to realise that the Sahara covers a large area (approximately 9.4×10 6 km 2 ), making it plausible that changes in its surface albedo, caused by the large-scale shift in vegetation during the Holocene, could have profound effects on global climate.For instance, previous studies have shown that drastic vegetation changes Published by Copernicus Publications on behalf of the European Geosciences Union.
in the Amazon can influence winter precipitation in the North Atlantic and Europe, which are brought about by large-scale circulation changes in the mid-and high latitudes (Gedney and Valdes, 2000).Accordingly, in this study, we have designed a series of experiments to assess the first-order impact of long-term radiative forcing of the Sahara during the Holocene.This involved performing a large suite of equilibrium simulations using LOVECLIM, an Earth model of intermediate complexity (EMIC).Given the large number of simulations performed (29) and the long time integration involved in performing such an array of simulations (equivalent to 83 000 model years of simulation), we feel that this type of model is ideally suited to address such an issue.Our discussion focuses on how Sahara vegetation changes during the Holocene have contributed to cooling in the Arctic (defined herein as north of 66.5 • N).

Model
We applied LOVECLIM, an Earth model of intermediate complexity (Goosse et al., 2010), which is comprised of a coupled atmosphere, ocean, sea-ice and vegetation model.The atmospheric component (ECBilt) is a quasi-geostrophic model of the atmosphere with a horizontal T21 truncation and three vertical layers, 800, 500 and 200 hPa (Opsteegh et al., 1998) and includes a full hydrological cycle.Clouds are prescribed based on the ISCCP D2 data set (1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995) from Rossow (1996), with the total upward and downward radiative fluxes computed as a function of this data set (Schaeffer et al., 1998;Goosse et al., 2010).Clouds and precipitation are decoupled from one another within the model, with the precipitation dependent on the total precipitable water content between the surface and 500 hPa, which is a prognostic variable within the model (Goosse et al., 2010).This variable is transported horizontally using a fraction (60 %) of the sum of the geostrophic and ageostrophic winds at 800 hPa to account for the fact that humidity is generally higher closer to the surface, where wind speeds are lower.Above 500 hPa, it is assumed that the atmosphere is completely dry, and thus, all water above this level is converted to precipitable water.Precipitation also occurs below 500 hPa, if the total precipitable water passes a threshold.This threshold is defined by 0.83 times the vertically integrated saturated specific humidity, assuming constant relative humidity within the layer (Goosse et al., 2010).
The oceanic component (CLIO) is a primitive-equation, free-surface general ocean circulation model (Deleersnijder and Campin, 1995;Deleersnijder et al., 1997) coupled to a thermodynamic-dynamic sea-ice model (Fichefet andMorales Maqueda, 1997, 1999), and has a resolution of 3 • × 3 • latitude-longitude and realistic bathymetry.The vegetation component (VECODE) is a reduced-form dynamic global vegetation model, and is capable of simulating the dynamics of two plant functional types, trees and grasses, as well as desert as a dummy type (Brovkin et al., 2002).These vegetation types have an effect on the surface albedo, soil moisture content and net precipitation.
LOVECLIM simulates a modern climate that is in reasonable agreement with observations, including large-scale distributions of temperature and precipitation, and sea-ice cover in both hemispheres (Goosse et al., 2010).It has a climate sensitivity to a doubling of the atmospheric CO 2 level of 1.9 • C, which is in the lower end of the range reported for general circulation models (GCMs) (Flato et al., 2013).It has successfully been applied to simulate various past climates, including the last millennium (Goosse et al., 2005), the last glacial maximum (LGM) (Roche et al., 2007), the Holocene (Renssen et al., 2009), the 8.2 ka event (Wiersma and Renssen, 2006), and the last interglacial (Bakker et al., 2013).The results of these simulations are consistent with those of comprehensive GCMs (Goosse et al., 2010;Bakker et al., 2013;Nikolova et al., 2013).

Experimental design
In order to calculate the contribution of Sahara desertification to cooling in the Arctic during the Holocene, a series of experiments were designed.The experiments can be broken into two categories.Firstly, we performed a series of transient Holocene simulations (9-0 ka) that allowed us to capture the natural evolution of the Sahara within our model.Secondly, we applied the vegetation fractions of the Sahara at 9, 6 and 0 ka from the transient simulations in a series of equilibrium experiments.Further details of each suite of experiments are given below.

Defining the Sahara
An initial transient simulation from 9 to 0 ka was performed, forced with appropriate orbital parameter settings (Berger, 1978) and greenhouse gas concentrations (Loulergue et al., 2008;Schilt et al., 2010), whilst the solar and volcanic forcings were fixed at pre-industrial levels.This simulation is referred to as OG.In addition, another simulation, which included the same forcings as OG, plus additional Laurentide (LIS) and Greenland (GIS) Ice Sheet meltwater fluxes and topography changes, was performed.This simulation is referred to as OGGIS.The LIS meltwater fluxes were based on the reconstructions of Licciardi et al. (1999), and those for the GIS on Blaschek and Renssen (2013).The associated topographic and surface albedo changes of the LIS were based on reconstructions by Peltier (2004) and applied at 50-year time steps.However, GIS topographic changes were not accounted for because the changes are only minor at the spatial resolution of our model.For a more detailed description of the experimental setup of OGGIS, the reader is referred to Blaschek and Renssen (2013).In both OG and 573 OGGIS, global vegetation changes were calculated interactively using VECODE.

Equilibrium simulations
The transient experiments OG and OGGIS provided simulated estimates of the evolution of the Sahara over the Holocene.Using these transient simulations, we could thus define the vegetation composition of the Sahara at 9, 6 and 0 ka under both OG and OGGIS forcings.Following this, a series of 18 equilibrium experiments were designed (see Appendix A) to assess the impacts of the different forcings on the temperature change from 9 to 6, 9 to 0 and 6 to 0 ka (see Appendix B).In these simulations, the relative percentages of Saharan vegetation cover at 9, 6 and 0 ka (see Appendix C for OG vegetation fractions), taken from both the OG and OGGIS simulations, were combined with orbital and trace gas levels for 9, 6 and 0 ka, which followed the guidelines of PMIP3 (http://pmip3.lsce.ipsl.fr/).The name of each experiment reflects the forcings that were used.For example, 0k6kEQ_ OG refers to the orbital and greenhouse gas forcings being set to 0 ka levels, vegetation was set to a 6 ka state, EQ indicates that this was an equilibrium simulation, and OG means that vegetation state was derived from the OG simulation.This nomenclature is used throughout the manuscript.
All equilibrium simulations were run for 2500 years, allowing the model, in particular the deep oceans, to reach a quasi-equilibrium state (Renssen et al., 2006).The initial conditions for these simulations were derived from a control experiment with pre-industrial forcings.For the OGGIS equilibrium simulations, LIS and GIS topography were kept constant.However, LIS and GIS melt fluxes were not included, as this would have resulted in continuous freshening of the oceans during the 2500 years of model integration, implying that the simulated climate would not reach a quasi-equilibrium, in particular in the deep oceans.Since our analysis is focused on a comparison of the equilibrium responses to forcings (see Appendix B), such analysis of LIS and GIS melt fluxes would not be particularly meaningful.Although we recognise this as a potential weakness of our 9 ka equilibrium simulations, it must be noted that with an earlier version of LOVECLIM, experiments performed including LIS melt fluxes showed that the impact of these on the Nordic seas was a modest 0.5 • C cooling between 9 and 8 ka (Blaschek and Renssen, 2013).Therefore, we feel that omitting this potentially important source of freshwater forcing to the high northern latitudes will not adversely affect our results and we are confident that our series of sensitivity experiments are suitable in assessing the first-order impact of radiative forcings on the Holocene climate.
From the OG and OGGIS equilibrium experiments we were able to deduce the contribution of Sahara desertification to Arctic cooling during the Holocene (see Appendix B).Data presented are averages over the last 500 years of each simulation.To investigate an "extreme" example of Sahara desertification between the mid-and late Holocene, we also performed six equilibrium simulations that had 100 % grass and desert in the Sahara at 9, 6 and 0 ka (see Appendix D).In LOVECLIM, the albedo of grass and desert are 0.2 and 0.4 respectively.The Sahara was defined as 15 • W-35 • E and 11-33 • N.

Results and discussion
We separate our results into several distinct sections.Firstly we describe the results from our sensitivity experiments (Sect.3.1).Following this we explain the mechanism that is responsible for the reduction in poleward heat transport during the Holocene (Sect.3.2), and the role that sea-ice feedbacks play in the Arctic cooling we simulate (Sect.3.3).We then introduce the results of additional sensitivity experiments (Sect.3.4), which we performed to test the robustness of our initial findings.We follow this with a discussion of the uncertainties of our results (Sect.3.5).To finish, we discuss these findings in the context of the limitations of our model and how these could be overcome in future work (Sect.3.6)

Results
Our results for both the OG and OGGIS equilibrium experiments can be separated into two distinct phases, 9-6 and 6-0 ka (Fig. 1a and b).In the OG experiment there is a total cooling in the Arctic from 9 to 0 ka of 2.9 • C (Figs. 1a  and 2a), with 0.7 • C occurring between 9 and 6 ka, and 2.1 • C between 6 and 0 ka.Of this total cooling, 2.1 • C is due to direct orbital and greenhouse gas forcing (Figs.1b  and 2b), whilst 0.5 • C (17 %) is due to Sahara desertification alone (Figs. 1a and 2c).The overall decrease in Arctic mean annual temperatures in our OG experiments, from 9 to 0 ka, is consistent with those observed in a transient experiment performed with a previous version of LOVECLIM (Renssen et al., 2005), in which a 1-3 • C decrease in mean annual temperatures in the high northern latitudes was observed.In addition, our observed decrease from 6 to 0 ka, essentially the Holocene thermal maximum (HTM) to present day, of 2.1 • C, is in line with the proxy-based temperature reconstructions from Greenland ice cores (Dahl-Jensen et al., 1998), in which there is an observed decrease of 2.5 • C between the mid-Holocene climate optimum and present day.Our results also lie within the range of proxy-based temperature reconstructions for 16 terrestrial sites from the western Arctic, where it has been shown that temperatures during the Holocene Thermal Maximum (HTM) were on average 1.6 ± 0.8 • C warmer than 20th century temperatures (Kaufman et al., 2004), which is consistent with proxy-based temperature reconstructions for northern Eurasia (Renssen et al., 2009;Sundqvist et al., 2010).
In our OGGIS transient experiment, the 9 ka climate is relatively cold due to the cooling effect of the LIS topogra- phy and albedo, and the GIS albedo, which leads to a delayed thermal maximum over most of the Arctic domain.This delayed onset of the thermal maximum for the early to mid-Holocene has been observed in numerous proxy records (Kaufman et al., 2004) and is consistent with previous experiments performed with LOVECLIM.In particular Renssen et al. (2009), who simulated a ∼ 1 • C increase in mean annual temperatures north of 60 • between 9 and 6 ka in a transient experiment that included the same forcings as our 6 and 9 ka OGGIS equilibrium experiments.This delayed response of high-latitude warming is truly representative of our observation of a 1 • C increase in mean annual temperature from 9 to 6 ka (Fig. 1b).The simulated temperature increase in OGGIS between 9 and 6 ka is also present in temperature reconstructions from lakes in northern Iceland, in which a temperature increase of ∼ 0.75 • C was observed (Axford et al., 2007).Additionally, the change in Sahara vegetation had a cooling effect of 0.2 • C between 9 and 6 ka and without this moderating effect, the warming would have even been higher (1.3 • C).
The second phase, 6-0 ka, of the OGGIS equilibrium experiment was identical to the OG equilibrium experiment, with a total decrease in mean annual Arctic temperature of 2.1 • C, with 0.  records from the Arctic (Kaufman et al., 2004;Renssen et al., 2009;Sundqvist et al., 2010).
For 9 to 0 ka, the cooling due to Saharan desertification is 0.5 • C in OG, which is 17 % of the total cooling of 2.9 • C simulated with all forcings (Fig. 1a).In the OGGIS experiment, the Saharan desertification suppressed the warming by 15 % between 9 and 6 ka (i.e. by 1.0 • C instead of 1.3 • C), and from 6 to 0 ka, it was responsible for 19 % of the observed cooling in the Arctic (−0.4 • C vs. −2.1 • C, Fig. 1b).

Mechanisms connecting the Sahara to the Arctic
The model experiments indicate that part of the cooling in the Arctic over the Holocene is a direct consequence of Sahara desertification, which invokes a land-atmosphere teleconnection.Due to the prescribed desertification in the Sahara in our equilibrium simulations, over the course of the Holocene net albedo and radiative heat loss increases, leading to a de-crease in surface temperatures and an increase in surface pressure in the Sahara.The decreasing temperatures in the Sahara cause an increase in the surface pressure (Fig. 3) with an extension over the Tropical Atlantic, leading to an easterly zonal shift, plus an expansion, of the Azores high.Additionally, we observe a weakening of the low-latitude trade winds and a net overall decrease in the atmospheric meridional heat flux (Fig. 4).In particular we observe a weakening of the mid-latitude westerlies.This overall weakening of the winds and atmospheric heat transport over the Atlantic Ocean is consistent with a decrease in the meridional temperature gradient due the relatively strong cooling over the Sahara.As a result, the Icelandic Low stabilises, which in turns results in a weakening of the polar easterlies.The Bjerknes compensation theory suggests that weakened atmospheric heat transport would be compensated by an increase in oceanic heat transport (Bjerknes, 1964).This is indeed what we observe (Fig. 4), however the increase in oceanic heat transport (0.058 PW at 7 • N, before gradually reducing to 0 PW in the Arctic) is less than the decrease in total atmospheric heat transport (VQ + VT), resulting in an overall decrease of poleward heat transport.Mid-latitude storms are responsible for the majority of the heat transport from the equator to the Arctic (Peixoto and Oort, 1992;Zhang and Rossow et al., 1997).In our simulations, the weakening of the Northern Hemisphere winds are shown to be robust (Fig. 4), in particular the westerlies, which results in a reduced meridional atmospheric heat transport from the low latitudes to the Arctic (Fig. 4), leading to widespread cooling north of 60 • N.
In the OGGIS simulation, we see the same long-range land-atmosphere-ocean teleconnection present.However, the localised effects of the increase in albedo between 9 and 6 ka, due to the diminishing LIS and the localised effects of insolation changes, results in a warming over North America.This warming emanates eastwards over the North Atlantic and the rest of northern Europe and Eurasia, as well as penetrating the Arctic, accounting for the observed warming between 9 and 6 ka in the OGGIS experiment.

Contribution of sea-ice feedbacks to Arctic cooling
Since the peak in the Holocene thermal maximum (∼ 6-8 ka), the high northern latitudes cooled, sea-ice and snow cover increased, leading to an increase in surface albedo, which thus induced further cooling (Kerwin et al., 1999;Crucifix et al., 2002).In addition, Arctic cooling also resulted in thickening of the sea-ice cover during the Holocene, leading to a reduction in the ocean-to-atmosphere heat flux and relatively cooler conditions in the lower atmosphere through the ice-insulation feedback, especially in winter (Renssen et al., 2005).It is important to see to what extent these important positive feedbacks of the climate system contributed to the Arctic cooling initiated by the Saharan desertification.
Our results for the OG equilibrium simulations show that from 9 to 0 ka there is a reduction in sea-ice coverage of 25 %  in the Beaufort Sea, 15 % across the rest of the Canadian Basin and 12 % across the majority of the eastern Arctic and its shelves (Fig. 5a).However, this is a result of all of the forcings, orbital, greenhouse gases and vegetation changes in the Sahara.This indicates that LOVECLIM has captured this important feedback and is in line with previous results (Renssen et al., 2005).Of this observed reduction, it can be seen that the majority is a direct result of the combination of changes in orbital and greenhouse gases (Fig. 5b) and that a minor amount of sea-ice reduction (3 %) can be accredited to the perturbed Sahara vegetation fractions within our experiments (Fig. 5c).Additionally, we see a similar response in the sea-ice thickness (not shown).This clearly shows that the positive feedbacks associated with sea-ice play a modest role in the Arctic cooling that we accredited to the change in the vegetation in the Sahara (17 %).

Sensitivity experiments
To evaluate the full potential of the impact of Saharan vegetation changes on the Arctic climate, we performed seven additional sensitivity experiments in which we prescribed either 100 % grassland or 100 % desert cover in the Sahara (see Appendix D).Biome reconstructions of pollen and plant macrofossils show that during the mid-Holocene the Sahara was fully covered by grass and shrubs (Hoelzmann et al., 1998;Jolly et al., 1998).Hoelzmann et al. (1998) produced estimated percentages of various tree types across northern Africa and the Arabian Peninsula for 6 ka, showing steppe and savannah covering much of northern Africa (> 90 %).Claussen et al. (1999) were the first to replicate the results of Hoelzmann et al. (1998) and Jolly et al. (1998) in a climate model, successfully reproducing the drastic decrease in Sahara vegetation fraction over the Holocene from 90 to 0 %.However, in our OG simulation from 9 to 0 ka, the vegetation fraction decreases from 65 to 20 %, clearly showing that whilst LOVECLIM is capable of simulating the general pattern of vegetation change observed (Claussen et al., 1999), it is unable to adequately capture the full amplitude of this highly unstable period that characterised the termination of the AHP.
Therefore, this suggests that the impact of Sahara vegetation on Arctic cooling could be greater than the 17 % we estimated from our initial sensitivity experiments.To account for this, and to constrain our results, we simulated extreme, early (9 ka) and late (0 ka) Holocene environments.By "extreme" we imply that for 9 ka we replicated a 100 % grass-vegetated environment, and for 0 ka we replicated a 100 % desert environment within our model.This was done because if we had decided to follow the results of Jolly et al. (1998) and Claussen et al. (1999), that of a 90 % vegetated Sahara region at 9 ka and a 0 % vegetated Sahara region at 0 ka, we would have been left with several questions.Firstly, how would we allocate the 90 % vegetation at 9 ka, given that we would have to make the decision of doing so per grid cell or over the entire Sahara domain?If over the Sahara domain, which regions do we leave arid?Compounding the difficulty of this, it must be noted that there is no proxy data available for the spatial extent of vegetation at 9 ka.Therefore, we opted for 100 % vegetated cover at 9 ka.The results from these experiments enable us to place an upper limit of the potential impact of Sahara desertification has upon Arctic cooling over the Holocene.
The results show that from a 9 ka, 100 % "green" Sahara to a 0 ka, 100 % "desert" Sahara, temperatures decrease by 4.0 • C (Fig. 2d), of which 1.6 • C (40 %) is attributable to the change in vegetation.However, given the extreme vegetation scenarios we prescribed at 0 and 9 ka, care must be taken when interpreting this result.

Uncertainty in our results
One important requirement for our results is that the Sahara region was warmer at 9 and 6 ka than at 0 ka, as suggested by our simulations.However, it is not possible to confirm or refute this on the basis of field-based evidence due to the absence of proxy-based temperature reconstructions for the Sahara region during the early and mid-Holocene.Although our LOVECLIM result is consistent with GCM simulations for the mid-Holocene (Braconnot et al., 2007), it is important to critically address this aspect of our simulations, especially because the prescribed cloud cover in our model has likely significantly affected the surface energy balance over the Sahara.
In LOVECLIM, the total upward and downward shortwave fluxes are calculated as a function of the ISCCP D2 cloud cover data set consisting of modern day observations.Whilst this is a reasonable first-order approximation for the modern day Sahara, it is not so during periods when the Sahara was vegetated, i.e. 9 and 6 ka.This is because having a cloud cover that is essentially associated with a desert environment, as opposed to a vegetated environment, results in less clouds than there was in reality for 9 and 6 ka.Hence, the simulations at 9 and 6 ka ignore changes in cloud albedo.These changes in cloud albedo would approximately balance impacts from changes in surface albedo, which we do include within our original equilibrium simulations.
For instance, if we assume that desert albedo is 0.4, and we have a 100 % cloudless desert environment, then the downward solar radiation reaching the surface of the Sahara is 342 Wm −2 (Kiehl and Trenberth, 1997).Thus, the total solar radiation absorbed at the surface is 205 Wm −2 [(1 − 0.4) • 342].In a cloud covered, vegetated region, with surface albedo of 0.2, and assuming that 23 % of incoming solar radiation is reflected by clouds (Kiehl and Trenberth, 1997) then the downward solar radiation reaching the surface of the Sahara is 265 Wm −2 [342 * (1 − 0.23)].Therefore, the solar radiation absorbed at the surface is 211 Wm Hence, in theory and as a first-order approximation, the surface of a cloudy vegetated region and that of a cloudless desert environment absorb approximately similar amounts of incoming solar radiation (211 and 205 Wm −2 respectively).As a result, temperature differences seen between these states would depend on differences between the latent and sensible heat fluxes, with less moisture available, lower latent heat flux and high sensible heat flux.
To test this we devised a series of sensitivity experiments that would allow for the addition of cloud cover over the Sahara at 9 and 6 ka.To achieve this we took the modern cloud cover that is prescribed in our model for the Amazon region and applied it directly to the Sahara in both the 9 and 6 ka sensitivity simulations.Hence, artificially induced cloud cover over the Sahara at these time periods is based on the modern day cloud cover for the Amazon.Given that the Amazon region consists of trees as opposed to a vegetated Sahara that more than likely consisted of grasses and shrubs during its vegetated period in the mid-and late Holocene, this prescribed cloud cover is probably an over estimation, but still serves the purpose of assessing the impact of cloud cover over a vegetated Sahara on surface temperature.
As can be seen (Fig. 6), the Sahara was warmer during the mid-Holocene in both the simulation with the original modern day cloud data set at 6 ka (Fig. 6a) and also in the simulation with the appropriate prescribed cloud cover for a vegetated Sahara at 6 ka (Fig. 6b).The difference between these two simulations, in essence the impact of the difference between the two cloud cover set ups we employed, is shown in Fig. 6c.Here, it can be seen that with the modern day cloud cover, the surface temperature in the Sahara is warmer by 0.7-1.3• C.This results in the Arctic being 0.3-0.7 • C warmer due to the modern day cloud data set over the Sahara.
As a final test of the sensitivity experiments, when we compare the temperature changes in the Arctic at 9 ka between simulations 9k9kEQ_OG and 9k9kEQ_OG_clouds, we see that the mean annual Arctic temperatures are −11.0±0.8 • and −11.1 ± 0.8 • C respectively (see Appendix D).Therefore, based on these sensitivity experiments, we can say that a major premise of our original experiments, that the Sahara was warmer at 9 and 6 ka, compared to 0 ka holds, and thus, the land-atmosphere teleconnection between the Arctic and Sahara that we previously introduced is robust.

Limitations of LOVECLIM
The results presented here provide a first-order approach to understanding the impact of long-term radiative forcing of the Sahara over the course of the Holocene (9-0 ka).Whilst we have good confidence in the mechanisms driving the land-atmosphere teleconnection we observe over the Atlantic region, it is important to highlight the limitations associated with our approach.
The land-atmosphere teleconnection we have outlined is initiated by a change in vegetation cover in the Sahara region and then carried through the atmosphere from the equatorial region to the Arctic.However, this draws close scrutiny of the atmospheric component of LOVECLIM.Given the relative simplicity of the ECBilt model compared to GCMs, it would be easy to dismiss our findings.However, as the model uses relatively simple physics, any findings in such a model should also be detectable in a more complex model.To highlight this, we shall compare the response of LOVECLIM in other studies to more complex models.Specifically we shall highlight studies that have focused on interglacial climate only as this will allow us to show the response of LOVE-CLIM to different perturbations, over a variety of temporal and spatial extents.The primary focus of our study is the cooling of the Arctic over the Holocene.Therefore, it would be apt for us to investigate how the temperature profile of LOVECLIM compares with more complex models.In this study and previous studies (Renssen et al., 2005), it is shown that LOVECLIM shows a gradual cooling over the Holocene as a response to reduced summer insolation forcing and decreasing greenhouse gases.In a recent study (Liu et al., 2014), the primary focus of which was trying to understand the apparent discrepancy between a reconstructed global temperature study (Marcott et al., 2013) and models, it was shown that the three models employed in the study, LOVECLIM, CCSM3 and FAMOUS (the latter two being GCMs) all simulated a robust annual mean global warming (0.5 • C) throughout the Holocene.
In a study of the last interglacial, Nikolova et al. ( 2013) provide a detailed analysis of a variety of climatic parameters, as simulated in both LOVECLIM and CCSM3.In these experiments they perturbed the climate for forcings equivalent to MIS-5e (127 ka).The results of these experiments show that the JJA surface temperature evolution during MIS-5e for both models is an exceptionally good match.This is especially true in the tropics over northern Africa, Greenland and North America.In the Arctic (defined as 60-90 • N), the simulated warming is 3 and 2.2 • C for LOVECLIM and CCSM3, respectively, which is in agreement with another model-data comparison (Otto-Bliesner et al., 2006).In addition, they show that the JJA precipitation between the two models to also be of a reasonable good fit, with both models simulating an intensification of both the African and Asian monsoons in line with proxy reconstructions.In addition, both models are consistent in replicating increased low-level wind speed and moisture transport over Africa during this period, as well a stronger upper-level (200 hPa) tropical easterly jet over Africa (Nikolova et al., 2013).
In a multimodel-data study (Bakker and Renssen, 2014), which included six GCMs and three EMICs, the data showed that when comparing the simulated period of maximum warmth during the last interglacial period against proxybased temperature reconstructions, LOVECLIM, along with the two other EMICs in the study, all overestimated the maximum temperature the least, across all latitudes.This is probably on account of models with lower resolution displaying less internal variability, when compared to GCMs (Gregory et al., 2005), leading to a more spatially coherent evolution of the climate.
In another series of sensitivity experiments performed for the period MIS-13 (∼ 500 ka), LOVECLIM was used to in-vestigate the East Asian monsoon activity during this period (Qiuzhen Yin et al., 2008).The results of this study were conferred by Muri et al. (2013), who replicated these results with HadCM3, a fully coupled atmosphere-ocean general circulation model.
As we have outlined above when compared to GCMs, LOVECLIM has shown itself to be able to replicate similar results across a range of interglacial climates.Although just a selection are included here, it serves to demonstrate that despite its apparent shortcomings, LOVECLIM is an effective tool for climate investigations over a variety of spatial and temporal scales.
Another limitation in our study is due to the simplistic cloud parameterisation scheme that is employed within our model.As discussed in Sect.3.4, we designed a series of sensitivity experiments to test this; although this is not a perfect solution, we were able to show that our original landatmosphere teleconnection holds.However, without the evolution of clouds in response to changes in the climate system, particularly in response to surface temperature and sea-ice coverage changes, we are still missing potentially important feedbacks.However, this issue is not easily solved as it is a problem that persists through the entire spectrum of the modelling community, and it is widely accepted that this is a major source of uncertainty in climate model sensitivity experiments (Randall et al., 2007).
Another potential source of error in our results lies in the vegetation fractions simulated in our transient simulations, which were then used to define the Sahara vegetation fractions in our equilibrium experiments.As previously noted (Sect.3.3), the vegetation fraction in our model decreases from 65 to 20 % from 9 to 0 ka, as opposed to reconstructions and modelling studies which indicate an approximate 90-0 % decrease over the same period (Claussen et al., 1999;Hoelzmann et al., 1998;Jolly et al., 1998).
As a result of this, our 9 ka simulations (with a 65 % tree fraction at 9 ka, as opposed to a 90 % tree fraction) have a higher surface albedo than proxy reconstructions and modelling studies indicate, and thus simulate a lower surface temperature at 9 ka.Additionally, if we consider our 0 ka simulations (with a 20 % tree fraction, as opposed to a 0 % tree fraction) results in a lower surface albedo, and thus a higher surface temperature at 0 ka.
As a result, the temperature gradient between the early and late Holocene in our experiments is not as steep as it possibly would be.Therefore, it would be interesting to see in future experiments, with a model that captures a more realistic temperature gradient between the early and late Holocene, how the polewards transported heat flux over the course of the Holocene is impacted.
To summarise, compared to comprehensive GCMs, LOVECLIM has several simplifications, of which the fixed cloud cover and the relatively simple vegetation scheme are likely to be the most important for this study.However, model inter-comparison studies have shown that LOVE-

579
CLIM produces similar interglacial climate responses to changes in forcings as GCMs, making it likely that our results will be confirmed in future studies with more comprehensive models.

Concluding remarks
In conclusion, our simulation experiments indicate that over the course of the Holocene, the observed cooling in the Arctic region is not only driven by local changes in insolation, but that the effects of desertification in the Sahara initiate a long-range land-atmosphere teleconnection.In our model experiments, this teleconnection accounts for between 17 and 40 % of the observed Arctic cooling between 9 and 0 ka, and it is likely that the actual effect is nearer the upper end of this range.We also show through a series of sensitivity experiments that despite an apparent weakness of our model, with its prescribed modern day cloud data set, our overall findings are conclusive and robust, withstanding the sensitivity experiments we performed.Therefore, our results indicate that high-latitude interglacial climate is sensitive to large-scale changes in Saharan vegetation cover.0 ka 0 ka EQ 0k0kEQ_OGGIS −13.9 0.9 0 ka 6 ka EQ 0k6kEQ_OGGIS −13.5 0.9 0 ka 9 ka EQ 0k9kEQ_OGGIS −13.2 1.0 6 ka 0 ka EQ 6k0kEQ_OGGIS −12.2 0.8 6 ka 6 ka EQ 6k6kEQ_OGGIS −11.8 0.8 6 ka 9 ka EQ 6k9kEQ_OGGIS −11.5 0.8 9 ka 0 ka EQ 9k0kEQ_OGGIS −13.5 0.9 9 ka 6 ka EQ 9k6kEQ_OGGIS −13.0 0.9 9 ka 9 ka EQ 9k9kEQ_OGGIS −12.8 0.9

Figure 1 .
Figure 1.Simulated temperature change in the Arctic for (a) OG and (b) OGGIS equilibrium simulations, showing the relative contributions of different forcings; ORB + GHG (orbital and greenhouse gases), VEG S (vegetation changes in the Sahara), OTHER (other factors outside the Sahara region such as vegetation changes), ORB + GHG + LISTOPO (orbital and greenhouse gases and Laurentide Icesheet Topography, which is only relevant for the period 9-6 ka) and ALL (for OG this includes GHG, ORB and prescribed Sahara vegetation; for OGGIS this includes GHG, ORB, prescribed Sahara vegetation, LIS melt, GIS melt and LIS topography changes).Error bars represent ±1σ .
4 • C due to desertification in the Sahara and the remainder of the cooling due to the localised effects of insolation changes.As previously discussed, the decrease in mean annual temperatures from the mid-to late Holocene are representative of a proxy temperature reconstruction from Greenland ice cores (Dahl-Jensen et al., 1998) and terrestrial

Figure 3 .
Figure 3. 9k0kEQ-9k9kEQ_OG 800 hPa geopotential height, allowing the effects of changing vegetation in the Sahara, on pressure, to be visualised.