The challenge of simulating the warmth of the mid-Miocene climatic optimum in CESM1

The mid-Miocene climatic optimum (MMCO) is an intriguing climatic period due to its above-modern temperatures in mid-to-high latitudes in the presence of closeto-modern CO2 concentrations. We use the recently released Community Earth System Model (CESM1.0) with a slab ocean to simulate this warm period, incorporating recent Miocene CO2 reconstructions of 400 ppm (parts per million). We simulate a global mean annual temperature (MAT) of 18C, ∼ 4C above the preindustrial value, but 4 C colder than the global Miocene MAT we calculate from climate proxies. Sensitivity tests reveal that the inclusion of a reduced Antarctic ice sheet, an equatorial Pacific temperature gradient characteristic of a permanent El Niño, increased CO2 to 560 ppm, and variations in obliquity only marginally improve model–data agreement. All MMCO simulations have an Equator to pole temperature gradient that is at least∼ 10C larger than that reconstructed from proxies. The MMCO simulation most comparable to the proxy records requires a CO 2 concentration of 800 ppm. Our results illustrate that MMCO warmth is not reproducible using the CESM1.0 forced with CO 2 concentrations reconstructed for the Miocene or including various proposed Earth system feedbacks; the remaining discrepancy in the MAT is comparable to that introduced by a CO 2 doubling. The model’s tendency to underestimate proxy derived global MAT and overestimate the Equator to pole temperature gradient suggests a major climate problem in the MMCO akin to those in the Eocene. Our results imply that this latest model, as with previous generations of climate models, is either not sensitive enough or additional forcings remain missing that explain half of the anomalous warmth and pronounced polar amplification of the MMCO.


Introduction
The mid-Miocene climatic optimum (MMCO 17-14.50Ma) (Zachos et al., 2008) is a period in Earth's history in which temperatures were significantly warmer in the deep ocean and in mid-to-high latitudes (Böhme et al., 2007;Pound et al., 2012;Zachos et al., 2008;Shevenell et al., 2008).These warm extratropical temperatures have been hard to reconcile with reconstructed below-modern tropical sea surface temperature (SST) records and boron and alkenone CO 2 reconstructions of 200-280 ppm (parts per million) levels (Pagani et al., 2005;Pearson and Palmer, 2000).
Recent re-evaluation of the proxy records has led to advancement in our understanding of MMCO warmth.First, the MMCO tropical SST records showing below-modern levels (Savin, 1977;Nikolaev et al., 1998;Bojar et al., 2005) are now understood to have a cool diagenetic bias (Stewart et al., 2004).Excluding these records indicates that tropical SSTs in the Miocene were above modern climate (Shevenell et al., 2004;You et al., 2009;LaRiviere et al., 2012;Zhang et al., 2013).Second, recent leaf stomatal studies reconstruct CO 2 concentrations at the MMCO to be 400-500 ppm (Kürschner et al., 2008) and these results have been confirmed in boron isotope-based reconstructions (Foster et al., 2012) and updated alkenone reconstructions (Zhang et al., 2013).
Published by Copernicus Publications on behalf of the European Geosciences Union.

A. Goldner et al.: Simulating warmth of the mid-Miocene climate optimum in CESM1
Nevertheless, even with higher CO 2 concentrations MMCO warming has been difficult to reproduce in an intermediate complexity Earth system model (Henrot et al., 2010), atmosphere and slab-ocean models (Tong et al., 2009;You et al., 2009), and fully coupled atmosphere-ocean models (Herold et al., 2011;Krapp and Jungclaus, 2011).For example, Herold et al. (2011) found that the Community Climate System Model (CCSM3.0)was ∼ 10 • C too cold compared to proxy records in high latitude regions like Alaska and Antarctica.In this study, we implement boundary conditions from Herold et al. (2011) within the National Center for Atmospheric Research (NCAR) Community Earth System Model (CESM1.0)using the Community Atmosphere Model (CAM4) framework to simulate the MMCO.This allows for a clean comparison with previous simulations done with CCSM3.0, using a latest generation model included in the Coupled Model Intercomparison Project (CMIP5).
To explore if the modeling framework is able to match MMCO warmth we conduct a pointwise model-data comparison using proxy records compiled for the MMCO (Tables S1  and S2).The MMCO is a good choice for climate model validation because the continental configuration is relatively close to modern (Herold et al., 2008) although differences exist (Potter and Szatmari, 2009).Additionally, the CO 2 levels during the MMCO are in the range of values for the next century, and paleoclimate records are better constrained compared to earlier warm periods such as the Eocene (∼ 56-33.9Ma) where there is large uncertainty in the CO 2 (Pagani, 2002;Pearson and Palmer, 2000;Royer et al., 2012) and temperature records (Huber and Caballero, 2011).
This makes exploring the MMCO especially important for understanding how the climate system responds to different boundary condition changes in a world with similar CO 2 and temperatures to those projected for the future.Additionally, because there is a possible decoupling of atmospheric CO 2 (LaRiviere et al., 2012) with the pronounced warmth during the MMCO, this study aims to explore many of the non-CO 2 forcings that could have operated during this period.Explaining the warmth is complicated and will likely involve many different alterations to boundary conditions as well as understanding nonlinear feedbacks to the imposed changes (Holbourn et al., 2005;Shevenell et al., 2008;Lyle et al., 2008).
CESM1.0 improvements in comparison to CCSM3.0 include the revised Zhang and McFarlane deep convection scheme in CAM4 allowing for parcel mixing (Zhang and Mc-Farlane, 1995;Neale et al., 2008Neale et al., , 2010)).The atmospheric component is also able to resolve tropical circulation better, including improvements in the Hadley circulation and tropical convection (Richter and Rasch, 2008).It also now includes improved modeling of the seasonal cycle (Bitz et al., 2012) and the model has a reduction in high latitude low cloud biases (Vavrus and Waliser, 2008).The slabocean configuration includes fully dynamic sea ice, whereas CCSM3.0 had only thermodynamic sea ice, and the sea-ice model includes updates in the scattering parameterizations, enhancing the realism of snow albedo (Briegleb and Light, 2007).Other CESM1.0 model improvements in comparison to older generation models are described in previous work for modern (Bitz et al., 2012;Neale et al., 2010) and for paleoclimate (Shields et al., 2012;Goldner et al., 2013).

