Deglacial climate changes as forced by different ice sheet reconstructions

. During the Last Deglaciation, the climate evolves from a cold state at the Last Glacial Maximum at 21 ka with large 15 ice sheets, to the warm Holocene at ~9 ka with reduced ice sheets. The deglacial ice sheet melt can impact the climate through multiple ways: changes of topography and albedo, bathymetry and coastlines, as well as fresh water fluxes. In the PMIP4 protocol for deglacial simulations, these changes can be accounted or not depending on the modelling group choices. In addition, two ice sheet reconstructions are available (ICE-6G_C and GLAC-1D). In this study, we evaluate all these effects related to ice sheet changes on the climate using the iLOVECLIM model of intermediate complexity. We show that the two 20 reconstructions yield the same warming to a first order, but with a different amplitude (global-mean temperature of 3.9°C with ICE-6G_C and 3.8°C with GLAC-1D) and evolution. We obtain a stalling of temperature rise during the Antarctic Cold Reversal (from ~14 ka to ~12 ka) similar to proxy data only with the GLAC-1D


Introduction
The Last Deglaciation is a time of large climate transition from the cold Last Glacial Maximum (LGM) with large ice sheets at ~21 thousand years ago (ka), to the warmer Holocene with reduced ice sheets at ~9 ka (Clarke et al., 2012). During that time, the global mean surface temperature likely rises by 4.5±0.9°C (Annan et al., 2022), and possibly 6.1±0.4°C (Tierney et al., 2020). Deglacial surface changes are neither smooth nor uniform. Rapid regional changes occur at different times 35 depending on location. In particular, the Northern and Southern high latitudes display a markedly different behaviour. The temperature rise in the Southern Hemisphere starts earlier than in the Northern Hemisphere on average (Shakun et al., 2012).
In Antarctica, the warming then halts between ~14.8 ka and 13 ka at the time of the Antarctic cold reversal (ACR, Jouzel et al., 2007). In Greenland, the temperature abruptly increases at ~14.7 ka at the beginning of the Bolling-Allerod, then drops from 12.8 ka to 11.6 ka during the Younger Dryas (Buizert et al., 2014). As the ice sheets shrink, especially in the Northern 40 Hemisphere, large fluxes of fresh water go to the ocean and lead to a considerable increase of sea level by around 120-130 m (Lambeck et al., 2014;Gowan et al., 2021). Besides temperature, proxy records indicate large changes of ocean circulation during the transition (McManus et al., 2004;Ng et al., 2018). They suggest a strong weakening of the A tlantic Meridional Overturning Circulation (AMOC) at the beginning of the deglaciation from ~19 ka followed by a resumption at ~15 ka. Ocean circulation changes result in a modification of heat transport and variations of temperature in the Northern and So uthern 45 Hemispheres. Sea ice changes are also likely to play a crucial role in oceanic circulation changes, because convection and deep-water formation are linked to sea ice formation. The sea ice formation rejects salt and thus makes surrounding waters saltier and denser. However, the links between ice sheets, ocean circulation and inter-hemispheric temperature evolution during the transition are difficult to disentangle. Ice-sheets, ocean circulation, atmospheric and oceanic temperature are thus linked, but disentangling the timing and causal links between all of their changes during the transition is difficult. 50 To better understand the drivers of the deglaciation and the response of the different components of the Earth system, modell ing groups have endeavoured to run transient simulations from the Last Glacial Maximum to the Holocene. To do so, models need to account for changes in orbital forcing greenhouse gases and ice sheets. Previous studies have focused on the evolution of climate in simulations with prescribed ice sheets from reconstructions (e. g. Menviel et al., 2011;He et al., 2013;Kapsch et al., 2022) or with interactive ice sheets (e. g. Bonelli et al., 2009;Ganopolski and Calov, 2010;Quiquet et al., 2021), which is 3 technically more challenging. The warming and ice sheet melt are largely driven by the insolation change due to the evolution of orbital parameters, especially the insolation increase in the high Northern latitudes in summer (Berger et al., 1978) as w ell as the atmospheric CO2 concentration increase (Bereiter et al., 2015;Gregoire et al., 2015).
The Last Deglaciation is an ideal test bed to evaluate models during a period of large warming and to better understand the underlying processes, as many proxy data exist for that period. As such, more and more models have undergone simulating 60 the transition and a new PMIP4 protocol has been set up for the Last Deglaciation (Ivanovic et al., 2016) to facilitate the comparison of model results in a common framework. In the protocol, the mai n prescribed boundary conditions are orbital parameters, atmospheric greenhouse gases and ice sheets. While the insolation and greenhouse gas concentration evolution are relatively well known, the ice sheet evolution is more uncertain, hence two different reconstructions are included in the PMIP4 protocol: ICE-6G_C and GLAC-1D ( Figure 1, Argus et al., 2014;Peltier et al., 2015;Briggs et al., 2014;Tarasov et al., 2012;65 Tarasov and Peltier 2002;Gowan et al., 2021). These reconstructions are obtained from inverse modelling based on GPS data and local sea level records. Large uncertainties remain, mostly due to the viscosity model use for solid Earth. Besides the ice sheet reconstruction itself, the PMIP4 protocol lets to some extent modelling group decide how to implement ice sheet related changes such as land-sea mask (i.e. coastlines), bathymetry and fresh water fluxes. The purpose of this 75 flexibility in the protocol is to allow any willing modelling group to join the effort, whatever the complexity of their model, from intermediate complexity to state-of-the art General Circulation Models (GCMs). However, not all models are able to run the entire range of possible simulations, it is thus crucial to evaluate the impact on climate of the different choices in setting up the experiments. More specifically, model groups have three choices related to ice sheets: (a) the choice of ice sheet reconstruction: ICE-6G_C or GLAC-1D; (b) to only account for topography/ice mask (i.e. albedo changes), or also land-sea 80 mask and bathymetry changes and (c) to account or not for the flux of meltwater from ice sheets to the ocean, and if included to do so in a uniform way or a more realistic way such as by routing melted water towards the ocean. The relative effect of all these individual modelling choices on the transient climate simulated is unclear has not been evaluated.
The effect of ice sheet related changes has been analysed to some extent in past studies, but not all of them, usually not within 85 the PMIP4 deglacial framework, and not in a systematic way. For example, an evolving bathymetry (including changes of coastlines) has recently been implemented in a GCM (Meccia and Mikolajewicz, 2018), and while it has been used in transient runs of the Last Deglaciation (Kapsch et al., 2022), the impact of evolving bathymetry compared to fixed bathymetry on climate has not been evaluated yet. Concerning fresh water fluxes, they were usually not included in early studies of deglacial changes, (Timm and Timmermann, 2007;Roche et al., 2011;Smith and Gregory, 2012), possibly explaining the lack of rapid 90 regional changes in these simulations. In contrast, in simulations with fresh water fluxes, rapid events were observed in addition to the global deglacial warming (Liu et al., 2009;Menviel et al., 2011;Bethke et al., 2012;He et al., 2013;Obase and Abe -Ouchi, 2019). However, the simulations that show temperature and circulation changes in best agreement with data were usually run with tuned fresh-water fluxes (Liu et al., 2009;Menviel et al., 2011;He et al., 2013;Obase and Abe-Ouchi, 2019).
On the one hand, such prescribed fresh water fluxes are usually inconsistent with the volume change inferred from ice sheet 95 reconstructions. On the other hand, simulations run with the fresh water flux evolution derived from ice sheet reconstructions show disagreement with observed climate changes (e.g. Bethke et al., 2012;Kapsch et al., 2022).
Here we evaluate the impact of ice sheet led changes, i.e. topography and albedo, bathymetry and coastlines, as well as fresh water fluxes, on climate, especially temperature and ocean circulation. For this, we use iLOVECLIM, a climate model of 100 intermediate complexity, fast enough to run several simulations and test these different effects. This work will help other modelling groups decide which choices to make in terms of ice sheet related changes.

