Ocean Response in Transient Simulations of the Last Deglaciation Dominated by Underlying Ice‐Sheet Reconstruction and Method of Meltwater Distribution

The last deglaciation was characterized by drastic climate changes, most prominently melting ice sheets. Melting ice sheets have a significant impact on the atmospheric and oceanic circulation, due to changes in the topography and meltwater release into the ocean. In a set of transient simulations of the last deglaciation with the Max Planck Institute for Meteorology Earth System Model we explore differences in the climate response that arise from different boundary conditions and implementations suggested within the Paleoclimate Modeling Intercomparison Project ‐ Phase 4 (PMIP4) deglaciation protocol. The underlying ice‐sheet reconstruction dominates the simulated deglacial millennial‐scale climate variability in terms of timing and occurrence of observed climate events. Sensitivity experiments indicate that the location and timing of meltwater release from the ice sheets into the ocean are crucial for the ocean response. The results will allow a better interpretation of inter‐model differences that arise from different implementations proposed within the PMIP4 protocol.

forcings, trigger mechanisms and dynamic feedbacks the protocol specifies options for different ice-sheet boundary forcing and methods for distributing meltwater from retreating ice sheets. While the latter allows the full range of models to participate in the PMIP4 deglaciation exercise, it is important to disentangle the differences in the model response due to the model setup.
Past studies have investigated the model sensitivity to different ice-sheet boundary conditions (e.g., Bakker et al., 2020;Gong et al., 2015;Hofer et al., 2012;Justino et al., 2005;Liakka & Nilsson, 2010;Liakka et al., 2016;Roe & Lindzen, 2001;Ullman et al., 2014). Using a variety of different model systems, ranging from simplified linear models of the atmosphere to complex fully coupled atmosphere-ocean models, these studies found a strong sensitivity of the modeled climate of the last glacial maximum (LGM; about 21 ka) to the ice-sheet topography. High LGM ice sheets significantly impact the atmospheric and oceanic circulation (see Kageyama & Valdes, 2000;Löfverström & Lora, 2017;Löfverström et al., 2014;Sherriff-Tadano et al., 2018;Zhu et al., 2014). A high Laurentide ice sheet leads to a zonalization and enhancement of the eddy-driven jet stream over the North Atlantic, associated with an increase in surface winds and enhanced wind-driven gyres in the North Atlantic, favoring a stronger AMOC. However, sensitivity experiments by Ullman et al. (2014) and Bakker et al. (2020) showed that the climate response is dependent on the ice-sheet reconstruction used as boundary condition; small deviations in the height of the North American ice sheet between reconstructions can lead to a significantly different atmospheric and oceanic circulation. This indicates that uncertainties in the climate response are largely influenced by the glacial boundary conditions. The inflow of meltwater from retreating ice sheets into the ocean can cause large changes in the AMOC (e.g., McManus et al., 2004), as it potentially reduces the overturning strength. Changes in the AMOC in turn are associated with significant climate changes (Broeker, 1991;Keigwin et al., 1991). Several modeling studies have investigated the response of the AMOC to freshwater forcing in form of so-called hosing experiments (e.g., Kageyama et al., 2010Kageyama et al., , 2013Lohmann et al., 2020;Maier-Reimer & Mikolajewicz, 1989;Schiller et al., 1997;Stouffer & Manabe, 2004;Stouffer et al., 2006;Smith & Gregory, 2009). Different methods have been applied to release the freshwater into the ocean, for example, freshwater distribution over latitude belts (Stouffer & Manabe, 2004;Stouffer et al., 2006), certain areas of the North Atlantic or Arctic Smith & Gregory, 2009), or point hosing (Condron & Winsor, 2012;Liu et al., 2018;Lohmann et al., 2020;Maier-Reimer & Mikolajewicz, 1989;Schiller et al., 1997). These studies showed that the AMOC is very sensitive to the location of the freshwater hosing. Specifically, it was shown that the AMOC responds more strongly to freshwater if it is released at higher latitudes of the North Atlantic or Arctic, close to or dynamically upstream of the deep water formation sites (Lohmann et al., 2020;Manabe & Stouffer, 1997;Rahmstorf, 1996;Roche et al., 2010;Smith & Gregory, 2009). However, Smith and Gregory (2009) found that if a certain magnitude of freshwater forcing is exceeded the location becomes only of secondary importance and the AMOC responds similarly, independent of where freshwater is added in the North Atlantic.
Using an ensemble of eight transient deglacial simulations with the Max Planck Institute for Meteorology Earth System Model (MPI-ESM), we extend previous studies that investigated the climate response to multiple icesheet boundary conditions in time-slice (e.g., Bakker et al., 2020;Liakka et al., 2016;Ullman et al., 2014) and hosing (e.g., Maier-Reimer & Mikolajewicz, 1989;McManus et al., 2004;Stouffer & Manabe, 2004) experiments. We aim to answer how sensitive the transient model climate of the last deglaciation is to the underlying ice-sheet boundary conditions, the implementation of meltwater distribution and the model version. By following the PMIP4 deglaciation protocol, we explore how PMIP4 models might respond to differences in the proposed experimental design. The intention is not to evaluate the simulations in regard to paleo-oceanographic evidence but to highlight differences that arise from the methods proposed in the protocol. This will help to interpret differences in the response of the models participating in the PMIP4 deglaciation exercise. We focus on the Northern Hemisphere, where ice-sheet changes were largest throughout the last deglaciation.

