{'fraction_non_alphanumeric': 0.031377551020408165, 'fraction_numerical': 0.020884353741496598, 'mean_word_length': 4.641465988678883, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling the Antarctic Ice Sheet across the mid-Pleistocene transition – implications for Oldest Ice

The international endeavour to retrieve a continuous ice core, which spans the middle Pleistocene climate transition ca. 1.2–0.9 Myr ago, encompasses a multitude of field and model-based pre-site surveys. We expand on the current efforts to locate a suitable drilling site for the oldest Antarctic ice core by means of 3-D continental ice-sheet modelling. To this end, we present an ensemble of ice-sheet simulations spanning the last 2 Myr, employing transient boundary conditions derived from climate modelling and climate proxy records. We discuss the imprint of changing climate conditions, sea level and geothermal heat flux on the ice thickness, and basal conditions around previously identified sites with continuous records of old ice. Our modelling results show a range of configurational ice-sheet changes across the middle Pleistocene transition, suggesting a potential shift of the West Antarctic Ice Sheet to a marine-based configuration. Despite the middle Pleistocene climate reorganisation and associated ice-dynamic changes, we identify several regions conducive to conditions maintaining 1.5 Myr (million years) old ice, particularly around Dome Fuji, Dome C and Ridge B, which is in agreement with previous studies. This finding strengthens the notion that continuous records with such old ice do exist in previously identified regions, while we are also providing a dynamic continental ice-sheet context.

lines. Eight ice core locations are highlighted, which are used as tuning targets with respect to ice core thickness and analysed in Figure 9 (West Antarctica) and 10 (East Antarctica), respectively.

2.1 Ice Sheet Model
We employ the 3D thermomechanical Parallel Ice Sheet Model (PISM) (Bueler and Brown, 2009;Winkelmann et al., 2011) in the hybrid shallow-ice/shallow-shelfshallow-shelf/shallow-ice mode (SIA+SSA)(ssa+sia) with a subgrid grounding line parameterisation (Gladstone et al., 2010;Feldmann et al., 2014) to allow for reversible grounding line migration despite using 5 a relatively coarse resolution. Basal sliding is calculated with a pseudo-plastic sliding law (Schoof, 2010) in which the yield stress (τ c ) is determined by the pore water content and the strength of the sediment which is set by a linear piecewise function dependent on the ice-bedrock interface depth relative to sea level. The relevant parameter for this approach is introduced in PISM via the till-friction angle (see Winkelmann et al. (2011) eq. 12) which is scaled linearly between till min and till max , depending on the bedrock elevation (see Table 2). Through this heuristic parameterisation marine-based ice has a more slip- 10 pery base as compared to ice above sea level, allowing for faster flowing marine outlet glaciers. The parameter space used here yields a basal friction coefficient (e.g. underneath Thwaites glacier) on the lower end compared to Yu et al. (2018). Since the simulations presented here span a long time periodframe (of 2 Myr2 MA), we abstain from the derivation of basal friction by inversion (optimization problem) as we want to prevent over-tuning of present-day flow patterns. All simulations are carried out on a 16x16 km 2 grid and 81 vertical levels with refined resolution near the base (≈ 18 m at the ice-bedrock interface). 15 The grid resolution resolves the major ice streams while allowing for reasonable computation times (ca. 100 model years per processor hour on 144 cores, i.e. ca 5-7 days for 2 Myr depending on the supercomputer load), yet small outlet glaciers such as in the Antarctic Peninsula cannot be simulated adequately on this resolution.
The initial topography used for the simulations consists of a 200 kyr thermal spinup of the BEDMAP2 (Fretwell et al., 2013, 20 see Figure 1) data set (present day steady state simulation with fixed ice-sheetice sheet geometry), refined around LDC and Dome Fuji by the new radar derived topographies published in Young et al. (2017) and Karlsson et al. (2018). As basal heat flux is crucial for the existence of 1.5 Myr old ice (Fischer et al., 2013;Van Liefferinge and Pattyn, 2013) as well as for ice dynamics, especially in the interior of the ice sheet (Larour et al., 2012), we consider four different geothermal heat flux (GHF) data sets (Shapiro and Ritzwoller, 2004;Purucker, 2013;An et al., 2015;Martos et al., 2017) in our simulations to account for 25 uncertainties in GHF and to illustrate their impact on ice dynamics and potential Oldest Ice candidate sites.
Sea level plays an important role in the stability of marine ice sheets as it affects the position of the grounding line via the floatation criterion. We employ three different sea level reconstructions (see section 2.3) to account for different glaciation patterns in the northern hemisphere and different sea level highstands in interglacials. PISM does not account for self-gravitational 30 effects yet, which can have a stabilising effect on the ice sheet locally in interglacials (e.g. Konrad et al., 2014). Ice-shelf melt rates are calculated based on the parameterisation in Beckmann and Goosse (2003) (eq. 1), with a square dependency on the temperature difference between the pressure dependent freezing point and the ambient ocean temperature as used in e.g. Pollard and DeConto (2012), where M is the melt rate in m/s, m b is a scaling factor, ρ w and ρ i are ocean water and ice-shelf density, respectively, c pw is ocean water heat capacitiy, γ T heat transfer coefficient, L i latent heat, T f freezing point at depth of ice and T 3D ocean ambient ocean temperature. The ambient ocean temperature T 3D ocean is derived from simple extrapolation of the 3D ocean temperature into the 5 ice-shelf cavity. Recently, there have been developments towards more realistic representations of basal shelf melt in standalone continental ice-sheetice sheet models, incorporating sub-shelf ocean circulation (e.g. Reese et al., 2018;Lazeroms et al., 2018) which improve the representation of basal ice-shelf melt rates, but they have not been included in this study. To better match present day observed sub ice-shelf melt rates (Rignot et al., 2013;Depoorter et al., 2013), we had to multiplylinearly scale the computed present day melt rates in the Amundsen and Bellinghausen Sea by a factor of m b = 10, around the Antarctic at the calving front. If the ice thickness drops below this threshold, ice in the respective grid node is purged; 2. independently 15 of 1. we additionally employ Eigencalving (cE), which calculates a calving rate from the ice-shelf strain rates (Albrecht and Levermann, 2014). Both calving parameterisations are active simultaneously throughout the simulations.

Climate Forcing
To adequately capture continental ice-sheetice sheet dynamics on multi-milleniallong timescales (i.e. millennia and more), in principle, a coupled modelling approach is required which resolvesto resolve climate-ice-sheetice sheet interactions is required. 20 First efforts to tackle multi-millennial timescales via a fully coupled modelling approach are promising and currently being developed (e.g. Ganopolski and Brovkin, 2017). However, coupled climate-ice-sheetice sheet models which resolve ice-shelfocean interactions are mostly limited to applications on the centennial time scale due to computational limitations. To bridge this shortcoming, we construct a transient climate forcing over the last 2 Myr by expanding time-slice snapshots from the Earth system model (ESM) COSMOS (Lunt et al., 2013) with a climate indexthe glacial index method as applied in Sutter 25 et al. (2016). The climate snapshots are based on Pliocene (Stepanek and Lohmann, 2012), LIG (Pfeiffer and Lohmann, 2016), LGM and Pre-Industrial orbital, atmospheric and topographic conditions. For each climate snapshot, anomaly fields with respect to the pre-industrial control run are calculated and added to a mean Antarctic climatology , derived from the regional climate model RACMO (van Wessem et al., 2014), or the extrapolated World Ocean Atlas 2009 (Locarnini et al., 2010)  calculated by interpolating the anomaly fields with a climate indexglacial index approach, utilizing either of two climateglacial indices (see Figure 6 middle panel) derived from the Dome C deuterium record from Jouzel et al. (2007) which is expanded to the last 2 Myr by a transfer function (Michel et al., 2016) using the benthic oxygen isotope stackmarine sediment core from Lisiecki and Raymo (2005) (LR04) or the global surface temperature data set from Snyder (2016). To obtain an "Antarctic" surface temperature record from the far field benthic oxygen isotope stack, the LR04 isotope values are scaled via: LR04 T is the new surface temperature record, LR04 is the benthic isotope stack data (time corrected to match the AICC2012 time scale (Bazin et al., 2013;Veres et al., 2013)) and LR04 810 is the mean LR04 isotope data for the last 810 kyr, standard 5 deviations of the EDC and LR04 810 record are denoted by σ(EDC2007) and σ(LR04 810 ), respectively, mean Dome C surface temperature record is denoted by EDC2007. The forcing variables (surface temperature T s , ocean temperature T o ) can then be calculated at every grid point in time by: where indices i,j denotes the grid point, z denotes the depth of the ice ocean interface at grid point (i, j), T i,j s (t) is the surface temperature at grid point i, j at time t, T i,j spd is the surface temperature at present day (mean climatology from 1979-2016), ∆T s i,j ig , ∆T s i,j g and ∆T s i,j p , are the climate anomalies for the LIG, the LGM and the Pliocene, respectively. Ocean temperatures are derived in the same way (eq. 4). The linear scaling factors ω x (t) are derived from the climate indexglacial index (CI) which interpolates the climate forcing at any given time between the respective climate states. The scaling ω x (t) is 15 computed by : where the subscripts g, ig and p stand for Glacial, Interglacial and Pliocene respectively and CI pd refers to the present day  . Times colder than present are shaded in cyan and times warmer than present in red. d) same as as a) but for the climate indexglacial index derived from the Snyder global surface temperature record and scaling factors ωx in c). Times warmer than the Last Interglacial are shaded in dark red.
The lines in b) and c) vanish when scaling factors are zero.