iLOVECLIM model 105
We use the iLOVECLIM model, which is a code fork and evolution of the LOVECLIM model (Goosse et al., 2010).
iLOVECLIM is an intermediate complexity model which shares the same atmosphere, ocean, sea ice and terrestrial biosphere components as LOVECLIM. The ocean grid has a 3°x3° resolution with 20 irregular vertical levels while the atmosphere is on a T21 grid with three vertical levels. It can simulate ~1000 years / day, making it very suitable for computing long term climate changes such as the work undertaken here. 110 To run simulations of the last deglaciation following the PMIP4 protocol, we have added modifications to account for the evolution of the ice sheets as reconstructed in ICE-6G_C and GLAC-1D. We consider the change of ice sheet topography and ice mask, the change of bathymetry and coastlines, as well as the routing of fresh water from ice sheet melt to the ocean.
To run simulations of the last deglaciation following the PMIP4 protocol, we have modified the code to be able to change the bathymetry and coastline regularly, and we have added the routing of fresh water from ice sheet melt to the ocean. We have 115 also generated the ice sheet topography, ice sheet mask, bathymetry and land-sea mask files on the iLOVECLIM grid using the ICE-6G_C or GLAC-1D reconstructions. This is detailed in the sections below.

Modification of ice sheet topography and ice mask
The ice sheet geometry change is accounted for in the atmospheric component of the model through two variables: the surface 120 elevation (topography) and the ice sheet mask (albedo). These two variables are updated interactively in the course of the simulation to follow the ICE-6G_C or GLAC-1D reconstructions, using anomalies with respect to ETOPO1 (NOAA, 2009).
The updates are done abruptly, i.e. without temporal interpolation, every 500 years for the ICE-6G_C reconstruction and every 100 years for the GLAC-1D reconstruction. The difference in frequency update is due to the difference in the frequency of available ice sheet reconstructions (500 years for ICE-6G_C and every 100 years for GLAC-1D). The ice sheet albedo in the 125 iLOVECLIM model is set to the constant and homogeneous value of 0.85. The vegetation is dynamically computed in the model by the terrestrial biosphere model (VECODE; Brovkin et al., 1997).

