Net effect of ice-sheet–atmosphere interactions reduces simulated transient Miocene Antarctic ice-sheet variability

. Benthic δ 18 O levels vary strongly during the warmer-than-modern early and mid-Miocene (23 to 14 Myr ago), suggesting a dynamic Antarctic ice sheet (AIS). So far, however, realistic simulations of the Miocene AIS have been limited to equilibrium states under different CO 2 levels and orbital settings. Earlier transient simulations lacked ice-sheet–atmosphere interactions and used a present-day rather than Miocene Antarctic bedrock topography. Here, we quantify the effect of ice-sheet–atmosphere interactions, running the ice-sheet model IMAU-ICE using climate forcing from Miocene simulations by the general circulation model GENESIS. Utilising a recently developed matrix interpolation method enables us to interpolate the climate forcing based on CO 2 levels (between 280 and 840 ppm), as well as varying ice-sheet conﬁgurations (between no ice and a large East Antarctic Ice Sheet). We furthermore implement recent reconstructions of Miocene Antarctic bedrock topography. We ﬁnd that the positive albedo–temperature


Introduction
The dynamics of the Antarctic ice sheet (AIS) during the early and mid-Miocene (23 to 14 Myr ago) are a topical issue in palaeoclimate science. This issue is of great interest because over the course of this period the climatic conditions, such as the atmospheric CO 2 concentration (Kürschner et al., 2008;Foster et al., 2012;Badger et al., 2013;Greenop et al., 2014;Super et al., 2018;Steinthorsdottir et al., 2021b), ranged from those of the recent past to those predicted for the future (∼ 200 to > 800 ppm; see Steinthorsdottir et al., 2021a, for a review). Benthic δ 18 O records show sizeable orbital-timescale fluctuations during this time, an indication of potential vigorous Antarctic ice volume variability (e.g. Zachos et al., 2008;Liebrand et al., 2011Liebrand et al., , 2017Holbourn et al., 2013;Levy et al., 2019). The δ 18 O signal is, however, also shaped by deep-sea temperature changes. How much of the signal is committed by each component, and thus how the AIS evolved during the Miocene, is currently still under de-bate. On the one hand, geological evidence supports strong AIS dynamism (Pekar and DeConto, 2006;Shevenell et al., 2008), with ice periodically advancing and retreating, e.g. over Wilkes Land  and the Ross Sea sector (Hauptvogel and Passchier, 2012;Levy et al., 2016;Pérez et al., 2021). Some studies, on the other hand, ad-L. B. Stap et al.: Miocene Antarctic ice-sheet variability vocate a more stable small AIS based on ice-proximal seasurface temperature records from dinoflagellate cyst assemblages  and TEX 86 data . In addition, a recent analysis of clumped isotope and Mg/Ca records has found high bottom-water temperatures, which combined with the benthic δ 18 O record paradoxically suggest a Miocene AIS volume persistently larger than in the present day (Modestou et al., 2020). Finally, using climate modelling and data comparisons, it has been suggested that the AIS varied mostly in spatial extent, rather than volume, which imprinted on Miocene δ 18 O variability through a strong effect on deep-water temperatures (Bradshaw et al., 2021).
Ice-sheet models can be used to understand the physical processes driving Miocene AIS variability. Broadly stated, earlier modelling efforts have followed two approaches to this task. One approach is to use ice-sheet models bidirectionally coupled to climate models with as much physical realism in terms of processes included and resolution as is computationally feasible for long-term simulations. Gasson et al. (2016) adopted this approach, using an isotope-enabled set-up of a 3D ice-sheet model coupled to a regional climate model nested in a general circulation model (GCM). They showed the importance of including ice-sheet-climate feedbacks, ice-calving parameterisations, and changes in the ice sheet's oxygen isotopic composition for simulating strong AIS-induced benthic δ 18 O changes. Halberstadt et al. (2021) used the same model set-up to determine the influence of CO 2 changes and uplift of the Transantarctic Mountains on glacial variability during the Miocene. These studies were so far limited to steady-state simulations, in which the climate forcing is kept constant until the ice sheet has equilibrated. Transient ice-sheet behaviour can be very different, as was, for example, demonstrated for the North American and Eurasian ice sheets during the Late Pleistocene (Abe-Ouchi et al., 2013). The second approach, which is the one taken in this study, is therefore to focus on transient AIS changes, applying gradually changing climatic conditions. Earlier studies employed 1D ice-sheet models in combination with parameterised temperatures (De Boer et al., 2010) or relatively simple climate models (Langebroek et al., 2009(Langebroek et al., , 2010Stap et al., 2017). Recently, Stap et al. (2019) showed that orbitaltimescale CO 2 changes cause smaller transient AIS variability than equilibrium simulations suggest by using the 3D thermodynamical Parallel Ice Sheet Model (PISM) forced by Miocene climates obtained from GCM (COSMOS) simulations (Stärz et al., 2017). Conceptual considerations indicate that the disequilibrium between ice volume and CO 2 also leads to out-of-phase changes in these quantities (Stap et al., 2020). However, the COSMOS climate simulations were performed using only one ice-sheet configuration. Hence, the effect of surface whitening (darkening) by ice-sheet advance (retreat) on temperature, the so-called albedo-temperature feedback, was not included, nor was the ice-sheet desertification effect, i.e. the reduction of snow deposition on a growing ice sheet (e.g. Oerlemans, 2004). Furthermore, the effect of mass balance changes due to surface height changes (surfaceheight-mass-balance feedback) could only be crudely captured using a scalar lapse-rate correction for temperature, which is known to be a simplification (Fortuin and Oerlemans, 1990). Moreover, a present-day Antarctic bedrock topography was used, while during the Miocene the topography was different, which also significantly affects the surface mass balance and ice flow (Stap et al., 2016;Colleoni et al., 2018;Paxman et al., 2020).
In this study, we expand the scope of previous studies by performing transient simulations (cf. Gasson et al., 2016;Colleoni et al., 2018;Paxman et al., 2020;Halberstadt et al., 2021), now including ice-sheet-atmosphere interaction (cf. De Boer et al., 2010;Stap et al., 2016Stap et al., , 2017. We run the 3D thermodynamical ice-sheet model IMAU-ICE (De Boer et al., 2013;Berends et al., 2018), which is based on the commonly used shallow ice and shallow shelf approximations and is therefore comparable to PISM. Climate forcing is obtained from Miocene simulations by the GCM GENESIS, in which different CO 2 levels and ice-sheet configurations ranging from zero ice to a substantial East Antarctic Ice Sheet are applied (Burls et al., 2021). Utilising a recently developed matrix interpolation routine (Berends et al., 2018), this set-up enables us to interpolate the climate forcing based on land ice changes in addition to CO 2 changes. In this manner, the effects of ice-sheet changes on local temperatures and precipitation are represented. Because surface height changes are incorporated in the forcing, the surface-height-mass-balance feedback is captured more realistically as well. We perform both equilibrium simulations, keeping the forcing CO 2 level constant, and idealised transient experiments, in which we vary the CO 2 level on long and short quasi-orbital timescales (40, 100, and 400 kyr). We compare experiments using the matrix interpolation method to ones in which we use a simpler index method to obtain climate forcing variability only based on CO 2 changes. In this way, we quantify the influence of the ice-sheet-atmosphere feedbacks on the variability in the AIS during the Miocene. In addition, we investigate the effect of bedrock topography changes over the Miocene by using topographies pertaining to the early (23-24 Myr ago) and mid-Miocene (14 Myr ago) (Hochmuth et al., 2020a;Paxman et al., 2019). We also carry out simulations using Eocene (34 Myr ago) and present-day bedrock topographies. Lastly, we test the sensitivity of our results to the parameterisation of basal melt underneath the floating ice shelves. This is achieved by performing identical simulations to our reference experiment but imposing relatively mild as well as harsh sub-shelf melt rates, easing and prohibiting ice shelf growth respectively. We quantify where and under which conditions ice shelves form and perish in our simulations, as well as how they influence the variability in the grounded ice volume.  Boer et al., 2013;Berends et al., 2018) to simulate the Antarctic ice sheet using atmospheric input data obtained from Miocene climate simulations by the GCM GENESIS (Burls et al., 2021).

