Anti‐Phased Miocene Ice Volume and CO2 Changes by Transient Antarctic Ice Sheet Variability

Geological evidence indicates large continental‐scale Antarctic ice volume variations during the early and mid‐Miocene. On million‐year timescales, these variations can largely be explained by equilibrium Antarctic ice sheet (AIS) simulations. In contrast, on shorter orbital timescales, the AIS needs not be in equilibrium with the forcing and ice volume variations may be substantially different. Here, we introduce a conceptual model, based on ice dynamical model results, to investigate the difference between transient variability and equilibrium differences of the Miocene AIS. In our model, an ice sheet will grow (shrink) by a specific rate when it is smaller (larger) than its equilibrium size. We show that phases of concurrent ice volume increase and rising CO2 levels are possible, even though the equilibrium ice volume decreases monotonically with CO2. When the AIS volume is out of equilibrium with the forcing climate, the ice sheet can still be adapting to a relatively large equilibrium size, although CO2 is rising after a phase of decrease. A delayed response of Antarctic ice volume to (covarying) solar insolation and CO2 concentrations can cause discrepancies between Miocene solar insolation and benthic δ18O variability. Increasing forcing frequency leads to a larger disequilibrium and consequently larger CO2‐ice volume phase differences. Furthermore, an amplified forcing amplitude causes larger amplitude ice volume variability, because the growth and decay rates depend on the forcing. It also leads to a reduced average ice volume, resulting from the growth rates generally being smaller than the decay rates.


