Competing influences of the ocean, atmosphere and solid earth on transient Miocene Antarctic ice sheet variability

Benthic δO levels vary strongly during the warmer-than-modern earlyand 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 CO2 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 icesheet-atmosphere interactions, running IMAU-ICE using climate forcing from Miocene simulations by the general circulation 5 model GENESIS. Utilising a recently developed matrix interpolation method enables us to interpolate the climate forcing based on CO2 levels (between 280 and 840 ppm) as well as ice sheet configurations (between no ice and a large ice sheet). We furthermore implement recent reconstructions of Miocene Antarctic bedrock topography. We find that the positive albedo-temperature feedback, partly compensated by the negative ice-volume-precipitation feedback, increases hysteresis in the relation between CO2 and ice volume (V). Together, these ice-sheet-atmosphere interactions decrease the amplitude of AIS variability caused 10 by 40-kyr forcing CO2 cycles by 21% in transient simulations. Thereby, they also diminish the contribution of AIS variability to benthic δO fluctuations. Furthermore, we show that under equal atmospheric and oceanic forcing, the amplitude of 40-kyr transient AIS variability becomes 10% smaller during the earlyand mid-Miocene, due to the evolving bedrock topography. Lastly, we quantify the influence of ice shelf formation around the Antarctic margins, by comparing simulations with Last Glacial Maximum (LGM) basal melt conditions, to ones in which ice shelf growth is prevented. Ice shelf formation increases 15 hysteresis in the CO2-V relation, and amplifies 40-kyr AIS variability by 19% using LGM basal melt rates, and by 5% in our reference setting.


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 (paleo)climate science. This issue is of great interest, because the climatic conditions of this period, such as the atmospheric 20 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, 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 ice volume 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 albedo-temperature and ice-volume-precipitation feedbacks are represented, and because surface height changes are incorporated in the forcing, the surface-height-mass-balance feedback is captured more 70 realistically as well. Performing several experiments that include both equilibrium and transient simulations, we test the influence of these feedbacks under varying CO 2 levels. Furthermore, we quantify the effect on transient Miocene AIS volume evolution of using bedrock topographies (Hochmuth et al., 2020a;Paxman et al., 2019) pertaining to the early-(23-24 Myr ago) and mid-Miocene (14 Myr ago), as well as the Eocene (34 Myr ago) and present day. Lastly, we investigate where and under which conditions ice shelves form, and what their influence on the grounded ice volume is. This is achieved by performing and 75 comparing simulations, in which relatively mild, and severe basal melt rates underneath the ice shelves are applied.

Models and methods
For this study, we employ the ice-sheet model IMAU-ICE (v1.1.1-MIO) (De 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). 80

Ice-sheet model
IMAU-ICE is a threedimensional 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 . We use a 40 x 40-km resolution 85 grid covering the Antarctic continent. This resolution ensures the feasibility of performing a large numbers of simulations, while still capturing the Antarctic ice flow in the amount of detail we deem appropriate for our purpose. For further settings, please see Appendix A

Basal mass balance
Basal melt underneath the ice shelves (M) is calculated using a combination of the parameterisation of Pollard and DeConto 90 (2009) and a linear relation to ocean temperature change (Beckmann and Goosse, 2003;Martin et al., 2011): with Here . The melt rate for the exposed shelves (M expo ) varies linearly between 3 m/yr for an applied CO 2 concentration of 280 ppm to 6 m/yr 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. Melt rate S m is dependent on the water density (ρ w ), specific heat (c p ), thermal exchange velocity (γ T ), water temperature (T water ), pressure-adjusted freezing point (T f reeze ), latent heat of fusion (H f us ) and shelf ice density (ρ i ). Similar to melt 100 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. 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. Therefore, in a separate sensitivity experiment (Table 1), we consider an extreme case using Last Glacial Maximum (LGM)-like melt parameter settings: M expo = 0 m/yr,

105
M deep = 2 m/yr, 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 similar to the ABUM experiment of the ice-sheet model intercomparison project ABUMIP (Sun et al., 2020).

Surface mass balance
The surface mass balance (SMB) calculation uses precipitation (P applied ) and 2-m air temperature (T applied ), obtained from 110 Miocene climate simulations by GENESIS (Sect. 2.2) utilising a matrix interpolation routine (Appendix B). 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) is calculated as:
For our reference simulations, we use the 2-m air temperature, precipitation, and surface height fields of simulations 1fumebi (cold) and 3nomebi (warm), bilinearly interpolated to the ice-sheet model grid (Fig. S1). Simulation 1fumebi has 280 ppm CO 2