Ice-sheet model
IMAU-ICE is a 3D thermodynamical ice-sheet model, which uses a superposition of the shallow ice approximation (SIA) and shallow shelf approximation (SSA) to simulate the dynamics of grounded and floating ice. In the interior of the ice sheet, the ice flow is dominated by the SIA, while towards the fast-flowing ice streams on the fringes the SSA gains importance, becoming the governing component on the ice shelf . No additional grounding-line parameterisations are included. We use a 40 × 40 km resolution grid covering the Antarctic continent. This resolution ensures the feasibility of performing a large number of simulations while still capturing the Antarctic ice flow in the amount of detail we deem appropriate for our purpose. Longitude and latitude coordinates are assigned to the grid using an inverse oblique stereographic projection (Reerink et al., 2010). For further settings, please see Appendix A.

Basal mass balance
Basal melt underneath the ice shelves (M(x, y)) is calculated using a combination of the parameterisation of Pollard and DeConto (2009) and a linear relation to ocean temperature change (Beckmann and Goosse, 2003;Martin et al., 2011): Weighing factor z deep grows with water depths (h w (x, y)) greater than 1800 m: and z expo (x, y) depends on the widest subtended angle to the open ocean (α sub (x, y)) and the shortest linear distance to the open ocean (r open (x, y)): Both factors are clipped at 0 and 1. The melt rate for the exposed shelves (M expo ) varies linearly between 3 m yr −1 for an applied CO 2 concentration of 280 ppm to 6 m yr −1 for CO 2 concentrations of 400 ppm and higher. The melt rate for the deep-water shelves (M deep ) similarly varies between 5 and 10 m yr −1 . Melt rate S m (x, y) is dependent on the water density (ρ w ), specific heat (c p ), thermal exchange velocity (γ T ), water temperature (T water ), pressure-adjusted freezing point (T freeze (x, y)), latent heat of fusion (H fus ), and shelf ice density (ρ i ). Because GENESIS has a slab ocean component (Sect. 2.2), realistic water temperatures cannot be taken from the GCM results Halberstadt et al., 2021). Instead, we use an ad hoc global ocean temperature parameterisation adapted from De Boer et al. (2013). Similar to melt rates M expo and M deep , T w varies linearly as a function of CO 2 between −1.7 • C (280 ppm) and 2 • C (≥ 400 ppm). Currently, calving is neglected. The only calving routine implemented in IMAU-ICE is thickness calving, i.e. calving floating ice with a thickness below a certain threshold. We exclude this in our simulations because it could hinder the growth of new ice shelves (e.g. Ritz et al., 2001). In combination with the heavily parameterised sub-shelf melt calculation, this makes the mass balance of the ice shelves arguably a weak point of our modelling effort. We therefore perform sensitivity experiments (Table 1) to investigate how our results depend on ocean forcing. We consider an extreme case using Last Glacial Maximum (LGM)-like melt parameter settings: M expo = 0 m yr −1 , M deep = 2 m yr −1 , and T w = −5 • C (De Boer et al., 2013). In another experiment, ice shelves are effectively inhibited from growing by applying a constant melt rate of 400 m yr −1 similar to the ABUM experiment of the ice-sheet model intercomparison project ABU-MIP . By comparing the results of these extreme scenarios to our reference experiment, we can quantify where and under which conditions ice shelves form, as well as how they subsequently influence the variability in the grounded ice volume.

Surface mass balance
The surface mass balance (SMB) calculation uses precipitation (P applied (x, y)) and 2 m air temperature (T applied (x, y)) obtained from Miocene climate simulations by GENESIS (Sect. 2.2) utilising a matrix interpolation routine (Sect. 2.3, Appendix B). The GENESIS fields are remapped from their native T31 grid (Sect. 2.2) to the ice-sheet model grid using bilinear interpolation based on longitude and latitude. Monthly SMB is the sum of accumulation and refreezing minus ablation (Berends et al., 2018). Refreezing is dependent on the available liquid water content, and the superimposed water content (Huybrechts and de Wolde, 1999;Janssens and Huybrechts, 2000). Accumulation (Acc(x, y)) is calculated as follows: Acc(x, y) = P applied (x, y) × f snow (x, y), with f snow (x, y) being the temperature-dependent snow fraction (Ohmura, 1999). Ablation (Abl(x, y)) follows from  , isostatically rebounded after ice-sheet removal TOPO Wilson_max Bedrock topography pertaining to late Eocene (Wilson et al., 2012), max reconstruction TOPO Wilson_mean Bedrock topography pertaining to late Eocene (Wilson et al., 2012), mean of min and max reconstructions

