Climate and ice sheet evolutions from the last glacial maximum to the pre-industrial period with an ice sheet – climate coupled model

The last deglaciation offers an unique opportunity to understand the climate – ice sheet interactions in a global warming context. In this paper, to tackle this question, we use an Earth system model of intermediate complexity coupled to an ice sheet model covering the Northern Hemisphere to simulate the last deglaciation and the Holocene (26-0 ka BP). We use a synchronous coupling every year between the ice sheet and the rest of the climate system and we ensure a closed water cycle 5 considering the release of freshwater flux to the ocean due to ice sheet melting. Our reference experiment displays a gradual warming in response to the forcings, with no abrupt changes. In this case, while the amplitude of the freshwater flux to the ocean induced by ice sheet retreat is realistic, it is sufficient to shut down the Atlantic meridional overturning from which the model does not recover within the time period simulated. However, with reduced freshwater flux we are nonetheless able to obtain different oceanic circulation evolutions, including some abrupt transitions between shut-down and active circulation 10 states in the course of the deglaciation. The fast oceanic circulation recoveries lead to abrupt warming phases in Greenland. Our simulated ice sheet geometry evolution is in overall good agreement with available global reconstructions, even though the abrupt sea level rise at 14.6 kaBP is underestimated, possibly because the climate model underestimates the millenialscale temperature variability. In the course of the deglaciation, large-scale grounding line instabilities are simulated both for the Eurasian and North American ice sheets. The first instability occurs in the Barents-Kara seas for the Eurasian ice sheet 15 at 14.5 kaBP. A second grounding line instability occurs circa 12 kaBP in the proglacial lake that formed at the southern margin of the North American ice sheet. With additional asynchronously coupled experiments, we assess the sensitivity of our results to different ice sheet model choices related to surface and sub-shelf mass balance, ice deformation and grounding line representation. While the ice sheet evolutions differ within this ensemble, the global climate trajectory is only weakly affected by these choices. In our experiments, only the abrupt shifts in the oceanic circulation due to freshwater fluxes are able 20 to produce some millenial-scale variability since no self-generating abrupt transitions are simulated without these fluxes.