Experimental design
The control preindustrial (PI) simulation employs the modeling components described above in standard configuration and with CO 2 concentrations set at 287 ppm (Table 1).The slab-ocean forcing file for the PI case has heat fluxes, salinity, and density inputs from a fully coupled atmosphere, ocean, ice, and land simulation (Bitz et al., 2012).Additionally, we run a PI simulation at 400 ppm CO 2 (PI400) to compare with our MMCO simulation (also at 400 ppm CO 2 ).This high CO 2 PI configuration allows us to isolate the temperature effect of including MMCO boundary conditions at constant CO 2 .
The MMCO simulation has vegetation cover (Herold et al., 2010) and topography described in Herold et al. (2011).In comparison to modern vegetation, the prescribed MMCO vegetation has reduced ice coverage over Antarctica (Pekar and DeConto, 2006), while Greenland is ice free.In the tropics and midlatitudes the desert regions are replaced with savannah and all C4 grasses are replaced with C3 grasses (Herold et al., 2010).Additionally, mid-to-high latitude regions have an abundance of temperate, broadleaf evergreen biomes in agreement with previous reconstructions that depict warmer and wetter conditions (Wolfe, 1985).
Previous slab-ocean and atmosphere MMCO simulations have been conducted within the CCSM3 framework (Tong et al., 2009;You et al., 2009), but here we improve upon their methodology by using ocean heat fluxes derived from a coupled ocean-atmosphere simulation of the Miocene.To create the Miocene slab-ocean forcing file we use a previous CCSM3.0Miocene simulation conducted by Herold et al. (2012).We include fluxes averaged from the last 100 yr of  (Herold et al., 2012).We admit that this simulation has CO 2 levels above reconstructed values, but this is a better background climate state for conducting our simulations because of CCSM3's known low sensitivity to CO 2 forcing (Otto-Bliesner et al., 2006).
The slab-ocean model employed in this study is a mixedlayer thermodynamic model (i.e., excludes dynamics) and is considered to be thoroughly mixed at all the ocean depths (Bitz et al., 2012).The slab-ocean model is forced with ocean heat convergence and mixed layer depths from a fully coupled, equilibrated Miocene simulation and its purpose is to recreate the equilibrated state including SST and sea ice of the fully coupled climatology.Total heat into the slab ocean is the combination of net long, short, sensible and latent heat from the atmosphere, upwelling heat from the deep ocean (prescribed in this study), and sensible and latent heat fluxes from ice and snow (Bitz et al., 2012).
Our use of the slab-ocean model in this study, as opposed to a fully dynamic ocean model, is justified given that (1) we are interested in simulating a large number of sensitivity experiments that demand already intensive computational resources.(2) Experience from modern and Eocene studies show that this slab-ocean approach produces very similar answers to those from coupled models (Gettelman et al., 2012;Bitz et al., 2012).(3) We can run the slab-ocean simulations with higher resolution in the atmosphere (1.9 • × 2.5 • ) than is standard for most paleoclimate studies because of the reduced computational requirement.(4) We explore imposing a permanent El Niño, which has been difficult to reproduce in fully coupled models (Haywood et al., 2007;Galeotti et al., 2010;Huber and Caballero, 2003).
Using slab fluxes from CCSM3.0 is not an issue because we find no substantial differences in SST (Fig. S1) or climate between CCSM3.0 and CESM1.0 for deep paleoclimate simulations such as the Eocene, as the ocean component biases are very similar between the two modeling frameworks (Danabasoglu et al., 2012).
The permanent El Niño imposed in some of our simulations is comparable to proxy based SST reconstructions for the equatorial Pacific (Wara et al., 2005;Dekens et al., 2008) and similar to previous studies which have prescribed this type of SST distribution (Shukla et al., 2009;Vizcaíno et al., 2010).The imposed permanent El Niño has a similar structure to the "canonical" El Niño; unlike the "Modoki-type" which has the largest temperature increases in the central Pacific (Ashok et al., 2007).In the paper we will refer to permanent El Niño as El Padre (Shukla et al., 2009).
The entire suite of simulations included in our analysis are listed in Table 1.The simulations conducted are run for over 60 yr with the last 20 used for analysis and are well equilibrated as evidenced by the radiative balance statistics found in Table S3.Additional MMCO CO 2 sensitivity experiments were run at 560 ppm CO 2 (MMCO560) to account for the uncertainty in Miocene CO 2 reconstructions and the model-data comparisons for this experiment is described in Fig. 3b.We also run a simulation at 800 ppm (MMCO800) CO 2 (Fig. 3c) (Table 1) to explore a wide range of CO 2 values although we note that this is well outside the range of the reconstructed CO 2 levels for the MMCO.

