Glacial Maximum world ocean simulations at eddy-permitting and coarse resolutions : do eddies contribute to a better consistency between models and palaeoproxies ?

Most state-of-the-art climate models include a coarsely resolved oceanic component, which hardly captures detailed dynamics, whereas eddy-permitting and eddyresolving simulations are developed to reproduce the observed ocean. In this study, an eddy-permitting and a coarse resolution numerical experiment are conducted to simulate the global ocean state for the period of the Last Glacial Maximum (LGM,∼ 26 500 to 19 000 yr ago) and to investigate the improvements due to taking into account the smaller spatial scales. The ocean state from each simulation is confronted with a data set from the Multiproxy Approach for the Reconstruction of the Glacial Ocean (MARGO) sea surface temperatures (SSTs), some reconstructions of the palaeocirculations and a number of sea-ice reconstructions. The western boundary currents and the Southern Ocean dynamics are better resolved in the high-resolution experiment than in the coarse simulation, but, although these more detailed SST structures yield a locally improved consistency between model predictions and proxies, they do not contribute significantly to the global statistical score. The SSTs in the tropical coastal upwelling zones are also not significantly improved by the eddy-permitting regime. The models perform in the mid-latitudes but as in the majority of the Paleoclimate Modelling Intercomparison Project simulations, the modelled sea-ice conditions are inconsistent with the palaeoreconstructions. The effects of observation locations on the comparison between observed and simulated SST suggest that more sediment cores may be required to draw reliable conclusions about the improvements introduced by the high resolution model for reproducing the global SSTs. One has to be careful with the interpretation of the deep ocean state which has not reached statistical equilibrium in our simulations. However, the results indicate that the meridional overturning circulations are different between the two regimes, suggesting that the model parametrizations might also play a key role for simulating past climate states.


Introduction
The Last Glacial Maximum (LGM) was a cold-climate event with a duration of around 6500 yr, centred approximately 23 000 yr ago (Clark et al., 2009), viz. at the end of the last glacial cycle and before the present warm phase.It is described, from palaeo-climate records, as the most recent maximum ice sheet extent over the continents, especially in the Northern Hemisphere with the large Laurentide and Fennoscandian ice caps over the Northern American and the Northern European continents, respectively (Peltier, 1994;Clark and Mix, 2002;Peltier, 2004).As a result of these large cryospheric changes, the sea level was around 120 m lower than today, exposing the continental shelves to the atmosphere and hereby modifying the present-day world ocean basins.The combination of an altered bathymetry, a changed hydrosphere and an atmosphere with lower greenhouse gas concentrations during this glacial phase may have led to modifications of the ocean state, for example, the temperature and salinity distributions, the tidal mixing and dissipation (Green et al., 2009), the transports of heat, mass and sediments (Seidov and Haupt, 1997) as well as the meridional overturning circulation (MOC).Consequently, this ocean state may have generated feedbacks to the global climate.The LGM hence constitutes a uniquely fascinating time slice of the earth's climate history, which can be used for understanding climate change, for testing general circulation models under different boundary conditions, and for reconstructing past scenarios on the basis of comparisons with the palaeo-record (Mix, 2001).
Reconstructions of the earth's climate variations are based on analyzing geological and biological samples (e.g.ice-and sediment cores, pollen, corals, tree rings or speleothems) and utilizing models.Palaeo-proxy data were first used to define and prescribe boundary conditions for atmospheric general circulation models (AGCMs) (Gates, 1976;Richard Toracinta et al., 2004), for AGCMs coupled with mixedlayer ocean models (Broccoli, 2000;Hewitt et al., 2003), and also for high-resolution atmospheric models (Kim et al., 2007).Subsequently it has proved possible to simulate climate variations with fully coupled ocean-atmosphere models (Braconnot et al., 2007a, b).The Paleoclimate Modelling Intercomparison Project (PMIP) was initiated in the early 1990s to evaluate and compare the response of numerical climate models under palaeoclimate conditions.Due to computational limitations, the evaluations were undertaken using coarse-resolution models.These simulate the large-scale structures of the ocean, but usually parameterize the subgridscale physics (unresolved structures) such as turbulence (see e.g.Gent and McWilliams, 1990).In eddy-permitting ocean models, the spatial resolution has been increased and the amount of subgrid-scale parametrization has been reduced.It has been demonstrated that for the present-day climate, eddy-permitting oceanic simulations improve the quality of the representations of the western boundary currents as well as those of the sea-ice conditions and the meridional heat transports in the North Atlantic and Southern oceans (FRAM Group, 1991;Tréguier et al., 2005;Hallberg and Gnanadesikan, 2006;Spence, 2010).Until now, this type of simulation has only been conducted over regional scales for the LGM climate (e.g.Yang et al., 2006;Mikolajewicz, 2011).
Since these high-resolution simulations are more realistic (due to small diffusive coefficients in the model, better transport of heat and salt in narrow passages or currents), they will become more and more important for testing the plausibility of specific past oceanic scenarios (Beal et al., 2011;Condron and Winsor, 2011) and for comparisons with palaeo-reconstructions (Otto-Bliesner et al., 2009).As pointed out by Hargreaves et al. (2011), proxy-data pertain to discrete source points and are valid over spatial scales which are smaller than that of the coarse-resolution model grid.Comparisons between the coarse-resolution climate simulations and reconstructed LGM sea-surface state have been performed.On one hand, they indicated that the ensemble of models designed under the PMIP can be regarded as globally reliable with respect to the Multiproxy Approach for the Reconstruction of the Glacial Ocean surface (MARGO) sea surface temperature (SST) data synthesis by Waelbroeck et al. (2009) (see, Hargreaves et al., 2011Hargreaves et al., , 2012)).On the other hand, it has been reported that although these models reproduce the strong SST meridional gradients and the cooling in the North Atlantic, they sometimes do not place the gradients at the right location or fail in estimating the magnitude of the regional cooling (Kageyama et al., 2006;Otto-Bliesner et al., 2009;Braconnot et al., 2012).The realism of model simulations of the LGM conducted during the PMIP2 project has also been discussed and summarized in the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (Jansen et al., 2007).There it is concluded that even if "the PMIP2 LGM simulations confirm that current AOGCMs are able to simulate the broad-scale spatial patterns of regional climate change recorded by palaeo data in response to the radiative forcing and continental ice sheets of the LGM", still, "Regional variations in simulated tropical cooling are much smaller than indicated by MARGO data, partly related to models at current resolutions being unable to simulate the intensity of coastal upwelling and eastern boundary currents".
The aim of the present investigation is to evaluate whether a LGM eddy-permitting oceanic simulation improves the results with regard to coarse-resolution models and palaeoproxy reconstructions, and to understand how a global ocean general circulation model (OGCM) behaves subject to glacial forcing and at the expected resolution of the next generation of global climate models.To achieve this, an eddy-permitting and a coarse resolution ocean simulations of the LGM period have been conducted.By applying a statistical analysis to the data sets, we estimate the accuracy of each type of simulation in representing the reconstructed surface state.A description of the model and boundary conditions for the two model regimes is provided hereafter.The LGM results are then analyzed and discussed with a focus on how the models behave compared with the available PMIP results, and how close the models are to the reconstructions of the LGM surface state.

