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.

Abstract. 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.

Introduction
The middle Pleistocene transition (MPT) is characterised by a shift from obliquity-driven climate cycles ( ∼ 41 000 years, 41 kyr) to the signature sawtooth ∼ 100 kyr cycles typical for the late Pleistocene. The drivers behind the MPT are still under debate and touch on the basic understanding of the climate system. The absence of any clear disruptive change during the MPT in orbital forcing makes the transition especially puzzling. Several theories have been put forth, striving to explain the enigmatic MPT (Raymo and Huybers, 2008). They include a shift in subglacial conditions underneath the Laurentide Ice Sheet (regolith hypothesis by Clark and Pollard, 1998), the inception of a large North American Ice Sheet (Bintanja and van de Wal, 2008) or marine East Antarctic Ice Sheet by Raymo et al. (2006), ice bedrock climate feedbacks , the buildup of large ice sheets between MIS24 and 22 identified by Elderfield et al. (2012), or the combination of changes in ice-sheet dynamics and the carbon cycle (Chalk et al., 2017). Ultimately, it seems likely that an interplay of the various proposed processes culminated in the MPT. To illuminate the potential role of these different processes and thus to solve one of the grand challenges of climate research, the recovery of a continuous ice core spanning at least beyond the MPT (in the following termed as "Oldest Ice") is crucial. An expansion of the currently longest ice core record from the EPICA Dome C project (Jouzel et al., 2007) to and beyond the MPT would provide the necessary atmospheric boundary conditions (i.e. atmospheric greenhouse gases and surface temperature) to revisit the current theories . It would provide a direct record of global atmospheric CO 2 and CH 4 concentrations and local climate during the MPT and beyond. A transient record of CO 2 concentrations would provide a key piece of the puzzle in answering the question of whether greenhouse gases were the main culprit behind the MPT, while proxies of climate conditions in Antarctica would illuminate the evolution of the Antarctic Ice Sheet leading to the MPT. However, retrieval of such an ice core is a challenging task, as a multitude of prerequisites must be met to recover an undisturbed ice core reaching more than a million years into the past Van Liefferinge and Pattyn, 2013;Parrenin et al., 2017). The European ice core community has identified a most promising target for an Oldest Ice drill site to be close to a secondary dome in the vicinity of Dome C, usually referred to as "Little Dome C" (LDC) . Also, other potential locations are targeted around Dome Fuji. The selection of sites is motivated by a series of recent studies based on radar observations of the internal ice-sheet stratigraphy and underlying bedrock topography Karlsson et al., 2018), local paleo-climate conditions , and 1-D and 3-D ice flow modelling Parrenin et al., 2017;Passalacqua et al., 2017). These studies provide a detailed view on the regional properties such as ice flow, thermal conditions and bedrock topography, enabling a localised assessment of promising drill sites. The only component missing so far in the analysis is the transient, continental paleo-ice-sheet dynamics perspective, which allows for the assessment of large-scale reorganisations of ice-sheet flow and geometry during glacial and interglacial cycles, their impact on divide migration, ice thickness changes along the East Antarctic ice divide and basal melt. There are many studies focussing on the transient evolution of the Antarctic Ice Sheet (AIS) during specific climate episodes in the past such as the Last Interglacial (LIG) (Sutter et al., 2016;DeConto and Pollard, 2016) or the Last Glacial Maximum (LGM) (e.g. Golledge et al., 2014). However, so far only a few studies cover the waxing and waning of the AIS during the MPT (Pollard and DeConto, 2009;de Boer et al., 2014) or late Quaternary (Tigchelaar et al., 2018) in transient model simulations with an evolving climate forcing. We build on these efforts by carrying out ensemble simulations of the Antarctic Ice Sheet across the last 2 Myr to investigate the MPT and the effect of glacial-interglacial variations in ice thickness and basal melting on potential Oldest Ice drill sites.