Introduction
The early and mid-Miocene (23 to 14 Myr ago) is an essential period for studying the dynamics of the Antarctic ice sheet (AIS) in a warmer climate, since it was the last period in Earth's history when extensive continental-scale (tens of meters sea level equivalent) AIS changes occurred. Indeed, geological evidence from sea level records (John et al., 2011;Kominz et al., 2008;Miller et al., 2005) and ice-proximal ocean sediment cores (e.g., Levy et al., 2016) indicates that the AIS ranged from being almost vanished to its modern-day marine-terminating configuration during this time. Moreover, benthic 18 O records show large fluctuations on shorter, predominantly orbital timescales (e.g., Cramer et al., 2009;De Vleeschouwer et al., 2017;Holbourn et al., 2013;Liebrand et al., 2011Liebrand et al., , 2017Levy et al., 2019;Zachos et al., 2008). However, different factors cause significant uncertainty in the interpretation of these benthic 18 O fluctuations.
One factor was the focus of an earlier study (Stap et al., 2019). It concerns to what extent these benthic 18 O fluctuations are caused by deep-sea temperature changes and to what extent by ice volume changes  Zachos et al. (2008) interpolated to 100 years (red), over the period 21 to 20 Myr ago. Discrepancies between insolation and benthic δ 18 O, that is, times when both simultaneously increase or decrease, are marked yellow.
On the basis of equilibrium simulations, AIS modeling studies have demonstrated the importance of including Antarctic bedrock topography changes (Colleoni et al., 2018) and ice sheet-climate feedbacks, ice calving parameterizations, and changes in the ice sheet's oxygen isotopic composition (Gasson et al., 2016), for simulating strong ice sheet-induced benthic 18 O variations. However, we found that the slow response time of the AIS to orbital timescale forcing changes significantly reduces the magnitude of ice volume variability and hence its contribution to benthic 18 O fluctuations, in transient simulations (Stap et al., 2019).
Here, we consider another factor causing uncertainty, which concerns the phase relation between Miocene benthic 18 O fluctuations and solar insolation forcing changes caused by orbital variations (Laskar et al., 2004). We expect an anti-phase relation between the effect of orbital variations on Southern Hemisphere summer insolation and Antarctic ice volume. Ignoring the potentially considerable influence of deep-sea temperature changes (e.g., Hartman et al., 2018;Hutchinson et al., 2019;Stap et al., 2019), this would be reflected in an anti-phase relation between Southern Hemisphere insolation and benthic 18 O. In Figure  1, we show that for the stacked record of Zachos et al. (2008), this not always the case during the period 21 to 20 Myr ago. Although their duration depends on the age model used, such anti-phased fluctuations of benthic 18 O and Southern Hemisphere insolation are a robust feature in different global benthic 18 O records throughout the early and mid-Miocene (Cramer et al., 2009;De Vleeschouwer et al., 2017;Zachos et al., 2008). A fundamental consequence of equilibrium simulations is that simultaneous growth of ice volume and CO 2 level or solar insolation increase has to be caused by a larger increase in accumulation than in ablation due to the induced warmer temperatures (Oerlemans, 1991). When this is not the case, assuming equilibrium, discrepancies between solar insolation and benthic 18 O variations can only be explained as a consequence of counterbalancing insolation and CO 2 effects on temperature due to anti-phased changes.
As an alternative explanation, we here explore discrepancies between solar insolation and benthic 18 O variability arising from disequilibrium between transient covarying insolation and CO 2 changes and ice volume changes. We take this disequilibrium perspective to study the relation between atmospheric CO 2 and ice volume and consequently the transient variability of the Miocene AIS, by using a conceptual model.
In previous research, conceptual models have been used to investigate transient ice volume variability during the Pleistocene, which was mainly constituted by the waxing and waning of large ice sheets on the North American and Eurasian continents (e.g., Calder, 1974;Huybers, 2009;Huybers & Wunsch, 2005;Imbrie & Imbrie, 1980;Parrenin & Paillard, 2003). These models contained a difference in timescale between growth and decay of ice sheets, based on physical knowledge from early dynamical ice sheet modeling (Weertman, 1964). Growth and decay rates were fixed (Huybers & Wunsch, 2005;Imbrie & Imbrie, 1980), stochastic (Huybers & Wunsch, 2005), or depending, for example, on the existing ice volume (Huybers, 10.1029/2020PA003971 2009Parrenin & Paillard, 2003). Criteria based on "control" parameters, mostly orbital parameters regulating solar insolation, were used to determine whether an ice sheet is in growth or in decay phase. These criteria were typically tuned to match ice sheet volumes with benthic 18 O data.
Here, our aim is to analyze the difference in tendency between the transient ice volume variability and equilibrium differences, due to CO 2 -driven climate change. In principle, this is possible using any conceptual model with a stable equilibrium solution and a finite response scale. However, we introduce a similar model to address this difference more directly. Our model is based on the notion that an ice sheet will grow when it is smaller and shrink when it is larger, than its explicitly prescribed equilibrium size.
First, we analyze a simple linear relation between the control parameter (C) and equilibrium ice volume with specified growth and decay rates. We investigate the influence on transient ice volume variability of the relationship between the change rate of the equilibrium ice volume and these growth and decay rates.
In particular, we focus on the possibility of ice sheet growth during times of rising C, even though the equilibrium ice volume decreases monotonically with C. Furthermore, we analyze the sensitivity of ice volume variability to the amplitude and frequency of the forcing variability. Thereafter, we apply our model to the Miocene AIS, using an input CO 2 -equilibrium ice volume relation derived from results of the 3D thermodynamical ice sheet model-the Parallel Ice Sheet Model (PISM) (Stap et al., 2019). By using these results, we circumvent the aforementioned problem that the precise influence of ice volume and deep-sea temperature changes on Miocene 18 O variations is still highly uncertain. We match our conceptual growth and decay rates to transient PISM results. Based on these rates, we explain the occurrence of phases of concurrent AIS growth (decay) and CO 2 level rise (fall), as well as the dependence of the duration of these phases on the forcing frequency. Furthermore, we predict the influence of reducing the amplitude of CO 2 variability on mean ice volume and ice volume variability. Lastly, we discuss the implications of our results, by analyzing simulations over the period 21 to 20 Myr ago forced by CO 2 variability based on Southern Hemisphere summer insolation. We show how, in absence of other effects, the delayed and nonlinear response of Antarctic ice volume to (covarying) solar insolation and CO 2 concentrations alone can cause discrepancies between solar insolation and benthic 18 O variability during the Miocene.