Modification of bathymetry and coastlines 130
We have extended the computation of bathymetry developed in Lhardy et al. (2021) for the deglaciation. Following the same procedure, a new bathymetry has been adapted to the iLOVECLIM ocean grid based on those provided in the ICE-6G_C and GLAC-1D reconstructions at regular intervals starting at 21ka. The bathymetry anomaly (with respect to the pre-industrial) at each time step is added to the pre-industrial bathymetry based on ETOPO1. The bathymetry associated with each ice sheet reconstruction is interpolated into the ocean grid in an automated way, with a few limited manual adjustments in straits or key 135 passages. We decide to maintain passages open at our grid resolution when they were open in the reconstructi ons. Since the ice sheet reconstruction (and thus bathymetry) is provided every 100 years for the GLAC-1D reconstruction and every 500 years for the ICE-6G_C reconstruction, we keep the same frequency for the bathymetry update in iLOVECLIM. To test the effect of the frequency update of bathymetry and land-sea mask, we have also run a simulation with GLAC-1D, updating the files every 500 years instead of 100 years. With the reconstructed bathymetry, the ocean volume increases across the 140 deglaciation due to the ice sheet volume decrease ( Figure 2). The initial volume of the ocean at the LGM is lower with ICE-6G_C (1.291 10 18 m 3 ) than with GLAC-1D (1.296 10 18 m 3 ) as more ice is stored in the ice sheets reconstructed with ICE-6G_C. Apart from this relatively constant shift, the evolution is rather similar in the two reconstructions until 14.5ka. From around 14.5ka, the ocean volume increase accelerates, and this acceleration is larger with ICE-6G_C so that the shift between both reconstructions is reduced. From around 12ka and as we get nearer to the Pre-Industrial, the ocean volume becomes 145 roughly the same in both reconstructions: as the ice sheet topography is prescribed as an anomaly with respect to ETOPO1 in the model, when the additional ice sheet volume disappears both ice sheet forcings converge.
There is more variability in the GLAC-1D evolution due to the update of bathymetry every 100 years compared to every 500 years for ICE-6G. On top of that, from 12.5 ka the evolution is more variable in GLAC-1D than ICE-6G, with accelerations and slow-downs. 150

Figure 2. Evolution of a) the ocean volume (m 3 ) with the ICE-6G_C and GLAC-1D bathymetry reconstructions and b) the surface area (10 6 km 2 ). The dashed lines indicate the values for the simulations with fixed bathymetry.
With the evolving bathymetry, we generate evolving land-sea masks (Figure 3). At the LGM, the continental area is larger than at the Pre-Industrial due to the lower sea level (Figure 2b). It then decreases throughout the deglaciation while the sea level rises and new ocean grids appear. At the LGM, the new land-sea mask implies a closing of the Bering strait and Canadian archipelago. The strait of Gibraltar is maintained open throughout the simulation. Every time the bathymetry is updated in the model, conservation of salt is ensured by redistributing the excess (or loss) of salt 170 in the global ocean. In the new oceanic cells, the variables (such as temperature) are initialised using the mean val ue of neighbour cells. When new continent grids appear, the vegetation model computes a new albedo given the vegetation distribution computed interactively by the model.

