Interactive comment on “ On the effect of orbital forcing on mid-Pliocene climate , vegetation and ice sheets ”

Abstract. We present results from modelling of the mid-Pliocene warm period (3.3–3 million years ago) using the Earth system model of intermediate complexity CLIMBER-2 analysing the effect of changes in boundary conditions as well as of orbital forcing on climate. First we performed equilibrium experiments following the PlioMIP (Pliocene Model Intercomparison Project) protocol with a CO2 concentration of 405 ppm, reconstructed mid-Pliocene orography and vegetation and a present-day orbital configuration. Simulated global Pliocene warming is about 2.5 °C, fully consistent with results of atmosphere–ocean general circulation model simulations performed for the same modelling setup. A factor separation analysis attributes 1.5 °C warming to CO2, 0.3 °C to orography, 0.2 °C to ice sheets and 0.4 °C to vegetation. Transient simulations for the entire mid-Pliocene warm period with time-dependent orbital forcing as well as interactive ice sheets and vegetation give a global warming varying within the range 1.9–2.8 °C. Ice sheet and vegetation feedbacks in synergy act as amplifiers of the orbital forcing, transforming seasonal insolation variations into an annual mean temperature signal. The effect of orbital forcing is more significant at high latitudes, especially during boreal summer, when the warming over land varies in the wide range from 0 to 10 °C. The modelled ice-sheet extent and vegetation distribution also show significant temporal variations. Modelled and reconstructed data for Northern Hemisphere sea-surface temperatures and vegetation distribution show the best agreement if the reconstructions are assumed to be representative for the warmest periods during the orbital cycles. This suggests that low-resolution Pliocene palaeoclimate reconstructions can reflect not only the impact of increased CO2 concentrations and topography changes but also the effect of orbital forcing. Therefore, the climate (Earth system) sensitivity estimates from Pliocene reconstructions which do not account for the effect of orbital forcing can be biased toward high values.


