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

Benthic δ¹⁸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₂ 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₂ levels (between 280 and 840 ppm), as well as varying ice-sheet configurations (between no ice and a large East Antarctic Ice Sheet). We furthermore implement recent reconstructions of Miocene Antarctic bedrock topography. We find that the positive albedo–temperature feedback, partly compensated for by a negative feedback between ice volume and precipitation, increases hysteresis in the relation between CO₂ and ice volume. Together, these ice-sheet–atmosphere interactions decrease the amplitude of Miocene AIS variability in idealised transient simulations. Forced by quasi-orbital 40 kyr forcing CO₂ cycles, the ice volume variability reduces by 21 % when ice-sheet–atmosphere interactions are included compared to when forcing variability is only based on CO₂ changes. Thereby, these interactions also diminish the contribution of AIS variability to benthic δ¹⁸O fluctuations. Evolving bedrock topography during the early and mid-Miocene also reduces ice volume variability by 10 % under equal 40 kyr cycles of atmosphere and ocean forcing.

1 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₂ 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 δ¹⁸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., 2011, 2017; Holbourn et al., 2013; Levy et al., 2019). The δ¹⁸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 debate. 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 (Sangiorgi et al., 2018) and the Ross Sea sector (Hauptvogel and Passchier, 2012; Levy et al., 2016; Pérez et al., 2021). Some studies, on the other hand, advocate a more stable small AIS based on ice-proximal sea-surface temperature records from dinoflagellate cyst assemblages (Bijl et al., 2018) and TEX₈₆ data (Hartman et al., 2018). In addition, a recent analysis of clumped isotope and $\mathrm{Mg}/\mathrm{Ca}$ records has found high bottom-water temperatures, which combined with the benthic δ¹⁸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 δ¹⁸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 δ¹⁸O changes. Halberstadt et al. (2021) used the same model set-up to determine the influence of CO₂ 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, 2010; Stap et al., 2017). Recently, Stap et al. (2019) showed that orbital-timescale CO₂ 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₂ 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 (surface-height–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., 2016, 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₂ 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₂ 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₂ level constant, and idealised transient experiments, in which we vary the CO₂ 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₂ 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.

2 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).

2.1 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 (Winkelmann et al., 2011). 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.

2.1.1 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):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}M(x,y)& =\left(z_{\mathrm{deep}}\right(x,y)-\mathrm{1})\left[\right(\mathrm{1}-z_{\mathrm{expo}}(x,y)\left)S_{\mathrm{m}}\right(x,y)\\ & +z_{\mathrm{expo}}(x,y)M_{\mathrm{expo}}]-(z_{\mathrm{deep}}(x,y)M_{\mathrm{deep}}),\end{array}\end{array}$$

with

$$\begin{array}{}\text{(2)}& S_{\mathrm{m}}(x,y)={\mathit{\rho }}_{\mathrm{w}}{c}_p{\mathit{\gamma }}_{\mathrm{T}}{\displaystyle \frac{T_{\mathrm{w}}-T_{\mathrm{freeze}}(x,y)}{H_{\mathrm{fus}}{\mathit{\rho }}_{\mathrm{i}}}}.\end{array}$$

Weighing factor z_deep grows with water depths (h_w(x,y)) greater than 1800 m:

$$\begin{array}{}\text{(3)}& z_{\mathrm{deep}}(x,y)={\displaystyle \frac{h_{\mathrm{w}}(x,y)-\mathrm{1800}\mathrm{m}}{\mathrm{200}\mathrm{m}}},\end{array}$$

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)):

$$\begin{array}{}\text{(4)}& z_{\mathrm{expo}}(x,y)={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{sub}}(x,y)-\mathrm{80}{}^{\circ }}{\mathrm{30}{}^{\circ }}}{e}^{-\mathrm{0.4}\times r_{\mathrm{open}}(x,y)}.\end{array}$$