Ice-sheet model
We employ the 3-D thermomechanical Parallel Ice Sheet Model (PISM) (Bueler and Brown, 2009;Winkelmann et al., 2011) in the hybrid shallow-ice/shallow-shelf mode (SIA+SSA) with a subgrid grounding line parameterisation (Gladstone et al., 2010;Feldmann et al., 2014) to allow for reversible grounding line migration despite using a relatively coarse resolution. Basal sliding is calculated with a pseudoplastic 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 slippery 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 period (of 2 Myr), we abstain from the derivation of basal friction by inversion (optimisation problem) as we want to prevent over-tuning of present-day flow patterns. All simulations are carried out on a 16 × 16 km 2 grid and 81 vertical levels with refined resolution near the base (≈ 18 m at the ice-bedrock interface). 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 d 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, see Fig. 1) data set (present-day steady-state simulation with fixed ice-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 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 geothermalheat-flux (GHF) data sets (Shapiro and Ritzwoller, 2004;Purucker, 2013;An et al., 2015;Martos et al., 2017) in our simulations to account for 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 Sect. 2.3) to account for different glacia-J. Sutter et al.: Modelling the AIS across the mid-Pleistocene transition 2025 Figure 1. Antarctic bedrock topography overlain by surface contours (grey lines). The present-day (PD) grounding line (grl) from BEDMAP2 (Fretwell et al., 2013) depicted by the dashed black line. The Last Glacial Maximum (LGM) grounding line reconstruction from Bentley et al. (2014) (thick black lines) is compared to simulated grounding line retreat in one of the ensemble members for the Last Interglacial (LIG, red line). Regions previously identified as potentially viable sites for Oldest Ice (Van Liefferinge and Pattyn, 2013) are outlined by thick black lines. Eight ice core locations are highlighted, which are used as tuning targets with respect to ice core thickness and analysed in Figs. 9 (West Antarctica) and 10 (East Antarctica), respectively. tion patterns in the Northern Hemisphere and different sea level highstands in interglacials. PISM does not account for self-gravitational 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 Pollard and DeConto (2012), where M is the melt rate in metres per second (m s −1 ); m b is a scaling factor; ρ w and ρ i are ocean water and ice-shelf density, respectively; c pw is ocean water heat capacity; γ T is the heat transfer coefficient; L i is latent heat; T f is the freezing point at depth of ice; and T 3-D ocean is the ambient ocean temperature. The ambient ocean temperature T 3-D ocean is derived from simple extrapolation of the 3-D ocean temperature into the ice-shelf cavity. Recently, there have been developments towards more realistic representations of basal shelf melt in stand-alone continental ice-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 Depoorter et al., 2013), we had to multiply the computed present-day melt rates in the Amundsen and Bellingshausen 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 Adélie and George V Land in East Antarctica are also multiplied by factor of 10 in a www.the-cryosphere.net/13/2023/2019/ The Cryosphere, 13, 2023-2041, 2019 subset of the simulations. 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 as follows.
(1) Thickness calving (cH) sets a minimum spatially uniform ice thickness (75 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 . Both calving parameterisations are active simultaneously throughout the simulations.

Climate forcing
To adequately capture continental ice-sheet dynamics on multi-millennial timescales, in principle, a coupled modelling approach which resolves climate-ice-sheet interactions is required. 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-icesheet models which resolve ice-shelf-ocean interactions are mostly limited to applications on the centennial timescale 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 index method as applied in Sutter et al. (2016). The climate snapshots are based on Pliocene (Stepanek and Lohmann, 2012), LIG (Pfeiffer and Lohmann, 2016), LGM and preindustrial 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) to provide the climate forcing for the individual climate epoch. The intermediate climate states between the snapshots are calculated by interpolating the anomaly fields with a climate index approach, utilising either of two climate indices 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 stack 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 timescale; Bazin et al., 2013;Veres et al., 2013) and LR04 810 is the mean LR04 isotope data for the last 810 kyr; standard 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); and T sig i,j , T sg i,j and T sp i,j 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 ω x (t) derived from the climate index (CI) interpolates the climate forcing at any given time between the respective climate states. The scaling ω x (t) is computed by where the subscripts g, ig and p stand for glacial, interglacial and Pliocene respectively and CI pd refers to the present-day climate index. The respective values of the climate indices are CI lgm = 0.0, CI pd = 0.70, CI lig = 1.0 and CI p = 1.13. The climate index is normalised 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 Fig. 2b). The present-day forcing derived from van Wessem et al. (2014) matches the presentday climatology in Antarctica (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): The precipitation is linearly dependent on the temperature anomaly with respect to present-day temperature. Note that this scaling underestimates the sensitivity of coastal mass balance to temperature changes. In this stand-alone approach, the AIS is responding to the external forcing; thus no feedbacks are acting between the ice sheet and the climate system. However, the climate index approach implicitly incorporates the 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 air temperature (Snyder, 2016) record. This allows one to investigate the dynamical response of the AIS to a shifting climate regime across the MPT, with the caveat that ice-sheet-climate interactions which are not included in the  (2016) forming ensemble branch B1 and B2. Prescribed input data consist of sea level (SL) data and geothermal-heat-flux (GHF) data sets. Both ensemble B1 and B2 are constructed with 12 different 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.
GCM time-slice approach might have also played a significant role in the evolution of the AIS during the MPT.