Introduction
The time period between 3.3 and 3 million years ago (Ma), traditionally named the mid-Pliocene warm period (MPWP), is the most recent interval in Earth's history which is believed to have been significantly warmer than the pre-industrial climate.Most palaeoproxies suggest that atmospheric CO 2 concentration during the MPWP was in the range from 300 to 400 ppm (Raymo et al., 1996;Pagani et al., 2009;Seki et al., 2010;Bartoli et al., 2011), the global surface air temperature was higher by 2-3 • C (estimate based on model simulations) (Haywood and Valdes, 2004;Jansen et al., 2007;Lunt et al., 2012b) and sea level was up to 25 m above present level during the warmest periods of the MPWP (Miller et al., 2012;Foster and Rohling, 2013).If these estimates are sufficiently accurate, this makes the MPWP a useful (although by no means direct) analogue for the future greenhouse world (Salzmann et al., 2009;Haywood et al., 2011b).
A comprehensive set of reconstructions of climate and boundary conditions for the MPWP is available from the PRISM3D (Pliocene Research, Interpretation and Synoptic Mapping) dataset (Dowsett et al., 2010).PRISM3D includes global reconstructions for sea-surface temperature (SST), vegetation distribution, ice-sheet extent and orography.Most of the data are the outcome of a combination of reconstructed data and modelling results.SST reconstructions show little change in the tropics relative to present day, while significant Published by Copernicus Publications on behalf of the European Geosciences Union.
warming is observed at mid-to high latitudes, up to 3-4 • C in the Northern Hemisphere (NH) with major warming occurring in the North Atlantic (Dowsett et al., 2009;Dowsett and Robinson, 2009).
There are numerous pieces of evidence from palaeodata indicating that the spatial distribution of vegetation during the MPWP was substantially different from today.The main difference was a northward expansion of taiga replacing tundra in the NH.The coast of Antarctica was at least locally covered by some vegetation.The vegetation distribution in mid and low latitudes was essentially similar to present day, although the moister climate favoured forest and woodland growth in some regions now covered by steppe or grassland (Thompson and Fleming, 1996).The PRISM3D vegetation reconstruction is based on a combined model-data approach and is described in Salzmann et al. (2008).It combines an internally consistent dataset of reconstructions at 202 palaeobotanical sites with results of simulations of the BIOME4 vegetation model driven by modelled mid-Pliocene climatology.Where data coverage is poor, the vegetation reconstruction is mostly based on the model output.Some important regions with low data availability include high northern Eurasia and North American mid-latitudes.This should be kept in mind when comparing model output with reconstructions.More applications of the BIOME4 model to the MPWP can be found in e.g.Haywood et al. (2002b) and Salzmann et al. (2009).Simulations with the BIOME4 model in offline mode do not account for the feedback of the vegetation on climate.MPWP simulations with a coupled climate-vegetation model generally show a reasonably good agreement with reconstructed vegetation (PRISM2), but regionally important differences exist, with the model underestimating the northward expansion of coniferous forest and overestimating the extension of arid areas (Haywood and Valdes, 2006).
MPWP sea level has been estimated to be up to 25 m higher than present during the warmest periods of the MPWP corresponding to the warm orbital configurations, which implies significant reduction of the Greenland Ice Sheet and the Antarctic Ice Sheet (Pollard and DeConto, 2009;Miller et al., 2012;Foster and Rohling, 2013) 1 .PRISM3D spatially explicit estimation of ice-sheet extent is primarily based on modelling approaches (Hill et al., 2007).Time-slice modelling experiments with different MPWP orbital configurations suggest that the Greenland Ice Sheet varied from almost ice-free to nearly present-day state, while Antarctic ice volume showed variations up to 30 % (Dolan et al., 2011).PRISM3D topography, which we use in this study, is more similar to present day than that in PRISM2, with small elevation reductions in the Andes and the Rocky Mountains (Rowley and Garzione, 2007).Because of the coarse spatial resolution of the CLIMBER-2 model, we first aggregate anomalies of MPWP orography from the PRISM3D dataset 1 We will loosely refer to orbital configurations with higher global temperatures as "warm".
to the CLIMBER grid and then add them to present-day model orography.
Some attempts were made to use the MPWP to derive the so-called Earth system sensitivity (ESS), i.e. the long-term quasi-equilibrium response of global temperature to radiative forcing.In contrast to the climate sensitivity (CS, also called Charney sensitivity) which involves only fast feedbacks, ESS includes also the slow feedbacks related to ice sheets, vegetation and atmospheric aerosols.The ESS estimates from MPWP range from modest 4.5 • C (Lunt et al., 2009) to very high 7.1-8.7 • C (Pagani et al., 2009) values.This large range can be explained by uncertainties in reconstructed CO 2 concentrations and model-derived global temperatures.Similar attempts have been made to estimate the dependence of global mean sea level on temperature or CO 2 (e.g.Foster and Rohling, 2013).Obviously, any attempt to employ MPWP data to constrain climate models and climate sensitivity depends not only on the accuracy of palaeoclimatic reconstructions but also on the validity of the implicit assumption that reconstructed MPWP climate represents quasi-equilibrium response to elevated level of CO 2 and modified topography.Testing this assumption has important implications for the interpretation of MPWP reconstructions.Even if climatic conditions during the MPWP were relatively stable compared to the more recent past, significant variability is evident in the benthic oxygen isotope (δ 18 O) record (Lisiecki and Raymo, 2005), suggesting a nonnegligible response of the climate system to orbital forcing.PRISM3D SSTs are well known to represent warmest conditions during the MPWP because of the warm peak averaging procedure involved in the reconstruction (Dowsett, 2007).The temporal resolution of vegetation proxies, mostly fossil pollen, is very coarse, and consequently vegetation reconstructions do not resolve orbital variations (Salzmann et al., 2008;Thompson and Fleming, 1996).Where information on vegetation variability in time is available, the biome representing the warmest climatic conditions has been included in the PRISM3D reconstructions (Haywood et al., 2013b).Furthermore, in regions with poor data coverage PRISM3D vegetation is heavily based on the BIOME4 simulation driven by climatology from an atmospheric general circulation model (AGCM) with prescribed PRISM2 boundary conditions, including SSTs (Salzmann et al., 2008).Since reconstructed SSTs represent warm MPWP conditions, it is possible that vegetation reconstructions are also warm-biased in regions with scarce data availability.
Up to now most simulations of MPWP climate represent equilibrium simulations with atmosphere general circulation models with prescribed SSTs (e.g.Chandler et al., 1994;Sloan et al., 1996;Haywood et al., 2000;Jiang et al., 2005;Haywood et al., 2009) or with coupled atmosphere-ocean climate models (e.g.Haywood and Valdes, 2004;Haywood et al., 2007;Lunt et al., 2007).In both cases, vegetation distribution and ice-sheet extent were prescribed.The most recent coordinated efforts as a part of PMIP3/CMIP5 (Paleoclimate Model Intercomparison Project/Coupled Model Intercomparison Project) is PlioMIP (Pliocene Model Intercomparison Project) (Haywood et al., 2010(Haywood et al., , 2011a)).In the case of PlioMIP, the protocol includes present-day orbital parameters and CO 2 concentration of 405 ppm with other greenhouse gases (GHGs) at the preindustrial level, as well as orography, ice sheets and vegetation from PRISM3D reconstructions.Results of eight models participating in PlioMIP reveal a wide range of global temperature responses to identical forcings and boundary conditions (1.8-3.6 • C) and a similarly wide range (3.5-6.8 • C) for the implied ESS, with some models having an ESS almost identical to the CS, while others have an ESS twice as large as the CS (Haywood et al., 2013a).A typical feature of most simulations of the MPWP with coupled climate models is an overestimation of warming at low latitudes and an underestimation of warming at high latitudes compared to available palaeoclimate reconstructions.The discrepancies between data and models in tropical temperatures are present not only in simulations of the MPWP but also other palaeoclimate simulations both for warm and cold climates (Lunt et al., 2012a).The reasons for such systematic difference remain debatable.Weaker simulated MPWP warming at high latitudes may have several reasons.One is that climate models underestimate the strength of climate feedbacks in the high latitudes.One such feedback is the climate-vegetation feedback, which in coupled climate models is accounted via prescribed MPWP vegetation.Another possibility is that the standard PlioMIP setup does not account for some orographic changes such as changes in ocean bathymetry which can affect ocean circulation (Robinson et al., 2011).Lastly, the low temporal resolution of most of Pliocene palaeoclimate records makes it possible that some of the reconstructions are warm-biased and represent climate conditions corresponding to warm orbital configurations rather than mean conditions over the entire 3.3-3 Ma period.
The effect of orbital forcing has previously been analysed with time-slice experiments using an AGCM coupled to a slab ocean (Haywood et al., 2002a).For different orbital configurations, they find small annual changes in SSTs and land surface temperatures, but seasonally more significant variability due to insolation distribution changes.Further, they show only minor variations in vegetation cover in offline experiments with the BIOME4 model.However, their model does not include interactive vegetation and ice sheets, which can amplify the direct effect of orbital forcing.
So far the role of orbital forcing on MPWP climate, vegetation and ice sheets has not been analysed with a fully coupled model.Here we present results from transient simulations using an Earth system model driven by time-dependent orbital forcing including experiments with interactive vegetation and ice sheets.