The ocean model
The OGCM NEMO (Nucleus for European Modelling of the Ocean) (Madec, 2008) is used for designing the experiment.This ocean model is based on the primitive equations under the spherical earth approximation, the thin-shell approximation, the turbulent closure hypothesis, the Boussinesq, the hydrostatic, and the incompressibility approximations.Various horizontal grid mesh resolutions and parametrizations are available for this ocean model.A resolution can be directly referred by its grid's name called ORCA.Thus, ORCA1 and ORCA025, the two configurations used in the study, correspond to approximately the 1 • horizontal resolution and the 0.25 • resolution, respectively.The vertical resolution is based on depth-coordinates levels.The ORCA1 configuration is set up with 64 vertical levels whereas the ORCA025 contains 46 levels.The vertical mesh is refined near the surface and the partial step method (i.e.adaptive bottom boxes) is selected for a better representation of the bathymetry (Barnier et al., 2006).The bottom boundary conditions are modelled by a non-linear formulation of the bottom friction.The ocean tracers, such as the temperature and the salinity, are linked to the density via a non-linear equation of state (Jackett and McDougall, 1995) and the sub-grid scale physics in the coarse resolution ORCA1 configuration is based on the Gent and McWilliams parametrization (Gent and McWilliams, 1990).Thus, the zonal, meridional and vertical velocities are made of one relative velocity plus an eddyinduced velocity.In the ORCA1 (ORCA025) configuration, the lateral diffusivity is parameterized by an iso-neutral (geopotential) Laplacian operator with an eddy-diffusivity coefficient of 1000 m 2 s −1 (300 m 2 s −1 for ORCA025).Momentum and tracers are mixed vertically using the turbulent kinetic energy scheme of Gaspar et al. (1990).The ocean model is coupled every two model hours with the multilayer thermodynamic-dynamic LIM sea-ice model version 2 (Fichefet and Maqueda, 1997).The sea-ice model resolves the thermodynamic growth and decay of the ice, the ice dynamic and transport.The sea-ice is considered as a 2-D viscous-plastic body and the model takes into account the sub-grid scale effect of snow and ice thickness.The model is integrated for a period of 150 yr by periodically repeating the surface forcing.Note that the runs do not reach a statistical equilibrium and firm conclusions can only be drawn for near-surface quantities.

The boundary conditions
For modelling the LGM climate, a reconstructed topography is required.The ICE-5G reconstruction (Peltier, 2004), that we used, includes bathymetry, altimetry and ice sheet reconstructions.The latter are based on geological insights as well as a sea-level model.The geomorphology of the continental plates during the LGM is similar to those of the present day.The major difference is the emergence of continental shelves due to a sea level approximately 120 m lower than today.In the ocean simulations, ice sheets and closed basins are considered as land points.
Different techniques can be used to initialize the ocean model in the LGM simulations.The PMIP2 protocol recommends starting the integrations either by using a spin-up procedure or from a previous LGM state generated by other simulations.The former procedure is based on integrations made from pre-industrial initial conditions and glacial boundary conditions, this in order to reach a cold equilibrium.Because the integration time is too long when starting from the recent-past ocean state, we have chosen the cold-state initialization technique and an ocean at rest.The temperature and salinity fields are interpolated onto the ocean grid mesh from a quasi-equilibrated integration carried out with the Community Climate System Model version 3 (CCSM3) (Brandefelt and Otto-Bliesner, 2009).The ocean state in this quasiequilibrated climate model integration differs substantially from the previous integration with CCSM3 used in the PMIP investigations (Brandefelt and Otto-Bliesner, 2009).In this equilibrated state, the annual mean surface temperatures are colder than previously reported in the Nordic Seas, north Atlantic Ocean and northern-most Pacific Ocean, mainly associated with the increase in Northern Hemisphere sea-ice extent, the reduction in the strength of the Atlantic overturning and the northward heat transport.The LGM global-ocean averaged salinity and temperature in this quasi-equilibrium is 36.59PSU and 0.60 • C.
The surface boundary conditions between the ocean, the sea ice, and the atmosphere are determined using the NCAR bulk formulae (Large and Yeager, 2004).This is the most popular method and has been used as the reference-surfacefluxes computational method for the numerical-model evaluations in, for example, the Drakkar experiments (Barnier et al., 2007;Brodeau, 2007;Brodeau et al., 2010).The LGM surface atmospheric variables originate from the CCSM3 model integration (Brandefelt and Otto-Bliesner, 2009), previously mentioned.The horizontal resolution of its atmospheric component is 128 longitudinal by 64 latitudinal points (T42).The horizontal resolution of the ocean component is approximately 1 • .The forcing is based on a 49 yr data set from the quasi-equilibrium LGM2 period 1412-1460 in Brandefelt and Otto-Bliesner (2009).The model is integrated for 150 yr by repeating the atmospheric forcing three times and no restoring term in temperature and salinity is applied at the sea surface since the main goal of our experiment is to investigate the impact of the ocean grid resolution on the representation of the surface state.Consequently, the salinity and temperature feedbacks on the atmosphere are not modelled.