Model and Experimental Design
For the transient simulations of the last deglaciation ice sheets were prescribed from two different ice-sheet reconstructions. Meltwater release was calculated from changes in the ice sheet thickness in the reconstructions. To test the sensitivity to the meltwater distribution from melting ice sheets and the sensitivity of the transient climate response to different model parameters, experiments with different methods of meltwater release and model versions were performed, respectively. We apply three model versions, hereafter referred to as P1, P2 and P3. As P2 and P3 include not only a different tuning of parameters than P1 but also different parameterizations and bug fixes, we focus our detailed analysis on the simulations with model versions P2 and P3. Details about the model versions are described in Text S1 of Supporting Information S1. The modeled climate differs between versions. P2 and P3 are significantly colder during the glacial than P1, indicating a stronger sensitivity to greenhouse gas changes ( Figure 1c). The glacial cooling over the North Atlantic, which is associated with a reduced meridional heat transport due to a weaker AMOC, is more realistic in P2 and P3 than in P1 (Figures 1b and 1d; Tierney et al., 2020). The colder model versions P2 and P3 are also more sensitive to small changes in the meltwater forcing than P1, as they are closer to the temperature threshold at which the AMOC becomes more sensitive to perturbations (see discussion in Section 3.1; Klockmann et al., 2018;Oka et al., 2012). All model experiments are summarized in Table 1.

Experimental Design
For each experiment, the model was integrated from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). Ice sheets and surface topographies were prescribed from either GLAC-1D (Briggs et al., 2014;Tarasov et al., 2012) or ICE-6G (Peltier et al., 2015) reconstructions. All forcing fields are updated every 10 years of the simulations and initiate changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask (Meccia & Mikolajewicz, 2018;Riddick et al., 2018). For this, all forcing fields are interpolated to a 10-year resolution. For each ice-sheet reconstruction, meltwater from ice sheets is calculated as the temporal derivative of ice thickness at grid points covered by grounded ice sheets (see Figure 1a). The derived meltwater is then distributed by the hydrological discharge model and finally released into the ocean as freshwater (for details, see Meccia & Mikolajewicz, 2018). Thereby salt is conserved and changes in water content within the ESM are consistent with the prescribed ice-sheet changes. In sensitivity experiments, we follow other suggestions within the PMIP4 deglaciation protocol and either distribute meltwater equally over all grid cells (land and ocean) or remove meltwater from the system. For the latter experiment, salt and water are not conserved throughout the simulation.

Results and Discussion
Using MPI-ESM-CR, we performed a comprehensive set of sensitivity experiments following the PMIP4 last deglaciation protocol (Ivanovic et al., 2016). We find that the ice-sheet boundary conditions, the method of meltwater distribution and the background climate (realized through different model tuning and parameterizations) play a crucial role for the deglacial climate. While MPI-ESM-CR did not participate in CMIP6/PMIP4 activities, the simulated climate of all experiments presented here fits well within the range of the PMIP4 models for pre-industrial and LGM climate conditions (see Text S2 and Figures S1 and S2 in Supporting Information S1).

Ice-Sheet Boundary Conditions Crucial for Glacial Climate
Meltwater forcing derived from GLAC-1D and ICE-6G differs significantly in terms of its millennial-scale variability and the timing of major melt events ( Figure 1a). These differences partly arise from a different native temporal resolution, which is 100 and 500 years for GLAC-1D and ICE-6G, respectively. For example, the major meltwater pulses (MWP1a and MWP1b at about 14.5 ka and 11.5 ka, respectively) are more confined in GLAC-1D than ICE-6G and the peak discharge occurs at different times. Further, the GLAC-1D meltwater release shows larger variability throughout the last deglaciation. Accordingly, the simulations with GLAC-1D boundary conditions show an enhanced millennial-scale variability over the North Atlantic and globally (Figures 1b-1d). In the following, we focus in more detail on the differences between the simulations with ICE-6G and GLAC-1D reconstructions.

