Orbital control on late Miocene climate and the North African monsoon: Insight from an ensemble of sub-precessional simulations. Past

. Orbital forcing is a key climate driver over multi-millennial timescales. In particular, monsoon systems are thought to be driven by orbital cyclicity, especially by precession. Here, we analyse the impact of orbital forcing on global climate with a particular focus on the North African monsoon, by carrying out an ensemble of 22 equally spaced (one every 1000 years) atmosphere–ocean–vegetation simulations using the HadCM3L model, covering one full late Miocene precession-driven insolation cycle with varying obliquity (between 6.568 and 6.589 Ma). The simulations only differ in their prescribed orbital parameters, which vary realistically for the selected time period. We have also carried out two modern-orbit control experiments, one with late Miocene and one with present-day palaeogeography, and two additional sensitivity experiments for the orbital extremes with varying CO 2 forcing. Our results highlight the high sensitivity of the North African summer monsoon to orbital forcing, with strongly intensiﬁed precipitation during the precession minimum, leading to a northward penetration of vegetation up to ∼ 21 ◦ N. The modelled summer monsoon is also moderately sensitive to palaeogeography changes, but it has a low sensitivity to atmospheric CO 2 concentration between 280 and 400 ppm. Our simulations allow us to explore the climatic response to orbital forcing not only for the precession extremes but also on sub-precessional timescales. We demonstrate the importance of including orbital variability in model–data comparison studies, because doing so partially reduces the mismatch between the late Miocene terrestrial proxy record and model results. Failure to include orbital variability could also lead to signiﬁcant miscorrelations in late