Statistical analysis
The simulated surface states from the models are compared with the MARGO data set (Waelbroeck et al., 2009).This is a compilation of almost 700 sediment samples located especially in the North Atlantic, the Southern Ocean and the tropical regions.From these data, reconstructions of the annual mean (hereafter ANN), the boreal winter (January-February-March, JFM) and the boreal summer (July-August-September, JAS) SSTs have been produced.
The annual-mean, the boreal winter and boreal summer SSTs are examined for each sediment core location.The comparison is made by taking the results at the model gridbox coordinates closest to the location of the proxy datapoint.The performance of the models is evaluated quantitatively by their skill score S and illustrated using Taylor diagrams (Taylor, 2001).The Taylor diagram includes the correlation coefficient, the standard deviations and the "centred" root mean square error (RMSE) between the two fields.The skill score S (Taylor, 2001) is a measure of the correlation between the simulated and reconstructed SSTs, given their respective variability: where R is the Pearson, product-moment correlation coefficient between simulated and reconstructed SSTs; σf the ratio between the model variance and the MARGO data set variance; and R 0 = 1 is the maximum (positive) correlation attainable.In our cases, the reference fields are the MARGO annual-mean, boreal winter (JFM) and boreal summer (JAS) SSTs, while the "test" fields are the model outputs.The skill scores are defined with R 0 set equal to 1 (i.e. when the model results exactly fit the reconstruction).A skill score S = 1 means that the model performs well while a skill score S close to zero is a bad score.The Taylor diagrams are evaluated over four latitudinal bands (50-25 • S, 25 • S-25 • N, 25-50 • N and 50-90 • N) to isolate the regional changes due to the introduction of the permitted eddies.
In order to investigate the local impact introduced by the high-resolution simulation, statistical analyses are applied on the regional scale in the Agulhas current (30-50 • E; 50-10 • S), the Gulf Stream (70-20 • W; 30-45 • N), the Kuroshio (120-160 • E; 20-35 • N) and regions in the Southern Ocean.The mesoscale eddies are relatively important in these regions and thus eddy-permitting and non-eddypermitting experiments could show particularly large differences in the SSTs.For each region, linear regression tests between MARGO and model's SSTs are performed and the corresponding p value, correlation coefficient, slope (between SST MODEL vs. SST MARGO ) and intercept of the regression line are given.We test the null hypothesis that the slope is equal to zero and consider that the model and proxy data are significantly correlated, if the p value of the test is less than 0.05.

Results and discussions
In this section, the model behaviours under glacial forcing are firstly described and discussed with regards to other palaeo-simulations.Then, the SSTs and the sea-ice cover and the meridional circulations from each simulation are compared with the palaeo-reconstructions.The analyses cover the last 50 yr of the model integrations.