Simulation name
Reconstruction Meltwater forcing location Note. GLAC-1D and ICE-6G indicate the underlying ice-sheet reconstructions, P1-P3 different parameter tuning within MPI-ESM-CR, and Glob and noMW the meltwater implementation. See Section 2 for details on the experimental design.

LGM
Although the meltwater history of the two reconstructions before about 15 ka is very similar, except for its millennial-scale variability, the ice-sheet height differs significantly (Figures 1a and 2b; see also Ivanovic et al., 2016). During the LGM, the Laurentide ice sheet is significantly higher in ICE-6G than GLAC-1D and the Fennoscandian/British ice sheet is lower, regionally by more than 950 and 650 m, respectively (Figure 2b). The relatively high Laurentide ice sheet in ICE-6G leads to changes in the zonal 250 hPa wind speed over the North Atlantic, associated with a zonalization and enhancement of the eddy-driven jet stream in the simulations with ICE-6G boundary conditions (Ice6G_P2 and Ice6G_P3; Figure S3 in Supporting Information S1). Similar changes in the 250 hPa winds also occur in Ice6G_P1 but the jet shift is less distinct ( Figure S4 in Supporting Information S1). Merz et al. (2015) found that non-topographic forcings, such as orbital forcings, greenhouse gases and sea-surface temperatures can play a crucial role for the jet position. Further, Andres and Tarasov (2019) showed that the jet position is dependent on changes in the background climate through changes in the position of the polar front and sea-ice margin. Our results confirm such sensitivity: Model version P1 shows only a weak cooling over the North Atlantic during the LGM and a less distinct shift of the jet position than P2 and P3 ( Figure S4 in Supporting Information S1). The weak cooling of North Atlantic temperatures during the LGM in P1 allows it to maintain a sea-ice margin that is similar to present day (not shown) and thereby a strong and stable AMOC, with a significant North Atlantic deep water cell and deep water formation that is extended toward the eastern North Atlantic ( Figure S5 in Supporting Information S1). Hence, the AMOC response to the ice sheet boundary conditions in model version P1 is significantly weaker than in P2 and P3. These findings are in line with Oka et al. (2012) and Klockmann et al. (2018), who showed that a temperature threshold exists at which the AMOC becomes more sensitive to perturbations. The colder model versions P2 and P3 are closer to such temperature threshold and the AMOC is more sensitive to small changes in the forcing than the warmer model. These findings indicate that the model climate affects the atmospheric and oceanic response in the simulations.
The enhancement of the jet in the simulations with ICE-6G ice sheets leads to stronger surface winds and an increase in the wind-driven gyres ( Figure S6 in Supporting Information S1). This is associated with an increased northward transport of salty waters and results in a stronger AMOC in the ICE-6G simulations. Hence, a higher Laurentide LGM ice sheet in ICE-6G results in a stronger AMOC (Figure 1), which is in line with Ullman et al. (2014) and other studies (e.g., Kageyama & Valdes, 2000;Löfverström & Lora, 2017;Sherriff-Tadano et al., 2018). The simulations with the GLAC-1D topography do not show a zonalization of the jet stream during the LGM and an enhancement of the wind driven ocean transport. The AMOC in Glac1D_P2 and Glac1D_P3 is more than 4 Sv weaker and North Atlantic surface temperatures are more than 2.6 K lower than in Ice6G_P2 and Ice6G_P3 during the LGM.
Interestingly, we find the largest shifts in the jet over the eastern North Atlantic. Andres and Tarasov (2019) have shown that the jet over the eastern North Atlantic is not very sensitive to the marginal position of the Laurentide ice sheet. This indicates that the shifts of the eastern jet in this study are likely triggered by changes in the stationary waves (e.g., Löfverström et al., 2014), transient eddies (e.g., Merz et al., 2015), and/or wave breaking (Li & Battisti, 2008).