7
The climate index is normalized with respect to the warmest climate period in the Dome C temperature record, therefore the LIG has index 1.0 in ensemble B1 and B2. The climate is linearly scaled between present day and LIG if the climate index surpasses the mean present day climate index CI pd , and between LIG and the Pliocene if the index is larger than 1.0. The major difference between the two glacial indices is the warmer overall climate state recorded in Snyder (2016) before the MPT (see (compared to in-situ measurements) very well, with biases in the high Antarctic plateaus of less than 5%.
We apply a temperature dependent scaling of precipitation (P ), using a scaling factor (percent precipitation change per degree Celsius) of α P of 3 and 5 %, respectively, motivated by central East Antarctic paleo precipitation changes (Frieler et al., 2015;Werner et al., 2018): 10 The precipitation is linearly dependent on the temperature anomaly with respect to present day temperature. Note that this scaling but underestimatesunderestimating the sensitivity of coastal mass balance to temperature changes.
In this standalone approach, the Antarctic Ice Sheet (AIS is responding to the external forcing, thus no feedbacks are acting between the ice sheet and the climate system. However, the climate indexglacial index approach implicitly incorporates the 15 integrated climate response to changes in orbital configuration and atmospheric CO 2 archived in the Dome C (Jouzel et al., 2007), the marine sediment core (Lisiecki and Raymo, 2005) or the global surface airair temperature (Snyder, 2016) record.
This allows to investigate the dynamical response of the AIS to a shifting climate regime across the MPT, with the caveat, that ice-sheetice sheet-climate interactions which are not included in the GCM time slice approach might have also played a significant role in the evolution of the AIS during the MPT transition.  (Locarnini et al., 2010) ocean temperatures) trying to fit the current sea level equivalent ice-sheetice sheet volume, geometry, ice flow, ice thickness at selected ice core locations (Figure 1), as well as the Antarctic sea level contribution during the last two glacial cycles.

30
The ensemble is built around two main branches of ensemble runs consisting of a set of boundary conditions (Table 1) and ice-sheetice sheet model parameterisations (Table 2). forcing combinations (three SL data sets and four geothermal heat flux fields). The parameter suite is derived from sensitivity studies in which the present day Antarctic Ice Sheet and its sea level contribution during the last two glacial cycles were the main tuning targets. In the first ensemble branch (B1) the climate indexglacial index is derived from an extrapolation of the EPICA Dome C temperature record (Jouzel et al., 2007) via correlation to the Lisiecki and Raymo (2005) time series to span the last two million years (Michel et al., 2016). In the second branch (B2) the climate indexglacial index is derived from the global air temperature record in Snyder (2016 (Stepanek and Lohmann, 2012) (see Figure 2).We explore two main parameter sets (P1 and P2) highlighted in Table 2. While we do take into account all sea level variations for ensemble B1, we only look at the sea level forcing derived from Lisiecki and Raymo (2005) (LR05) in ensemble B2. We also experimented with other parameter choices based on Table 2 (VP) but do not covering all individual forcing sets for these, thus they are not . The four panels illustrate the four GHF input data sets (Shapiro and Ritzwoller, 2004;Purucker, 2013;An et al., 2015;Martos et al., 2017) used in this study. As a reference the rounded GHF (in mW/m2) at selected ice core locations is provided. discussed in this study. In total we carried out 186 individual simulations. The ensemble members discussed in this manuscript consist of 8 experiments for each ensemble B1 and B2 with sea level forcing from LR05.

Results & Discussion
The main objective of this work is to assess the existence of 1.5 Myr old ice along the East Antarctic ice divide. We simulate 5 the evolution of the AISice dynamics throughout the last 2 Myr, focussing on ice volume changes specifically across the MPT (see Figure 5 and 6), as well as on ice-sheetice sheet configurations in glacials (focussing on marine isotope stage 2) and interglacials (with a focus on marine isotope stage 11 and 5, see Figure 7 and 8). We further investigate ice thickness changes at ice core locations in West and East Antarctica (see Figure 9 and 10). We conclude with a map of promising sites providing suitable conditions for an Oldest Ice ice core around Dome Fuji, Dome C and Ridge B, following the approach of Fischer et al.

Antarctic ice volume changes
We divide our discussion of the evolution of the AIS volume into three time frames: 1. pre-MPT (2 -1.  Figure 6, we present two clusters of the model ensemble from branch B1 (ice/marine sediment core climate indexglacial index) and the branch B2 (surface airglobal air temperature climate indexglacial index). Depicted are three simulations from both B1 and B2 with two identical model parameterisations (P1 and P2, respectively, see Table 2) using three different GHF data sets (Shapiro and Ritzwoller, 2004;Purucker, 2013;Martos et al., 2017).

15
Simulated ice volume changes before the MPT are characterised by a strong obliquity (≈ 41 kyr) cycle resembling the climate indexglacial index forcing which is formed by the integrated planetary response to orbital variations. The two clusters in the upper panel of Figure 6 show an a present day ice-sheetice sheet configuration similar to present day (B1-branch) and a strong Interglacial configuration (B2-branch) in which the West Antarctic Ice Sheet (WAIS) has collapsed. Note that most of the ensemble members do not allow a significantly increased glaciation as encountered during the last 800 kyr. The increase in 20 glacial AIS volume increase during the pre-MPT phase is limited to less than 2-4 m sea level equivalent ice volume throughout temperatures at present day and +2 • at peak interglacial). Additionally, we increase the sensitivity of the basal melt rate to ocean temperature changes in certain ocean basins (see method section). Pollard and DeConto (2009) prescribe basal melt rates directly, scaling them via the far field benthic isotope record (Lisiecki and Raymo, 2005) and austral summer insolation.  Table 2) and sea level forcing LR04. Geothermal heat flux depicted with blue line (Shapiro and Ritzwoller, 2004) upper light blue dashed line (Purucker, 2013) and lower light blue line (Martos et al., 2017). The red curves show B2 with parameter set P1 and using the climate indexglacial index from Snyder (2016). Red line (Shapiro and Ritzwoller, 2004), dashed maroon line (Martos et al., 2017), maroon line (Purucker, 2013)). The simulated ice volume from Pollard and DeConto (2009) Ultimately, this scaling leads to larger bulk ice-shelf melt rates and smaller melt rates close to the grounding line compared to the ones calculated in our approach. Overall, this leads to a more muted ice-loss response to warmer interglacial conditions during the pre-MPT in our ensemble and a generally lower variability in ice volume (sea level equivalent of ca. 4 − 8 m), while the growth and retreat phases are more or less synchronous to the variations in DeConto (2009) andde Boer et al. (2014). Interestingly, the differences in interglacial AIS volume between the three studies are largest in pre-MPT times, while 5 they are rather similar for the last four interglacials (see Figure 6). Evidently, the strongest interglacial AIS retreat is found in MIS7 for both DeConto (2009) andde Boer et al. (2014), while in our ensemble it is MIS11 and MIS5, with MIS5 producing a slightly stronger response. A coherent result ofOne main similarity between our study and the results from throughout the pre-MPT.