Routing of melted water from ice sheets 175
In iLOVECLIM, contrary to the original LOVECLIM model, the routing is no longer done in fixed routing basins. Rather, water routing is computed along the greater slope of the model topography and is updated accordingly whenever the topography is updated. This allows for a portable evolution of the first order changes in river routing along the deglaciation.
The routing thus computed is used for the precipitation (interactive) and for the anomalous freshwater flux arising from the ice-sheet melt. The latter is computed as the change in ice-sheet thickness between two reconstruction snapshots, averaged as 180 an annual flux and applied homogeneously over the year. The water fluxes are added as a water flux into the ocean model hence changing the salinity. With the new routing scheme, the fresh water flux from ice sheet melt is thus not added homogeneously in the ocean but routed towards the closest ocean grid cell following the topography.

Simulations following the PMIP4 protocol 185
All transient simulations are forced with evolving orbital parameters (Berger et al., 1978), as well as the increased concentration of greenhouse gases, notably CO2 (Bereiter et al., 2015), CH4 (Loulergue et al., 2008) and N2O (Schilt et al., 2010), following the PMIP4 deglacial protocol (Ivanovic et al., 2016). To evaluate the impact of ice sheet led changes we have run a set of simulations (Table 1)

Impact of different ice sheet reconstructions (topography and albedo)
The effects on climate of the two ice sheet reconstructions (ICE-6G_C and GLAC-1D) are first compared when only considering the changes of topography and albedo (no bathymetry nor fresh water changes). These simulations will serve as a 205 reference for the other ones with the added effects of bathymetry and fresh water fluxes.
In both simulations, the global mean temperature is ~4°C colder at the LGM compared to the pre-industrial (Figure 4), in the range of the PMIP4 models temperature change (3.3° to 7.2°C; Kageyama et al., 2021) and reconstructions (Annan et al., 2022). The global mean temperature is slightly warmer in the simulation with the GLAC-1D reconstruction compared to the simulation with ICE-6G_C due to slightly smaller and lower ice sheets with GLAC-1D (Ivanovic et al., 2018; Figure 1). This 210 temperature difference is observed until ~11.5 ka, when temperatures from both simulations have reach similar values, due to the ice sheet volume evolution. Indeed, the ice sheet height and extent become similar for both reconstructions from ~11.5 ka.
The global mean temperature evolution is relatively similar to proxy data reconstruction (Shakun et al., 2012), with a progressive warming until ~14ka. At this time, while the data show a slow-down and then a decrease in temperature, the temperature with ICE-6G continues to increase while the temperature with GLAC-1D shows a stalling more similar to data, 215 but at a later date. This can be linked to the Northern hemisphere ice sheet reconstructions: in GLAC-1D, following a strong ice volume decrease until 14 ka, the volume stays constant from 14 ka to 12 ka, while it continues to decrease, although at a lower rate, with ICE-6G_C. This can be linked to the Northern hemisphere ice sheet reconstructions. In GLAC-1D the ice volume decreases strongly until 14ka, then stays constant from 14 ka to 12 ka. With ICE-6G_C it continues to decrease, although at a lower rate. From ~11 ka, both data and model simulations show an increase in temperature followed by a plateau 220 at the beginning of the Holocene.

230
The sea surface temperature distribution at the LGM and the warming pattern at the beginning of the Holocene at ~8 ka compared to the LGM are also similar in the two simulations (Figure 5a). Most of the warming at 10ka has occurred in the North Atlantic and Southern Ocean, where deep water forms and penetrates in the ocean interior. In summary, the first order change of temperature is similar in the two simulations with the two reconstructions, and in 240 relatively good agreement with global mean data. The main differences are a temperature shift from 21 ka to 12ka, and a difference in the evolution between 14 ka and 12 ka, at the time of the Antarctic Cold Reversal and Younger Dryas.

Impact of evolving bathymetry and land-sea mask (vs fixed ones)
On top of the changes of topography and albedo, we have run simulations where we also modify the bathymetry and land-sea 245 mask (later referred to as "with bathymetry") for the two ice sheet reconstructions. The evolving bathymetry and land -sea mask lead to an increase of the ocean volume (Figure 2), and hence a decrease of the global salinity as the global salt content is conserved (Figure 6). The salinity change at the LGM that is computed directly from the volume change (1.4 psu for ICE-6G_C and 1.3 psu for GLAC-1D)) is larger than the one prescribed in the PMIP4 standard protocol without bathymetry change where the bathymetry is fixed to the pre-industrial one (1 psu). 250