glacial terminations precede interglacial periods that show reduced ice sheets. The study of glacial terminations can help us 25 to understand the mechanisms behind large scale ice sheet retreat but also the key role of ice sheets within the global climate system.
During the last deglaciation (~21-7 kaBP), the sea level rose by about 120 metres to reach approximately its present-day level (Waelbroeck et al., 2002;Lambeck et al., 2014). This rise is mostly explained by the disintegration of the North American and 30 Eurasian ice sheets while Greenland and Antarctica together probably contributed to less than 20 metres (Whitehouse et al., 2012;Briggs et al., 2014;Lecavalier et al., 2014). The extent of the Northern Hemisphere ice sheets across the deglaciation is relatively well known although it can sometimes present large (>1 kyrs) dating uncertainties (Hughes et al., 2016;Dalton et al., 2020). However, the volume evolution of the individual ice sheets remains weakly constrained. In particular, sea level archives have suggested the presence of abrupt sea level rises standing out from the gradual sea level rise of the deglaciation 35 (Deschamps et al., 2012;Abdul et al., 2016;Harrison et al., 2019). These so-called melt water pulses suggest large-scale ice sheet instabilities but the contribution of the different ice sheets to these events remains debated (e.g. Liu et al., 2016).
Parallel to the non-linear ice sheet retreat, the atmosphere and the ocean have also undergone some large and abrupt variations.
For example, while atmospheric temperatures above Greenland rise gradually since the last glacial maximum (LGM), at the 40 onset of the Bølling-Allerød period at 14.7 kaBP, they rise abruptly by more than 10 • C in a few decades (Severinghaus and Brook, 1999;Buizert et al., 2014). After 500 years of interglacial conditions, the climate abruptly returns to a cold state during the Younger Dryas (Alley, 2000a) from which the temperatures rise again steadily to reach their Holocene values. The evolution of the oceanic conditions are more uncertain. It seems nonetheless that the North Atlantic Deep Water (NADW) was shallower at the LGM compared to today (Curry and Oppo, 2005). The 3D evolution of the water masses across the deglacia-et al., 1997;Opsteegh et al., 1998); a free surface oceanic general circulation model on a 3 • x3 • spherical grid which includes a thermodynamic sea ice model (CLIO, Goosse and Fichefet, 1999); and a dynamic vegetation and carbon allocation model (VECODE, Brovkin et al., 1997). iLOVECLIM has been extensively used to study millenial climate change during the Qua-95 ternary. For example, it has proven able to reproduce the glacial-interglacial variability of the hydrological cycle in the Tropics (Caley et al., 2014). It has also been used to study Heinrich events during the last glacial period (Roche et al., 2014b) or to investigate the processes responsible for changes in the carbon cycle during the last eight interglacial periods (Bouttes et al., 2018). With a similar model configuration to the one used in this work, iLOVECLIM results were included in the fourth phase Since Roche et al. (2014a), the model also includes a 3D thermomechanically coupled ice sheet model (GRISLI Ritz et al., 2001;Quiquet et al., 2018a). GRISLI solves the ice sheet mass conservation equation on a Cartesian grid. Like most ice sheet model, deformation is computed with a Glen flow law in which anisotropy is artificially accounted for using an flow 105 enhancement factor (E f ) that facilitates deformation induced by vertical shear. For the entire domain, the velocity field is the sum of velocity driven by vertical shearing (Shallow Ice Approximation, SIA) and the velocity driven by horizontal shearing (Shallow Shelf Approximation, SSA). In doing so, the SSA is used as a sliding law (Bueler and Brown, 2009;Winkelmann et al., 2011). Basal dragging τ b is assumed to follow a linear friction law: (1) 110 where β is the basal drag coefficient and u b is the velocity at the base. Cold based grid points have a virtually infinite friction at the base (5×10 5 Pa yr m -1 ), while floating ice shelves have no friction. For grid points at the pressure melting point, we use a friction computed from the effective pressure at the base of the ice sheet N : with c f a parameter that has to be calibrated. For the experiments shown here, we impose an ice flux at the grounding line 115 that follows the analytical solution of Tsai et al. (2015). Calving at the ice shelf edge occurs if the ice thickness falls below a critical threshold and if the upstream Lagrangian ice flux does not allow to maintain an ice thickness above this threshold.

Ice sheet model coupling
The inclusion of GRISLI into iLOVECLIM has been presented in Roche et al. (2014a). However, the coupling procedure has 125 been largely modified from this work. In particular, we have substantially improved the computation of surface and sub-shelf mass balance. Water conservation between GRISLI and the rest of the climate model has also been considerably improved.
Details on this coupling is given in the following while its schematic representation is shown in Fig. 1. It is important to mention that only the Northern Hemisphere ice sheets are interactively simulated, while the Antarctic ice sheet topography and ice mask remains prescribed.

Surface mass balance
In Roche et al. (2014a), the ice sheet surface mass balance (SMB) was computed from the annual mean precipitation and the annual and July mean near-surface air temperature using a positive degree day method (Reeh, 1989). Although computationally inexpensive and easy to implement in a model, this method does not account for some important physical quantities that 135 influence the SMB. In particular, the surface shortwave radiation is only implicitly taken into account, through the temperature.
Instead, we use here the insolation temperature melt method (ITM) following Pollard (1980) and van den Berg et al. (2008).
The amount of melt M s over one time step dt is in this case: With T s is the near surface air temperature, SW s is the shortwave radiation at the surface, α is the surface albedo, ρ w is the 140 density of liquid water and L m is the specific latent heat of fusion. λ and c rad are empirical parameters that need calibration.
In the literature, this calibration has been performed on observations of present-day glaciers. The λ parameter is generally set to 10 W m -2 K -1 (Pollard, 1980;van den Berg et al., 2008;Robinson et al., 2010). The parameter c rad is less constrained and is adjusted for the region considered (van den Berg et al., 2008). It is set to -50 W m -2 in Pollard (1980), it ranges from -40 to -60 W m -2 in Robinson et al. (2010), whilst it is equal to -117 W m -2 in van den Berg et al. (2008). We used λ =10 W m -2 K -1 145 and c rad =-40 W m -2 . However, iLOVECLIM presents an important warm bias in Eastern North America and a cold bias in Northern Europe that lead to unrealistic simulated ice sheet under glacial forcing. A problem also identified in Heinemann et al. (2014). To account for this, we use a local modification of the melt parameter c rad to partially correct these temperature biases. To this aim, we compute the annual mean temperature bias with respect to ERA-interim (Dee et al., 2011) and use a linear correction in which a +10 • C bias leads to c rad =-80 W m -2 (instead of the reference value of -40 W m -2 ).
at the resolution of the ice sheet model with the ITM method (Eq. 3). The near-surface air temperature and the SMB are accumulated along the course of the year to generate the yearly forcing fields required by the ice sheet model. We made a few adjustments compared to the downscaling procedure presented in Quiquet et al. (2018b). In particular, some large-scale climate fields are now bi-linearly interpolated onto the high-resolution grid before the energy and moisture computation. This 160 prevents the strong discontinuities that could exist between two sub-grid points belonging to two different large-scale grid cells.