Sea surface temperature and salinity
The simulated ORCA1 and ORCA025 time-averaged SST distributions are found to be almost symmetric around the equator with the strongest meridional gradients in the midlatitude regions (Fig. 1a).In comparison with the modern state, the North Atlantic, the North Pacific, the Arctic and the Antarctic show a tendency toward cold surface temperatures (less than 1 • C) due to the sea-ice cover.The highest SSTs are found in the equatorial region of the western Pacific and the eastern Atlantic, and the eastern and equatorial Pacific cold tongue is also captured by the models.These latter features are consistent with those found in the PMIP2 simulations (Otto-Bliesner et al., 2009).In contrast, the low SSTs simulated in the North Atlantic differ and are caused by the cold atmospheric conditions extracted from the numerical experiment by Brandefelt and Otto-Bliesner (2009).The simulated LGM annual mean Sea Surface Salinity (SSS) shows that the most saline surface waters are found in the Southern Ocean (brought about by brine rejection), the south tropical band and with maxima in the tropical Atlantic and Mediterranean Sea (due to evaporation) (Fig. 2a).Two surface waters face each other in the North Atlantic: the warm, saline tropical water and the cold, fresh mid-latitude water created by sea-ice melting.Fresher surface waters are also found at river mouths (e.g. the Congo and Mackenzie Rivers) and over the littoral where the ice-sheets melt.The SSS difference in the Sea of Japan originates from a larger inflow of North Pacific saline water via the Tsugaru Strait simulated in the coarse resolution experiment.
In order to identify the regions where the mesoscale eddies play a role in the LGM thermohaline surface state, the maps of the time-averaged "ORCA1 minus ORCA025" SSTs and SSSs are represented in Figs.1b and 2b.The permitted mesoscale eddies mainly modify the distribution of temperature and salinity in the 30 to 60 • latitudinal bands, i.e. in the regions where the strongest meridional gradient and the strongest surface currents occur, such as the Gulf Stream, the Kuroshio, the Antarctic Circumpolar Current (ACC), and the Agulhas current.The Gulf Stream is about 1 • C warmer in the eddy-permitting simulation than in the coarse resolution experiment, whereas the central North Atlantic basin is 1 • C colder and 1 PSU fresher due to about 10 % larger formation of sea-ice in the ORCA025 experiment (see Supplement).The Agulhas current and Agulhas retroflection are warmer in the ORCA1 simulation probably due to larger diffusive coefficients.The ORCA025 simulation is about 1 PSU saltier along the 45 • S circle and the Kuroshio extension transports more saline (1 PSU) and warm (3 • C warmer) waters than in the ORCA1 experiment, associated with the presence of mesoscale eddies and reduced sea-ice cover (Supplement).The surface waters in the Arctic basin are fresher in the ORCA025 simulation than in the coarse resolution simulation, especially at the river mouth and might be the results of the different mixing parametrizations between the two regimes.Finally, in the inter-tropical region, the SST differences are weak suggesting that the eddy permitting may not improve the regional cooling near the upwelling regions.

Sea ice cover
The large sea-ice cover during the LGM can modify the distribution of the surface temperatures and salinities.During the boreal winter season, the models simulate sea-ice cover in the Nordic Seas, the northwestern North Atlantic, the Labrador Sea, and the North Atlantic sea-ice extent reaches almost 40 • N (Fig. 3a).The central Arctic basin and western Fram Strait have perennial sea ice.In the Southern Ocean, the austral summer sea-ice edge is located near 50 • S between 60 • W and 120 • E and near 55 • S in the rest of the circumpolar region (Fig. 3b).For the boreal summer season, the LGM sea-ice fraction is reduced in the Labrador Sea, the central North Atlantic and the Norwegian Sea, suggesting possibly ice-free conditions during the glacial period (Fig. 4a).In the Southern Ocean, the sea-ice extent during austral winter is relatively similar to the austral summer, showing a weak seasonality (Fig. 4b).The simulated sea-ice cover is however consistent with the CCSM3 simulation carried out by Brandefelt and Otto-Bliesner (2009), which showed an increased sea-ice fraction (especially in the Northern Hemisphere) at their quasi-equilibrated stage.Increasing trends of sea-ice fraction have been observed in almost all PMIP2 models (Braconnot et al., 2007b), but their magnitudes differed.It has also been reported that the CCSM3 model simulates an approximately 50 % larger sea-ice area in the Southern Ocean than the other PMIP models (Murakami et al., 2008).Since both simulations are forced by an atmospheric forcing state originating from a CCSM3 integration and that sea ice is strongly controlled by the atmospheric state, the sea-ice covers are also particularly large in our simulations.
The comparison of the sea-ice seasonality between the eddy-permitting and the coarse experiments shows that the sea-ice area is weakly affected by the model resolution in the Northern Hemisphere, whereas in the Southern Ocean the area is about 2 × 10 6 km 2 smaller in an eddy-permitting simulation (Fig. 5a and b).This difference is the results of the warmer SSTs in the ORCA025 simulation due to presence of eddies in the ACC and larger meridional heat transport discussed hereafter (Fig. 8).Larger differences between the two model integrations appear for the sea-ice volume seasonality (Fig. 5c and d).In the Northern Hemisphere, the seaice volume is about 20 × 10 3 km 3 larger in the ORCA025 simulation.It mainly occurs in the Arctic Ocean and might come from the different parametrizations between ORCA1 and ORCA025 and from the fact that the Arctic Ocean is  almost closed, constraining a weak export of sea ice and a large thickness growth.The Southern Ocean sea-ice volume is about 10 × 10 3 km 3 larger in the ORCA1 simulation than in the ORCA025 experiment, as a result of a larger southward heat transport in the eddy-permitting experiment.