Model and Methods
In this study, we employ a conceptual model of transient ice volume variability. It is based on the notion that at any level of a control parameter, which can for instance be CO 2 or temperature, an ice sheet will grow or shrink toward an equilibrium state. In our model, the control parameter-equilibrium ice volume (C-V eq ) relation of the ice sheet, as well as the growth and decay rates, have to be prescribed. The C-V eq relation can contain hysteresis, so that it has a top branch and a bottom branch. If the ice volume (V) is smaller than the specified bottom branch equilibrium solution (V eq_bot ) for the applied value of the control parameter (C) at a time step, the ice volume increases by the growth rate (g). Oppositely, it decreases by the decay rate (d) when it is larger than the top branch equilibrium solution (V eq_top ). The growth and decay rates can be constant but may also depend on the control parameter, the ice volume, and the difference between the ice volume and the equilibrium ice volume. If V is in between the top and bottom branch solutions, it is not clear what would happen to an ice sheet in reality. In our model, by default, the volume is kept unchanged in this case, but we discuss the influence of this choice on our results.

The Difference Between Transient Ice Volume Variability and Equilibrium Differences
To explore the difference between transient ice volume variability and equilibrium differences, we first prescribe a simple linear C-V eq relation without any hysteresis:  (Table 2), used as input for the conceptual model. The blue shaded area marks where the ice sheet is smaller than equilibrium and the red shaded area where it is larger. (c) Ice thickness from steady-state simulations of the Miocene Antarctic ice sheet using PISM (Stap et al., 2019), started with a present-day topography (Bedmap2) after deglaciation and isostatic rebound, for the lowest CO 2 level and (d) for the highest CO 2 level. Intermediate level results are not shown here. The gray line indicates the edge of the grounded ice sheet, the black line shows the edge of the continent or ice shelf, and gray areas indicate land that is not glaciated.
Note. Five forcing cycles are performed per experiment.
are unitless. We linearly change the control parameter (C) from 1 to 0 to 1 with a period of 200 time steps and run the model for 1,000 time steps (Table 1; Figure 3a, gray). Initial ice volume is zero.
If we set both the growth and the decay rate equal to the equilibrium change rate (g 0.01 d 0.01 ), the ice sheets remain in equilibrium throughout the run and the ice volume simply follows C (not shown). When the growth and decay rates are reduced (g 0.004 d 0.004 ), however, the ice sheet is smaller than equilibrium at the lowest level of C ( Figure 3, green solid lines). It continues to grow when C is increased again, to the point where equilibrium is reached. Likewise, the ice sheet continues to shrink after a reversal from an increasing to a decreasing C. Hence, phase differences between ice volume and the forcing are introduced. Furthermore, the amplitude of ice volume variability increases with decreasing forcing frequency (g 0.004 d 0.004 _slow; Figures 3c and 3d, green dashed lines). This can be attributed to the fact that the ice sheet has a longer time to build up and decay, which means a larger ice body can wax and wane. In case of lower forcing frequency, the ice sheet remains closer to its equilibrium size. The ice sheet then attains its equilibrium size relatively faster after a reversal of the forcing change. Therefore, the periods of concurrent ice volume increase and CO 2 rise and ice volume decrease and CO 2 decline are relatively shorter.
In the equilibrium cycle of simulation g 0.004 d 0.004 , the growth and decay phases of ice volume take an equal amount of time. When we implement a growth rate that is smaller than the decay rate (g 0.002 d 0.004 ), however, the growth phase lasts longer than the decay phase in the equilibrium cycle, even though the forcing cycle is still symmetrical (Figures 3e and 3f, black lines). This constitutes a sawtooth pattern in the evolution of ice volume over time. The duration of the phases is determined by the ratio between the growth and decay rates and the amplitude of ice volume variability by their absolute magnitudes. Decreasing the growth rate prolongs the growth phase and reduces ice volume variability. Increasing the decay rate also prolongs the growth phase but increases ice volume variability.
A 50% reduction of forcing variability amplitude (g 0.002 d 0.004 _half) leads to an increase of the mean ice volume but does not affect the variability (Figures 4a and 4b). When the growth and decay rates depend on the difference between the ice volume and the equilibrium ice volume (g nc d nc and g nc d nc _half), however, we obtain reduced amplitude ice volume variability with smaller forcing amplitude (Figures 4c and 4d). A reduction of the forcing amplitude leads to smaller differences between the ice volume and the equilibrium ice volume. This now causes smaller growth and decay rates and consequently smaller amplitude ice volume variability.
In the analysis so far, we have considered a simple linear C-V eq relation without any hysteresis. In Appendix A, we show that most of our results remain qualitatively unchanged when we use a C-V eq relation with appreciable hysteresis. This holds true for different ways of handling the situation when an ice sheet has a volume in between the bottom and top branch equilibrium ice volume (i.e., letting the ice volume grow, keeping it unchanged, and letting it decay). It is important to note, however, that when the ice sheet is kept unchanged between the bottom and top branches (our default option), a reduction of the forcing variability amplitude leads to smaller ice volume variability even in the case of constant growth and decay rates for the C-V eq relation with hysteresis. Note. In this study, these ice volumes are used to construct an CO 2 -index-V eq relation serving as input for the conceptual model ( Figure 2b).

