Strongly eddying ocean simulations required to resolve Eocene model-data mismatch

Model simulations of past climates are increasingly found to compare well with proxy data at a global scale, but regional discrepancies remain. A persistent issue in modeling past greenhouse climates has been the temperature diﬀerence between equatorial and (sub-)polar regions, which is typically much larger in simulations than proxy data suggest. Particularly in the Eocene, multiple temperature proxies suggest extreme warmth in the southwest Paciﬁc Ocean, where model simulations consistently suggest temperate conditions. Here we present new global ocean model simulations at 0.1 ° horizontal resolution for the middle-late Eocene. The eddies in the high-resolution model aﬀect poleward heat transport and local time-mean ﬂow in critical regions compared to the non-eddying ﬂow in the standard low-resolution simulations. As a result, the high-resolution simulations produce higher surface temperatures near Antarctica and lower surface temperatures near the equator compared to the low-resolution simulations, leading to better correspondence with proxy reconstructions. Crucially, the high-resolution simulations are also much more consistent with biogeographic patterns in endemic-Antarctic and low-latitude-derived plankton, and thus resolve the long-standing discrepancy of warm subpolar ocean temperatures and isolating polar gyre circulation. The results imply that strongly eddying model simulations are required to reconcile discrepancies between regional proxy data and models, and demonstrate the importance of accurate regional paleobathymetry for proxy-model comparisons.


Introduction
Model-data comparisons for warm periods in the geological past can be used to test the performance of climate models under greenhouse conditions (Tierney et al., 2020; -2-manuscript submitted to AGU Advances Tabor et al., 2016;Hutchinson et al., 2021;Lunt et al., 2021;Kennedy-Asser et al., 2020;Braconnot et al., 2012;Zhu et al., 2020;Liu et al., 2009;Schmidt et al., 2014;Dowsett et al., 2013;Cramwinckel et al., 2018).Some fully coupled climate models using stateof-the-art Eocene geographic boundary conditions (Baatsen et al., 2016) and greenhouse gas forcing simulate climates that correspond well to reconstructions of tropical sea surface temperature (SST) and deep ocean temperature (Cramwinckel et al., 2018).However, in such simulations these models regionally simulate much cooler conditions in extratropical regions than proxy data suggest, particularly in the southwest Pacific (Lunt et al., 2021;Cramwinckel et al., 2018;Baatsen et al., 2020;Lunt et al., 2012;Huber & Caballero, 2011).Consequently, depending on the radiative forcing, models either produce SSTs near the equator that are higher than proxy data indicate or SSTs at midto-high latitudes that are much lower than proxy reconstructions, leading to stronger meridional SST gradients.
One challenge in paleoclimate model-data comparisons is the scale difference between proxies and models.The proxies capture a regional environment and effects of smallscale regional setting (e.g.geography, bathymetry, and oceanography), while general circulation models have difficulties capturing regional climate correctly due to the coarse resolution that is typically used (1 • horizontally or coarser for the ocean) (Nooteboom et al., 2020;Dowsett et al., 2013;Harrison et al., 2016;Eyring et al., 2019;Kennedy-Asser et al., 2020;Tabor et al., 2016).The quality of ocean models improves considerably at a higher horizontal resolution (0.1 • ) (Griffies et al., 2015;Dong et al., 2014;Sun et al., 2019;Müller et al., 2019;Dong et al., 2014;Viebahn et al., 2016;Hewitt et al., 2016;Mc-Clean et al., 2006), especially their regional flow (Delworth et al., 2012;Marzocchi et al., 2015;Nooteboom et al., 2020).This is not only due to higher level of detail, but also because of the smaller scale interactions resolved (including mesoscale eddies of 10-30km size) that influence the large-scale flow properties (Porta Mana & Zanna, 2014) and increase the importance of the local setting (i.e. the paleogeography and bathymetry) in the resulting regional ocean flow.
Biogeographic patterns of microplankton (e.g.dinoflagellate cysts; dinocysts) in Southern Ocean marine sediments have been used as tracer of past surface oceanography (Huber et al., 2004).For instance, Eocene sediments deposited near Antarctica contain dinocyst species that are endemic to circum-Antarctic locations (Bijl et al., 2011).
Hence, Southern Ocean regions with many of these endemic species, as opposed to those -3-manuscript submitted to AGU Advances with abundant cosmopolitan species, must be oceanographically connected.This implies that these biogeographic patterns of dinocysts provide a direct proxy of the flow direction itself (Bijl et al., 2011).So far, climate models were broadly able to match the circulation patterns deduced from microplankton endemism in the Southern Ocean, sometimes after adaptations of the model paleobathymetry (Houben et al., 2019;Bijl et al., 2013;Huber et al., 2004) or details of the configuration of critical Southern Ocean gateways (Sijp et al., 2016).However, these model simulations cannot explain the occurrence or absence of endemic dinocysts at some sites.In addition, state-of-the-art fully coupled climate model simulations did come close to the proxy-based warmth in the southwest Pacific Ocean, but this required a flow through the Tasmanian Gateway which was incompatible with microplankton-based evidence of surface ocean flow (Stickley et al., 2004;Cramwinckel et al., 2020).Consequently, no model simulation exists that can reconcile southwest Pacific Ocean warmth with ocean flow that is compatible with the plankton records (Baatsen et al., 2020).
Here we show that high resolution ocean model simulations partly solve this mismatch, using sinking Lagrangian particles to represent biogeographic patterns of microplankton in the ocean model simulations (Nooteboom et al., 2019;Huber et al., 2004).We present the first simulations of a global eddying Eocene ocean model with a 0.1 • horizontal resolution (HR2 and HR4; Table 1).These simulations are initialized and forced with atmospheric fields from an equilibrium state of a coarser (1 • ) resolution model with a fully coupled ocean and atmosphere (LR2 and LR4; Table 1) (Baatsen et al., 2020).Hence, the high-and low-resolution simulations have a similar atmospheric forcing and bathymetry.
The new high-resolution simulations are run for a few decades (42 and 27 years for HR2 and HR4 respectively), sufficient for upper-ocean circulation to equilibrate.