MPT Antarctic Ice-Sheet evolutiondynamics
Commonly, the onset of the MPT is put at 1.2 Myr BP and the MPT ends about 0.9 Myr BP culminating in the extended cold conditions between marine isotope stages 24 and 22 (ca. 940-880 kyr BP). Two pronounced interglacials (MIS31 and MIS25) 5 and several "colder-than-present" interglacials"lukewarm" events (between 1050 kyr and 950 kyr BP) separated by moderate glacial conditions throughout the MPT can be identified in the climate indexglacial index forcing ( Figure 6 panel b). The AIS response during the MPT is dominated by two proto-glacial states between 1.1 -1.0 Myr BP separated by interglacial MIS31, which can be interpreted as a first expression of a 100 kyr cycle. However, obliquity pacing still dominates ice ice volume changesdynamics at this stage. The second proto-glacial state after MIS31 is interupted by MIS25 followed by an extended glaciation/deglaciation before the MPT. This finding is robust across all tested sea level forcing data sets. We find no evidence of large changes in the EAIS margin during the MPT as hypothesized byas in Raymo et al. (2006), where a transition from a mostly land based EAIS to a marine EAIS similar to todays configuration is proposed. However, most simulations from B2, which include warm Pliocene climate conditions, show a major re-organisation of West Antarctica into a "present day" ice-sheetice sheet configuration at the end of the MPT (see Figure 9). This might represent a West Antarctic analoguependant 20 to the theory that the EAIS transitioned to a marine configuration during the MPT (Raymo et al., 2006), which does not require significant changes in the EAIS margin during the MPT. Such a configurational WAIS-shift would potentially implicate strong climate feedback mechanisms due to the formation of an ocean gateway between the Weddell, Ross and Amundsen Sea (Sutter et al., 2016) affecting climate dynamics across the MPT. This transition is not simulated in B1 and calls for a more crucial analysis outside the scope of this publication, e.g. incorporating a fully coupled ESM with a dynamical ice-sheetice sheet 25 component. Accordingly, the climate state in B1 does not allow a waxing and waning of the WAIS for pre-MPT interglacial conditions. We note, that other modelling studies either focussing on warmer Pliocene stages (DeConto and Pollard, 2016) or regional sensitivity studies (Mengel and Levermann, 2014) show large scale retreat of the grounding line into the Wilkes and Aurora subglacial basins, therefore a potential re-organization of the EAIS across the MPT cannot be excluded using all model experiments. The end of the MPT is marked by a pronounced glacial state at MIS 22 akin to the Last Glacial Maximum 30 reflecting a strong growth of the AIS at the end of the MPT. This result is robust across all ensemble members for both branch B1 and B2. It is interesting to note that the glaciations in the MPT interval become progressively stronger and reach a full late Quaternary glaciation state in MIS 22.