255
Accounting for the progressive change of bathymetry over time due to the ice sheet changes has a limited impact on the global mean temperature: its evolution is relatively similar with or without accounting for bathymetry changes (Figure 7a, b) apart from small changes, such as between 13 ka and 12 ka with ICE-6G, and a constant shift towards slightly colder temperatures from ~12 ka onwards with ICE-6G_C and from ~14 ka onwards with GLAC-1D. The shift of temperature towards colder values is also visible in the temperature evolution in Greenland (Figure 7c and d) and Antarctica (Figure 7e and f), with a larger amplitude than at global scale. The relative cooling in the simulations with bathymetry changes compared to the simulations with fixed bathymetry is higher locally than at a global scale, especially in 270 (1.1°C with ICE-6G_C and 0.7°C with GLAC-1D) compared to ~0.4°C globally (0.5°C with ICE-6G_C and 0.3°C with GLAC-1D). With both reconstructions, this shift is visible mostly from ~14 ka onwards. With ICE-6G_C, Figure 7e also shows a rapid increase then decrease of temperature between 14 ka and 12 ka in Antarctica. This is due to a change in the bathymetry and sea-land mask around Antarctica, where a grid cell switches from land to ocean. Forcing this grid cell to remain as land 275 suppresses this unrealistic response (Figure 7 simulation "ICE-6G_C evolving bathy mod"). In these simulations, there is no large abrupt climate change in Greenland contrary to data, as the AMOC strength remains strong throughout the simulation in the absence of fresh water flux (see later discussion and figure 14).
The shift towards colder temperature with evolving bathymetry is mainly due to the albedo change associated with the landsea mask modification. As shown on Figure 8, except for a few time periods, the albedo becomes higher in all simulations 280 with evolving bathymetry compared to the corresponding simulations with fixed bathymetry, a difference that becomes more visible from ~15-14 ka onward, at the same time as the shift in temperature observed on Figure 7. In these simulations, the continental surface area decreases with time while the ocean expands due to the ice sheet volume decrease. While the ocean volume area starts changing earlier, the change in surface area only starts to be significant from ~15 ka (Figure 2). At high latitudes where it is cold enough to have sea ice, the continental surface whose albedo can vary from ~0.1 without snow to 285 ~0.8 with snow is replaced by surfaces with sea ice, which has a high albedo of ~0.8 (which can decrease depending on the sea ice type, down to 0.1 for very thin ice), leading to the cold temperatures.

290
In agreement with the albedo change, the sea ice area is also significantly impacted by the evolving bathymetry (Figure 9), because the ocean surface area, especially in the North Atlantic and Arctic, but also in the Southern Ocean, is very differen t when accounting for bathymetry changes. Due to the sea level increase from ice sheet melt, the land-sea mask is modified, and the ocean surface in the high latitudes of the North Hemisphere increases with time, leaving more surface for sea ice to form 295 where it is cold enough. From 15 ka onwards, this leads to more sea ice in the simulation with evolving bathymetry compared to the simulation with fixed bathymetry. With fixed bathymetry, the ocean surface is kept the same as the LGM one throughout the simulation (Figure 2b), and the primary effect governing sea ice changes is the warming resulting in less area covered by sea ice over time. With evolving bathymetry, this effect is counteracted by the increase of ocean surface at high latitudes so that the sea ice area shows a more limited evolution reduction. 300 Figure 9. Evolution of sea ice area (10 12 km 2 ) in the Northern and Southern Hemispheres.