Effect of model resolution on Eocene flow
The resulting ocean circulation is different between the eddying and non-eddying configurations (Fig. 1).In the eddying simulations, the time-mean flow strength has a higher spatial variability, the bathymetry has a larger influence on the flow strength and direction (especially in the Southern Ocean; see Supporting Information Fig. S3 for the bathymetry), and local scale features are much more pronounced, compared to the lowresolution model.All western boundary and equatorial currents increase in strength, except in the North Atlantic.The spatial structure and separation locations of the western boundary currents are also shifted.For instance, the eastward Agulhas separation  The EAC flow provides an example of the stronger influence of the paleobathymetry on the flow in HR4 compared to LR4, even though the bathymetry is the same in both configurations.Eddies are responsible for the downward transfer of momentum input at the ocean surface by winds that is eventually balanced by bottom form stresses (Munday et al., 2015).As a consequence, the flow is strongly determined by isobaths (i.e.lines of constant bathymetry) (Rintoul, 2018;Marshall, 1994).Hence, the bathymetry has a much larger influence on the flow if the ocean is eddying (in HR4 and HR2) than if it is not (LR4).In HR4, the EAC is steered further southeastward than in LR4 along the submerged continental block of Lord Howe Rise (see Supporting Information Fig. S3 for the bathymetry).Moreover, jets like the EAC have a narrower structure in the eddying flow, due to interactions between eddies and the time-mean flow (Waterman et al., 2011), which has profound impacts on the regional oceanography.