BMB LGM
LGM ocean conditions (M expo = 0 m yr −1 , M deep = 2 m yr −1 , and T w = −5 • C) BMB no_shelves Constant sub-shelf melt rate M = 400 m yr −1 an insolation-temperature melt parameterisation (Bintanja et al., 2002): The parameters are set to c 1 = 0.788 m yr −1 K −1 , c 2 = 0.004 m 3 J −1 , and c 3 = −0.50 m yr −1 . The internally calculated monthly albedo (α(x, y)) is derived as follows. First, a background albedo (α bg (x, y)) is determined: 0.1 for oceancovered grid cells, 0.2 for soil, and 0.5 for ice. Thereafter, the effect of snow, with an albedo (α snow ) of 0.85, is added: where Abl yr−1 (x, y) is the ablation calculated for the previous year. The firn layer depth (z firn (x, y), clipped at 0 and 10 m) is changed over time by the addition of snowfall minus melt. The albedo α(x, y) is not allowed to be lower than the background albedo or higher than the albedo of snow. In Eq. (6), incoming solar radiation at the top of the atmosphere Q TOA (x, y) is taken from Laskar et al. (2004) and kept constant at the present-day values so as to isolate climate variability caused by CO 2 changes. This surface mass balance parameterisation has recently been shown to yield realistic results over the present-day Greenland ice sheet (Fettweis et al., 2020).

Climate model
GENESIS (version 3.0) (Thompson and Pollard, 1997;DeConto et al., 2012) has a spectral resolution of T31 (∼ 3.75 • × 3.75 • ). The atmospheric model is coupled to 2 • × 2 • surface models including a non-dynamical 50 m slab diffusive mixed-layer ocean with dynamic sea ice and dynamic vegetation. The applied palaeogeography comes from Herold et al. (2008), with dedicated modifications to Antarctica based on simulations by DeConto and Pollard (2003). Different simulations were carried out as described in Burls et al. (2021). Here, we use simulation 1fumebi for the cold climate and simulation 3nomebi for the warm climate ( Fig. S1 in the Supplement). Simulation 1fumebi has a CO 2 level of 280 ppm (1 × pre-industrial (PI)) and an AIS configuration with a full East Antarctic Ice Sheet (18.975 × 10 6 km 3 volume, area shown in Fig. S2). Simulation 3nomebi is forced by a CO 2 level of 840 ppm (3 × PI) with no ice cover. Both simulations use present-day orbital settings. The vegetation on non-glaciated land, which has an influence on the regional climate, is determined by the dynamic vegetation model BIOME4.

Matrix interpolation method
The applied climate forcing is interpolated between the cold (1fumebi) and warm (3nomebi) GENESIS fields using a matrix interpolation method (Berends et al., 2018), described in detail in Appendix B. Briefly, the applied temperature interpolation is based on a prescribed value for CO 2 and on the internally modelled ice sheets. The higher (lower) the CO 2 level and the smaller (larger) the ice coverage, the more the climate forcing is like the warm (cold) climate. The feedback of the ice sheets on the climate is calculated via the effect on absorbed insolation through changes in surface albedo. The precipitation interpolation is determined solely from the ice volume. Because the cold climate is generally drier, this constitutes a negative ice-volume-precipitation feedback, representing the ice-sheet desertification effect in our simulations.

Simulations
We perform an experiment using our default settings (REF), as well as 11 sensitivity experiments with settings altered as indicated in Table 1. Each experiment consists of 15 simulations, 11 of which are steady state and 4 of which are transient simulations ( Table 2). The steady-state simulations are run for 150 kyr with a constant forcing CO 2 level -enough time for the ice sheet to equilibrate with the forcing climate. Mind that the CO 2 level, rather than the climate itself, is kept fixed, which means the climate can still change due to land ice changes. First, six of these steady-state simulations are run using CO 2 levels bracketing the 280 to 840 ppm range: 280, 392, 504, 616, 728, and 840 ppm. In our reference case, the initial conditions are obtained from recent reconstructions of the Antarctic bathymetry (Hochmuth et al., 2020a) and bedrock topography (Paxman et al., 2019) pertaining to 23 to 24 Myr ago (dataset: Hochmuth et al., 2020b). No ice is present at the start of these runs. These simulations constitute the lower, ascending branch of the hysteresis curve in the equilibrium CO 2 -ice-volume relation (blue connected diamonds in figures, e.g. Fig. 1b). Thereafter, the remaining five steadystate simulations are continued from results of the 280 ppm IMAU-ICE run, with CO 2 levels of 392, 504, 616, 728, and 840 ppm. These simulations constitute the upper, descending branch of the hysteresis curve in the equilibrium CO 2ice-volume relation (red connected diamonds in figures, e.g. Fig. 1b).
Next, we perform the idealised transient simulations (Table 2). Three of these are started from the final state of the equilibrium 840 ppm IMAU-ICE run. The evolution of CO 2 over time has a triangular shape (e.g. pink line in Fig. 1a). The CO 2 is first linearly reduced from 840 to 280 ppm and then increased back to 840 ppm. In the three simulations, this triangular CO 2 evolution is implemented over the course of 100 and 400 kyr and five consecutive times over 40 kyr (200 kyr in total; e.g. pink line in Fig. 1c). The final transient simulation is continued from the 280 ppm IMAU-ICE run and forced with an inverted triangular 40 kyr CO 2 evolution. That means the CO 2 is first linearly increased from 280 to 840 ppm and then decreased back to 280 ppm (e.g. purple line in Fig. 1c).