Deep-water formation, volume and heat transports
The differences in the surface conditions can alter the formation of the deep waters.The mixed layer depth (MLD) is an indicator for the thermocline ventilation and the deep-water formation.It can be calculated as the maximum depth h where the potential density difference ρ between the surface and h is smaller than 0.01 kg m −3 .In the present-day climate, the maximum MLDs are found in the North Atlantic (in the Nordic Seas and the Labrador Sea) and are associated with the formation of North Atlantic Deep Water (NADW) (de Boyer et al., 2004).In the Southern Hemisphere, the deep mixed layers are located in the southeast Pacific ACC region where Subantarctic Mode Water and Antarctic Intermediate Water are formed.In the LGM eddy-permitting simulations, the deep mixed layers are found in the region where thick sea-ice is formed and where the warm and saline tropical waters encounter the cold and fresh conditions near the sea-ice edge.In the Northern Hemisphere, deep-water formation hence takes place in the Arctic Ocean and near 30 • N (Fig. 6a), suggesting that the LGM Atlantic thermohaline circulation may differ from that of the present-day conditions.In the Southern Hemisphere, this process takes place in the northern branch of the ACC and along the coast of Antarctica (adjacent to the Weddell and Ross Seas), also here associated with the sea-ice dynamics (Fig. 6b).These changes are also diagnosed in the equilibrated simulation by Brandefelt and Otto-Bliesner (2009) with a MLD that is reduced (increased) in the North Atlantic (in the Weddell Sea) and the convection sites shifted equatorward.In comparison with the ORCA025 experiment, the mixed layer depths are slightly deeper in the ORCA1 simulation probably due to a larger vertical diffusion coefficient, and in the Arctic Ocean because the basin is relatively homogenous in temperature and salinity (see Supplement).The locations of these convection sites control the overturning circulations.The Atlantic meridional overturning circulation (AMOC) is characterized by a northward surfacewater transport, a sinking of dense waters at high latitudes and a deep return flow towards the Southern Ocean (Fig. 7a  and b).The sinking of dense water in the North Atlantic takes place between 30 and 65 • N, where the maximum deep mixed layers are found.The return flow is found at a depth of ∼ 1500 m in the ORCA025 simulation and 2000 m in the ORCA1 simulation, which is in line with the larger MLD diagnosed for the ORCA1 simulation.Moreover, the AMOC in the coarse resolution simulation penetrates up to 80 • N.These differences between the eddy-permitting and non-eddy-permitting configurations can originate from the different parametrizations (e.g.larger diapycnal and lateral diffusivity in the ORCA1 configuration) which impose the differences in the advection of salinities, as shown in Fig. 7a  and b.The maximum AMOC is about 10 Sv at 600 m near 25 • N in the ORCA025 simulation and about 18 Sv at 650 m near 35 • N in the ORCA1 simulation.This maximum in the coarse resolution simulation is within the range of values found in the PMIP2 models (between 13.8 and 20.8 Sv) (Otto-Bliesner et al., 2007) whereas it is particularly weak in the eddy-permitting simulation.Hence, it seems that the strength and geometry of the LGM AMOC can vary between climate models (Otto-Bliesner et al., 2007;Weber et al., 2007) and also between different regimes or choice of the resolution and parametrization in the models.Weak AMOC (similar to the ORCA025 simulation) has also been observed in some LGM climate model simulations (Schmittner, 2003;Shin et al., 2003;Brandefelt and Otto-Bliesner, 2009) and eddy-permitting regional simulations (Yang et al., 2006).It has been attributable to the expansion of the Antarctic Bottom Water (AABW) brought about by the increased rates of sea-ice formation (saltier and denser AABW).The large sea-ice covers in our simulation as well as the atmospheric forcing fields from the climate simulation by Brandefelt and Otto-Bliesner (2009) could also explain the weak surface thermohaline circulation and an abyssal ocean filled with extremely cold and saline waters.The larger value diagnosed in the ORCA1 simulation suggests that the parametrization also plays a key role in the structure of the LGM AMOC.
The northward heat transport in the Atlantic Ocean is controlled by the strength of the overturnings (NADW and AABW) and the temperature difference between the northward and southward flows.The northward heat transport (NHT) in the eddy and the non-eddy-permitting experiments is relatively similar in the Atlantic Ocean (Fig. 8b) which is explained by the combination of weak surface and deep circulations in the ORCA025 simulation and stronger NADW and AABW in the ORCA1 simulation.The transport is close to zero in the South Atlantic Ocean, suggesting that the global tropical cooling is achieved by a stronger meridional heat transport in the Indo-Pacific Ocean.Therefore, the meridional heat transport for the global ocean attains a maximum of 1.5 PW near 15 • N and 15 • S.These features are consistent with previous model diagnostics (Ganopolski et al., 1998;Weaver et al., 1998;Hewitt et al., 2001), but not with recent PMIP2 model results, which show larger transports between 40 and 60 • N for the global ocean, a larger transport in the Atlantic Ocean (south of 30 • N) and a smaller transport south of 40 • N in the Pacific (Murakami et al., 2008).These large discrepancies are connected with the differences in the meridional overturning circulations and the reorganization of the temperatures and water masses (due to the sea-ice formation) between our simulations and those from the PMIP2 experiment.The main differences in the NHT between the ORCA025 and ORCA1 experiments occur in the Indo-Pacific Ocean and the Southern Ocean, in the area where the eddy field is particularly important (Fig. 8a and c).The eddies participate in about 0.5 PW for the meridional transport of heat between 20 and 60 • S and in the Indo-Pacific tropical band, suggesting that the comparison between models and the proxy might differ in these regions.