Deglaciation
Due to the reduction in the differences between the ice sheets' height in ICE-6G and GLAC-1D, the differences in temperatures and the AMOC between the simulations with GLAC-1D and ICE-6G boundary conditions get smaller throughout the deglaciation (Figure 1; see also Ullman et al., 2014).
In the reconstructions, an increase in meltwater release mark the onset of the deglaciation (see Figure 1a). However, as the amount of meltwater injected into the ocean is small until the occurrence of MWP1a in both reconstructions, the AMOC remains relatively strong in all simulations until about 14.5 ka (Figure 1b). Interestingly, none of the simulations shows a BA warm period, characterized by a strong AMOC and warm North Atlantic temperatures (Figures 1b and 1d; e.g., Clark et al., 2002;Weaver et al., 2003). Weaver et al. (2003) argued that a MWP1a from Antarctica triggered the BA, as it led to an increase in the strength of North Atlantic deep water formation and thereby an increase in the AMOC. While the source of the meltwater is still debated in the literature (Bentley et al., 2010;Harrison et al., 2019), both reconstructions used in the present study prescribe the largest ice volume changes during MWP1a in the Northern Hemisphere (Figure 1a). A northern hemispheric meltwater release of more than 0.4 Sv is associated with a significant freshening of the North Atlantic and Arctic ocean (Peltier, 2005), a reduction of the deep water formation and the AMOC in all simulations. Hence, the North Atlantic realms cool significantly and prevents a BA warming in our simulations (Figure 1 and Figure S7 in Supporting Information S1). These findings indicate that the prescribed location of the meltwater pulse in the underlying reconstruction is crucial for the modeled climate response and determines whether a BA warming may occur in deglacial simulations. After the end of MWP1a, the AMOC recovers in the model. In the simulations with GLAC-1D boundary conditions, this AMOC recovery is accompanied with an overshooting of the AMOC, peaking after 200-300 years. While the physical mechanisms of the AMOC slowdowns due to MWP1a in the different simulations are identical, the timing and magnitude of the response is dependent on the underlying reconstruction (Figures 1a-1d; Figure S7 in Supporting Information S1). The peak freshwater discharge of MWP1a occurs about 100 years earlier in GLAC-1D and is about 0.15 Sv weaker than in ICE-6G.
After the AMOC recovery in response to MWP1a and before MWP1b, three additional AMOC slowdowns occur. However, they either occur in simulations with GLAC-1D (about 13.45 ka) or simulations with ICE-6G boundary conditions (about 12.85 ka and 11.75 ka). Note, that the AMOC does not slow down in Glac1D_P1 at 13.45 ka and Ice6G_P1 at 11.75 ka, as model P1 is less sensitive to freshwater changes (see earlier discussion about instability). The AMOC slowdowns are associated with enhanced freshwater discharge into the Arctic due to changes in river directions (Figures 1, 2e, and 2f). Small changes in the topography prescribed from the reconstructions allow a rerouting of meltwater that stems from the Laurentide ice sheet through the MacKenzie instead of the Mississippi river (Figures 2e and 2f). The rerouting leads to a freshening of the Arctic, suppressed overturning and a decrease in AMOC strength. The 11.75 ka event in the simulations with ICE-6G show only small changes in the river directions, but ice-sheet melt occurs mainly at the northeastern side of the Laurentide ice sheet (not shown). This leads to discharge of meltwater into the Labrador Sea and Arctic, hence, a decrease in the AMOC strength. The 12.85 ka event in the simulations with ICE-6G coincides with the YD. Our simulations confirm findings by Condron and Winsor (2012), who suggested that a rerouting of freshwater into the MacKenzie valley led to the AMOC weakening during the YD, as coastal boundary currents effectively transport freshwater from the MacKenzie river toward the deep water formation sites of the subpolar North Atlantic. As no significant increase in meltwater release is evident in ICE-6G during the YD, the results are also in line with previous studies, indicating that no meltwater pulses initiated the YD (Abdul et al., 2016;Stanford et al., 2006;Tarasov & Peltier, 2005).
A meltwater discharge of about 0.3 Sv at 11.45 and 11.25 ka in ICE-6G and GLAC-1D, respectively, characterizes MWP1b (approx. 11.45 ka; Abdul et al., 2016) in the simulations. However, uncertainties of the origin of MWP1b exist (e.g., Bentley et al., 2010;Harrison et al., 2019). In GLAC-1D the largest ice melt during MWP1b occurs in the Northern Hemisphere, while in ICE-6G melt is dominant in the Southern Hemisphere ( Figure 1a). In the simulations with GLAC-1D, meltwater discharge into the North Atlantic reduces the North Atlantic deep water formation and leads to a slowdown of the AMOC around the peak of the meltwater discharge ( Figure 1). This results in a significant cooling of the Northern Hemisphere, specifically over the North Atlantic (Figures 1a-1d; Figure S7 in Supporting Information S1). In the simulations with ICE-6G, meltwater is released into the Southern Ocean and associated with a cooling of the Southern Hemisphere (Figure 1a and Figure S7 in Supporting Information S1). Interestingly, meltwater release into the Southern Ocean does not significantly affect the AMOC (Figure 1b), which is in line with findings by Stouffer et al. (2007). It also confirms findings by Stanford et al. (2006), who argued that the location of the meltwater injection is more important than the rate and magnitude of the meltwater injection for the North Atlantic deep water formation. The significantly different location of the freshwater release prescribed by the reconstructions explains the fundamental differences in the AMOC response to MWP1b in the simulations with GLAC-1D and ICE-6G reconstructions.
The results emphasize that the climate variability in transient deglacial simulations is largely dependent on the underlying ice-sheet boundary conditions. The height of the reconstructed ice sheets has an influence on the atmospheric circulation during the LGM and early in the deglaciation until about 14 ka, when northern hemispheric ice sheets were still extensive, affecting the wind-driven overturning and the AMOC strength. However, the results also indicate that the meltwater history and changes in river directions dominate the millennial-scale variability in the simulations during the last deglaciation. The reconstructions determine the magnitude, rate and location of meltwater injections, hence, significantly impact the deep water formation in the North Atlantic and the AMOC. These in turn have a significant influence on the northern hemispheric climate (see Figure S7 in Supporting Information S1). Additionally, the model version significantly controls how large the model responds to changes in the boundary conditions.