MMCO terrestrial and sea surface temperature compilation
For the model-data comparison we update the compilation of terrestrial and SST proxy records described in Pound et al. (2012), Herold et al. (2011), and others (Tables S1, S2).
We present the longitudinal and spatial distribution of the proxy records in Fig. 1.The proxy reconstruction spans over the MMCO (17-14.50Ma), however, because of the sparseness of data over this period we include records that have an average age between 20 and 13.65 Ma, where they fill spatial gaps (i.e., Southern Hemisphere).This data compilation can be used as a reference data set for future MMCO model-data comparisons.
We update the minimum error in our compiled terrestrial proxy records for a number of reasons.Firstly, recent work suggests that for physiognomic leaf-climate methods there should be a minimum error of ±5 • C (Royer, 2012).Secondly, studies have suggested that there is large uncertainty in estimating MAT (Grimm and Denk, 2012) using the coexistence approach (Mosbrugger and Utescher, 2007).For  S1,S2.b) The spatial distribution of the terrestrial and SST proxy records used in the model data comparisons overlain onto the Miocene topography (Herold et al., 2008).S1 and S2.
(b) The spatial distribution of the terrestrial and SST proxy records used in the model-data comparisons overlain onto the Miocene topography (Herold et al., 2008).
our intended purposes increasing the minimum proxy record uncertainty should make matching the simulations more obtainable.If our model still fails to match proxy data even with generous error bars this merely proves our main results further.
To address the temporal variability of including records across a long timescale we calculate approximate global mean temperatures, relative to the MMCO, based on the known temperature relationship to oxygen isotopes and applying a factor of 0.5 to account for polar amplification using a stacked oxygen isotope record for the early-to-middle Miocene (Zachos et al., 2008).We find that the difference in temperature before and after the MMCO is −0.84 and −0.99 • C, respectively.Thus, our error bars encompass the likely variation in global mean temperature that may be present in the records used here.This interpretation is buttressed by more sophisticated studies like that of Liebrand et al. (2011).
The SST records are compiled from available published data in the literature and we describe these records in detail in Table S2.We separate out some tropical SST records that may have a diagenetic bias (Sexton et al., 2006;Huber, 2008).Tropical SSTs are few and far between for the MMCO, but more common in the mid-to-late Miocene, thus we may omit proxy records from over almost half the surface area of the planet (30 • N and 30 • S) or utilize data from intervals slightly outside the MMCO.Because there is a lack of tropical SST data points for the MMCO we compile SSTs from the late Miocene and justify this based off the minimal change between middle and late Miocene SSTs at other locations (LaRiviere et al., 2012).Given that the Pliocene tropical SSTs were ∼ 4-6 • C (Brierley et al., 2009;Dekens et al., 2007;Ravelo et al., 2006;Fedorov et al., 2013) above modern SSTs in the tropical and sub-tropical upwelling zones and the late Miocene was ∼ 7-9 • C above modern SSTs (LaRiviere et al., 2012) it is reasonable to conjecture MMCO tropical SSTs were this warm or warmer.This characterization is supported by the recently published SST data set of Zhang et al. (2013).Either approach introduces potential errors in interpretation and here we choose to utilize SST estimates in data sparse regions that lie generally within the early-tomiddle Miocene, but may be outside the MMCO.Our updated minimum error bars are large enough to encompass the temporal variation in these records.
Previous work has discussed the importance of including orbital variations when quantifying uncertainty in modeldata comparisons (Haywood et al., 2013).To quantify the possible error introduced by aliasing of orbital variability in our interpretation of model-data mismatch, we conduct two sensitivity experiments varying obliquity to minimum and maximum Miocene values (22 and 25 • respectively).We run each of the sensitivity tests for 60 yr and we then calculate using an average of the last 20 yr the maximum and minimum model-derived temperatures at each proxy location from both extreme orbit simulations and use this absolute anomaly as an estimate of orbitally induced variance.These maximum and minimum values are plotted as vertical error bars on the modeled MAT in our pointwise model-data comparisons (Figs.2-6).