CO 2 -Index-V eq Relation
To apply our conceptual model to the Miocene AIS, we prescribe a C-V eq (in this case CO 2 -index-V eq ) relation based on ice dynamical model results we presented in an earlier study (Stap et al., 2019) (Figures  2c and 2d). These were obtained by running the 3D thermodynamical PISM version 0.7.3 (Bueler & Brown, 2009;Winkelmann et al., 2011) into equilibrium using Miocene forcing conditions (precipitation, near-surface atmospheric, and 3D ocean temperatures). PISM uses a combination of the shallow ice approximation (SIA) and shallow shelf approximation (SSA) to solve the grounded ice velocities and only the SSA for the floating ice velocities. The surface mass balance is calculated using a positive degree day scheme and subshelf melt using a quadratic temperature relation. The forcing conditions were taken from Stärz et al. (2017), who simulated the early and mid-Miocene climate using the coupled atmosphere-ocean general circulation model (GCM) COSMOS (Jungclaus et al., 2006) under three CO 2 levels (278, 450, and 600 ppm). As detailed in Stap et al. (2019), we refer to the CO 2 levels in a relative (low, medium, and high) rather than absolute sense. While retaining the shape, this serves to eliminate model and model setup dependency of the absolute values of the CO 2 -ice volume relation from our results. Anomalies of the COSMOS simulations with respect to a preindustrial control run were added to a base climate provided by the 1970-2000 average of ERA-40 reanalysis data (Uppala et al., 2005) and the 1955-2006 average of the World Ocean Database 2009 data set (Locarnini et al., 2010). PISM simulations were conducted using climate forcing with low (CO 2 -index i = 0), medium (i = 0.5342), and high (i = 1) CO 2 levels. Additional simulations were performed using forcing conditions in between the low and medium (i = 0.2671) and medium and high (i = 0.7671) CO 2 levels, by taking the average of the two enclosing climates. The simulations were started without any ice, and with Antarctic bedrock topography taken from the Bedmap2 reconstruction (Fretwell et al., 2013), isostatically rebounded after removal of the ice. We supplement these simulations here by previously unpublished PISM runs under the same climate forcing conditions, but started from ice and Antarctic bedrock topography conditions taken from the results of the low-CO 2 PISM simulation, to complete the hysteresis loop (Table 2). In order to construct a continuous CO 2 -index-V eq relation, we assume that the equilibrium ice volume varies linearly between these equilibrated states (Figure 2b).
The growth and decay rates are dependent on the difference between the ice volume and equilibrium, to simulate the observed increase of the rates when this difference become larger. For the decay rate, we find that this dependence becomes less strong at larger differences in the PISM results, as reflected in the square root dependence. Furthermore, a cut-off factor of value 9 is invoked, because otherwise, the decay rate would increase too quickly after a reversal from a growing to a diminishing ice volume. We use the default option of keeping the ice volume unchanged in between the bottom and top branches of the CO 2 -index-V eq relation. This simplification does not impede correspondence between the transient conceptual model simulations and the PISM simulations, because this domain is rather small in the equilibrium PISM simulations (see also section 4).