Reference simulations
As a minimal check for model performance and parameter settings, we first perform steady-state control simulations as in de Boer et al. (2015). These are executed using ERA40 present-day (PD;1957 precipitation and temperature forcing (Uppala et al., 2005). Constant presentday-like melt parameter settings are used to calculate the basal mass balance: M expo = 3 m yr −1 , M deep = 5 m yr −1 , and T w = −1.7 • C. Depending on whether the run is started from BedMachine PD ice and topography  or without ice and an isostatically rebounded topography, we simulate a total ice volume above flotation of 56.6 m sea level equivalent (m s.l.e.; 23.1 × 10 15 m 3 ) or 55.0 m s.l.e. (22.4 × 10 15 m 3 ) respectively. This is very close to the 55.0 m s.l.e. obtained from remapping the BedMachine data to our grid. Like other SIA-and SSA-based ice-sheet models (see de Boer et al., 2015), IMAU-ICE generally simulates slightly thinner ice in the interior and thicker ice along the margins (Fig. S3). The grounding line is simulated quite well, although it is slightly more advanced in the Filchner-Ronne and the Amery embayments. Due to lack of calving, the simulated ice shelves are extended much farther than the observations. However, since the increase in back pressure due to these overestimated shelf areas is very small (Reese et al., 2018b), this is unlikely to affect the evolution of the grounded ice. In this study, we will therefore focus on grounded ice.
In the reference equilibrium simulations using Miocene forcing (experiment REF), the AIS consists only of disconnected glaciers in high mountainous regions in the case of a 840 ppm CO 2 level (Fig. 2a). It grows into a full East AIS (EAIS) and small West AIS (WAIS) in the case of 280 ppm forcing (Fig. 2b). In the latter case, the WAIS grounds on Marie Byrd Land and the peninsula, while it does not grow into the Ross and Filchner-Ronne embayments. Intermediate states for different forcing CO 2 levels are shown in Fig. S4.
We show the reference CO 2 and equilibrium ice volume (above flotation) (V eq ) relation in Fig. 1b (Pollard and De-Conto, 2005). It decreases monotonically and exhibits considerable hysteresis, as displayed by the difference between the ascending (blue) and descending (red) branches. Transient simulations with different imposed timescales for the forcing variability are indicated by the dashed (400 kyr) and dotted (100 kyr) black lines. Peak ice volumes of 45.5 m s.l.e. (400 kyr) and 37.8 m s.l.e. (100 kyr) are yielded. In agreement with Stap et al. (2020), this peak is only reached later in time than the CO 2 minimum (Fig. 1a). The peak is attained relatively earlier in time (Fig. 1a), i.e. at a lower CO 2 level (Fig. 1b), in the 400 kyr simulation than in the shorter 100 kyr simulation. In the case of 40 kyr forcing cycles, the equilibrium cycle is attained very quickly and does not depend on the initial conditions (Fig. 1c, d; black lines partly overlapping). In the equilibrium cycle, the ice volume varies between 3.6 m s.l.e. at 728 ppm and 30.7 m s.l.e. at 532 ppm, a 27.1 m s.l.e. difference. Since the maximum ice volume is reached longer after the minimum CO 2 level (252 ppm difference) than the minimum ice volume after the maximum CO 2 level (112 ppm difference), the growth phase lasts longer than the decay phase. That also means that ice-sheet decay is on average 40 % faster than growth, but both rates, in particular the growth rate, are certainly not constant. In the transient simulations, the AIS initially builds up slowly in  the high mountainous regions at ∼ 0.1 m s.l.e. kyr −1 , speeding up to a maximum of 3 m s.l.e. kyr −1 when the surface mass balance (SMB) becomes positive on a large part of the continent as separate ice domes start to merge (gradients in Fig. 1c). When the transient ice volume is in between the ascending and descending branches of the equilibrium CO 2 -V eq relation, i.e. larger than the ascending (blue) volume but smaller than the descending (red) volume (Fig. 1d), the rate of transient ice volume change is relatively small (< 1 m s.l.e. kyr −1 ). In earlier research using IMAU-ICE, different values for the SMB parameter c 3 in Eq. (6) have been used (e.g. De Boer et al., 2013Boer et al., , 2014Berends et al., 2018). Our reference setting C 3 = 0.50 m yr −1 leads to satisfactory PD results, but using C 3 = 0.25 m yr −1 as De Boer et al. (2014) did, a similar PD ice volume of 55.5 m s.l.e. is obtained (not shown). Using this setting (experiment SMB; Table 1) increases surface melt, with the onset of AIS growth and decay at lower CO 2 levels (Fig. S5) as a result. Our REF settings yield an ice sheet similar in extent to the results of Stap et al. (2019) and DeConto and Pollard (2003) at 280 ppm, which disappears almost completely at 840 ppm. This ensures that we can most effectively study large-scale variability in the East Antarctic Ice Sheet. A different tuning would shift the CO 2 at which these ice-sheet sizes are simulated but does not (qualitatively) affect the shape of the CO 2 -ice-volume relation.

Influence of ice-sheet-atmosphere interaction
Using a matrix interpolation routine enables us to interpolate the climate forcing for our ice-sheet model based on both CO 2 and ice-sheet size. Here, we make a comparison to re- sults of our model set-up using a simpler index routine, as, for example, used in Stap et al. (2019), in which the interpolation is solely based on the CO 2 level. In experiment NOFEEDB (Table 1), we implement w tot = CO 2 −CO 2,cold CO 2,warm −CO 2,cold instead of Eqs. (B7) (for temperature) and (B11) (for precipitation) to interpolate between the cold (1fumebi) and warm (3nomebi) climates.
Because the positive albedo-temperature feedback is not invoked in this case, hysteresis in the CO 2 -V eq relation is smaller (Fig. 3b, solid and dashed red and blue lines). The temperature declines more quickly towards lower CO 2 values, causing the equilibrium simulations to yield larger ice volumes (Fig. 3b, blue lines). At 280 ppm, V eq in the NOFEEDB case is 55.8 m s.l.e. compared to 48.4 m s.l.e. in the REF case. This difference in ice volume is mainly manifested in regions with low bedrock topography through a thicker WAIS in Marie Byrd Land, and EAIS in Mac. Robertson Land and Coats Land (not shown). Because the AIS deglaciates at similar levels as in the REF case (Fig. 3b, red lines), the difference between the CO 2 levels of glaciation and deglaciation is smaller in the NOFEEDB case. That reduces the hysteresis in the CO 2 -V eq relation.
Transient ice volume variability is larger in the NOFEEDB case than in the REF case (Fig. 3a, b, black and orange lines). In the NOFEEDB case, peak ice volumes are increased to 54.3 m s.l.e. (400 kyr forcing CO 2 cycle) and 47.4 m s.l.e. (100 kyr). Forced by 40 kyr CO 2 cycles, the amplitude of variability is 34.1 m s.l.e. in the equilibrium cycle. That means the ice-sheet-atmosphere interactions that are included in the REF case decrease the amplitude of AIS variability by 21 %.
As mentioned above, the difference between the REF and NOFEEDB results is mainly governed by the albedotemperature feedback. However, in the REF case, precipitation decreases with increasing ice volume, via Eq. (B10) combined with a cold climate that is generally drier. This constitutes a negative ice-volume-precipitation feedback.
The effect of this feedback alone is shown by a comparison between the results of experiment NOFEEDB and those of NOFEEDB-P (Table S1, Fig. 3c, d). In the latter experiment, the index method is only used for precipitation, and temperature is interpolated in the normal sense (Eq. B7). The hysteresis in the CO 2 -V eq relation is now larger, and the amplitude of transient variability in the equilibrium cycle of the 40 kyr simulation is 19.3 m s.l.e., 43 % smaller than in the NOFEEDB case. In contrast, using the index method only for temperature (experiment NOFEEDB-T), the amplitude of 40 kyr transient variability becomes 6 % larger (36.1 m s.l.e.) than in the regular NOFEEDB case. The results of experiment NOFEEDB-T are much closer to those of experiment NOFEEDB. Therefore, of the two ice-sheetatmosphere feedbacks, the albedo-temperature feedback is dominant in our experiments.

Effect of bedrock topography changes
Recent geological reconstructions have indicated that Antarctic bedrock topography has evolved substantially since the Eocene-Oligocene boundary (34 Myr ago) mostly due to thermal subsidence and glacial erosion (Paxman et al., 2019). In the REF simulations, we have used the bedrock topography reconstruction pertaining to the Oligocene-Miocene boundary (23-24 Myr ago). In experiment TOPO 14 Ma (Table 1), we instead use the one pertaining to the mid-Miocene (14 Myr ago). A larger part of West Antarctica (Marie Byrd Land and Palmer Land) has subsided below sea level in this reconstruction (Fig. S6a, b). This leads to a smaller simulated WAIS at low CO 2 levels in the TOPO 14 Ma case (Fig. 4). The bedrock subsidence in Victoria Land and the Amery embayment has a smaller effect on equilibrium ice volume. The total equilibrium ice volume above flotation is 44.6 m s.l.e. at 280 ppm compared to 48.4 m s.l.e. in the REF case (−8 %). The amplitude of transient ice volume variability is affected as well, even on relatively short orbital (obliquity) timescales. It decreases to 24.4 m s.l.e (−10 % with respect to REF) when 40 kyr forcing cycles are imposed (Fig. 5).
The effect of long-term Antarctic landscape evolution on ice volume becomes more apparent when we compare simulations using bedrock conditions of the Eocene-Oligocene boundary (34 Myr ago; TOPO 34 Ma) and relaxed PD topography (TOPO PD) (Fig. S6c, d). In the latter case, less ice builds up in Victoria Land and Coats Land where the initial bedrock topography has subsided below sea level compared to the TOPO 14 Ma case (Fig. S7a) S7b) and is centred around a generally larger ice volume in experiment TOPO 34 Ma (Fig. S8a, b, long-dashed lines). This is also  due to the fact that the EAIS is more stable, as shown by larger ice volumes already at high CO 2 levels in the CO 2 -V eq relation, particularly in the upper, descending branch. Hysteresis is indeed larger than in the REF case because the positive surface-height-mass-balance feedback is stronger due to the build-up of thicker ice. At 280 ppm, the equilibrium ice volume is 57.4 m s.l.e. in this experiment (+19 % with respect to REF). A much thicker ice sheet is built up in Coats Land, the Amery embayment, and Marie Byrd Land (Fig. S7b). Still in this case, the ice sheet does not significantly advance into the Ross and Filchner-Ronne embayments. This is different when we use the earlier geological reconstructions of Antarctic palaeotopography at the Eocene-Oligocene boundary of Wilson et al. (2012) (Fig. S6e, f). At 280 ppm, the ice then grounds in the Ross embayment (TOPO Wilson_mean, using the mean of the minimum and maximum reconstructions; Fig. S7c) or in both embayments (TOPO Wilson_max, using the maximum reconstruction; Fig. S7d). Due to the higher CO 2 levels at which deglaciation takes place, a relatively large ice volume persists throughout 40 kyr transient simulations (Fig. S8c, d). In contrast to the other experiments, the initial conditions of these transient simulations now affect the ice volume during the equilibrium cycle. The ice sheet remains generally larger when the equilibrium 280 ppm simulation, rather than the 840 ppm simulation, is used as starting point (not shown).