Proxy derived MAT value
To determine the difference in global MAT between Miocene and preindustrial climate we take the proxy records and perform a pointwise anomaly of proxy-derived MAT compared to modern observed MAT at paleolatitudes and paleolongitudes.We split the resulting anomalies into tropical (30 • N-30 • S) midlatitude (30 • N/S-60 • N/S) and polar (60 • N/S-90 • N/S) regions and conduct a weighted average anomaly over each latitudinal region.This latitudinal binning and area weighting addresses issues of having more proxy records in certain regions (i.e., the midlatitudes).The error bars for the proxy derived MAT were calculated to include two types of error.First, we assume a ±5 • C random error on each proxy estimate.Second, a random error is introduced across the spatial domain calculated from the standard deviation over samples within the three binning regions.The two errors are combined and normalized by the square root of the total number of proxy records.Using proxy records for the MMCO (Table S1 and S2) we calculate a global MAT change of ∼ 7.6 • C ± 2.30 (we report two standard errors from the mean) compared to PI.The proxyderived temperatures compared against modern observations (ECMWF 40 Year Re-analysis Project) is 6.8 • C ± 2.20 as there is ∼ 1.0 • C of warming between modern observations and PI climate.
To validate our approach for estimating proxy derived MAT we calculate a resampled MAT using our methodology and compare against a globally weighted MAT (we will call this true MAT) from both model runs and modern observational data sets.The globally weighted true MAT value of the MMCO simulation is 18.00 • C (Table 2) whereas our calculation for MAT resampled over the proxy record regions using the methodology from above is 17.12 • C. The calculated standard error from the mean including proxy record uncertainty is 1.33 • C, which illustrates that our resampled MAT value is well within the calculated standard error.We also calculate the resampled MAT using modern observations and with other Miocene simulations and find that all the resampled MAT estimates fall within two standard errors of the true MAT.For all intended purposes we are confident that our approach for reconstructing global MAT from our proxy record compilation is a valid estimate.

MMCO simulation compared against the proxy records
The MMCO simulation is 4.04 • C warmer than the control PI simulation, but the simulation is about 4 • C cooler than globally averaged MMCO proxy temperature reconstructions (Table 2).The MMCO simulation generally captures the tropical and midlatitude temperature distribution of the proxy records, but fails to achieve above-freezing temperatures in the high latitudes (Fig. 2b, Table 3).The nature of this discrepancy can be clarified by examining the Equator to pole surface temperature gradient.It is 17 • C larger in the MMCO simulation than in the proxy records (Table 2).
Using the methods described in Lunt et al. (2012), the Equator to pole temperature gradient is calculated by averaging the mean annual temperatures over the absolute latitudes of 60-80 • minus 0-30 • ; except here we use 80 • because this is the maximum latitudinal extent of proxy records.Addition-  S1, S2).Horizontal error bars are the modeled pointwise maximum and minimum temperatures from the extreme obliquity simulations (see Methods Sect.2.3) and methodological error is plotted as the vertical error bars.The best fit line (black dashed) is weighted to include proxy uncertainty and is fitted across all points.The weighting for each proxy record is calculated by 1/(error 2 ).The y intercept and slope are reported in Table 2.
ally, an error weighted best fit line for the pointwise comparison reveals a root mean square error (RMSE) of ∼ 6 • C and y intercept of 9.61 • C, and the slope of the regression line is 0.62 (Table 2).In summary, the MMCO simulation (at 400 ppm CO 2 ) is unable to produce high latitude warmth or a sufficiently warm global mean temperature compared to the paleotemperature records.

Effect of MMCO boundary conditions and CO 2 sensitivity experiments
We find that our MMCO simulation is 2.43   S1, S2).2.43 • C of the temperature difference between our MMCO and PI simulations are a result of changes in continental positions, topography, and vegetation.This change is consistent with late Miocene modeling, which finds 3.0 • C of warming due to changes in vegetation and topography (Knorr et al., 2011).
A CO 2 sensitivity experiment run at 560 ppm CO 2 (above most reconstructed CO 2 records) is also too cold at high latitudes compared to proxy records (Fig. 3b) and the Equator to pole temperature difference is still too large by ∼ 13 • C (Table 2).This simulation has a global MAT 5.89 • C higher than the control PI simulation, and is ∼ 2 • C colder than the proxy-derived global MAT.The error weighted best fit line for the MMCO560 pointwise comparison gives a y intercept of ∼ 7.93 • C, but the calculated RMSE is still 5.7 • C (Table 2).The MMCO800 simulation has a MAT 7.26 • C above PI (Table 2), which is our best comparison with the proxy derived MAT value.The error weighted best fit line is closer to the one-to-one line and has our best y intercept of 6.62.(Fig. 3d).Overall, MMCO800 matches the proxy  2.
compilation the best and we use this comparison to prove that matching global MMCO warmth can be accomplished, but at CO 2 concentrations approximately twice that reconstructed from proxies.These results are very similar to those found in the Eocene (Huber and Caballero, 2011;Lunt et al., 2012) Below, we test hypotheses that have been proposed to explain Miocene warmth, with the goal of improving the model-data comparison without having to increase CO 2 above reconstructed levels.

