Modelling the enigmatic Late Pliocene Glacial Event — Marine Isotope Stage M2

for the Pliocene signi ﬁ variability and intervals apparently a Marine Isotope Ma) a recognisable cooling event that disturbs an otherwise relatively (compared to present-day) warm background climate state. It remains unclear whether this event corresponds to signi ﬁ cant ice sheet build-up in the Northern and Southern Hemisphere. Estimates of sea level for this interval vary, and range from modern values to estimates of 65 test the hypothesis larger-than-modern ice sheet con ﬁ gurations have existed Climate model results are compared with proxy climate data available for M2 to assess the plausibility of each ice sheet con ﬁ guration. Whilst the outcomes of our data/model comparisons are not in all cases straight forward to interpret, there is little indication that results from model simulations in which signi ﬁ cant ice masses have been prescribed in the Northern Hemisphere are incompatible with proxy data from the North Atlantic, Northeast Arctic Russia, North Africa and the Southern Ocean. Therefore, our model results do notprecludethepossibilityoftheexistenceoflargericemassesduringM2intheNorthernorSouthernHemisphere. Speci ﬁ callytheyarenotabletodiscountthepossibilityofsigni ﬁ canticemassesintheNorthernHemisphereduring the M2 event, consistent with a global sea-level fall of between 40 m and 60 m. This study highlights the general need for more focused and coordinated data generation in the future to improve the coverage and consistency in proxy records for M2, which will allow these and future M2 sensitivity tests to be interrogated further.


Introduction
The Late Pliocene (3.6-2.58 Ma; Piacenzian Stage; Fig. 1) represents the most recent period in Earth history where average global temperatures were warmer than present-day (Haywood and Valdes, 2004;Dowsett et al., 2010) and peak carbon dioxide (CO 2 ) levels were be-tween~50 and 125 ppmv higher than pre-industrial (e.g. Bartoli et al., 2011). The interval between 3.264 and 3.025 Ma (the mid-Pliocene Warm Period; Dowsett et al., 2010) has been widely studied, in terms of both modelling and geological data collection, in order to understand the causes of Pliocene warmth , and to establish the long-term effect (Earth System Sensitivity) of increased levels of CO 2 in the atmosphere (Lunt et al., 2010).
However, although not as dramatic as the glacial and interglacial cycles that typified the Pleistocene, the Late Pliocene also exhibited climate variability and periods which were apparently cooler than modern (Lisiecki and Raymo, 2005; Fig. 1). Of particular interest is the major cooling event which occurred at 3.3 Ma in Marine Isotope Stage (MIS) M2 (during the Mammoth reversed polarity subchron; Fig. 1). This "Pliocene glacial" represents the largest positive isotope excursion during the Pliocene prior to the Plio-Pleistocene boundary (Lisiecki and Raymo, 2005; Fig. 1). A sharp peak in oxygen isotopes occurs over a 20,000 year interval within the longer 3.312-3.264 Ma MIS M2 (De Schepper et al., 2013). It has been referred to as a failed attempt of the climate to reach a full glacial state, not reached until the latest Pliocene (~2.75 Ma), due to unfavourable conditions in some aspect of the climate system (Haug and Tiedemann, 1998).
The period of sea level fall culminating in the M2 event begins at around 3.4 Ma and is part of a longer term cooling trend (~100 kyr). However, the M2 isotope excursion steepens deeply between 3.305 and 3.285 Ma (De Schepper et al., 2013;Fig. 1) and the magnitude of this event in isolation would potentially require ice sheets to advance and retreat over a geologically and glaciologically short time period (perhaps as short as 20,000 years). The total Antarctic contribution to post-LGM sea level rise was probably b 10 m of equivalent eustatic sea level (Bentley et al., 2014). Using this as a basis implies that the M2 event could be potentially entirely explainable by changes in Antarctica or additionally require a substantial build-up of ice in the NH (dependent on which sea level reconstruction is considered). Local evidence for glaciation prior to the large scale initiation of NHG can be found across the high latitude Northern Hemisphere (McDougall and Wensink, 1966;Geirsdóttir and Eiríksson, 1994;Krissek, 1995;Lagoe and Zellers, 1996;Kleiven et al., 2002;St. John and Krissek, 2002;Gao et al., 2012;Menzies et al., 2013;Knies et al., 2014;see also De Schepper et al. (2014)), but there is no clear indication of large ice sheets sufficient to explain M2 sea level falls. The isolation of the M2 glacial event represents an intriguing and rapid climate shift in an otherwise comparatively warm Late Pliocene (Fig. 1). This study presents the first global climate model simulations specifically targeting the M2 glacial event and uses different potentially analogous Plio-Pliestocene ice sheet configurations that encompass a range of plausible M2 ice sheet volumes.