Sensitivity to ocean forcing
We perform and compare two additional sensitivity experiments that are end-members with respect to basal melt rate settings to access the sensitivity of our results to the parameterisation of ocean forcing. In experiment BMB LGM (Table 1), we use Last Glacial Maximum (LGM)-like basal melt parameter settings in Eq. (1): M expo = 0 m yr −1 , M deep = 2 m yr −1 , and T w = −5 • C (De Boer et al., 2013). The mild basal melt rates in this case create an environment in which ice shelves can grow relatively easily. At the other extreme, in experiment BMB no_shelves, ice shelves are effectively inhibited from growing by applying a constant melt rate of 400 m yr −1 . The effect of ice shelf formation becomes notable from 504 ppm in the equilibrium simulations (Fig. 6b, solid and dashed blue lines). Ice shelves form along the margins at Coats Land and Oates Land (Fig. 7a). While in Coats Land this leads to a modest thickening of the grounded ice, this happens to a much lesser degree in Oates Land. In the long 400 kyr transient simulation, these ice shelves form after around 110 kyr at ∼ 530 ppm (Fig. 6b). At 392 ppm in the equilibrium simulations, and after around 170 kyr at ∼ 360 ppm in the 400 kyr transient simulation, ice shelves appear approximately simultaneously along almost the entire remainder of the (east and west) Antarctic margin (Fig. 7b). The effect on grounded ice is particularly noticeable in West Antarctica, as well as Wilkes Land, George V Land, and Mac. Robertson Land, where the ice sheet thickens by kilometres. In the deglaciation phase these ice shelves remain influential until 728 ppm in the equilibrium simulations (Fig. 6b, solid and dashed blue lines) and until around 360 kyr at ∼ 728 ppm in the 400 kyr transient simulation (Fig. 6a, b, green and orange lines). In the transient simulations, the amplitude of variability is increased. The hysteresis in the CO 2 -V relation is increased as well, although only slightly. In the 40 kyr transient simulations, the ice volume variability amounts to 30.7 m s.l.e. in case LGM BMB and 25.7 m s.l.e. in case LGM no_shelves, a 5 m s.l.e. (19 %) difference (Fig. 6c, d).