Reducing Antarctic ice-sheet volume
Recent work estimates the volume of the middle Miocene Antarctic Ice Sheet (AIS) to be ∼ 30-50 % less than modern (Shevenell et al., 2008).Consequently the Herold et al. (2008) reconstruction for AIS elevation and extent is likely too large (Fig. 4a).To correct this, we utilize a new AIS reconstruction derived from a fully interactive terres-trial ice and atmosphere model (D.Pollard, personal communications, 12 April 2012) (Fig. 4b).This AIS is half the volume of that used in Herold et al. (2011) (Fig. 4a) and is within the range of estimates from proxy records (Pekar and DeConto, 2006;Billups and Schrag, 2003).We also reduce the area of glacier albedo over Antarctica by half and replace it with a combination of unvegetated and tundra-like land cover.We introduce this new AIS topography and vegetation cover (Fig. 4b) into the MMCO boundary conditions described in Herold et al. (2008) and denote this simulation LOW AIS.The difference in surface albedo over the AIS between these two simulations ends up being similar as snow (also with a high albedo) ends up covering the areas that were once glacier because Antarctica stays below freezing year round.
The LOW AIS simulation is 4.15 • C warmer than PI and 0.10 • C warmer than the previously described MMCO simulation with a high AIS (Fig. 4d).The global temperature response to adding or removing the AIS in the Miocene has not been explicitly calculated and we find that there is no substantial global mean temperature impact from decreasing  S1,S2).The best fit line (black dashed) is weighted to include error uncertainty is fitted across all points and the y-intercept and slope are reported in Table 1.S1, S2).The best fit line (black dashed) is weighted to include error uncertainty and is fitted across all points and the y intercept and slope are reported in Table 2.
the size of the AIS.This result is consistent with previous work, which found a small global temperature response to adding and removing the AIS in the Eocene (Goldner et al., 2013).Although recent coupled MMCO simulations have found warmer and wetter conditions regionally over Europe due to the reducing ice extent in Antarctica highlighting the importance of including ocean feedbacks for resolving regional temperature distributions (Hamon et al., 2012).The temperature difference between LOW AIS and the MMCO simulation is largest over Antarctica (Fig. 4c) because of the imposed elevation and surface albedo changes.Although lowering the AIS warms the Antarctic continent, the Miocene LOW AIS simulation results in negligible improvement in matching proxy records elsewhere in the high latitudes (Table 3).A slight warming occurs in the Ross Sea between the LOW AIS simulation and the MMCO simulation, but overall there is minor improvement in the model-data comparison (Fig. 4d) by lowering the height and reducing glacier extent of the AIS (Fig. 4b).

El Padre
It has been hypothesized that pre-Quaternary climates were characterized by a reorganization of tropical oceanatmosphere circulation inducing a permanent El Niño SST distribution (Philander and Fedorov, 2003;Lyle et al., 2008;Ravelo et al., 2006), which has been called El Padre (Shukla et al., 2009).A reduced temperature gradient in the eastern equatorial Pacific (EEP) should induce high latitude warming in Alaska and other high latitude regions, because this is a standard teleconnected response during modern El Niños (Molnar and Cane, 2007).Prior modeling studies have demonstrated the effectiveness of this mechanism (Barreiro et al., 2006;Vizcaíno et al., 2010;Bonham et al., 2009;Haywood et al., 2007;Goldner et al., 2011), although no modeling study has explicitly studied its impacts with realistic MMCO boundary conditions.To explore the impacts of an El Padre SST anomaly in our simulations, we take the heat convergence and mixed Pointwise EP case global mean MAT compared against the proxy record MAT (˚C).These are the same terrestrial and SST records and error bars described in Figure 1.The best fit line (black dashed) is weighted to include error uncertainty is fitted across all points and the y-intercept and slope reported in Table 1.  2.
layer depths derived from a fully coupled Miocene simulation (Herold et al., 2012) and zonally average these quantities across the equatorial Pacific (10 • N and 10 • S of the Equator).We introduce the zonally averaged ocean heat convergence and mixed layer depths into a new slab-ocean forcing file and simulate the MMCO with a low AIS at 400 ppm CO 2 .The El Padre SST anomaly is presented in Fig. 5a and will be referred to as EP.We are confident the CAM4 CESM1.0 framework reproduces modern day observational teleconnection patterns induced by El Niño forcing, as described in detail in other studies (Wang et al., 2013;Shields et al., 2012).
Although an interesting question for past warm periods like the MMCO is how these global and regional responses to ENSO (El Niño-Southern Oscillation) have varied throughout geologic time, as modeling of the late Miocene has shown that ENSO teleconnections can be modified from modern teleconnections (Galeotti et al., 2010;Tang et al., 2014).
In the EP simulation high latitude regions warm, especially Alaska and Antarctica (Fig. 5a).The pointwise model-data comparison for the EP simulation is plotted in Fig. 5b.This simulation is ∼ 4.6 • C warmer in global mean than the PI simulation and ∼ 0.5 • C warmer than the MMCO and LOW AIS simulations.Warming due to adding El Padre is largest in regions where the model previously performed the worst (Fig. 5a).Roughly 2 • C of warming occurs in Alaska, but the simulation is still ∼ 8.5 • C too cold in this region (Table 3) and still has a ∼ 13 • C larger Equator to pole surface temperature gradient compared to the proxy records (Table 2).Imposing an El Padre illustrates a mechanism capable of warming the high latitudes without elevating CO 2 , consistent with the results of LaRiviere et al. ( 2012); Sriver and Huber (2010); Brierley et al. (2009).Nevertheless this change does not reconcile the warmth of the MMCO, as temperatures are still ∼ 2 • C too cool globally and ∼ 8.5 • C too cool in the high latitudes.
Adding EP and increasing obliquity to 25 • results in a simulation that is 5.64 • C warmer than PI (Fig. 6).These boundary condition changes in conjunction with one another are not clearly detected in proxy reconstructions and we present this simulation as an example of including a suite of forcings aimed at increasing temperature in high latitude regions.Although higher than modern obliquity is reconstructed in orbital configurations calculated for the MMCO (Laskar et al., 2004), this MAT anomaly compared to PI is similar to the warming found in the MMCO560 simulation.The MMCO560 simulation does not include any of the boundary condition changes aimed at increasing high latitude warmth.Interestingly the EP, AIS, and obliquity forcing results in a 4 • C improvement in simulating the Equator to pole temperature gradient compared to MMCO560 (Table 2).Both comparisons are too cold compared to the proxy derived global MAT value as matching the proxy records in high latitudes requires a CO 2 concentration of double what is predicted in the reconstructions.