Model-data comparison: plankton biogeography
The new Eocene ocean model velocity fields enable us to use sinking Lagrangian particles (Nooteboom et al., 2020) to reveal biogeographic provinces of endemic microplank- The sedimentary endemism of the data is the percentage of measured endemic species at the site (Bijl et al., 2011), representative of 41-39 Ma.
Labeled sites are named in the main text.
-7-manuscript submitted to AGU Advances ton in the Eocene Southern Ocean.In this way, we can test how representative the modeled flow is compared to the reconstructed ocean flow from sediment records.In this approach, it is determined where sedimentary particles originated from at the ocean surface, while taking into account how the particles were advected by ocean currents during their sinking journey.If these virtual particles originate from an environment with a temperature below a threshold value indicated by SST , the particle is assumed to originate close to Antarctica, and flagged as representing Antarctic-endemic dinocyst species (see Materials and Methods section and Supporting Information Fig. S4 for an illustration).
Due to the circulation differences between eddying and non-eddying simulations, the model-derived occurrence of Antarctic-endemic sedimentary dinocysts is clearly different between both configurations (Fig. 2).While the endemism is more strongly dependent on latitude and a sharper boundary exists between low-endemism and high-endemism in LR4, sinking particles are transported further away from Antarctica in specific areas (especially near western boundary currents) in HR4.As a consequence, the occurrence of several recorded endemic species can be explained in HR4, while it cannot in LR4 (see e.g.site SanB).Moreover, the modeled endemism in the non-eddying LR4 cannot match with both DSDP277 and MH at the same time, because these sites contain an opposite signal (i.e.MH contains endemic species and DSDP277 does not) while being located closely to each other.In HR4 on the other hand, the sedimentary particles in site DSDP277 (Fig. 2) originate only from the warm waters of the southeastward flowing EAC, while the closely located site MH also contains particles originating from cold waters in the east, in agreement with the occurrence of endemic species at MH.
Overall we find that only the eddying simulations produce circulation patterns consistent with plankton biogeographic patterns.As a result, the model-data comparison has a better overall fit in HR4 compared to LR4 (Fig. 2a,c).The model-data fit improvement in HR4 compared to LR4 highlights the need for accurate reconstructions of the geographic boundary conditions (Baatsen et al., 2016)   3).However, the transient response of the upper ocean equilibrates similarly in the 2× and 4×pre-industrial CO 2 cases in a few decades, which also results in a similar timemean surface flow (Supporting Information Fig. S2).This implies that plankton biogeographic patterns and surface ocean circulation are to a large extent affected by bathymetry, rather than the climate boundary conditions (e.g.atmospheric CO 2 ) of the model.
At the beginning of the HR2 and HR4 simulations, much of the energy input at the surface is used to set up the circulation and the development of eddies, as can be seen from a reduction of Southern Ocean gateway transports in the first 5 years (similar in both HR2 and HR4), after which they recover (Fig. 3c-f).After 9 years, the Drake Passage transport (through the gateway between South America and Antarctica) exceeds the transport in the low-resolution simulations and equilibrates at a higher level.The increased Drake Passage transport is mainly caused by the lower (more realistic) viscosity that the high-resolution models allows compared to the low-resolution model (which becomes numerically unstable at this low viscosity value).Interestingly, the volume transport through the Tasman Gateway in HR2 and HR4 does not exceed the volume transport in LR2 and LR4.Instead, a larger fraction of the water is transported north of Australia, resulting in the stronger southeastward East Australian Current (EAC) in the South Pacific (Fig. 1).