Discussion
The CO 2 levels at which East Antarctic glaciation is simulated are known to be dependent on the climate model used to construct the climate forcing (Gasson et al., 2014). This could (partly) explain why in our study EAIS inception occurs at higher CO 2 levels than in Stap et al. (2019). They used climate forcing from COSMOS which has a higher climate sensitivity (reported as 4.1 K per CO 2 doubling) than GENESIS (2.9 K). Note that in neither study are greenhouse gases other than CO 2 varied, probably biassing the CO 2 concentrations high. Another factor that has a large influence on the glaciation threshold is the ablation factor that we use in our surface mass balance calculation (C 3 in Eq. 6), as demonstrated by our SMB sensitivity experiment (Fig. S5). Arguably the best way to further constrain this parameter in the future is by performing a comprehensive comparison of our mass balance calculation, and other methods and observations, as was done previously for the Greenland ice sheet (Fettweis et al., 2020). Differences in ice dynamics, e.g. in the treatment of ice flow enhancement and basal sliding, might cause further disparity. Probably for the same reasons, we obtain larger transient AIS volume variability than Stap et al. (2019) when using a similar index method (experiment NOFEEDB). We nonetheless corroborate that disequilibrium between the AIS and the forcing climate causes smaller transient ice volume variability than steady-state simulations suggest (Stap et al., 2019), as well as out-of-phase changes between these quantities (Stap et al., 2020).
Here, we include ice-sheet-atmosphere feedbacks in our reference case. We have shown that these feedbacks tend to increase hysteresis in the CO 2 and equilibrium ice volume relation and decrease the amplitude of AIS variability under equal 40 kyr forcing CO 2 cycles by 21 % (Fig. 3). Also on longer timescales AIS variability is suppressed by the ice-  sheet-atmosphere feedbacks. This is partly due to the smaller equilibrium ice volumes at low CO 2 levels and partly to decreased growth rates of the transient ice volume relative to the change in equilibrium ice volume. In short, the ice-sheetatmosphere feedbacks cause a slower build-up of the AIS to smaller peak ice volumes. Correspondingly, the contribution of AIS volume changes to benthic δ 18 O fluctuations gets smaller. This calls for a larger contribution of deep-sea tem-peratures to the amplitude of Miocene benthic δ 18 O fluctuations. The actual strength of the impact will, however, depend on the isotopic composition of Antarctic snow as well Rohling et al., 2021).
Furthermore, we use recent reconstructions of Antarctic bedrock topography during the Miocene (Hochmuth et al., 2020a;Paxman et al., 2019), rather than the present day as in Stap et al. (2019), as a boundary condition in our AIS simulations. In agreement with Colleoni et al. (2018) and Paxman et al. (2020), we find that the WAIS becomes more vulnerable over the early and mid-Miocene (Fig. 5). This is demonstrated by the lower equilibrium ice volumes at equal CO 2 levels in our experiment TOPO 14 Ma, which has a bedrock topography pertaining to the mid-Miocene (14 Myr ago) compared to REF, in which the topography is from the Oligocene-Miocene boundary (23-24 Myr ago). We also corroborate the finding of Paxman et al. (2020) that from 34 to 24 Myr ago bedrock changes primarily concern the EAIS, while from 24 Myr ago onwards mostly the WAIS is affected. Here, we have additionally examined the effect of the evolving Antarctic landscape on transient Miocene AIS variability, which is 10 % smaller when subjected to 40 kyr forcing CO 2 cycles in TOPO 14 Ma than in REF due to impeded WAIS growth. Hence, CO 2 minima must have dropped increasingly further over the course of the early and mid-Miocene to obtain AIS variability with similar amplitude and consequently a similar contribution of the AIS to Miocene benthic δ 18 O variability.
In this study, we focus on the variability in a generally small AIS, which is located primarily on the east side of the continent, while the WAIS size remains relatively limited. In the Ross and Filchner-Ronne embayments growth of extended ice shelves never occurs in our Miocene simulations, not even when the basal melt rates are set to LGM values. On the one hand, it is certainly possible that our basal melt parameterisation is not suited for the simulation of (largescale) Antarctic ice shelf inception. On the other hand, the fact that the Ross and Filchner-Ronne ice shelves do grow when we impose steady present-day climate forcing points out that West Antarctic ice shelves have to be sustained by a sufficiently large ice flow coming from the WAIS. That means that in addition to the oceanic conditions, suitable atmospheric conditions are very important for ice shelf formation as well. Therefore, WAIS growth is hampered by the air temperatures remaining too warm in our forcing climate simulations because in the cold case (simulation 1fumebi) there is only a small WAIS present (Fig. S2) in the forcing climate simulation. In the future, the forcing steady-state climate simulations used in the matrix interpolation should preferably be constructed by a coupled set-up of the same ice-sheet and climate models that are deployed in the transient simulations. This would ensure congruence between the results of the equilibrium and transient simulations so that the transient ice sheet will not grow beyond its largest size in the forcing equilibrium results.
Although the AIS remains mostly grounded, the formation of floating ice shelves does have an impact on our result. This is demonstrated by comparing simulations in which LGM basal melt rates are applied to simulations in which ice shelf growth is prevented (Figs. 6, 7). The 40 kyr transient AIS variability is amplified by 19 % in the case of LGM basal melt rates -and 5 % in the reference case -compared to the case without ice shelves. In our simulations, ice shelves first appear along the coasts of Oates Land and Coats Land when the CO 2 concentration drops below ∼ 530 ppm (T w = 2 • C), but the effect of ice shelves on grounded ice volume only becomes significant below ∼ 360 ppm CO 2 (T w = 0.8 • C) when the EAIS -with a volume of almost 30 m s.l.e. -is already well-established. This transition, from a land-based ice sheet to one that is in contact with the ocean, at CO 2 levels around 400 ppm is in general agreement with other model results (Halberstadt et al., 2021) and geological data inferences (Naish et al., 2009;Levy et al., 2016Levy et al., , 2019. Below ∼ 360 ppm CO 2 , ice shelves fringe almost the entire continent in our simulations. This includes the Wilkes Subglacial Basin and the Ross Sea sector, where records from IODP (Integrated Ocean Drilling Program) site U1356 and ANDRILL AND-2A show evidence of periodic ice retreat and readvance (Hauptvogel and Passchier, 2012;Levy et al., 2016;Pérez et al., 2021;Sangiorgi et al., 2018). When we increase CO 2 levels from cold conditions in our simulations, the ice shelves keep having an effect on grounded AIS ice volume until relatively high CO 2 levels of ∼728 ppm (T w = 2 • C) are reached. Only then does the AIS retreat into the Wilkes Subglacial Basin, suggesting that in our model set-up periodic AIS waxing and waning in this region have to be governed by substantial CO 2 fluctuations. Growth of ice shelves hence increases hysteresis in the CO 2 -ice-volume relation, in general agreement with the findings of Garbe et al. (2020). In line with this result, we find a different effect on final ice volume (above flotation) in steady-state PD simulations with a 400 m yr −1 melt rate depending on whether the simulation is started with or without ice. Starting from an isostatically rebounded PD topography without ice, we obtain a volume of 39.5 m s.l.e., 28 % smaller than the PD run with normal basal melt rates (Fig. S9a, c). When the simulation is started from a PD topography including grounded and floating ice, the final volume is 48.4 m s.l.e., which is only 14 % smaller than the normal PD run (Fig. S9b, d). Hence, in our model, the effect of ice shelves on grounded ice volume is smaller when existing ice shelves are melted away than when ice shelves cannot form in the first place. Despite having a longer model runtime and different surface mass balance forcing, these simulations are similar to the ABUMIP ABUM experiment . Our simulated 8.2 m s.l.e. loss of ice volume above flotation after ice-sheet removal under PD forcing conditions (starting with ice) is on the large side of the wide range of model responses (∼ 1 to 11 m s.l.e. loss) reported there. Showing complete collapse of the WAIS and strong retreat in the Amery embayment but limited ice loss in the Recovery, Wilkes, and Aurora subglacial basins, our spatial response is most similar to that of ABUMIP models PISM1, PSU3D1, and BISICLES. However, our simulated response, using a Coulomb sliding law (q = 0.3), is weaker than that of Elmer/Ice and SICOPOLIS (Simulation Code for Polythermal Ice Sheets) which used a Weertman sliding law (m = 3) and that of a different version of IMAU-ICE that used the Schoof flux boundary condition (Schoof, 2007) at the grounding line. A stronger response of grounded ice volume to ocean warming and ice shelf break-up might lead to earlier retreat, i.e. at a lower CO 2 level, of the AIS in, for example, the Wilkes Subglacial Basin, hence reducing hysteresis in the CO 2 -ice-volume relation.
Our climate input comes from the same model (GENE-SIS) as Halberstadt et al. (2021), who obtain relatively large EAIS volumes at CO 2 levels up to 1140 ppm. Differences between their results and ours can be explained by the different ice-sheet models and surface mass balance schemes used. We use an insolation temperature melt (ITM) scheme (Bintanja et al., 2002), accounting for the direct effect of insolation on surface melt, whereas they deploy a positive degree day (PDD) scheme. Additionally, they use a regional climate model for downscaling the forcing, which allows for capturing orographic precipitation in more detail. We use a bias correction to obtain a more detailed reference climate. This may, however, lead to overestimated temperatures be-cause we use a present-day (ERA-40) base climate, which is slightly warmer than our PI reference (Stap et al., 2019). In model tests, we simulate higher AIS ice volumes at high CO 2 levels when we refrain from using this bias correction, although this is dependent on the ablation factor (c 3 in Eq. 6; not shown). All in all, even with the bias correction included, we need very large CO 2 variations on relatively short orbital timescales (40 kyr) to obtain substantial East Antarctic Ice Sheet variability. This could point to a relatively stable EAIS and a consequent reduced contribution to benthic δ 18 O fluctuations during the Miocene Hartman et al., 2018). Alternatively, the WAIS could play a more major role in establishing ice-sheet variability. As explained above, WAIS growth is hampered in our simulations by the climate forcing. However, geological evidence was recently found for a dynamic Miocene WAIS in the Ross Sea section , in agreement with the modelling results of, for example, Gasson et al. (2016) and Halberstadt et al. (2021). The generally low benthic δ 18 O values during the Miocene in that case would require persistent high ocean temperatures, which were indeed recently found by Modestou et al. (2020). Further modelling efforts are still needed to reconcile these warm ocean temperatures with the existence of a large marine AIS.