Model description -HadCM3
All of the General Circulation Model (GCM) simulations described in this paper were carried out using the UK Met Office coupled Fig. 1. Marine isotope stage M2 in relation to the long-term climate evolution of the Late Pliocene as shown in the benthic oxygen isotope record of Lisiecki and Raymo (2005). The PRISM warm interval (3.264 to 3.025 Ma) is shown by the shaded grey bar. Modern day benthic values are represented by the horizontal dashed line in the top panel. Obliquity (°), eccentricity, precession and the July insolation at 65°N (W m −2 ) is also shown for reference (Laskar et al., 2004).
atmosphere-ocean GCM HadCM3 (Version 4.5; Gordon et al., 2000). The horizontal resolution of the atmosphere is 2.5°in latitude by 3.75°i n longitude and consists of 19 layers in the vertical (Pope et al., 2000). The atmospheric model has a time step of 30 min and includes a radiation scheme that can represent the effects of major and minor trace gases (Edwards and Slingo, 1996). A parameterisation of simple background aerosol climatology is also included (Cusack et al., 1998). The convection scheme is that of Gregory and Mitchell (1997). This version of HadCM3 uses the MOSES2.1 land-surface scheme, which includes the representation of the freezing and melting of soil moisture. The representation of evaporation includes the dependence of stomatal resistance on temperature, vapour pressure and CO 2 concentration (Cox et al., 1999). In all of the simulations presented here the dynamic vegetation model TRIFFID (Cox, 2001) was also employed.
The spatial resolution over the ocean in HadCM3 is 1.25°by 1.25°a nd the model has 20 vertical layers. The horizontal resolution allows the use of a smaller coefficient of horizontal momentum viscosity leading to an improved simulation of ocean velocities compared to earlier versions of the model (Gordon et al., 2000). The ocean model includes the use of the Gent-McWilliams mixing scheme (Gent and McWilliams, 1990). The sea-ice model uses a simple thermodynamic scheme and contains parameterisations of ice drift and leads (Cattle and Crossley, 1995). Gordon et al. (2000) illustrated that HadCM3 is capable of reproducing many aspects of the observed heat budget of the Earth and the model has been used for future climate predictions in the IPCC AR4 (Solomon et al., 2007). HadCM3 has also been used in the Palaeoclimate Modelling Intercomparison Project to simulate Last-Glacial Maximum and Mid-Holocene climates (Braconnot et al., 2007) and also the mid-Pliocene Warm Period (e.g. Haywood and Valdes, 2004;Dolan et al., 2011;Haywood et al., 2013).