Model-data comparison: sea surface temperatures
Now that the high-resolution Parallel Ocean Program (POP) model simulates an Eocene ocean flow, which is consistent with proxy data for ocean circulation, we compare the results of these simulations to proxy data for SST.SST distributions however, are also influenced by the model background state and sensitive to their global-scale equilibration.Moreover, the background flow affects the distribution of heat differently in the eddying versus non-eddying simulations.Meso-scale eddies are important for the distribution of heat, and eddying ocean models do a better job in representing heat transport compared to non-eddying models that use parameterizations for eddy-induced heat transport (Viebahn et al., 2016;Griffies et al., 2015;Dong et al., 2014).
Indeed, heat is distributed differently in the top km of the eddying compared to the non-eddying simulations (Fig. 3a and 3b).Eddies efficiently transport heat to the subsurface (Delworth et al., 2012), which leads to subsurface warming in both eddying simulations (HR2 and HR4) and a lower vertical temperature gradient compared to LR2 -10-manuscript submitted to AGU Advances and LR4.However, in HR2 the surface cools more, while the subsurface warms less compared to HR4.
Much of the heat transport change from LR to HR is related to the Southern and Northern Meridional Overturning Circulation (SMOC and NMOC respectively).In both HR2 and HR4, North Pacific sinking develops (in a few decades) next to existing South Pacific sinking, while in the low-resolution simulations there is only Southern Hemisphere sinking (see Supporting Information Fig. S8).Overall, the North Pacific sinking leads to an increase in the NMOC and a decrease in the SMOC.These changes in the MOC are stronger in HR2 compared to HR4, and both the NMOC and SMOC are still increasing in magnitude at the end of the HR2 simulation.
The SMOC also differs in structure between the high-and low-resolution simulations (see the mixed layer depth in Supporting Information Figure S8).In HR2 and HR4, more volume transport through the Drake Passage increases the surface salinity in the South Atlantic resulting in denser surface water in the Weddell Sea (Tumoulin et al., 2020).
Therefore, the main deepwater formation location is the South Atlantic in HR2 and HR4, while it is the South Pacific in LR2 and LR4.
These results imply that HR2 and HR4 are run long enough for the upper-ocean circulation to equilibrate, while the deep ocean is not in equilibrium yet, as can be seen from the MOC in HR2 (Fig. 3g).Although the transient evolution of the deep ocean circulation differs between HR2 and HR4, we can nevertheless investigate their impact on SST distributions and compare those to proxy-data.
Both the tropical and Arctic Ocean cool significantly in HR2 compared to LR2, while in HR4 the equatorial regions cool less and high-latitude (north and south) regions warm more as compared to LR4 (see Fig. 4).For both atmospheric CO 2 levels, local SST differences between the high-and low-resolution simulations mostly occur near western boundary currents of which the location shifts in the eddying simulation (Fig. 4a and     d).These shifts have an effect on the model-data comparison at sites near western boundary currents.In fact, the EAC transports warm waters southeastwards in the southwest Pacific, which (partly) explains why sites in the southwest Pacific are found to be warmer compared to model simulations with a coarse resolution.Notably, similar SST changes occur near the Kuroshio and Agulhas currents.The Weddell Sea warms up in HR2 and HR4 compared to LR2 and LR4 respectively, which is related to the South Atlantic sinking that occurs in HR2 and HR4.
-11-manuscript submitted to AGU Advances  Hinsbergen et al., 2015), each site is compared to the model SST value from up to 3 • distance of the site that minimizes the RMSE of the scatter plot (similar to (Baatsen et al., 2020); see Supporting Information Fig. S7 for a point-to-point comparison).The dashed black line is the one-to-one line representing the perfect match between model and proxy data.
-12-manuscript submitted to AGU Advances Climate models generally do not produce the low meridional temperature gradients of warm climates as inferred from proxy data (Huber & Caballero, 2011;Sijp et al., 2014).While the simulations LR2 and LR4 were found to generate a lower meridional SST gradient compared to other models of 1 • horizontal resolution or coarser (Baatsen et al., 2020), this gradient reduces further in HR2 and HR4.The tropics are cooler in