135
(1 x pre-industrial (PI) level), a full ice sheet, and medium Antarctic summers, while simulation 3nomebi has 840 ppm CO 2 (3 x PI), no ice sheet, and medium Antarctic summers. The suffix 'bi' in the names of the GENESIS simulations indicates that dynamic vegetation model BIOME4 was used, which is of no direct concern for the current study.

140
We perform an experiment using our default settings (REF), as well as 12 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 equili-145 brate 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. 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 million years (Myr) ago (dataset: Hochmuth et al., 2020b). No ice is present at the start of these runs. The CO 2 levels are set to different levels bracketing the 280 to 840 ppm range: 280, 392, 504, 616, volume relation (blue lines in figures, e.g. Fig. 2b). Thereafter, the remaining five steady-state 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 2 -ice volume relation (red lines in figures, e.g. Fig. 2b).

155
Next, we perform the transient simulations. Three of these are continued from the 840-ppm IMAU-ICE run. The CO 2 is first linearly reduced from 840 to 280 ppm, and then increased back to 840 ppm, over the course of 100, 400 kyr, and five consecutive times over 40 kyr (200 kyr total). The final transient simulation is continued from the 280-ppm IMAU-ICE run.
In this simulation, CO 2 is first gradually increased from 280 to 840 ppm, and then reduced back to 280 ppm five consecutive times over 40 kyr. This is essentially the normal transient 40-kyr run mirrored.

Reference simulations
As a minimal check for our model and parameter settings, we first perform steady-state simulations using ERA40 present-day (PD;1957-2002 precipitation and temperature forcing (Uppala et al., 2005). Depending on whether the run is started from 165 Bedmachine PD ice and topography (Morlighem et al., 2020) or without ice and an isostatically rebounded topography, we sim- 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 the extra shelf area is very small (Reese et al., 2018a), 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 case of a 840-ppm CO 2 level. It grows into a full East AIS (EAIS) and small West AIS (WAIS) in case of 280-ppm forcing (Fig. 1). In the latter case, the WAIS grounds on Marie Byrd Land and the Peninsula, while not growing into the Ross and Filchner-Ronne embayments. The reference CO 2 -equilibrium ice volume (above floatation) (V eq ) relation decreases monotonically and shows considerable hysteresis (Fig. 2b). Depending on the timescale of the kyr). In agreement with Stap et al. (2020), this peak is only reached after the CO 2 minimum (Fig. 2a), and relatively earlier, i.e. at a lower CO 2 level, in the 400-kyr simulation than in the shorter 100-kyr simulation. In case of 40-kyr forcing cycles, the equilibrium cycle is attained very quickly and does not depend on the initial conditions (Fig. 2c,d). In the equilibrium cycle, the ice volume varies between 3.6 m.s.l.e. at 728 ppm to 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