Sub-shelf melt rate
The sub-shelf melt rate in Roche et al. (2014a) was imposed arbitrarily to an homogeneous and constant value for the entire Northern Hemisphere. Instead, we use here a physically based computation of the sub-shelf melt rate following Beckmann and 165 Goosse (2003). For each vertical oceanic layer, z, we estimate the potential sub-shelf melt rate as: wit c p specific heat capacity of sea water, ρ i is the density of ice, γ T is the thermal exchange velocity and T F (z) is the thermal forcing at depth z, defined as the difference between the ambient temperature and the temperature of the salinity dependent freezing point. F g is a weakly constrained dimensionless parameter and can be changed to explore the response of the ice 170 sheet to different sub-shelf melt sensitivities to oceanic temperature change. In our reference experiment, we chose a parameter value (15×10 -3 ) that produce about 0.1 m yr -1 in the Arctic, a value similar to what is experiencing the Ross ice shelf today in Antarctica. In addition, in order to avoid unrealistic ice shelf expansion over the deep ocean we also impose a high sub-shelf melt rate of 20 m yr -1 where the bathymetry is greater than 1500 m. Also, to mimic the fact that observed melt rates are greater in the vicinity of the grounding line, we double the value of the inferred melt rate in Eq. 4 for the floating points that are in 175 contact with the grounding line. Eq. 4 is computed for each oceanic timestep (one day) and integrated over the year in order to provide the yearly forcing needed by the ice sheet model. There is no downscaling of the sub-shelft melt rate to the highresolution ice sheet model grid, except that the depth of the ice shelf draft is used to determine the vertical layer z in Eq. 4 that produces the melt.
On the other hand, freshwater fluxes resulting from the ice sheet melting are transferred to the oceanic model. In Roche et al. (2014a), the total ice sheet volume variation was transferred to the continental routing scheme assuming a uniform distribution 190 over the ice sheet. Only the calving flux was separated from the total volume variation to eventually feed an iceberg model (Bügelmayer et al., 2015). This method has the advantage to ensure a closed water budget within the model but the spatial information about ice sheet runoff is lost.  The climate model uses time-varying information of greenhouse gases (Lüthi et al., 2008) and insolation (Berger, 1978). The carbon dioxide mixing ratio evolution and the 65 • N insolation in June is depicted in Fig. 3. We use a recent implementation 205 of the last glacial maximum bathymetry at 21 kaBP (Lhardy et al., 2020), which is left unchanged for the duration of the experiments. Topography and ice mask are both provided by the ice sheet model.
To define our initial state, we run uncoupled ice sheet and climate experiments. First, we run the climate model using the last glacial maximum boundary conditions for 3000 years. In this case, the ice sheet topography and ice mask correspond to the one 210 of the GLAC-1D reconstruction (Tarasov et al., 2012;Tarasov and Peltier, 2002;Briggs et al., 2014)   Our reference experiment (DGL) is an ice sheet -climate experiment, synchronously coupled. This experiment starts at 26 kaBP and uses the initial conditions presented in Sect. 2.3.1. The climate and the ice sheets used as initial conditions are not fully consistent between each other since they have been obtained with uncoupled long-term equilibriums. As such, the first thousand years or so of our experiments have to be discussed with care since part of the response can arise from artefacts due to the start of the coupling.