Conclusions
Using the ice-sheet model IMAU-ICE, we have performed steady-state and transient simulations of the Antarctic ice sheet (AIS) during the early and mid-Miocene (23 to 14 Myr ago). Climate forcing was obtained from Miocene simulations by GENESIS (Burls et al., 2021), with different CO 2 levels and ice-sheet configurations, using a recently developed matrix interpolation routine (Berends et al., 2018).
Compared to a simpler index method, as used, for example, in Stap et al. (2019), our matrix interpolation method incorporates two important ice-sheet-atmosphere feedbacks: the positive albedo-temperature feedback and a negative icevolume-precipitation feedback. The former tends to increase hysteresis in the CO 2 and ice volume above flotation (V ) relation and decrease the amplitude of transient AIS variability. The negative ice-volume-precipitation feedback, which represents the ice-sheet desertification effect, has the opposite effect (Fig. 3). The albedo-temperature feedback is dominant in our simulations. In the 40 kyr transient reference simulation, the ice volume variability amplitude is 21 % smaller than when an index method is used. Hence, the net effect of accounting for ice-sheet-atmosphere interactions in the interpolation of climate forcing reduces simulated transient Miocene Antarctic ice-sheet variability.
Differences between the geological reconstructions of Antarctic bedrock topography at the Oligocene-Miocene boundary (23-24 Myr ago) and the mid-Miocene (14 Myr ago) are mainly located in West Antarctica. The subsidence of land below sea level during the early and mid-Miocene impedes WAIS build-up. Consequently, it diminishes hysteresis in the CO 2 -V relation and reduces 40 kyr AIS variability by 10 % in amplitude (Figs. 4,5). Hence, CO 2 minima must have dropped increasingly further over the course of the early and mid-Miocene to obtain AIS variability with similar amplitude.
Here, we have conducted idealised transient simulations, addressing the effect of CO 2 variability on Miocene AIS evolution. Moving forward towards realistic transient simulations of AIS evolution, orbital forcing variability will need to be accounted for in the model set-up, for instance by integrating it in the matrix interpolation method (e.g. Ladant et al., 2014;Tan et al., 2018). A further improvement to be made is the implementation of a more advanced scheme to calculate the basal mass balance underneath the ice shelves; several alternatives have been developed in recent years (e.g. Reese et al., 2018a;Lazeroms et al., 2019;Pelle et al., 2019). This should be accompanied by ocean forcing supplied by sophisticated (3D) ocean models. As a final prospect, inclusion of a coupled sediment model will ultimately open up the possibility of quantifying the transient effect of changing bedrock topography on AIS dynamics over the Miocene (Pollard andDeConto, 2003, 2020).

Appendix A: Ice-sheet model settings
In IMAU-ICE, Glen's flow law exponent is set to n = 3, and the SIA-and SSA-based velocities are enhanced by factors of 5 and 1 respectively. The basal stress is calculated using a Coulomb power law model for pseudo-plastic till (exponent q = 0.3 and threshold velocity u threshold = 100 m s −1 ) (Bueler and Brown, 2009). The till friction angle, which determines the yield stress, varies linearly between 5 • for bedrock elevations of −1000 m and lower and 20 • for those of 0 m and above. The geothermal heat flux is taken from Shapiro and Ritzwoller (2004). A basic elastic lithosphererelaxing asthenosphere (ELRA) model (Lingle and Clark, 1985;Le Meur and Huybrechts, 1996) is used to calculate bedrock adjustment, using a mantle density of 3300 kg m −3 , a lithospheric flexural rigidity of 1.0 × 10 25 kg m 2 s −2 , and a relaxation time of 3000 years. The sea level is kept constant at present-day level (0 m) because local Antarctic sea level changes are very uncertain during the Miocene as they are mainly due to Antarctic ice-sheet changes. Local sea level changes can differ substantially from the global mean (e.g. Stocchi et al., 2013), and their calculation would require solving the sea level equation which is not yet available in IMAU-ICE. We have performed three tests applying a constant 7 m higher sea level (mimicking the absence of the Greenland ice sheet), as well as applying internally calculated global mean and minus the global mean sea level changes. These tests have shown negligible differences in comparison to our reference results (not shown).