Model ensemble approach
We choose a model ensemble approach to address the multitude of uncertainties regarding the paleo-climate state during the last 2 million years, the applied boundary conditions and the physics of ice flow. The aim of the ensemble design ( Fig. 3) is to investigate the impact of different climate forcings, the response of the AIS to different geothermal-heatflux signatures (Fig. 4) and the impact of sea level on the transient configuration of marine ice sheets. Ultimately, different manifestations of ice-sheet flow and climate response are investigated via a set of ice-sheet model parameterisations. The model parameters are preselected in equilibrium simulations under present-day forcing (1979-2011 climatology from van Wessem et al., 2014, andWorld Ocean Atlas Locarnini et al., 2010, ocean temperatures) trying to fit the current sea level equivalent ice-sheet volume, geometry, ice flow, ice thickness at selected ice core locations (Fig. 1) and the Antarctic sea level contribution during the last two glacial cycles. The ensemble is built around two main branches of ensemble runs consisting of a set of boundary conditions (Table 1) and ice-sheet model parameterisations (Table 2).
In the first ensemble branch (B1) the climate 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 2 million years (Michel et al., 2016). In the second branch (B2) the climate index is derived from the global temperature record in Sny- Table 1. Available choices of selected forcing fields for the model ensemble. B1 and B2 stand for the two glacial indices derived from Michel et al. (2016) and Snyder (2016); SL data from Lisiecki and Raymo (2005) . 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 cover all individual forcing sets for these; thus they are not discussed in this study. In total we carried out 186 individual simulations. The ensemble members discussed in this paper consist of eight experiments for each ensemble B1 and B2 with sea level forcing from LR05.

Results and 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 the evolution of the AIS throughout the last 2 Myr, focussing on ice volume changes specifically across the MPT (see Figs. 5 and 6), as well as on ice-sheet configurations in glacials (focussing on marine isotope stage 2) and interglacials (with a focus on marine isotope stage 11 and 5; see Figs. 7 and 8). We further investigate ice thickness changes at ice core locations in West and East Antarctica (see Figs. 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. (2013) and Van Liefferinge and Pattyn (2013) (Fig. 11).

Antarctic ice volume changes
We divide our discussion of the evolution of the AIS volume into three time frames: pre-MPT (2-1.2 Myr BP), MPT (1.2-0.9 Myr BP) and post-MPT (0.8-0 Myr BP). To put our results into perspective, we compare them to two published Figure 4. 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 m 2 ) at selected ice core locations is provided.  Tigchelaar et al. (2018) and with respect to different choices of GHF and climate index. In Fig. 6, we present two clusters of the model ensemble from branch B1 (ice/marine sediment core climate index) and branch B2 (surface air temperature climate index). Depicted are three simulations from both B1 and B2 with two 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).

