Effects of Last Glacial Maximum (LGM) sea surface temperature and sea ice extent on the isotope–temperature slope at polar ice core sites

. Stable water isotopes in polar ice cores are widely used to reconstruct past temperature variations over several orbital climatic cycles. One way to calibrate the isotope– temperature relationship is to apply the present-day spatial relationship as a surrogate for the temporal one. However, this method leads to large uncertainties because several factors like the sea surface conditions or the origin and transport of water vapor inﬂuence the isotope–temperature temporal slope. In this study, we investigate how the sea surface temperature (SST), the sea ice extent, and the strength of the Atlantic Meridional Overturning Circulation (AMOC) affect these temporal slopes in Greenland and Antarctica for Last Glacial Maximum (LGM, ∼ 21 000 years ago) to preindus-trial climate change. For that, we use the isotope-enabled


Introduction
Stable isotopologues of water (H 16 2 O, H 18 2 O, and HD 16 O, hereafter called stable water isotopes) are integrated tracers of climate processes occurring in diverse parts of the hydrological cycle (Craig and Gordon, 1965;Dansgaard, 1964). Because of their differences in mass and symmetries, an isotopic fractionation happens at each phase change of water. This process is reflected by a change in the water isotope ratio values, expressed hereafter in the usual δ notation (as δ 18 O and δ 2 H with respect to the Vienna Standard Mean Ocean Water, V-SMOW, if not stated otherwise). As a result, water isotopes have been widely used to describe past variations of the Earth's climate. For example, their measurements in polar ice cores made it possible to reconstruct the temperature variations over several glacial-interglacial cycles Jouzel, 2013, and references therein;NEEM Community Members, 2013).
For such a reconstruction, the present-day isotopetemperature spatial slope can be taken as a surrogate for the temporal gradient at a given site. For example, a spatial slope of 0.80 ‰ • C −1 for δ 18 O in Antarctica was calculated based on a compilation of measured surface temperatures and δ 18 O of snow at various locations on the continent (Masson-Delmotte et al., 2008). However, this method often leads to a large error in the temperature reconstructions because the temporal isotope-temperature slope depends on many factors like the sea surface temperature (SST) (Risi et al., 2010), the sea ice extent (Noone and Simmonds, 2004), the ice sheet elevation (Werner et al., 2018), and the origin and transport of water vapor (Casado et al., 2018). For example, it has been suggested that the relationship between temperature and the isotopic signature for warmer interglacial periods in East Antarctica can vary among ice core sites, with an error in the temperature reconstruction that can reach up to 100 % (Sime et al., 2009;Cauquoin et al., 2015). In Greenland, the use of the spatial relationship between the δ 18 O in Greenland ice core records and surface temperature to evaluate the local temperature variations during the last deglaciation leads to a large uncertainty of a factor of 2 (Jouzel, 1999;Buizert et al., 2014). Recently, Buizert et al. (2021) proposed a reconstruction of surface cooling in Antarctica during the Last Glacial Maximum (LGM, ∼ 21 000 years ago) using borehole thermometry and firn properties of different ice cores. Based on these results, they proposed new estimates of temporal δ 18 Otemperature slopes at these ice core stations, varying from 0.8 to 1.45 ‰ • C −1 .
The LGM is a period with full glacial conditions and represents the beginning of the last deglaciation. It is one of the key climate periods chosen by the Paleoclimate Modeling Intercomparison Project (PMIP, Kageyama et al., 2018Kageyama et al., , 2021 because it allows evaluating how well state-of-the-art models are able to simulate climate changes as large as those expected in the future. In addition to being very different from the preindustrial climate (PI), the LGM period also offers a wealth of isotope proxy data, including stable water isotopes in polar ice cores, for an in-depth comparison with outputs from isotope-enabled models (Lee et al., 2008;Risi et al., 2010;Werner et al., 2016Werner et al., , 2018. One way to capture the physical processes influencing the temporal isotope-temperature slope in polar regions is the use of atmospheric general circulation models (AGCMs) equipped with prognostic stable water isotopes. Such models can simulate different climate conditions, like the LGM and PI periods. Moreover, the use of isotope-enabled AGCMs in combination with isotopic observations allows us to investigate the physical processes controlling the variations of isotopic delta values at a given site. This method makes it possible to estimate the temporal isotope-temperature slope for LGM-to-preindustrial climate change (Lee et al., 2008;Risi et al., 2010;Werner et al., 2018). Even if such models simulate various temporal isotope-temperature slopes, implying that processes like water vapor transport, post-depositional effects, or the polar atmospheric boundary layer are poorly or not represented (Krinner et al., 1997;Werner et al., 2000;Casado et al., 2018), these models are very useful for evaluating the sensitivity of the temporal slopes to parameters like the change in elevation (Werner et al., 2018).
Ocean surface conditions represent one of the factors that influence LGM-PI isotope changes (Risi et al., 2010;Noone and Simmonds, 2004). Two reconstructions of SST and one of sea ice extent during the LGM period have been released recently. Paul et al. (2021) reconstructed both the SST and the sea ice extent fields based on faunal and floral assemblage data from the Multiproxy Approach for the Reconstruction of the Glacial Ocean Surface (MARGO) project and several recent estimates of the LGM sea ice extent. The Data-Interpolation Variational Analysis (DIVA) software was used to optimally interpolate sparse SST reconstruction data. The resulting reconstruction was called GLOMAP (Glacial Ocean Map). Tierney et al. (2020a) reconstructed the LGM SST field with a different method by combining a large collection of geochemical proxies for sea surface temperature with simulation outputs from the isotope-enabled model iCESM1.2 (Brady et al., 2019) using an offline data assimilation technique to produce a field reconstruction of LGM temperatures. Tierney et al. (2020a) LGM cooling is globally larger than in GLOMAP (3.6 and 1.7 • C, respectively), with possible impacts on LGM-to-PI isotope changes and their temporal relationship with nearsurface air temperature. In addition, other SST and sea ice fields, with different characteristics compared to the reconstructions of LGM sea surface conditions described above, can be extracted from atmosphere-ocean coupled model simulations like MIROC 4m (Obase and Abe-Ouchi, 2019).
Are air temperatures near the surface and the isotopic composition of precipitation in the polar regions influenced by LGM-to-PI changes in SST and sea ice distribution in the same way? What are the underlying dynamics, for example, in terms of changes in concentrations and transport of water vapor? To answer to these questions, we performed multiple simulations with the isotope-enabled AGCM ECHAM6wiso driven by different LGM SST and sea ice boundary conditions. We evaluate the modeled LGM-PI δ 18 O anomalies with available observations, and we investigate how the SST and sea ice extent patterns influence the model-data agreement on a global scale and at polar ice core stations. The influence of ocean circulation, particularly the strength of the Atlantic Meridional Overturning Circulation (AMOC), on sea surface conditions and by extension on our modeled δ 18 O of meteoric water is also investigated. Finally, the impacts of the sea surface boundary conditions on the δ 18 Otemperature slopes for LGM-to-preindustrial climate change are evaluated and discussed for Greenland and Antarctic ice core stations.

ECHAM6-wiso
ECHAM6  is the sixth generation of the atmospheric general circulation model ECHAM, developed at the Max Planck Institute for Meteorology. It consists of a dry spectral-transform dynamical core, a transport model for scalar quantities other than temperature and surface pressure, a suite of physical parameterizations for the representation of diabatic processes, and boundary datasets for externalized parameters (trace gas and aerosol distributions, land surface properties, etc.). ECHAM6 forms the atmospheric component of the fully coupled Earth system model MPI-ESM Mauritsen et al., 2019). The implementation of the water isotopes in ECHAM6 as part of MPI-ESM has been described and evaluated in detail by Cauquoin et al. (2019b), and this model version has been labeled ECHAM6-wiso. At a later stage, Cauquoin and Werner (2021) updated the water isotope module of ECHAM6-wiso in several respects. The supersaturation has been slightly retuned, the kinetic fractionation factors for the evaporation over the ocean are now assumed to be independent of wind speed, and the isotopic content of snow on sea ice is taken into account for sublimation processes in sea-ice-covered regions. The latter leads to a stronger depletion of surface water vapor over such sea-ice-covered areas (while the surface temperature remains the same). As a consequence, this change is expected to contribute to a steeper temporal isotope-temperature slope over sea-ice-covered areas.
2.2 Sea surface temperature and sea ice extent boundary conditions for LGM conditions

SST
The Tierney et al. (2020a) SST reconstruction has a larger and more homogeneous cooling than GLOMAP, except for the high southern latitudes at which the Pacific sector cools more than the Atlantic sector ( Fig. 1). On the other hand, the LGM cooling in the northern North Atlantic Ocean is stronger in GLOMAP than in the Tierney et al. (2020a) reconstruction (−5.4 and −4.8 • C, respectively, see Table 4 in Paul et al., 2021). These differences between the two SST reconstructions are due to the use of different proxy datasets for the reconstructions (geochemical proxies only for Tierney et al., 2020a, MARGO dataset for GLOMAP) and the methods applied to produce SST gridded maps from scattered observations (see Sect. 1). For their offline data assimilation technique, Tierney et al. (2020a) used results from the coupled climate model iCESM1.2, which shows one of the largest coolings among the PMIP4 models ( Fig. 1b of Kageyama et al., 2021). In addition to these two reconstructions, we used SST and sea ice extent outputs from a MIROC 4m LGM simulation (Obase and Abe-Ouchi, 2019) with oscillating AMOC strength. The global LGM cooling is between −2.3 and −2.7 • C according to the considered simulations ( Fig. 1), i.e., higher than GLOMAP and lower than the Tierney et al. (2020a) reconstruction. The main specificity of MIROC 4m LGM SST is a very strong cooling in the North Atlantic (more than 10 • C, Fig. 1) and more uniform temperature anomalies between −2 and −4 • C in the other areas, including off the coast of Greenland. We extracted the MIROC 4m SST outputs, averaged over a 100-year period, at two different times of the LGM simulation depending on the AMOC strength: during a weak AMOC phase (average AMOC index was equal to 8.44 Sv) and a strong AMOC phase (19.95 Sv). A weaker AMOC during LGM implies larger cooling in the North Atlantic ( Fig. 1) and more extended sea ice (Fig. 2), with less cooling in the Southern Ocean. The strong AMOC phase period in the MIROC 4m simulation was selected in the middle of the AMOC peak ( Fig. S1   On the other hand, sea ice is more extensive in the Northern Hemisphere for MIROC 4m in the weak AMOC phase. For the stronger AMOC case, a decline of the sea ice in the Northern Hemisphere is seen, accompanied by weaker cooling (see Sect. 2.2.1). In its parameterization, MIROC 4m uses a threshold of 95 % for the sea ice fraction to allow sub-grid "sea ice leads". This threshold is not rigid, but it is difficult to exceed sea ice concentrations of 95 % unless there is significant convergence of sea ice. Consequently, while the sea ice is, on average, more extensive in the north in MIROC 4m for the weak AMOC phase compared to the GLOMAP reconstruction, the sea ice area fraction in grid cells near coastal areas like Greenland is lower in MIROC 4m than in GLOMAP (95 %-98 % against 100 %, respectively).

Model setup and experiments
We performed an ensemble of LGM simulations with ECHAM6-wiso, forced with different combinations of SST and sea ice boundary forcings presented in Sect. 2.2. All our PI and LGM simulations are designed with respect to the PMIP4 protocol (e.g., greenhouse gas and orbital conditions). The LGM SST boundary fields are expressed relative to the Atmospheric Model Intercomparison Project (AMIP, Eyring et al., 2016) mean SST (averaged over the period 1870 to 1899) used for the preindustrial simulations ( Table 1). The GLOMAP reconstruction has the advantage of providing a monthly climatology of LGM SST and sea ice extent, while only annual mean SST is available from the reconstruction by Tierney et al. (2020a), without a sea ice map distribution. So, the Tierney et al. (2020a) LGM SST for ECHAM6-wiso was produced by taking the annual mean SST anomaly from Tierney et al. (2020a) and adopting the monthly climatology temperature variability from GLOMAP. Since we also used GLOMAP sea ice extent data in this case, the SST was adjusted slightly to maintain consistency (e.g., SST set to freezing temperature where there is sea ice). In order to investigate the impact of sea ice extent on our isotope results and the related isotope-temperature slopes for LGM-to-PI climate change, we used LGM SST outputs from MIROC 4m simulations combined with sea ice extent data from the same MIROC 4m simulations or from the GLOMAP dataset. To evaluate the LGM-PI anomalies, we performed several PI simulations with different sea ice boundary conditions depending on the setup of LGM experiments using climatological monthly mean sea ice area fractions from AMIP or MIROC 4m coupled simulations (Table 1). The prescribed LGM ice sheet is GLAC-1D (Tarasov and Peltier, 2002;Tarasov et al., 2012Tarasov et al., , 2014Abe-Ouchi et al., 2013;Briggs et al., 2014) for all LGM simulations. As with SST and sea ice distribution, mean δ 18 O of surface seawater needs to be prescribed. For the PI simulations, we used the δ 18 O reconstruction from the global gridded dataset of LeGrande and Schmidt (2006). As no equivalent dataset of the δ 2 H composition of seawater exists, the deuterium isotopic composition of the seawater in any grid cell has been set equal to the related δ 18 O composition, multiplied by a factor of 8, in accordance with the observed relation for meteoric water on a global scale (Craig, 1961). As in Werner et al. (2018), a prescribed glacial seawater enrichment of +1 ‰ and +8 ‰ is assumed for δ 18 O and δ 2 H in the LGM simulations, respectively. The PI and LGM simulations were run for 60 and 120 model years, respectively, and we used the last 30 model years for our analyses. The simulation characteristics are summarized in Table 1. Two additional sensitivity simulations have been performed to evaluate the impacts of lower MIROC 4m sea ice area fraction in coastal grid cells (Sect. 2.2.2) and the consideration of the isotopic composition of snow on sea ice in ECHAM6-wiso (Sect. 2.1) on the modeled δ 18 O p -temperature temporal slopes between LGM and PI (see Supplement Sect. S1). Also, an LGM simulation using the PMIP3 ice sheet reconstruction instead of GLAC-1D (see Fig. 3b and d of Werner et al., 2018, respectively) has been performed to evaluate the impact of ice sheet topography on the isotopically enriched bias in Antarctica (Sect. S1).

Observational data
To evaluate the modeled δ 18 O of precipitation and snow values at ice core stations, we use a selection of 6 Greenland and 11 Antarctic ice cores for the preindustrial and LGM climates (Fig. 3). The observed δ 18 O values were defined as averages over the last 200 years for the preindustrial pe- LGM and PI sea ice area fractions used as boundary conditions for ECHAM6-wiso simulations. LGM_miroc4m_strong_AMOC_sst_glomap_sic MIROC 4m GLOMAP Mean PI SST and sea ice from AMIP AMOC oscillation: strong phase LGM_miroc4m_strong_AMOC_sst_and_sic MIROC 4m MIROC 4m Mean PI SST from AMIP and PI sea ice from MIROC 4m AMOC oscillation: strong phase riod and in the 21 ± 1 ka period for the LGM. We also use LGM-PI δ 18 O anomalies from five (sub)tropical ice cores that are reported in Table 2 of Risi et al. (2010). The ice core data used in this study are summarized in Table 2. Similarly, we use temperature reconstructions based on borehole reconstructions or δ 15 N for seven Antarctic cores and one Greenland ice core (Buizert et al., 2021;Dahl-Jensen et al., 1998)

Evaluation of ECHAM6-wiso under LGM conditions
We evaluate the global distribution of δ 18 O of precipitation (δ 18 O p ) changes between LGM and PI ( δ 18 O p ) from our different ECHAM6-wiso simulations. Figure 4 shows the comparison of modeled δ 18 O p anomalies with isotope measurements from ice cores and speleothems for the simulation LGM_miroc4m_sst_and_sic (i.e., SST and sea ice boundary conditions from the MIROC 4m simulation in a weak AMOC phase). Well-known patterns of global δ 18 O p distribution are found in ECHAM6-wiso, like the negative anomalies across Canada, Greenland, and northern Europe due to the presence of glaciers in these areas during LGM period (Fig. 4c). Generally, negative δ 18 O p anomalies are also simulated over Antarctica and the Southern Ocean, where the LGM cooling is stronger compared to lower southern  latitudes (Fig. 4a). In the middle to low latitudes, δ 18 O p is mainly controlled by precipitation anomalies (Fig. 4b).
For example, lower modeled precipitation in the Amazonian area, over parts of Southeast Asia, and in the western Pacific Ocean during the LGM leads to positive modeled δ 18 O p anomalies. Despite some biases in modeled δ 18 O p , like in southern Amazonia ( Fig. 4c and d) where negative anomalies are measured in ice cores (between −2 ‰ and −6 ‰, see green dots in Fig. 4d) while positive anomalies are simulated (between 0 ‰ and 4 ‰), modeled δ 18 O p anomalies are in rather good agreement with observations from ice cores and speleothems (Fig. 4d). The isotope fractionation is mainly controlled by changes in temperature and in the water cycle. Even though all the ECHAM6-wiso simulations show a similar global distribution of 2 m air temperature (T 2 m ) and precipitation responses to the various SST and sea ice boundary fields, we find some differences too (Figs. S2 and S3). As expected, the modeled global cooling using SST from GLOMAP is lower, while it is stronger when using SST from Tierney et al. (2020a) (cooling of −4 and −6.3 • C, respectively). Average T 2 m anomalies in the middle range are obtained when using the MIROC 4m SST fields (between −4.4 and −5.3 • C depending on the MIROC 4m data used). The temperatures over sea-icecovered areas are largely impacted by the sea ice forcings used (i.e., GLOMAP or MIROC 4m). The modeled T 2 m anomalies over the Southern Ocean vary between −10 and −15 • C with an extensive LGM sea ice (GLOMAP), while the cooling is only between −4 and −10 • C when ECHAM6wiso is forced by a less extensive one (MIROC 4m sea ice). For the Arctic region, a strong cooling is simulated with the very extensive sea ice from MIROC 4m with a weak AMOC phase (a cooling of more than 20 • C), more than with the sea ice from GLOMAP (between −20 and −10 • C). The different SST boundary conditions have a strong influence on the precipitation anomalies, especially at middle to low latitudes including the western Pacific area and the East Asian monsoon region (Fig. S3). All these differences in T 2 m and precipitation responses have profound impacts on modeled δ 18 O p anomalies (Fig. S4) and their agreements with observations ( Fig. 4c and d).
The statistics of δ 18 O p model-data agreements are shown for our different ECHAM6-wiso simulations in Table 3. The best model-data agreement in terms of model-data slope (1 is a perfect match) is found when using SST and sea ice from GLOMAP (slope = 0.70) as boundary conditions for ECHAM6-wiso, but a better coefficient of determination (r 2 ) and root mean square error (RMSE) are obtained with LGM SST from Tierney et al. (2020a) (r 2 = 0.58 and RMSE = 3.5 ‰). We notice worse model-data agreements in δ 18 O p changes when both SST and sea ice changes from MIROC 4m simulations are provided as sea surface boundary conditions (slopes lower than 0.582 and RMSE higher than 4.1 ‰). This is in agreement with Werner et al. (2018), who showed worse model-data agreement when using SST and sea ice boundary conditions from a coupled model instead of reconstructed ones. The substitution of MIROC 4m sea ice changes by GLOMAP ones, which are more extensive in the south and less extensive in the north, improves the δ 18 O p model-data agreement for all cases (i.e., weak or strong AMOC phase). For example, the model-data slope when using SST changes from the MIROC 4m simulation during a strong AMOC phase is similar to the one for the simulation with the Tierney et al. (2020a) SST reconstruction (0.66 and 0.65, respectively).

Impacts of SST boundary conditions on the ∆δ 18 O model-data agreement at polar ice core stations
The modeled values of δ 18 O of snow (δ 18 O sn ) and T 2 m at polar ice cores stations for different boundary conditions in LGM-PI SST changes are compared to isotopic observations and temperature reconstructions in Fig. 5a. Only simulations using sea ice from GLOMAP are selected here. Except for Renland station in the north and Taylor Dome in the south (both coastal sites), ECHAM6-wiso generally underestimates δ 18 O sn changes. We find generally good agreement with the reconstructed temperature data, except for GRIP station where the cooling is 6-10 • C too weak (dots and crosses in Fig. 5a). The substitution of GLAC-1D reconstruction by the PMIP3 one strongly improves the model-data agreement of δ 18 O in Antarctica (Fig. S5), leading to better model-data agreement at the global scale (slope = 0.87, r 2 = 0.62 and RMSE = 3.2 ‰) compared to the LGM_GLOMAP experiment. This agrees with the findings of Werner et al. (2018), who showed that the isotopic model-data correlation for Antarctic ice core stations is weaker when using GLAC-1D instead of thePMIP3 ice sheet reconstruction (RMSE = 2.1 ‰ and 1.1 ‰ for 11 Antarctic stations, respectively). However, this better δ 18 O modeldata agreement is more likely due to a bias compensation than a more realistic ice sheet because the simulation of Antarctic temperatures by ECHAM6-wiso is degraded at the same time (markers in Fig. S5). Except for the Taylor Dome station, all modeled δ 18 O sn values at polar ice core stations are in better agreement with measurements (blue bars in Fig. 5a) when reconstructed SST fields (i.e., GLOMAP or Tierney et al., 2020a) are used (orange and green bars in Fig. 5a, respectively), confirming the results of Werner et al. (2018) about the worse model-data agreement when using sea surface boundary conditions from a coupled model instead of reconstructed ones. The change from one MIROC 4m SST field to another one (i.e., weak or strong AMOC phase) as input for ECHAM6-wiso does not modify the modeled δ 18 O sn values much (red and purple bars in Fig. 5a).
For a stronger SST cooling in the Southern Ocean (GLOMAP and Tierney et al., 2020a), ECHAM6-wiso simulates higher δ 18 O p changes (right maps of Fig. 5) that are in better agreement with the observations. A stronger cooling in the Atlantic sector of the Southern Ocean (GLOMAP SST compared to the Tierney et al., 2020a, one) has the consequence of enhancing the δ 18 O p depletion in the Atlantic-Indian Ocean sector of Antarctica (right map of Fig. 5c) despite similar temperatures (left map of Fig. 5c). This area includes the ice core stations Dome Fuji and EDML, where better δ 18 O and similar temperature model-data agreements are found (orange markers and bars in Fig. 5a). This is also true for other stations further to the east and west, like WDC and EDC (green bars and markers in Fig. 5a). When using an SST field with stronger cooling in the northern North Atlantic Ocean (i.e., GLOMAP), a stronger cooling is simu-  Table 3. Values of δ 18 O model-data slope (1 is better), coefficient of determination r 2 , and root mean square error (RMSE) for our ECHAM6-wiso simulations using different SST and sea ice boundary fields. For each column, worst to best model-data agreements are shown with a yellow-to-green color map. lated there to the south of Greenland. A stronger cooling is simulated in the southern and central parts of inner Greenland too (left maps of Fig. 5) with consequently higher δ 18 O changes between LGM and PI in Greenland (right maps of Fig. 5) except in the northern part like at Camp Century station (Fig. 5a). Better agreement is obtained with the Greenland δ 18 O observations under this configuration (orange bars in Fig. 5a), without improving the overly weak cooling in the model (colored markers for GRIP station in Fig. 5a).

Impacts of sea ice change boundary conditions on
the ∆δ 18 O model-data agreement at polar ice core stations To analyze the effects of sea ice boundary conditions on the modeled δ 18 O changes in polar regions between LGM and PI, we compare the results from the simulations using the same SST (here from the MIROC 4m simulation with the weak AMOC phase) but different sea ice area fraction fields: GLOMAP and MIROC 4m (i.e., LGM_miroc4m_sst_ glomap_ sic and LGM_miroc4m_sst_and_sic simulations, respectively). For all Antarctic ice core stations, a stronger depletion in δ 18 O sn between LGM and PI is simulated when the LGM sea ice in the Southern Ocean is more extensive (GLOMAP, orange bars in Fig. 6a). Except for Taylor Dome, better agreement with isotopic observations is then found as is better agreement with temperature reconstructions in West Antarctica. It has a huge impact on modeled T 2 m anomalies over the Southern Ocean (between 2 and 10 • C), and the simulated cooling is higher by 1 to 4 • C in western Antarctica and in coastal regions of the continent (left map of Fig. 6c). As a consequence, higher LGM-PI anomalies in δ 18 O of precipitation and snow are simulated: more than 5 ‰ over the Southern Ocean and around 1 ‰-2 ‰ on the continent, especially in the western part (right map of Fig. 6c). More extensive northern sea ice (i.e., MIROC 4m) increases the cooling by 5 to 10 • C in the Arctic Ocean and from 0.5 to 5 • C in Greenland (left map of Fig. 6c), giving higher δ 18 O p anomalies of up to 2 ‰ (right map of Fig. 6c) generally in better agreement with the temperature and δ 18 O observations (green bars and markers in Fig. 6a). The opposite is true with the less extensive sea ice distribution from MIROC 4m under a strong AMOC phase, weakening the model-data agreement for Greenland (Fig. S6). A smaller sea ice area fraction in grid cells near coastal areas in MIROC 4m compared to GLOMAP (95 %-98 % against 100 %, respectively. See Sect. 2.2.2), resulting in less important LGM-PI sea ice change, makes the cooling near the Greenland coast slightly lower. This lower sea ice change in MIROC 4m combined with the isotopic content of snow on sea ice taken into account for sublimation processes in sea-ice-covered regions leads to a reduction of the LGM-PI δ 18 O p changes in Baffin Bay (right map of Fig. 6c). This aspect is investigated in detail in Sect. 4.2.

Impacts of LGM AMOC strength on the ∆δ 18 O model-data agreement at polar ice core stations
Here, we investigate the impacts of AMOC strength on the modeled δ 18 O in polar regions. For that, sea surface outputs (i.e., both SST and sea ice spatial distribution) from the MIROC 4m simulation with different LGM AMOC strengths are used as boundary conditions for ECHAM6-wiso. We focus first on the North Pole region because the AMOC strength mainly influences the climate of the Northern Hemisphere, as shown in SST and sea ice distributions used in this study (Figs. 1 and 2). A weaker AMOC during the LGM involves less heat transported in the north and thus lower LGM temperatures (i.e., larger cooling relative to PI), as shown in the left map of Fig. 7c. A difference in T 2 m of up to 10 • C in the North Atlantic and Arctic oceans is seen in the LGM_miroc4m_strong_AMOC_sst_glomap_sic and LGM_miroc4m_ strong_AMOC_sst_and_sic simulations. Cooling in Greenland is reduced by 2 to 6 • C when the AMOC is increased, increasing the warm bias in GRIP (markers in Fig. 7a).
LGM-to-PI changes in δ 18 O in Greenland are mainly controlled by this change in mean temperature with an increase in LGM δ 18 O sn of between 1.2 and 2.5 ‰ at Greenland ice core stations for a stronger LGM AMOC (orange and green bar in Fig. 7a). As ECHAM6-wiso generally underestimates the LGM-PI δ 18 O changes at the poles, a weaker AMOC generally improves the model-data agreement (blue and orange bars in Fig. 7a). In the Southern Ocean and Antarctic regions, only small T 2 m changes are simulated by ECHAM6-wiso due to a change in AMOC strength during LGM (left map of Fig. 7c), as shown by the absence of any effect on model-data temperature agreement (markers in Fig. 7a). As a consequence, modeled δ 18 O sn values are very similar between the two simulations (orange and green bars in Fig. 7a). These small differences are due to the selection of the strong AMOC phase period in the middle of the peak in the MIROC 4m simulation (see Sect. 2.2.1). The impact of the period selection for the strong AMOC phase (e.g., the start or the end of the interstadial) on surface temperature and δ 18 O in Antarctica will be investigated in more detail in a future study. Finally, the SST values alone due to AMOC strength variations change by only less than 1 ‰ of the modeled δ 18 O sn (red and purple bars in Fig. 5a). This shows that the LGM-to-PI changes in sea ice distribution, related to the AMOC strength variations, have a large impact on modeled T 2 m anomalies and consequently on the isotopic signals in the North Pole region.

Impacts of sea surface boundary conditions on δ 18 O-T 2 m temporal slope for LGM-PI climate change
We have analyzed the effects of LGM-to-PI changes in SST and sea ice distribution on modeled T 2 m and δ 18 O of precipitation and snow in the polar regions, as well as the impacts of the LGM AMOC strength. Next, we investigate the repercussions for modeled δ 18 O p -T 2 m temporal slopes and the underlying causes in terms of changes in moisture transport. In other words, are T 2 m and δ 18 O signals in the polar regions influenced by LGM-to-PI changes in SST and sea ice distribution in the same way? What are the underlying dynamics, for example, in terms of changes in concentrations and transport of water vapor? A correction for the prescribed glacial seawater change of 1 ‰ has been applied to LGM δ 18 O values before temporal slope calculation, according to Eq. (1) of Stenni et al. (2010). As in Cauquoin et al. (2019b), the calculation of temporal slopes was restricted to grid cells wherein simulated annual mean temperatures are below +20 • C for both PI and LGM. Moreover, we selected only the grid cells showing an absolute LGM-PI annual mean T 2 m difference of at least of 0. ous modeling studies (Schmidt et al., 2007;Werner et al., 2018;Cauquoin et al., 2019b). For Greenland ice core stations, we find a modeled spatial slope of 0.71 ‰ • C −1 , also in agreement with previous model results (Schmidt et al., 2007;Cauquoin et al., 2019b).

Antarctica
Water vapor δ 18 O at coastal and western low-elevation sites is controlled by nearby local sources, while the evaporative moisture source of high-elevation East Antarctic ice cores is typically further north at around 40-45 • S (Sodemann and Stohl, 2009). So, the values of the δ 18 O p -T 2 m slope in Antarctica are expected to be influenced by sea surface boundary conditions in different ways depending on the ice core site location and elevation. Greater SST cooling in the Southern Ocean influences both the LGM T 2 m and δ 18 O p in the same direction (i.e., toward lower values) but with different magnitudes. Stronger cooling in the eastern part of the ocean (Tierney et al., 2020a, LGM SST is colder than the GLOMAP one) increases the temporal slopes at EDC and Talos Dome by 0.23 ‰ • C −1 and 0.07 ‰ • C −1 , respectively. A higher cooling in the Atlantic sector of the Southern Ocean (GLOMAP SST in Fig. 8a) makes the Antarctic temporal slope values higher between 0 and 90 • E longitude compared to the other simulations (meaning that δ 18 O p is more impacted than temperature). It especially impacts the Dome Fuji and EDML ice core sites, where values of temporal δ 18 O p -T 2 m slopes reach 0.66 ‰ • C −1 and 0.69 ‰ • C −1 , representing an increase of at least 25 % compared to simulations using SST with less cooling. This can be explained by a change in moisture transport. With a stronger SST cooling in the Southern Ocean (GLOMAP against MIROC 4m), the westerlies around Antarctica are enhanced and the atmosphere is wetter in the midlatitudes, while a drier belt appears around the continent (Fig. 9b). More water vapor from the Atlantic sector at 30-40 • S, where the cooling is less strong (compared to MIROC 4m SST) and the water vapor more depleted in δ 18 O (Fig. S4a) because of enhanced evaporation, is transported to eastern Antarctica. As a consequence, the larger cooling in GLOMAP Southern Ocean SST increases the temporal slope at Dome Fuji and EDML stations. A larger SST cooling near the Amundsen Sea (i.e., Tierney et al., 2020a, SST compared to GLOMAP, Fig. 5c) impacts the temperature from this region to western Antarctic sites (2 to 4 • C, left map of Fig. 5c). On the other hand, the δ 18 O of water vapor and precipitation in the Amundsen Sea area is not so impacted by imposed stronger SST cooling (by 2 ‰ at maximum, right map of Fig. 5c). The decrease in the transport of this not-so-depleted water vapor to western Antarctic sites (Fig. S7) increases the temporal slopes by ∼ 0.1 % • C −1 at WDC and Byrd stations (orange marker in Fig. 12). At the same time, this water vapor influences the δ 18 O p of nearby coastal regions like the Antarctic peninsula, decreasing their temporal slopes (Fig. 8b and orange marker in Fig. S9).
LGM-to-PI changes in sea ice area fractions have a strong impact on the slopes in coastal regions, as shown by the comparison between plots (c)-(d) and (e)-(f) in Fig. 8. Law Dome ice core is representative of this impact in coastal areas, with a slope of 0.35 ‰ • C −1 and 0.64 ‰ • C −1 depending on whether less extensive sea ice (LGM_miroc4m_sst_and_sic, Fig. 8d) or more extensive sea ice (LGM_miroc4m_sst_glomap_sic, Fig. 8c) is used, respectively. More extensive sea ice in the Atlantic and Indian sectors of the Southern Ocean changes the transport of vapor only slightly (Fig. 9c). On the other hand, the presence of more sea-ice-covered area close to the continent increases the amount of local water vapor depleted in δ 18 O (compared to the δ 18 O of open seawater) that influences the isotopic composition of snow in coastal sites, increasing their temporal slopes. The more extensive sea ice in the Indian sector of the Southern Ocean increases the slope in a geographical band from Law Dome to Vostok and EDC stations (Fig. 8c  and d). This can be explained by a decrease in transport of enriched water vapor from the Indian Ocean and marine region south of Australia, especially in austral winter (Fig. S9). On average, the modeled temporal δ 18 O p -T 2 m slopes of East Antarctic ice core stations are increased by more than 25 % when more extensive sea ice is used (Fig. 12, dark and light purple markers against the red and grey ones, respectively). This conclusion remains the same when using the average slope across the entire East Antarctic area but with an enhanced difference because the entire coastal area is under consideration (Fig. S9). No clear influence of sea ice extent on temporal slopes at West Antarctica ice core stations is found (Fig. 12). On the other hand, effects on the Antarctic peninsula, Ronne Ice Shelf, and the coast of the Amundsen Sea (Fig. 8), influencing the average slope values in the entire western part of the continent (Fig. S9), are simulated. These are due to changes in the nature of the water source nearby (i.e., open water or sea ice) depending on whether the sea ice used is less extensive or not.

Greenland
Using an SST field with more LGM cooling off the coast of Greenland and in the northern North Atlantic Ocean (GLOMAP) as forcing for ECHAM6-wiso, we model higher δ 18 O p -T 2 m temporal slope values at all Greenland ice core stations (Fig. 10a) compared to all other simulations. The average of the temporal slope values at ice core stations is 0.50 ‰ • C −1 with GLOMAP sea surface boundary forcing (blue marker in Fig. 12) and less than 0.42 ‰ • C −1 in the other ECHAM6-wiso simulations. This difference is enhanced by considering the entire Greenland area (blue marker in Fig. S9). For Greenland, most of the moisture comes from northern North Atlantic Ocean at latitudes 30-40 • N (Drumond et al., 2016), south of the ice sheet (Fig. 11a). A stronger SST cooling in the Arctic-North Atlantic enhances the transport of isotopically depleted water vapor from North America (Fig. 11b) and strengthens the water vapor transport belt around Iceland. Also, ECHAM6-wiso simulates a slight decrease in local water vapor content in the Greenland Sea (with lower LGM δ 18 O because the cooling is stronger there) and its transport to land. It increases the isotope-temperature temporal slope from the sea to the inland area, passing through the eastern Greenland coast where Renland station is located (Fig. 10b and c compared to a).
More extensive sea ice (i.e., MIROC 4m compared to GLOMAP) makes the Arctic Ocean area drier, especially at 50 • N, and it slightly slows down the transport of water vapor from the North Atlantic to Greenland area (Fig. 11c). On the other hand, all this area is covered by sea ice in the "more extensive sea ice" simulation (i.e., MIROC 4m sea ice). It decreases the δ 18 O of water vapor above this surface, increasing the isotope-temperature temporal slope on the eastern Greenland coast and in the Greenland Sea ( Fig. 10d and f compared to c and e, respectively). In the latter, more extensive sea ice, especially in summer, decreases the LGM δ 18 O too, while the effect on temperature is small, again increasing the local temporal slope. It has a slight effect on modeled temporal slopes (∼ 0.1 ‰ • C −1 ) over the eastern coastal regions of Greenland, including Renland station. In ECHAM6wiso, the isotopic composition of sea ice surfaces also reflects the isotopic composition of snow deposited on this surface. Then the formation of new sea ice from PI to LGM acts as a positive feedback effect in the decrease in surrounding δ 18 O p , leading to steeper modeled δ 18 O p -T 2 m temporal slopes (see Sect. S1 and Fig. S10). Finally, ECHAM6-wiso forced with MIROC 4m sea ice, whose fractional areas are artificially lower (i.e., not 100 % sea-ice-covered) in coastal grid cells, simulates lower δ 18 O p -T 2 m temporal slope values over Baffin Bay (between 0.3 ‰ • C −1 and 0.9 ‰ • C −1 , Fig. 9d and f) compared to when the model is forced with GLOMAP sea ice (between 0.6 and 1.25 ‰ • C −1 , Fig. 10c and e). If the MIROC 4m sea ice is corrected by setting the sea ice fraction as 100 ‰ as in GLOMAP (see Sect. S1 and Fig. S11), we obtain temporal slope values similar to those in the simulations forced by GLOMAP sea ice (Fig. S12).
A stronger AMOC increases the amount of water vapor and enhances its transport from the North Atlantic to European coastal areas because of the less extensive sea ice (Fig. 11d). More water vapor with higher δ 18 O is present southeast of Greenland because sea ice is replaced by open water. However, there is only a slight increase in the transport of this water vapor toward north in the Greenland interior (Fig. 11d), while the cooling inland is largely reduced (Fig. 7c). So, isotope-temperature temporal slopes are slightly increased over the interior of Greenland for stronger  On the other hand, temporal slopes are decreased over the Greenland Sea (Fig. 10f) because of the presence of open water instead of sea ice, enhancing the LGM δ 18 O p locally. Note that stronger AMOC worsens the model-data agreement for both T 2 m and δ 18 O (Fig. 7a).

Conclusions and perspectives
In this study, we elevated the importance of sea surface boundary conditions for the relationship between nearsurface air temperature and δ 18 O p for LGM-to-PI climate change. Figure 13 illustrates the main findings of our study. In Antarctica, the temporal slopes in coastal areas like the Antarctic peninsula, the Ronne Ice Shelf, near the Amundsen Sea, and the eastern coast are influenced by sea ice extent. The more extensive the sea ice, the steeper the slope in these areas (including Law Dome). The transport of mois-ture to Antarctica is generally only slightly changed with variations in the sea ice extent. On the other hand, the distribution of open water and sea-ice-covered areas around the continent mainly controls the nearby δ 18 O p and therefore the isotope-temperature temporal slopes. We found that temporal slopes further inland, at the EDC and Vostok sites, are influenced by changes in sea ice extent through a weakening of moisture transport from the Indian Ocean and marine region south of Australia when sea ice is more extensive (left map Fig. 13). The values of δ 18 O p -T 2 m temporal slopes in inland East Antarctic ice core stations like Dome F, EDML, EDC, and Talos are mainly controlled by the change in SST in our ECHAM6-wiso simulations. Stronger cooling in the Atlantic sector of the Southern Ocean (GLOMAP) leads to steeper temporal slopes at Dome F and EDML due to enhanced water vapor transport from lower latitudes. Strong cooling in the Amundsen Sea weakens the transport of relatively less depleted water vapor (compared to the large cooling) in inland West Antarctica. It slightly increases the temporal slopes at the WDC and Byrd sites. At the same time, this water vapor contributes to the nearby coastal region, decreasing the temporal slopes there (left map of Fig. 13). These various effects on the temporal slopes can complicate the interpretation of δ 18 O p signals in West Antarctica, depending on the sites considered. Our modeled results clearly demonstrate that the δ 18 O p -T 2 m temporal slopes in Greenland are influenced by the sea surface temperatures very near the coasts. The greater the LGM cooling off the coast of southeastern Greenland, the larger the δ 18 O p -T 2 m temporal slopes (right map of Fig. 13) because of an enhancement in the transport of isotopically depleted moisture from the western North Atlantic Ocean. Similarly, the presence or absence of sea ice very near the coast can impact the modeled temporal slopes at some Greenland ice core stations. It has a large influence in Baffin Bay and the Greenland Sea, too. A more or less large  southern extent of the sea ice (MIROC 4m in weak AMOC phase or GLOMAP, respectively) has only a slight impact on modeled temporal slope values (light purple and red markers in Fig. 12, respectively). The seasonal variation of the sea ice distribution, especially its retreat in summer, changes the origin type of the water source (from evaporation of open water or sublimation of snow on sea ice) and thereby influences the δ 18 O of local vapor contributing to the δ 18 O p -T 2 m temporal slopes in the eastern part of inner Greenland (right map of Fig. 11). Finally, a stronger LGM AMOC slightly impacts the temporal slopes modeled by ECHAM6-wiso over Greenland, mainly because of the changes in the sea ice dis-tribution. Higher temporal slopes are modeled by ECHAM6wiso, mainly because the cooling in Greenland is largely reduced while changes in δ 18 O anomalies remain localized in the Greenland Sea. For Antarctica, only small changes in surface temperature and δ 18 O are modeled by ECHAM6-wiso because the strong phase period was selected from within the middle of the AMOC peak. The impact of the precise period chosen to represent the strong AMOC phase, for example the start or the end of the interstadial, will be investigated in more detail in a future study.
In Greenland, ECHAM6-wiso simulates δ 18 O p -T 2 m temporal slopes oscillating between 0.2 ‰ • C −1 and 0.7 ‰ • C −1 inland and at northwestern coastal sites, respectively, which are lower than the spatial ones (0.71 ‰ • C −1 , Fig. 12), as already reported in previous studies (Buizert et al., 2014;Cauquoin et al., 2019b;Jouzel, 1999;Werner et al., 2000). Our modeled temporal slope values for stations NEEM (around 0.7 ‰ • C −1 ) and NGRIP (between 0.37 ‰ • C −1 and 0.57 ‰ • C −1 ) are in agreement with previous reconstructions (Buizert et al., 2014), too. In Antarctica, the ECHAM6wiso modeled δ 18 O p -T 2 m temporal slopes for LGM-to-PI climate change are on average lower than the PI spatial slopes of the same model by at least 0.20 ‰ • C −1 and 0.48 ‰ • C −1 for eastern and western ice core locations, respectively (Fig. 12), regardless of the simulation being considered. By extension, we found much lower temporal slope values than the ones estimated by Buizert et al. (2021). We simulate a maximum temporal slope value of 0.9 ‰ • C −1 for the South Pole, while Buizert et al. (2021) found temporal slopes in Antarctic ice core stations ranging from 0.9 to 1.4 ‰ • C −1 , which are higher than the observed spatial δ 18 O p -T 2 m slope of 0.8 ‰ • C −1 (Masson-Delmotte et al., 2008). The use of the thicker PMIP3 ice sheet reconstruction compared to GLAC-1D increases the resulting modeled δ 18 O p -T 2 m temporal slopes in ECHAM6-wiso ( Fig. S13) with mean values for East and West Antarctic ice core stations equal to 0.68 ‰ • C −1 and 0.92 ‰ • C −1 , respectively, by decreasing the isotopically enriched bias in the model for LGM (Fig. S5). However, the temperature model-data agreement is reduced in this case. Even if ECHAM6-wiso shows biases in the representation of the isotopic content of Antarctic precipitation, we emphasize that the purpose of our study was to investigate the relative effects of sea surface conditions and AMOC strength on the links between δ 18 O p and near-surface air temperature, regardless of agreement or disagreement with other slope reconstructions.
In addition to effects of orography, fractionation during the sublimation of surface ice is not taken into account in ECHAM6-wiso as in many isotope-enabled AGCMs. This process would lead to a further decrease in the δ 18 O of water vapor in the polar regions, contributing to steeper modeled δ 18 O p -T 2 m temporal slopes in regions with low temperature. The mismatch between our model slopes and the reconstructed ones from Buizert et al. (2021) could be related to the representation of the atmospheric boundary layer and the related inversion temperature (Krinner et al., 1997;Masson-Delmotte et al., 2006;Cauquoin et al., 2019a), too. Still, despite these biases that potentially affect our modeled δ 18 O p -T 2 m temporal slopes for LGM-to-PI climate change, our ensemble of simulations provides information on how sea surface conditions partially control the links between δ 18 O p and near-surface air temperature in the polar regions through changes in the sources or transport of moisture arriving at the poles.
Because only ECHAM6-wiso is used in this study, we cannot exclude the model dependency of our results. So, the use of isotope-enabled AGCMs other than ECHAM6-wiso would be beneficial to confirm or refute our findings. A set of SST reconstructions for the LGM, based on both model re-sults and observations, is now available. We stress the importance of providing sea ice cover reconstruction, which contributes not only to the δ 18 O p of coastal sites but also of some inland ice core stations, for this period too. The sea ice cover simulated by coupled GCMs for the LGM period takes various forms. An alternative reconstruction to the GLOMAP one, also based on observations, would help to better assess the impact of sea ice cover on the δ 18 O p -T 2 m relationship for LGM-to-PI climate change. We also emphasize that more proxy measurements of temperature and sea ice are necessary for the Southern Ocean. Relatively large uncertainties remain in the reconstruction of the climatology in this area, while the water vapor from this region contributes largely to δ 18 O p in Antarctica. Finally, by showing the sensitivity of δ 18 O p -T 2 m temporal slopes to sea surface boundary conditions, the potential uncertainties of the latter could have an impact on the reconstruction of the former (Jouzel et al., 1997(Jouzel et al., , 2003Markle and Steig, 2022).
As a first step, the focus of this study was to identify and quantify the important factors influencing the isotopetemperature relationship in the polar areas for the LGM-to-PI climate change. Future studies will investigate the evolution of this relationship along the whole of the last deglaciation. For that, an ensemble of equilibrium isotopic simulations using the sea surface and ice sheet boundary conditions from MIROC 4m for different time periods between the LGM and PI will be performed.
Code availability. The ECHAM model code is available under a version of the MPI-M software license agreement (https:// www.mpimet.mpg.de/en/science/models/license/, Mauritsen et al., 2019). The code of the isotopic version ECHAM6-wiso is available upon request on the AWI's GitLab repository (https://gitlab. awi.de/mwerner/mpi-esm-wiso, Cauquoin et al., 2019b;Cauquoin and Werner, 2021  Competing interests. At least one of the (co-)authors is a member of the editorial board of Climate of the Past. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement.
This article is part of the special issue "Ice core science at the three poles (CP/TC inter-journal SI)". It is not associated with a conference.