305
In the Southern Ocean, the land-sea mask changes less (Figure 3), and later than in the North Atlantic, yielding apparent changes between simulations with and without bathymetry at a later date than in the Northern Hemisphere. The difference between simulations with and without bathymetry starts to be visible from 14 ka with GLAC-1D, and from 13 ka with ICE-6G. As can be seen on Figure 1, the Antarctic ice sheet starts melting earlier with GLAC-1D compared to ICE-6G, we thus observe a stronger and earlier change in sea ice with the two GLAC-1D simulations than the ICE-6G simulations. With both 310 reconstructions, apart from between 13 ka and 12 ka with ICE-6G, having an evolving bathymetry leads to a larger sea ice area compared to the simulation with fixed bathymetry, as the ocean surface where sea ice can form increases with evolving bathymetry.
The sea ice formation is tightly linked to convection and AMOC changes: the sea ice forms where it is cold enough for water 315 to freeze, which corresponds to cold places where the water becomes denser. In the North Atlantic, deep convection occurs where the heat loss from the ocean to the atmosphere is large enough for the salty and warm water coming from the south to significantly cool and make the water dense enough to convect. The dominant effect of sea ice is to isolate the ocean from the atmosphere, preventing heat loss and hence convection, so that convection takes place at the sea ice edge, where it is very c old and heat loss can take place. In the Southern Ocean, the dominant effect of sea ice is the release of brines (very salty water) 320 during sea ice formation. This makes the water denser and favours convection.
At the LGM, due to the colder climate, the simulated Northern Hemisphere sea ice edge is shifted to the south with respect to the pre-industrial with both ice sheet reconstructions (Figure 10a and b), and the mixed layer depth is deep, showing strong convection in the Iceland, Norwegian and Irminger Basins. At 10 ka, with the warmer conditions the sea ice edge retreats 325 toward the north and the mixed layer depth is reduced in all simulations ( Figure 10). However, the convection sites differ between simulations. In particular, while there is a convection site in the Labrador Sea in the ICE-6G simulations with fixed bathymetry, this convection site is not visible anymore in the simulation with evolving bathymetry, where the sea ice edge is shifted to the south. This is due to the colder conditions in the simulation with evolving bathymetry. The sea ice cover isolates the ocean from the atmosphere preventing cooling and deep convection from taking place in the Labrador Sea. Instead, another 330 convection site is observed south of Greenland. This difference in convection sites observed in the model (Renssen et al., 2005) could explain the strong change in AMOC strength in ICE-6G_C (Figure 11a): the AMOC strength decreases from 12 ka when accounting for bathymetry change while 340 it stays relatively constant with a fixed bathymetry, leading to a ~10 Sv difference between the two simulations. This is not the case with GLAC-1D where only small changes between the two simulations can be seen. This can also be seen in the meridional streamfunction shown on Figure 12. With evolving bathymetry and ICE-6G_C, the AMOC is weaker and shallower at 10 ka compared to the simulation with fixed bathymetry (with ICE-6G). In the GLAC-1D simulations, the sea ice edge and convection sites are also shifted southward in the simulation with evolving bathymetry compared to the one with fixed 345 bathymetry (Figure 10d and f), but the latitude change is smaller, and not sufficient to shut down the Labrador Sea convectio n site, hence it has a limited effect compared to the ICE-6G_C simulations on the AMOC evolution (Figures 11 and 12).  The update frequency for bathymetry and land-sea mask is different for the two reconstructions: 500 years for ICE-6G_C and 100 years for GLAC-1D. To test the impact of different update frequency, we have run an additional simulation with the 365 GLAC-1D reconstruction, with an update frequency of 500 years, similar to ICE-6G_C. As shown on Figure 13, the frequency update modifies the global temperature evolution. It often results in a delayed response compared to the simulation with 100 year frequency, as the bathymetry and land-sea mask change happens later, for example at ~19 ka, ~16.5 ka, ~11 ka or ~10 ka.
Yet the effect is limited, and the difference between the two simulations with the same reconstruction but different update frequency is smaller than the difference between the simulations with the two reconstructions and the same update frequency. 370

375
In summary, accounting for bathymetry changes has limited impact on the global mean temperature change, but can lead to large changes in regional climate, such as the convection sites, sea ice, AMOC and temperature in the North Atlantic. This depends on the reconstruction, the model grid resolution, and probably also on the model sensitivity.

Impact of fresh water flux from me lting ice sheets
In addition to topography and bathymetry, changes of ice sheets can also have an impact on ocean circulation and climate through the input of fresh water coming from the melting of the ice. We now consider this impact by accounting for fresh water being routed from the site of ice sheet melting to the nearest ocean grid cell. The total water flux from ice sheet melting is relatively small (less than 0.1 Sv) until around 15 ka with both reconstructions (Figure 13a). At 15 ka, the flux becomes very 385 large, up to more than 0.5 Sv, corresponding to meltwater pulse 1A (Deschamps et al., 2012). After 14 ka, the fluxes from the two simulations differ more, with the flux from GLAC-1D showing more variability (due partly to the higher frequency update).