Model
For our analysis we use the Earth system model of intermediate complexity CLIMBER-2 (Petoukhov et al., 2000;Ganopolski et al., 2001).CLIMBER-2 includes a 2.5dimensional statistical-dynamical atmosphere model, a 3basin zonally averaged ocean model including sea ice, a dynamic model of the terrestrial vegetation (Brovkin et al., 1997(Brovkin et al., , 2002) ) and the three-dimensional polythermal ice-sheet model SICOPOLIS based on the shallow-ice approximation (Greve, 1997).The latter is only applied for the Northern Hemisphere because it does not include an ice shelf model that is crucial for simulation of the Antarctic Ice Sheet.Although this is a limitation of our modelling approach, it is obvious that variations of the Antarctic Ice Sheet on the orbital timescales had a rather localised effect on climate.The most relevant processes contributing to the response of climate to orbital forcing take place in the NH because of the larger land area where ice sheets can expand and vegetation can change.Different components of the model have different spatial resolution.Atmospheric and land components have latitudinal resolution of 10 • and longitudinal of approximately 51 • .The ocean model has a latitudinal resolution of 2.5 • and 21 levels vertically.The ice-sheet model has a latitudinal resolution of 0.75 • , longitudinal resolution of 1.5 • and vertically 20 levels in sigma coordinates.Atmosphere and ice sheets are coupled bi-directionally using a physically based energy balance approach (Calov et al., 2005).Surface runoff and ice calving simulated by the ice-sheet model are applied to the ocean as surface freshwater flux.Since CO 2 concentration was prescribed, the marine carbon cycle model was not used in this study.The CLIMBER-2 model in different configurations has been used for numerous studies of past and future climates (e.g.Ganopolski et al., 1998Ganopolski et al., , 2001;;Claussen et al., 2006;Schneider von Deimling et al., 2006;Ganopolski and Calov, 2011;Brovkin et al., 2012;Eby et al., 2013).