Influence of interaction between ice sheets and atmosphere
Using a matrix interpolation routine enables us to interpolate the climate forcing for our ice-sheet model based on both CO 2 and ice volume. Here, we make a comparison to results of our model set-up using a simpler glacial index routine, as e.g. used in Stap et al. (2019), in which the interpolation is solely based on the CO 2 level. In experiment FEEDB (Table 1), we implement 215 w tot = CO2−CO 2,cold CO2,warm−CO 2,cold instead of Eqs. B7 (for temperature) and B11 (for precipitation) to interpolate between the cold and warm climates. Because the positive albedo-temperature feedback is not invoked in this case, hysteresis in the CO 2 -V eq relation is smaller (Fig. 4b). The temperature declines more quickly towards lower CO 2 values, causing the equilibrium simulations to yield larger ice volumes. At 280 ppm, the difference in ice thickness with the REF simulation is again mainly manifested in regions with low bedrock topography (not shown). This also results in larger transient ice volume variability (Fig. 4a,b). Forced 220 by 40-kyr CO 2 cycles, the amplitude of variability is 34.1 m.s.l.e. in the equilibrium cycle, 26% larger than in the REF case.
Stated the other way around, the ice-sheet-atmosphere interactions decrease the amplitude of AIS variability by 21%.
As mentioned before, the difference between the REF and FEEDB results is mainly governed by the albedo-temperature feedback. However, in the REF case, precipitation decreases with increasing ice volume, via Eq. B10 combined with a drier climate at lower CO 2 . This constitutes a negative ice-volume-precipitation feedback. The effect of this feedback alone is shown 225 by a comparison to the results of experiment FEEDB onlyP (Table S1, Fig. 4c,d). In this experiment, the glacial index method is only used for precipitation, and temperature is interpolated in the normal manner (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., 29% smaller than in the REF case. In contrast, using the glacial index method only for temperature (experiment FEEDB onlyT), the amplitude of 40-kyr transient variability becomes 6% larger (36.1 m.s.l.e.) than in the regular FEEDB case. In conclusion, 230 the albedo-temperature feedback is the strongest of these two ice sheet-atmosphere feedbacks in our experiments.

Influence of bedrock topography
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 235 ago). In experiment TOPO 14Ma (Table 1) 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 34Ma) and relaxed PD topography (TOPO PD). In the latter case, less ice builds up in Victoria Land and Coats Land where the initial bedrock topography has subsided 245 below sea level, compared to the TOPO 14Ma case (Fig. S3a). The equilibrium ice volume at 280 ppm reduces to 42.8 m.s.l.e. Conversely, it increases to 28.4 m.s.l.e. (+5% w.r.t. REF), and is centered around a generally larger ice volume, in experiment TOPO 34Ma (Fig. 7a,b). 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 in Coats Land, the Amery embayment, and Marie Byrd Land (Fig. S3b). 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 older geological reconstructions of Antarctic paleotopography at the Eocene-Oligocene boundary of Wilson et al. (2012). At 280 ppm, the ice then grounds in the 255 Ross embayment (TOPO Wilson_mean, using the mean of the minimum and maximum reconstructions; Fig. (Fig. S3c)), or in both embayments (TOPO Wilson_max, using the maximum reconstruction; Fig. S3d). Due to the higher CO 2 levels at which deglaciation takes place, a relatively large ice volume persists throughout 40-kyr transient simulations (Fig. 7c,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, 260 is used as starting point (not shown).

Influence of ice shelf formation
We perform and compare two additional sensitivity experiments that are end-members with respect to basal melt rate settings, to assess the influence of ice shelf formation on AIS variability. In experiment BMB LGM (Table 1) (Fig. 8a). Ice shelf form along the margins at Coates Land and Oates Land (Fig. 9a). While in Coates 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 270 simulation, these ice shelves form after around 110 kyr, at ∼530 ppm (Fig. 8b). 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. 9b). The effect on grounded ice is particularly noticeable in West-Antarctica, and Wilkes Land, George V Land and MacRobertson Land, where the ice sheet thickens by kilometers.
In the deglaciation phase these ice shelves remain influential until 728 ppm in the equilibrium simulations (upper, descending  (Fig. 8c,d).

280
The CO 2 levels at which East-Antarctic glaciation is simulated, are known to be dependent on the climate model used to construct the climatic forcing (Gasson et al., 2014). This could (partly) explain why 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 greenhouse gasses other than CO 2 are varied, probably biasing the CO 2 concentrations high. In contrast, our climate input comes from the same model as Halberstadt et al. (2021), who obtain relatively large EAIS volumes at CO 2 levels up to 1140 ppm. In absolute sense, this CO 2 difference can simply be caused by the value of the ablation factor that we use in our surface mass balance calculation (C 3 in Eq. 4), as demonstrated by our SMB sensitivity experiment (Fig. 3). 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). Furthermore, they additionally use a regional climate model for downscaling 290 and refrain from using bias corrections for the temperature and precipitation forcing. Nevertheless, these factors likely do not explain the narrower CO 2 range between inception of the East-and West-Antarctic ice sheet in our simulations compared to Halberstadt et al. (2021). At the root of this discrepancy can be the different schemes we use to calculate precipitation and ablation, which may be more sensitive to temperature changes. 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 295 volume variability than Stap et al. (2019), when using a similar glacial index method (experiment FEEDB). Conducting an intermodel comparison of early AIS dynamics would be advisable for a more comprehensive investigation of the different results yielded by ice sheet models. We nonetheless corroborate earlier findings that disequilibrium between the AIS and the forcing climate causes smaller transient ice volume variability then 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 300 reference case. We have shown that these feedbacks tend to increase hysteresis in the CO 2 -equilibrium ice volume relation, and decrease the amplitude of AIS variability under equal 40-kyr forcing CO 2 cycles by 21%. This implies that correspondingly, the contribution of AIS volume changes to benthic δ 18 O fluctuations gets smaller, although the strength of the impact will depend on the isotopic composition of Antarctic snow as well (Gasson et al., 2016;Rohling et al., 2021). 305 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. 6). This is demonstrated by the lower equilibrium ice volumes at equal CO 2 levels in our experiment TOPO 14Ma, which has a bedrock topography pertaining to the mid-Miocene (14 Myr ago), compared to REF, in which 310 the topography is from the Miocene-Oligocene boundary (23-24 Myr ago). We also corroborate the finding of Paxman et al.
(2020), that from 34 to 24 million years ago bedrock changes primarily concern the EAIS ( Fig. 7), 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 14Ma than in REF, due to impeded WAIS growth. Hence, CO 2 minima must have dropped increasingly further over the course of the early-315 and mid-Miocene to obtain AIS variability with similar amplitude.
We simulate the inception of ice shelves only at relatively low CO 2 levels. Ice shelves first appear along the coasts of Oates Land and Coates Land when the CO 2 concentration drops below ∼530 ppm, but the effect of ice shelves on grounded ice volume only becomes significant below ∼360 ppm CO 2 when the EAIS -with a volume of almost 30 m.s.l.e. -is already well-established. This transition from a land-based to a marine ice sheet 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 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). In the 325 Ross and Filchner-Ronne embayments, however, 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 (large-scale) Antarctic ice shelf inception. On the other hand, the fact that the Ross and Filchner-Ronne ice shelf do grow when we impose steady present-day climatic 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 330 oceanic conditions, suitable atmospheric conditions are very important for ice shelf formation as well. Possibly, WAIS growth is therefore 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 in the forcing climate simulation. In the future, the forcing steady-state climate simulations used in the matrix interpolation could 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 335 and transient simulations, so that the transient ice sheet will not grow beyond its largest size in the forcing equilibrium results.
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 are reached. Only then does the AIS retreat into the Wilkes Subglacial Basin, suggesting periodic AIS waxing and waning in this region has to be governed by substantial CO 2 fluctuations. Growth 340 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 floatation) in steady-state PD simulations with a 400-m/yr 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. S4a,c). When the simulation is started from a PD topography including grounded and floating ice, the final 345 volume is 48.4 m.s.l.e., which is only 14% smaller than the normal PD run (Fig. S4b,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 (Sun et al., 2020). Our simulated 8.2 m.s.l.e. loss of ice volume above floatation after ice sheet removal under PD forcing conditions (starting with ice) is on the large side of the wide range of model responses (∼1 350 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, that 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 355 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 e.g. the Wilkes Subglacial Basin, hence reducing hysteresis in the CO 2 -ice volume relation.
Here, we have conducted idealised transient simulations, addressing the effects of the ocean, atmosphere and solid earth on Miocene AIS evolution. Moving forward towards realistic transient simulations of AIS evolution, our results of experiment 360 INSOL indicate that orbital forcing variability will need to be accounted for in the model set-up (Fig. 3), 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., 2018b;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 365 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).

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 370 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). We find competing influences of the atmosphere, ocean and solid earth on transient