Implementation of Meltwater Release Determines Millennial-Scale Variability
The PMIP4 protocol also permits different implementations of meltwater distribution in the experiments (Ivanovic et al., 2016). Sensitivity experiments with a global distribution of meltwater and no meltwater show that the multi-millennial climate variability during the deglaciation is largely controlled by the choice of meltwater implementation (Table 1 and Figures 1e-1h). If meltwater is distributed globally, most AMOC slowdowns throughout the deglaciation still occur (MWP1a, MWP1b, YD). However, a global distribution of meltwater reduces the amount of freshwater that reaches the deep water formation sites in the North Atlantic, resulting in a weaker AMOC response in Ice6G_P2_Glob as compared to Ice6G_P2 ( Figure S8 in Supporting Information S1). This indicates that the prescribed meltwater pulses of the aforementioned events are large enough to reduce the overturning strength independent of where freshwater is added in Ice6G_P2_Glob (e.g., Smith & Gregory, 2009;Stanford et al., 2006). It is noteworthy, that the magnitude of the AMOC slowdown during the YD is similar in Ice6G_P2_Glob and Ice6G_P2 but the timing of the AMOC slowdown differs (Figure 1). This indicates that changes in the river directions are likely contributing to an early initiation of the YD AMOC slowdown in Ice6G_P2. However, the persistence of a meltwater peak of more than 0.1 Sv over more than 500 years prescribed from ICE-6G seems large enough to affect the deep water formation in Ice6G_P2_Glob. For other events the AMOC response is fully absent in Ice6G_P2_Glob. For example, the 11.75 ka event is not associated with a longlived meltwater peak in ICE-6G and the AMOC response in Ice6G_P2 is mainly a result of the discharge into the Labrador Sea and Arctic, hence, the discharge location (see earlier discussion). This emphasizes that a global distribution of relatively small amounts of meltwater does not significantly affect the AMOC in Ice6G_P2_Glob (Figure 1; see also Manabe & Stouffer, 1997;Smith & Gregory, 2009;Rahmstorf, 1996). Additional differences in the AMOC response throughout the deglaciation may arise from the high non-linearity of the AMOC response to changes in greenhouse gases and temperatures (Stouffer & Manabe, 2004;Wang et al., 2002), hence, at times when the CO 2 concentration increases rapidly (see Figure S9 in Supporting Information S1).
In Ice6G_P2_noMW, multi-millennial scale variability is largely absent and the global mean sea-surface temperature is warmer than in Ice6G_P2 from about 16.5 ka onward ( Figure 1). As no meltwater pulses are prescribed in Ice6G_P2_noMW that could potentially lead to AMOC slowdowns, an increase in North Atlantic sea-surface temperatures is evident around the BA. The North Atlantic warming is associated with a reorganization of the deep water formation sites throughout the deglaciation as well as a northward shift of the jet stream (see Section 3.1).