Post-MPT Antarctic Ice-Sheet evolutiondynamics
The simulated Quaternary AIS volume evolution can be roughly divided into two parts, the first spanning the window from 900 kyr BP to 420 kyr BP (MIS11) and the second from MIS11 to today. After MIS11, ice volume variability increases with smaller interglacial and bigger glacial ice sheets compared to the preceeding 500 kyr. This pattern mostly reflects the stronger interglacial atmospheric and ocean temperature forcing from MIS11 onwards. MIS11 is the first late Quaternary interglacial in 5 which the WAIS recedes to a land based ice-sheetice sheet configuration in B1 with the second major interglacial being MIS5e.
The ensemble mean sea level contribution in MIS5e amounts to ≈ 2.5 − 3 m with a full ensemble range between 1 and 4.5 m (see Figure 8). Glacial ice volume in the late Quaternary grows by ca. −8 to −10 m sea level equivalent ice volume (see Figure  8) with strongest glacials represented by MIS16 and 2. In general, the glacial extent of the AIS matches reconstructed LGM grounding margins rather well, with the notable exception of the Ammundsen and Bellinghausen sea sectors. In this region, 10 the ocean forcing seems to be too warm to allow for an advance of the ice margin to the continental shelf edge in the model.
LGM ice growth in the whole ensemble is strongly dependent on the SIA enhancement factor sia e , with values larger than 1.5 leading to an underestimation of ice thickness, albeit not necessarily ice extent. In the Ross Sea, ice thickness calving exerts Quaternary ice-sheetice sheet volume never drops below the present AIS volume. Note also, that the stability of the marine WAIS is crucially dependent on the choice of the GHF, as all simulations with the GHF field from Purucker (2013) exhibit a collapse of the WAIS in the LIG with a much smaller loss for both Martos et al. (2017) and Shapiro and Ritzwoller (2004). In our simulations, the applied sea level forcing plays a minor role in the stability of the WAIS in lnterglacials (not shown).

The role of geothermal heat flux
It is a well established fact, that the heat flux at the ice-bedrock interface can be a major driver of ice-sheetice sheet evolutiondynamics.
Albeit, the GHF for the AISAntarctic Ice Sheet is poorly constrained (Martos et al., 2017) and the few published continental 5 data sets available differ substantially (see Figure 4 and Martos et al. (2017)). We can analyse the impact of GHF in our model ensemble by gauging the fit of diagnostic variables (ice thickness, basal melt, basal temperature, ice volume) in comparison to observed data. Both data sets from Purucker (2013) and An et al. (2015) show relatively low heat flux along the East Antarctic ice divide and overall for the WAIS. In consequence, potential "Oldest Ice" candidate sites indicated by the model ensemble and using those two data sets are unrealistically large, due to the absence of basal melting (see Figure 11), specifically for the 10 Dome Fuji region. This is the case regardless of the model parameterisation and climate forcing or sea level input data. The choice of geothermal heat flux imprints on both East and West Antarctic ice thicknessdynamics., While mostly Modulating ice thickness (leading to thickness changes of up to 20%) in East Antarctica itbut also impactsimpacting the simulated ice divide along the transect of Dome A-Ridge B-Vostok-Dome C (see Figure 11). The impact of heat flux in West Antarctica can be drastic, as it acts as a major control on the marine ice-sheetice sheet instability. All simulations with the GHF field from Purucker 15 (2013) exhibit a collapse of the WAIS in the LIG with a much smaller percentage for both Martos et al. (2017) and Shapiro and Ritzwoller (2004). A thorough analysis of this result is beyond the scope of this manuscript, but it could be caused by larger glacial ice cover caused by the colder basal conditions in Purucker (2013). This would lead to an overdeepened bedrock and larger surface gradients along the coast at the onset of interglacials and therefore to favourable conditions for the marine ice sheet instability. Overall ice-sheetice sheet variability between ensemble members (for identical parameter settings) due to 20 different choices of GHF is consistently larger than due to the choice of different sea level forcing. This emphasises the strong role of GHF in modulating Antarctic Ice Sheet dynamics and overall evolution of ice volume.

Ice thickness variability in WAIS and EAIS
Ice thickness is an important parameter controlling the availability of a continuous 1.5 million year old ice core record at a given location. If the ice is too thick, its insulating properties lead to melting at the ice bedrock interface, if it is too thin, either 25 old ice is transported away too quickly or the layering might be too thin to decipher a meaningful transient climate signal.   thickness for the selected ensemble members in B1 and B2). A notable exception is the Talos Dome ice thickness, which is too thin in all discussed B1 and B2 ensemble members except for the runs simulated with the relatively cold (Purucker, 2013) or intermediate (Shapiro and Ritzwoller, 2004) GHF data set. The applied model resolution of 16 km is generally too coarse to accurately reconstruct smaller outlet glaciers and therefore might overestimate advection away from coastal ice domes. also has a substantial impact on overall ice thickness (up to 20% change between forcing sets, see Figure 10, comparing the difference between simulations with different GHF in the same ensemble).

Mapping potential Oldest Ice sites
Our findings indicate, that the existence of Oldest Ice is not only dependent on the choice of GHF but can also be influenced by the applied climate forcing or ice flow parameterisation modulating the regions illustrated in Figure 11. Despite regional

