Coupled climate-ice sheet modelling of MIS-13 reveals a sensitive Cordilleran Ice Sheet

Previous modelling efforts have investigated climate responses to different Milankovitch forcing during Marine Isotope Stage (MIS) 13. During this time the climate has been highly variable at atmospheric CO2 concentrations of ~240 ppm. As yet, ice sheet-climate feedbacks were missing in previous studies. Therefore we use the state-ofthe-art coupled climate-ice sheet model, AWI-ESM-1.2, to investigate the MIS-13 climate and corresponding Northern Hemisphere ice sheet (NHIS) evolution by performing simulations under three different astronomical configurations representing 495, 506 and 517 kyr BP. The simulated excess ice compared to present-day is mainly over the Cordillera, Arctic islands and Tibet. The global mean surface air temperature for the MIS-13 experiments have the same magnitude. At 506 kyr BP with boreal summer at perihelion, the Northern Hemisphere continents are warmer during summer than the other experiments, which could potentially inhibit the development of the ice sheets. The Cordilleran Ice Sheet is found to be especially sensitive to orbital (precession) forcing, at an intermediate CO2 level. This is probably due to its high elevation where the freezing point could be easily maintained. The other ice sheets over northeast America and Eurasia, however, are absent in our simulations. We propose that the alpine-based Cordilleran Ice Sheet is more sensitive and easier to build up than other NHISs in response to the astronomical controlled summer insolation. Dynamic surges are simulated for the Cordilleran Ice Sheet under fixed low orbital forcing. These surges due to internal ice sheet-climate feedbacks could potentially be the mechanism for the millennial scale H-like events.