Miocene AIS variability:
-Compared to a simpler glacial index method, as used e.g. in Stap et al. (2019), our matrix interpolation method incorporates two important ice-sheet-atmosphere feedbacks: the positive albedo-temperature feedback, and the negative 375 ice-volume-precipitation feedback. The former tends to increase hysteresis in the CO 2 -V relation and decrease the amplitude of transient AIS variability, while the latter has the opposite effect (Fig. 5). The albedo-temperature feedback is the strongest of the two, as in the 40-kyr transient reference simulation, the ice volume variability amplitude is 21% smaller than when a glacial index method is used.
-Differences between the geological reconstructions of Antarctic bedrock topography at the Oligocene-Miocene boundary 380 (23-24 Myr ago) and the mid-Miocene (14 Myr ago) are mainly located in West-Antartica. The subsidence of land below sea level during the early-and mid-Miocene impedes WAIS build-up. Consequently, it diminishes hysteresis in the CO 2ice volume above floatation (V) relation and reduces 40-kyr AIS variability by 10% in strength (Figs. 5,6). Hence, under equal atmospheric and oceanic forcing conditions, the amplitude of transient AIS variability becomes progressively smaller during the early-and mid-Miocene, due to glacial erosion.
-The importance of ice shelf formation is demonstrated by comparing simulations in which LGM basal melt rates are applied, to simulations in which ice shelf growth is prevented (Figs. 8,9). Ice shelves start to form around CO 2 levels of ∼530 ppm at the margins of Coates Land and Oates Land and at ∼360 ppm approximately simultaneously around all other margins. They continue to affect grounded ice volume until the forcing CO 2 concentration is increased to ∼728 ppm in our 400-kyr transient simulation. The formation of ice shelves leads to a thickening of the hinterlying 390 grounded ice, particularly of the WAIS, and consequently increases hysteresis in the CO 2 -V relation. 40-kyr transient AIS variability is amplified by 19% in case of LGM basal melt rates -and 5% in the reference case -compared to the case without ice shelves.
In our simulations, the ice-sheet-atmosphere feedbacks and the evolution of bedrock topography, that is due mostly to glacial erosion, hence tend to reduce the contribution of AIS changes to orbital timescale benthic δ 18 O variability during the Miocene.  Data availability. All relevant data related to this article will be made openly accessible from the PANGAEA database (https://pangaea.de/).

400
Appendix A: Ice-sheet model settings In IMAU-ICE, Glenn's flow law exponent is set to n = 3, and the SIA-and SSA-based velocities are enhanced by factors 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) (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 405 geothermal heat flux is taken from Shapiro and Ritzwoller (2004). A basic elastic lithosphere-relaxing 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 x 10 25 kg m 2 s −2 , and a relaxation time of 3000 yr. The sea level is kept constant at present-day level (0 m).
Appendix B: Climate matrix forcing 410 Using a so-called matrix method (Pollard, 2010;Pollard et al., 2013), a time-continuous climate forcing for the ice-sheet model can be calculated. This is achieved by interpolating between the GENESIS time-slice climate simulations, based on the varying 13 https://doi.org/10.5194/tc-2021-309 Preprint. Discussion started: 1 November 2021 c Author(s) 2021. CC BY 4.0 License. CO 2 levels and ice coverage. The matrix interpolation routine used in this study is described in detail in Berends et al. (2018), and 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(Berends et al., , 2021.

B1 Surface air temperature
For the applied surface air temperature, first two weighing factors for CO 2 (w CO2 ) and ice-sheet coverage (w ice ) are calculated.
where CO 2 is the applied CO 2 level, and CO 2,cold and CO 2,warm the low (280 ppm) and high (840 ppm) CO 2 levels. Weight factor w ice 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 calculate w ice , first the absorbed insolation I abs is calculated from the incoming insolation (Q T OA ) and internally calculated surface albedo (α) for the modeled field (I abs,mod ), and the fields interpolated between (I abs,cold , I abs,warm ): Then, the weighing field for insolation (w ins ) is calculated: w ins = I abs,mod − I abs,cold I abs,warm − I abs,cold Finally, the weighing field for ice-sheet coverage is calculated from the smoothed (using a Gaussian smoothing filter with a radius of 200 km) and average (over the ice-sheet model grid) w ins , to account for local as well as regional effects: For each grid cell in the domain the final weighing factor w tot is calculated as: CO2 is the ratio between the effects of CO 2 and ice-sheet coverage, predetermined (in our reference simulation) on the basis of GENESIS simulations 3fumebi (840 ppm CO 2 , full ice sheet) and 1nomebi (280 ppm CO 2 , no ice sheet): where T is the average temperature over the Antarctic domain. In the reference case CO2 = 0.62. The weighing factors are clipped at -0.25 and 1.25, to avoid spurious extreme temperatures and precipitation. and warm (T warm , h warm ) forcing: and likewise The applied 2-m air temperature field T applied is obtained by first correcting for the height differences between the modeled 445 (h mod ) and reference (h ref ) surface height using a uniform lapse rate λ = 8 K/km. Finally, a spatially-varying bias correction based on the difference between the modeled preindustrial temperature (T P I ) and the 1957-2002 ERA40 reanalysis data (Uppala et al., 2005) (T ERA40 ) is applied: Precipitation is weighed based on the total ice volumes modeled (V mod ) and in the cold (V cold ) and warm climates (V warm ): A reference precipitation field is obtained, by interpolating between the warm (P warm ) and cold (P cold ) states: Similarly as for the temperature, a bias correction is applied: The term 0.008 is added to avoid problems caused by division by very small numbers. And finally, the precipitation field is downscaled to the finer ice-sheet model grid, using a Clausius-Clapeyron relation (Lorius et al., 1985;Huybrechts, 2002;De Boer et al., 2013):