Experimental design and boundary conditions
A suite of ten simulations were carried out with HadCM3 (Table 1). All simulations have been run for at least 500 years in order for the climate to reach a satisfactory state of equilibrium. The surface boundary conditions, which are the prescribed two-dimensional fields that remain fixed for the duration of these HadCM3 simulations, are detailed in Table 1. Vegetation, SSTs and sea-ice distribution are dynamically predicted by HadCM3 during the simulation.
We have carried out two control simulations that are equivalent to the HadCM3 pre-industrial and Pliocene (Bragg et al., 2012) standards used within the Pliocene Model Intercomparison Project (PlioMIP; Haywood et al., 2011). The Pliocene control (Pliocene Ctrl ) uses boundary conditions (e.g. topography and ice sheet extent) derived from the Pliocene Research Interpretation and Synoptic Mapping dataset (PRISM3: Dowsett et al., 2010 and Fig. 1). The key difference between this Pliocene Ctrl and the standard HadCM3 Pliocene simulation submitted to PlioMIP (e.g. Bragg et al., 2012) is the use of a dynamic vegetation model.
The remaining eight simulations include changes to (i) CO 2 , (ii) orbit, (iii) ice and (iv) topography. Reconstructions of Pliocene CO 2 have tended to focus on estimating high end-members rather than potential glacial values. However, some records do provide suggestions as to the atmospheric CO 2 concentration at MIS M2. Bartoli et al. (2011) show a lower limit of~220 ppmv at MIS M2. DeConto et al. (2008) additionally demonstrate that initiation of Northern Hemisphere Glaciation (NHG) is unlikely at CO 2 levels above 280 ppmv and therefore we have chosen to run all of our glacial simulations with both 280 ppmv and 220 ppmv CO 2 . Orbital forcing appropriate to 3.3 Ma has also been applied in all M2 glacial sensitivity experiments. The astronomical solution of Laskar et al. (2004) was used to derive the parameters of eccentricity, obliquity and precession required by HadCM3 for 3.3 Ma.
Deciding upon the climate model boundary conditions for ice extent and topography is less straightforward than CO 2 or orbit. In contrast to the mid-Pliocene Warm Period (3.264-3.025 Ma), where detailed reconstructions are available , there are no explicit boundary conditions designed for use when simulating a Pliocene glacial climate. Ice sheet extent and elevation, and associated changes in sea level are poorly constrained at MIS M2. However, to perform 'snap-shot' GCM modelling, some information regarding the configuration of ice masses needs to be given the model.
We chose to implement Quaternary ice sheet reconstructions based on ice sheets presented in Singarayer and Valdes (2010) as an approximation to different MIS M2 ice sheet scenarios. The ice sheet reconstructions were originally based on the work of Peltier and Fairbanks (2006), which used the Martinson et al. (1987) SPECMAP record of δ 18 O history to constrain the evolution of land ice volume from the Last Interglacial (LIG) up to the Last Glacial Maximum (LGM). The specific time points used for the Pliocene glacial scenarios were 76 ka (Large-M2 M2Orbit ) and 116 ka (Medium-M2 M2Orbit ) and represent large and mediumsized M2 ice sheets respectively (Fig. 2). The total ice volume prescribed in these boundary conditions gives an approximate sea level change of 60 m (~8 m from Antarctica) for the large scenario and 40 m (~6 m from Antarctica) for the medium scenario. Changes in the global ice volume also necessitate alterations to the land-sea mask configuration (see Fig. 2). The main differences in land-sea configuration due to a drop in sea level are seen in shallow shelf areas such as the Barents Sea and around the coast of North America.
We also have a small MIS M2 scenario, which uses the pre-industrial ice sheet distribution (Small-M2 M2Orbit ; Fig. 2a) and finally a very small ice sheet scenario that utilises the PRISM3 ice sheet reconstruction (PRISM-M2 M2Orbit ; Fig. 2c; Table 1), which is compatible with a sea level high-stand of less than 25 m . Topography outside the ice sheet regions for this simulation has been kept as present day (Fig. 2) for consistency between the experiments. PRISM-M2 M2Orbit therefore, has a different topography than the PRISM3 reconstruction Table 1 Boundary condition details used in Pliocene Glacial experiments in terms of ice sheet configuration, CO 2 concentration, prescribed orbit and topography outside of the ice sheet regions. Where more than one value is displayed (e.g. for CO 2 and topography) all of these scenarios have been undertaken in separate experiments. There are ten simulations in total.  Fig. 2d) and the HadCM3 PlioMIP experiment submission (Bragg et al., 2012). Given the limited constraints from geological data as to the nature of the ice sheets during MIS M2, the GCM boundary conditions used here should be viewed as potential scenarios rather than complete reconstructions. We have not explicitly considered how the altering the distribution of ice between the hemipsheres in any of our scenarios might affect the results. Our strategy is to simply explore, as far as current geological constraints allow, what may be a more or less likely ice sheet volume and distribution scenario for the M2 event using plausible global ice sheet configurations. Specifically, we assess whether our specified ice sheet reconstructions lead to a simulation of climate that is consistent or inconsistent with the available palaeoenvironmental data from MIS M2.  Table 2). Prescribing larger ice sheets in HadCM3, also strengthens the equator-to-pole temperature gradient (Table 2), with the high northern latitudes becoming up to 3°C cooler (or up to 5.6°C cooler if CO 2 is lowered to 220 ppmv; Fig. 3g and Table 2). This cooling principally reflects the larger (and topographically higher) land ice coverage, causing a substantial change in surface albedo and decreasing the amount of solar radiation absorbed at the Earth's surface. Some cooling can also be seen downstream of the ice sheets and this is related to the increased strength of NH westerly flow. In regions where the land-sea mask has been altered to account for exposed shelf areas as sea level falls (e.g. the Barents Sea; Fig. 2) there is also a significant local temperature decrease as the land becomes colder than the replaced sea surface ( Fig. 3e and g). In all but the large (Large-M2 M2Orbit 280 ) ice sheet scenario, HadCM3 simulates a MAT warming (of up to 3°C) over areas of South America, southern Africa and North Australia (Fig. 3a, c and e), which might be associated with reduced precipitation rates in these regions (Fig. 3b, d and f). The warming and reduction in precipitation over the Amazon (especially in Medium-M2 M2Orbit 280 ) is also linked to enhanced trade wind flow and changes in vegetation in this region.

Mean annual surface temperature
In spite of the M2 orbit and no increase in CO 2 , the simulation using the PRISM3 distribution of ice sheets (PRISM-M2 M2Orbit 280 ) remains warmer than pre-industrial (0.64°C) but with relatively small overall polar amplification (Table 2). However, the warming seen in PRISM-M2 M2Orbit 280 is primarily due to the warmer than pre-industrial temperatures simulated over ice sheets regions where ice has been removed in the PRISM3 reconstruction ( Fig. 3a; Dowsett et al., 2010). M2 ice sheet sensitivity experiments run using 220 ppmv CO 2 were between 1.4°C and 1.6°C cooler than those using 280 ppmv, however similar regional patterns of temperature change were simulated.