Both factors are clipped at 0 and 1. The melt rate for the exposed shelves (M_expo) varies linearly between 3 m yr⁻¹ for an applied CO₂ concentration of 280 ppm to 6 m yr⁻¹ for CO₂ 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(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 (Gasson et al., 2016; 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₂ 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⁻¹, M_deep=2 m yr⁻¹, and $T_{\mathrm{w}}=-\mathrm{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). 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.

(Hochmuth et al., 2020a; Paxman et al., 2019)(Hochmuth et al., 2020a; Paxman et al., 2019)(Morlighem et al., 2020)(Wilson et al., 2012)(Wilson et al., 2012)

Table 1 Name and description of the experiments carried out. Each experiment comprises 15 simulations, as explained in the main text and listed in Table 2.

2.1.2 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:

$$\begin{array}{}\text{(5)}& \mathrm{Acc}(x,y)=P_{\mathrm{applied}}(x,y)\times f_{\mathrm{snow}}(x,y),\end{array}$$

with f_snow(x,y) being the temperature-dependent snow fraction (Ohmura, 1999). Ablation (Abl(x,y)) follows from an insolation–temperature melt parameterisation (Bintanja et al., 2002):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathrm{Abl}(x,y)& =c_{\mathrm{1}}\left(T_{\mathrm{applied}}\right(x,y)-\mathrm{273.15}\mathrm{K})\\ & +c_{\mathrm{2}}\left(Q_{\mathrm{TOA}}\right(x,y)\times (\mathrm{1}-\mathit{\alpha }(x,y)\left)\right)+c_{\mathrm{3}}.\end{array}\end{array}$$

The parameters are set to c₁=0.788 m yr⁻¹ K⁻¹, c₂=0.004 m³ J⁻¹, and $c_{\mathrm{3}}=-\mathrm{0.50}$ m yr⁻¹. The internally calculated monthly albedo (α(x,y)) is derived as follows. First, a background albedo (α_bg(x,y)) is determined: 0.1 for ocean-covered 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:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathit{\alpha }(x,y)& ={\mathit{\alpha }}_{\mathrm{snow}}-({\mathit{\alpha }}_{\mathrm{snow}}-{\mathit{\alpha }}_{\mathrm{bg}}(x,y\left)\right)\\ & \times \exp [-\mathrm{15}\times z_{\mathrm{firn}}(x,y\left)\right]-\mathrm{0.015}\\ & \times {\mathrm{Abl}}_{\mathrm{yr}-\mathrm{1}}(x,y),\end{array}\end{array}$$

where Abl${}_{\mathrm{yr}-\mathrm{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₂ 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).

2.2 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₂ level of 280 ppm (1 × pre-industrial (PI)) and an AIS configuration with a full East Antarctic Ice Sheet (18.975×10⁶ km³ volume, area shown in Fig. S2). Simulation 3nomebi is forced by a CO₂ 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.

2.3 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₂ and on the internally modelled ice sheets. The higher (lower) the CO₂ 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.

2.4 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).

Table 2 Description of the 15 simulations spanning each experiment.

The steady-state simulations are run for 150 kyr with a constant forcing CO₂ level – enough time for the ice sheet to equilibrate with the forcing climate. Mind that the CO₂ 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₂ 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₂–ice-volume relation (blue connected diamonds in figures, e.g. Fig. 1b). Thereafter, the remaining five steady-state simulations are continued from results of the 280 ppm IMAU-ICE run, with CO₂ levels of 392, 504, 616, 728, and 840 ppm. These simulations constitute the upper, descending branch of the hysteresis curve in the equilibrium CO₂–ice-volume relation (red connected diamonds in figures, e.g. Fig. 1b).

Figure 1 (a) Evolution of forcing CO₂ levels (pink) and ice volume above flotation (in metres sea level equivalent; black) over time (relative to the length of the simulation) for the transient 100 kyr (dotted) and 400 kyr (long-dashed) REF simulations. (b) Relation between CO₂ and equilibrium ice volume (ascending branch, blue; descending branch, red) and transient ice volume (black as in a). (c, d) Same for the transient 40 kyr simulations starting from warm (CO₂ pink, ice volume solid black) and cold (CO₂ purple, ice volume dashed black) conditions. Arrows in (b) and (d) indicate the progression direction of ice volume. In (d), the transient ice volume evolutions of both simulations partly overlap.

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₂ over time has a triangular shape (e.g. pink line in Fig. 1a). The CO₂ is first linearly reduced from 840 to 280 ppm and then increased back to 840 ppm. In the three simulations, this triangular CO₂ 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₂ evolution. That means the CO₂ is first linearly increased from 280 to 840 ppm and then decreased back to 280 ppm (e.g. purple line in Fig. 1c).

3 Results

3.1 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–2002) precipitation and temperature forcing (Uppala et al., 2005). Constant present-day-like melt parameter settings are used to calculate the basal mass balance: M_expo=3 m yr⁻¹, M_deep=5 m yr⁻¹, and $T_{\mathrm{w}}=-\mathrm{1.7}$ ^∘C. Depending on whether the run is started from BedMachine PD ice and topography (Morlighem et al., 2020) 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¹⁵ m³) or 55.0 m s.l.e. (22.4×10¹⁵ m³) 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₂ 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₂ levels are shown in Fig. S4.