Conclusion and outlook
We have shown that an eddying Eocene ocean simulation provides a more detailed ocean flow compared to a non-eddying version of the same model.As a result, modeldata mismatches in the geologic past (Lunt et al., 2021;Hutchinson et al., 2021;Baatsen et al., 2020;Houben et al., 2019;Bijl et al., 2011;Huber et al., 2004) can at least partly be explained by the lack of eddies in the ocean models used.Our eddying simulations of the late Eocene are better able to explain the occurrence or absence of endemic dinocyst species near Antarctica compared to non-eddying simulations.The SST model-data comparison also improved in the eddying compared to non-eddying simulations.
The explicit representation of eddies in ocean models may have implications for comparison of models with other proxy types than considered here.For instance, pollen-based temperature reconstructions imply that it did not freeze at the Antarctic coast during winter in the early Eocene (globally ∼ 6 • C warmer than the late Eocene), despite polar darkness (Pross et al., 2012).Eddy-induced flow, and its impact on ocean heat transport, could in part explain such conditions.
The simulations in this paper are computationally expensive.However, other types of model set-ups may be interesting if not limited by computational capabilities.First, the strong influence of bathymetry on the eddying flow implies that the uncertainty of paleogeography reconstructions has a major impact on model-data comparisons.Future studies could make adaptations on the bathymetry within uncertainty of paleogeographic reconstructions, to find its impact on the modeled ocean circulation and model-data comparison.Moreover, since the eddying flow has a direct response to bottom topography, it seems suitable for a downscaling, or eddy parameterization type of approach to obtain this influence of bathymetry on the flow with reduced computational costs.However, these type of approaches are found to be challenging in present-day configurations (Fox-Kemper et al., 2019;Nooteboom et al., 2020;Lanzante et al., 2018).
Second, we used the model equilibrium of the non-eddying climate model simulations (which are in radiative equilibrium (Baatsen et al., 2020)) to start and force the eddying model.However, this switch induces a drift of the deep ocean circulation, which is not equilibrated yet in the high-resolution simulations of this paper.Hence, the background state of the model will change further if the model is run for longer time periods (a few millennia).Future simulations may have the capabilities to perform longer simulations.These changes of the model background state on long time scales might have implications for the regional flow and the quality of the model-data comparisons.
Finally, atmospheric feedbacks greatly influence the ocean model background state on long time scales, such as the meridional overturning circulation (den Toom et al., 2012;Rahmstorf & Willebrand, 1995;Arzel et al., 2011;Zhang et al., 2010).Hence, the highresolution ocean should be coupled to a high-resolution atmosphere, which could further enhance the meridional transport of heat and lead to an improved model-data comparison.
4 Materials and Methods

Data
We used two datasets in this paper.The first includes the SST proxies from U k 37 , TEX H 86 , Mg/Ca, ∆ 47 and δ 18 O, which are described in detail in (Baatsen et al., 2020).
Proxy-based SST reconstructions come with uncertainties, limitations and biases (Hollis et al., 2019), related to the depth, or season they represent.The second dataset are sediment samples with dinocysts from (Bijl et al., 2011), combined with the samples described in (Houben et al., 2019;Cramwinckel et al., 2020;Bijl et al., 2021).We averaged dinocyst abundance of Endemic-Antarctic, cosmopolitan and low-latitude-derived for the respective time slices.