Seasonal surface temperature
Fig. 4 displays the JJA (June, July, August) and DJF (December, January, February) temperature responses for the ice sheet sensitivity experiments. In general seasonal changes in surface air temperature track mean annual changes over large parts of the Earth in all four M2 ice sheet scenarios. An exception to this occurs primarily in coastal regions of Northern Europe and the Barents Sea. Due to the different heat capacities of land and ocean, areas which have become land due to sea level falls given the proposed ice sheet scenarios (Fig. 2) exhibit a greater range in seasonal temperatures than in the Pre-Ind Ctrl . In the Large-M2 M2Orbit 280 simulation for example, where the ice sheet loading on Antarctica, North America and Eurasia is equivalent to a sea level change of~60 m, coastal regions of Alaska, the Baltic and North Sea and the Barents and Kara Sea became subaerial (Fig. 2g). In the boreal summer months, simulated temperatures in the Arctic coastal regions were over 10°C higher than Pre-Ind Ctrl , with temperature increases of 2-4°C also propagating further into the continent (Fig. 4).

Mean annual precipitation
For the simulations in which globally ice volume is greater than present day, the mean annual average precipitation is reduced by between 0.03 and 0.15 mm day −1 (Fig. 3 and Table 2). The simulation, Small-M2 M2Orbit 280 , which essentially displays the effect on the pre-industrial climate of the orbital conditions at 3.3 Ma, shows minimal changes in MAP (ΔP = 0.005 mm day −1 ; Fig. 3d). In the tropical regions shifts in the positioning of the Intertropical Convergence Zone (ITCZ) and also a strengthening/weakening of the South Asian summer Monsoon cause considerable difference between the pre-industrial control and the M2 sensitivity experiments. Broadly, precipitation changes in the Tropical regions are consistent throughout all M2 experiments, suggesting this is a robust feature of the M2 climate (according to HadCM3).
Generally, changes in mean annual precipitation outside of the tropical regions reflect alterations in the ice sheet regions in all M2 ice sensitivity simulations (Fig. 3). Areas where ice is removed in the climate model become significantly wetter in the PRISM-M2 M2Orbit 280 scenario (e.g. margins of East Antarctica, West Antarctica and Greenland; Fig. 3b). The ice sheet topography for ice scenarios Medium-M2 M2Orbit 280 and Large-M2 M2Orbit 280 resulted in strong localised cooling due to albedo and elevation feedbacks, which in turn reduced MAP over the ice sheet regions (~1 mm d − 1 ; Fig. 3e-h). In Large-M2 M2Orbit 280 pronounced decreases in precipitation are also evident downstream of the implemented ice sheets, with the atmosphere becoming drier especially in the northern hemisphere (NH; Fig. 3h).

Table 2
Summary of climate diagnostics for all simulations including global mean annual temperature (°C), boreal summer (JJA) and winter (DJF) global mean temperature, the change in the equator-to-pole temperature gradient from Pre-Ind Ctrl (°C), global mean annual precipitation (MAP; mm d −1 ), the strength of the Atlantic Meridional Overturning Circulation (AMOC; Sv) and the hemispheric maximum and minimum sea ice areal extent (and changes from Pre-Ind Ctrl as a percentage).

Ocean surface temperature
Over the oceans sea surface temperatures (SSTs; Fig. 3) generally become globally cooler relative to the pre-industrial as the planetary ice volume increases ( Fig. 3e and g). Given the prescription of PRISM3 ice sheets in our M2 glacial experiments, SSTs still increase in large regions of the world and especially in areas where there are changes in the sea-ice extent ( Fig. 3a; Table 2). This pattern of change is similar to that seen in other standard Pliocene experiments (e.g. Haywood et al., 2013), although the magnitude of warming is reduced here due  to the reduction in CO 2 and the prescription of an M2 orbit. The greatest cooling in the SSTs occurs in Large-M2 M2Orbit with decreases in temperature of over 3°C at high latitudes (Fig. 3g). The model predicts the propagation of cooler temperatures into the tropics and also the Southern Ocean in this experiment and there are few regions where any warming exceeds 0.5°C.