Figure 2 (a) Simulated equilibrated grounded ice thickness from experiment REF with 840 ppm and (b) 280 ppm CO₂ forcing, starting without ice. Cyan areas indicate ice shelf extent and grey areas ice-free land.

We show the reference CO₂ and equilibrium ice volume (above flotation) (V_eq) relation in Fig. 1b (Pollard and DeConto, 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₂ minimum (Fig. 1a). The peak is attained relatively earlier in time (Fig. 1a), i.e. at a lower CO₂ 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₂ level (252 ppm difference) than the minimum ice volume after the maximum CO₂ 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⁻¹, speeding up to a maximum of 3 m s.l.e. kyr⁻¹ 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₂–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⁻¹).

In earlier research using IMAU-ICE, different values for the SMB parameter c₃ in Eq. (6) have been used (e.g. De Boer et al., 2013, 2014; Berends et al., 2018). Our reference setting C₃=0.50 m yr⁻¹ leads to satisfactory PD results, but using C₃=0.25 m yr⁻¹ 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₂ 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₂ at which these ice-sheet sizes are simulated but does not (qualitatively) affect the shape of the CO₂–ice-volume relation.

3.2 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₂ and ice-sheet size. Here, we make a comparison to results 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₂ level. In experiment NOFEEDB (Table 1), we implement $w_{\mathrm{tot}}=\frac{{\mathrm{CO}}_{\mathrm{2}}-{\mathrm{CO}}_{\mathrm{2},\mathrm{cold}}}{{\mathrm{CO}}_{\mathrm{2},\mathrm{warm}}-{\mathrm{CO}}_{\mathrm{2},\mathrm{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₂–V_eq relation is smaller (Fig. 3b, solid and dashed red and blue lines). The temperature declines more quickly towards lower CO₂ 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₂ levels of glaciation and deglaciation is smaller in the NOFEEDB case. That reduces the hysteresis in the CO₂–V_eq relation.

Figure 3 (a) Evolution of forcing CO₂ levels (pink) and ice volume over relative time for the transient 40 kyr REF (black) and NOFEEDB (orange) simulations. (b) Relation between CO₂ and equilibrium ice volume (REF, solid red and blue; NOFEEDB, dashed red and blue) and transient ice volume (as in a). Only the equilibrium cycle (final 40 kyr) is shown. (c, d) Same for the NOFEEDB-P (dotted red, blue, and orange) and NOFEEDB-T (long-dashed red, blue, and orange) equilibrium and transient 40 kyr simulations. The amplitudes of transient ice volume variability (in m s.l.e.) in the equilibrium cycle are indicated in panels (a) and (c).

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₂ cycle) and 47.4 m s.l.e. (100 kyr). Forced by 40 kyr CO₂ 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 albedo–temperature 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₂–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-sheet–atmosphere feedbacks, the albedo–temperature feedback is dominant in our experiments.

3.3 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₂ 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).

Figure 4 (a) Simulated equilibrated grounded ice thickness from experiment TOPO 14 Ma with 280 ppm CO₂ forcing. Cyan areas indicate ice shelf extent and grey areas ice-free land. (b) Difference in ice thickness between TOPO 14 Ma and REF with 280 ppm CO₂ forcing (grounding line and continental edge as in a).

Figure 5 (a) Evolution of forcing CO₂ levels (pink) and ice volume over relative time for the transient 40 kyr REF (black) and TOPO 14 Ma (orange) simulations. (b) Relation between CO₂ and equilibrium ice volume (REF, solid red and blue; TOPO 14 Ma, dashed red and blue) and transient ice volume (as in a). Only the equilibrium cycle (final 40 kyr) is shown.

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). The equilibrium ice volume at 280 ppm reduces to 42.8 m s.l.e. (−12 % with respect to REF) and the amplitude of transient 40 kyr ice volume variability to 23.3 m s.l.e. (−14 % with respect to REF; Fig. S8a, b dotted lines). Conversely, it increases to 28.4 m s.l.e. (+5 % with respect to REF; Fig. 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₂ levels in the CO₂–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₂ 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).