Difference Between Transient and Equilibrium Ice Volume Variability
Our conceptualization explains the main differences between the transient and equilibrium PISM results of Stap et al. (2019) (Figure 5). Using the PISM-derived CO 2 -index-V eq relation with these growth and decay rates, we perform experiment 100kyr (Table 3) in which we linearly change the CO 2 -index (i) from 1 to 0 to 1 five successive times. Starting from i = 1 and zero ice volume, we apply a forcing period of 100 kyr for one CO 2 cycle, so 500 kyr for the complete simulation ( Figures  6a and 6b). Because the growth and decay rates are smaller than the equilibrium change rate, we obtain phases of concurrent CO 2 increase (decrease) ice sheet growth (decay), even though the equilibrium ice volume decreases monotonically with CO 2 . The onset of ice volume decay trails CO 2 increase by 24.3 kyr or 24% of the forcing period in the equilibrium cycle. Ice volume increase follows 11.0 kyr (11%) after the forcing starts to decrease. The growth phase thus lasts 14.3 kyr longer than the decay phase, as a consequence of the growth rate generally being smaller than the decay rate.
We additionally perform experiments with different forcing periods: 40kyr (Figures 6c and 6d, orange lines) and 400kyr (Figures 6c and 6d, green lines), which show that a lower forcing frequency leads to larger amplitude ice volume variability. The difference between the maximum and minimum ice volumes in the equilibrium cycle increases from 3.2 × 10 6 km 3 in experiment 40kyr, to 6.2 × 10 6 km 3 in 100kyr, and 10.1 × 10 6 km 3 in 400kyr. Meanwhile, the relative phase differences between ice volume and CO 2 are larger in 40kyr: 16% (6.2 kyr) for the onset of the growth phase and 28% (11.1 kyr) for the decay phase. Conversely, with 6% (25.6 kyr) and 20% (78.6 kyr), respectively, the growth and decay phase differences are smaller in 400kyr.

Response to Smaller Forcing Amplitude
Our conceptualization also predicts transient PISM behavior. We here investigate the sensitivity of ice volume to a 50% reduction of the forcing amplitude (Figures 6e and 6f). In experiment 100kyr_half, we cycle the CO 2 -index between 0.25 and 0.75. Because the growth rates are generally smaller than the decay rates, this forcing amplitude reduction leads to a slight increase of the equilibrium-cycle mean ice volume from 4.9 to 5.1 × 10 6 km 3 . Meanwhile, the difference between the maximum and minimum ice volumes now reduces by 52% to 3.0 × 10 6 km 3 , most probably as a consequence of the growth and decay rates being dependent on the difference between the ice volume and the equilibrium ice volume. Finally, the phase delays for the growth phase decrease to 10.6 kyr, while it increases to 27.8 kyr for the decay phase. Qualitatively similar results are obtained in the case of 40-kyr and 400-kyr forcing cycles (not shown).