Appendix B: Climate matrix forcing
Using a so-called matrix method (Pollard, 2010;Pollard et al., 2013), a time-continuous climate forcing for the icesheet model can be calculated. This is achieved by interpolating between the GENESIS time-slice climate simulations every 10 model years, based on the varying CO 2 levels and ice coverage. The matrix interpolation routine used in this study is described in detail in Berends et al. (2018), and it was previously deployed to simulate the evolution of the North American, Eurasian, Greenland, and Antarctic ice sheets during the past 3.6 Myr (Berends et al., 2019.

B1 Surface air temperature
For the applied surface air temperature, first we calculate a scalar weighing factor for CO 2 (w CO 2 ) and a spatially varying one for ice-sheet coverage (w ice (x, y)).
where CO 2 is the applied CO 2 level, and CO 2,cold and CO 2,warm are the low (280 ppm) and high (840 ppm) CO 2 levels applied in the climate forcing simulations. Weight factor w ice (x, y) is based on scaling between local absorbed insolation of the forcing climates, accounting for the effect of surface whitening (darkening) by ice-sheet advance (retreat) on temperature. In order to obtain w ice (x, y), first the absorbed insolation I abs (x, y) is calculated from the incoming insolation (Q TOA (x, y)) and internally calculated surface albedo (α(x, y)). This is done every 10-year coupling time step for the modelled field (I abs,mod ). For the cold (1fumebi) and warm (3nomebi) GCM snapshot fields (I abs,cold , I abs,warm ), this is done at the start of each simulation by running the SMB model (Sect. 2.1.2), without ice dynamics, for 10 years to determine α(x, y).
Incoming insolation Q TOA (x, y) (Laskar et al., 2004) is kept constant at the present-day values so that differences between the cold, warm, and modelled I abs (x, y) fields are solely determined by the albedo. Then, the weighing field for insolation (w ins (x, y)) is calculated: w ins (x, y) = I abs,mod (x, y) − I abs,cold (x, y) I abs,warm (x, y) − I abs,cold (x, y) .
Weighing factors w CO 2 and w ins (x, y) are clipped at −0.25 and 1.25 to avoid spurious extreme temperatures and precipitation. Finally, the weighing field for ice-sheet coverage is calculated from the smoothed (using a Shepard smoothing filter with a radius of 200 km) and average (over the ice-sheet model domain) w ins to account for local as well as regional effects: w ice (x, y) = 1 7 w ins,smoothed (x, y) + 6 7 w ins,avg .
The viability of the partitioning w ins,smoothed (x, y) and w ins,avg was demonstrated by Berends et al. (2018), who simulated ice-sheet evolution over the last glacial cycle, showing that model results agree well with available data in terms of ice-sheet extent, sea level contribution, ice-sheet surface temperature, and contribution to benthic δ 18 O. For each grid cell in the domain the final weighing factor w tot (x, y) is calculated as follows: w tot (x, y) = × w CO 2 + (1 − ) × w ice (x, y), where is the ratio between the effects of CO 2 and ice-sheet coverage. To determine for our simulations, we use GENE-SIS simulations 3fumebi (840 ppm CO 2 , large ice sheet) and 1nomebi (280 ppm CO 2 , no ice sheet). These climate simulations isolate the separate effects of CO 2 and the ice sheet respectively. The (spatially averaged) relative influence of CO 2 on temperature is estimated as the isolated effect of CO 2 divided by the combined influence of CO 2 and ice-sheet coverage.
where T is the average temperature over the Antarctic domain. Synergy between the effects of CO 2 and icesheet coverage is relatively small so that ≈ (T 3fumebi − T 1fumebi )/(T 3nomebi − T 1fumebi ). Next, the reference temperature (T ref (x, y)) and surface height (h ref (x, y)) are obtained by interpolating between the cold (T cold (x, y), h cold (x, y)) and warm (T warm (x, y), h warm (x, y)) forcing climate fields: T ref (x, y) = w tot (x, y) × T warm (x, y) + (1 − w tot (x, y)) × T cold (x, y), and likewise h ref (x, y) = w tot (x, y) × h warm (x, y) The applied 2 m air temperature field T applied (x, y) is obtained by first correcting for the height differences between the modelled (h mod (x, y)) and reference (h ref (x, y)) surface height using a uniform lapse rate λ = 8 K km −1 . Finally, a spatially varying bias correction based on the difference between the modelled pre-industrial temperature (T PI ) and the 1957-2002 ERA40 reanalysis data (Uppala et al., 2005) (T ERA40 ) is applied: T applied (x, y) = T ref (x, y) − λ × (h mod (x, y) − h ref (x, y)) + (T ERA40 (x, y) − T PI (x, y)).

B2 Precipitation
Precipitation is weighed based on the total ice volumes modelled (V mod ) and in the cold (1fumebi) and warm (3nomebi) GENESIS forcing climates (V cold , V warm ): A reference precipitation field is obtained by interpolating between the cold (P warm ) and warm (P cold ) forcing fields: P ref (x, y) = exp w tot × log[P warm (x, y)] Similarly as for the temperature, a bias correction is applied: P ref,corr (x, y) = P ref (x, y) × P ERA40 (x, y) + 0.008 P PI (x, y) + 0.008 .
Author contributions. LBS designed the research and performed the experiments, with technical assistance from CJB. LBS, CJB, MDWS, and RSWvdW analysed the results. EGWG provided the GENESIS results used as climate forcing in this study. LBS drafted the paper, with input from all co-authors.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Acknowledgements. We thank our colleague Peter Bijl for commenting on an earlier draft of this article. We further thank the researchers that reconstructed the Antarctic palaeotopographies and palaeobathymetries and making them publicly available, in particular Isabel Sauermilch for providing them to us. Simulations were performed on the Gemini computing cluster of the Faculty of Science, Utrecht University. Last but not least, we would like to thank two anonymous reviewers, whose comments have helped to improve the quality of this article.
Financial support. Lennert B. Stap is funded by the Dutch Research Council (NWO) through VENI grant VI.Veni.202.031. Constantijn J. Berends is supported by PROTECT, which has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 869304 (PROTECT contribution number 31).
Review statement. This paper was edited by Alexander Robinson and reviewed by two anonymous referees.