3.4 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⁻¹, M_deep=2 m yr⁻¹, and $T_{\mathrm{w}}=-\mathrm{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⁻¹. 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₂–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).

Figure 6 (a) Evolution of forcing CO₂ levels (pink) and ice volume over relative time for the transient 400 kyr BMB no_shelves (solid orange) and BMB LGM (dashed green) simulations. (b) Relation between CO₂ and equilibrium ice volume (BMB no_shelves, solid red and blue; BMB LGM, dashed red and blue) and transient ice volume. (c, d) Same for the transient 40 kyr simulations. Only the equilibrium cycle (final 40 kyr) is shown.

4 Discussion

The CO₂ 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₂ 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₂ doubling) than GENESIS (2.9 K). Note that in neither study are greenhouse gases other than CO₂ varied, probably biassing the CO₂ 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₃ 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).

Figure 7 (a) Simulated equilibrated grounded ice thickness from experiment BMB LGM with 504 ppm and (b) 392 ppm CO₂ forcing, starting without ice. Cyan areas indicate ice shelf extent and grey areas ice-free land.

Here, we include ice-sheet–atmosphere feedbacks in our reference case. We have shown that these feedbacks tend to increase hysteresis in the CO₂ and equilibrium ice volume relation and decrease the amplitude of AIS variability under equal 40 kyr forcing CO₂ 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₂ levels and partly to decreased growth rates of the transient ice volume relative to the change in equilibrium ice volume. In short, the ice-sheet–atmosphere feedbacks cause a slower build-up of the AIS to smaller peak ice volumes. Correspondingly, the contribution of AIS volume changes to benthic δ¹⁸O fluctuations gets smaller. This calls for a larger contribution of deep-sea temperatures to the amplitude of Miocene benthic δ¹⁸O fluctuations. The actual strength of the impact will, however, depend on the isotopic composition of Antarctic snow as well (Gasson et al., 2016; 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₂ 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₂ cycles in TOPO 14 Ma than in REF due to impeded WAIS growth. Hence, CO₂ 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 δ¹⁸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 (large-scale) 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₂ 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₂ (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₂ 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., 2016, 2019). Below ∼ 360 ppm CO₂, 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₂ levels from cold conditions in our simulations, the ice shelves keep having an effect on grounded AIS ice volume until relatively high CO₂ 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₂ fluctuations. Growth of ice shelves hence increases hysteresis in the CO₂–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⁻¹ 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 (Sun et al., 2020). 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₂ level, of the AIS in, for example, the Wilkes Subglacial Basin, hence reducing hysteresis in the CO₂–ice-volume relation.

Our climate input comes from the same model (GENESIS) as Halberstadt et al. (2021), who obtain relatively large EAIS volumes at CO₂ 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 because 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₂ levels when we refrain from using this bias correction, although this is dependent on the ablation factor (c₃ in Eq. 6; not shown). All in all, even with the bias correction included, we need very large CO₂ 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 δ¹⁸O fluctuations during the Miocene (Bijl et al., 2018; 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 (Marschalek et al., 2021), in agreement with the modelling results of, for example, Gasson et al. (2016) and Halberstadt et al. (2021). The generally low benthic δ¹⁸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.

5 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₂ 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 ice-volume–precipitation feedback. The former tends to increase hysteresis in the CO₂ 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₂–V relation and reduces 40 kyr AIS variability by 10 % in amplitude (Figs. 4, 5). Hence, CO₂ 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₂ 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 and DeConto, 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⁻¹) (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 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⁻³, a lithospheric flexural rigidity of 1.0×10²⁵ kg m² s⁻², 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 ice-sheet 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₂ 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, 2021).

B1 Surface air temperature

For the applied surface air temperature, first we calculate a scalar weighing factor for CO₂ ($w_{{\mathrm{CO}}_{\mathrm{2}}}$) and a spatially varying one for ice-sheet coverage (w_ice(x,y)).

$$\begin{array}{}\text{(B1)}& w_{{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2}}-{\mathrm{CO}}_{\mathrm{2},\mathrm{cold}}}{{\mathrm{CO}}_{\mathrm{2},\mathrm{warm}}-{\mathrm{CO}}_{\mathrm{2},\mathrm{cold}}}},\end{array}$$