MPWP boundary conditions
In all MPWP experiments (equilibrium and transient) we prescribed global orography according to the PRISM3D dataset.Since our ice-sheet model is only applied to the Northern Hemisphere, in all experiments the Antarctic ice-sheet extent and elevation were also prescribed according to PRISM3D.Following the PlioMIP protocol (Haywood et al., 2010), the CO 2 concentration was set to 405 ppm in all experiments.The 405 ppm value is at the high end of concentration estimates (Pagani et al., 2009;Seki et al., 2010), but is chosen to account for unknown additional contributions from non-CO 2 greenhouse gases such as methane.In the set of equilibrium simulations which followed the PlioMIP protocol, we also prescribed the Greenland Ice Sheet and vegetation according to PRISM3D reconstructions (Dowsett et al., 2010).PRISM3D data are adjusted and aggregated to the CLIMBER-2 resolution.For vegetation we used the PRISM3D product in its megabiome version.Because in CLIMBER-2 only forest, grass and deserts are distinguished, we aggregated the megabiome data to these three classes at the model resolution (Appendix A).As a final result we obtained the fraction of forest, grass and desert cover in each grid cell.Additionally, annual mean values of the leaf area index are derived from the vegetation cover in a semiinteractive vegetation model setup with fixed cover fractions.Orbital parameters in PlioMIP experiments were set to their present-day values.In addition to the PlioMIP-type of experiment, we performed a set of transient experiments through the entire MPWP period with time-dependent orbital forcing computed according to Laskar et al. (2004).In these experiments the vegetation cover and the NH ice sheets were either interactive or prescribed according to PRISM3D (see below).

PlioMIP experiment and factor separation analysis
Since CLIMBER-2 is a model of intermediate complexity, to enhance credibility of modelling results, it is useful to compare its results with available results of more-complex climate models for identical experiments.Since transient MPWP experiments have not yet been undertaken with coupled AOGCMs, we performed CLIMBER-2 experiments using the PlioMIP protocol.In addition, to disentangle the effects of CO 2 concentration, orography, vegetation and ice sheets on the mid-Pliocene climate, a full factor separation is performed following the technique described in Lunt et al. (2012b).For four factors this requires running 2 4 = 16 different experiments, with all possible combinations of changes in the boundary conditions.The total global annual MPWP temperature change will be a combination of the contribution from the single factors.The way the single contributions are computed is described in Appendix B.

Transient experiments
To estimate the effect of orbital forcing on the MPWP climate, ice sheets and vegetation we perform a suite of transient experiments with time-dependent orbital forcing.With respect to both ice sheets and vegetation we performed three types of experiments: the first one with prescribed present-day ice sheets (vegetation), the second with prescribed PRISM3D ice sheets (vegetation) and the third with interactive ice sheets (vegetation).

Results and Discussion
First the results from the factor separation analysis are presented, then those from the transient experiments.