Introduction
The Northern Hemisphere summer insolation is recognized as a major control of the glacial-interglacial cycles during the Quaternary Period (Milankovitch, 1941;Imbrie et al., 1993;Huybers and Wunsch, 2005). The sawtooth pattern of ice volume change indicates that internal nonlinear feedbacks between ice sheets and other climate components are non-negligible (Lisiecki and Raymo, 2005). Contributions from different processes (like the ice-albedo feedback, elevation-mass balance feedback, ice-solid Earth feedback, etc.) to ice sheet evolution likely vary between different glacial-interglacial stages (e.g. Abe-Ouchi et al., 2007;Kubatzki et al., 2006). However, the individual or combined effects from these feedbacks are still not fully understood, and are difficult to quantify in a stand-alone ice sheet or climate model. This is due to the circularity between the modelled climate and ice sheets, where the climate forcing is generated from a prescribed ice sheet configuration and vice versa. A coupled climate-ice sheet model needs to be used to properly represent these internal feedback correctly.
Modelling interactive ice sheets with other Earth system components has been implemented to date mostly with Earth systems Models of Intermediated Complexity (EMICs) due to their efficiency in numerical computation, especially on paleo timescales Charbit et al., 2005;Fyke et al., 2011;Stap et al., 2014;Heinemann et al., 2014;Choudhury et al., 2020). The computational costs limit the applications of more advanced General Circulation Models (GCMs) to millennial timescales or to equilibrium simulations at certain time-slices (Gregory et al., 2012;Lipscomb et al., 2013;Vizcaino et al., 2015;Ziemen et al., 2019). Most of those studies focused on either glacial inception and termination within a glacial cycle, or future scenarios with anthropogenic forcing. Older interglacials like Marine Isotope Stage 13 (MIS 13, ~530 to ~480 kyr BP, Railsback et al., 2015) received little attention for ice sheet-climate modelling.
Extraordinary features during MIS-13 are exhibited in the benthic δ 18 O records from deep-sea sediments, and temperature reconstruction and Greenhouse gas (GHG) records from Antarctic ice cores (Lisiecki and Raymo, 2005;EPICA Community Members, 2006;Jouzel et al., 2007;Siegenthaler et al., 2005). The amplitude of the glacial cycles seems to be smaller before ~420 kyr BP (including MIS-13) than after. The MIS-13 climate was likely less warm than more recent interglacials, and thus the period is considered either as a long-lasting cold interglacial (e.g. Yin et al., 2009;Rachmayani et al., 2016), or one interglacial and several interstadials (e.g. Past Interglacials Working Group of PAGES, 2016; Tzedakis et al., 2017).
Regardless of the definition, MIS-13 can be split into three sub-stages (MIS-13 a-c), with corresponding extreme marine isotopic events . (Local maxima and minima in δ 18 O values are notated with decimals; Railsback et al., 2015). Furthermore, the MIS-13a can also be divided into three sub-stages, with related marine isotopic subevents MIS-13.11, 13.12, 13.13 (Fig. 1;Yin et al., 2009). The astronomical configurations during different sub-stages were different, causing direct or indirect effects onto local climate (Yin et al., 2009;Rachmayani et al., 2016). The Greenhouse gas (GHG) values during MIS-13 were relatively low compared to the preindustrial, where CO 2 was around 240 ppm (Siegenthaler et al., 2005, Fig . 1). Meanwhile, the sizes or locations of the ice sheets are still unknown, since larger ice sheets than present day may have existed due to the cooler climate. The presence of excess ice could potentially induce additional feedbacks within the climate system, e.g. modify the Earth's energy balance via changes in albedo, and also the atmospheric circulation via changes in topography (Abe-Ouchi et al., 2007).
Various modelling studies but with prescribed ice sheets have been conducted to infer the background climate during MIS-13, with either EMICs (e.g. ) or GCMs (e.g. Muri et al., 2012Rachmayani et al., 2016). In Rachmayani et al. (2016), the intrainterglacial climate variability and associated astronomical forcing is studies during interglacials including MIS-13 with the GCM CCSM3 and fixed preindustrial ice sheet condition. Larger Northern Hemisphere ice sheets (NHIS) are proposed to have existed at that time according to some studies (Yin et al., , 2009Sundaram et al., 2012;Muri et al., 2012Muri et al., , 2013Yin et al., 2014;Karami et al., 2015), and could have nonnegligible influence on monsoons and precipitation distribution. By using the EMIC LOVECLIM,  found that in addition to the astronomical forcing at perihelion, the existence of Eurasian ice sheets could also reinforce the East Asian Summer Monsoon (EASM) through an atmospheric wave train. Yin et al. (2009) tested the individual and combined effects of ice sheets and precession on MIS-13 climate, by prescribing 5 different ice sheet configurations and 2 opposing precessional angles. Similar experiments were also conducted with more advanced GCMs (Muri et al., , 2013, with more focus on the exceptionally strong EASM (Chen et al., 1999;Yin and Guo, 2008;Yin et al., 2014;Karami et al., 2015). However, Shi et al. (2020a) recently show that a large Antarctic ice sheet but removed Greenland ice sheet could also induce a strong EASM. The removed Greenland ice sheet could induce additional atmospheric wave train. The larger Antarctic ice sheet could result in stronger upwelling of circumpolar deep water, warmer SST in the Southern Ocean which propagates to the North and influences the EASM.
The ice sheet sizes or extents during MIS-13 is difficult to infer with geomorphological reconstruction, since the repeated glaciations erode the evidence. Instead, some studies assume a simple linear relationship between total ice volume change on the Earth and deep-sea benthic δ 18 O values (e.g. Yin et al., 2009). Under this premise, the size and locations of the excess ice compared to present day were inferred based on the ice sheet configuration at the Last Glacial Maximum (LGM), i.e. the ice volumes of different ice sheets at MIS-13 are proportional to that at the LGM. However, the nonlinear response of ice sheets to climate change and additional variations from the deep-sea temperature on the δ 18 O are neglected.
In light of the important role of ice sheets on MIS-13 climate, Earth system modelling with interactive ice sheets should be applied for testing the possible configuration of NHISs, as well as a re-evaluation of the related climate conditions. In this paper, we make a first attempt of simulating MIS-13 climate with a sophisticated coupled GCM-ice sheet model, AWI-ESM-1.2 with high spatial resolution. Three time slices during MIS-13 with different orbital configurations and one preindustrial (PI) control run are simulated with asynchronous climate-ice sheet coupling. The Northern Hemisphere ice sheet evolution and the effects from different astronomical configurations and lower GHG concentrations under these circumstances are discussed. Since the NHIS are mainly land-based ice sheets, the atmosphere and ice sheet components are the primary focus of this study.

The AWI-ESM with interactive ice sheets
The earth system model we use here is the state-of-the-art model AWI-ESM-1.2, with the GCM AWI Climate Model (AWI-CM) asynchronously coupled to the Parallel Ice Sheet Model PISM (version 1.1) . The coupling between the ice sheets and the climate are only applied to the Northern Hemisphere in this study. The Antarctic ice sheet is kept constant at present-day configuration.
AWI-CM is a coupled atmosphere-ocean model developed at Alfred Wegener Institute (Rackow et al., 2018;Sidorenko et al., 2015), consisting of the atmosphere model ECHAM6 and the ocean model FESOM1.4. ECHAM6 is the sixth generation of the ECHAM (short for ECMWF and Hamburg) atmospheric general circulation models (Stevens et al., 2013). It incorporates diabatic processes with large-scale  Berger (1978): Obliquity (light blue), Eccentricity (brown), and Precession parameter (green). Benthic δ 18 O from Lisiecki and Raymo (2005, black). δD from Antarctic Dome C (Jouzel et al., 2007, blue). CO 2 record from Hönisch et al. (2009, purple). The round dots indicate the respective values for preindustrial (PI). The shaded yellow area indicates the period of MIS-13. The crosses on the δ 18 O curve are the selected peaks (marine isotopic events 13.12, 13.13, 13.2); The vertical dashed lines are the times used in the MIS-13 experiments. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) circulations of the atmosphere, and contains the JSBACH land vegetation model that accounts for the land processes, and a hydrological model. The resolution is T63 (~1.9 ∘ ) with 47 vertical levels, resolving the atmosphere at up to 0.01 hPa. FESOM1.4 (Finite Element Sea ice-Ocean Model) is a global sea-ice ocean model on an unstructured mesh with multi-resolution modelling functionality (Danilov et al., 2004;Wang et al., 2014). With the unstructured mesh, high spatial resolution can be applied in dynamically active regions while a relatively coarse resolution can be used elsewhere to reduce computational demands. The resolution ranges spatially from 10 km to 180 km, and 46 vertical z-level grids (Fig. A1). The exchange between ECHAM6 and FESOM1.4 is through the coupler OASIS-MCT. AWI-CM has comparable performance to other GCMs, when comparing the present day state (Danilov et al., 2004;Wang et al., 2014), paleoclimate scenarios like the LGM and the Holocene and future scenarios (Brown et al., 2020;Lohmann et al., 2020;Shi and Lohmann, 2016;Shi et al., 2020b;Ackermann et al., 2020).
PISM is a three-dimensional thermo-mechanically coupled shallow ice sheet model (the PISM authors, 2015;Bueler and Brown, 2009;Winkelmann et al., 2011). The shallow ice approximation and shallow shelf approximation are used for computing the ice sheet stress balance. The ice deformation is a function of temperature, pressure and liquid water fraction (called the Glen-Paterson-Budd-Lliboutry-Duval flow law; Aschwanden et al., 2012). Enhancement factors of 2.0 and 1.0 are applied for SIA and SSA respectively (Cuffey and Paterson, 2010). Glacial isostatic adjustment (GIA) is calculated using the Lingle and Clark method (Lingle and Clark, 1985). The spatial domain is defined on a Northern Hemisphere polar stereographic grid. The horizontal resolution is 20 km. The other model setup parameters of the ice sheet dynamics are similar to Niu et al. (2019).
Asynchronous iterative coupling is applied between AWI-CM and PISM due to the computional time limits, with 1-year of AWI-CM for every 100-years of PISM. (The resolution of the AWI-CM is higher than other asynchronous coupled studies and the simulations are computationally expensive. It takes ca. 1 h per climate model year.) The exchange frequency is reasonable due to the different response times among Earth system components (Ruddiman, 2001). The atmosphere, the land and the ocean surface have rapid response. The response timescales to perturbations are within tens of years, while it is thousand of years for deep ocean and tens of thousands of years for large-scale ice sheets. As a result, the climate model with fast components can be integrated for a shorter period, and applied for a longer period in the ice sheet model. We note here that the 1-year climate interval may introduce large interannual variability (climate noise) that may hamper or accelerate the ice sheet build-up. However, Pollard et al. (1990) found the effects to be minor when applied to simple flowline models for climate-ice sheet time step ratios of 1:100 or greater. We tested initialization with a ratio of 5:100, and the initial build-up patterns of the ice sheets are similar. The simulations are intended to integrate until quasiequilibrium, and the influence to the equilibrated patterns is minor. The AWI-CM is firstly integrated for 1 year with prescribed ice sheet boundaries. Then the ice sheet model PISM is switched on and integrated for 100 years, forced by the output from AWI-CM. Afterwards, the simulated ice sheet state is used as a new boundary condition to AWI-CM, and the AWI-CM is integrated for another 1 year. As such, the simulations between AWI-CM and PISM switch back and forth until they reach relatively stable states.
Variable exchanges between the models are illustrated in Fig. 2. For PISM, monthly mean surface air temperature and precipitation fields from ECHAM6 are transfered as atmospheric forcing to PISM, with surface mass balance calculated by a positive-degree-day (PDD) scheme (Reeh, 1989;. Ocean temperature and salinity from FESOM1.4 (average from 0 m to 500 m in depth) are transfered as ocean forcing for deriving the basal melt rate and sub ice-shelf basal temperature (Holland and Jenkins, 1999). For AWI-CM, the updated albedo, orography and freshwater discharge from PISM are used as new boundary conditions in AWI-CM. The ice volume changes due to ablation, calving, and basal melting under grounded ice or ice shelves are all considered as freshwater discharge into the ocean via the hydrological model in ECHAM6. The change of bathymetry in FESOM1.4 due to the migration of the grounding line or GIA is not technically applicable in the current model. Because FESOM1.4 has an unstructured mesh, changing in bathymetry would require a change of the ocean mesh which is not yet feasible. Bilinear interpolation is used to account for the mismatch of the resolution between AWI-CM and PISM. The coupling scheme is described by  for the PI condition and has been applied in future climate scenarios (Ackermann et al., 2020).

Experimental design
We conduct four experiments. A preindustrial scenario is used as a reference (PI-ref), and the other experiments focus on MIS-13 climate with three fixed orbital configurations (495, 506, 517 kyr BP; Table 1). The summer insolation at 65 ∘ N was at a local maximum at 506 kyr BP, with boreal summer at perihelion (Fig. 1, red line). At 495 kyr BP and 517 kyr BP, summer insolation at 65 ∘ N reached two local minimums, when boreal summer were both at aphelion. The obliquity was the highest at 495 kyr BP, the lowest at 517 kyr BP, and was intermediate at 506 kyr BP (similar to PI). The GHG values are kept the same among the MIS-13 experiments and are the same as that in previous studies (Table 1, Yin et al., 2009;Muri et al., 2012Muri et al., , 2013, where CO 2 concentration is lower than that in PI-ref (240 ppm versus 280 ppm). The overall experiment settings for MIS13-495, MIS13-506 and MIS13-517 are corresponding to marine isotopic events at MIS 13.12, MIS 13.13 and MIS 13.2 respectively (Fig. 1, black line). The land-ocean map used in the MIS-13 experiments is the same as the PI condition.
The experiment PI-ref is initialized from an equilibrium PI simulation without ice sheet dynamics (AWI-CM alone, integrated for 2961 years), then it is integrated with dynamic ice sheets. The MIS-13 experiments are first initialized from an AWI-CM standalone simulation under fixed boundary conditions at 506 kyr BP (745 model years). Then, the

Table 1
Orbital parameters and greenhouse gas values used in the experiments, global annual mean surface air temperature (GSAT) and the simulated ice volume (iV, in sea level equivalent, m SLE, calculated with 20,000-year mean) for all Northern Hemisphere ice sheets (All), and the separated Greenland ice sheet (GIS), Cordilleran ice sheet (CIS), Tibet, and other regions.

The MIS-13 climate
3.1.1. The MIS-13 climate relative to the preindustrial Although forced with different orbital configurations, the simulated global annual mean surface air temperature (GSAT) for the MIS-13 experiments are around the same magnitude (~11.8 ∘ C, Table 1), and are about 1 ∘ C cooler than that of PI-ref (12.795 ∘ C). The differences in GSAT could largely attribute to the difference in GHG concentrations (Arrhenius, 1896; Yin and Berger, 2012).
However, the spatial and seasonal patterns from different MIS-13 experiments differ in response to different astronomical forcing (Fig. 3). For MIS13-495, the averaged zonal mean temperature deviations relative to PI-ref in different seasons is the smallest among the MIS-13 experiments (Fig. 3a-c). In experiment MIS13-506, significant warming in summer (JJA) and cooling in winter (DJF) are exhibited over most of the continents (Fig. 3d-f). The summer monsoon in tropical Africa and India is intensified due to higher land-sea thermal contrast ( Fig. A2e-f). These features are likely due to the precessional effect, since significant warming over the Northern Hemisphere continents occurs in experiment MIS13-506 when boreal summer is at perihelion (but not in MIS13-495 with high obliquity). Lastly, the obliquity is the lowest in experiment MIS13-517. The cooling over high-latitudes is found to be the largest, which can be up to 8 ∘ C in both winter and summer (Fig. 3g-i).

Impact of orbital changes at 506 and 495 ka with significant precession differences
For the climatic effects induced by varied orbital forcing during MIS-13, we compared the differences between the MIS-13 experiments. For experiment MIS13-506 and MIS13-495, the difference in obliquity is relatively small (comparing to that between 517 kyr BP and 495 kyr BP). However, the precessional angles for MIS13-506 and MIS13-495 are at opposite phases (boreal summer at 506 kyr BP is at perihelion). This results in a significant increase in boreal summer insolation, and reduction in austral summer insolation (Fig. 4).
In boreal summer, the SAT increases over the Northern Hemisphere ocean and most of the continents response directly to the higher  insolation. The land-sea temperature contrast is strengthened, which leads to intensified summer monsoon. This could explain the significant cooling over tropical Africa and India, where increased monsoon rainfall leads to increased cloud cover and decreased short-wave radiation received by the surface (Rachmayani et al., 2016). The warming over northern mid-latitude continents is up to 8 ∘ C (Fig. 5b). During boreal winter, there is cooling over the mid-and low-latitudes and Antarctica (Fig. 5c). Warmer signals are present over northwestern North America, western Arctic ocean and some of the Southern Ocean in MIS13-506 compared to MIS13-495. The warming over the ocean is mainly due to reductions of sea-ice extent in JJA that persist through DJF. The warming over northwestern North America is mainly a local response from the absence of Cordilleran Ice Sheet (CIS, Section 3.2.1).

Impact of orbital changes at 517 and 495 ka with significant obliquity differences
The boreal summer at 517 kyr BP and 495 kyr BP are both at aphelion, while obliquity is the lowest and highest respectively. Lower obliquity reduces the insolation received over polar regions during both hemispheric summers, while there are increases in the insolation received during winter especially for the Southern Hemisphere (Fig. 4b).
The SATs over land respond faster than the ocean due to lower heat capacity. In other words, the continents are in general cooler (warmer) in summer, and warmer (cooler) in winter over high (low) latitudes for MIS13-517 than in MIS13-495 ( Fig. 5d-f). However, consistent cooling (warming) is shown in high (low) latitudes of the ocean all year round. This is probably due to the larger heat capacity of the ocean and the residual effect of sea ice that lasts from previous summer. In addition, no significant change in precipitation is shown especially over highlatitudes, which is slightly different from the results in Rachmayani et al. (2016).

The different MIS-13 climate states under different orbital forcing
From our coupled climate-ice sheet model experiments, increased ice compared to present day are mainly simulated over the Cordillera, Arctic islands or Tibet. The global patterns simulated for the three MIS-13 climates are remarkably similar as in previous studies with different levels of complexity of the climate models (Yin et al., 2009;Yin and Berger, 2012;Muri et al., 2013;Rachmayani et al., 2016). Particularly Rachmayani et al. (2016) investigated different interglacials with varies time slices including MIS-13 with a GCM, and the different effects of orbital configurations and GHG were studied. Our results support their findings: GHG plays a dominant role on the variations of the global annual mean temperature. The surface air temperature pattern can largely be interpreted by direct radiative forcing. When boreal summer is at perihelion (aphelion), a global scale warming (cooling) especially over the Northern Hemisphere mid-latitude continents is shown in boreal summer. Lower obliquity are generally characterized by cooling over the Northern Hemisphere high-latitudes and slightly warming over low-latitudes in boreal summer. Feedbacks between different earth system components also modify the regional climate significantly (e.g. monsoon regions or polar oceans).

The equilibrium response under different orbital forcing
Significant increases in the simulated ice sheets are shown in experiment MIS13-517 and MIS13-495 (Fig. A3, where large differences than MIS13-506 and PI-ref are exhibited after 5000-year integration). The response of ice sheets in MIS13-517 is slightly faster than MIS13-495. Although the GSAT is lower in experiment MIS13-506 than in PIref, the simulated ice volume of NHIS is slightly smaller (Table 1). It is likely due to the precessional effect that significant warming in summer occurs over the Northern Hemisphere continents in response to the higher summer insolation than the other MIS-13 experiments. Comparing to PI-ref, more ice is built up over the North American Cordillera, polar islands (e.g. Arctic Archipelago, Svalbard) and Tibet in experiments MIS13-495 and MIS13-517 (Fig. 6) in response to the cooler summer temperature.
The ice volume over Greenland, Cordillera, Tibet and other polar regions are calculated separately (Table 1). The simulated changes of Greenland ice volume in these 4 experiments are small (~8.7 m SLE). Most differences are in the Cordillera, Tibet and the polar islands. In the PI-ref experiment, there is a small but significant growth of the Cordilleran Ice Sheet while only glaciers exist today (The locations of the CIS matches with that of the glaciers). One possible reason of the larger ice volume could due to a cold bias in the GCM for that region. On the other hand, the ice caps were found expanding during the PI period (Luckman, 2000), so the climate is likely very close to the tipping point towards the buildup of ice. In the absence of global warming, it is possible that an ice sheet would have eventually developed. Nevertheless, significant deviations of simulated CIS are still shown under different MIS-13 orbital conditions. The simulated CIS is the smallest in MIS13-506 (1.27 m SLE), which is around half the size as PI-ref. A much larger CIS is simulated in MIS13-495 and MIS13-517 (more than 8 m SLE), consistent with lower summer insolation and GHG. A similar feature is also shown for ice sheets over the polar islands. In addition, no ice sheet is simulated over Tibet in experiments PI-ref and MIS13-506, while large amounts of ice sheets are simulated in MIS13-495 and MIS13-517. There are different hypotheses regarding whether or not there was a coalescence of one large Tibet ice sheet during glaciations (e.g. Kuhle, 2005;Derbyshire et al., 1991). It is also possible that the simulated ice sheet over Tibet results from a cold bias of that region in the GCM, or because the spatial resolution is not high enough for resolving the steep and mountainous topography there. However, we will not expand this here since Tibet is not the focus of the current study. In addition, direct short-wave radiation effect on the computed surface mass balance over ice sheets might not be fully represented by the PDD scheme, which could also induce bias on the simulated results.

The transient response of the Cordilleran Ice Sheet
Apart from the mean state, periodic surge oscillations are found over northwest of CIS in MIS13-495 and MIS13-517 when it reaches to full glacial states (Fig. 7a). In the simulations, the north of the CIS is slowly built up over southern Alaska, advancing northwest to the ocean in around 20,000 years ( Fig. 7b-c). When it reaches the ocean, a surge process is triggered and the CIS retreats rapidly back to a minimum state within 5000 years (Fig. 7d). The overall buildup and retreat process resembles a saw-tooth pattern with a period of around 25,000 years. The difference in ice volume between the maximum and minimum states can be up to 5.2 m SLE. In order to uncover the physical mechanism of the dynamic surge, the mean values of the atmosphere, ocean and ice sheet conditions over the northwest of the Cordilleran Ice Sheet (hatched area in Fig. 7b-d) are calculated (Fig. 8). From Fig. 7b-d, the CIS advances slowly towards the northwest and covers southern Alaska before a rapid retreat is triggered at the land-ocean boundary. From variables in Fig. 8, the different climate components couple closely with each other and respond relatively synchronously, except that the bedrock topography declines monotonically before the retreat. Thus, one possible mechanism is suggested as follows: When the bedrock is depressed and reaches to a critical threshold (e.g. ~75 m at the start of the first retreat around 30,000 years), part of the bedrock goes below sea level. This allows the ocean to interact with the northwest marine edge of the ice sheet, causing calving or oceanic melting of ice shelves, leading to a positive feedback for accelerating retreat. In the end, an ice sheet state is reached with minimized ocean influence, after which the bedrock is uplifted slowly towards its original state. Subsequently, another ice sheet buildup processes proceeds until another surge process is triggered.
The temperature-elevation change acting as a positive feedback could also speed up the retreat process. Specifically, the surface temperature (especially summer) increases due to the thinning of the ice, and more melting is induced. Conversely, during the ice sheet retreat, more freshwater goes into the ocean and causes decreases of the ocean temperature and salinity. This could act as a negative feedback and slow down the ice sheet melting, while is likely a secondary factor.

The sensitivity of the Cordilleran Ice Sheet to different orbital forcing
Summer surface air temperature is considered as a primary factor for the ice sheet evolution of the Northern Hemisphere (e.g. Gallée et al., 1992;Niu et al., 2019). Correspondingly, the ice volume of the CIS in MIS13-506 is less than half of that in PI-ref (Table 1), showing that the insolation effect here is larger than the GHG effect. In the MIS-13 experiments, MIS13-495 is with the highest obliquity and MIS13-517 is with the lowest obliquity, while the boreal summer are both at aphelion. The simulated CIS in both cases reach to a full glacial state. However, for MIS13-506, with obliquity in the middle but boreal summer is at perihelion, the simulated CIS is much smaller. This shows that precession plays an important role on the significant reduced CIS. The changes in ice volume of the Greenland Ice Sheet due to insolation change is limited at current GHG level. From this results, we propose that the CIS is very sensitive and has rapid response to summer insolation.
Geological reconstruction of the CIS during MIS-13 is difficult, due to the repeated glaciations eroding the evidence. Evidence of sub-glacial volcanic eruptions in southwest British Columbia suggests multiple glaciations that do not align with the MIS based glacial-interglacial timeline, including MIS 13 and 15 (Wilson and Russell, 2018;Wilson et al., 2019). At least eight glacial advances are recognized around Gulf of Alaska at the end of the mid-Pleistocene transition, indicating a more dynamic CIS than recent glaical cycles (Montelli et al., 2017). Although fragmentary, abundant evidence suggest the CIS formed and decayed many times during the Pleistocene (Clague and Ward, 2011;Booth et al., 2003). In addition, various glacial deposits indicate a MIS-4 glaciation over central, northern or southern sector of the CIS domain and a deglaciated MIS-3 (Clague and Ward, 2011;Ward et al., 2007;Stroeven et al., 2010Stroeven et al., , 2014Batchelor et al., 2019;Wilson et al., 2019), if the last glacial cycle is treated as an analogy. During MIS-3 and MIS-4, the GHGs level were both low, while the insolation values differ. These support our results that a full glacial CIS could be reached when summer insolation (mainly precession induced) is low and the GHG is at intermediate values.
We propose that the alpine-based Cordilleran Ice Sheet is more sensitive and easier to build up than other NHISs in response to the astronomical controlled summer insolation. The initial expansion of ice sheet first requires a strong snow-albedo feedback. For the CIS, a freezing point could be easily reached due to its high topography. With lowering snowline elevations intersecting the high topography, ice caps grow and are accelerated by a temperature-elevation feedback. However, for the Laurentide or Eurasian Ice Sheets, their topographic seed regions (Baffin Island ranges, Norwegian highlands) are lower (and probably warmer) than that for CIS. Thus, stronger and more persistent cooling (e.g. lower greenhouse gas) is required for the initiation of these ice sheets. Once massive ice has been built, the temperature-elevation feedback could then play a role and help stabilize the ice sheets.

The absence of the Laurentide ice sheet or Eurasian ice sheets
In our simulations, the initial conditions (with Greenland ice sheet only) favor smaller ice sheets, while the simulated states leads to larger ice sheets. From the benthic δ 18 O records, the excess ice from presentday condition is estimated from 13 to 40 m SLE (Yin et al., 2009) (Also note that this is based on a simple linear assumption, and variations from the deep ocean temperature is not taken into account). In our experiments, the simulated ice sheets are mainly on Greenland (~8.7 m SLE), Cordillera (1.27 to 8.53 m SLE), Tibet (0 to 8.22 m SLE) and Arctic islands (0 to 1.55 m SLE). The simulated excess ice from Northern Hemisphere ice sheets compared to Present day potentially ranges from 4 to 21 m SLE. In addition, variations from the Antarctic ice sheet, though it responds in opposite precessional phases to Northern Hemisphere, could also potentially add up to 10 m or more SLE (on cancelling the precessional timescales in the benthic δ 18 O records; Raymo et al., 2006). Our results exhibit one possible climate-ice sheet configuration during the enigmatic MIS-13, with absent Laurentide ice sheet (LIS) and Eurasian ice sheets (EIS) but more extensive and dynamic CIS.
However, there is the possibility of the presence of LIS/EISs during this period. From the ocean sediment records that are located off the Newfoundland continental margin from Channell et al. (2012), detrital layers including signals of ice-rafting events are detected during not only glacial periods, but also interglacials such as MIS-5(a-d) and MIS-13. However, they also pointed that these signals are not necessarily only associated with Laurentide Ice Sheet surges, but also Northwest Atlantic Mid-Ocean Channel debris flows or glacial-lake drainage events. In addition, no ice-rafting event is identified farther in Atlantic ocean (Naafs et al., 2013); No significant ice sheet surges, i.e. no large-scale LIS/EIS, are recognized during MIS-13 (McManus et al., 1999;Bailey et al., 2010). In this case, there might be remnant ice of LIS/EIS transitioning from the previous deglaciation, but its size is probably very limited. Thus, we propose that our conclusion with sensitive CIS to insolation changes is still robust.
Furthermore, from the aspect of Earth system modelling, we speculate that a multi-stability of the climate-ice sheet system might exist which depends on the initial ice sheet configuration. From  with EMIC, the equilibrium response of the climateice sheet system has pronounced multi-stability and hysteresis behavior to a certain isolation forcing. Studies focusing on the last glacial transient simulations found that the bifurcation transitions from interglacial to glacial states are triggered by strong snow-albedo feedback . When an interglacial ocean-vegetation is prescribed, the initial ice sheet expansion over North America (LIS) could be suppressed (Kubatzki et al., 2006). As is shown from our simulations, if the LIS/EIS were completely deglaciated from MIS-14, the LIS/EIS is hard to be resumed during MIS-13. While if there is significant LIS/EIS, the climate-ice sheet system could possibly stabilize at a different state (due to the high topography and high albedo of the ice covered areas). The multi-stability of the climate-ice sheet system resulting from different ice sheet configurations could be tested by prescribing different sizes or locations of the ice sheets, orbital configurations and CO 2 values with our complex Earth system model. This requires substantial computing resources. This is interesting and worthy for a followup study. A complete picture of the climate-ice sheet condition during MIS-13 require further investigation from both geological reconstructions and model simulations. A transient simulation starting from the previous glaciation of MIS-13 would be necessary for determining the sizes, locations and thermodynamic structures of the ice sheets.

The transient response of the Cordilleran Ice Sheet
From Section 3.2.2, periodic surging occurs over the northeast CIS. This is simulated when the CIS reaches full glacial states, which happens in experiment MIS13-495 and MIS13-517. The analysis of the underlying mechanism points to the important role of the ocean, as well as the glacial isostatic adjustment process. Interactions between lagged bedrock rebound, ocean grounding-line depths/fluxes or oceanic melt can lead to internal oscillations (Pollard, 1983;Bassis et al., 2017;Kingslake et al., 2018). From Bassis et al. (2017) for the Laurentide Ice Sheet, the Heinrich events could be triggered by subsurface ocean warming and causing massive iceberg discharge from the Hudson Strait. Extensive West Antarctic Ice Sheet retreat during Holocene has also been detected, and the reverse of ice loss and re-advance of the ice sheet could result from the bedrock rebound-driven stabilizing processes (Kingslake et al., 2018). Increased oceanic influence could speed up the ice discharge process and induce a catastrophic ice sheet collapse and increased freshwater discharge. Meanwhile, the depressed bedrock favors the flow of the warm seawater towards the grounding line, which maintains the retreat, together with other feedback (e.g. elevationsurface mass balance feedback). Based on current analysis, though, another ice sheet oscillatory mechanism with free oscillations induced by basal melting below grounded ice cannot be ruled out (the binge/ purge oscillations; MacAyeal, 1993). The identification of the triggering mechanism requires more investigation that is beyond the scope of the current study. This H-like events has been recently observed in geological records in the North Pacific, but for other time intervals (Maier et al., 2018). The oscillating mechanism from our simulations could be applied to different timescales including the H-like millennial variability in the Pacific Ocean during MIS-3.
The transient response of the CIS is simulated under fixed orbital forcing. It is independent of the insolation control, and are mainly attributed to internal feedbacks. With this mechanism, the CIS is destabilized when it reached to a maximum state, and causes ice volume change by up to 5.2 m SLE. Also, the surging process resembles a sawtooth pattern, and has a period of around 20,000 to 25,000 years.

Conclusions
With the coupled climate-ice sheet model AWI-ESM-1.2, the climate during MIS-13 is simulated with three orbital configurations (495 ka, 506 ka, 517 ka). The simulations are initialized from the present-day ice sheet configuration, and the internal ice sheet-climate feedback is included.
The excess ice compared to present day from MIS-13 simulations are mainly from North American Cordillera, Arctic islands or Tibet, but the extent and ice volume of the ice sheets are different in response to different orbital configurations. In experiment MIS13-506, significant , ocean temperature (average from 0 to 500 m depth, same for ocean salinity), ocean salinity, ice thickness, ice sheet surface velocity, bedrock topography, bedrock uplift rate, ice sheet basal melt, and the total freshwater into the ocean over the northwest of the Cordilleran ice sheet (hatched area in Fig. 7b-d).
warming is found especially over the Northern Hemisphere continents with respect to both the preindustrial control run and the two MIS-13 simulations at 517 and 495 ka. Correspondingly, significant reduction of the Cordilleran ice sheet is simulated, since summer insolation is the main controlling factor for the Northern Hemisphere ice sheet changes. For experiments at 495 ka and 517 ka, no significant warming over the Northern Hemisphere continents in summer is simulated, and the CIS both reach full glacial conditions. Comparing the orbital configurations among the three MIS-13 experiments, the obliquity in MIS13-495 and MIS13-517 are at the highest and lowest points respectively, while boreal summers are both at aphelion. In contrast, boreal summer in MIS13-506 is at perihelion. We conclude that the precessional effect plays an important role on the evolution of the CIS. However, quantifying the individual effects of the precession or obliquity onto the climate-ice sheet system requires further idealized experiments with fully separations of the orbital parameters.
Our simulations with the highly sensitive Cordilleran Ice Sheet supports geologic evidence that the CIS has been glaciated many times, e.g. MIS-4 and earlier interglacials. This is probably due to the high elevations of these areas, that means the freezing point could be reached easily. However, no ice sheet is simulated over northeastern North America or Eurasia which may be related to a multi-stability of the ice sheets that requires model initialization from a glacial configuration. On the other hand, this also shows that the CIS is easier to build up than the other ice sheets. The initial expansion of the ice sheets requires a strong snow-albedo feedback. This can be easily reached for the CIS, while the other ice sheets might require stronger and more persistent cooling due to their warmer or lower elevations.
Periodic surges over the northwest of the CIS, which is independent of the insolation change is also simulated when the CIS reaches a full glacial states. This internal climate feedback points to the important role of the ocean and the glacial isostatic adjustment process. The oscillating mechanism from our simulations could be applied to different timescales including the H-like millennial variability.
For future studies, the multi-stability of the climate-ice sheet system should be further investigated with a complex Earth system model. A more advanced surface mass balance model should be applied, replacing the PDD scheme. In addition, only the dynamic Northern Hemisphere ice sheets are included in the current study, a dynamic Antarctica should also be included as a next step.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.