where CO₂ is the applied CO₂ level, and CO_2,cold and CO_2,warm are the low (280 ppm) and high (840 ppm) CO₂ 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).

$$\begin{array}{}\text{(B2)}& I_{\mathrm{abs}}(x,y)=(\mathrm{1}-\mathit{\alpha }(x,y\left)\right)\times Q_{\mathrm{TOA}}(x,y).\end{array}$$

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:

$$\begin{array}{}\text{(B3)}& w_{\mathrm{ins}}(x,y)={\displaystyle \frac{I_{\mathrm{abs},\mathrm{mod}}(x,y)-I_{\mathrm{abs},\mathrm{cold}}(x,y)}{I_{\mathrm{abs},\mathrm{warm}}(x,y)-I_{\mathrm{abs},\mathrm{cold}}(x,y)}}.\end{array}$$

Weighing factors $w_{{\mathrm{CO}}_{\mathrm{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:

$$\begin{array}{}\text{(B4)}& w_{\mathrm{ice}}(x,y)={\displaystyle \frac{\mathrm{1}}{\mathrm{7}}}{w}_{\mathrm{ins},\mathrm{smoothed}}(x,y)+{\displaystyle \frac{\mathrm{6}}{\mathrm{7}}}{w}_{\mathrm{ins},\mathrm{avg}}.\end{array}$$

The viability of the partitioning $w_{\mathrm{ins},\mathrm{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 δ¹⁸O. For each grid cell in the domain the final weighing factor w_tot(x,y) is calculated as follows:

$$\begin{array}{}\text{(B5)}& w_{\mathrm{tot}}(x,y)=\mathit{ϵ}\times w_{{\mathrm{CO}}_{\mathrm{2}}}+(\mathrm{1}-\mathit{ϵ})\times w_{\mathrm{ice}}(x,y),\end{array}$$

where ϵ is the ratio between the effects of CO₂ and ice-sheet coverage. To determine ϵ for our simulations, we use GENESIS simulations 3fumebi (840 ppm CO₂, large ice sheet) and 1nomebi (280 ppm CO₂, no ice sheet). These climate simulations isolate the separate effects of CO₂ and the ice sheet respectively. The (spatially averaged) relative influence of CO₂ on temperature is estimated as the isolated effect of CO₂ divided by the combined influence of CO₂ and ice-sheet coverage.

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}\mathit{ϵ}& ={\displaystyle \frac{{\stackrel{\mathrm{‾}}{T}}_{\mathrm{3}\mathrm{fumebi}}-{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}\mathrm{fumebi}}}{({\stackrel{\mathrm{‾}}{T}}_{\mathrm{3}\mathrm{fumebi}}-{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}\mathrm{fumebi}})+({\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}\mathrm{nomebi}}-{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}\mathrm{fumebi}})}}\\ & =\mathrm{0.62},\end{array}\end{array}$$