PlioMIP experiment and factor separation
The simulated total global annual mean mid-Pliocene warming compared to pre-industrial climate is 2.5 • C (see also Fig. 1).This value is fully consistent with the 1.8-3.6 • C range obtained in PlioMIP (Haywood et al., 2013a).The main contribution, 1.5 • C, comes from the higher CO 2 concentration.Additional significant warming comes from changes in vegetation (0.4 • C), orography (0.3 • C) and ice sheets (0.2 • C) (Table 1).The synergy between the factors is small, 0.02 • C.This implies that the total response of the system is to a good approximation a linear superposition of the contributions of the single factors.While the effects of CO 2 and vegetation are similar to those found by Lunt et al. (2012b), the effects of orography and ice sheets in our experiments are lower.At least some of the discrepancies might arise from the fact that Lunt et al. (2012b) used PRISM2, while we use PRISM3D reconstructions.In some regions differences in orography between PRISM2 and PRISM3D are significant with the latter more similar to present day.Together with the much coarser CLIMBER-2 resolution, this could at least partly explain the smaller effect of orography and ice sheets in our study compared to Lunt et al. (2012b).In general, some differences are possibly due to different climate sensitivities in the two models (3.04 • C for HadCM3, 2.7 • C for CLIMBER-2).
Total annual global precipitation increases by 0.27 mm day −1 in the PlioMIP simulation compared to the pre-industrial climate.This is above the range from PlioMIP models, which give a precipitation increase in the range 0.09-0.18mm day −1 .Of the total precipitation change, 0.16 mm day −1 are attributable to CO 2 , 0.06 to vegetation, 0.03 to orography and 0.01 to ice sheets (Table 1).As expected, the relative contribution from the single factors is comparable to the proportions found for global temperature, except for vegetation, which has a larger impact, and ice sheets, which have a smaller impact on precipitation than on temperature.The global annual temperature and precipitation factor separation is summarised in Table 1.
Temperature and precipitation anomalies simulated for the PlioMIP modelling setup are illustrated in Figs.1a and  2a, respectively.Our spatial patterns for temperature compare favourably well with the model mean from PlioMIP (Haywood et al., 2013a).In particular, the largest warming, around 10 • C, is modelled around Greenland and over parts of Antarctica.Warming above 5 • C is simulated in the northern parts of Eurasia and North America.Also, temperature anomalies in CLIMBER-2 are consistent with PlioMIP results also in areas affected by changes in orography, such as central Asia and South America.The Sahara is 3-4 • C warmer than in the pre-industrial climate.Annual precipitation shows a large increase over the tropics, particularly over South America and around Australia.Drier climate is simulated over continental Eurasia and over parts of the South Atlantic and the central Pacific.In general, compared to the  PlioMIP models, CLIMBER-2 simulates a larger increase in annual precipitation.In particular, no reduction in precipitation is simulated over parts of South and Central America and the dipole in precipitation anomalies over the tropical Pacific is reversed compared to mean PlioMIP results (Haywood et al., 2013a).However, the PlioMIP precipitation anomalies show a large spread between the models, particularly in the tropics (Haywood et al., 2013a).
The geographical distribution of the separate effects of CO 2 , orography, ice sheets and vegetation on annual mean temperature and precipitation are shown in Figs. 1 and 2 together with the synergy between the factors.The main contribution to the total signal is from the higher CO 2 concentration, while the other factors can be equally important but only on a regional scale.Ice sheets have a very local impact limited to Greenland and Antarctica.Regarding temperature, vegetation has a significant effect at high latitudes and over an extended region around northern Africa, Europe and central Asia.All factors are important for the warming in polar regions.Orography changes are also important in Southeast Asia and South America.The largest changes in precipitation occur in the tropics and have significant contributions from CO 2 , orography and vegetation.The synergy term is small for both temperature and precipitation.Our spatial patterns for the annual temperature factor separation compare favourably well with those from Lunt et al. (2012b) except for orography, where discrepancies are attributable to different versions of PRISM datasets used.As already pointed out for global mean values, also regional differences exist regarding the amplitude of the warming, with results from Lunt et al. (2012b) showing higher temperature anomalies.
We additionally performed an experiment with PRISM3D boundary conditions but interactive vegetation.Simulated global temperature in the experiment with interactive vegetation is 0.2 • C lower compared to the experiment with PRISM3D vegetation, i.e. with interactive vegetation and present-day orbital configuration the model can explain only half of the warming attributed to vegetation by the factor separation analysis.A comparison of the modelled forest cover with its PRISM3D counterpart is presented in Fig. 6c,d.Overall the model shows lower forest fractions compared to the reconstructions.The discrepancy is particularly evident at high northern latitudes, where modelled forest area is significantly lower than suggested by MPWP reconstruction.Also, in South America and in central Asia, PRISM3D shows more widespread forest cover than our model.Over central South America, although CO 2 and orographic changes would lead to an increase in precipitation (Fig. 2c,e), it is not enough to promote forest growth in the model.As will be discussed in the next section, some of the discrepancies between the model and reconstructions can be reconciled if we consider the effect of orbital forcing, particularly at high latitudes.