Comparison with previous MMCO CCSM3.0 simulations
The most comparative study to the experiments presented here are the CCSM3.0MMCO simulations described in Herold et al. (2011) (Table 2).In the present study, the Miocene simulations are ∼ 2.0 • C warmer than the Miocene CCSM3.0 simulations (Herold et al., 2011) at the same CO 2 levels.This temperature difference is explained in large part because CAM4 CESM1.0 is a more sensitive model to background CO 2 concentrations.The older CCSM3.0 has a 2.5 • C change in global mean surface temperature to a doubling of CO 2 (Kiehl et al., 2006), while CSEM1.0CAM4 has a 3.5 • C temperature change to a doubling of CO 2 (Gettelman et al., 2012), roughly a 1   2.
less warming in the upper troposphere compared to the surface warming at varying latitudes and a reduction in low level cloud fraction (Bitz et al., 2012).

Comparison with other fully coupled MMCO simulations
Krapp and Jungclaus (2011) simulated the MMCO and found a MAT of 17.1 and vegetation across the MMCO and these simulations are too cold in the midlatitudes compared to the records.They simulate warming of 2.9 and 3.4 • C above PI when CO 2 concentrations are increased to 500 ppm and vegetation is altered respectively, which is half of the temperature change needed to explain the proxy derived MAT.

Further implications of the model-data mismatch
Additional MMCO simulations exploring precession (Sloan and Huber, 2001;Lawrence et al., 2003), eccentricity (Westerhold et al., 2005), and vegetation (Knorr et al., 2011) could be important in our understanding of MMCO warmth.Also, because we have fixed ocean heat transport (OHT), using a dynamic ocean could help explain missing warmth.However, on the contrary, every study that has been performed studying the sensitivity of zonal mean OHT to changes in a wide range of boundary condition alterations occurring throughout the Cenozoic robustly shows small changes in OHT, none of which go as far as to explain high latitude warmth (Huber and Sloan , 2001;Krapp and Jungclaus, 2011;Zhang et al., 2013;Haywood et al., 2013).We point out that previous fully coupled oceanatmosphere simulations at 560 ppm CO 2 were unable to reproduce MMCO warmth (Herold et al., 2011).In fact, the CCSM3.0MMCO simulation at 560 ppm CO 2 performs worse against the proxy records than our MMCO CAM4 CESM1.0 simulation at 400 ppm CO 2 (Table 2).We also reiterate that the temperature effect of including MMCO boundary conditions induces 2.43 • C of warming compared to the PI400 simulation.This is roughly a third of the warming needed to explain the MMCO warmth of ∼ 7.6 • C ± 2.30.
The issue of some terrestrial records, especially in mid-tohigh latitude regions, having difficulty in matching modeling simulations has been pointed out more extensively in a model-data comparison of the Eocene conducted by Huber and Caballero (2011).Possible reasons behind the modeldata mismatch include undersampling for leaf physiognomic techniques (Wilf, 1997), skewness of high latitude temperatures due to the "toothiness" (Boyd et al., 1994), and other systematic biases mentioned below (Burnham et al., 1989;Spicer et al., 2005;Peppe et al., 2010).Additionally, due to the length of the period we explore there is the question of whether the record is recording a seasonal or MAT signal.All of this uncertainty is additional justification for why we have chosen to include larger error bars on all our terrestrial proxy records (Royer, 2012;Grimm and Denk, 2012).
Other factors related to the model boundary conditions could be contributing to the modeling simulations being too cold compared to the terrestrial records.The modeling framework is limited to the fact that paleoelevation in some areas is simplified due to resolution and also poorly constrained in other regions (Herold et al., 2008).This could significantly affect MAT temperatures in the midlatitude regions depending on the imposed elevation especially if specific locations are too high.Additionally, further altering orbital variations (Haywood et al., 2013) and Miocene vegetation cover (Dutton and Barron, 1997) in the high latitudes can have a significant effect on the regional temperature over the land surface.Finally, we point out that a generally good match for terrestrial records occurs at 800 ppm CO 2 (Fig. 3c,  d) as there is enhanced high latitude warmth compared to the tropical regions.This again highlights that we are missing a forcing equal to a doubling of CO 2 to explain MMCO warmth.Thus, finding the right combination of boundary condition changes and the right model framework will continue to improve the model-data mismatch regionally.