Pre-MPT Antarctic ice-sheet evolution
Simulated ice volume changes before the MPT are characterised by a strong obliquity (≈ 41 kyr) cycle resembling the climate index forcing which is formed by the integrated planetary response to orbital variations. The two clusters in the upper panel of Fig. 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. Note that most of the ensemble members do not allow a significantly increased glaciation as encountered during the last 800 kyr. The increase in glacial AIS volume during the pre-MPT phase is limited to less than 2-4 m sea level equivalent ice volume throughout the model ensemble. Variability in the B2 branch (red) is higher than in B1 (blue), due to the waxing and waning of the marine WAIS and stronger glacial-interglacial surface mass balance variability. Ice volume in the B1 branch is predominantly driven by surface mass balance and sea level. The comparison to Pollard and DeConto (2009) 2014), with smaller differences between assumed peak interglacial and present-day uniform ocean temperature. Peak interglacial ocean temperatures for ensemble B1 are approximately 2 • C (Last Interglacial) warmer than present day and 3 • C warmer (Pliocene) in B2 (3.7 • C in de Boer et al., 2014, with −1.7 • circum-Antarctic ocean 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 De-Conto (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. Overall, this leads to a more muted iceloss 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 they are rather similar for the last four interglacials (see Fig. 6). Evidently, the strongest interglacial AIS retreat is found in MIS7 for both Pollard and DeConto (2009)

MPT Antarctic Ice Sheet evolution
Commonly, the onset of the MPT is put at 1.2 Myr BP and the MPT ends at about 0.9 Myr BP, culminating in extended cold conditions between marine isotope stages 24 and 22 (ca. 940-880 kyr BP). Two pronounced interglacials (MIS31 and MIS25) and several colder-than-present interglacials (between 1050 kyr and 950 kyr BP) separated by moderate glacial conditions throughout the MPT can be identified in the climate index forcing (Fig. 6b). The AIS response during the MPT is dominated by two proto-glacial states between 1.1 and 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 volume changes at this stage. The second proto-glacial state after MIS31 is interrupted by MIS25 followed by an extended cold period, allowing for the formation of marine ice sheets in the Weddell and Ross Sea similar to what is observed during late Quaternary glacials, marking the end of the MPT and the onset of unperturbed 80-120 kyr cycles in AIS volume. Our icesheet model results are in line with the notion of a 900 kyr event by Elderfield et al. (2012), which is centred around MIS25-22, manifesting itself in a qualitative difference in the formation of glacials: a long buildup phase ended by a sharp decline of the ice volume into interglacials (late Quaternary sawtooth pattern) during the last 800 kyr and a more symmetric 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 hypothesised by Raymo et al. (2006), where a transition from a mostly land based EAIS to a marine EAIS similar to today's configuration is proposed. However, most simulations from B2, which include warm Pliocene climate conditions, show a major reorganisation of West Antarctica into a present-day ice-sheet configuration at the end of the MPT (see Fig. 9). This might represent a West Antarctic analogue 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 po-   Table 2) and sea level forcing LR04. Geothermal heat flux is depicted with the 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 index from Snyder (2016). The red line refers to (Shapiro and Ritzwoller, 2004), the dashed maroon line refers to (Martos et al., 2017) and the maroon line refers to (Purucker, 2013). The simulated ice volume from DeConto (2009), de Boer et al. (2014) and Tigchelaar et al. (2018) is shown for comparison (black, grey and black dashed line). (b) shows the climate index used in B1(dark grey) and B2 (light grey) with the horizontal grey dashed line depicting the average Holocene index. (c) shows the sea level reconstructions used in the model ensembles (Lisiecki and Raymo, 2005, LR04;de Boer et al., 2014, dB14;and Rohling et al., 2014, R14).
tentially 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-sheet 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 largescale retreat of the grounding line into the Wilkes and Aurora subglacial basins; therefore a potential reorganisation 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 MIS22 akin to the Last Glacial Maximum 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 MIS22.

Post-MPT Antarctic Ice Sheet evolution
The simulated Quaternary AIS volume evolution can be roughly divided into two parts, the first spanning the window from 900 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 preceding 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 which the WAIS recedes to a landbased ice-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 Fig. 8).
Glacial ice volume in the late Quaternary grows by ca. −8 to −10 m sea level equivalent ice volume (see Fig. 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 Amundsen and Bellingshausen 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 an 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 m generally leads to a good representation of LGM ice margin reconstructions (Bentley et al., 2014) (Fig. 7b), while simulations with a thickness limit of 150 m 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. The different forcing approaches of our study and DeConto (2009), de Boer et al. (2014) and Tigchelaar et al. (2018) are apparent, for example, in the largest ice-sheet retreat occurring in our ensemble at MIS5e while occurring during MIS7 (ca. 210 kyr BP) for both DeConto (2009) andde Boer et al. (2014), whereas in Tigchelaar et al. (2018) Quaternary ice-sheet volume never drops below the present AIS volume.

The role of geothermal heat flux
It is a well established fact that the heat flux at the icebedrock interface can be a major driver of ice-sheet evolution. However, the GHF for the AIS is poorly constrained (Martos et al., 2017) and the few published continental data sets available differ substantially (see Fig. 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 Fig. 11), specifically for the 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 thickness. Modulating ice thickness (leading to thickness changes of up to 20 %) in East Antarctica also impacts the simulated ice divide along the transect of Dome A-Ridge B-Vostok-Dome C (see Fig. 11). The impact of heat flux in West Antarctica can be drastic, as it acts as a major control on the marine ice-sheet instability. 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). A thorough analysis of this result is beyond the scope of this paper, 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 icesheet variability between ensemble members (for identical parameter settings) due to 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, the layering might be too thin to decipher a meaningful transient climate signal. Figures 9 and 10 show the ice thickness change at the four East and West Antarctic ice core locations which are depicted in Fig. 1.
In the simulations with a collapsed late Pliocene WAIS (ensemble B2), the MPT leads to the advance of the WAIS into a present-day configuration going along with a closure of the open ocean connection between the Weddell, Amundsen/Bellingshausen and Ross Sea. However, in those simulations with no early Pleistocene WAIS, the WAIS remains relatively small throughout the Quaternary, which is likely to be an indication that the climate forcing is too warm. The ensemble members with a stable interglacial pre-MPT WAIS transition to a higher variability during the MPT with no major reorganisation of ice flow and grounding line dynamics.
Changes in the EAIS manifest in an increase in ice thickness across the MPT 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 in glacial-interglacial ice volume (standard deviation of ice thickness of individual ensemble mem-bers) of ca 50 %. The shift to larger glacial-interglacial ice thickness changes in individual ensemble members is evident at all ice core locations and mostly due to the more pronounced climate cycles following the MPT. The temporal evolution of ice thickness changes for the different East Antarctic ice core locations follows a similar pattern with muted variability during the pre-MPT and a gradual increase completed by ca. 900 kyr BP. Note that we find different points in time regarding ice thickness maxima for the central dome positions Dome C and Dome Fuji compared to positions away from the central domes such as around the EDML and Talos Dome ice core sites. While the former generally show ice-sheet maxima during mid-interglacials, the latter show maxima generally at the onset of interglacials and declining ice thicknesses during the interglacial. This and the higher glacial-interglacial variability for the two locations (EDML and Talos) can be explained by the impact of grounding line migration and larger glacial-interglacial surface mass balance differences. Mean ice thickness variability for Dome Fuji and Dome C during the late Quaternary is 165 and 195 m,. Overall, the simulated present-day ice thickness after 2 Myr at the highlighted ice core locations in East Antarctica is in good agreement with the ice thickness derived from the BEDMAP2 (Fretwell et al., 2013) data set (within ≈ 5 % discrepancy in ice 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.

Mapping potential Oldest Ice sites
We apply the conditions for the existence of 1.5 Myr old ice derived in Fischer et al. (2013) (ice thickness larger than 2000 m, basal melting zero and surface ice velocity slower than 1 m a −1 ) to our simulations in order to investigate the impact of the transient paleo-climate forcing, GHF and different model parameterisations on the sites Dome Fuji, Dome C, Dome A and Ridge B (Fig. 11). Overall, we identify similar Oldest Ice regions as in Van Liefferinge and Pattyn (2013), Van Liefferinge et al. (2018), Parrenin et al. (2017) and Passalacqua et al. (2017). This gives us confidence in the robustness of our model results and shows that the transient forcing and continental setup used here does not change the general conclusions in Van Liefferinge et al. (2018). The regions with major overlaps to Van Liefferinge et al. (2018) are thus promising sites from a paleo-ice-climate evolution viewpoint as well as from the detailed dissection of the present-day conditions. Overall, applying the GHF from the data set by Shapiro and Ritzwoller (2004) Parrenin et al. (2017) and Passalacqua et al. (2017). However, the geothermal-heatflux data set from Martos et al. (2017) matches the observed/derived heat flux at Dome C to a better degree. In the simulations with the Martos et al. (2017) data set, the present-day ice divide around Dome C is shifted strongly (ca. 160 km) in the direction of the Ross Sea/Belgica Subglacial Highlands compared to the other three data sets due to the low geothermal heat flux south-west of the Dome C region contrasting with the high heat flux around Dome C (see Fig. 11). Simulations using the Martos et al. (2017) data set Figure 9. Ice thickness evolution during the last 2 Myr as simulated in our model ensemble for the four West Antarctic ice core locations (see Fig. 1). Blue lines are from ensemble B1 (as in Fig. 6) and red lines from B2. Sampling rate is 1 kyr. For comparison, the ice thickness evolution simulated in Pollard and DeConto (2009) is plotted in black (sampling rate 10 kyr). The observed present-day ice thickness, derived from BEDMAP2 (Fretwell et al., 2013), is shown by the horizontal black dashed line, with 5 % variation illustrated by grey dashed lines. The ice thickness simulated in Pollard and DeConto (2009) is based on a 5 Myr transient simulation using BEDMAP1 (Lythe et al., 2001) as the initial ice-sheet configuration.
show almost no viable conditions for the existence of Oldest Ice in East Antarctica (with the exception of Oldest Ice patches around Dome A, Ridge B and in the Belgica Subglacial Highlands west of Dome C), as basal temperatures are relatively high due to the large heat flux at the base of the ice leading to sustained basal melting. This is the case despite the relatively low ice thickness (and therefore insulation) simulated with the Martos et al. (2017) data set. Simulations with the Purucker (2013) and An et al. (2015) GHF reconstruction overestimate the Oldest Ice area substantially (showing viable conditions also for the Dome Fuji and Dome C ice core location where subglacial melting is observed). Interestingly, however, the Purucker (2013) data sets yields the best agreement between the simulated and observed presentday ice divide for Dome Fuji as well as Dome C. The effect of the geothermal-heat-flux forcing is not limited to the basal temperature of the ice but also has a substantial impact on overall ice thickness (up to 20 % change between forcing sets; see Fig. 10, comparing the difference between simulations with different GHF in the same ensemble).
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 Fig. 11. Despite regional differences, all three sites (Dome Fuji, Dome C, Vostok/Ridge B) show suitable conditions for an Oldest Ice core throughout the last 1.5 Myr. This accounts for ensemble members using either the GHF from Shapiro and Ritzwoller (2004) or Purucker (2013) and An et al. (2015) (see Fig. 11). However the latter two show unreasonably large Oldest Ice patches. While previously identified regions of Oldest Ice (e.g. in the Dome C area; Parrenin et al., 2017) are also identified in our model results, the regions of Oldest Ice at Dome Fuji are somewhat shifted in comparison to the findings in Van Liefferinge and Pattyn (2013) and Van Liefferinge et al. (2018). This might be due to the different employed input data or ice dynamics simulated, or af- Figure 10. Same as Fig. 9 but for the four East Antarctic ice cores.
fected, by differences in model resolution or the thermal state at the base of the ice. However, the generally robust agreement between high-resolution studies Parrenin et al., 2017;Passalacqua et al., 2017) and our coarse-resolution paleo-dynamics approach strengthens the notion of viable conditions for Oldest Ice both at Dome Fuji and Dome C as well as Ridge B.

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 rearrangement. 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 expansion of the Antarctic Ice Sheet around 0.9 Myr BP by Elderfield et al. (2012). However, we do not find a major reorganisation 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 ma-rine configuration in this interval (Raymo et al., 2006). However, such a reorganisation cannot be excluded at this point, as several other modelling studies show potential phases of major EAIS grounding line retreat in the Pliocene or under strong interglacial conditions (DeConto and Pollard, 2016;Mengel and Levermann, 2014), and proxy reconstructions indicate large-scale EAIS grounding line retreat in the late Pleistocene (Wilson et al., 2018). While we do not simulate a transition to a marine-based EAIS in the Wilkes and Aurora basins, such a process would certainly imprint on ice flow around LDC due to proximity alone. We do find a clear transition of the WAIS configuration around 0.9 Myr in model runs where warm boundary conditions led to a collapse of the WAIS in the late Pliocene but allowed for glaciation during the colder interglacials in the MPT. We argue that such a transition between an ice-free and a present-day WAIS configuration around 900 kyr BP would have a potentially stronger influence on the global climate system (compared to an advance of the EAIS grounded ice margin) via the closing of the gateway between the Ross, Amundsen and Weddell Sea. Additionally, climate conditions favourable for a retreated configuration of the EAIS ice margin (Raymo et al., 2006) would imply a collapsed WAIS according to ice-sheet modelling studies (Golledge et al., 2015;DeConto and Pollard, 2016;Sutter et al., 2016). Our study confirms a strong contribution of the Antarctic Ice Sheet to the LIG sea level highstand, with a mean contribution of 2.5-3 m (depending on the climate index record used and the applied boundary conditions) and a maximum contribution of ca. 4.5 m, which is in line with previous studies by Golledge et al. (2015), Sutter et al. (2016) and DeConto and Pollard (2016). This corroborates the major impact of the Antarctic Ice Sheet on LIG global sea level (Dutton et al., 2015). The spatial pattern of GHF can be decisive in modelling dynamics of the WAIS. 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 paleo-climate ice-sheet simulations and therefore the assessment of Oldest Ice sites. Accordingly, the successful retrieval of an ice core spanning the last 1.5 Myr would provide a transient data benchmark and proxy horizons against which ice-sheet models can be calibrated.
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 are available upon request.
Author contributions. JS devised the experiments, ran the simulations and carried out the analysis. JS wrote the manuscript with contributions and proofreading from all co-authors. The forcing scheme was implemented in PISM by TK and JS. Radar data for Dome Fuji were contributed by NBK and OI. OE and KG coordinated the project.
Competing interests. OE is CEIC of The Cryosphere.
Special issue statement. This article is part of the special issue "Oldest Ice: finding and interpreting climate proxies in ice older than 700 000 years (TC/CP/ESSD inter-journal SI)". It is not associated with a conference.
Acknowledgements. We would like to thank Bas de Boer, Jorge Alvarez-Solas 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 express our gratitude to Christian Stepanek, Madlene Pfeiffer and Martin Werner for providing climate snapshot data from the ESM COSMOS. We are very grateful to Adrien Michel for providing the temperature reconstruction for ensemble branch B1. Hubertus Fischer gratefully acknowledges the long-term support by the Swiss National Science Foundation (SNSF). This publication was generated in the frame of Beyond EPICA -Oldest Ice (BE-OI). It is furthermore supported by national partners and funding agencies in Belgium, Denmark, France, Germany, Italy, Norway, Sweden, Switzerland, the Netherlands and the United Kingdom. Logistic support is mainly provided by AWI, BAS, ENEA and IPEV. The opinions expressed and arguments employed herein do not necessarily reflect the official views of the European Union funding agency, the Swiss Government or other national funding bodies.
Financial support. This publication was generated in the frame of Beyond EPICA -Oldest Ice (BE-OI). The project has received funding from the European Union's Horizon 2020 research and innovation programme (grant no. 730258) (BE-OI CSA). It has received funding from the Swiss State Secretariat for Education, Research and Innovation (SERI) (contract no. 16.0144). The work of Thomas Kleiner has been conducted in the framework of the PalMod project (grant no. 01LP1511B). Development of PISM is supported by NASA (grant no. NNX17AG65G) and NSF (grant nos. PLR-1603799 and PLR-1644277).
The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association.
Review statement. This paper was edited by Alexander Robinson and reviewed by Jorge Alvarez-Solas and Bas de Boer.