where $\stackrel{\mathrm{‾}}{T}$ is the average temperature over the Antarctic domain. Synergy between the effects of CO₂ and ice-sheet coverage is relatively small so that $\mathit{ϵ}\approx ({\stackrel{\mathrm{‾}}{T}}_{\mathrm{3}\mathrm{fumebi}}-{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}\mathrm{fumebi}})/({\stackrel{\mathrm{‾}}{T}}_{\mathrm{3}\mathrm{nomebi}}-{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}\mathrm{fumebi}})$.

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:

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}{T}_{\mathrm{ref}}(x,y)& =w_{\mathrm{tot}}(x,y)\times T_{\mathrm{warm}}(x,y)\\ & +(\mathrm{1}-w_{\mathrm{tot}}(x,y\left)\right)\times T_{\mathrm{cold}}(x,y),\end{array}\end{array}$$

and likewise

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}{h}_{\mathrm{ref}}(x,y)& =w_{\mathrm{tot}}(x,y)\times h_{\mathrm{warm}}(x,y)\\ & +(\mathrm{1}-w_{\mathrm{tot}}(x,y\left)\right)\times h_{\mathrm{cold}}(x,y).\end{array}\end{array}$$

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⁻¹. 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:

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}{T}_{\mathrm{applied}}(x,y)& =T_{\mathrm{ref}}(x,y)-\mathit{\lambda }\times \left(h_{\mathrm{mod}}\right(x,y)-h_{\mathrm{ref}}(x,y\left)\right)\\ & +\left(T_{\mathrm{ERA}\mathrm{40}}\right(x,y)-T_{\mathrm{PI}}(x,y\left)\right).\end{array}\end{array}$$

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):

$$\begin{array}{}\text{(B10)}& w_{\mathrm{tot}}=\mathrm{1}-{\displaystyle \frac{V_{\mathrm{mod}}-V_{\mathrm{warm}}}{V_{\mathrm{cold}}-V_{\mathrm{warm}}}}.\end{array}$$

A reference precipitation field is obtained by interpolating between the cold (P_warm) and warm (P_cold) forcing fields:

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}{P}_{\mathrm{ref}}(x,y)& =\exp \left[w_{\mathrm{tot}}\times \log \left[P_{\mathrm{warm}}\right(x,y\left)\right]\right.\\ & \left.+(\mathrm{1}-w_{\mathrm{tot}})\times \log \left[P_{\mathrm{cold}}\right(x,y\left)\right]\right].\end{array}\end{array}$$

Similarly as for the temperature, a bias correction is applied:

$$\begin{array}{}\text{(B12)}& P_{\mathrm{ref},\mathrm{corr}}(x,y)=P_{\mathrm{ref}}(x,y)\times {\displaystyle \frac{P_{\mathrm{ERA}\mathrm{40}}(x,y)+\mathrm{0.008}}{P_{\mathrm{PI}}(x,y)+\mathrm{0.008}}}.\end{array}$$

The term 0.008 is added to avoid problems caused by division by very small numbers. 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):

$$\begin{array}{}\text{(B13)}& \begin{array}{rl}& P_{\mathrm{applied}}(x,y)=P_{\mathrm{ref},\mathrm{corr}}(x,y)\times {\left({\displaystyle \frac{T_{\mathrm{I},\mathrm{ref}}(x,y)}{T_{\mathrm{I}}(x,y)}}\right)}^{\mathrm{2}}\\ & \times \exp \left[\mathrm{22.47}\times \left({\displaystyle \frac{\mathrm{273.16}K}{T_{\mathrm{I},\mathrm{ref}}(x,y)}}-{\displaystyle \frac{\mathrm{273.16}K}{T_{\mathrm{I}}(x,y)}}\right)\right],\end{array}\end{array}$$

with

$$\begin{array}{}\text{(B14)}& {\displaystyle }{T}_{\mathrm{I}}(x,y)=\mathrm{88.9}+\mathrm{0.67}\times T_{\mathrm{ref}}(x,y),\text{(B15)}& {\displaystyle }\begin{array}{rl}{T}_{\mathrm{I},\mathrm{ref}}(x,y)& =\mathrm{88.9}+\mathrm{0.67}\times \left(T_{\mathrm{ref}}\right(x,y)\\ & -\mathit{\lambda }\times \left(h_{\mathrm{mod}}\right(x,y)-h_{\mathrm{ref}}(x,y\left)\right)).\end{array}\end{array}$$