Conclusions
Paleoclimate modeling studies need to conduct a pointwise model-data comparison to be confident that their modeling results match proxy records and consequently we will make the presented MMCO temperature data set available for these types of comparisons.Simulating the MMCO at 400 ppm CO 2 using the CAM4 CESM1.0 framework produces a significant model-data mismatch in global MAT and in high latitudes.The discrepancy in the MAT comparison is equal to that introduced by a full doubling of CO 2 , as the model matches the data best at 800 ppm CO 2 .A similar conclusion about climate model sensitivity to background CO 2 forcing was reached based on fully coupled ocean-atmosphere Eocene simulations where a CO 2 levels nearly double the reconstructions was required to match the proxy records (Huber and Caballero, 2011).It is interesting to note that the reconstructed CO 2 used in this study of 400 ppm is equivalent to the concentration used in simulations of the Pliocene, where global temperatures were not as warm as the Miocene.
Including two of the most discussed Earth system feedbacks (El Padre and reduced ice volume) had small impacts on improving the model predictions even when we included uncertainty associated with time varying and possible aliasing of orbital forcing.Like previous fully coupled atmosphere ocean efforts (Herold et al., 2011;Krapp and Jungclaus, 2011), matching proxy records at the MMCO is challenging even in the latest generation of models and using a model with a climate sensitivity near the median of Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (AR5) estimates (Andrews et al., 2012).Given the variety of methods used for reconstructing Miocene climate (Tables S1, S2), we are confident in the broad trends reflected in the proxy record.Thus, explaining the warming will require additional incremental changes in boundary conditions (such as an even higher CO 2 ), a more sensitive model to background CO 2 concentrations, and/or identification of some -as yet unknown -process or forcing that accounts for almost half of the difference in temperature between today and the MMCO.We understand that we are also limited in this study to the use of one model configuration and our conclusions should be understood within this context.However, we have tested the most common (and not so common) boundary condition changes understood to have a likely impact on MAT and even with all of these changes there is a negative bias in the model.Furthermore, the CESM exhibits biases identified in numerous other paleoclimate models, adding confidence that the climate simulated here is not fundamentally different to what would be expected from other models.
Although some terrestrial CO 2 proxies suggest CO 2 was higher than 500 ppm, this would not solve the data-model mismatch, as increasing CO 2 past 560 ppm would likely make the tropics too warm (e.g., Fig. 3b, d).Ultimately, our inability either to identify a missing paleoclimate forcing or formulate models with sufficient positive feedbacks to recreate substantial increases in global mean temperature with strong polar amplification represents a persistent weakness of climate models.

Fig
Fig. XX a) MMCO terrestrial temperatures (red diamonds) and SST (blue crosses) with methodological error plotted as the vertical bars and described in TableS1,S2.b) The spatial distribution of the terrestrial and SST proxy records used in the model data comparisons overlain onto the Miocene topography(Herold et al., 2008).
Fig. 1.(a) Longitudinal distribution of MMCO terrestrial temperatures (red diamonds) and SST (blue crosses) with proxy record error plotted as the vertical bars and described in TablesS1 and S2.(b) The spatial distribution of the terrestrial and SST proxy records used in the model-data comparisons overlain onto the Miocene topography(Herold et al., 2008).