Transient runs
The time series of global mean annual temperature differences from pre-industrial climate for different transient experiments of the MPWP are presented in Fig. 3b the orbital timescales is dominated by the precessional component of orbital forcing, as was shown already in Claussen et al. (2006).Interestingly enough, the precessional component is only clearly seen during the last 100 kyr (1 kyr: 1000 yr) of the MPWP.To clarify the origin of such nonlinear response to the precessional component of the orbital forcing, we performed an additional experiment in which the fraction of vegetation was fixed in the Sahara/Sahel region (between 10 and 30 • N) but interactive elsewhere.Fixing Sahara/Sahel vegetation leads to significant weakening of precessional variability, which unambiguously points to surface albedo in Sahara/Sahel region as the main cause for the precessional component of global temperature variability.Since parameterisation of surface albedo in the CLIMBER-2 model is rather crude, it is possible that the effect of vegetation on surface albedo is overestimated.On the other hand, the "greening" of Sahara would reduce the atmospheric dust load and, as a result, lower the planetary albedo and warm the NH.This effect is not accounted for in PlioMIP-type simulations.
In addition to the precessional cycle, global mean temperature also contains a pronounced 100 kyr cycle.This cycle is not related to the vegetation, because the same variability is observed also in a transient simulation with fixed vegetation (not shown).This 100 kyr cyclicity is explained by corresponding eccentricity variations (Fig. 3a), which causes global annual mean insolation variations of up to 0.5 W m −2 (Paillard, 2001).
Global temperature in simulations with interactive ice sheets but fixed vegetation is also dominated by the precessional component.This is related to a strong response of the NH ice sheets to the summer 2 insolation (Fig. 4).During periods of high insolation, large portions of the Greenland Ice 2 When we refer to summer, we hereinafter always mean boreal summer.Sheet disappeared, while during cold summers Greenland Ice Sheet occupies the entire island, and a relatively small ice sheet grows over the Arctic Archipelago.As it was shown in Ganopolski and Calov (2011), high levels of CO 2 prevent Earth from falling into the full glacial state.Previous studies have already addressed the effect of orbital forcing on ice sheets in offline time-slice experiments of the MPWP showing very different ice-sheet states under different orbital configurations (Dolan et al., 2011).Considerable variations in NH ice volume are also modelled in our coupled transient experiments (Fig. 4).Most of the variations come from the Greenland Ice Sheet, but an expansion of ice caps over North America and Scandinavia is simulated during the coldest orbital intervals.The magnitude of NH ice-sheet volume variations corresponds to approximately 15 m sea level equivalent, which is not in conflict with magnitude of benthic δ 18 O variability (Lisiecki and Raymo, 2005).However, modelled icevolume variability is dominated by precession, while the benthic δ 18 O record is dominated by obliquity.This fact can be at least partly explained by a compensation of precessional variability in the ice volume between the two hemispheres (Raymo et al., 2006;Dolan et al., 2011).
The combined effect of interactive ice sheets and vegetation results in further amplification of the precessional forcing.Global mean modelled warming varies now within the even larger range 1.9-2.8• C. Again, the best agreement with T PlioMIP is found in warm orbital intervals.More detailed information can be inferred from zonal temperature anomalies (Fig. 5).The most striking feature is the large temperature variability in the fully coupled transient experiment both over land and ocean.The variability is particularly significant in summer over land where the temperature anomalies span a 10 • C interval, with temperatures ranging from values similar to or lower than present day to up to 12 • C higher.This has fundamental implications for the interpretation of reconstructions over both ocean and land, since dating accuracy and temporal resolution of MPWP proxy records are generally not high enough to resolve orbital timescales.Again, the zonal temperatures for the PlioMIP-like experiment are shown for comparison.The mean MPWP zonal temperatures of the fully coupled transient run are generally lower than in the PlioMIP simulation, particularly in the NH, with annual mean differences up to 0.5-1 • C.However, the PlioMIP temperatures fall well into the variability range of the transient experiment.
The SST anomalies are also compared with the PRISM3D data.The model generally underestimates the warming in mid-to high northern latitudes, while it overestimates the temperature increase in the tropics, where PRISM3D SSTs are essentially close to present-day values.Including orbital forcing improves agreement between model and data assuming that PRISM3D SSTs represents the warmest annual conditions during the entire MPWP.The agreement is even better if reconstructed temperatures are biased towards a warm season.In fact, the summer SST increase in the NH estimated from proxy data falls well into the wide range of SSTs modelled throughout the MPWP, except above 80 • N.
Finally, we analyse the impact of orbital forcing on simulated vegetation distribution, and compare these with PRISM3D reconstructions.Orbital forcing has a significant effect on simulated vegetation.It causes large temporal variations in vegetation cover, particularly at high northern latitudes.Substantial forest area variations are simulated north of 60 • N (Fig. 4) with oscillations in time closely following the maximum summer insolation at 65 • N. The corresponding value derived from PRISM3D is also plotted.The best agreement between modelled and reconstructed NH vegetation is found during warm orbital intervals, although the simulated forest area always stays below the PRISM3D value.Figure 6a and b illustrate the geographical distribution of maximum and minimum simulated forest fraction during the whole MPWP.Apart from the large variability at high northern latitudes, significant variations are observed also over Asia.The modelled equilibrium forest cover with present-day orbital parameters is very similar to the time-averaged forest fraction from the transient experiment (not shown).Simulated (CLIMBER-2) and reconstructed (PRISM3D) zonal mean forest fraction anomalies relative to present day are compared in Fig. 7.The main feature in the PRISM3D data is a large increase in forest at high northern latitudes and in the subtropics, while minor reductions are found in NH mid-latitudes.CLIMBER-2 under Pliocene boundary conditions but with present-day orbital configuration simulates small changes in forest cover compared to the pre-industrial state everywhere except above 60 • N, where forest expands at the expense of tundra.However, the large increase in forest fraction between 70 and 80 • N suggested by palaeoclimate reconstruction cannot be explained by the higher CO 2 concentration alone.Here, a generally better agreement between model and data is found in the transient runs during warm orbital periods (Fig. 7).These results are apparently in contrast with those from Haywood et al. (2002b), who show small variations in vegetation cover in time-slice experiments with different orbital configurations.