Sea surface temperature
The model SSTs are firstly compared to the MARGO palaeoproxy reconstructions (Waelbroeck et al., 2009) in 4 latitudinal bands by using the Taylor diagram representation.The models are distinguished by a number (1 or 2, for ORCA1 and ORCA025 respectively) and the MARGO data is symbolized by a "star" on the x axis (as the standard deviation of the data set).The model data standard deviation is given by the radial distance from the origin of the diagram.The correlation coefficient is represented by the azimuthal angle between the simulated SST field and the MARGO data, and the centred RMSE between the simulated data and the reference MARGO data is given by the distance between the number and the symbol (in other words, the closer the model is to the reconstructed field, the lower is its centred RMSE).
The model SSTs are uncorrelated with the MARGO reconstruction between 50 and 90 • N (Fig. 9), (R < 0.25 for each period).Moreover, the standard deviations of the simulated SSTs are close to zero, suggesting that the simulated SSTs in this region are relatively homogenous (viz.a small variability of the SSTs).The North Atlantic zonal SST gradients diagnosed in the MARGO reconstruction are thus not resolved by the models.This can be explained by the presence of large North Atlantic sea-ice covers in our simulations which prevent SST gradients in this area.Consequently, the models' RMSE are between 3.0 and 4.5 • C and the skill scores are close to zero (Table 1), i.e. the models fail to represent the SSTs in this region.Between 25 and 50 • N (Fig. 10), the models have similar behaviour.They capture the variability and the large-scale pattern of the MARGO SSTs.The correlation coefficients R are between 0.7 and 0.8 for the ANN and JAS reconstructions and between 0.8 and 0.9 for the JFM reconstruction.Note that the coarse-resolution simulation has the best correlation coefficients R and the lowest RMSE.The skill scores of the models are between 0.83 and 0.95, with the best score reached by the ORCA1 experiment.Between 25 • S and 25 • N (Fig. 11), the models also have a similar behaviour, which, in view of the available proxy-data, is not in favour of the high-resolution simulation.For each periods, they underestimate the variability of the MARGO SSTs.For the ANN and JFM reconstructions, the correlation coefficients R are near 0.5, the models' RMSE slightly above 3 • C and the skill scores S between 0.55 and 0.66, suggesting that the inter-tropical SSTs are not fully consistent with the reconstructed SSTs.Hence, the tropical upwellings are poorly resolved by the models.The ORCA1 model has also better skill scores than the eddy-permitting simulation.Between 25 and 50 • S (Fig. 12), the models are highly correlated with regards to the palaeo-reconstructions (the correlation coefficients R are between 0.95 and 0.99) and the RMSE is near 2 • C.Although the ORCA025 simulation is slightly closer (better R) to the MARGO reconstruction in this region for the JFM and JAS periods, the models have the same skill scores: 0.97 for the annual mean reconstruction, 0.98 for JFM and 0.97 for JAS.Overall, this analysis shows that for each latitudinal band the coarse-resolution and eddypermitting simulations have similar skills, therefore, the permitted eddies do not contribute significantly to a better consistency with regards to the large-scale features proposed by palaeo-proxy reconstructions.
By performing similar statistical analyses on hindcast simulations of the present-day climate (see Figs. S5 to S9 in Supplement), it is shown that, relatively to the coarse resolution, the eddy-permitting regime improves the statistical scores (better correlation and variability) for replicating the global World Ocean Atlas 1998 SST observations, particularly in the 25-50 • N band and in the 25 • S-25 • N band for the annual mean data set.The improvements are weaker in the other regions.For each period (ANN, JFM and JAS), the predicted SSTs are highly correlated to the observed SST state (i.e.R > 0.9) and the models capture the observed SST variability in the latitudinal bands.In considering only the observations at the MARGO core locations (Figs.S10 to S13 0.0 0.1 0. 2 0 .30 .40 .5 0 .6 0 .7 0 .8 0 .90 .9 5 0.9 in Supplement), it is shown that the improvements introduced by the higher spatial scales are less obvious.Although, the correlations coefficients and standard deviations are improved due to the reduced size of the sample, the ORCA1 and ORCA025 show similar behaviour for replicating the observed SSTs.This suggests that the distribution of the LGM data has an impact on the statistical scores and that more sediment cores may thus be required to draw reliable conclusions about the improvements introduced by the higher spatial scales for reproducing the modern and the past global SSTs.
The statistical analysis at regional scales (Table 2) provides insight on the plausible local improvements where the eddy activity is particularly large and could play an important role on the SST structures.The models and data are significantly correlated in all the areas except in the Southern Ocean Indian sector for the ANN and JAS periods, and in the Southern Ocean Pacific sector for the JFM period.For the ANN and JFM periods, the ORCA1 and ORCA025 simulations are highly correlated with the MARGO data (p values < 0.001) in the Agulhas and the Gulf Stream regions.The statistical analysis reveals that for the regions where models and data are significantly correlated, the high-resolution 0.0 0.1 0. 2 0 .30 .40 .5 0 .6 0 .7 0 .8 0 .90 .9 5 0.9 9 C o r r e l a t i o n 0.0 0.8 1.6 2.4 3.2 4.0 4.8 5.6 6. simulation often improves the statistical scores (higher correlation coefficient, slope closer to 1 and intercept closer to 0 • C), particularly for the Agulhas current and the Agulhas retroflection.For the JFM period, the correlations between models and proxies are statistically highly significant in all regions except for the Pacific sector in the Southern Ocean (Table 2).It is however less obvious to diagnose improvement induced by the permitted eddies for this period.For the JAS period, the correlations of the ORCA1 and ORCA025 simulations with the MARGO data are statistically highly significant only in the Agulhas and the Gulf Stream region.Hence, on the basis of the available Southern Hemisphere proxy data in the mid-to high latitudes, and the two statistical approaches presented here, it seems that the eddy-permitting simulation mainly improves the regional representation of the Southern Ocean SSTs.Similarly, modern SSTs are better replicated by the eddy-permitting simulation in the Gulf Stream and Agulhas regions whereas the improvements are weaker in the Southern Ocean and in the region where the number of LGM SST estimates is too small (Table S1 in Supplement).