Conclusions
The search for the major drivers of the mid Pleistocene transition is ongoing, with ice sheets acting as a crucial component of climate system re-arrangement. Our model ensemble simulations indicate both rapid transitions and gradual growth of the Antarctic Ice Sheet during the MPT, featuring the buildup of a large glacial ice mass in West Antarctica driven by extended glacial conditions and muted interglacials between 1.2 and 0.9 Myr BP. These findings fit well to the notion of a significant expansionexpanse of the Antarctic Ice Sheet around 0.9 Myr BP by Elderfield et al. (2012). However, we do not find a major re-organisation of the EAIS grounding line across the MPT which is in contrast to the theory that the EAIS transitioned from a mostly land based ice sheet to a marine configuration in this interval (Raymo et al., 2006). However, such a re-organisation Despite this uncertainty, we identify promising candidate sites at Dome Fuji, Dome C and Ridge B which provide favourable conditions for the existence of old ice throughout the last 2 Myr. This study illustrates that uncertainties in climate forcing and boundary conditions have a large impact on paleoclimate ice-sheetice sheet simulations and therefore the assessment of Oldest Ice sites on long time scales. Accordingly, the successful retrieval of an ice core spanning the last 1.5 Myr would provide a 25 transient data benchmark and proxy horizons against which ice-sheetice sheet models can be calibrated.tuned and providing insights into a potential major re-arrangement of the East Antarctic coastal margin during the MPT.
Code and data availability. PISM is freely available via github, we use PISM v0.73 with modifications. Further details on the code modifications can be provided upon request. Diagnostic ice sheet model output data is available upon request.

Competing interests. O. Eisen is CEIC of The Cryosphere
Acknowledgements. We would like to thank Bas de Boer, Jorge Bernales and the editor Alexander Robinson for their very constructive comments in the review process. We thank the regional climate initiative REKLIM and the AWI Strategy Fund for funding of the project.  We thank Jorge Alvarez-Solas for his thorough review on our paper and for his constructive comments and remarks. We hope that we are able to address them in a satisfactory manner and discuss the modifications we introduced motivated by his remarks below. Please find all changes and the new or modified figures in the attached manuscript tracked in red. We also included the main changes below each comment in this document in red and our response to the reviewer's comments in blue.
1) Main goal of the study: I think this study nicely focuses on providing new information on where to look for the oldest Antarctic ice within the framework of continentally wide ice-sheet modeling. And it is also appreciable to subscribe such an approach within the context of transient simulations covering the MPT.
Nevertheless, several parts of the abstract and the introduction can be misleading because they are somehow suggesting that an important part of the paper will be devoted to analyze the effects of parametric uncertainty on the Antarctic ice dynamics during the MPT.
For example, the abstract reads: "We discuss the effects of changing climate conditions, sea level and geothermal heat flux boundary conditions on the mass balance and ice dynamics of the Antarctic Ice Sheet." This is not strictly true, particularly so with respect to ice dynamics (there is not a single panel or comment devoted to ice velocities or changes in the ice flow dynamics as a result of the changing climate). Changes in the grounding line position (also very superficially tackled), basal temperature and thickness evolution are related to changes in ice velocities, but this relationship is not addressed in the paper.
Similar misleading sentences to that of the abstract mentioned above can be found throughout the introduction (e.g. page 2, lines 27-29) and methods (page 7, lines 5-7).
I suggest two different ways of mitigating this issue: a) Lower the reader's expectations concerning ice dynamics. b) Expand the analysis on the effects that different parameters controlling the ice flow have on the conclusions of the study.
We concur with the reviewer, that the use of the term "dynamics" can potentially misguide the reader's expectations. While we do discuss e.g. grounding line migration in the Last Interglacial we do not dive into the intricacies of e.g. ice flow changes and individual rerouting of outlet glaciers. We rather put the focus on ice thickness changes around Oldest Ice candidate sites and analyse ice volume throughout the last 2 million years as an integrative measure of configurational changes of the Antarctic Ice Sheet. We appreciate the reviewer's suggestion to rephrase certain parts of the manuscript to lower the readers expectations with regard to the term "dynamics". While rewriting the manuscript with a detailed focus on ice dynamical changes during the last 2 Million years would be certainly worthwhile, we think that this would be better suited for a separate manuscript and would divert too much attention from the main focus of this study: to identify potential Oldest Ice sites for the IPICS endeavour.
We rephrased all potentially misleading sentences containing the term "dynamics" in the manuscript.

2) Reproducibility and description of the methods:
The manuscript contains a large amount of small inaccuracies, typos and erratum, particularly so in the methods section. These, taken individually, do not constitute a major problem, but taken together have a double negative effect: i) do not allow for the suitable reproducibility of this study, and ii) make the reading of the manuscript frustrating.

1) Ice-sheet model and ensembles
Page 4 This is indeed not self explanatory, we therefore add a short sentence clarifying, that the scaling factors multiply the computed melt rates.

P5, l8 -
To better match present day observed sub ice-shelf melt rates (Rignot et al., 2013;Depoorter et al., 2013), we had to multiply the computed present day melt rates in the Amundsen and Bellinghausen Sea by a factor of m_b=10, around the Antarctic Peninsula by 5, and underneath the Filchner Ice Shelf by a factor of 1.5. Shelf melt rates adjacent to Wilkes, Terre Adelie and George V Land in East Antarctica are also scaled by factor of 10 in a subset of the simulations.
Furthermore, the above quoted sentence seems incongruent with having a parameter between 1 and 10 as a part of the ensemble controlling the basal melt values (gamma_EAIS [1;10]). Caption of table 2 refers to George V and Wilkes Lands, while the above scaling factors refer to other basins as well. Thus, are the values of the shelf melt being changed as a part of the ensemble? If yes, both basins simultaneously? On the other hand, if these scaling factors reflect the uncertainty associated to the processes that determine shelf basal melt in time, why not exploring the factors of the rest of the basins as well? Please define and clarify.
We agree with the reviewer that this is unclear in the manuscript. Table 2 is now revised according to the suggestions of both reviewers and provides the two main parameter sets discussed in the manuscript, while additionally providing further parameter combinations used in the ensemble but not discussed in the results. We now address specifically which ensemble members are mainly analysed in the study (p10, l4-11) and changed table 2 (p10) accordingly. Table 2 now provides the two main parameter sets used for all sea level and geothermal heat flux combinations in ensemble B1 and for all geothermal heat flux combinations in B2. We further provide additional parameter combinations simulated in the ensemble but not discussed in the manuscript in the third line of table 2. We hope this clarifies the main scope of the ensemble.
P10, l4-11 The constituting forcing set for the ensemble consists of four different geothermal heat flux and three sea level data sets, i.e. twelve individual experimental settings. We further explore two main parameter sets (P1 and P2) highlighted in table 2. While we do take into account all sea level variations for ensemble B1 (48 individual experiments), we only look at the sea level forcing derived from Lisiecki and Raymo (2005) (LR05) in ensemble B2. We also experimented with other parameter choices based on table 2 (VP) but not covering all individual forcing sets, thus these are not discussed in this study. The ensemble members discussed in this manuscript consist of 8 experiments for each ensemble B1 and B2 with sea level forcing from LR05.
How is the shelf basal melting evolution in time achieved in the model? Is it dependent on the temperatures of the ocean given by equation 2? If yes, how is it done? Perhaps an anomaly method with respect to the scaled PD basal melt defined in page 4, line 29? Please define.
Basal melt rates are calculated according to the formulation in Beckmann and Goosse (2003) with a square dependency on the temperature difference between the pressure dependent freezing point and the ambient ocean temperature as used in e.g. Pollard and DeConto (2012). We added the melt rate equation on p5. The ocean temperature over time is derived as shown in equation 4.
P4, L32, calving parameterisations: Are these two parameterisations (threshold and eigen) exclusive to each other? If yes, then specify in the ensembles description that you explore the use of two different calving "laws", the threshold one with two values of its parameter and the eigen one with only one value. If no, then please remove the eigen parameter value to the table summarizing the explored values of the parameters for the ensemble.