240
In addition to this reference experiment, we have performed additional synchronously coupled experiments for the first set of experiments which aims at investigating the importance of oceanic changes in shaping the last deglaciation.
First, we have run experiments in which the amount of freshwater is reduced in order to gradually limit their influence.
In DGL_FWF/2 and DGL_FWF/3 we divide the flux resulting from ice sheet melting by two, respectively three, while in Here, we discard completely the role of freshwater flux to the ocean in the accelerated experiments. The ADGL experiments is the accelerated counterpart of the DGL experiments, and as such will define the new reference for the accelerated experiments.

270
The other accelerated experiments are used to assess the sensitivity of our simulated deglaciation to important processes related to ice sheet dynamics: modelling choices for ice dynamics, for the surface mass balance and for the sub-shelf melt rate.
First, we explore two aspects related to ice dynamics: grounding line dynamics and ice deformation. The ice sheet model GRISLI accounts for two formulations of the flux at the grounding line. For the Antarctic ice sheet, the use of Schoof (2007) 275 instead of Tsai et al. (2015) leads to slower grounding line retreat during deglaciation phases (Quiquet et al., 2018a). For this reason, in the ADGL_schoof experiment we use the Schoof (2007) formulation of the flux at the grounding instead of Tsai et al. (2015). A second aspect for ice dynamics is the choice of the flow enhancement factor E f which is a tuned parameter that has consequences on the ice velocity. In the ADGL_ef experiment we use a larger flow enhancement factor (larger velocities) since the simulated North American ice thickness at the LGM is overestimated (Fig. 2).

280
Then, to explore the sensitivity of our results to the surface mass balance we have performed two experiments in which the weakly constrained melt parameter c rad (Eq. 3) is changed. In ADGL_accplus we use a smaller value for this parameter in order to reduce surface melt to delay the deglaciation. In the ADGL_nocor experiment we use an homogeneous value of c rad instead of using the spatial heterogeneous value defined from the temperature bias.
Finally, to assess the sensitivity of our results to the sub-shelf melt rate, in the ADGL_bmbplus we enhance the sub-shelf melt 285 rate to increase the relative importance of oceanic changes with respect to atmospheric changes.
The list of the different experiments is available in Tab. 1.

Results
In this section we first describe the general evolution of the simulated climate in the synchronously coupled experiments before 290 examining the ice sheet changes. Then we examine the results for the accelerated asynchronously coupled experiments to infer the sensitivity of our results to different ice sheet evolutions.