Sea ice
The implementation of larger-than-modern ice sheets in the NH cause the maximum winter and minimum summer sea ice extents to increase southward into the Bering Sea and the North Pacific as well as reaching considerably further down the coast of Canada towards the Gulf of St. Lawrence (Fig. 5 and Fig. S1). The overall change in NH maximum sea ice extent from Pre-Ind Ctrl is + 1.03 × 10 6 km 2 in the Large-M2 M2Orbit 280 simulation, − 2.09 × 10 6 km 2 in the Medium-M2 M2Orbit 280 simulation, + 0.47 × 10 6 km 2 in the Small-M2 M2Orbit 280 simulation and −0.22 × 10 6 km 2 in the PRISM-M2 M2Orbit 280 simulation (Table 2; Fig. 5). For the simulations with larger-than-modern ice sheets, the southward sea ice advance is partially offset by reduced ocean extent in the Arctic associated with changes in the land-sea mask ( Fig. 2; Table 2). SH maximum sea ice extent increases in the larger-than-modern ice sheet scenarios, with the greatest increase shown in the Large-M2 M2Orbit 220 simulation (+ 7.84 × 10 6 km 2 (+ 33.7%) relative to Pre-Ind Ctrl ; Table 2; Fig. 5g). In the Large and Medium M2 simulations, the SH increase in both maximum ( Fig. 5f and h) and minimum (Fig. S1) sea ice extent is also associated with intensification and a northward shift of the South Polar westerlies relative to the Pliocene control (Fig. S2).

Mixed layer depth and salinity
The depth of the mixed layer as an annual mean is strongly affected in the Large-M2 M2Orbit 280 and Medium-M2 M2Orbit 280 ice sheet scenarios. Associated with the southward deflection of the westerlies over the Atlantic Ocean due to the larger-than-modern ice masses, there is deepening (over 15 m) of the mixed layer depth in the southern North Atlantic with a southward shift of the North Atlantic surface gyre ( Fig. S3c and S3d). In comparison to the Pre-Ind Ctrl , the effect of prescribing an M2 orbit (e.g. Small-M2 M2Orbit 280 ) is negligible (Fig. S3b), however when PRISM3 ice sheets are prescribed (PRISM-M2 M2Orbit 280 ) the mixed layer depth increases (by 5 to 15 m) in the North Atlantic in areas south of Iceland (Fig. S3a). This deepening of the mixed layer is related to surface temperature increases and enhanced ocean salinity values in this region, in addition to an overall reduction in sea ice relative to the Pre-Ind Ctrl (Fig. S4a) and an increase in evaporation.

Ocean circulation
Changes in ocean circulation are also predicted by HadCM3 for the glacial sensitivity scenarios (e.g. Fig. 6), although it should be noted that the response, for example of the Atlantic Meridional Overturning Circulation (AMOC), exhibited by climate models is not consistent between models both for the LGM (Otto-Bliesner et al., 2006;Weber et al., 2007) and the Pliocene . Nevertheless, in the larger-than-modern ice sheet scenarios presented here, the AMOC response is similar to that suggested by data and most models for the LGM . There is a reduction in overturning strength in the main cell and a slight southward shift of the downward branch of the AMOC relative to the Pre-Ind Ctrl (Fig. 6e and f). Changes in the top 500 m of the Atlantic ocean in Large-M2 M2Orbit 280 and Medium-M2 M2Orbit 280 reflect wind-driven changes in the mixed-layer depth ( Fig. 6; Fig. S3c and S3d).

Regional data-model comparison
In order to determine the plausibility of the glacial scenarios presented, it is important to compare the resulting climatologies to available palaeoenvironmental estimates. We have chosen four regions where high resolution data exists which gives some information about wider climatic conditions at 3.3 Ma; (i) the North Atlantic, (ii) Africa, (iii) Northeast Russia and (iv) the Southern Ocean.