Sea-ice cover
In complement to the SST reconstructions, sea-ice reconstructions are available for the LGM interval (Pflaumann et al., 2003;De Vernal et al., 2006;Nørgaard-Pedersen et al., 2003;Gersonde et al., 2005).For the JFM period (Fig. 3a  and b), the models simulate sea-ice cover in the northwestern North Atlantic and in the Labrador Sea as well as Nordic Seas.However, this North Atlantic extent reaches almost 40 • N, whereas the margin is located between 50 and 60 • N in the reconstructions (Pflaumann et al., 2003;De Vernal et al., 2006).The ice-free conditions in the Irminger current up to Iceland are also less pronounced in our simulations.These large sea-ice covers in our simulations might be inferred from the atmospheric condition applied at the ocean surface.The near-surface forcing fields originate from a quasi-equilibrated simulation which showed particularly cold conditions over the North Atlantic region (Brandefelt and Otto-Bliesner, 2009).In the Southern Ocean, our results are also not fully consistent with the palaeoreconstructions.As the majority of the PMIP simulations (Roche et al., 2012), the models simulate a larger sea-ice extent than in the proxy reconstructions by Gersonde et al. (2005).It is also known that ocean/sea-ice models have difficulties in simulating the variance of the Southern Ocean seaice extent, thickness and volume for the present-day period (Arzel et al., 2006;Zunz et al., 2013).Consequently, it is also difficult to draw reliable conclusions for the LGM Southern Ocean sea-ice dynamics.For the JAS period (Fig. 4a and b), the LGM sea-ice fractions are reduced in the Labrador Sea, the central North Atlantic and the Norwegian Sea, suggesting possibly ice-free conditions during the glacial period as  (2005).However large uncertainties in the proxy-data exist in the Nordic Seas (De Vernal et al., 2006), which make the interpretation and confrontation between model and proxies more complex.The central Arctic Ocean and the western part of Fram Strait have perennial sea-ice in accordance with the reconstructions of Nørgaard-Pedersen et al. (2003).In the Southern Ocean, the modelled maximal seaice area is around 39 × 10 6 km 2 (Fig. 5b), which is similar to a reconstruction (Gersonde et al., 2005) with sea-ice extending to 45 • S in the Atlantic and Indian oceans and 55 • S in the Pacific.

Atlantic meridional overturning circulation
The reconstruction of the structure of the LGM Atlantic thermohaline circulation is challenging and palaeo-proxy data suggest various geometries.Some palaeo-proxy reconstructions based on geochemical tracers suggest that the LGM Atlantic overturning was as strong as it is today and the return flow occurred at about 2000 m (Yu et al., 1996;Lynch-Stieglitz et al., 2007;Gherardi et al., 2009;Lippold et al., 2012).Other suggest that the overturning rate was weaker than today (Lynch-Stieglitz et al., 1999;Hesse et al., 2011).Reconstructions of the abyssal transport suggest that the circulation was not sluggish and was comparable in strength with present-day transports (Yu et al., 1996;Lynch-Stieglitz et al., 2007), whereas others suggest weak renewable or poorly ventilated water mass (Gherardi et al., 2009;Lippold et al., 2012).Although the strength of the AMOC varies between the eddy-permitting and non-eddypermitting regime, both simulations capture an AMOC in the upper 2000 m and a large intrusion of an AABW below, as proposed by the reconstructions.The difference in the strength of the overturning between the two regimes may be due to the different diffusive parametrizations (e.g.Nilsson et al., 2003;Prange et al., 2003).In our simulations, the deep circulation is particularly stagnant.This can be inferred from the non-equilibrated solution in the deep ocean and from the highly dense and homogenous waters that filled the abyssal basins and prevent density contrast between north and south for maintaining a strong abyssal overturning.In the models of Butzin et al. (2005) and Hesse et al. (2011), a shallowed and weakened AMOC as it is simulated in the eddypermitting experiment, but an intensified AABW give the best agreement for capturing the general δ 13 C distribution (i.e.water masses geometry) derived from sediment analysis.These results were also obtained by an enhanced northward sea-ice export in the Southern Ocean which can maintain an expanded and dense AABW.