Conclusions
Results  signal.Regionally, and especially seasonally, the effect of orbital forcing is comparable or even more important than higher CO 2 and orographic changes.The most significant effect of orbital forcing is found in summer over NH land, where temperatures vary by up to 10 • C during the 300 kyr interval.Simulated response of vegetation and ice sheets to the orbital forcing is much larger than to change in CO 2 concentration alone.Only a combination of elevated CO 2 and "optimal" orbital configuration leads to SST and vegetation cover being more consistent with reconstructions.
In summary, we find that orbital forcing plays a fundamental role in controlling modelled MPWP climate and consequently vegetation distribution and ice-sheet extent.If the MPWP reconstructions are warm-biased and therefore contain orbital fingerprints, then a large portion of discrepancies between data and the model can be reconciled.This is the case for NH SSTs, boreal vegetation cover and ice-sheet extent, suggesting that at least some of the proxy data may reflect not only the effect of increased CO 2 concentration and changed topography but also the effect of orbital forcing.This has wider implications for the ongoing attempts to estimate climate/Earth system sensitivity with a combined model-reconstructions approach.Our results imply that Earth system sensitivity estimates based on MPWP reconstructions are likely to be biased to high values.As has already been proposed elsewhere (Haywood et al., 2013b), a possible solution to this problem is to use reconstructions from a time slice with close to present-day orbital configuration instead of the whole MPWP slab, although it will be a challenge to obtain enough data to reconstruct a mid-Pliocene time slice, especially for vegetation.
Our modelling results do not support a notion that the Earth system sensitivity is significantly higher than the classical climate (or Charney) sensitivity.Indeed, results of model simulations (Dolan et al., 2011) in conjunction with the existence of a strong temporal variability of the benthic δ 18 O on orbital timescales suggest that a significant fraction of the ice sheet's reduction during the warmest intervals of MPWP was related to orbital forcing rather than elevated CO 2 .As far as vegetation feedback is concerned, it amplifies CO 2 -induced temperature change by only 0.2 • C in our model (less than 10 % of the climate sensitivity), which is consistent with previous studies (e.g.O'ishi et al., 2009).Admittedly, there are climate feedbacks which are not included in our model and which can be important for amplification of CO 2 induced warming, such as aerosols (e.g.aeolian dust), atmospheric chemistry, etc.However, the fact that our model is able to simulate the MPWP climate anomaly in reasonable agreement with reconstructions can be interpreted as evidence against the existence of very strong and unaccounted for climate feedbacks.