Introduction
Late Miocene (11.61-5.33 Ma;Hilgen et al., 2005;Gradstein et al., 2004) climate is thought to have been globally warmer and wetter than the present day, as indicated by the available proxy reconstructions and modelling studies (e.g. Bradshaw et al., 2012aBradshaw et al., , 2015Bruch et al., , 2011Eronen et al., 2010Eronen et al., , 2011Pound et al., 2011Utescher et al., 2011). It is suggested that the Antarctic Ice Sheet was already present throughout this time period (e.g. Lewis et al., 2008;Shackleton and Kennett, 1975) while the presence of a much reduced Greenland Ice Sheet means that the Northern Hemisphere was nearly ice-free (Kamikuri et al., 2007;Moran et al., 2006). This period was also characterised by significant tectonic reorganisation, such as the gradual closure of the Panama Gateway (e.g. Duque-Caro, 1990;Keigwin, 1982) and the major uplift of modern mountain chains (e.g. Himalayas, Molnar et al., 1993;Andes, Garzione et al., 2000;Alps, Kuhlemann, 2007;East African Rift, Yemane et al., 1985; Rockies, Morgan and Swanberg, 1985). There is considerable uncertainty in the reconstructed CO 2 concen-trations for the late Miocene, with values ranging from about 140 to 1400 ppm, but with most of the data converging between the preindustrial (280 ppm) and present-day (400 ppm) concentrations (see Fig. 1 in Bradshaw et al., 2012a, and references therein). Several recent studies suggest that CO 2 concentrations at this time were closer to the upper end of this range (e.g. Zhang et al., 2013;Bolton and Stoll, 2013;LaRiviere et al., 2012). This is, however, still a matter of debate; in fact, model-data comparisons based on vegetation reconstructions argue that the presence of European seasonal temperate forests in the late Miocene (as indicated by the fossil record) is consistent with preindustrial CO 2 concentrations (Forrest et al., 2015).
The quantitative global proxy record for the late Miocene is spatially biased, with the majority of the terrestrial data coming from the Mediterranean region. From a global perspective, the available biome reconstructions suggest a very different vegetation distribution from that of the present day, for instance with the presence of boreal and temperate forests at higher northern latitudes (Pound et al., 2011. A previous modelling study highlighted the mismatch between the available proxy reconstructions and climate simulations for this period (Bradshaw et al., 2012a), exhibiting significantly lower surface air temperatures in the model than in the palaeo-reconstructions. These authors carried out an extensive terrestrial model-data comparison and demonstrated that part of the proxy-derived temperature and precipitation differences between the late Miocene and preindustrial reference climates could be explained by changes in the palaeogeography (a combination of marine gateways, continental position and ice extent). In addition, reconstructed surface air temperatures could only be accounted for by assuming that late Miocene CO 2 concentrations were towards the higher end of the range of estimates (400 ppm). The palaeosimulations by Bradshaw et al. (2012a) were carried out using present-day orbital forcing, whereas the compilation of proxy reconstructions used in their study is likely to have been generated under a range of different stages of the orbital cycle and climate states. Here, we partially explain the extent of the remaining model-data mismatch for the late Miocene by carrying out the same analysis as Bradshaw et al. (2012a) while taking into account not only changes in palaeogeography and CO 2 concentrations but also orbital variability on sub-precessional timescales.
As well as holding the majority of the terrestrial proxy data for the late Miocene, the Mediterranean region also contains proxy data with an exceptionally high temporal resolution, commonly at precessional resolution or higher. The Mediterranean's geological record throughout the Neogene is characterised by regular alternations which have been interpreted as a sedimentary response to orbital forcing (e.g. Krijgsman et al., 2001;Sierro et al., 2001;Hilgen et al., 1999). The mechanism is thought to be increased freshwater input as a consequence of enhanced runoff into the basin, causing both stratification of the water column and  (2004) orbital solution. Obliquity (green), eccentricity (blue) and climatic precession (red), which is defined as e sin , where is the longitude of perihelion and e is eccentricity. Bottom: experimental design for the set of 22 late Miocene orbital simulations with 280 ppm atmospheric CO 2 concentrations. Simulations are spanned 1 kyr apart throughout this precession cycle. Each simulation is indicated by a number (1 to 22) and all simulations are designed based on the precession cycle, but orbital parameters all vary at the same time, as shown by the dotted line for experiment 22. The precession maximum experiment is indicated as pMAX and the precession minimum as pMIN. Obliquity is expressed both in radians and degrees, and insolation in W m −2 . enhanced surface productivity at times of high summer insolation. This would lead to the deposition of organic-rich sediments known as sapropels (Kidd et al., 1978), which are preserved by the anoxic conditions on the seafloor. Today the main source of fresh water to the Mediterranean Sea is the Nile River, whose discharge is driven by summer monsoonal rainfall. During the late Miocene, the structure of the North African catchment area is thought to have been different, with the extensive north-central catchment draining into the eastern Mediterranean Sea via the Chad-Eosahabi River (Griffin, 2002(Griffin, , 2006 rather than into the Niger River as it does today. With both the Chad and Nile catchments draining a large area of North Africa that is influenced by the summer monsoon, the seasonal variability of the freshwater input into the Mediterranean basin may have been greater than today. This could have caused significant shifts in the Mediterranean Sea's hydrologic budget at times of intensified monsoon. In addition, the Mediterranean region underwent significant palaeogeographic changes during the late Miocene, driven by the motion of the African and Eurasian plates. In particular, the tectonic evolution of the Gibraltar Arc during the Messinian (7. 25-5.33 Ma) led to reducedand at times absent -exchange between the Mediterranean Sea and the Atlantic Ocean (Hsu et al., 1973), which triggered widespread changes in Mediterranean sea level and salinity . The late Miocene, therefore, represents an ideal time period in terms of considering sub-precessional scale processes and North African monsoon dynamics, given their implications for the Mediterranean environment. The North African region is of particular interest as it has also been previously identified as an optimal target for late Miocene palaeodata acquisition (Bradshaw et al., 2015;Lunt et al., 2008b).
A correlation between changes in Earth's orbital parameters and the North African monsoons has previously been suggested in several modelling studies (e.g. Bosmans et al., 2015a;Tuenter et al., 2003;Prell and Kutzbach, 1987) and sedimentary analyses (e.g. Larrasoaña et al., 2003;Lourens et al., 2001;Rossignol-Strick and Planchais, 1989). According to both models and observations, the strength of the African summer monsoon is enhanced at times of minimum climatic precession (hereafter, climatic precession will be referred to simply as precession, which is defined in the caption of Fig. 1) and, to a lesser extent, of maximum obliquity. This enhanced strength is a result of the increased amplitude of the seasonal cycle of solar radiation in the Northern Hemisphere, which increases the land-ocean temperature contrast. The variability of the African monsoon as a result of changes in orbital forcing is commonly studied using idealised numerical experiments, mainly considering the extremes of the orbital cycle (precession/insolation maximum and minimum). In this study, we analyse changes in the intensity and seasonality of the North African summer monsoon on sub-precessional timescales throughout a full late Miocene precession-driven insolation cycle with varying obliquity (hereafter, this insolation cycle will be simply referred to as a late Miocene precession cycle, in line with the experimental design of Sect. 2.2) and explore the phasing between orbital forcing and climatic responses not only for the precession extremes but throughout the entire precession cycle.

Numerical model
The experiments presented in this study were carried out using a global general circulation model (GCM), the UK Hadley Centre Coupled Model (HadCM3L, version 4.5), a coupled atmosphere-ocean GCM with a horizontal resolution of 3.75 • longitude by 2.5 • latitude for both the atmosphere and ocean components. There are 19 vertical levels in the atmosphere and 20 levels throughout the ocean (Cox et al., 2000). The resolution of this GCM is typical for palaeoclimate studies because it allows the computation of long integrations, from centuries to millennia, and of numerous ensemble members. The model has also been used in several other pre-Quaternary palaeoclimate studies, both for the late Miocene (Bradshaw et al., 2012a(Bradshaw et al., , 2015 and Eocene (Loptson et al., 2014;Lunt et al., 2010;Tindall et al., 2010). Other late Miocene simulations have also been carried out running the higher-resolution-ocean (1.25 • latitude and longitude) version of the model (Ivanovic et al., 2013(Ivanovic et al., , 2014a. However, here we used the lower-resolution and more computationally efficient HadCM3L because of the availability of an existing 2000-year spin-up, for consistency with the Bradshaw et al. (2012a) study, and because of the number of simulations conducted in the ensemble.
HadCM3L is coupled to the dynamic global vegetation model TRIFFID (Top-down Representation of Interactive Foliage and Flora Including Dynamics; Hughes et al., 2004;Cox, 2001), which can simulate five plant functional types (PFTs): broadleaf and needleleaf trees, C3 and C4 grasses, and shrubs. Land surface processes are simulated by the MOSES-2.1 (Met Office Surface Exchange Scheme) land surface scheme (Essery and Clark, 2003), which includes nine surface types (the five PFTs plus those representing bare soil, lakes, ice and urban surfaces). Previous studies highlighted the importance of including land surface processes and vegetation to simulate the warm conditions inferred from the late Miocene palaeorecord, especially with relatively low CO 2 concentrations (e.g. Bradshaw et al., 2015Bradshaw et al., , 2012aKnorr et al., 2011). Here, TRIFFID is run in equilibrium mode, which has a ∼ 5-year coupling period and efficiently produces equilibrium states even for the slowest variables (Cox, 2001).
The late Miocene palaeogeography used in our experiments is the same as Bradshaw et al. (2012a), which is characterised by significant reductions in the elevation of most of the world's highest mountain chains compared to the present day (e.g. lower Tibetan Plateau and Andes) and by a much smaller extent of the Greenland Ice Sheet. These late Miocene orography and boundary conditions are based on the reconstructions by Markwick (2007); the full technique is described in Markwick (2007) and Markwick and Valdes (2004). Other significant differences from the present-day continental configuration in our late Miocene simulations are the more southerly position of Australia, the closed Bering Strait, the open Panama Gateway and non-restricted Indonesian Seaway, and the presence of the Barents/Kara Sea landmass (Supplement Fig. S1). A full description of the model setup and results of preindustrial and late Miocene control experiments with both 280 and 400 ppm CO 2 concentrations are provided by Bradshaw et al. (2012a), including an assessment of the model's performance with respect to modern observations (see Bradshaw et al., 2012a; Appendix B, Sect. 1.1).

Experimental design
Modelling studies exploring the impact of orbital forcing on climate and monsoon systems are traditionally performed as idealised sensitivity experiments, mainly simulating the most extreme configurations of the orbital cycle (e.g. precession maxima and minima; Bosmans et al., 2012, 2015a, Braconnot et al., 2008Tuenter et al., 2003;Kutzbach,1981). In contrast, our ensemble spans a full orbital cycle so that no assumptions about which phase of the orbit will be more (or less) extreme for a particular variable or region are made. This setup allows us to analyse the evolution of the global climate system (e.g. monsoons) on sub-precessional timescales and to investigate the importance of orbital variability when evaluating the mismatch between proxy reconstructions and model results for the late Miocene. We carried out 22 equally spaced simulations (one every 1000 years) throughout a full late Miocene precession cycle, between 6.568 and 6.589 Ma. This time period was chosen because of its relatively high eccentricity values, which enhance the precession-induced climatic signal (Fig. 1a). Another study which used realistically varying orbital parameters to assess the impact of orbital forcing on climate variability has recently been carried out by Prescott et al. (2014) for the mid-Pliocene warm period (∼ 3.3 to 3 million years ago). Prescott et al. (2014) used a lower temporal resolution (one simulation every 2 or 4 kyr) but a conceptually similar experimental design. The high temporal resolution of our ensemble for an entire late Miocene precession cycle also allows direct comparison with micropalaeontological and geochemical data from the Mediterranean Sea for the same time period (e.g. Perez-Folgado et al., 2003;Sierro et al., 2001). The Mediterranean model-data comparison on sub-precessional timescales will be explored in a future study.
The initial model integration for the core orbital ensemble is taken from Bradshaw et al. (2012a). Each one of the orbital simulations begins from the end of their 2000-year integration at 280 ppm CO 2 with a present-day orbital configuration and a late Miocene palaeogeography. The trend in the global mean temperature for this simulation is very small: < 8 × 10 −4 • C per century (Bradshaw et al., 2012a). Choosing 280 ppm as the baseline rather than 400 ppm means that a comparison can be made between the effect of varying orbital parameters and increasing CO 2 , to address the cold temperature bias in late Miocene simulations with respect to proxy reconstructions (e.g. Bradshaw et al., 2012aBradshaw et al., , 2015Knorr et al., 2011;Micheels et al., 2007;Steppuhn et al., 2006). All orbital parameters were changed for each simulation, and they were derived from the Laskar et al. (2004) orbital solution. Each ensemble member has been run for 200 years, and here we analyse the climatological means of the last 50 years of simulation. The deep and intermediate ocean has not reached equilibrium by the end of our simulations, but as we investigate relatively short-term atmospheric processes, this is not expected to influence our analysis greatly. This approach is consistent with that used by Bosmans et al. (2015a), who ran their experiments for 100 model years and did not find strong trends in surface air temperatures and precipitation. In addition, the climate system was found to be in equilibrium for the discussed atmospheric variables in the transient orbital experiments performed with an earth system model of intermediate complexity, which also justifies the use of snap-shot simulations from more complex models (Tuenter et al., 2005). Trends for global annual surface air temperatures and summer precipitation over North Africa, after 200 years of simulation, are shown for four of our experiments in the Supplement (Fig. S2). The complete experimental design for the main orbital ensemble is shown in Fig. 1b.
For the presentation of our results we use a modern-day calendar. This does not take into account the changes in the length of the seasons determined by variations in the date of perihelion along a precession cycle (Kutzbach and Gallimore, 1988;Joussaume and Braconnot, 1997). This socalled "calendar effect" has the potential to introduce biases in the seasonal interpretation of the analysed variables. The use of a modern-day calendar is not uncommon for this kind of climate simulations, but its impact has been evaluated. In their assessment of this issue, Prescott et al. (2014) found that seasonal surface air temperatures were not influenced by the calendar effect. Chen et al. (2011) also found only minor seasonal biases in the responses of temperature and precipitation to precessional forcing of the main global summer monsoon systems. The use of a variable celestial calendar is more important when carrying out seasonal comparisons between model results and the proxy record (Chen et al., 2011). Here, however, the calendar effect will not significantly influence the terrestrial model-data comparison (Fig. 8) and the analysis of the phase relationship between surface air temperatures and precession forcing (Fig. 6), in which we only consider maximum and minimum values of the analysed variables, regardless of the season or month in which these occur. After testing the effect of a calendar-effect correction on our results (see Figs. S3 and S4), we choose to use a modernday calendar in our analysis to facilitate comparison with other modelling studies (e.g. Bosmans et al., 2012Bosmans et al., , 2015aPrescott et al., 2014;Braconnot et al., 2008;Tuenter et al., 2003Tuenter et al., , 2005. Two 200-year control simulations have also been run, one with a preindustrial setup (PIctrl) and one for the late Miocene period (LMctrl) but using present-day orbital configurations, both using the same setup as the simulations in Bradshaw et al. (2012a). In the main ensemble of 22 orbital simulations, the Strait of Gibraltar has been kept closed in order to simulate the significantly reduced or missing Mediterranean-Atlantic exchange of the latest Messinian. Given the shortcomings in the parameterisation of Mediterranean-Atlantic exchange for the HadCM3 model (Ivanovic et al., 2013), we decided to close the gateway as neither solution (open or close) would have been entirely realistic. The Mediterranean-Atlantic connection was not yet fully restricted during our simulated time period, but nonetheless it was different compared to today's exchange through the Strait of Gibraltar. Due to the uncertainty in late Miocene CO 2 concentrations, two additional sensitivity experiments have also been run. These are equivalent to the precession extreme experiments at 280 ppm CO 2 (pMIN and pMAX), but the CO 2 concentration was increased to 400 ppm (pMIN400 and pMAX400). The CO 2 sensitivity experiments were run for 200 years and begin from a late Miocene 1500-year integration, also run at 400 ppm CO 2 , which was spun off from a 500-year late Miocene 280 ppm CO 2 integration (Bradshaw et al., 2012a). Running the full ensemble at both 280 and 400 CO 2 , or other intermediate concentrations, would have been too computationally expensive.
For the model-data comparison we only consider the Messinian part of the data set compiled by Bradshaw et al. (2012a), as our time period is within this time period. We have incorporated additional precipitation data reconstructed by Eronen et al. (2012) for the North American continent, which were not included in the original data set. An updated version of the full data set used in this study for the terrestrial model-data comparison, including the additional palaeo-precipitation reconstructions, is provided in the Supplement (Table S1). Note that the correct database queries are used in this work, following the correction applied in Bradshaw et al. (2012b), but the data used are the same as in Bradshaw et al. (2012a).
The output of these simulations can be accessed online: www.bridge.bris.ac.uk/resources/simulations.

Mean climate state
In this section we explore the mean climate state, both throughout the full precession cycle and in more detail for the two precessional extremes, both with 280 and 400 ppm CO 2 concentrations. Note that the occurrence of precession minimum relatively close to the obliquity maximum (as seen in Fig. 1b) is likely to enhance seasonality in the Northern Hemisphere (Tuenter et al., 2005).

Surface air temperature response through a full precession cycle
The evolution of the modelled global annual mean surface air temperature (SAT) throughout the full precession cycle is not significantly influenced by the varying orbital parameters in either of the two hemispheres ( Fig. 2a), with changes within ∼ 0.3 • C. It also does not appear to be related to the evolution of the annual mean incoming shortwave radiation, perpendicular to the Earth's surface, at the top of the atmosphere (from now on referred to as insolation), which is the same in the Northern and Southern hemispheres (Fig. 2a). Variations in global mean annual insolation are only a result of eccentricity, while precession and obliquity have no impact on it. In the model, the mean global SAT is the result of a combination of the two hemispheres balancing each other out, with neither of the two clearly dominating the trend (Fig. 2a). Finally, global mean SATs in the Southern Hemisphere are generally higher (by ∼ 2 • C) than in the Northern Hemisphere in our late Miocene simulations. The present-day configuration is the opposite of this, with the Northern Hemisphere on average 1.5 • C warmer than the Southern Hemisphere (Feulner et al., 2013, and references therein). This difference is caused by the open Panama Seaway in our late Miocene simulations (Lunt et al., 2008a), leading to a weaker Atlantic Meridional Overturning Circulation in the late Miocene, compared with that of the present day. This maintains warmer temperatures in the Southern Hemisphere and colder in the Northern Hemisphere in our palaeosimulations, as opposed to modern temperatures.
In the model, seasonal temperature variations are driven by orbital forcing and related to changes in insolation, which in turn exhibit opposite phasing between the two hemispheres in every season . In addition, in both hemispheres the seasonal cycles cancel each other out in pairs and therefore produce only small variations in the annual mean (as seen in Fig. 2a). In winter, SAT in the Northern Hemisphere is in phase with local insolation (in anti-phase with precession) but with a lead of ∼ 2 kyr (Fig. 2b). The same lead can be seen in the Southern Hemisphere but with the difference that SAT is in anti-phase with insolation (in phase with precession; Fig. 2b). In summer, modelled SAT is in phase with insolation with a lead of ∼ 3 kyr, both in the Northern and Southern hemispheres (Fig. 2c). In the model, the difference between maximum summer temperatures and minimum winter temperatures in the Southern Hemisphere is much lower (∼ 7 • C) than in the Northern Hemisphere (∼ 20 • C), due to the more extended presence of land in the Northern Hemisphere. A warmer winter season results in a higher annual mean temperature in the Southern Hemisphere than in the Northern Hemisphere (as seen in Fig. 2a).
In spring, modelled SAT is in phase with insolation in both hemispheres, with a lead of ∼ 2 kyr in the Northern Hemisphere and of ∼ 2.5 kyr in the Southern Hemisphere (Fig. 2d). In autumn, the Northern Hemisphere SAT is in Table 1. Number of overlaps and total number of Messinian palaeodata considered. Numbers are shown for both the experiments described in this study and those from Bradshaw et al. (2012a). Note that in Bradshaw et al. (2012a) the five precipitation data points in North America from Eronen et al. (2012) are not included because they were published at a later date. phase with insolation and it leads it by 3 to 4 kyr (Fig. 2e).
The same phasing can be observed in the Southern Hemisphere, with the difference that SAT leads insolation by about 3 kyr (Fig. 2e). Note, however, that in our simulations SAT is influenced not only by precession but also by the varying obliquity, which changes from a maximum of ∼ 23.9 • to a minimum of ∼ 22.8 • . In contrast to precession, obliquity will have the same effect on both hemispheres, with maximum obliquity resulting in stronger summer insolation and weaker winter insolation, thus leading to an enhanced seasonal contrast.
Using idealised orbital transient simulations with an earth system model of intermediate complexity, Tuenter et al. (2005) have shown that a simple mechanism can explain how the leads and lags in the climatic response (e.g. surface air temperatures) to insolation within a year can result in leads or lags in time with respect to orbital parameters. However, the same simulations performed with interactive vegetation indicated how quickly this mechanism can be extensively modified, because of vegetation feedbacks and changes in the albedo, as induced by different sea ice distributions (Tuenter et al., 2005). The simulations analysed in this study are carried out using a more complex general circulation model and using interactive vegetation. As such, understanding the causes of leads and lags in the climate system is sometimes challenging and the simple mechanism described by Tuenter et al. (2005) is no longer applicable. However, we investigate in more detail the vegetation responses to changes in orbital forcing for the North African monsoon region (Sect. 3.4.1), as well as the spatio-temporal phasing of surface air temperatures globally (Sect. 3.2).

Global climate response to orbital forcing: precession extremes
In the model, the DJF SAT anomalies between precession minimum and maximum (pMIN-pMAX) are generally negative (i.e. cooler; Fig. 3a), especially in north-central Asia and India (except adjacent to the Tibetan Plateau), central-east North America, and most of Australia. Maximum cooling is found around the Sea of Japan, which is probably a model artifact caused by the enclosed nature of the basin in the model, which intensifies the signal. The Nordic Seas exhibit much higher DJF temperatures during the precession minimum. Yin and Berger (2012) interpreted the substantial winter warming of the Arctic as a response to insolation forcing, resulting from the "summer remnant effect". This effect has also been discussed by Lunt et al. (2013) and , who attributed the warming to delays in sea ice formation during the winter season as a result of excess solar radiation during the summer months. In this study, however, the location of the warm anomaly is shifted further south, localised in the Nordic Seas area rather than in the Arctic. This is possibly a consequence of using a late Miocene palaeogeography, rather than the modern palaeogeographies applied by Lunt et al. (2013) and , and of the different sea ice distribution in our experiments. In these regions, the main differences are the presence of the Kara/Barents Sea landmass in our late Miocene simulations (Fig. S1) and the higher concentration of sea ice in the Nordic Seas and in the subpolar North Atlantic in general, compared to the preindustrial (Fig. S5). The modelled JJA mean SAT differences exhibit globally warmer temperatures during precession minimum (Fig. 3b), especially on land in the Northern Hemisphere, with maximum warming (positive anomalies) over central Eurasia. The exception is the cooling (negative anomalies) that occurs over the monsoon regions in North Africa and India, due to the intensified cloud cover as a consequence of enhanced monsoonal precipitation (Braconnot et al., 2007). Maximum warming is generally centred over the main land masses rather than the ocean, because of the ocean's greater heat capacity and potential for latent cooling. However, the North Atlantic also shows significantly higher SATs (up to 5.5 • C) at times of precession minimum.
The cold and warm month means (which at each model grid square represent the SAT for the coldest and warmest months, respectively) exhibit clear differences in the sign of the anomaly in each hemisphere (Fig. 3c, d). This represents the opposite effect of precession on both hemispheres. In the Northern Hemisphere, the cold month mean (Fig. 3c) mirrors DJF values (Fig. 3a), with reduced warming in the Nordic Seas. The warm month mean (Fig. 3d) largely mirrors JJA values (Fig. 3b), except the intensified warming in the North Atlantic and the differences in the Northern Hemisphere monsoon regions, where cooling is no longer visible.
The Southern Hemisphere's cold month mean (Fig. 3c) anomalies are mainly positive during precession minimum, especially in northern Australia and in the Southern Ocean along the Antarctic coast. Negative anomalies dominate the warm month mean (Fig. 3d), with the exception of northerncentral South America and part of central-southern Africa, as a result of vegetation changes modifying the albedo feedback. The mean annual temperature difference (Fig. 3e) is characterised by maximum warming in the Nordic Seas and part of the Arctic regions during the precession minimum, whereas cooling is found in the African and Indian monsoon belts (and off the coast of South America, around 10 • S), with the most negative values in the sensitive area around the Sea of Japan (which could be the result of a model artifact, as previously discussed). Changes are generally small elsewhere (mostly within ∼ 1.5 • C) and this lack of a clear signal is mainly due to the positive and negative forcing during summer and winter seasons balancing each other out. The colder (negative anomalies) mean annual temperatures over the monsoon regions are caused by the dominant JJA cooling, and the mean annual warmer (positive anomalies) values in the Arctic and Nordic Seas are a result of the overall warming in both seasons.
Absolute plots of modelled JJA precipitation at precession maximum ( Fig. 4a) and minimum (Fig. 4b) portray the distribution of enhanced precipitation in the equatorial regions and clearly highlight the northward shift of the inter-tropical convergence zone (ITCZ) during precession minimum, which is most clearly seen over the monsoon regions. Coupled climate models are typically affected by a split of the ITCZ over the tropical West Pacific Ocean (e.g. Bellucci et al., 2010;Oueslati and Bellon, 2015), which leads to large disagreement between models and observations for present-day simulations. This is also clearly visible in our simulations, especially in the absolute plots of precipitation (Fig. 4a, b).
Modelled precipitation in DJF shows small differences between the two precession extremes at high latitudes in both hemispheres (Fig. 4d). Prominent features include changes in the North Atlantic storm tracks which take a more southerly route during precession minimum, leading to the widespread spatial precipitation anomaly which extends over the Mediterranean Sea and south-west Europe (Fig. 4d). Most of the significant changes in precipitation patterns between the two precession extremes, both in DJF and JJA, are found around the location of the ITCZ, depicting its migration between the two hemispheres in response to changes in orbital forcing (Fig. 4c, d). In JJA the ITCZ shifts northward, towards the warmer Northern Hemisphere as a result of the higher insolation forcing in summer. This can be clearly identified in the monsoon regions, especially in Africa and Asia, which experience much higher summer precipitation (more than 3.5 mm day −1 increase) during precession minimum (Fig. 4c). In JJA wetter (positive anomalies; up to 1.5 mm day −1 ) conditions during precession minimum are also found north of ∼ 50 • N in the Northern Hemisphere, as well as across the Southern Ocean and over most of Australia in the Southern Hemisphere (Fig. 4c). In contrast, sig- nificant negative anomalies (up to 3.5 mm day −1 ) dominate the North Pacific, North America and the North Atlantic between ∼ 10 and 40 • N. Finally, precipitation anomalies are small in Antarctica and across the Arctic regions because of the reduced amount of precipitation over these areas.

Global climate sensitivity to atmospheric CO 2 concentrations
In addition to the full set of 22 simulations with preindustrial CO 2 concentrations (280 ppm), two sensitivity experiments were carried out at 400 ppm for the two precessional extremes (pMIN400 and pMAX400), in order to explore the global and local climatic response to varying CO 2 concen-trations and to assess the impact of CO 2 on the model-data comparison. This was necessary because of the uncertainty in reconstructed CO 2 concentrations for the late Miocene, as discussed in Sect. 1. In line with previous modelling studies (e.g. Bradshaw et al., 2012aBradshaw et al., , 2015, late Miocene climate warms significantly as CO 2 concentrations increase, especially at high latitudes, and the greatest warming is found on land in the Northern Hemisphere (not shown). For the precession minimum simulation the mean annual global SAT is 14.6 • C at 280 ppm while at 400 ppm it is 17 • C. Modelled SATs in DJF (Fig. 5a) are more sensitive to orbital changes (where the "orbital sensitivity" described here is the difference between precession minimum and maximum) at 400 ppm CO 2 (over 5 • C) in the  subpolar North Atlantic south-east of Greenland, in the regions around the Sea of Japan and in the north-west Pacific Ocean. In contrast, some areas reveal DJF SATs with increased sensitivity at 280 ppm (Fig. 5a), including the regions north of India, Canada and other parts of North America, central Africa and most of South America.
The most extended significant SAT anomalies in JJA (Fig. 5b) are found in the subpolar North Atlantic, Nordic Seas and the Arctic Ocean, where modelled temperatures are more sensitive to orbital changes at 400 ppm. This is also true for the North Pacific and the area around the Sea of Japan. Higher sensitivity at 280 ppm for JJA SATs is found locally in the Southern Ocean, especially around the Ross Sea, and in the Northern Hemisphere over Greenland (up to 4 • C). In both seasons, differences in and around the polar regions are most likely to be linked to changes in sea ice distribution. In fact, relatively warm initial conditions for the late Miocene simulations lead to enhanced sea ice loss during precession minima, triggering a strong positive sea ice feedback mechanism as CO 2 concentrations increase (Bradshaw et al., 2015).
In the model, the precipitation response to increasing CO 2 is a moderate increase at mid to high latitudes in both hemispheres (not shown), as illustrated by Bradshaw et al. (2012a). The most significant differences in orbital sensitivity of precipitation patterns are found across the equatorial regions, largely driven by shifts in the ITCZ, both in DJF and JJA (Fig. 5c, d). This is generally most pronounced over the ocean, but nonetheless precipitation over central northern Asia, eastern North America, and western Africa is significantly more sensitive to orbital changes at 280 ppm in JJA. In contrast, over central Greenland, western Europe, and central North Africa JJA precipitation is more sensitive at 400 ppm (Fig. 5d). In DJF the most significant changes in precipitation sensitivity over land are found in the Southern Hemisphere, especially in South America and central and southern Africa, in both cases with some regions exhibiting higher sensitivity at 280 ppm but dominantly at 400 ppm (Fig. 5b).

Spatio-temporal phasing of surface air temperatures
While comparison of orbital extremes is probably adequate to investigate the links between climate and orbital forcing, we argue that it may not capture the full variability and leads and lags between the orbital forcing and the climatic response. Our results through a full late Miocene precession cycle show that maximum warming and cooling are not spatially synchronous and strongly vary in time across different regions (Fig. 6). Consequently, the warmest or coldest SATs do not necessarily correspond to precession minima and maxima, respectively.
Our experiments only capture the full variability of a single precession cycle. Obliquity (and eccentricity) values also realistically vary in the ensemble together with precession, but our simulations are not designed to fully capture the variability of an entire eccentricity or obliquity cycle. A detailed separation of the effect of precession and obliquity forcing is beyond the scope of this work, but this could be addressed in a future study. The effect of obliquity on seasonal insolation is especially significant at high latitudes. Nonetheless, an obliquity-driven signal has been found in some low-latitude proxy record for the late Miocene in the Mediterranean region (e.g. Hilgen et al., 1995Hilgen et al., , 2000, despite the small influence of obliquity on low-latitude insolation. The influence of obliquity may explain some of the leads and lags between modelled SAT and precession discussed in this section. For example, there are regions showing largely synchronous warming or cooling, especially in the Northern Hemisphere, but in other areas (even neighbouring ones and in the same hemisphere) maxima and minima can be out of phase with the precessional maximum or minimum by as much as 6 kyr (Fig. 6a). This might be expected in the monsoon regions because of the intensified cloud cover reached at times of minimum precession, but it is less understandable for the other locations. Modelled SATs are more out of phase over the ocean than on land, which may relate to the more direct link between solar forcing and temperature over land than over the ocean. Maximum SATs are consistently not synchronous (4-6 kyr out of phase) with precession minimum/maximum in the eastern North Pacific Ocean, in the region of the Indonesian Throughflow, and in the Southern Ocean (Fig. 6a). Given the location and latitudinal extension across the Southern Ocean, here the lag could be associated with changes in ocean circulation and linked to the pathway of the Antarctic Circumpolar Current. Moderate out-of-phase behaviour (2-3 kyr) is also found in northern and southern Asia, over central North America, part of Greenland, in the Arctic regions, the Indian Ocean, the South Atlantic and over several parts of the Pacific Ocean. In the Southern Hemisphere, the monsoon regions in South America, southern Africa and northern Australia are out of phase by 5 kyr or more.
The obliquity maximum reached in experiment 2 will tend to shift SAT maxima away from the precession minimum towards the obliquity maximum, when the system is sensitive to obliquity and responds directly to summer insolation. In that case, a maximum lead of 5 kyr with respect to precession (minima) can be explained. Note, however, that 65 • N summer insolation varies in anti-phase with precession (see Fig. 1) and is not shifted in the direction of the obliquity extremes. The different response in the two hemispheres, with stronger out-of-phase behaviour in the Southern Hemisphere, might also be partially explained by the use of a modern calendar in these simulations (Chen et al., 2011).
Minimum modelled SATs (Fig. 6b) are mostly not synchronous (4-6 kyr out of phase) with precession minimum/maximum in the North Atlantic Ocean and Nordic Seas, as well as part of the South Atlantic. Strong out-ofphase behaviour is also found over Greenland, northern and central Asia, South America, south of Africa and at several locations in the equatorial regions and in the North Pacific. Moderately out-of-phase (2-3 kyr) temperatures extend over North America, north and central Asia, part of the Arctic and in several locations over the ocean both in the Northern and Southern hemispheres. Because of their location, we suggest that the patterns observed across the North Atlantic and North Pacific, with areas out of phase by up to 4 kyr, are associated with the winter storm tracks. Overall, minimum temperatures exhibit an even more complicated mosaic of patterns than the maximum ones. The response of the climate system at high latitudes is more complex due to vegetation, snow, and sea-ice albedo feedbacks (Tuenter et al., 2005). This could therefore exacerbate leads and lags with the orbital forcing in these regions.
These results further demonstrate the importance of considering orbital variability in order to capture the entire magnitude of the warming/cooling (or wettest/driest periods), especially locally and when considering model-data comparisons. Prescott et al. (2014) also found significant out-ofphase responses when investigating peak warming around two Pliocene interglacials. These authors argued that proxybased reconstruction of temperature time series that rely on cold/warm peak alignment and averaging (e.g. Dowsett and Poore, 1991;Dowsett et al., 2012) could potentially result in significant temporal miscorrelations. This is confirmed by our results from a single late Miocene precession cycle. The bias is relevant for all pre-Quaternary model-data comparison studies, which require a methodology incorporating the effect of orbital variability on climate (Prescott et al., 2014). The more traditional time-average approach must be avoided in order to compare model results with proxy reconstructions robustly (Prescott et al., 2014;Dowsett et al., 2013;Haywood et al., 2013;Salzmann et al., 2013). Bradshaw et al. (2012a) carried out a quantitative terrestrial model-data comparison using a late Miocene data set, which incorporated a conservative estimate of uncertainties associated with both the model output and the data reconstructions. As well as calibration uncertainties in the proxies, model bias and interannual variability, their methodology also considered the potential impacts of poor temporal constraint on determination of the data palaeolocation (see Bradshaw et al., 2012a, for full details of the model-data comparison method- (f) Model-data mismatch is defined as the minimum possible distance to overlap, but here we show that the maximum possible differences could be much greater if the true values for both the model and the data were to lie at the extremes of the uncertainty ranges (Bradshaw et al., 2012a). Note that the relative contributions of model and data uncertainties will vary depending on the variable analysed and for each experiment. The real values are not indicated here as this figure is schematic. ology). The available late Miocene terrestrial proxy record is biased by a sparse and patchy distribution, and low temporal resolution. Despite these large uncertainties, Bradshaw et al. (2012a) found significant discrepancies between the climate model output and the available late Miocene terrestrial proxy record.

Global terrestrial model-data comparison
These authors applied a modern-day orbital configuration to their simulations. Here, as described in Sect. 2, we use the same numerical model and initial setup, but we take into account the full range of variability through the analysed late Miocene precession cycle when undertaking the model-data comparison. This is achieved by selecting the maximum and minimum value through the orbital cycle from the 22 simulations, for every analysed variable in each grid cell. Our definitions and estimates of model-data agreement or mismatch (Fig. 7) and the uncertainties in the model and data are the same as those described by Bradshaw et al. (2012a) but with an extension to the envelope of model uncertainty to include orbital changes.
The methodology developed by Bradshaw et al. (2012a) includes a bias correction which corrects for the offset between the model's simulated preindustrial climate and preindustrial observations. This assumes that even if simulated temperatures and precipitation are not necessarily accurate in an absolute sense, there is a robust relationship between the late Miocene climate and that of the present day (Bradshaw et al., 2012a). For the model, the uncertainty associated with the natural interannual variability within the simulation is also included. For each value this is calculated as 1 standard deviation of the interannual variability of the last 50 years of the model simulation. In addition, given that the observational data sets are characterised by a higher spatial resolution than the model, in the model-data comparison all the model grid cells adjacent to the ones containing the proxy data are considered, where the minimum and maximum values from all of the eight adjacent cells, rather than only the value on the specific grid cell, are used (nine grid cells in total). This is a way to account for the poorly constrained age control on the data, plate rotation uncertainties, and the location of the climate signal recorded by the proxy record (Bradshaw et al., 2012a). Finally, the calibration error for each proxy type is also included and calculated based on modern proxies. Any remaining error must be due to model structural or parametric uncertainties which could only be addressed through multi-model inter-comparison studies.
Overlap or mismatch (Fig. 7) depends on whether the range between the maximum possible model value (M max ) and minimum possible model value (M min ) overlaps with the range between the maximum and minimum data values (D max , D min ). In our case, for each variable and in each grid cell, M max is the maximum value out of 198 (22 × 9, where 9 are the 8 grid cells surrounding the data location plus the grid cell itself, and 22 is the number of orbital simulations) grid cells, plus 1 standard deviation of the interannual variability, and similarly for M min . Finally, the bias correction is applied.
In this way we are able to capture the entire range of variability simulated by the model throughout the full precession cycle for each variable, allowing us to check whether the proxy reconstructions would overlap with model results at any point during the precessional cycle. We can therefore test whether part of the mismatch obtained by Bradshaw et al. (2012a) may be explained by orbital variability. Our model results are compared to mean annual SATs and precipitation, and warm and cold month SATs from the Messinian, reconstructed from proxy data. The data set used is the same compilation of terrestrial proxy reconstructions as Bradshaw et al. (2012a) but with the addition of palaeo-precipitation data for North America by Eronen et al. (2012) Table S1).
The comparison of our results at 280 ppm CO 2 including orbital variability with those of Bradshaw et al. (2012a) demonstrates an overall reduction of the model-data mismatch almost everywhere, both for the mean annual temperature and precipitation records (Fig. 8ia, ib). The only exception is a single data point in the South American continent, showing slight deterioration (see Fig. 7d, depicting the mismatch because of the added orbital variability). We find 766 overlaps (Table 2) from a total of 1193 data points between our orbital ensemble and the Messinian terrestrial proxy data, as opposed to the 610 overlaps found in the Bradshaw et al. (2012a) simulation for the Messinian part of the data set. The cold month temperature and annual precipitation over the Asian, African and North American continents are well matched between our simulations and the data. However, over most of the European continent the proxy record is still warmer and wetter than the climate reproduced by our simulations, both in the warm month and the annual means (Fig. 8iia, iib, iic). The Mediterranean region, where there is the greatest density of observations, gives both the highest match and mismatch between the model and the data. Even when considering the wider envelope of model variability, the simulations still largely fail to capture the magnitude of warming found in the Messinian data, exhibiting mostly colder temperatures (both annual and warm month mean) especially in the Mediterranean region (Fig. 8iia, iic), but with a good match (187 overlaps out of 238 data points) for the cold month mean (Fig. 8iid).
As previously discussed, it is important to remember that in our realistic late Miocene simulations we are only considering one specific precession cycle. We, therefore, cannot capture the full variability of obliquity. In addition, there are other higher-amplitude precession cycles in the Messinian (higher eccentricity values), which means that we are not able to fully capture the maximum amplitude of precessional variations either. Running idealised simulations including the full orbital variability throughout the late Miocene period (using absolute maximum and minimum values for all orbital parameters) may result in an even better match with the proxy record. Bradshaw et al. (2012a) also investigated the impact of using different CO 2 concentrations when modelling the late Miocene climate and obtained a better match with the proxy record using 400 ppm (719 overlaps for the Messinian part of the data set) rather than preindustrial values of 280 ppm (610 overlaps). We have therefore carried out an additional model-data comparison, taking into account the variability between the two precession extreme experiments at 400 ppm www.clim-past.net/11/1271/2015/ for each analysed variable (Fig. 8iii, iv). Our simulations with higher CO 2 concentrations and including orbital variability also exhibit a significantly better match with the Messinian observational record (Fig. 8iv) than the orbital ensemble carried out at 280 ppm, both for mean annual temperature (MAT) and warm month mean temperature (WMT). This is indicated by the presence of 172 overlaps for the MAT (Table  2) and 183 overlaps for the WMT, compared to the 86 MAT and 121 WMT overlaps obtained in the 280 ppm ensemble (Fig. 8iiia, iiic). WMTs in the model at 400 ppm show a good agreement with the proxy data (Fig. 8ivc), which is more improved than the match achieved in the 280 ppm simulations, except for the north-east Asian region. All Messinian WMT reconstructions overlap with the model results in the Mediterranean region (Fig. 8ivc), and there is an almost complete overlap in this region also for the cold month temperature (CMT). Modelled MATs (Fig. 8iva) exhibit both some warmer and colder data points compared to the Messinian observational record, despite generally good agreement over the European continent. There are no major differences in the comparison between CMT and mean annual precipitation (MAP) at 280 and 400 ppm CO 2 concentrations. However, a slight deterioration is found in the CMT (Fig. 8ivb) and in the MAT (Fig. 8ivd). This is indicated by the presence of 185 overlaps for the CMT (Table 2) and 370 overlaps for the MAP at 400 ppm, compared to the 187 CMT and 372 MAP overlaps obtained in the 280 ppm ensemble (Fig. 8iiib, iiid). Bradshaw et al. (2015) also discussed the reasons for modeldata comparison deterioration with higher CO 2 concentrations in certain areas and found that the best fit for mean annual precipitation occurred at 180 ppm CO 2 , despite the best match for SATs resulting at 400 ppm. The reasons for these discrepancies are still not clear, and our results show that these cannot be reconciled by including orbital variability.
As the warmest or coldest temperatures do not necessarily correspond to precession minimum and maximum, the 400 ppm precessional extreme sensitivity experiments do not necessarily capture the full variability of the precession cycle (refer to Fig. 6). At 280 ppm CO 2 , the model-data comparison outputs for the true minimum and maximum resulting from the full ensemble of simulations covering the whole precession cycle are almost identical to the model-data comparison results for just the precession minimum and maximum. In fact, there is a difference of only five overlaps (Ta-ble 2), because the differences in the simulations are smaller than the uncertainties in the proxy reconstructions. However, this may not be the case for regions where well-constrained data are available, such as the Mediterranean Sea.
To summarise, our results imply that accounting for orbital variability, when combined with higher CO 2 concentrations, reduces model-data mismatch by more than 25 % as compared to previous experiments for the late Miocene using a modern orbital configuration (Bradshaw et al., 2012a). In regions where good agreement is obtained between model and data and where in addition high-resolution and more precisely dated proxy records are available for our specific Messinian modelled time period, it would also be possible to estimate during which part of the precessional cycle the proxy reconstruction has been generated (assuming that the model realistically simulates orbital and seasonal variability). For instance, this can be applied to Messinian micropalaeontological data from the Mediterranean Sea that have been sampled on sub-precessional timescales.

African summer monsoon variability between precession extremes
The majority of the late Miocene terrestrial proxy data are concentrated around the margins of the Mediterranean Sea. River discharge into the Mediterranean today is dominated by the Nile. In the late Miocene another North African river which is now dry, the Eosahabi, may also have drained from Lake Chad into the eastern Mediterranean (Griffin, 2002(Griffin, , 2006. Changes in the discharge of these rivers is driven by the summer North African monsoon, which is in turn influenced by precession (e.g. Rossignol-Strick and Planchais, 1989;Lourens et al., 1996;Larrasoaña et al., 2003) and, to a lesser extent, by obliquity (Bosmans et al., 2012(Bosmans et al., , 2015a. We therefore analyse the dynamics of the North African monsoon and its seasonal precipitation and SAT changes throughout our full simulated precession cycle. Here, we consider the North African monsoon system as the combination of both the present-day West African and Central African monsoon dynamics, predominantly controlled by the overriding northsouth large-scale Hadley circulation. Our model results highlight the prominent effect that different orbital configurations have on the African summer monsoon. For instance, the minimum precession simulation exhibits significantly higher SATs over Europe (and gen- Figure 9. (a) SAT and (b) precipitation difference between minimum (pMIN) and maximum (pMAX) precession during the monsoon season (JJAS). The dashed boxes in (a) illustrate how North Africa is split in two areas, northern "box" and southern "box", for analysis (where only the land component is considered these are defined as northern region and southern region). Latitudes and longitudes for the southern "box" are defined according to Thorncroft and Lamb (2005) for the West African monsoon. Absolute values for the monsoon season (JJAS) precipitation at (c) precession minimum and (d) precession maximum.
erally the Northern Hemisphere, as shown in Fig. 3b), but lower values over part of North Africa, ∼ 10-20 • N (Fig. 9a), as a result of increased cloud cover caused by major changes in precipitation patterns over this area (Fig. 9b). The northward shift of the ITCZ is clearly visible in the absolute changes in precipitation between the two precession extremes (Fig. 9c, d). During precession minimum, precipitation > 10 mm day −1 reaches as far north as ∼ 18 • N and intensifies over land (Fig. 9d). By contrast, during precession maximum higher precipitation (positive anomalies) occurs over the Atlantic and only reaches ∼ 10 • N (Fig. 9c).
As well as this land-sea contrast, the northernmost part of the northern African continent exhibits very different patterns from the more southerly area. These two regions (as defined in Fig. 9a and b) are therefore analysed separately.
In the northern region (land-only component of the northern "box" in Fig. 9a), modelled SATs show a very similar seasonal distribution and a single seasonal peak around the month of July in all the simulations presented in Fig. 10a: both extreme precession experiments with 280 ppm CO 2 , in the two control runs (late Miocene and preindustrial palaeogeography with present day orbital forcing), and in the minimum precession simulation with 400 ppm. Considerably higher temperatures are reached during precession minimum (over 35 • C at 280 ppm and close to 40 • C at 400 ppm) while the lowest summer temperatures (< 30 • C) occur in the precession maximum simulation (Fig. 10a). Precipitation ex-hibits a bi-modal distribution, which is most pronounced in the precession minimum simulation, when even in the drier parts of North Africa precipitation reaches 0.8-1 mm day −1 in August (Fig. 10b). Precipitation generally peaks around the months of June and September, but during precession minimum this second and most pronounced peak occurs about 1 month earlier in the season (August) and later (October) in the late Miocene control. The winter months are characterised by extremely dry conditions in all simulations, with precipitation consistently below 0.10 mm day −1 (Fig. 10b).
Modelled SATs in the southern region (land-only component of the southern "box", as defined in Fig. 9a, where latitudes and longitudes are defined as in Thorncroft and Lamb (2005) for the present-day West African monsoon) of North Africa show a weak bi-modal distribution with peaks in April-May and September-October, and this second peak is most pronounced in both precession minimum simulations (Fig. 10c). These summer temperatures of ∼ 28 • C are considerably lower than those in the northern region and are caused by the increased cloud cover during the monsoon season, with little variation in the seasonal distribution between the different simulations (Fig. 10d). However, considerably higher precipitation values are reached in the precession minimum experiments (over 8 mm day −1 ) irrespective of which CO 2 concentrations are used (Fig. 10d). This may reflect a non-linear relationship between changes in African monsoonal precipitation and an increase in CO 2 concentrations, Figure 10. SAT and precipitation seasonal distribution over North Africa (averaged over land in the northern and southern "boxes", as indicated in Fig. 8) for the two precession extremes (pMIN and pMAX) at 280 ppm, precession minimum at 400 ppm and the two control experiments (late Miocene and preindustrial at 280 ppm). Differences due to orbits, palaeogeography and CO 2 concentrations are highlighted by the vertical bars relative to the month of August when the seasonal distribution is not varying. Note that the scales in panel similarly to what has been demonstrated for present-day simulations (Cherchi et al., 2011).
The only difference between the two control experiments (LMctrl and PIctrl) is the palaeogeography, and this results in significantly different precipitation values. For instance, in the late Miocene control experiment, precipitation rates are up to ∼ 2 mm day −1 lower in southern North Africa than they are in the preindustrial control run (Fig. 10d). There are also smaller differences (< 0.3 mm day −1 ) in the northern region (Fig. 10b). Across the whole of North Africa, the preindustrial control experiment is on average warmer than the late Miocene control. This is most pronounced in the northern region, where SATs are up to ∼ 2 • C greater (Fig. 10a), as a result of stronger local summer insolation. Analysis of the different simulations demonstrates that in the northern region the biggest influence on temperature range is orbital variability (∼ 7 • C; Fig. 10a), while CO 2 results in ∼ 3 • C temperature difference and palaeogeography ∼ 1 • C. The striking differences in precipitation in the southern region are again most strongly influenced by orbital variability, which contributes up to ∼ 2.5 mm day −1 to the August peak (Fig. 10d). The June-July-August-September average of SAT and pre-cipitation for each of the experiments summarised in Fig. 10 can also be found in Table 1. This highlights the extended length of the monsoon season during precession minimum at 400 ppm CO 2 , resulting in increased modelled precipitation in the month of September (but no change with respect to the 280 ppm simulation in the month of August) and therefore for the entire period.
Note that, to a lesser extent, obliquity forcing also has an impact on the North African summer monsoon (e.g. Bosmans et al., 2015a, b). Given our experimental design, it is possible to compare experiments 1 and 22, where obliquity is a maximum and minimum, respectively, and precession has very similar values in both simulations (see Fig. 1). The seasonal distribution of precipitation in the North African monsoon region (southern "box") is, however, very similar between the two extreme obliquity simulations and so are the mean annual values, only showing significant differences in the summer months (below ∼ 1 mm day −1 ; Supplement, Fig. S7). This results in much smaller hydrologic changes compared to the extreme precession simulation and these are, therefore, not discussed in more detail in this study.  In the model, the variability in the African summer monsoon between the two precession extremes can largely be explained by changes to the regional circulation -for instance, in the strength of the African westerly jet, which transports moisture into North Africa during precession minima. Because of a greater land-sea temperature differential, low level winds are stronger (> 10 m s −1 ) in the precession minimum simulation (Fig. 11a) and weaker (< 4 m s −1 ) during precession maxima (Fig. 11b), relative to the modern orbit late Miocene control experiment (Fig. 11c). The importance of perturbations to the large-scale atmospheric circulation is also shown by the differences in the strength of the Hadley circulation between these three simulations (Supple-ment, Fig. S11). During precession minimum, the ascending branch is much stronger than in the late Miocene control run and it shows a northward propagation of ∼ 4 • . During precession maximum, the ascending branch is significantly weaker than in the control and located ∼ 3 • further south. This clearly indicates the shifts in the position of the ITCZ during these three simulations.

Impact on vegetation
In our experiments, the substantially increased precipitation at times of precession minimum (Fig. 9b, d) results in a greening of the areas south of the Sahel region (Fig. 12).
During the precession minimum C4 grasses shift to the north (Fig. 12a), colonising areas around 15-20 • N, which are instead covered by the desert fraction (bare soil) during the precession maximum (Fig. 12b). Further south, between ∼ 5 and 15 • N, bare soil is also partially substituted by broadleaf trees in the precession minimum simulation (Fig. 12b, c). A similar amplified precession signal in the monsoon and an extended seasonality within a year when interactive vegetation is included has also been found in both transient (Tuenter et al., 2005) and time-slice simulations (e.g. Doherty et al., 2000;Brostrom et al., 1998). A greening around the Sahel region during this time period is also consistent with geochemical and mineralogical studies (Colin et al., 2014), and a northward displacement of the tree line during precession minima has also been observed in an idealised modelling study (Tuenter et al., 2005). The permanent presence of an extensive desert area in North Africa throughout the entire precession cycle also appears realistic, since both observational (Schuster et al., 2006) and modelling studies (Zhang et al., 2014) suggest that the formation of the Sahara may have been initiated as early as the late Miocene. Vegetation reconstructions for the Late Miocene are also consistent with this hypothesis, indicating the presence of arid conditions starting at around 7 Ma .
We have also investigated the sensitivity of changes in vegetation distribution to varying CO 2 concentrations. However, since the precipitation simulated by precession minima experiments with both 280 and 400 ppm CO 2 are nearly identical over southern North Africa, where the significant vegetation changes are found at 280 ppm (Fig. 12a-c), the small difference in vegetation across this area is unsurprising ( Fig. 12d-f). No major changes are found over North Africa between the two experiments in the expansion of the tree fraction (Fig. 12f), and the differences further south are unrelated to the North African summer monsoon. Patchy differences in C4 grasses distribution increase with CO 2 in the central part of North Africa, where they cover areas that are desert at 280 ppm. C4 grasses decrease to the western side, where they are substituted by the desert fraction (Fig. 12d, e). The less predictable distribution of these changes is also perhaps not unexpected, since CO 2 and vegetation feedbacks do not necessary combine linearly (Bradshaw et al., 2012a). More significant differences in sensitivity may be found in regions outside North Africa, where vegetation productivity is higher (Supplement, Figs. S8, S9 and S10). Exploring this further is, however, beyond the scope of this work.
The recurrence of the so-called African humid periods has been intensively studied both in observational (e.g. Larrasoaña et al., 2003, and references therein) and modelling (e.g. Hely et al., 2009;Liu et al., 2007;Renssen et al., 2006;Joussaume et al., 1999) investigations, especially for the Quaternary period. The proxy record indicates that these periods were characterised by a northward shift in precipitation as a result of a stronger African summer monsoon, paced by astronomically forced insolation changes. To date, mod-elling studies have largely failed to simulate the northward penetration of the African summer monsoon beyond 21 • N and increase precipitation sufficiently to simulate the mid-Holocene "Green Sahara" (Brovkin et al., 1998;Claussen et al., 1999) conditions (e.g. Harrison et al., 2015;Bosmans et al., 2012;Braconnot et al., 2007;de Noblet-Ducoudre et al., 2000). These conditions would allow savanna-like vegetation to expand northward, beyond the central Saharan watershed (Larrasoaña et al., 2003). Bosmans et al. (2012) hypothesised that the lack of interactive vegetation could be the main reason for the insufficient precipitation over the Sahara in mid-Holocene simulations. However, in our simulations which are coupled with a vegetation model, the summer precipitation increase during precession minimum is still confined south of 21 • N in North Africa. Assuming that the monsoon system in the late Miocene was similar to that of the Quaternary, this indicates that even our fully coupled model still fails to represent relevant processes driving precipitation in the Sahel regions. This is perhaps suggesting the lack of relevant teleconnections in the model, such as those found with North Atlantic dynamics (e.g. Barandiaran and Wang, 2014;Zhang and Delworth, 2006). In addition, this could also be due to the low sensitivity of the land-atmosphere coupling which characterises models of the HadAM3 family (Koster et al., 2006).

Seasonality of the African summer monsoon on sub-precessional timescales
Our experimental design allows us to analyse the seasonal distribution of SATs and precipitation patterns over North Africa not only for the two precessional extremes but also throughout the different stages of the orbital cycle (Fig. 13). The highest SATs (up to 35 • C) are reached in the northern region during the summer months (Fig. 13a). In the southern region, SATs remain below 30 • C throughout the entire cycle (Fig. 13b). The highest quantity of precipitation (up to 2500 mm day −1 ) is found in the southern region during the summer months and especially around the precession minimum (Fig. 13d). In the northern region, which is outside the area influenced by the summer monsoon, drier conditions persist throughout the entire cycle and in all seasons (Fig. 13c), with values consistently below the driest periods experienced in the southern region (maximum 250 mm day −1 ). The mean annual values show the correlation with the precession forcing, which is positive for modelled precipitation and SAT in the southern region, and negative for SAT in the northern region. However, some lags between orbital forcing and the climate response can also be seen. For instance, maximum precipitation in the southern region occurs at the same time as the precession minimum, but minimum precipitation lags the precession maximum by about 2 kyr (Fig. 13d). In the northern region, where precipitation rates are an order of magnitude smaller than in the southern region, the phas- Also note that the four annual mean panels are not on the same scale, as their aim is to show the phasing with orbital forcing rather than to compare the actual values.
ing with precession is less clear; maximum annual precipitation corresponds to the precession minimum, but the signal flattens out in the remaining part of the cycle around the precession maximum, and minimum values are reached around simulations 13 and 20 (Fig. 13c). Note that part of the changes in the northern region, characterised by low precipitation value, may be due to interdecadal/centennial variability in the model.
Minimum annual SAT in the southern region occurs with a 1 kyr lag after the precession minimum, while maximum SAT lags the precession maximum by 2-3 kyr (Fig. 13b). In this area, the SAT response to orbital forcing is linked to the increased cloud cover at times of precession minimum (maximum monsoon strength), as discussed in Sect. 3.4. In the northern region, maximum annual SAT leads the precession minimum by ∼ 5 kyr, while minimum SAT leads the precession maximum by ∼ 4 kyr.
Although maximum values for modelled SATs in the northern region and for precipitation in the southern region can be correlated with the precession minimum, the seasonal response is not "symmetrical", but rather it exhibits an elongated and slightly tilted structure (Fig. 13a, d). This asymmetrical response around the precession minimum has also been observed in transient idealised orbital simulations with a model of intermediate complexity (Tuenter et al., 2005) and has been explained by the extended length of the monsoon season around the precession minimum. At this stage in the orbital cycle, the North African monsoon can start up to 1 month in advance and end a month later than average monsoon timing (Fig. 13d), in agreement with the results of Tuenter et al. (2005). One possible explanation for this phe-nomenon is the presence of a larger vegetated area during precession minimum, which modifies the albedo feedback and results in a longer monsoon season (Tuenter et al., 2005).
Finally, modelled SATs in the northern region consistently peak during the month of July (Fig. 13a), with the exception of simulations 11 to 15 (August) and 20 to 22 (June). In the southern region, the month of August consistently exhibits the highest values of precipitation (Fig. 13d), apart from simulations 2, 3, and 22, where July is the wettest month, and simulation 21, in which June has the highest precipitation rates. This differs from the results of Tuenter et al. (2005), whose idealised simulations consistently showed maximum precipitation rates during the month of July, in all the precession minimum experiments throughout their entire simulated time period. This difference is likely due to the fact that our simulations use realistically varying orbital parameters throughout one precession cycle and that interactive vegetation is also included.

Synthesis and conclusions
Using a fully coupled ocean-atmosphere-vegetation model (HadCM3L) we have investigated the impact of orbital forcing on global climate, with a focus on the North African monsoon, throughout a full late Miocene precession-driven insolation cycle with varying obliquity. In the model, the evolution of global mean annual SATs through our simulated orbital cycle is not directly influenced by changes in insolation. Annually, the Southern Hemisphere generally exhibits higher SATs than the Northern Hemisphere (Fig. 2), in contrast with the present-day configuration. Seasonal SAT changes show a more complex response throughout one full late Miocene precession cycle, exhibiting leads and lags with insolation. This response is in part caused by the use of interactive vegetation in our simulations, in agreement with idealised transient simulations carried out with an earth system model of intermediate complexity (Tuenter et al., 2005). In addition, the influence of obliquity in our simulations may also play a role.
The difference in mean annual modelled SAT between the two precession extremes is generally small on a global scale (Fig. 3e). There are, however, significant local changes resulting from the dominance of strong seasonal signals over the annual mean. Examples are the Nordic Seas, which are over 1.5 • C warmer during precession minimum, and the African and Indian monsoon regions, which are more than 1.5 • C colder. Seasonal differences in SAT are much larger. In line with previous studies (e.g. Lunt et al., 2013;Yin and Berger, 2012), widespread global cooling can be observed during DJF (Fig. 3a), especially on land, with the exception of warming (up to 4.5 • C warmer during precession minimum) in the Arctic and the Nordic Seas. Widespread warming is instead found in JJA (Fig. 3b), especially in the Northern Hemisphere, with lo-cal cooling (up to 3.5 • C colder) over the African and Indian monsoon regions, as a result of increased cloud cover. In the model, precipitation differences between the two orbital extremes (Fig. 4) are largest during JJA (> 3 mm day −1 ) and are mostly centred on the equatorial African and Asian monsoon regions. DJF anomalies are smaller in the Northern Hemisphere and reflect changes in the storm track systems, especially across the North Atlantic. The sensitivity of SATs and precipitation to CO 2 in combination with orbital forcing is moderate on a global scale (Fig. 5), with locally enhanced SAT sensitivity to orbital forcing at both 400 and 280 ppm.
Globally, the warmest and coldest modelled SATs are not necessarily reached during precession minimum and maximum, because maximum warming and cooling are not spatially synchronous and vary in time across different regions (Fig. 6). This is especially relevant in the Southern Hemisphere, where maximum SAT can be out of phase with the precession minimum by as much as 6 kyr. This has implications for the correlation of proxy-based temperature reconstructions with warm/cold peaks (e.g. Dowsett and Poore, 1991) and could result in significant temporal miscorrelations, as discussed by Prescott et al. (2014) for the Pliocene warm period. Bradshaw et al. (2012a) found that the available proxy data for the late Miocene exhibit significantly warmer mean annual temperatures than at the present day for all locations where the palaeodata are different from their modern-day climate estimates. The same was observed for mean annual precipitation, with the palaeorecord appearing wetter at all locations where the climate of the late Miocene climate is different from a potential natural modern one (see Figs. 7 and 11 in Bradshaw et al., 2012a). Our simulations through one full precessional cycle demonstrate that some of the modeldata mismatch for the late Miocene found by Bradshaw et al. (2012a) can be explained by orbital variability (Fig. 8). However, our simulations at 280 ppm still largely fail to reproduce the same magnitude of warming indicated by the Messinian proxy reconstructions. A better match is achieved including the variability between the two precession extremes and 400 ppm CO 2 . This demonstrates that, in addition to using an appropriate palaeogeography and higher CO 2 concentrations (Bradshaw et al., 2012a), accounting for orbital variability can reduce the model-data mismatch for the late Miocene. However, some disagreement between the model output and the data is still present across some areas.
In our simulations, the North African monsoon is highly sensitive to orbital forcing (Fig. 9), which strengthens the African westerly jet during precession minimum (Fig. 11), intensifying precipitation over North Africa significantly and leading to a greening of the region south of the Sahel (Fig. 12). The African monsoon is also sensitive to palaeogeographic changes, but largely insensitive to varying CO 2 concentrations between 280 and 400 ppm (Fig. 10d). Nonlinear behaviour with respect to CO 2 forcing for the late Miocene is consistent with modern-day climate simulations of the North African monsoon (Cherchi et al., 2011). Our ensemble of simulations demonstrates that both SATs and precipitation over the North African monsoon regions exhibit significant differences in their seasonal distribution through a full precession cycle. SAT is significantly influenced by the amount of cloud cover during the monsoon season, while precipitation is enhanced between June and September during precession minimum (Fig. 13). The evolution of these two variables is, however, not symmetrical around precession minimum and maximum because of the extended length of the monsoon season as a result of vegetation feedbacks (Tuenter et al., 2005).
In conclusion, we suggest that future studies comparing model and proxy data will need to take into account not only differences in palaeogeography and CO 2 concentrations but also orbital variability. This is not only relevant for the late Miocene but more generally for all pre-Quaternary modeldata comparison studies, where the proxy reconstructions largely rely on time-averaged palaeoenvironmental syntheses (Prescott et al., 2014).
The Supplement related to this article is available online at doi:10.5194/cp-11-1271-2015-supplement.