Figure 1
Figure1.a) High AIS topography(Herold et al., 2011), c)  Pointwise MAT comparison between the MMCO simulation and proxy records (TableS1,S2).Vertical error bars are the modelled pointwise maximum and minimum temperatures from the extreme obliquity simulations (see Methods section methodological error is plotted as the horizontal error bars.The best fit line (black dashed) is weigh to include error uncertainty and is fitted across all points.The y-intercept and slope are reported in T 1.

Fig. 3 .
Fig. 3. (a) Modeled temperature anomaly for the MMCO560 (560 ppm CO 2 ) simulation minus the MMCO simulation ( • C).(b) Pointwise MMCO560 simulated global MAT compared against the proxy record MAT ( • C).(c) Modeled temperature anomaly for the MMCO800 (800 ppm CO 2 ) simulation minus the MMCO simulation ( • C).(d) Pointwise MMCO800 simulated global MAT compared against the proxy record MAT ( • C).These are the same terrestrial and SST records described in Fig. 1.Horizontal error bars indicate the uncertainty recorded by maximum and minimum temperatures of extreme orbital obliquity parameters (see Methods Sect.2.3).The best fit line (black dashed) is weighted to include error uncertainty and is fitted across all points and the y intercept and slope reported in Table2.

Figure 2
Figure 2. a) LOW AIS topography based on offline ice-sheet modeling (David Pollard, personal comms), b) modelled temperature anomaly between the LOW AIS simulation and the MMCO simulation with the high AIS.c) Pointwise MAT comparison between the LOW AIS simulation and proxy records (TableS1,S2).The best fit line (black dashed) is weighted to include error uncertainty is fitted across all points and the y-intercept and slope are reported in Table1.
Fig. 4. (a) High AIS topography (meters) used in Herold et al. (2011), (b) LOW AIS topography based on offline ice-sheet modeling (David Pollard, personal communications, 12 April 2012), (c) modeled temperature anomaly ( • C) between the LOW AIS simulation and the MMCO simulation with the high AIS.(d) Pointwise MAT comparison between the LOW AIS simulation and proxy records (TablesS1, S2).The best fit line (black dashed) is weighted to include error uncertainty and is fitted across all points and the y intercept and slope are reported in Table2.

Figure 3
Figure 3. a) Modelled temperature anomaly for the EP simulation minus the LOW AIS simulation (˚C), b)Pointwise EP case global mean MAT compared against the proxy record MAT (˚C).These are the same terrestrial and SST records and error bars described in Figure1.The best fit line (black dashed) is weighted to include error uncertainty is fitted across all points and the y-intercept and slope reported in Table1.
Fig. 5. (a) Modeled temperature anomaly for the EP simulation minus the LOW AIS simulation ( • C), (b) pointwise EP case global mean MAT compared against the proxy record MAT ( • C).These are the same terrestrial and SST records and error bars described in Fig. 1.The best fit line (black dashed) is weighted to include error uncertainty and is fitted across all points and the y intercept and slope reported in Table2.
Fig. 6.(a) Modeled temperature anomaly for the EP+ORBITAL25 simulation minus the LOW AIS simulation ( • C), (b) pointwise EP+ORBITAL25 case global mean MAT compared against the proxy record MAT ( • C).These are the same terrestrial and SST records described in Fig. 1.Horizontal error bars indicate the uncertainty recorded by maximum and minimum temperatures of extreme orbital obliquity as in Fig. 1.The best fit line (black dashed) is weighted to include error uncertainty fitted across all points and the y intercept and slope reported in Table2.

Table 1 .
Simulation names and boundary conditions for the PI and MMCO simulations.

Table 2 .
Compilation of model simulation names, imposed boundary condition changes, model and proxy MAT values, Equator to pole temperature gradient values, and model-data pointwise comparison statistics.The Equator to pole surface temperature gradient is calculated by averaging the mean annual temperatures over the absolute latitudes of 60-80 • minus 0-30 • ; 80 • is the maximum latitudinal extent of proxy records.b The slope and y intercept of the best fit line for the pointwise model and proxy comparisons in Figs.2-6.The best fit line is weighted to include the error uncertainty found in the proxy records (Tables a

523/2014/ Clim. Past, 10, 523-536, 2014 gure
DR2. a) Modelled temperature anomaly for the EP+ORBITAL25 simulation minus the LOW AIS mulation (˚C), b) Pointwise EP+ORBITAL25 case global mean MAT compared against the proxy cord MAT (˚C).These are the same terrestrial and SST records described in Figure1.Vertical error bars dicate the uncertainty recorded by maximum and minimum temperatures of extreme orbital obliquity rameters (22˚ and 25˚ respectively).The best fit line (black dashed) is weighted to include error certainty is fitted across all points and the y-intercept and slope reported in Table1.
• C higher climate sensitivity.The higher temperature sensitivity to CO 2 has been attributed to www.clim-past.net/10/ Henrot et al. (2010)and 19.2 • C at 720 ppm CO 2 .These simulations are roughly 4 and 2 • C colder compared to the MAT calculated from the proxy records presented here.This study also comes to similar conclusions about their model's inability to reproduce reconstructed warmth in the high latitude regions especially in the Southern Hemisphere.Hamon et al. (2012)also conducted fully coupled MMCO simulations under a variety of different changes in boundary conditions.Comparison to this study is difficult because results focused on regional temperature changes to AIS forcing and they did not report global MAT values.Henrot et al. (2010), using an intermediate complexity planet simulator explored changes in topography, seaways, CO 2 ,