P5, l12-17
These scaling factors are kept constant throughout the paleo simulations. Ice shelf calving and therefore the dynamic calving front is derived via two heuristic calving parameterisations: 1. thickness calving (cH) sets a minimum spatially uniform ice thickness (75 m or 150 m) at the calving front, if the ice thickness drops below this threshold, ice in the respective grid node is purged; 2. independently of 1. we additionally employ Eigencalving (cE), which calculates a calving rate from the ice shelf strain rates (Albrecht and Levermann, 2014). Both calving parameterizations are active simultaneously throughout the simulations. Table 2 shows the explored values of two parameters that are not described nor mentioned at all in the text: "sia" and "ssa". By looking at the values, I assume these refer to the enhancement factors of the SIA and the SSA parts of the simulated ice sheet. They need to be defined and properly named (sia and ssa are not parameters per se but approximations).
See changes applied to table 2 and caption table 2. Table 2 shows the parameters till_min and till_max. These are not defined nor described in the text. The reader can only speculate about the possibility that these parameters have something to do with the description of page 3, last line: "[...] the yield stress (tau_c) is determined by the pore water content and the strength of the sediment which is set by a linear piecewise function dependent on the ice-bedrock interface depth relative to sea level". If the reader keeps digging and tries to identify these parameters in the references of the model given here, will still not succeed because no mention to "till_max" and "till_min" can be found in Brown 2009, nor in Winkelmann et al, 2011. I had to go to Martien et al, 2011 (The Cryosphere) and assume that "till_min" and "till_max" correspond to the upper and lower numbers given by their equation 10. Is this correct?
Please define, clarify and cite accordingly.
The role of till_min and till_max was indeed not evident from the methods section, we added a short reference to the respective equation in Winkelmann et al 2011.
P4, l8-10 The relevant parameter for this approach is introduced in PISM via the till-friction angle (see Winkelmann et al. (2011) eq. 12) which is scaled linearly between till min and till max , depending on the bedrock elevation.

1) Glacial index description.
Having the glacial index in Figure 5 and knowing that you use 3 climate snapshots weighted in time by such an index, one can have an idea of how you are forcing the ice sheet model. However, the description of the method (including equations 1 to 5) is a bit odd. Perhaps it would be simply solved by providing the values of GI_pd and GI_max (as far as I saw the value of these parameters is not given in the manuscript).
We agree, that the description of the climate forcing should be expanded. To this end we introduce a plot which illustrates how the glacial index is implemented (p.7 new Figure 2). We also correct several typos in equations 5-7 (p6). Figure 2. climate index derived from Dome C deuterium record a) and corresponding scaling factors ω x in b). Times colder than present are shaded in cyan and times warmer than present in red. d) same as as a) but for climate index derived from the Snyder global surface temperature record and scaling factors ω x in c) . Times warmer than the last Interglacial are shaded in dark red.
Otherwise the reader could wonder: i) why is the value of the glacial index of ca. 0.7 during the Holocene? Does that mean you are taken a 30% of the LGM anomaly? (I guess not, and assume you put GI_pd to be approximately 0.7, so w_g in equation 3 goes to 0) The climate index (called glacial index before, see ii)) is normalized with respect to the warmest climate period in the Dome C temperature record, therefore the LIG has index 1.0 in ensemble B1 and B2. The climate is linearly scaled between present day and LIG if the climate index surpasses the mean Holocene climate index which is ca. 0.7 (see inclusion starting on p6, l20).
P6, l20- The respective values of the climate indices are CI lgm = 0.0, CI pd = 0.7022776, CI lig = 1.0,CI p = 1.130952. The climate index is normalized with respect to the warmest climate period in the Dome C temperature record, therefore the LIG has index 1.0 in ensemble B1 and B2. The climate is linearly scaled between present day and LIG if the climate index surpasses the mean present day climate index CI pd , and between LIG and the Pliocene if the index is larger than 1.0.
ii) Is "glacial index" the right term for a curve that goes to 0 during glacial times? Would not be more appropriate to define it as (1-GI) or simply call it "climate index"?
We agree that the term "glacial index" could lead to the impression that the index solely scales the climate between full glacial and present day conditions. We replace "glacial index" with "climate index" throughout the manuscript.
ii) What happens, for example when GI = 1.2 ? Do you take a linear interpolation of the Pliocene and the Last Interglacial climate fields? (I would assume so, because in equations 1 to 5 there is not any explicit differentiation of the time period, just a dependence on the values of the index, GI). If yes, can you justify it or elaborate?
Correct. Since the warmest climate period in the Snyder surface temperature record is not the LIG we interpolate between the LIG (climate index 1.0) and the warmest period in the Snyder climate index (ca. 1.2). Here, the climate is linearly scaled between peak LIG conditions (climate index 1.0) and peak climate conditions in the Pliocene (climate index ca 1.2). The justification for the manner we implement the climate index is, that we assume, that qualitatively the climate system fluctuated throughout the Quaternary between interglacials and glacials as it did during the last glacial interglacial cycle. We are aware that this is a simplification, but we only had access to paleoclimate time slices from the LIG, LGM and Pliocene.
Furthermore, equation 6 must contain a typo or be wrongly formulated. In its current form, the more you cool from PD the more you increase the precipitation. Is that correct?