Vegetation aggregation
The PRISM3D megabiome vegetation dataset distinguishes nine plant functional types (PFTs), while CLIMBER-2 differentiates only between forest and grassland.In CLIMBER-2 each grid cell on land is composed of a mosaic of forest, grass and desert with fractional areas computed by the vegetation module.However, in PRISM3D a unique PFT or desert is assigned to each grid cell.According to the PRISM3D -CLIMBER-2 transfer function between surface types in Table B1, we transform the PRISM3D from 10 to 3 surface types keeping the 2 • × 2 • resolution.Then for each 10 • × 51 • CLIMBER-2 cell we count the number of corresponding 2 • × 2 • cells covered by forest, grass and desert covering the same area, and from that compute the forest, grass and desert share which is then used as boundary condition for CLIMBER-2.

Appendix B Factor separation
Adopting the notation from Lunt et al. (2012b) we indicate with a subscript the boundary conditions that are prescribed from PRISM3D (c for CO 2 , o for orography, i for ice sheets and v for vegetation).E.g.T iv is the temperature from the experiment with prescribed PRISM3D ice and vegetation, but with present-day CO 2 and orography.Assuming that the single factors are independent from each other and including  The synergy between the factors is given by S = T ciov − T + dT CO 2 + dT oro + dT ice + dT veg . (B6) It is simply the portion of the total signal that cannot be explained by the linear superposition of the contributions from single factors.
Fig. 3. (a) Maximum summer insolation at 65 • N computed from Laskar et al. (2004) (black solid line) and present-day value (black, dashed) and eccentricity computed from Laskar et al. (2004) (grey solid line) and present-day value (grey, dashed).(b) Annual global temperature anomalies (MPWP pre-industrial) for the fully coupled transient experiment (solid red), transient experiment with fixed PRISM3D vegetation and interactive ice sheets (blue), transient experiment with fixed PRISM3D ice sheets and interactive vegetation with fixed (dashed green) and interactive Sahara/Sahel (solid green).The PlioMIP experiment equilibrium warming is shown for comparison (dashed black).
Fig. 4. (a) Maximum summer insolation at 65 • N computed from Laskar et al. (2004).(b) Total NH ice-sheet volume (blue) and forest area between 60 and 80 • N (green) for the fully coupled experiment (solid lines).Dashed green lines represent forest area between 60 and 80 • N computed from PRISM3D reconstructions.

Fig. 5 .
Fig. 5. Zonal mean temperature anomalies (MPWP-pre-industrial).Annual (a) and summer (c) surface air temperatures over land, and annual (b) and summer (d) sea-surface temperatures for the fully coupled transient simulation, showing the time average over the entire MPWP (solid red line) as well as the temperature variability range (shaded red) and the PlioMIP experiment equilibrium temperature anomalies (dashed black).The blue solid lines represent PRISM3D SST anomalies.

Fig. 6 .
Fig. 6.Forest fractions.Modelled maximum (a) and minimum (b) forest fraction (excluding glacier area) in each grid cell during the MPWP for the fully coupled transient simulation.(c) Modelled equilibrium forest fraction for the PlioMIP experiment with interactive vegetation.(d) Forest fraction derived from PRISM3D reconstructions (as described in Appendix A).

Fig. 7 .
Fig. 7. Zonal forest fraction anomalies (MPWP-present day) for PRISM3D (blue), for the fully coupled transient experiment at the time of maximum summer insolation at 65 • N (red) and for the PlioMIP equilibrium simulation with interactive vegetation (green).

Table 1 .
Annual global temperature and precipitation factor separation.The sum of the single contributions and the total changes including the synergy term are shown as well.