395
The fresh water input has a direct effect on convection as it reduces the surface salinity in convection zones, leading to re duced AMOC strength. This is the case with the FWF simulation showing a reduction of AMOC strength starting at around 15 ka (when the fresh water flux becomes very large), and a shut down from around 14 ka (Figure 14). In this simulation, the AMOC does not recover and the simulation crashes at year 11792 BP, i.e. after 9208 years of simulations, for an unknown reason. The 400 response to a fresh water flux is very model dependent (Gottschalk et al., 2019). Hence, to emulate the response of model les s sensitive to fresh water inputs, we have reduced the flux by a factor of 3 and 4. With the reduction of a factor 3 (FWF/3), the AMOC still collapses and recovers only around 8 ka. With a larger factor of 4 (FWF/4), the AMOC is weakened but quickly recovers. Similar changes of AMOC have been observed in an iLOVECLIM version coupled to an ice sheet model . As in other modelling studies using realistic fresh water fluxes computed from ice sheet melt (Kapsch et al., 405 2022), the timing of the AMOC weakening disagrees with data (McManus et al., 2004;Ng et al., 2018), taking place too late at around 15 ka when the fresh water flux becomes large, while proxy data indicate an AMOC weakening from ~17 ka, when the fresh water flux is too small to have an impact on the simulated AMOC. Because of the AMOC reduction and shutdown, less heat is brought to the North Hemisphere and the temperatures in Greenland are strongly reduced in the FWF simulation compared to the standard simulation without FWF (Figure 15b). This 415 is also the case when the water flux is divided by 3. With the fresh water flux divided by 4 (FWF/4) the AMOC strength is only slightly diminished, and the temperature decrease in Greenland is very limited, both in amplitude, and in time as the AMOC recovers quickly. In this case, we observe two AMOC drops (Figure 14, similar to those observed in Kapsch et al., 2022), and simultaneously two temperature drops in the Greenland temperature evolution (Figure 15b).

425
While the Bolling-Allerod warming in the Greenland record at ~14.5 ka cannot be simulated in GLAC-1D FWF as the fresh water flux from ice sheet melting becoming substantial only later at ~15 ka, the temperature evolution in Antarctica and in t he global mean (Figure 15a and c) is in better agreement in the model compared to data, with a stalling in the temperature increase at the time of the Antarctic cold reversal from ~14.5 ka. Yet the warming that starts again at ~12 ka in the data lags in the 430 simulations with FWF.

Discussion and conclusion
Using the iLOVECLIM model of intermediate complexity, we have evaluated the impact of different ice sheet reconstructions, separating their different effects: the change of topography (and albedo), the change of bathymetry (and land-sea mask) and 435 the change of fresh water fluxes from the ice sheet melt.

Different ice sheet reconstructions
In our simulations, the ice sheets are prescribed based on two reconstructions: ICE-6G_C and GLAC-1D, following the PMIP4 protocol. We show that the two reconstructions result in differences in terms of amplitude and timing of climate changes. In 440 particular the global mean temperature is ~0.3°C warmer with GLAC-1D compared to the simulation with ICE-6_C, from the LGM to ~12 ka. Between 14 ka and 12 ka the evolution between the two simulations differ. Only the simulation with GLAC-1D displays a stalling in the temperature warming similar to data, albeit with a lag compared to data.