Conclusions
This study presents the response of the first global eddypermitting ocean general circulation model forced with atmospheric fields representing a 49 yr sample of the Last Glacial Maximum climate and extracted from the climate model simulation carried out by Brandefelt and Otto-Bliesner (2009).In order to identify whether eddies contribute to a better consistency between models and palaeoproxies, a coarse resolution simulation is also designed and confronted to the palaeo-reconstructions.The results reported here are consistent with those from the CCSM3 quasi-equilibrated simulation (Brandefelt and Otto-Bliesner, 2009), but in other respects may differ from the PMIP model analyses due to the cold atmospheric conditions used to force our models.As pointed out by Brandefelt and Otto-Bliesner (2009), their "new equilibrium differs substantially from the first quasi-steady-state with 1.1 • C colder global mean temperature and regional differences of 5-15 • C in the North Atlantic region and a 30 % reduction of the strength of the AMOC", this associated with the equilibration of the abyssal ocean.The most significant discrepancies have been diagnosed for the seasonality of the Southern Ocean sea-ice fractions which may be inferred from these cold atmospheric conditions.As a consequence of the larger sea-ice areas simulated in the models, the global ocean surface temperature and salinity, as well as the deep-water formation are modified.Substantial differences in the ORCA1 and ORCA025 distribution of the water properties (temperature and salinity), the AMOC and the deep mixed layer may be also inferred from the different mixing parametrizations between the two experiments.These differences raise the question whether the tuning applied for the present-day climate is suitable for past climate simulations.The comparison between eddy-permitting and non-eddypermitting regimes shows that the permitted eddies modifies the sea surface states (SSTs, SSS, sea-ice cover) particularly in the mid-latitudes and in the Southern Ocean.Their impacts are weak in the northern North Atlantic and the Arctic Ocean, which is particularly isolated from the North Atlantic during the LGM period.The performances of the two model regimes for reproducing the SSTs are also summarized by using Taylor diagrams (a quantification of the skill score) as well as by regional statistical tests between the modelled SSTs and the proxy reconstructions.The models perform well between 25-50 • S and 25-50 • N, i.e. in the region where the mesoscale eddies are the most important.They are less adequate in the Northern Hemisphere high latitudes (due to large sea-ice cover) and in the tropical regions, implying that the equatorial dynamics and coastal upwelling zones are still poorly represented in the eddy-permitting model.The lack of vertical resolution as well as the choice of parametrization of the upwelling processes can also have an impact on the quality of the simulated tropical SSTs.The large covers of sea-ice in high latitudes is significantly different than the palaeo-oceanographic reconstructions.
It is hence not obvious that the eddy-permitting simulation contributes significantly to the improvement of the sea-surface state as estimated from the available palaeooceanographic reconstructions.Some differences in the structure of the coarse-resolution and eddy-permitting SSTs are noticeable in specific regions, e.g. the width of the Agulhas current along the coast of Mozambique and South Africa or the Kuroshio are narrower in the eddy-permitting simulation.These results may locally yield a better consistency between model and proxy data but do not contribute significantly to the global statistical score.Similar statistical tests applied to hindcast simulations of the present-day period  reveal that the number of sediment cores might also have an impact on the statistical score.Hence, to draw firm conclusions on the benefit of the high resolution simulation, more reconstructions are required in areas such as the ACC and the mid-latitudes (e.g. the Gulf Stream and the Agulhas current) where the data coverage is limited and the mesoscale dynamics play a key role for the surface temperature structure and the climate variability.
It is thus challenging to reconstruct past climate conditions such as the LGM climate but the synergy of numerical models and climate proxy reconstructions can provide insights into this climate at global and regional scales.The consistencies and discrepancies between models and climate proxies also reveal the complexity of the reconstruction techniques.Each method contains uncertainties or simplification of the reality.Our analysis does not consider the uncertainties in the palaeoproxies, as for example the conflicting conclusions remaining in the North Atlantic sea-ice reconstructions by various palaeo-sensors (De Vernal et al., 2006).In summary, this investigation indicates that a higher model resolution is not a panacea yielding the hoped-for better correspondence between simulations and proxy archives.In addition, it is shown that both local and global statistical analyses are required to identify model behaviour relative to the climate proxy reconstructions.Finally, until more experiments, either corroborating or invalidating the results of the present study, have been undertaken, only provisional conclusions may be drawn.One of these is that it may be more appropriate to aim at improving the parametrizations of small-scale processes (e.g.upwellings and vertical mixing) and Southern Ocean dynamics rather than to focus on increasing model resolution.

Fig. 1 .
Fig. 1.Map of (a) the annually averaged sea surface temperature (SST) in the LGM eddy-permitting simulation and (b) the difference between the coarse resolution and the eddy-permitting simulation.Units in • C.

Fig. 2 .
Fig. 2. Map of (a) the annually averaged Sea Surface Salinity (SSS) in the LGM eddy-permitting simulation and (b) the difference between the coarse resolution and the eddy-permitting simulation.Units in PSU.

Fig. 3 .
Fig. 3. Polar stereographic map of the boreal winter (JFM) sea-ice fraction (a) in the Northern Hemisphere and (b) in the Southern Hemisphere as simulated by the eddy-permitting model.Units in %.

Fig. 4 .Fig. 5 .
Fig. 4. Polar stereographic map of the boreal summer (JAS) sea-ice fraction (a) in the Northern Hemisphere and (b) in the Southern Hemisphere as simulated by the eddy-permitting model.Units in %.

Fig. 6 .
Fig. 6.Maps of (a) the boreal winter (JFM) and (b) the boreal summer (JAS) mixed layer depth (m) computed in the eddy-permitting simulation.

Fig. 8 .
Fig. 8. Meridional ocean heat transport in the eddy-permitting (ORCA025) and coarse (ORCA1) resolution simulations for (a) the Global Ocean, (b) the Atlantic basin and (c) the Indo-Pacific basin.The dashed line represents the contribution of the eddies in the meridional heat transport.Units in PW = 10 15 W.

Fig. 9 .
Fig. 9. Taylor diagrams summarizing the correlation coefficient (along the arc), the standard deviations (on the axis) and RMSE (grey solid arcs of circles) for the LGM annual, boreal winter (JFM) and boreal summer (JAS) modelled and reconstructed sea surface temperatures between 50 and 90 • N. The proxy-data locations are shown in the upper right-hand diagram.The NEMO-ORCA1 is symbolized with the number 1, NEMO-ORCA025 with the symbol 2.

Fig. 10 .Fig. 11 .Fig. 12 .
Fig. 10.Taylor diagrams summarizing the correlation coefficient, the standard deviations and RMSE for the LGM annual, boreal winter (JFM) and boreal summer (JAS) modelled and reconstructed sea surface temperatures between 25 and 50 • N. The proxy-data locations are shown in the upper right-hand diagram.The NEMO-ORCA1 is symbolized with the number 1, NEMO-ORCA025 with the symbol 2.

Table 1 .
Taylor (2001)as defined byTaylor (2001)for each latitudinal band in the ORCA025 and ORCA1 experiments.The best scores reached are in bold.

Table 2 .
Regional statistical analysis (linear regression test)between the simulated and reconstructed SSTs summarized by the p value, the correlation coefficient, the slope and the intercept (in • C) of the regression line.N is the sample size of the regional MARGO data for each period.In bold, the most significant values.