North Atlantic
Using a combination of oxygen isotopes, Mg/Ca and alkenone palaeothermometry, De Schepper et al. (2009Schepper et al. ( , 2013 characterise the North Atlantic Ocean during MIS M2. A conceptual model explaining the reduction in the North Atlantic Current (NAC) and the initiation of NH ice sheets at this time is put forward, invoking changes in the tectonic configuration of the Central American Seaway (CAS) and thus an increase in Pacific-to-Atlantic flow (De Schepper et al., 2013). Whilst we have not explicitly tested mechanisms for NH ice sheet initiation, the modelled scenarios can be assessed in terms of how well they display a reduction in the strength of the NAC. Ocean current velocities for DSDP Site 610, IODP Site U1308 and IODP Site U1313 have been calculated for the glacial sensitivity scenarios. At U1313 for both Large-M2 M2Orbit 280 and Medium-M2 M2Orbit 280 there is a reduction in the surface velocity of the ocean (−0.6 and −1.2 m s −1 respectively), which is consistent with palynological data and alkenone-based SSTs suggesting a southward shift and/or slowdown of the NAC at this site (Naafs et al., 2010;De Schepper et al., 2013). Large-M2 M2Orbit 280 also displays a weakening of the surface currents around M2 at U1308 (−0.47 m s −1 ). However  (Dowsett et al., unpublished;Lawrence et al., 2009;Naafs et al., 2010). Fig. 7 shows the SST anomalies from the modern value (derived from the World Ocean Atlas Database (Locarnini et al., 2005)) for each data site overlain on the HadCM3 predicted SST anomalies from Pre-Ind Ctrl . Benthic oxygen isotope records at M2 (Fig. 1) suggest a cooling in bottom water temperatures and/or an increase in ice volume at 3.3 Ma. However, data depicting the surface expression of M2 show no clear pattern of cooling when compared to modern temperatures at the four sites, although all M2 temperatures are cooler than the background Pliocene state (De Schepper et al., 2009;Lawrence et al., 2009;Naafs et al., 2010;De Schepper et al., 2013). The proxy-based SSTs considered here range from approximately modern (to within ±0.5°C) to warmer than present-day (up to 3.7°C at U1308, Dowsett et al., unpublished). Where multi-proxy data exists (Sites U1313 and U1308) there are notable differences in the SST anomalies predicted by alkenones, Mg/Ca and faunal-based methods (over 4°C at U1313; Fig. 7). When comparing the model results with proxy data, the warming exhibited at Sites 610 and 982 relative to modern is broadly simulated in PRISM-M2 M2Orbit 280 (Fig. 7a), but is also evident to a lesser degree in Medium-M2 M2Orbit 280 and Large-M2 M2Orbit 280 ( Fig. 7c and d). At Site U1308 the proxy data is difficult to interpret, however the larger-than-modern ice sheet scenarios would seem to overestimate cooling in this region as shown by the data. The performance of the model at Site 1313 is again difficult to determine given the conflicting nature of the proxy-record (Fig. 7). The coldest proxy record at U1313 (the faunal-based SSTs; Dowsett et al. unpublished) suggests that the cooling exhibited by HadCM3 in the Medium-M2 M2Orbit 280 and Large-M2 M2Orbit 280 is most consistent, however the temperature change shown in PRISM-M2 M2Orbit 280 agrees better with the Mg/Ca records (De Schepper et al., 2013).
It is very difficult to draw robust conclusions from this datamodel comparison (DMC). Moreover, it has been shown that the response of Pliocene climate models is inconsistent in the North Atlantic in particular  and therefore the results showed here are likely to be model dependent. Over time an increase in the understanding of differences between proxy results and further data collection may shed light on the plausibility of the M2 ice configurations considered here.

Africa
There are a number of high-resolution climate archives that cover MIS M2 and could potentially allow for DMC over the African continent (e.g. Leroy and Dupont, 1994;Bonnefille et al., 2004;deMenocal, 2004). Whilst some records display a significant change in North African climate (e.g. Leroy and Dupont, 1994;Bonnefille et al., 2004), other records suggest little change at this time (e.g. deMenocal, 2004). Evidence for climate variability from dust deposits suggest that the amplitude of wet-dry cycles in East and West Africa was low at the M2 event and no discernable climate shifts occurred (deMenocal, 2004).
However, based on a marine pollen record (ODP Site 658), Leroy and Dupont (1994) describe periods of marked aridity in northwestern African climate which occur during positive isotope excursions, including stage 132 (Tiedemann, 1991) which correlates to Stage M2 as proposed by Shackleton et al. (1995: S95). Increased accumulation of terrigenous matter in the East Atlantic also supports aridity on the northwest African continent at 3.3 Ma (Stein, 1985). Lisiecki and Raymo (2005) demonstrate that although the nature of the M2 excursion is not identical between the LR04 and S95 age models, the timing of M2 could be considered comparable between the two. Therefore, we can have confidence that aridity in northwestern Africa in any of our M2 scenarios would be consistent with the marine record.
When compared to modern, precipitation changes in northwest Africa are very small (b0.1 mm d −1 ; Fig. 3 and Fig. S5a) for the M2 glacial scenarios with CO 2 prescribed at 280 ppmv. In the PRISM-M2 M2Orbit 280 experiment, the implementation of PRISM3 ice sheets leads to a reduction in precipitation over land in northwest Africa, but an increase in coastal precipitation north of the Canary Islands. Whilst Small-M2 M2Orbit 280 produces widespread increases in precipitation (both coastal and continental), Medium-M2 M2Orbit 280 displays a reduction in coastal precipitation that begins to penetrate inland. Finally, Large-M2 M2Orbit 280 shows more widespread aridification over northwest Africa and the Mediterranean, which is consistent with the available pollen records (Fig. S5). Increasing aridification in northwest Africa alongside increases in NH ice sheet extent is also displayed in the 220 ppmv CO 2 experiments (Fig. S5b).
Climate data derived from fossil pollen assemblages from East Africa can be used to assess the plausibility of modelled ice scenarios. Bonnefille et al. (2004) record up to a 5°C mean annual temperature cooling and a 0.5-0.8 mm d − 1 increase in rainfall at Hadar, Ethiopia consistent with the global marine 18 O isotopic shift at 3.3 Ma. It should be noted that models of Pliocene climate often overestimate the temperature at Hadar and thus perform particularly poorly at this site . Additionally, there is disagreement regarding evidence for tectonic (elevation) changes in this region since the Pliocene (Aronson and Taieb, 1981;Redfield et al., 2003) which have not been explicitly taken into account in this study and therefore caution should be applied when judging the significance of the Hadar comparison.
Nevertheless, comparing HadCM3 predictions for the glacial sensitivity experiments to the standard Pliocene control (Pliocene Ctrl ) demonstrates the potential impact on temperature and precipitation in relation to the background warm Pliocene state (Fig. S6). The change in temperature at the Hadar site varies from − 2.37°C (PRISM-M2 M2Orbit 280 ) to − 5.39°C (Large-M2 M2Orbit 280 ), showing that even in the large ice sheet scenario, the temperature drop is reconcilable with the data. The pattern of precipitation change at Hadar is less consistent, with some regions displaying a precipitation reduction and others displaying an increase (Fig. S6e-f). At the Hadar site, the Medium-M2 M2Orbit 280 simulation produces the greatest precipitation increase (0.18 mm d − 1 ), although regional changes in the Large-M2 M2Orbit 280 ice sheet scenario are up to 0.7 mm d − 1 , which is more consistent with the Bonnefille et al. (2004) record. PRISM-M2 M2Orbit 280 displays a precipitation reduction (− 0.3 mm d − 1 ) when compared to Pliocene Ctrl .

Northeast Arctic Russia
High resolution pollen-based reconstructions from Lake El'gygytgyn in northeast Arctic Russia display a dramatic cooling event centred on MIS M2 . This cooling is combined with a reduction of detrital input to the lake and a termination of calcite formation (suggesting a reduction in weathering), which are both linked to the first development of permafrost in the region at MIS M2 . The significant decrease in the abundance of tree and shrub pollen around 3.3 Ma also indicates dryer conditions during this time Nowaczyk et al., 2013). Based on pollen spectra and the modern analogue approach (Melles et al., 2012), mean warm month temperature (MWMT) and mean annual precipitation (MAP) values detailed in Brigham-Grette et al. (2013) can be compared with the M2 glacial sensitivity simulations. Brigham-Grette et al. (2013) show a~5°C reduction in MWMT compared to the background Late Pliocene state. All M2 glacial sensitivity simulations display a reduction in summer (June, July, August; JJA) temperature at Lake El'gygytgyn relative to the Pliocene Ctrl of around 5°C (4.28 to 5.18°C), which shows that temperatures at this site are not significantly affected by the different prescribed ice sheets in the NH (Fig. S7). Actually, in the Large-M2 M2Orbit 280 simulation the changes in the land-sea mask in the Yakutia and Chukotka coastal regions promote increased continentality, displaying the warmest JJA temperatures at Lake El'gygytgyn of all of the M2 sensitivity experiments (Fig. S7d).
Estimated peaks of MAP at Lake El'gygytgyn are between 450 and 600 mm a −1 , but the record shows a prominent decrease towards modern levels of~280 mm a −1 coincident with cooling temperatures at M2 . Although HadCM3 predicts Pliocene precipitation rates which are broadly in line with the Lake El'gygytgyn record (500 mm a −1 ), the change shown in the glacial scenarios is at most 140 mm a − 1 (Large-M2 M2Orbit 280 ) and thus the precipitation levels at M2 remain wetter than the proxy-record might suggest.
Additional simulations presented in Brigham-Grette et al. (2013) using an atmosphere-slab-ocean climate model suggest that a full Last Glacial Maximum (LGM;~120 m equivalent sea level) is inconsistent with the proxy record at Lake El'gygytgyn, with the regional climate being too cold and too dry. The results we have presented here would suggest that an ice sheet half the size of the LGM, could still be consistent with the broad climate changes at MIS M2. The sensitivity seen in previous models may reflect the location of large ice masses in the ICE-4G reconstruction of the LGM (Justino et al., 2006), particularly in the Russian Arctic or Cordilleran ice masses.

Southern Ocean
Increased production of IRD (Passchier, 2011) suggests that ice sheet behaviour in the SH was also changing during MIS M2. Recent evidence from the AND-1B core has pointed towards an expansion of the Antarctic ice sheet on to the continental shelf at 3.3 Ma associated with reduced SSTs (− 2.5°C), increased sea-ice extent and altered Southern Ocean circulation (McKay et al., 2012). In the M2 glacial scenarios tested here, the areal extent of the Antarctic ice sheet has not been altered to exceed modern, however, there are parts of the ice sheet in the Medium-M2 M2Orbit 280 and Large-M2 M2Orbit 280 scenarios that have a higher elevation compared with the present day (Fig. 2). Both these experiments produce a 2-3°C cooling in mean annual coastal temperatures around the Ross Sea in comparison to the Pre-Industrial control, which is broadly consistent with the data (Fig. 3). An even greater cooling (N7°C) is predicted if the M2 model scenarios are compared with the Pliocene Ctrl (with no WAIS and a reduced EAIS; Dowsett et al., 2010). Predicted winter sea-ice extent around Antarctica (relative to present day) also increases in the Small-M2 M2Orbit 280 , Medium-M2 M2Orbit 280 and Large-M2 M2Orbit 280 ice sheet scenarios. Model-predicted sea ice increase is associated with a northward migration of the Southern Ocean westerlies and an increase in the overall westerly wind strength (Fig. S2). For example, there is up to a 3 m s − 1 increase in the mean wind strength in the Atlantic sector of the Southern Ocean in the Large-M2 M2Orbit 280 scenario (Fig. S2f). Such wind strength increases are consistent with the prescription of a larger Antarctic ice sheet which increases the temperature gradient (by up to 3.9°C) and the pressure gradient in the SH in Large-M2 M2Orbit 280 . Broadly, the simulations in which the Antarctic ice sheet is increased in volume (height) relative to modern are in agreement with the general conditions in the Southern Ocean at M2 according to McKay et al. (2012). It is probable that expanding the AIS further onto the continental shelf in HadCM3 would lead to an amplification of the trends displayed in the experiments here (e.g. greater increases in sea ice and more intense westerlies). However, it would be impossible to determine the most plausible ice configuration on Antarctica from a DMC focussed solely on information from one site.

Broader implications and future work
Both terrestrial and marine proxy records provide consistent evidence for significant changes in climate and continental glaciation around MIS M2. The simulations presented here provide the first assessment of the plausibility of larger-than-modern ice sheets in the NH and SH at MIS M2. Outcomes of our data/model comparison reinforce the possibility of larger ice sheets in the NH and SH during M2. Specifically, most available proxy data appear to be consistent with either a Medium (~34 m SLE) or a Large (~52 m SLE) ice extent in the NH.
De Schepper et al. (2013) invoke an increase in Pacific-to-Atlantic flow due to an open CAS to explain changes in reconstructed SSTs and currents at the M2. Therefore, a logical next step would be to test the hypothesis for the mechanism of ice sheet growth (the temporary opening of the CAS) within the modelling framework presented here. Climate and ice sheet modelling presented by Lunt et al. (2008aLunt et al. ( , 2008b investigating the initiation of NHG suggest that the closing of the CAS was not a sufficient trigger for ice growth on Greenland, however, an open CAS did lead to marginally more ice growth on North America (Lunt et al., 2008b).
Finally, the experimental design of this study was such that an increase in ice area in Antarctica was not explicitly considered. However in light of recent evidence from the Southern Ocean (McKay et al., 2012), a future study should test the compatibility of global data at M2 with the hypothesis that significant ice changes comes from Antarctica and not the NH.

Conclusions
This paper presents the results from the first coupled atmosphereocean climate model experiments specifically targeting the simulation of Marine Isotope Stage M2. Plausible ice sheet scenarios (based on Quaternary ice sheet reconstructions) have been implemented in HadCM3 to test the possibility of larger-than-modern ice sheets in the Northern Hemisphere at~3.3 Ma. We have found that: • Large and Medium ice sheet scenarios representing 60 m and 40 m of sea level rise compared with present day, induce a cooling of global annual mean climate of between 0.39°C and 0.29°C under CO 2 levels of 280 ppmv and an M2 orbit. • MIS M2 glacial scenarios are also globally drier than the pre-industrial control (Pre-Ind Ctrl ) and the Pliocene control (Pliocene Ctrl ). • Larger-than-modern ice sheets are not conclusively inconsistent with palaeoenvironmental proxy data relating to surface currents, salinity and sea surface temperatures in the North Atlantic. • An intensitfication of aridification in northwest Africa with larger ice sheets is consistent with marine pollen data from offshore Africa 3.3 Ma. • In Northeast Arctic Russia summer temperatures at Lake El'gygytgyn are insensitive to the range of ice sheet sizes we have imposed in the wider NH. This means that all ice sheet configurations tested here are potentially compatible with pollen-based reconstructions of mean warm month temperature. • Prescribing a larger-than-modern Antarctic ice sheet, sea surface temperatures, winter sea ice extent and the prediction of the Southern Ocean westerlies are all in agreement with proxy-based data.