Bathymetry 445
In most past studies of deglacial simulation, the bathymetry was fixed, either at the LGM or the PI. We have shown here that the impact of accounting for bathymetric changes is a shift towards colder values, which is relatively limited for global mea n temperature, but more important regionally. This effect is mainly due to the different albedo with and without bathymetry changes. With evolving bathymetry, some grid cells are progressively changed from continent to ocean throughout the deglaciation, increasing the space available in the Nordic seas especially. This results in a larger sea ice area with a higher 450 albedo, which cools the climate and modifies the convection zones in the North Atlantic. As we demonstrate that bathymetry changes can impact the local and global climate, we advise to account for them (when technically possible) in deglacial simulations. However, changes of bathymetry can also trigger unrealistic changes, such as observed here in the simulations with the ICE-6G_C reconstruction at ~13.5ka around Antarctica. Using a higher resolution might help, but this shows that accounting for realistic bathymetry changes remains challenging. 455 The evolution of bathymetry depends on the ice sheet reconstructions. As the reconstructions are improved over time, the bathymetry change will also become more constrained. However, the model grid change resulting from the bathymetry change is not always straightforward when seaways become small and especially smaller than the grid cells: should a passage stay open in the model even if it results in a larger seaway, or should it be closed, even though it is not completely closed in the reconstructions? In this study, we have chosen to keep seaways open when they were open in the reconstructions, even if it 460 was smaller than the grid cell. Yet, it would be interesting to test the impact of closing them. This issue is related to the model resolution: with higher resolution models it is possible to keep small seaways open, while in lower resolution models the gri d cells are sometimes too large. In addition, we have chosen to keep the frequency of bathymetry update from the original frequency of the ice sheet reconstructions. This results in two different frequencies for the two reconstructions (500 and 100 years). Testing the impact of different update frequencies for the same reconstructions show some limited impact: a less 465 frequent update leads to delays in the climate response as the changes take place later. Yet this effect is small compared to the change in climate from the two different reconstructions with the same update frequency.
In LGM simulations, a previous study showed that accounting for bathymetry and land-sea mask change was crucial for the carbon cycle (Lhardy et al., 2021b). A change of bathymetry, especially the one associated ocean volume change, modifies the ocean carbon storage capacity. As it is a very large carbon reservoir compared to the atmosphere, even changes of around 3%, 470 the order of magnitude between LGM and PI volume changes, can largely i mpact atmospheric CO2. This effect remains to be tested for the deglaciation.
In our study, the ice sheets are prescribed from reconstructions. In reality, ice sheets and climate interact. While the use of prescribed ice sheets is a limitation of this study, this topic has been addressed in Quiquet et al. (2021) with an interactive ice sheet model coupled to the iLOVECLIM model. These simulations show similar results in terms of temperature and ocean 475 circulation changes. However, the bathymetry was fixed in such simulations. As shown here, the regional temperature can be impacted by the changes of bathymetry, which influences oceanic circulation and hence the transport of heat. In the simulations with evolving bathymetry, the temperature is colder over Greenland and Antarctica compared to the simulation with fixed bathymetry (Figure 7). In simulations with interactive ice sheets, this would impact the mass balance and could lead to a different evolution of ice sheets. 480

Fresh water fluxes
In iLOVECLIM, the AMOC is very sensitive to the fresh water flux, resulting in an AMOC shut dow n when freshwater input from ice sheet melting is included. The AMOC then remains in this collapsed state for several thousand years. The sensitivity to fresh water flux is very model dependant (Gottshalk et al., 2019), and some other models display less s ensitivity. For 485 example, with the MPI-ESM model (Kapsch et al., 2022), the same fresh water fluxes computed from the ice sheet melt from the ICE-6G_C and GLAC-1D reconstructions result in weaker AMOC but no collapse, more similar to our FWF/4 simulation.
In a simulation with HadCM3, a larger sensitivity of AMOC to the fresh water flux is obtained, but only for the 19 ka to 16 ka period as the simulation was not continued after (Ivanovic et al., 2018).
In most past studies of deglacial simulations (apart from Kapsch et al., 2022 andIvanovic et al., 2018), fresh water fluxes were 490 prescribed and not computed from the ice sheet melting. In addition, in most previous studies the areas where fresh water flux is added are prescribed. In this study the fresh water flux is computed from the ice sheet volume change and routed towards the ocean following the topography. In figure 13, we compare the fresh water fluxes used in those simulations. In most model studies, the fresh water flux is much smaller and / or with a different timing compared to the ones computed from the ICE-6G_C and GLAC-1D reconstructions. This leads to a smaller evolution of sea level equivalent than in the reconstructions 495 ( Figure 16). These models obtained a similar evolution of AMOC compared to proxy data and accordingly a relatively good evolution of temperature in Greenland, but the prescribed fresh water flux is in disagreement with the currently available ic e sheet reconstructions. It is noteworthy that wWhile in most models, fresh water fluxes are necessary to trigger large ocean circulation changes, one model has displayed changes in ocean circulation without fresh water fluxes (Zhu et al., 2014). In t his simulation, the ocean circulation change was due to the orographic change of the ice sheet only (Zhu et al., 2014). However, 500 the prescribed fresh water fluxes are also smaller than in the reconstruction leading to smaller sea level change compared to data. Larger fresh water fluxes similar to the data may lead to different and possibly degraded results compared to data.
It is thus difficult for models to explain the large change of ocean circulation and the related change of temperature in the Northern and Southern Hemisphere between 18 and 15ka by freshwater fluxes. Either the ice sheet and sea level reconstructions should be revisited, the modelled AMOC sensitivity modified, or other processes explaining the AMOC change tested. 505 Figure 16. Equivalent sea level evolution (m) due to the fresh water input (displayed on Figure 13).

Data availability 510
The source data of the figures presented in the main text of the paper will be available on the Zenodo repository.

Author contribution
All authors have developed the model code. NB has performed the simulations. NB ran the analysis and prepared the manuscript with contributions from all co-authors.