2) Other specific/technical comments:
P5, L2 reads. "To adequately capture continental ice sheet dynamics on long timescales (i.e. millennia and more), in principle, a coupled modelling approach is required to resolve climate-ice sheet interactions." I think I understand what you mean here, but I also believe a sentence that says "To capture A, a given approach is required to resolve B" does not make sense.
Agreed and changed.
P5, l 19-20 To adequately capture continental ice-sheet dynamics on multi-millenial timescales, in principle, a coupled modelling approach which resolves climate-ice-sheet interactions is required.
P9. L24: "The two clusters in the upper panel of Figure 5 show a present day ice sheet configuration (B1-branch) and a strong interglacial configuration in which the WAIS is collapsed (B2-branch)". I think this sentence needs some rephrasing. By ".. show a present day and a strong interglacial" do you mean that their mean state is similar to the ones expected during present day and strong interglacial respectively?
Agreed and changed.
P12, l3-5 The two clusters in the upper panel of Figure 6 show an ice-sheet configuration similar to present day (B1-branch) and a strong interglacial configuration (B2-branch) in which the West Antarctic Ice Sheet (WAIS) has collapsed.
At the end of section 3.1 you state: " [...] all simulations with the GHF field from Purucker (2013) exhibit a collapse of the WAIS in the LIG with a much smaller percentage for both Martos et al. (2017) and Shapiro and Ritzwoller (2004)." Why is this? I see no obvious explanation since the West Antarctica GHF values from Purucker seem lower than those from Shapiro and way lower than Martos's.
Agreed. This sentence needs further explanation. However, we think that the role of geothermal heat flux on the stability of the WAIS and glaciers in general is a research field in its own right and therefore a detailed discussion would be beyond the scope of this manuscript. Regardless, we would like to keep the sentence in order to point out that GHF can have a decisive role for ice sheet systems which are close to a nonlinearity (i.e. MISI). We added a sentence, providing one potential reason for the outstanding effect of a very cold GHF underneath WAIS in Interglacials.
P16, l3-All simulations with the GHF field from Purucker (2013) exhibit a collapse of the WAIS in the LIG with a much smaller percentage of simulations for both Martos et al. (2017) and Shapiro and Ritzwoller (2004). A thorough analysis of this result is beyond the scope of this manuscript, but it could be caused by larger glacial ice cover caused by the cold Purucker (2013) basal conditions. This would lead to an over-deepened bedrock and larger surface gradients along the coast at the onset of interglacials and therefore favourable conditions for the marine ice sheet instability.  We also add a discussion of the LGM grounding line position and the spread of results regarding LGM sea level contribution in Figure 8 (formerly Figure 6).

P14, l29 -
In general, the glacial extent of the AIS matches reconstructed LGM grounding margins rather well, with the notable exception of the Ammundsen and Bellinghausen sea sectors. In this region, the ocean forcing seems to be too warm to allow for an advance of the ice margin to the continental shelf edge in the model. LGM ice growth in the whole ensemble is strongly dependent on the SIA enhancement factor sia e , with values larger than 1.5 leading to a underestimation of ice thickness, albeit not necessarily ice extent. In the Ross Sea, ice thickness calving exerts a strong influence on grounding line advance. A calving thickness of 75 generally leads to a good representation of LGM ice margin reconstructions (Bentley et al., 2014), while simulations with a thickness limit of 150m underestimate Ross Sea LGM grounding line advance. Furthermore, the parameterisation of ice shelf calving can play a preeminent role in interglacials, which underlines the dire need for a physical rather than heuristic representation of calving in ice-sheet models. We expanded Figure 7, to illustrate the effects of the different forcing approaches on the spatial configuration of the AIS during the LIG, LGM and PD.
"Middel" is wronlgy spelled. "Left" should be right. And invert the order of LGM and Last Interglacial.
Addressed (see new Figure 7). P15, L6 reads: "Mean ice thickness variability for Dome Fuji and Dome C during the late Quaternary is 165 and 195 m, respectively (105 and 140 during pre-MPT)." How is this "mean ice thickness variability" calculated? Is it the temporal standard deviation of the mean evolution of the ensemble, or has it something to do with the standard deviation of the ensemble itself?
Corrected. The increase is in mean ice thickness, not in variability. We added a sentence, pointing out the change in ice thickness variability (std) before and after the MPT.
P18, l8-10 Changes in the EAIS manifest in an increase in ice thickness by ca. 30% (mean ice thickness at ice core locations calculated for pre-MPT (1.8-1.2 Myr BP) and Quaternary (0.8-0.0 Myr BP)) time intervals alongside an increase in variability (standard deviation of ice thickness of individual ensemble members) of ca 50%. P15, L7 reads: "Overall, the simulated present day ice cover after 2 million years at the highlighted ice core locations is in good agreement (within ≈ 5%) with the BEDMAP2 (Fretwell et al., 2013) data set." This sentence does not seem precise enough. Does "present day ice cover" mean ice thickness? Here it is important to be precise, because it is not the same saying that the the error of the simulated region of Antarctica covered by ice falls within 5% with respect to BEDMAP2 than saying that the simulated thicknesses are within 5%. If the latter is what you meant, how did you calculate it?
We thank Bas de Boer for his thorough review on our paper and for his constructive comments and remarks. We hope that we were able to address them in a satisfactory manner and discuss the modifications we introduced motivated by the main remarks below. Please find all changes and the new or modified figures in the attached manuscript tracked in red. We also included the main changes below each comment in this document (response to reviewer comments in blue, changes to the manuscript in red).
1. Transfer function As pointed out on line 15, page 5, you have used a transfer function to describe 2 Myr of climate variability, combining the EDC ice core record and LR04 benthic oxygen isotope stack. Please use a couple of sentences to describe this function in the main text.
This information was indeed missing and we now include a short description including the transfer function used (p5 l27-32).