Conclusions and Implications for PMIP4
Underlying ice-sheet reconstructions and the implementation of meltwater distribution are crucial for the modeled climate response in transient simulations of the last deglaciation. These results are obtained by performing an ensemble of eight transient simulations with the state-of-the-art MPI-ESM in coarse resolution. MPI-ESM-CR is most suitable for transient simulations, as it accounts for changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask. The objective of this study is to compare the climatic impact from the different methods proposed in the PMIP4 deglaciation protocol (Ivanovic et al., 2016) and to highlight the uncertainties that arise from different boundary conditions. A more thorough analysis of the associated physical processes is essential and will be conducted in a future study.
Specifically, we find that differences in the topography between ICE-6G and GLAC-1D significantly impact the atmospheric circulation. With glacial ICE-6G ice sheets the Atlantic jet in the Northern Hemisphere shows higher wind speeds and a more zonal and equatorward position (e.g., Merz et al., 2015). Glacial GLAC-1D ice sheets are lower and do not lead to changes in the jet position. Hence, the ocean response significantly differs between the simulations with ICE-6G and GLAC-1D. Small changes in the topography throughout the deglaciation affect changes in river directions and explain substantial differences in the climate response between simulations. For example, during the YD a rerouting of meltwater from the Mississippi into the Arctic is evident in the simulations, but only ICE-6G reconstructions prescribe a sufficient amount of melt water to significantly affect the ocean circulation. Therefore, a YD cooling is only evident in simulations with ICE-6G. Other observed climate events, such as the BA warm period, are not simulated in the experiments with prescribed meltwater, as MWP1a is released in the Northern Hemisphere in both reconstructions, where it leads to a reduction in the AMOC strength and a significant cooling. Sensitivity experiments with different model versions show that the background climate is essential for the ocean response in the simulations, as it determines how close the AMOC is to the bifurcation point between a strong and weak AMOC (Klockmann et al., 2018;Oka et al., 2012). These findings are in line with previous studies, investigating the effect of ice-sheet boundary conditions on the atmospheric and oceanic circulation (e.g., Bakker et al., 2020;Löfverström et al., 2014;Merz et al., 2015;Pausata et al., 2011;Ullman et al., 2014) and the importance of meltwater injections for the AMOC stability (e.g., Stanford et al., 2006;Stouffer et al., 2007). Performing a first systematic ensemble of transient simulations with different ice-sheet boundary conditions and methods of meltwater distribution, the present study extends previous studies and shows that differences in the topography and meltwater history of ICE-6G and GLAC-1D dominate the simulated millennial-scale climate variability in simulations of the last deglaciation.
The PMIP4 deglaciation protocol is designed for a large range of models to participate in the exercise. Most modeling groups will likely only contribute to PMIP4 with a subset of the simulations presented here, due to computational or technical limitations. Our experiments point toward the challenges in interpreting results from simulations with the different methods proposed in the PMIP4 deglaciation protocol, as the differences in the climate response due to the implementation choices proposed in the protocol can be as large as the climate variability that the simulations try to capture. Further, a direct comparison of the model results with proxy evidence will be difficult, as the proposed boundary conditions do not allow for the modeling of all observed climate events (e.g., BA), and the climate effect arising from the uncertainties in the ice-sheet reconstructions is large. This points toward the necessity of a more process-based analysis and interpretation of the PMIP4 last deglaciation ensemble.

Data Availability Statement
Model data and scripts used for the analysis will be available online on the MPG.PuRe repository (http://hdl. handle.net/21.11116/0000-0009-128D-4) upon publication. Post-processed Paleoclimate Modeling Intercomparison Project -Phase four model data by LSCE is available at http://dods.lsce.ipsl.fr/pmip4/db/ (last access: 13 December 2021). Max Planck Institute Earth System Model model code is available upon request from the Max Planck Institute for Meteorology (reinhard.budich@mpimet.mpg.de) under the Software License Agreement version 2 (https://mpimet.mpg.de/en/science/modeling-with-icon/code-availability -last access: 25 August 2021).
Additionally, the full model output will be made available on the ESGF repository. Planck Institute for Meteorology Earth System Model simulations were performed at the German Climate Computing Center. The authors would also like to thank the Laboratoire des Sciences du Climate et de l'environment for providing post-processed Paleoclimate Modeling Intercomparison Project -Phase 4 data and acknowledge two anonymous reviewers as well as the editor Gudrun Magnusdottir, who have helped to improve our manuscript. Open access funding enabled and organized by Projekt DEAL.