Response to (Covarying) Miocene Solar Insolation and CO 2 Concentrations
We perform simulations using our conceptual model with the PISM-derived Miocene AIS CO 2 -index-V eq relation and CO 2 -index forcing based on solar insolation variability. We force the model over the arbitrarily chosen 1.5 million-year time frame from 21.5 to 20 Myr ago and use only the final million years simulated (21 to 20 Myr ago) for our analysis. We vary the CO 2 -index forcing, in sync with the 80 • S January mean insolation curve (Laskar et al., 2004) (Figure 1, blue), between 0 and 1 (experiment orbit, Table 3). Mind that solar insolation itself is not included in the PISM results underlying the conceptual model setup.
The resulting ice volume exhibits an asymmetric, sawtooth response to the forcing (Figure 7a), as is also observed in benthic 18 O records (Liebrand et al., 2017). The ice volume furthermore shows no clear relation to the forcing CO 2 -index (Figure 7b), in contrast to the equilibrium ice volume. This is due to the intertwined 10.1029/2020PA003971 frequencies with different, time-varying strengths comprising the forcing signal. The amplitude of ice volume variability is suppressed as a consequence of the dominant precession component in the forcing, which grants only short periods of growth to the ice sheet. Generally, ice volume follows the precession cycles, but some are skipped. This can happen when eccentricity decreases from one precession cycle to the next, and consequently, the amplitude of precession variability becomes smaller.
Similar to our findings in section 3.2, decreasing the amplitude of CO 2 -index variations from 1 to 0.5-now varying it between 0.25 and 0.75 (experiment orbit_half)-leads to a 7% larger mean ice volume of 5.92 compared to 5.52 × 10 6 km 3 (not shown). Similar features have been described for the influence of temperature variability on future Greenland (Mikkelsen et al., 2018) and Pleistocene Northern Hemisphere ice sheet evolution (Abe-Ouchi et al., 2013;Niu et al., 2019). We perform two more simulations with the CO 2 -index covarying with insolation between 0 and 0.5 (orbit_low; Figure 7d, orange lines) and between 0.5 and 1 (orbit_high; Figures 7c and 7d, green lines). In orbit_high, the ice volume remains very small throughout the entire simulated period. For the CO 2 range in orbit_low, however, the CO 2 -index-V eq relation is much steeper. This causes higher growth rates and consequently larger amplitude ice volume variability, because of the increased differences between ice volume and equilibrium ice volume. Particularly toward eccentricity minima, the ice volume grows substantially, hence constituting a 400-kyr eccentricity cycle, similar to the one observed in the stacked benthic 18 O record of Zachos et al. (2008) (Figure 1, red).
These simulations show how the delayed and nonlinear response of the AIS to (covarying) solar insolation and CO 2 concentrations alone can already modulate power in the benthic 18 O spectrum. Moreover, they corroborate a significant effect of the mean CO 2 and amplitude of CO 2 variability, on the relation between insolation and ice volume (Figures 7e and 7f) and hence benthic 18 O during the Miocene (Holbourn et al., 2005(Holbourn et al., , 2013(Holbourn et al., , 2015(Holbourn et al., , 2018. (e, f) Relation between 80 • S January mean insolation (Laskar et al., 2004) and ice volume, for the same experiments.

Discussion
Discrepancies between Southern Hemisphere solar insolation and benthic 18 O variability during the Miocene can be caused by the influence of deep-sea temperature changes and by CO 2 (and other greenhouse gas) concentrations not covarying with solar insolation. In section 3.3, we have shown that these discrepancies can also result from the delayed and nonlinear response of Antarctic ice volume to (covarying) solar insolation and CO 2 concentrations. In reality, however, CO 2 concentration changes do not simply follow insolation variability even on orbital timescales and can have counterbalancing effects on the AIS. In fact, AIS changes can affect the carbon cycle (Wadham et al., 2019), for example, by their effect on chemical weathering (Lear et al., 2004;Shevenell et al., 2008) and organic carbon burial rates (Badger et al., 2013;Kump & Arthur, 1999;Sosdian et al., 2020), as is reflected in benthic 13 C records. In the future, the carbon cycle should therefore be integrated in both the physically based setup (Ganopolski & Brovkin, 2017) and in the conceptual model (Paillard & Parrenin, 2004), to additionally capture interactions between ice sheets, the climate, and the carbon cycle.
Furthermore, the PISM results underlying our conceptual model setup do not include the effects of orbital variability, and they are climatically forced offline. Consequently, important climate-ice sheet feedbacks, such as the albedo-temperature feedback, are neglected (e.g., DeConto & Pollard, 2003;Gasson et al., 2016;Gregory et al., 2012;Herrington & Poulsen, 2011;Ladant et al., 2014;Stap et al., 2014;Tan et al., 2018;Verbitsky et al., 2018). In addition, GCM modeling has shown large model dependencies on the AIS temperature forcing at different CO 2 levels, stemming from the implemented physics (e.g., Gasson et al., 2014;Lunt et al., 2012). Hence, in reality, the range of CO 2 levels where significant AIS changes occur is most likely different from what the PISM simulations we have used here suggest. We have accounted for this by normalizing the CO 2 changes using an index. Moreover, hysteresis in the CO 2 -index-V eq relation might be larger (e.g., Gasson et al., 2016). Our analysis in Appendix A shows that increased hysteresis does not qualitatively affect our main results. However, how to conceptually model ice sheet changes in between the 10.1029/2020PA003971 bottom and top branches of the CO 2 -index-V eq relation remains elusive. Keeping the ice volume unchanged in between the bottom and top branches suffices to obtain reasonably close correspondence between the transient conceptual model simulations and the PISM simulations of Stap et al. (2019). With increased hysteresis, however, this will likely not be the case, necessitating a comprehensive investigation of equilibrium (coupled climate and) ice volume simulations started from "in between" initial CO 2 and ice volume states.
Lastly, inclusion of climate-ice sheet feedbacks can also affect the pace of ice volume changes and therefore the difference between the growth and decay rates and the equilibrium ice volume change rates of the Miocene AIS. Hence, our results should be re-evaluated quantitatively, once a sufficient number of equilibrium simulations from coupled ice sheet-climate models including orbital variability are available to serve as a basis for our conceptual model setup.

Conclusions
Usually, concurrent ice volume increase and CO 2 level rise are explained as a consequence of increased precipitation outweighing increased surface melt in a warmer climate (Oerlemans, 1991). This reasoning only needs to consider ice sheets in equilibrium with the forcing climate. However, here we have explored how it is possible to obtain ice volume growth during CO 2 increase without the necessity for this characteristic of the equilibrium mass balance. We have deployed a conceptual model in the vein of earlier models of transient ice volume variability (e.g., Calder, 1974;Huybers, 2009;Huybers & Wunsch, 2005;Imbrie & Imbrie, 1980;Parrenin & Paillard, 2003), to analyze the difference in tendency between the transient ice volume variability and equilibrium differences.
Our main finding is that in our simulations, the ice sheet continues to adapt to a relatively large equilibrium ice sheet after a CO 2 decrease, when the CO 2 level is already increasing again. This is purely a consequence of the disequilibrium of the ice sheet, the growth rate being smaller than the change rate of the equilibrium ice volume, in combination with the transient forcing. We have used Southern Hemisphere insolation variability as a basis for CO 2 -index forcing over the period 21 to 20 Myr ago. Our results indicate that the delayed and nonlinear response of the AIS to (covarying) solar insolation and CO 2 concentrations can lead to "skipped" forcing cycles and modulation of power in the benthic 18 O spectrum. Hence, we conclude that this response is in part the cause of discrepancies between Southern Hemisphere solar insolation and benthic 18 O variability during the Miocene.
We have further found that higher frequency CO 2 variability increases the relative phase lag between ice volume and CO 2 /temperature. Furthermore, because the growth rates are smaller than the decay rates of the Miocene AIS, the growth phase lasts longer than the decay phase under symmetrical forcing. Finally, the amplitude of CO 2 variability influences the mean ice volume, as a consequence of the growth rates being smaller than the decay rates. Decreased CO 2 variability also leads to a smaller amplitude of ice volume variability, because the growth and decay rates are not constant but depend on the forcing.
In light of our results, not only the value of CO 2 concentrations but also an accurate timing of their tendency is very important to constrain ice modeling studies. In particular, the phase relation with benthic 18 O variations is of great interest. This calls for the retrieval of high temporal resolution proxy-CO 2 and benthic 18 O records with good age control, over at least parts of the early and mid-Miocene. 10.1029/2020PA003971 (e) (f) Figure A1. Results of the conceptual model using a piecewise linear C-V eq relation with appreciable hysteresis. Equilibrium-cycle relation between the control parameter and ice volume for our reference experiment (solid lines), decreased forcing frequency (long dashed lines), and decreased forcing amplitude (dashed lines). Results for different ways of handling the situation when an ice sheet has a volume in between the bottom and top branch equilibrium ice volume: (a, b) Case 1, keeping the ice volume unchanged; (c, d) Case 2, letting it grow; and (e, f) Case 3, letting it decay.
in all investigated cases, although the duration of these phases differs among the cases. In all cases, these phases last relatively longer with increased forcing frequency ( Figure A1, solid and long dashed lines). Furthermore, the ratio between the growing and decaying phases and hence the amplitude of ice volume variability differs among the cases, but the phases of ice sheet growth always last longer than that of ice sheet decay in the equilibrium cycle. Lastly, a reduced amplitude forcing variability leads to an increase of the mean ice volume ( Figure A1, solid and dashed lines). Other than before, however, the amplitude of ice volume variability is reduced when the ice sheet is kept unchanged when its volume is in between the bottom and top branch equilibrium ice volume ( Figure A1b). In the other cases, ice volume variability remains unaffected ( Figures A1d and A1f).

Data Availability Statement
The data from the conceptual model simulations presented in this paper (https://doi.pangaea.de/10.1594/ PANGAEA.907521) and the code of the conceptual model (https://github.com/lbstap/Comics) are available and can be retrieved online.