Model set-up
We used the Parallel Ocean Program (POP) (Viebahn et al., 2016;den Toom et al., 2014;Smith et al., 2010) to perform eddying ocean model simulations for the middlelate Eocene (38Ma).To derive the forcing of this model, we made use of the fully-coupled (ocean and atmosphere) simulations with Community Earth System Model v1.0.5 (CESM) from (Baatsen et al., 2020), with a non-eddying ocean.We used both CESM simulations with 2× pre-industrial atmospheric CO 2 (LR2) and 4× pre-industrial CO 2 (LR4) configuration.
The high-resolution POP is forced at the surface by a fixed atmosphere of the CESM simulation.To construct the surface forcing, we interpolated the average (over the last 50 model years of LR2 and LR4) Sea Surface Temperature (SST), Sea Surface Salinity (SSS) and wind stress (zonal and meridional) of the CESM simulation for every month of the year (such that a seasonal cycle is included in the surface forcing).These SST and SSS fields were used as restoring boundary conditions at the surface.The restoring boundary conditions imply that POP is 'pushed' towards the SST and SSS output of the CESM at the surface with a specific timescale (30 and 10 20 days respectively).This means that differences between the SST and SSS at different model resolutions arise due to the internal transport (vertical and horizontal) of heat and salt in the ocean, not due to the surface forcing.The bathymetry that CESM uses was interpolated linearly on the highresolution grid that POP uses, making both bathymetries similar (see the code at https:// github.com/pdnooteboom/MCEocene).
For initialisation of the eddying model, the three-dimensional ocean output at the end of the CESM simulations (LR2 and LR4) is interpolated on the higher resolution grid that the POP (HR2 and HR4) uses.We simulated 42 and 27 years in total for HR2 and HR4, respectively.Since we investigate the response of the simulations to an increase in horizontal resolution, the same 5 model years of both HR2 and HR4 are used in most analyses in this paper: year 23 to 27.For the same analyses of the low-resolution simulations (LR2 and LR4), we used the last 5 years of these simulations.
Using this setup of POP, we can investigate the sensitivity of simulations to the studied resolution difference only, because the model is forced by the same atmosphere and their geographic boundary conditions are based on the same reconstruction of (Baatsen et al., 2016), and the three-dimensional eddying ocean is initialized by the equilibrated output of the CESM.As a result, the atmosphere is representative of the middle-late Eocene climate, but does not respond to changes in the ocean.We hence cannot investigate the effect of atmospheric feedbacks on the results (den Toom et al., 2012;Rahmstorf & Willebrand, 1995;Arzel et al., 2011;Zhang et al., 2010).
The model set-up is suited to study the effects of model resolution on Eocene ocean flows, but it is not suitable to study dynamics which involve atmospheric coupling, such as the El Niño Southern Oscillation.The model set-up can best be used to investigate the upper ocean circulation, as the deep ocean is not in equilibrium yet.Therefore, we can only use this setup to obtain a transient response of the deep meridional overturning, not its equilibrium.

Sinking Lagrangian particles
To quantify sedimentary dinocyst endemism in the model, we applied a similar backtracking analysis of virtual sinking Lagrangian particles as in (Nooteboom et al., 2019) (Supporting Information Fig. S4).This implies that we released these particles at the ocean bottom and tracked them back in time while sinking and being advected by the three-dimensional flow from POP, until they reached 10m depth.We released particles on a 2 • ×1 • grid of locations between 32-80 • S every day for a year and waited until all of the particles reached the near-surface (i.e.17,520 particles in total).This analysis requires a higher than monthly temporal resolution of model output (Nooteboom et al., 2020;Qin et al., 2014).Therefore we used daily fields for the years 35-42 (HR2) and years 20-27 (HR4) to perform this backtracking analysis.
-16-manuscript submitted to AGU Advances The used particle sinking speed of the Lagrangian particles in this paper is 6 m day −1 .
This represents a low sinking speed for single dinocysts (Anderson et al., 1985).We choose this low sinking speed, because it is considered as a lower bound of the realistic sinking speeds where most lateral transport occurs, which makes it easier to explain low abundances of dinocyst species.However, this sinking speed could in reality be different due to e.g.aggregation with other particles.We also applied a sinking speed of 25 m day −1 (see Supporting Information Figure S6), which represent small aggregates (Nooteboom et al., 2019).The main conclusions on the model-data comparison do not change if 25 m day −1 instead of 6 m day −1 sinking speed is used.
The percentage of dinocyst endemism in the model is determined by the percentage of particles that originated from an environment with a temperature below SST (which must be close to Antarctica; see Supporting Information Figure S4; similar approach as in (Huber et al., 2004)).The percentage of modeled dinocyst endemism is not expected to compare well with the percentage of measured endemic dinocyst, because this match is also sensitive to the species-specific susceptibility of dissolution during the sinking journey and their productivity at the ocean surface (Nooteboom et al., 2019).Therefore, we compare whether any endemic species occur in sites (0% or ¿0%) between model and data instead of the exact percentage.
We assume that the sinking Lagrangian particles are not greatly influenced by the fact that the deep circulation is not in full equilibrium yet in the eddying simulations.
Most of the lateral particle displacement occurs near the surface which is in equilibrium and where the currents are the strongest.Moreover, the eddying simulations are initialised with output from the non-eddying simulations, which are in reasonable equilibrium.The mechanistic development of the flow, given the heat and salt distribution from the initialisation, occurs in a few years (see also Fig. 3c-h).Hereafter, the flow changes slowly and may equilibrate after ∼1000 years due to the flow response to changing density distributions.The assumption that sinking Lagrangian particles are not greatly affected by the deep ocean equilibration, is supported by the results that use sinking Lagrangian particles in HR2 and HR4: These results are similar, even though the deep ocean circulation is different in HR2 and HR4.
to optimize model-data matches as in Fig. 2a,b: It are the details in the ocean flow that induce a better model-data fit in HR4 compared to LR4.The modeled dinocyst endemisms in the 2× and 4×pre-industrial atmospheric CO 2 configurations are similar (see Supporting Information Fig S1), even though HR2 and HR4 are forced by a different atmosphere and respond differently after initialisation (Fig.

Figure 4 .
Figure 4. Model-proxy data comparison: sea surface temperature (SST).The 2× and 4×preindustrial case are compared to SST proxy data of 38-34Ma and 42-38Ma respectively.(a), (d) SST difference of the high-compared to the low-resolution model with the site locations of the SST proxies for 2× and 4×pre-industrial carbon configuration, respectively and (b), (e) their zonal mean.(c), (f) the zonally averaged annual mean SST in the high resolution (black) and the low resolution (red) model for 2× and 4×pre-industrial carbon configuration respectively.The shaded areas show zonal spread (i.e.minimum and maximum) of the annual mean SST.Markers indicate SST proxy estimates with their uncertainty.(g-j) Scatter plots between proxy-derived and model-derived SST for all four configurations, with Root Mean Squared Errors (RMSE).Error bars represent proxy calibration errors.To consider the paleolocation uncertainty of sites (van HR2 and HR4 compared to LR2 and LR4, while in the zonal-mean the southern highlatitudes are only slightly warmer in HR4 (Fig.4d-f).Regionally, there is, however, significant warming of Southern Ocean SSTs in HR4.Overall, this improves consistency between the high-resolution model results and SST proxies in the tropics, while the modeled high-latitude SST values are often still lower than the proxy-derived SST values.The eddying simulations show stronger horizontal gradients in the time-mean SST field compared to the non-eddying simulations, which results in a higher time-mean SST variation in the model around the sediment sample sites.The model-data fit greatly improves in the eddying compared to non-eddying simulations (Fig.4g-j), although a mismatch with some sites remains (especially for the 2×pre-industrial CO 2 case) and the high-latitude temperatures are overall lower compared to the proxy data.Overall, the eddying ocean model improves the SST model-data match from the non-eddying model, because it alters the local transport of heat.However, the SST modeldata comparison is also sensitive to the model background state (i.e. the state of the ocean at a global scale), which depends on the used atmospheric forcing, paleogeography and long time scales phenomena, such as the deep meridional overturning circulation.Hence, the SST model-data mismatch could be reduced even further if better model boundary conditions are used which lead to a more realistic background state of the late Eocene.

Table 1 .
The ocean model simulations of the middle-late Eocene (38Ma) in this paper