Climate evolution in the synchronously coupled ice sheet -climate experiments
The simulated global mean surface temperature evolution for the synchronously coupled ice sheet climate experiments is shown 295 in Fig. 3, together with the strength of the AMOC. In response to the forcings, the different experiments produce a gradual warming from the last glacial maximum towards its maximum value during the Holocene. The glacial-interglacial temperature difference ranges from 3.1 to 3.8 • C and is in good agreement with palaeo-temperature stack (Shakun et al., 2012), even though iLOVECLIM is one of the warmest model at the LGM within the PMIP4 ensemble (Kageyama et al., 2020) . The glacial-interglacial temperature difference is mostly explained by the cold temperatures at the LGM resulting from the large 300 ice sheets that induce higher surface elevations and a strong albedo effect. A polar amplification is simulated since the northern and southern high latitudes both show a greater temperature difference from the pre-industrial compared to the Tropics (Fig. 4). This pattern is consistent with recent reconstructions (e.g. Tierney et al., 2020, shown in Fig. 4d), even though with a smaller amplitude in our model. However, our simulated glacial-interglacial temperature difference is within the range of other estimates (4 ± 0.8 • C, Annan and Hargreaves, 2013 erence DGL simulation shows a slight decrease in temperature for about 2 kyrs followed by a moderate warming until~7 kaBP. On the contrary, the experiment in which the freshwater flux are discarded (DGL_noFWF) displays a brief period during which the temperature ceases to increase followed by a sharp temperature increase. In this case, the maximal surface temperature is reached at 10 kaBP after which there is a slight decrease until 7 kaBP. DGL_FWF/3 shows a very similar temperature change as DGL_noFWF while DGL_FWF/2 presents similarities with both DGL and DGL_noFWF. This temperature evolution is in While some experiments show very abrupt shifts in the ocean, the atmospheric temperature evolution is nonetheless mostly 345 gradual. This is visible at the global scale (Fig. 3b), but also when examining the temperature change above the Greenland ice sheet (Fig. 6a).

Simulated ice sheet changes
The large-scale differences amongst the different experiments discussed in Sect. 3.1 are largely driven by differences in the 360 amount of the freshwater released to the ocean related to ice sheet melting. This freshwater flux is shown in Fig. 7a for the reference experiment DGL. Even though this flux displays some variability, its evolution is generally gradual and shows a maximum around 14 kaBP where it peaks above 0.3 Sv (1 Sv corresponds to 10 6 m 3 /s) with 100-yr mean values about 0.23 Sv.
In Fig. 7a we also show the meltwater flux computed from the ice thickness changes in the ICE-6G_C and GLAC-1D geologically constrained reconstructions. These fluxes have the same order of magnitude of the simulated flux in the DGL experiment. Instead, the model produces important fluxes (greater than 0.1 Sv) over a few thousand years. An other way to discuss these fluxes is to integrate them in time to have an idea of the total ice sheet volume evolution through the deglaciation (Fig. 7b). In doing so, we can see that the coupled iLOVECLIM-GRISLI model setup produces an ice volume evolution in general agree- This is mostly due to our methodology used to define the initial state for the coupled experiments. When the coupling starts, at the beginning of our experiments, there is an abrupt change in the climate model in terms of ice mask and surface elevation, from GLAC-1D to our spun-up ice sheets. Our spun-up ice sheets at 26 kaBP (Fig. 2a) show a higher North American ice sheet surface elevation than the GLAC-1D reconstruction used during the climatic spin-up, suggesting an overestimation of 390 the precipitation in this area. When the coupling starts, this precipitation bias is amplified due to higher surface elevation and related increased orographic precipitation. The iLOVECLIM climate model likely shows an underestimation of the elevation desertification effect over the ice sheets (Quiquet et al., 2018b). The simulated volume of the Eurasian ice sheet displays a similar evolution than the North American ice sheet with a maximum around the last glacial maximum. This agrees well with the GLAC-1D reconstruction. Given its smaller volume, the absolute rate of volume loss is smaller for the Eurasian ice sheet

395
(1.3 m per millenia) compared to the one of the North American ice sheet (5.7 m per millenia). However, the Eurasian ice sheet already lost half its volume by 14.5 kaBP where this occurs at 12.8 kaBP for the North American ice sheet. The Greenland ice sheet presents only a small volume reduction of 2.6 m of sea level equivalent, in good agreement with the reconstructions.
However, the Greenland ice sheet volume at the end of the simulation is largely overestimated compared to the present-day observations (about 40% volume overestimation). As for the total volume, the individual ice sheets deglaciate faster in the 400 DGL_noFWF experiment. This is particularly visible for the North American ice sheet for which there is a difference of one thousand years at about 11 kaBP.
A map of the simulated ice sheet configuration for selected snapshots is shown in Fig. 9. Greenland ice sheet expands considerably onto the continental shelf during the glacial period and retreats until about 10 kaBP.

425
It does not display any substantial change in the ice extent during the Holocene but it displays some ice elevation changes.
The ice elevation evolution near the summit shows a maximum at about 10 kaBP and decreases afterwards in agreement with palaeo-elevation reconstructions at the deep ice core drilling sites (Vinther et al., 2009).
In Fig. 10 we present the rate of total ice mass change and its individual components: surface mass balance, basal mass balance 430 and calving. The total mass change remains positive until 20.5 kaBP due to a positive integrated surface mass balance, not entirely compensated by the basal mass loss (mostly sub-shelf melt) and calving. After this date, the total mass change becomes negative for the rest of the duration of the experiment. The total mass loss peaks at -9.7 × 10 3 Gt yr -1 at 13.8 kaBP when surface ablation and loss by calving almost synchronously display a maximum (surface ablation slightly precedes the calving increase). At this date, the mass loss due to the ocean and lake represent more than half the loss by surface mass balance. In 435 fact, if both basal mass loss and calving remain almost constant until 14.5 kaBP (-1.4 × 10 3 Gt yr -1 ), they nonetheless show some variability after this date. These fluxes are maximal at the time of the grounding instabilities shown in Fig. 9 for both the Eurasian (14.5-13.5 kaBP) and the North American (12.8-10 kaBP) ice sheets. While the mass loss is primarily driven by surface ablation until 12.8 kaBP, after this date the oceanic and lake forcing become the major driver for the ice sheet retreat.
The total mass loss finally reaches zero (ice sheet equilibrium) at 6.5 kaBP.

Accelerated experiments to assess specific sensitivities
The aim of this section is to assess the sensitivity of the simulated climate evolution to the choice of critical ice sheet model parameters and assumptions. To do so, we have performed additional experiments in which the forcings are accelerated. Three mass balance and sub-shelf melting rates.
The evolution of some large-scale climate variables for these additional experiments are shown in Fig. 11. Since we do not feed back the freshwater related to ice sheet melting to the ocean in the accelerated experiments, they have to be compared to the DGL_noFWF experiment. The reference accelerated experiment ADGL (black in Fig. 11) is in fact similar to the DGL_noFWF 450 (grey): rapid temperature increase and active Atlantic circulation throughout the deglaciation. However, the accelerated experiment displays larger ice sheets than the non-accelerated (about 8 m of sea level equivalent at 14 kaBP) and as a result a colder climate (~0.4°C in global mean surface temperature at 14 kaBP). If the timing of the ice sheet retreat can be slightly different, the overall pattern of this retreat is only weakly affected by the acceleration factor.

455
The two experiments related to ice sheet dynamics (ADGL_ef and ADGL_schoof) do present some differences in their simulated ice sheet volume. The increased enhancement factor (ADGL_ef) leads to thinner ice sheets (smaller ice volume) and, as such, deglaciates faster than the reference accelerated experiment (ADGL). The experiment in which we use the formulation of Schoof (2007) instead of Tsai et al. (2015) (ADGL_schoof) also produces a lower ice sheet volume during the glacial period.
However, this experiment shows a slower ice sheet retreat during the deglaciation compared to the reference ADGL experi-460 ment. This is mostly related to the greater grounding line sensitivity in the formulation of Tsai et al. (2015), already shown in (Quiquet et al., 2018a) for the Antarctic ice sheet. These differences in terms of ice sheet evolution have nonetheless only a limited impact on the climate evolution. The ADGL_ef produces a slightly more rapid warming during the deglaciation (related to the smaller ice sheets) while it is the contrary for the ADGL_schoof (slower ice sheet retreat). The Atlantic circulation is also weakly impacted by the different ice sheet evolution. Only the ADGL_ef produces a decrease in the overturning circulation 465 strength slightly earlier than the ADGL experiment while the ADGL_schoof displays insignificant differences.
The two experiments related to modification of the surface mass balance parameters induce larger simulated ice sheet volume differences. Associated with a larger ice sheet surface mass balance, the ADGL_accplus produces larger ice sheet volumes throughout the whole simulated time period. For this experiment the maximal ice volume is reached circa 19 kaBP and it is 470 larger by about 15 m of sea level equivalent than the ADGL experiment. This excess ice also explains the delayed ice sheet retreat: at 10 kaBP the simulated ice sheets still represent about 45 m drop in eustatic sea level in the ADGL_accplus experiment compared to about 8 m in the ADGL experiment. This has consequences on the simulated climate: i-the global mean temperature rises more slowly, reaching eventually a comparable value to the reference simulation at 0 kaBP; ii-the phase of very active overturning in the middle of the deglaciation is extended by 2 ka. Even though the ADGL_accplus experiment dis-475 plays larger ice sheet volume, the pattern of the ice sheet retreat is similar to the one of the ADGL experiment. On the contrary, the ADGL_nocor experiment provides alternative ice sheet histories. In the ADGL experiment the Barents-Kara sector of the Eurasian ice sheet is almost fully deglaciated at 13.5 kaBP, while it is the case only after 8.5 kaBP in the ADGL_nocor experiment. Conversely, the North American ice sheet retreats faster in the ADGL_nocor experiment. This is a direct consequence of the cold temperature bias in Northern Europe and the warm bias in North America. If the climate evolution is not drastically changed as a result of these different ice sheet chronologies, it nonetheless shows some interesting differences. The overturning circulation remains moderate for a longer time period compared to the ADGL experiment since it increases only after 15 kaBP (w.r.t. 17.5 kaBP in ADGL). As a result the global mean temperature in ADGL_nocor is colder than in ADGL even though it shows smaller ice sheets at least until 15 kaBP. The oceanic circulation in the model seems largely affected by the Eurasian ice sheet size.

485
Finally, the experiment in which we increase the sub-shelf melting rate, ADGL_bmbplus, shows only negligible changes with respect to the ADGL experiment. This suggests that, in our model, the ice sheet retreat is mostly driven by surface ablation and not sub-shelf melt.

Discussion
We have shown that in our reference experiment the freshwater flux to the ocean resulting from ice sheet melting leads to a progressive weakening of the Atlantic overturning circulation from the last glacial maximum, eventually leading to a com- considering the depth of the freshwater release, its seasonality or the impact of the iceberg transport.

510
The simulated temperature change during the deglaciation is generally very gradual with no abrupt transitions. For example, in our experiments, over the Greenland ice sheet, the local temperature change is strongly correlated to the global mean tem-perature change and most of the time does not display abrupt events such as the one recorded in ice cores (Alley, 2000b). In fact only the abrupt AMOC recoveries in certain experiments (DGL_FWF/3 at 10.7 kaBP and DGL_brines at 3.8 kaBP) are able to produce abrupt temperature changes in Greenland comparable to the ice core record. Since these AMOC recoveries are 515 lacking in the majority of our experiments we generally largely underestimate the millenial scale variability observed at high latitudes. This variability could largely influence the ice sheet evolution. For example, since the Bølling-Allerød warming is not simulated in our model, we are not able to quantify its impact on the North American or the Eurasian ice sheets (Gregoire et al., 2016;Brendryen et al., 2020).

520
In addition, within the experiments presented here, only changes in the AMOC related to freshwater flux are able to produce some abrupt temperature changes. For example, all the accelerated experiments, in which this process is not considered, produce a smooth temperature increase since the LGM. However, these experiments show different ice sheet evolutions with some rapid ice sheet retreat at times. This suggest that in our model and for the time period simulated, external forcing and ice sheet changes alone are not able to produce millenial-scale climate variability without invoking freshwater hosing.

525
Finally, in our experiments we did not consider the potential changes of the Antarctic ice sheet since we use a constant topography and ice mask in the Southern Hemisphere. Similarly we do not take into account the freshwater flux resulting from Antarctic ice sheet retreat from the last glacial maximum. This simplification was motivated by the fact that an earlier study already identified that freshwater hosing around Antarctica with our model has a negligible impact on the simulated climate 530 . In this region, the circumpolar current tends to rapidly dilute the released freshwater leading to a very limited impact on vertical oceanic mixing. However, the gradual retreat of the ice sheet from the continental shelf margin can also facilitate the sinking of brines and as such enhance dense water formation. If the sinking of brines around Antarctica seems to play a moderate role in our experiments, it can nonetheless produce an abrupt AMOC recovery at 3.8 kaBP, not occurring in the reference experiment. As such, this process should be more thoroughly investigated with, e.g., interactive Antarctic topog-535 raphy and bathymetry.

Conclusions
In this paper, we have presented climate model experiments in which the Northern Hemisphere ice sheets are synchronously coupled to the rest of the system (atmosphere and ocean). For the majority of our experiments, the atmospheric changes are     for the freshwater release to the ocean due to ice sheet melting is shown in grey (DGL_noFWF). The other lines are accelerated simulations (factor of acceleration of 5) and they similarly do not account for the freshwater flux to the ocean. The accelerated reference experiment is in blue (ADGL). The dark green line is in experiment in which we use the Schoof (2007) formulation of the flux at the grounding line (ADGL_schoof) instead of Tsai et al. (2015). The light green line is an experiment for which we use a larger enchancement factor (3.5 instead of 1.8, ADGL_ef). The light orange line is a version of the model with a lower c rad coefficient (-50 instead of -40 W m -2 , ADGL_accplus) which induces a more positive surface mass balance. The dark orange is for an experiment in which with do not apply the spatial correction of the c rad parameter (ADGL_nocor). The blue line is for an experiment with increase sub-shelf melting rate (ADGL_bmbplus).