Description of temperature forcing (eqns. 1 and 2)
The description of the temperature forcing in equations (1) and (2) should be expanded a bit more and made clearer. For example, since you use the temperatures at the LIG, LGM and Pliocene as anomalies, use a _T in the equations and description of the variables in the text below. Also describe all three anomalies below the equation for both surface and ocean temperatures as (for example): "The ESM anomalies _Ts ig, _Ts g, and _Ts p, represent the LIG, LGM and Pliocene, respectively. Ocean temperatures (eqn. (2)) are defined in the same way." Secondly, when describing the calculation of the weighting factor w, using the glacial indices, it is unclear how you distinguish between the three states (wg, wig and wp). Also, it would be good to have a Figure 2 close by to refer to the GI values over time (instead of places in Figure 5 only). Also, what are the values of GIPD and GImax? Do these vary between the two GI records? Also, why did you not shifted the index to have GIPD = 1 or 0 for example? (since it is an index and can be shifted any way you want, as long as you coherently adapt equations 3-5). I do understand that since you have two GI records, based on Snyder and ice-core/_18O records, the differences between the two need to stay intact.
We agree, that the description of the climate forcing should be expanded. To this end we introduce a plot which illustrates how the glacial index is implemented (p.7 new Figure 2). We also correct a mistake in equations 4-6 (p6). Figure 2. climate index derived from Dome C deuterium record a) and corresponding scaling factors ω x in b). Times colder than present are shaded in cyan and times warmer than present in red. d) same as as a) but for climate index derived from the Snyder global air temperature record and scaling factors ω x in c) . Times warmer than the last Interglacial are shaded in dark red.

Comparison with PD09 and dB14
As mentioned on line 3, page 10, you shortly explain the type of forcing used in the two other studies. For a good comparison I do think you have to add a bit more explanation on how these two studies derived a transient climate forcing, and use this later on in your discussion as well. Both studies derive their long-term forcing in temperature and sea level from the LR04 benthic _18O record, although applied in different ways. Also, they both use a weighting index (eqn 6 in PD09, equation 7 in dB14) to prescribe the variation in sub-oceanic melting. This weighting index includes solar variations (through an anomaly at 80S) which mainly controls the waxing and waning of the WAIS (see dB14 supplement figure S4). The large peaks towards lower ice volume are caused by including this insolation anomaly, so take this into account in your comparison/discussion. For more description of the methods, please see Methods section at the end of both papers. Some more thoughts: It is a bit hard to distinguish individual simulations in Figure  5, but I guess your simulations of the B1-ensemble should be, in terms of timing, compare rather well to both PD09 and dB14 (depending on the sea-level forcing used). As you noted on line 6, page 10.
We agree, that the reader should be able to grasp major differences between the different studies approaches without doing a literature search first. We include a more detailed discussion of the approach in dB14 and P09 and a quick comparison of the main differences in the forcing approach with our study. We put the focus on basal shelf melt as we deem this to be the major difference affecting the simulated ice volume between the studies (p11 l 9-15, p12 l 1-2 & l 5-8.) One of the main differences in the approach here and the forcing applied in DeConto (2009) andde Boer et al. (2014), is the handling of basal melting underneath the ice shelves. This forcing component arguably exerts the strongest influence on grounding line migration of the AIS in interglacials. Our calculation of basal melt rates is very similar to de Boer et al. (2014), with lower differences between assumed peak interglacial and present day uniform ocean temperature. Peak interglacial temperature for ensemble B1 is approximately 2°C warmer than present day and 3°C warmer in B2 (3.7°C in de Boer et al (2014), -1.7°C at present day and +2°C at peak interglacial). Additionally, we increase the sensitivity of the basal melt rate to ocean temperature changes in certain ocean basins (see method section). Pollard and DeConto (2009) prescribe basal melt rates directly, scaling them via the far field benthic isotope record (Lisiecki and Raymo, 2005) and austral summer insolation. Ultimately this scaling leads to larger bulk ice shelf melt rates and smaller melt rates close to the grounding line compared to the ones calculated in our approach.

Structure of the ensemble
Based on the variables in Tables 1 and 2 you have a possible ensemble of 12 x 92 members. Please clarify in Section 2.3 how the two ensemble branches are built up. How many members does each have, and which specific experiments did you perform. All possible combinations? Or only a set? Also, in Table 2 the variable cE has only 1 value, so not a variable that you vary within the ensemble? If so, please remove form the table. What are the other settings of the model? Do you intend to include a parameter table? Or are these settings similar to previous simulations with PISM (please refer)? This is indeed confusingly presented in the manuscript. We now address specifically which ensemble members are mainly analysed in the study (p8, l 8-10 & p9, l 1-4) and changed table 2 (p9) accordingly. Table 2 now provides the two main parameter sets used for all sea level and geothermal heat flux combinations in ensemble B1 and for all geothermal heat flux combinations in B2. We further provide additional parameter combinations simulated in the ensemble but not discussed in the manuscript in the third line of table 2. We hope this clarifies the main scope of the ensemble.
The constituting forcing set for the ensemble consists of four different geothermal heat flux and three sea level data sets, i.e. twelve individual experimental settings. We explore two main parameter sets (P1 and P2) highlighted in table 2. While we do take into account all sea level variations for ensemble B1 (48 individual experiments), we only look at the sea level forcing derived from Lisiecki and Raymo (2005) (LR05) in ensemble B2. We also experimented with other parameter choices based on table 2 (VP) but do cover all individual forcing sets, thus these are not discussed in this study. The ensemble members discussed in this manuscript consist of 8 experiments for each ensemble B1 and B2 with sea level forcing from LR05. Figure 6 In general, I think your figures look really good. I do suggest you add panel labels to all figures (a,b,c, etc.) and use these in the captions. However, the left panel of figure  C4 6 is unclear, and I cannot distinguish between the dots or triangles at all. I suggest you create a separate figure from the left panel and the two other panels (the maps). In the new figure 6, make the symbols much bigger and perhaps put the time frames next to each other (LGM, LIG and you could also add the simulated PreInd or PD ice volumes in the middle perhaps?). Please also mention the total number of simulations shown in the plot. Nice to add the values, that should be included again. In the two panels with the map, please indicate which specific simulations you have used here. Is this a minimum/maximum within the ensemble, or a (sort of) reference simulation that represents the middle/median of the ensemble? Are the lines of LGM, PD and LIG shown in both panels?  LGM (c/d) for the full ensemble forced with climate index B1 and climate index B2. The ensemble members focused on in this paper are highlighted with colours (P1 red colors, P2 blue/green colors). Horizontal black line depicts the full ensemble mean and the dashed line the standard deviation.

Revise
We expanded Figure 7, to illustrate the effects of the different forcing approaches on the spatial configuration of the AIS during the LIG, LGM and PD.  7. Reference list I have added some minor comments in the reference list, seems to have some issues with the copy from Bibtex (in case you have used latex). Please carefully check your list of references on errors in mainly page numbers and hyperlinks via doi. This is now corrected.
We thank Bas de Boer for his positive review of our manuscript and hope that the revised manuscript satisfactorily addresses all critical points.