Geologically constrained 2-million-year-long simulations of Antarctic Ice Sheet retreat and expansion through the Pliocene

Domain: earth_science
Published: 2024-08-15
Editorial outcome: minor_revision

Abstract
Pliocene global temperatures periodically exceeded modern levels, offering insights into ice sheet sensitivity to warm climates. Ice-proximal geologic records from this period provide crucial but limited glimpses of Antarctic Ice Sheet behavior. We use an ice sheet model driven by climate model snapshots to simulate transient glacial cyclicity from 4.5 to 2.6 Ma, providing spatial and temporal context for geologic records. By evaluating model simulations against a comprehensive synthesis of geologic data, we translate the intermittent geologic record into a continuous reconstruction of Antarctic sea level contributions, revealing a dynamic ice sheet that contributed up to 25 m of glacial-interglacial sea level change. Model grounding line behavior across all major Antarctic catchments exhibits an extended period of receded ice during the mid-Pliocene, coincident with proximal geologic data around Antarctica but earlier than peak warmth in the Northern Hemisphere. Marine ice sheet collapse is triggered with 1.5 °C model subsurface ocean warming.

Introduction
Based on atmospheric CO2 concentrations and global temperatures, the warm Pliocene provides an analog for current and future climate and sea level1. However, large uncertainties hamper geologic estimates of Pliocene global sea level, and paleoshoreline reconstructions are limited in their ability to resolve the relative amplitudes and timing (hemispheric phasing) of sea-level contributions from the Antarctic and Greenland ice sheets. An upper limit on Pliocene sea level remains elusive, which propagates deep uncertainty in future sea-level projections2. During the Pliocene, Antarctic Ice Sheet (AIS) behavior dominated the global sea-level signal; therefore, reconstructing ice sheet dynamics during this key time period is crucial for providing context for global sea-level reconstructions, understanding glacial stability, and improving future sea level rise projections.

Previous model explorations of AIS contribution to Pliocene sea level have simulated stable ice sheet configurations under static boundary conditions2–6; while this approach is fairly computationally straightforward, a constant climate forcing may artificially build up (or melt) ice sheets compared to a time-evolving climate7, and the equilibrium snapshot method can introduce additional uncertainty due to hysteresis in initial conditions8. This approach is also restricted to a specific time period corresponding to the specified boundary conditions; most work has focused on the mid-Piacenzian Warm Period (MPWP, 3.264–3.025 Ma)9, but Southern Hemisphere maximum insolation occurred earlier in the Pliocene (4.23 Ma)6, and other geologic proxies indicate sea-level highstands or Antarctic-proximal temperature maxima during different time intervals than the MPWP10–12. Previous modeling of time-evolving Pliocene AIS dynamics spans only short intervals13 or has been inextricably tied to the benthic δ18O record using inversion methods14,15 which assumes a linear relationship between δ18O and CO2 as well as ice volume, though this relationship is known to be complex16.

Here we use numerical ice sheet and climate modeling to explore ice sheet dynamics throughout the mid- and late Pliocene (4.5–2.6 Ma). Model simulations evolve transiently, reproducing unique patterns of glacial cyclicity as the ice sheet responds to variable climatic forcing driven by astronomical orbits and CO2 fluctuation.

Results
Modeled ice sheet behavior ranges widely due to key parameter variation (Fig. 1 and Supplementary Fig. S1). We first synthesize the available geologic records from across the Antarctic continent, compiling a suite of different data types, proxies, and geologic settings to validate and constrain our model simulations. This compilation is used to evaluate geologic fidelity: below we summarize model-data comparison results for (a) ice advance and retreat patterns, (b) extent of ice retreat, and (c) ice thickness changes. See “Methods” for a comprehensive sector-by-sector description of these datasets and the specific model evaluation criteria for each sector.

Ice advance and retreat patterns are recorded by marine geophysical data and drill core sediments, illuminating the extent and frequency of Pliocene glacial expansions across the continental shelf. In the Amundsen Sea, the West Antarctic Ice Sheet (WAIS) grew out to the continental shelf break multiple times during the early and later Pliocene (≥8 and ≥3 times, respectively), with a prolonged period of mid-Pliocene ice sheet retreat from about 4.2–3.2 Ma (the Pliocene Amundsen Sea Warm Period, PAWP)22,23. In the Ross Sea, seismic stratigraphy reveals ≥7–10 episodes of widespread Plio-Pleistocene glacial advances from both the WAIS and East Antarctic Ice Sheet (EAIS)24–26. Although the exact ages of these unconformities remain relatively unconstrained, some of these events likely correspond to glacial erosional surfaces identified at the ANDRILL-1B site during the later Pliocene (~13 advances)12,27. Poor age control precludes the definite identification of a period of prolonged Ross Sea ice retreat at the same time as in the Amundsen Sea, although ANDRILL-1B paleoenvironmental reconstructions indicate an extended warm interval from 4.5 to 3.4 Ma12, slightly earlier than the PAWP. Reconstructions of WAIS and EAIS dynamics in the Weddell Sea are extremely limited due to persistent sea ice obstructing ship access; however, the accumulation of glacially triggered debris flows on the continental shelf slope during the Pliocene suggests that the ice sheet periodically advanced to the shelf break28,29. Offshore of the Wilkes Subglacial Basin, sediment and drill cores indicate ≥12 EAIS advances to the shelf break alternating with large-scale grounding line retreat30–34. Similarly, in Prydz Bay, glacial unconformities and core data reveal periodic EAIS advances across the continental shelf during the Pliocene35–38. In summary, these datasets reconstruct a dynamic marine ice sheet that reached the continental shelf edge during many, if not most, glacial expansions, and receded during interglacials.

These geologic criteria, with varied confidence levels based on the robustness of the geologic constraint (“Methods”), are used to evaluate our ensemble of model simulations (Fig. 2). The computational effort of performing an ensemble of multimillion-year simulations requires a relatively coarse (40 km) model spatial resolution, so our interpretation of the geologic record and model-data comparison efforts are correspondingly large-scale; however, higher-resolution nested simulations demonstrate similar patterns of grounding line fluctuation (Supplementary Fig. S2). Simulations generally reproduce orbitally paced dynamic EAIS and WAIS migration across the continental shelf during the Pliocene. However, some model members with the highest sensitivity to ocean temperature, or fastest ice cliff failure rates (maximum enhancement of MICI parameters), are not able to grow sufficiently far across the continental shelf to satisfy this set of geologic constraints. In the Ross Sea region, only those simulations with lower sensitivity to ocean temperatures advance all the way to the shelf break as indicated by the geologic record. Simulations with lower sensitivity to ocean temperature and less extreme MICI parameters are able to reproduce the observed patterns of ice sheet fluctuation.

Most model simulations reconstruct a long period of ice sheet retreat during the PAWP in all catchment regions (not just the Amundsen Sea). This modeled warm interval is therefore slightly offset from the Ross Sea and Prydz Bay geologic records, with prolonged ice recession from ~4.1 to 3.2 Ma (rather than 4.5–3.4 Ma as in ref. 12, or 4.6–4.0 Ma as in ref. 38). Also, model ensemble members that generally produce sufficient glacial expansions across the Ross Sea continental shelf during the later Pliocene also advance during the early Pliocene, although direct geologic evidence for glacial expansion during that time is absent12,27.

Constraining the extent of past ice sheet retreat beyond the modern configuration requires more indirect geologic datasets. During the Pliocene, large-scale ice sheet collapse events are recorded by iceberg-rafted debris accumulation rates and sediment provenance analyses, as well as inland outcrops of open-marine sediments. Specifically, far-traveled iceberg-rafted debris pulses are attributed to destabilization events of large-scale ice collapse in the Wilkes Subglacial Basin and Aurora Subglacial Basin under warmer-than-present conditions39,40, suggesting significant grounding line retreat into these subglacial basins. Marine diatoms in Transantarctic Mountain outcrops41 have also been interpreted as indicators of ice collapse over Aurora and Wilkes subglacial basins42,43. Further evidence for inland erosion is provided by offshore records of terrigenous sediments32,44. The inland extent of retreat across Wilkes Subglacial Basin during these collapse events can be constrained by (a) εNd measurements of sediments that were eroded from geochemically distinct regions of bedrock and transported offshore45, suggesting that grounding line retreat never entered an inland source region; and (b) low cosmogenic nuclide concentrations in ANDRILL-1B sediments which preclude land exposure of much of the Transantarctic Mountain region and the southernmost part of the Wilkes Subglacial Basin46. In Prydz Bay, large-scale EAIS retreat is also indicated by inland outcrops of open-marine sediments that were deposited during periods of grounding line and ice shelf retreat by hundreds of kilometers47,48. In Aurora Subglacial Basin, however, model-data comparison is complicated by directly conflicting data-based interpretations: geophysical evidence suggests that grounding line retreat across the Aurora Subglacial Basin was limited to ~150 km inland from its modern position49, but iceberg-rafted debris pulses likely originated from larger-scale retreat in this region39,40.

Model members with little or no ice cliff failure (MICI parameters set to zero or low) do not produce sufficient grounding line retreat to satisfy the geologic evidence for large-scale grounding line retreat across Wilkes Land continental shelf or into the Wilkes Subglacial Basin. Model ensemble members with zero or low MICI parameterizations also do not drive enough ice sheet and ice shelf recession in Prydz Bay to simulate periodic open-marine environments occurring upstream of the glacially reworked diatomaceous sediment outcrops. However, model simulations with very high parameterized MICI sensitivity produce frequent ice sheet retreat into a geologically contraindicated inland source region45,46. Only simulations with intermediate MICI parameterizations are consistent with the evidence for ice sheet retreat across Wilkes Subglacial Basin as well as the Prydz Bay geologic record.

Past ice thickness changes can be reconstructed from cosmogenic nuclides measured at exposed mountain peaks. Although these terrestrial data from the Pliocene are extremely limited, Yamane et al.50 report episodes of interior Pliocene ice sheet thickening at various Antarctic nunataks. Their observations of ice sheet thickening are consistent only with models that have lower parameterized sensitivity to ocean temperatures: in these simulations, inland thickening occurs during interglacials due to precipitation, while coastal thickening occurs during glacials due to marine ice growth. Halberstadt et al.51 also used cosmogenic nuclide exposure ages to characterize the fraction of time that each elevation along a mountain peak has spent ice-covered, thus reconstructing the frequency behavior of ice sheet thinning and thickening. This approach has only been employed at the Pirrit Hills; however, the pattern of cyclic bedrock exposure at this location is not consistent with any model simulations. Model-data comparison using exposure age datasets is hampered by the coarse model resolution, which does not resolve the mountain peaks where data were collected; additionally, at the Pirrit Hills site, the ice thickness frequency dataset is integrated across a different time period as the simulations in our model ensemble.

This synthesis of Pliocene geologic data reconstructs a dynamic marine ice sheet that grew across the Antarctic continental shelf during glacial periods and retreated beyond the modern configuration during interglacials (with evidence of periodic large-scale ice sheet collapse). These detailed datasets are leveraged as sector-by-sector model evaluation criteria (“Methods”; Supplementary Table S1) and used to narrow down the full model ensemble (Fig. 1) to identify the most geologically consistent simulations (Fig. 3). In accordance with the geologic record, best-fit model simulations reproduce glacial periods of ice sheet expansion to the continental shelf edge, with episodic ice sheet retreat deep into EAIS marine basins (Fig. 3d, e).

Methods
We conduct transient Pliocene model simulations (4.5–2.6 Ma) of the AIS using the PSU-ISM17, a hybrid shallow ice/shallow shelf approximation ice sheet model (ISM) with a grounding line ice-flux formulation18 that demonstrates grounding line reversibility and reproduces theoretical and full-Stokes modeled grounding line behavior in idealized model intercomparison studies81–83. Model simulations are conducted as in ref. 2, with the description of additional techniques below.

Three key model parameters are systematically varied within plausible ranges to produce an ensemble of simulations. (1) Ice sheet sensitivity to ocean temperature (OCFACMULT = 2, 3, 4) is a dimensionless coefficient multiplying the sub-ice basal melt rates (parameterized as a quadratic dependence on temperature following ref. 84). (2) Marine ice cliff instability parameters (VCLIFF = 0, 1, 2, 3, 12, and CALVLIQ = 0, 15, 50, 100, 180) describe the maximum rate of cliff wastage rate horizontally into the ice edge due to cliff failure (VCLIFF; km yr−1; for context, terminus retreat velocities at Jakobshavn Isbrae have been measured up to ∼12 km yr-1 (see ref. 85)), and the enhancement of surface crevassing with increasing rain and surface meltwater availability (CALVLIQ; m−1 yr2). VCLIFF and CALVLIQ parameter combinations range from inactive ‘HF off’ to the maximum value “HF max” tested in ref. 60. (3) The matrix methodology for extrapolating the ice-extent weight was based on either the total grounded volume or the grounded area of the ice sheet. In other words, the ice-extent weight applied to each snapshot climatology is evaluated based on how closely the grounded ice volume (area) of the preceding timestep matches the total grounded ice volume (area) of the ice sheet topography that was used to produce that climatology. Each parameter combination produces a reasonable ice sheet configuration under preindustrial (Supplementary Fig. S4) and Last Glacial Maximum boundary conditions. Surface mass balance is calculated using a positive-degree-day scheme, with a lapse-rate correction for topographic differences. Eustatic sea level is kept at zero; we find minimal model sensitivity to eustatic sea-level variations up to 60 m scaled by Northern Hemisphere insolation. For the computational efficiency necessary for multimillion-year simulations, we use a 40 km grid resolution; higher-resolution (15 km) nested simulations produce similar grounding line fluctuation, indicating that our results are not biased by the necessarily coarse resolution of our continental ensemble (Supplementary Fig. S2).

Initial conditions for our simulations are provided by a modern ice sheet configuration2. Pliocene paleotopographic reconstructions deviate only slightly from modern86,87; sensitivity tests conducted using paleotopography produce similar modeled ice sheet fluctuations, but with slightly reduced glacial/interglacial variability, potentially due to the slightly shallower Pliocene subglacial bathymetry and coarser resolution and smoother bed topography of the paleotopographic reconstruction.

At each timestep, unique temperature and precipitation fields are provided to the ISM from a matrix of 18 climate model simulations, following the matrix method19,20. Climate model “snapshots” in the matrix were produced at two levels of atmospheric CO2 (285 and 421 ppm; Supplementary Fig. S5a, based on ref. 88), three orbital configurations from Pliocene time periods with minimum, maximum, and median Southern Hemisphere summer insolation (2.967 Ma, 2.956 Ma, and 2.892 Ma, respectively; Supplementary Fig. S5b), and three ice sheet topographies (collapsed West Antarctic Ice Sheet and loss of East Antarctic marine basins; modern; and a representative Pliocene expanded glacial topography; Supplementary Fig. S5c). For each matrix model member, a global atmospheric circulation model with a slab ocean (GENESIS GCM)89 was equilibrated under a set of unique boundary conditions. In the GCM, the vertical heat flux from the ocean to the base of sea ice is iteratively set to be proportional to 50–60°S slab ocean temperatures within the range of previously validated values. GCM output was then dynamically downscaled to a 60 km resolution over the Antarctic Ice Sheet using a regional climate model (RegCM3)90 to provide the temperature and precipitation fields passed to the ice sheet model (with temperature and precipitation lapse-rate corrections as necessary). The GCM includes a dynamic vegetation module91 rather than prescribe Pliocene-specific paleovegatation.

We also leverage the model matrix climatologies to provide time-evolving ocean temperatures to the ISM. We establish an empirical relationship between water temperatures of the upper mixed layer of the GCM ocean and subsurface (400 m water-depth) temperatures that influence grounding lines and are used to force the ice sheet model. This relationship was established from fully coupled atmosphere/ocean climate models spanning the last deglaciation92, the last interglacial (lig127k CESM experiment)93 as well as the warm Pliocene (PlioMIP2 CCSM4-Utr experiments under 400 and 560 ppm)94. CCSM4-Utr has relatively low climate sensitivity and good data-model agreement compared to the PRISM4 dataset95, and uses the same model as the Liu et al.92 deglacial simulation. CCSM4-Utr is one of the warmest PlioMIP2 ensemble members95; we use this simulation as a warm end-member to establish the ocean temperature scaling methodology. Our higher-CO2 snapshot climatologies fall directly in the middle of the PlioMIP2 ensemble spread (Supplementary Fig. S6). Using these datasets, we relate mixed-layer temperatures in the 50–60°S latitude band with 400 m temperatures of the ocean grid cells nearest the ice sheet margin (Supplementary Fig. S7). A subsurface ocean temperature scaling can therefore be calculated for each GCM 50–60°S mixed-layer temperatures within the model matrix; as the model steps forward in time, the matrix weighting scheme determines a uniform anomaly correction to a modern ocean climatology96 for each timestep.

A drawback of this ocean scaling approach is that it preserves the spatial structure of modern ocean temperatures throughout the model simulation, despite the transient and dynamic nature of ocean structure through time and potential feedbacks with ice growth97. We mitigate this issue by assuming that modern ocean temperatures are representative of warmer-than-modern times throughout the Pliocene (following the scaling approach as described above), but for colder-than-present times, the Liu et al.92 glacial-state ocean is more representative. Therefore, when GCM matrix ocean temperatures fall below modern, the ice sheet model input ocean climatology is scaled between a modern ocean96 and a glacial-state ocean at 20 ka92, as shown in Supplementary Fig. S7. This approach avoids extrapolating the significant ocean warmth currently observed offshore the Amundsen/Bellingshausen region98 throughout Pliocene glacial periods.

The matrix method relies on a high-resolution CO2 timeseries in order to establish the appropriate climatology for each timestep. Although recent work has filled in many gaps in the Pliocene CO2 record99, studies that reconstruct CO2 fluctuations at orbital-scale resolution remain limited in time (e.g., dela Vega et al. (MPWP)88; Chalk et al.100, Hönisch et al. (MPT)101; Martínez-Botí et al. (Late Pliocene)102). Splicing together proxy records introduces possible CO2 variability arising from differences between sites and between proxies, further complicating the use of a continuous proxy-based CO2 record to drive the ice sheet model through time. We therefore employ an alternative method of reconstructing past CO2 variability using a benthic δ13C-based proxy for atmospheric CO2 (cf. Lisiecki103). Benthic δ13C records in sediment cores reveal changes in deep water ventilation and deep ocean carbon storage104, and therefore can be used to infer atmospheric CO2; specifically, Lisiecki103 found that a modified δ13C gradient between the deep Pacific and intermediate North Atlantic (Δδ13C P-NA/2) correlates well with CO2 measured in ice cores across the last 800 ka. In this work, we extend the Δδ13C P-NA/2 proxy further back in time to 4.5 Ma, producing a continuous orbital-scale CO2 timeseries spanning the Pliocene.

The relationship between ocean δ13C and atmospheric CO2 becomes increasingly uncertain as we extend this proxy further back in time; changes in ocean circulation, carbon burial, and paleoproductivity can alter these δ13C records, as well as long-term trends in weathering and tectonics that impact the global carbon cycle. Although these processes may have modified the absolute values of benthic δ13C on long-term timescales, we assume that the timing of glacial/interglacial cyclicity preserved in these records remains robust. Therefore, after stacking the individual δ13C records, we scale the resultant Δδ13C P-NA/2 curve based on the mean and amplitude of boron isotope-based CO2 reconstructions during the Pliocene99 (Supplementary Fig. S8). The timing of ice sheet grounding line fluctuations is sensitive to this highly uncertain paleo-CO2 formulation, though the amplitude is robust (model simulations across a time interval where orbital-scale reconstructions are available (3.3–2.6 Ma) produce similar amplitudes of glacial cyclicity as models forced by the Δδ13C P-NA/2 CO2 proxy).

Although our δ13C-derived CO2 proxy varies widely, model ice sheet behavior does not directly mirror the CO2 timeseries forcing (Fig. 3a, b). However, the simulated period of receded ice from ~3278 to 3142 ka was likely driven by elevated CO2 in the proxy timeseries which may be an artifact of deep ocean reorganization rather than a change in deep ocean carbon, indicating elevated atmospheric CO2.

Geologic records and model-data comparison
Below we review the currently available geologic records of Pliocene ice sheet behavior, organized by region, along with a description of model-data agreement. We synthesize these records into specific data-based criteria for model evaluation (Supplementary Table S1), which we then use to assess each simulation in the model ensemble (Fig. 2). Models simulations are evaluated with respect to each criterion, producing a model-data agreement score (‘1’ for simulations deemed consistent with the geologic record, ‘0’ for a poor fit or unclear model-data comparison, or ‘-1’ for a model that violates the geologic record). Each criterion is given a confidence weighting that reflects the strength of the data interpretation or the robustness of the model-data comparison; given the wide range of data quality and ambiguity or certainty around proxies and data interpretation, the weighting factor correspondingly varies exponentially, from ‘10’ (most confident) to ‘5’ (less confident) to ‘1’ (least confident). The model simulation score is calculated from the sum of the model-data agreement scores multiplied by the weighting factor for each criterion (Fig. 2).

The computational effort of performing an ensemble of multimillion-year simulations requires a 40 km model spatial resolution, so our interpretation of the geologic record and model-data comparison efforts are correspondingly coarse resolution. For example, we interpret modeled ice sheets that extend across most of the exposed continental shelf as being consistent with geologic records of grounded ice at or near the shelf break. We make allowances for these slight discrepancies because the position of the continental shelf edge was changing throughout the Pliocene in many regions, as the ice sheet actively eroded the bed, prograded the shelf, and constructed trough mouth fans.

We also note that our simulations do not produce a significant M2 glaciation (models simulate a glacial period at 3.32 Ma, but it is not significantly stronger than other glacials). Our model is forced by time-evolving CO2 and insolation; the CO2 proxy dataset (Fig. 3b) does not indicate particularly low CO2 at this time, and although CO2 is thought to play a secondary role in triggering the M288, our model does not produce an orbitally driven glaciation either.

Amundsen Sea Embayment: Integrated seismic stratigraphy, drill core physical properties, and sedimentological data reveal dynamic WAIS behavior across the modern Amundsen Sea continental shelf during the early Pliocene (≥8 glacial advances out to the continental shelf break) and late Pliocene (≥3 glacial advances), interrupted by a period of prolonged mid-Pliocene ice sheet retreat from about 4.2–3.2 Ma, dubbed the Pliocene Amundsen Sea Warm Period (PAWP)22,23. Buried continental shelf grounding zone wedges identified from seismic records indicate that at least four glacial advances occurred during this million-year warm period, but were separated by long periods of ice sheet retreat spanning multiple glacial/interglacial cycles that buried these grounding zone wedges in hemipelagic sedimentation. This extended warm period is associated with high diatom contents and low terrigenous sedimentation rates in a sediment core spanning these seismic packages, suggesting reduced glacial processes on the Amundsen Sea continental shelf that trigger downslope transport of sediments from the shelf to the core site22.

Most model results independently correlate with this interpretation, reproducing multiple dynamic WAIS advances across the continental shelf during the early and late Pliocene, with a long period of ice sheet retreat during the PAWP (Supplementary Fig. S9). End-member models with the highest sensitivity to ocean temperature or maximum MICI parameters are not able to grow sufficiently far across the continental shelf. Model-data comparison for this region is designated “most confident” given the detailed history of glacial expansion and retreat across the Pliocene.

Ross Sea Embayment: Large-scale unconformities along the Ross Sea continental shelf record periodic erosive ice sheet advances across the continental shelf during the Pliocene. Seismic stratigraphy mapping has identified at least 7–10 episodes of widespread glacial advance of the WAIS and EAIS into the Ross Sea during the Plio-Pleistocene24–26, although the exact ages of these unconformities remain relatively unconstrained. At the ANDRILL-1B drill core in the western Ross Sea, periodic open ocean conditions alternated with grounded ice advance across this region throughout the Pliocene12,27. In all, 13 glacial erosional surfaces are identified during the later Pliocene (2.6–3.4 Ma), with no evidence for glacial advance before 3.4 Ma12,27. Paleoenvironmental reconstructions at this site reveal an extended warm interval from 4.5 to 3.4 Ma characterized by open ocean conditions, increased sea surface temperatures, and minimal marine-based ice and summer sea ice, followed by cooling and glacial expansion at about 3.3 Ma12.

Most simulations produce frequent glacial expansions beyond the modern grounding line, but only models with low sensitivity to ocean temperatures advance all the way to the shelf break as indicated by the geologic record. The limited continental shelf expansions in most simulations may be due to model resolution issues or uncertainties in sub-ice-shelf bathymetry under the Ross Ice Shelf. In general, model ensemble members that produce glacial expansions across the shelf during the later Pliocene also advanced during the early Pliocene, despite the lack of geologic evidence for glacial expansion during that time. In the Ross Sea, the model ensemble generally reproduces an extended warm interval, but slightly offset from the geologic record (~4.1–3.2, rather than 4.5–3.4 Ma as in ref. 12). Model-data comparison for this region is designated as “most confident” given the detailed history of Pliocene glacial expansion and retreat.

Wilkes Subglacial Basin: Drill core data suggest that the ice sheet periodically advanced to the continental shelf break and then retreated inland hundreds of kilometers across the Wilkes Subglacial Basin during the Pliocene. On the continental shelf, alternating diamicts and open-marine sediments reveal dynamic glacial advance and retreat behavior30,31. Offshore turbidite deposits record periods of glacial advance to the continental shelf edge (≥12 advances from 4.5 to 2.6 Ma)32,33; ice sheet advance to the shelf edge and the onset of grounding line retreat is associated with pulses of iceberg-rafted debris33,34. During warm interglacial periods, diatom-rich/bearing muds accumulated offshore, with terrigenous components sourced from far inland suggesting large-scale grounding line retreat into the subglacial basin32,44. Large-scale glacial retreat across the Wilkes Subglacial Basin is also inferred from observations of iceberg-rafted debris pulses collected offshore Prydz Bay, which are attributed to destabilization and large-scale ice collapse in the Wilkes Subglacial Basin and Aurora Subglacial Basin under warmer-than-present conditions39,40. Recent work on marine sediment provenance offshore Wilkes Subglacial Basin provides a spatial constraint on ‘large-scale’ ice collapse events45: if the ice sheet margin periodically retreated by several hundred kilometers into the Wilkes Subglacial Basin, glacial erosion of the geochemically distinct Adelie Craton region would be detected in offshore Pliocene sediments (as in warm Miocene intervals105). Erosion of significant amounts of Adelie Craton material in times of peak Pliocene warmth are, however, not observed in Wilkes Subglacial Basin provenance records32,44, suggesting near total loss of the marine-based ice in the Wilkes Subglacial Basin did not occur.

Pliocene marine diatoms have been found in the Sirius Group formation, which outcrops along the Transantarctic Mountains. Here we do not use these data as an explicit model constraint because multiple interpretations explain the presence of these diatoms; they could indicate a shallow inland sea depositional environment following WAIS collapse41, but more recently have been hypothesized as wind blown grains from ice-free Wilkes or Aurora subglacial basins42,43.

Most model simulations produce periodic glacial advance and retreat across the continental shelf, except for ensemble members with a high sensitivity to ocean warming which never advance all the way to the shelf break. Model members with MICI physics set to zero or low values (“HF off”, or “HF low”) do not retreat enough across this region to satisfy the geologic evidence for large-scale collapse32,39,40,44; however, models with maximum parameterized MICI sensitivity (‘HF max’) produce frequent ice sheet collapse into the geologically contraindicated Adelie region45 (also Shakun et al.46, see “ANDRILL catchment region” below).

In addition, all models (barring “HF off”, which does not retreat sufficiently) show an episode of retreat at around 3.55 Ma, consistent with the pulse of Wilkes-sourced ice rafted debris reaching Prydz Bay38,39. Model-data comparison in this region is split into two categories: the compilation of numerous multi-proxy geologic evidence of grounding line advance and retreat; and the geochemical constraint on maximum extent of ice retreat, both designated a “most confident” model constraint.

ANDRILL catchment region: Cosmogenic nuclide concentrations in ANDRILL-1B sediments are extremely low46. These sediments are interpreted to reflect ice dynamics across the relatively large land catchment area delivering sediments to the AND-1B core site throughout the Pliocene; the low nuclide concentrations therefore preclude large-scale or periodic land exposure of much of the Transantarctic Mountain region or the southernmost part of the Wilkes Subglacial Basin.

This is consistent with model simulations, which never deglaciate the terrestrial Transantarctic Mountain region even during episodes of marine ice sheet collapse. The indirect nature of this evidence leads to a ‘less confident’ model-data comparison designation. This dataset provides a similar constraint to the geochemical evidence for Wilkes Subglacial Basin ice sheet retreat; model ensemble members where deglaciation extends across the southernmost Wilkes Subglacial Basin also erode into the Adelie craton region, which is precluded by sediment provenance analysis45 as discussed above.

Aurora Subglacial Basin: The Aurora Subglacial Basin is substantially less studied than Wilkes Subglacial Basin. Large-scale glacial retreat across this region is also inferred from observations of far-traveled iceberg-rafted debris39,40. However, geophysical interpretations directly conflict, suggesting that grounding line retreat across the Aurora Subglacial Basin did not exceed ~150 km inland from its modern position since the Miocene49 which precludes large-scale ice sheet collapse across this region during the Pliocene.

Only model members with zero (or low) parameterized sensitivity to MICI mechanisms prevent more than 150 km of grounding line retreat during past warm periods. These parameter combinations are inconsistent with the geologic record in other catchments, supporting the interpretation of Pliocene ice loss across the Aurora Subglacial Basin39,40. Since model-data comparison in this region is ambiguous, it is designated “least confident” due to the conflicting geologic reconstructions.

Prydz Bay: Periodic glacial advances eroded large-scale glacial unconformities across the Prydz Bay continental shelf during the Pliocene35, constructing a large trough mouth fan and progradational continental slope deposits36, and depositing ice-proximal diamictites across and beyond the continental shelf37. Inland outcrops of open-marine diatomaceous sediments in the Lambert Graben rift flank47,48 suggest that these large-scale ice advances were punctuated by periods of ice sheet and ice shelf retreat by hundreds of kilometers. Although the ages of these retreat events are only loosely constrained, glacially reworked marine fossils in these outcropping formations date between 5.8 and 3.6 Ma at the Bardin Bluffs Formation47, 4.9–4.1 Ma at the Larsemann Hills106,107, and 4.2–4.1 Ma in the Vestfold Hills47,108. Offshore records similarly reveal cyclic ice sheet advance and retreat. Detailed analysis of iceberg-rafted debris accumulation rates indicates a long period of receded ice, both grounded ice and ice shelf, between 4.6 and 4.0 Ma38. After 3.3 Ma, IRD accumulation rates dramatically increase in amplitude and variability, with provenance data indicating more contribution from distal Wilkes Subglacial Basin sources40. This region is designated “most confident” due to the multi-proxy agreement of dynamic grounding line behavior.

Model ensemble members with zero, low, or intermediate MICI parameters do not produce grounding line retreat beyond the modern configuration, or ice sheet and ice shelf recession upstream of the glacially reworked open-marine sediment outcrops during warm interglacials, as inferred from the geologic record. None of the model members show reduced ice extents between 4.6 and 4.0 Ma; the modeled warm period occurs later (4.1–3.2 Ma). No modeled change in ice behavior is observable after 3.3 Ma; this may be when outer shelf erosion intensified and the trough mouth fan was constructed (3.9–3.6 Ma)109,110, which likely influenced ice sheet dynamics across this time period but are not physical processes currently represented by the model.

Weddell Sea Embayment: Persistent sea ice conditions in the Weddell Sea have long hampered data collection efforts, so only limited information is available regarding ice sheet dynamics during the Pliocene. Seismic surveys of the eastern Antarctic Peninsula (northwestern Weddell Sea), correlated to the nearby SHALDRILL drill core, reveal 10 Pliocene ice sheet grounding events on the continental shelf111, but it is unclear if these are related to WAIS fluctuations in the Weddell Sea. Seismic surveys of the Crary Trough Mouth Fan (on the other side of Weddell Sea, correlated to site 693 ODP leg 113), reveal extensive mass transport events during the Pliocene28, associated with continental slope progradation as the ice sheet periodically advanced to the outer shelf edge29.

All model simulations produce periodic glacial advance and retreat across the Weddell continental shelf; only extreme end members with high sensitivity to ocean temperatures and MICI parameters do not expand sufficiently. Model-data comparison in this region is ‘least confident’ given the lack of data coverage across the Weddell shelf, and the relatively unconstrained timing of trough mouth fan formation.

Interior nunataks: Cosmogenic nuclides measured in situ at exposed mountain peaks record episodes of interior ice sheet thickening during the Pliocene50, attributed to increased precipitation during warm interglacials. Four locations reveal episodes of Pliocene ice thickening: Sør Rondane Mountains in Dronning Maud Land (+400 m between 2.5 and 3 Ma, possibly +600, before 4 Ma); Grove Mountains in Princess Elizabeth Land (+150–220 m at least once before 3.5 Ma); McMurdo Sound-Dry Valleys area in Victoria Land (+700 at least once before 2.8 Ma); and Shackleton range in Coats Land (+750 m at least once before 3 Ma).

All model members simulate ice sheet thickness changes of many hundreds of meters at each the specified nunatak locations; however, only models with lower parameterized sensitivity to ocean temperatures produce sufficient thickening (specifically, a difference arises at Dronning Maud Land and Victoria Land sites). This model-data comparison is designated ‘least confident’ for two main reasons: (a) the model grid resolution of 40 km does not resolve the mountain peaks where data were collected; and paleo-ice dynamics around nunataks are influenced by processes occurring at the sub-grid scale112; and (b) modeled interior ice thickening is heavily influenced by parameterizations that are not explored here, for example, relating to precipitation and accumulation.

Pirrit Hills nunatak: As an ice sheet grows and shrinks across multiple glacial cycles, the cyclic exposure of cosmogenic nuclides along a mountain peak can be measured to reconstruct ‘fraction of time spent ice-covered’ at each sampled elevation113,114. These measurements, along an elevation transect, can be directly compared to the frequency behavior of model ice sheet thickness fluctuations across millions of years51. The only location with sufficiently detailed cosmogenic nuclide data to construct an ice cover frequency curve is currently the Mt. Tidd Nunatak in the Pirrit Hills.

None of the model simulations produce a similar pattern of cyclic exposure and ice cover frequency as indicated by the transect of cosmogenic nuclide exposure ages at the Pirrit Hills. This model-data comparison is designated ‘least confident’, primarily given the coarse-resolution model grid size which complicates model-data comparisons, but also due to the fact that the Pirrit Hills ice thickness frequency dataset integrates glacial fluctuations across the last ~5 Myr while our modeled ice thickness frequency curves reflect ice sheet behavior during the Pliocene only.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper represents a significant step forward in the use of this particular formulation of ice sheet dynamics. The evidence for this parameterisation of ice sheet retreat is lacking and it is an unnecessary addition to produce both modern and past observations (Edwards et al., 2019). Therefore, much more groundwork is required to test and understand its behaviour and this manuscript represents a good step in this direction. However, some significant clarifications are required in this manuscript, so that an audience not immersed in ice sheet and climate physics can understand it clearly.\n\nThe Pliocene is an important test-bed of ice sheet models, as it is the last period of similar to modern greenhouse gases, sustained increased global temperatures and the best-attested period of Antarctic ice sheet retreat. So, this a good choice of time period for this study. However, the matrix of climate models that are used within this study do not include a Pliocene climate state. Even with the ocean temperatures modified using a Last Glacial Maximum (LGM) simulation (albeit a rather old simulation) when they fall below modern temperatures, the same is not extended to warmer than modern temperatures, which could have produced a much more Pliocene like climate state. The CCSM4-Utr model’s Pliocene simulation was used to relate surface and sub-surface temperatures, so not using it in the same way as the LGM simulation needs justifying.\n\nIt needs to be made clear in the main text of the paper exactly what climate forcings are going into these transient simulations. There are only 2 atmospheric carbon dioxide levels (285 and 421ppm), three orbital configurations (please also state these time periods) and three ice sheet configurations, none of which incorporate any East Antarctic Ice Sheet retreat. This means there is very limited ice sheet climate feedbacks in East Antarctica and using a slab ocean model means there is no ocean circulation changes in these simulations. The key evidence for the Pliocene being warmer than today is increases in sea surface temperature and the best explanation for this is changes in ocean circulation, brought about by geological changes in ocean gateways (Otto-Bliesner et al., 2016).\n\nGiven the above and the fact that the Charney sensitivity of this model is on the lowest end of modern estimates, I would expect the warming of the Pliocene climate in these simulations to be particularly low. A higher sensitivity climate model with ocean circulation feedbacks and set up with Pliocene feedbacks (palaeogeography, vegetation and more ice loss), could potentially see much more warming; close to double from non-CO2 and WAIS Pliocene feedbacks (Lunt et al., 2012) and climate sensitivity could be twice as much and still be within IPCC best estimates (IPCC AR6). As such, the modelled EAIS may not require additional physics to collapse and, indeed, in the PLISMIP experiments ice sheet models without these additions are able to simulate significant retreat in East Antarctica (Dolan et al., 2018).\n\nUsing the direct proximal evidence for the Pliocene Antarctic Ice Sheet seems like a good option for evaluating these simulations. However, the threshold between the simulations that are taken as geological constrained and those dismissed is extremely narrow. The two simulations considered as “best-fit” and “highest fidelity” are inconsistent with two of the 10 lines of evidence and a poor fit to two more and you use these to state that the simulations suggest that “the Antarctic Ice Sheet is highly sensitive to past and future warm conditions”. However, one of your simulations (OC3-HFlow-area), with only 5 msle of AIS variability, is inconsistent with two lines of evidence and a poor fit to three. It seems to me that the available evidence is not of sufficient quality to rule out scenarios based on one extra unclear misfit to the data, especially with the magnitude of the claims that are produced form this. It needs to be made clear in the manuscript that the available data and these simulations really do not rule out a very stable Antarctic Ice Sheet, at least given the assessment made here.\n\nThere are a couple of significant changes that need to be implemented for this manuscript to be suitable for publication, namely:\n• The close match between the geological record and simulations with very little ice sheet fluctuation needs to be recognised much more in the manuscript. Particularly, this needs to be included in the summary statement and the abstract. I would recommend finishing the final sentence of the summary with the clause “without ruling out very stable ice sheet parameter sets”. In the abstract, I would rewrite the sentence starting on line 11 to be “Geologically consistent simulations reveal a dynamic ice sheet that could have contributed up to 25m of sea level change across glacial-interglacial cycles, although very stable ice sheet configurations cannot be ruled out.”\n• The climate model simulations that make up the matrix driving the ice sheet model needs to be included in the main text, somewhere around line 48. This description should include the CO2 levels, orbital configurations (including times represented) and Antarctic Ice Sheets incorporated (maybe with an image of the latter). It should also mention the things that aren’t included, especially Pliocene palaeogeographic changes, climate feedbacks from EAIS loss and changes in ocean circulation.\n\nThere are also a number of changes that would improve the manuscript:\n• The opening sentence of the abstract makes it sound like the Pliocene is rarely as warm as modern, whereas temperature records suggest that it was warmer for very long periods of time, at least 100s of thousands of years, if not millions of years. Please rephrase to reflect this.\n• Line 5 of the abstract states that you are using coupled ice sheet and climate models. This is only true in the most generalised of senses and not in the normal scientific meaning of this phrase. The following would be a much clearer sentence to include here. “Here we employ an ice sheet model with forcing provided by a matrix of climate model simulations to simulate transient Antarctic ice sheet cyclicity from 4.5-2.6 Ma. This provides …”\n• Line 8: From this manuscript I cannot see how the data is used to “constrain the model” or “ensure geologic compatibility”. Better verbs would be to “evaluate” and “test” respectively.\n• Line 10: I do not see any point in the manuscript where the Antarctic geological record is extrapolated. Either remove this sentence or make clear what this extrapolation is.\n• Line 26: Do you mean the Pliocene sea level signal? It currently reads as the future sea level, since this is the last sea level record you mention, which would make no sense.\n• Line 31: using a constant climate could melt the ice sheet rather than build it up, just depending where in a climate cycle is simulated.\n• Line 32: Many modelling studies (e.g. PLISMIP) incorporate the uncertainty due to hysteresis (and it is relatively straightforward to do so). Maybe add “although modelling can take this into account” at the end of the sentence.\n• A better description of how subsurface temperatures are calculated should be included in line 49.\n• Line 139: I think it is important to point out here that although your model cannot retreat into the Wilkes Subglacial Basin without MICI that other models are able to reproduce this, albeit with different climate forcings (e.g. PLISMIP)\n• Section 2.3: Please could you explain why warming episodes in your model setup do not increase the height of the Antarctic Ice Sheet, through increased precipitation. Is this something the model is unable to do or does not reproduce? This relationship is attested by present day observations, future ice sheet simulations and Pliocene climate and ice sheet modelling, so why does this not hold for your modelling setup?\n• Line 196: It is fair to say that there are significant sources of uncertainty in oxygen isotope based sea level estimates, but I am a little uncomfortable dismissing them as simply unreliable. Maybe something like “though δ18O estimates have significant uncertainties, particularly when the multiple factors that can affect the signal are considered16.”\n• Line 306: I am confused about your comments on the mPWP, given that the results presented here show the highest CO2 concentration (and therefore presumably temperatures) and the longest period of retreated Antarctic Ice Sheets during the mid-Pliocene warm period. It is very plausible that the early Pliocene was warmer than the mid-Pliocene, which is suggested by long climate records, and especially so in the Southern Hemisphere when the Isthmus of Panama was open (although you wouldn’t capture this here).\n• Line 322: I don’t know of anyone who suggested that the mPWP was the warmest period of the Pliocene, globally or in Antarctica. So please give a reference here, if you are trying to refute something in particular. In fact many people have suggested the early Pliocene was significantly warmer than the mid-Pliocene (e.g. ANDRILL or Ohneiser et al., 2020).\n• Line 378: I don’t think it makes any difference, but it is worth noting that CCSM4-Utr is one of the warmer models in PlioMIP2 and by some estimates has the strongest response to non-CO2 Pliocene boundary conditions (Haywood et al., 2020).\n• Line 385: It should be noted that modern melting of Antarctica is driven by upwelling of deeper waters penetrating into glacial troughs on the Antarctic shelf (Walker et al., 2007; “Oceanic heat transport onto the Amundsen Sea shelf through a submarine glacial trough”). So I don’t think this formulation would do a good job of simulating observed melting of Antarctica, as it is not closely related to the mixed layer temperatures. Is this something you have considered?\n• Line 393: You have data from the CCSM4-Utr Pliocene model, so please justify why the same scaling was not applied when temperatures were above modern.\n• Line 559: Please note here that there are marine sediments above the modern ice surface, which suggests significant topographic change that is not included in these models.\n• Line 599: should say “model-data” rather than “mode-data”.\n• Table 1: “Sea level amplitude” is a good metric, but most attempts to produce a Pliocene sea level estimate aim at maximum sea level. It would be great to have another column with the maximum sea level in metres above present day sea level, so the simulations can be compared to data.\n• Figure M3 is an absolutely key figure, in that it shows the forcings that go into the model, so needs to be in the main paper. It could perhaps be added as panels in Figure 2. It would also be great to have the three orbital forcing configurations that are actually simulated added to panel c with a symbol.", "reviewer_2": "Review Halberstadt et al.\nGEOLOGICALLY CONSTRAINED ANTARCTIC ICE SHEET SIMULATIONS: 2 MILLION\nYEARS OF GLACIAL RETREAT AND EXPANSION THROUGH THE WARM PLIOCENE\nHalberstadt et al. present an ensemble of paleo ice sheet simulations spanning the time from 4.5 – 2.6 Ma BP. Employing a matrix/climate index offline forcing approach and a coarse resolution model setup (40 km) they asses the evolution of the AIS during glacial-interglacial cycles and find extended periods of retreated ice sheet configuration driven by MISI and MICI but also individual periods in which the AIS attains a full glacial state. They calibrate and discuss their model results in light of the available proxy data finding good agreement with i.e. marine sediment cores suggesting sustained ice free conditions in different ocean sectors. They conclude that an intermediate parameterization of MICI in PSU-ISM yields the best results with respect to the proxy records, showing that low MICI parameters do not produce sufficient ice retreat compared to proxies and high MICI does not allow for an expanded ice sheet margin. Both cases are not in line with current estimates from paleo-proxy reconstructions. This is an important step in calibrating a highly uncertain process such as MICI in an ice sheet model. The manuscript is very well written and easy to follow despite the complexity of the approach and the various proxy records taken into account. The method section is substantial and very well written shedding light on the more involved aspects of the model setup employed. The main flaw of the study is the model resolution (40 km) which is very coarse for current standards. While I realize that simulations covering such a long time-span can be prohibitive from a computational standpoint it should be possible to at least produce a higher resolution run of the best fit ensemble-members highlighted in table 1. In his highly cited paper from 2009 Pollard and DeConto (both co-authors here) carried out 5 Ma of model simulations on 40 km including nested 10 km res simulations as well. I am aware that model physics where probably not as complex as the current version of PSU-ISM but say 15-20 km res simulations for a subset (see above) should be manageable and would improve both the quality of the paper as well as the robustness of the results. The second aspect I was missing is a PD or PI calibration of the model setup. As the authors correctly state, proxy recs for the Pliocene are sparse and sometimes of qualitative nature while we have a very good coverage of the current AIS state. I am aware that it is almost customary to reference previous instalments of the model stating that calibration has been carried out there and the model performs acceptably under well constrained reference conditions. It would be nice to see some more information on that in the supplements i.e. equilibrium thickness change and grounding zone positions for e.g. PI-cntr conditions.\n\nThis being said I congratulate the authors on an interesting study and I would suggest publication in light of my points above and some minor comments outlined below.\nBest Regards,\nJohannes Sutter\nL70 .. Best-fit model simulations are used to extrapolate pinpoint geologic records, disparate in space and time, into a continuous and geologically constrained reconstruction of AIS contribution to Pliocene sea level.\nThis is a nice idea, what about uncertainty assessment due to ISM and boundary condtions?\nL93 “trough mouth fan development”\nmaybe provide a reference as this manuscript is intended to be understandable by a wider (non-expert) audience.\nL188 “GMSL ranges from +18m during an Antarctic interglacial/Greenland glacial period (+18m from the Antarctic Ice Sheet plus 0m from an approximately modern-size Greenland) to -2m during an Antarctic glacial/Greenland interglacial (-7m from Antarctica plus +5m from Greenland), resulting in GMSL amplitude of 20m.”\nI don’t quiet understand the reasoning here. During a period of NH glacial advance (I’m not sure whether the term Greenland glacial period is fortunate here as globally your either in an interglacial or glacial but not both at the same time) antiphased to SH retreat you can’t exlude a Laurentide or fennoscandian component growing so the assumption +0m is a bit unfounded. Likewise, why capping the max Greenland contr. to 5 m, PD GIS has about 7m SLE.\nL211 “Specifically, zero or low MICI parameter values cannot generate enough grounding-line retreat across Wilkes Subglacial Basin (to produce large pulses of far-traveled iceberg rafted debris) or Prydz Bay (to deposit inland open marine sediments),”\nCouldn’t an expanded Wilkes Subglacial Basin ice cover (grounding line adv. + ice thickenning) which quickly retreats to PD conditions also produce a large far traveled iceberg armada signal). At least this would be an alternative explanation which might be as conceivable as a large inland retreat of the grounding zone.\nL219 If MICI is a process which has driven ice sheet instability in the past, a valid parameterization for these past climate states (both on the warmer and colder spectrum than PD) should also be valid for projections such as done e.g. in DeConto et al. 2016. Otherwise, one would have to change parameterisations of MICI depending on the climate state which would not be very helpful I guess. Wouldn’t it be possible to at least inform the parameter range which DeConto et al. used and therefore reduce the associated uncertainty in hindsight?\nL226 “Ice sheet fluctuations are driven by variation in climate forcing” as well as internal dynamics/feedbacks?\nL285 “In other words, our simulations indicate that the ice sheet thickening at coastal locations reported in Yamane et al.50 occurred during Pliocene glacials, while thickening at interior locations occurred during Pliocene interglacials. During these interglacial warm periods, however, ice sheet contribution to GMSL during warm periods was overwhelmingly dominated by mass loss at marine margins despite increased surface accumulation.”\nI guess this is somewhat to be expected as it basically corresponds to late Qauternery glacial interglacial thickening/thinning patterns in the interior EAIS where both ice core proxies as well as modelling suggest thickening during interglacials and thinning during glacials due to SMB changes.\nL322 “Antarctica but earlier than the MPWP, suggesting that the globally recognized MPWP may not have been the warmest Pliocene period in Antarctica”\nDid you test your model setup on the MPWP, would be interesting to see your model-response in comparison.\nL324/325 You mentioned before that it is not straightforward to assess the DeConto et al 2016 parameterization with your approach here, but it would be still nice to contextualize this at least approximately. E.g. the MICI parametertization is on the lower/higher/mid range compared to what’s been used the 2016 paper.\nMethods:\nL336 I suggest to specify this a bit more, acceptably is a very subjective term in this context. To see caveats it would be good to know whether e.g. grounding line reversibility etc. is comparable to the mean model response or belongs to the better/weaker performing models.\nL340 what parameterization is used for sub-ice basal melt rates and what is the reasoning for the scaling parameters? I am aware that this is probably mentioned somewhere in ref 2 but a short sentence here would be helpful.\nL341 maybe I’m reading this wrong, but does that mean, ice cliffs retreat 100 km, 200 km, 300 km or 1200 km per year? I’m sure this is not the case but for readers which are not familiar with the MICI parameterization in PSU-ISM it would be nice to quickly specify what this means in real wastage rates and under which conditions these are switched on.\nL349 It would be great to see some delta plots (e.g. deltaH) for PD/PI conditions and e.g. ice extent for LGM in the supplements. This gives the reader a clearer idea about potential biases in the modelling approach.\nL352/353 it would be good to have some sense of resolution-sensitivity. I suggest to run a shorter timeframe (couple of glacial-interglacial cycles e.g 200 kyr) on e.g. 15-20 km to see how the model response will vary. This should be possible computationaly for a subset of parameterizations (e.g. min/max and middle of the road). Alternatively a nested approach as in Pollard & DeConto 2009 could be reasonable.\nL354 This means there’s no model spinup? Probably not a larger issue due to multi-million year period you are investigating but would be good to mention how the initialization shock looks like.\nL423-426 the geologic constrained compiled by you for this study could be a very valuable resource for the community in general if available from a public repository. If all the data is published somewhere and can be uploaded this would be a great bonus outcome of this study.\nL427-430 I think a similar sentence should appear in the main part of the manuscript so the reader can interpret the results with all biases taken into account.\nL508 How well can ice retreat be quantified from these sediment records as opposed to just saying there was retreat but we don’t really know the baseline state of the ice sheet neither how much grl retreated. What I am getting at is whether these records are robust enough to make a statement about MICI parameterization. I assume there is grl retreat into the Wilkes Subglacial Basin and Aurora basin also without MICI, what is the reasoning behind that this is too small of a change (maybe I have missed this in the main manuscript).\nL522 It is impossible to adequately resolve the transatlantic mountain range at the resolution considered here. There will always be too much ice accumulating as small outlet glacies are insufficiently resolved and ice drainage is too slow. This is a typical problem with coarser resolution continental scale simulations so I’d be very careful interpreting results from this region.\nTable 1 it would be interesting to disentangle in the matrix which element has low confidence/poor fit due to known model biases. As you mention in the manuscript, comparing ice thickness changes at Nunataks at this resolution is not very informative, and I’m wondering why one would include it at all in this assessment.\nFigure 1, 2, 5: the 2d AIS plots are cut of at a latitude which masks the Peninsula and parts of the EAIS. I suggest to make the cut at a lower latitude.\nFigure 1 and 2: from the time series it seems that some glacial-interglacial-glacial transitions don’t show hysteresis (i.e. fast retreat, slow regrowth). It would be nice to hear the authors view on why this is (forcing, model parameterization, combination of the two)?\nFigure 3: this is a nice way of visualizing the results. Maybe add cm/yr to the colorbar label.\nFigure 4: while this is a nice way of illustrating the results I think one could move it to the supplements. It is clear from a physical point of view that deglaciation rates in cases of ice sheet collapse are much faster than growth rates (the latter being mainly fueled by a cessation of high surface and basal melt rates and not by nonlinear effects such as MISI and MICI). Seeing this quantified in a figure is nice but does not really show anything new per se.\nFigure 5 I’m a little confused by the figure caption (d) your model results show interior EAIS thickening in interglacials and in the manuscript you attribute this to higher precip rates which makes sense and is what we see in late quaternary interglacials. However, in the caption of (d) you state, that your model produces higher interior precip rates for cold glacial orbits (keeping CO2 at the same rate). Is this representative of your glacial states in the matrix method or just a special case? If it is representative than I don’t understand the thickenning in interglacials even though integrated SMB would be smaller. Please clarify and modify the respective paragraphs if needed.\nFigure 6 I realize that this is supposed to be a qualitative figure, but it is difficult to disentangle what it illustrates. To begin with, the continental shelf edge is relatively close or far away for PD reference conditions. Is the vertical axis calibrated with respect to this? Is distance measured along latitude or nearest point? Does the lower limit of the y-axis represent the continental shelf edge and the upper limit is positioned according to scale? If so the lower end of the y-axis would be zero right? It would help to add y-axis ticks as well. From an aesthetic point of view, inlet (c) could be improved.", "reviewer_3": "Halberstadt et al. present model experiments to simulate the transient behaviour of the Antarctic ice sheet through the Pliocene (4.5-2.6 Ma). They also present a compilation of geological evidence to document advances and retreats during this time period. These geological evidence are used to constrain their model results. These simulations are clearly interesting are they intrinsically represent the transient nature of this time period while the majority of previous modelling works have focused on selected boundary conditions supposedly representative of typical “Pliocene” climate (e.g. warm/cold orbit, high/low co2, etc.). I generally enjoy reading this paper as I think that the chosen methodology is sound and the results provide a better understanding of the Antarctic ice sheet dynamics through the Pliocene. However, I have a few comments that might deserved some discussions.\n\nMajor comments\n1. Uncertainties and model constraints. I appreciate the authors' effort to compile available data to constrain their model results. That being said, I found the \"constraints\" offered by this data relatively weak as they don't offer precise indication on space nor time. Since these constraints are so weak I feel that many model experiments could in principle agree with them. The authors here have explored only a very restricted model parameter ensemble (sensitivity to oceanic temperature change, MICI and a methodological parameter). I feel that the ice sheet dynamics shown here is affected by some modelling choices: the ice sheet model used, the climate forcing used to generate the climate matrix, the parameter in the regional climate model used,... On top of that the greenhouse gas concentration reconstructions is far from being perfect and is associated with a considerable uncertainty. I do not expect the authors to quantify explicitly all these uncertainties as the simulations presented here are probably numerically expensive. However, after reading the paper, I have the feeling that one could have simulate different sea level amplitudes and different waxing/waning frequency using a slightly different setup. Perhaps at least a discussion on model biases and how these biases could affect the results could be added in the manuscript? And the same for the CO2 curve?\n\n2. Sensitivity to atmospheric and oceanic climate change. From the results presented it seems that without the MICI (HF off) the model is unable to simulate grounding line retreat in East Antarctica. From Fig. 1c it also seems that no ice volume lower than PI is simulated in HF off. For the forcings representative of the Pliocene we can expect some ice sheet retreat in East Antarctica without invoking the MICI. For example in recent intercomparison exercise for future warming scenario in ISMIP6 (Seroussi et al., 2020) or for the buttressing intercomparison exercise (Sun et al., 2020), some models produce East Antarctica grounding line retreat. The model used here might present a relatively weak sensitivity to the oceanic forcing. Also, I was wondering about the atmospheric sensitivity. There is no information on surface mass balance computation in the manuscript (PDD? lapserate?...), only precipitation changes are discussed. With the forcings used, do you simulate surface melt? Does it explains part of the grounding line retreat simulated?\n\n3. Early Pliocene ice sheet volume. I was surprised to see a systematically larger ice sheet than today for the earlier part of the record (Fig. 1 and Fig. 2). Is this linked with the CO2 reconstruction based on delta13C? For this part of the simulation, we can expect a role of Northern Hemisphere ice sheets seem the climate seems favourable to their growth. I do not think that Northern Hemisphere ice sheets for this time period are supported by geological evidence... More generally in the paper it could have been useful to discuss the fact that to sustain a large scale Antarctic glaciations it might requires a cool climate that can allow for Northern Hemisphere ice growth.\n\n4. Erosion and bedrock evolution. It is hard to imagine that the bedrock has remained the same as today over these long timescales. At present there are some deeply incised trough in some areas that might be less marked during the early Pliocene for example. The amount of subglacial sediment may have been different as well. Is the bedrock used the one of the pre-industrial? The bedrock topography is a key driver for ice sheet dynamics and as such these considerations could have been more thoroughly discussed. In the model it would have been possible to use alternative and/or evolving bedrocks to partially quantify this effect?\n\nMinor comments / questions\nL. 370 Which topography the regional climate model use? And how the simulated climate by the RCM is downscaled to the ice sheet model? Using a lapserate? Is there any precipitation correction?\nFig. M1 The GCM matrix lines for deltaTmixed layer are all below the CCSM4 runs. If I understand correctly it means that the CCSM4 runs are all warmer than the GENESIS runs for this particular variable. Does a direct scaling really make sense then?\nL. 437-441 There is a potential role of Northern Hemisphere ice sheets here that is disregarded in this paper.\n\nTechnical corrections\nL.45 and L.335 The Schoof (2007) formulation of the flux at the grounding line is not really a parametrization.\nL. 341-342 The value given for VCLIFF are in m/yr right?\nL. 376 \"glacial cycle\" -> deglaciation\nFig 1. What ESL means here? Volume above flotation in time with respect to the volume above flotation for PI?\nFig. 1 shows poor quality.\n\nReference\nSchoof, C. (2007), Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res., 112, F03S28, doi:10.1029/2006JF000664.\nSeroussi, H. et al.: ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070, https://doi.org/10.5194/tc-14-3033-2020, 2020.\nSun S, Pattyn F, Simon EG, et al. Antarctic ice sheet response to sudden and sustained ice-shelf collapse (ABUMIP). Journal of Glaciology. 2020;66(260):891-904. doi:10.1017/jog.2020.67"}, "author_response": "Response to reviewer comments on:\n“Geologically constrained Antarctic Ice Sheet simulations: 2 million years of\nglacial retreat and expansion through the warm Pliocene”\nby A. R. Halberstadt et al.\nWe would like to thank the reviewers for their constructive and insightful suggestions\nwhich have helped us to improve this manuscript. Reviewer comments are shown in\nblack font, and the author response is shown in red italics. Author’s comments\nreference line numbers in the revised manuscript with Track Changes turned on (we\nalso submit a revised version with the changes fully incorporated, as well as updated\nreference formatting).\n\nReviewer #1 (Remarks to the Author):\nThis paper represents a significant step forward in the use of this particular formulation\nof ice sheet dynamics. The evidence for this parameterisation of ice sheet retreat is\nlacking and it is an unnecessary addition to produce both modern and past observations\n(Edwards et al., 2019). Therefore, much more groundwork is required to test and\nunderstand its behaviour and this manuscript represents a good step in this direction.\nHowever, some significant clarifications are required in this manuscript, so that an\naudience not immersed in ice sheet and climate physics can understand it clearly.\nWe thank Reviewer 1 for their thoughtful consideration and constructive suggestions for\nmaking our work more widely accessible.\n\nThe Pliocene is an important test-bed of ice sheet models, as it is the last period of\nsimilar to modern greenhouse gases, sustained increased global temperatures and the\nbest-attested period of Antarctic ice sheet retreat. So, this a good choice of time period\nfor this study. However, the matrix of climate models that are used within this study do\nnot include a Pliocene climate state. Even with the ocean temperatures modified using a\nLast Glacial Maximum (LGM) simulation (albeit a rather old simulation) when they fall\nbelow modern temperatures, the same is not extended to warmer than modern\ntemperatures, which could have produced a much more Pliocene like climate state. The\nCCSM4-Utr model’s Pliocene simulation was used to relate surface and sub-surface\ntemperatures, so not using it in the same way as the LGM simulation needs justifying.\nThis paragraph raises two concerns, which we address below in turn: first, that our\nmatrix of climate models do not include a Pliocene climate state; and second, that we\nuse a different ocean temperature scaling methodology for cooler-than-present times\nversus warmer-than-present times.\n(1) We interpret the reviewer’s statement that our matrix of climate models “do not\ninclude a Pliocene climate state” to refer to the observation from Fig. M1 (now Fig S7)\nthat the PLIOMIP CCSM4 ensemble members are significantly warmer than our climate\nmodel snapshots (colored vertical lines).\nThese CCSM4 simulations are, however, some of the warmest model members from\nthe PLIOMIP ensemble. In the figure below, we compare surface air and ocean\ntemperatures (SAT and SST) from our climate model snapshots to the full PLIOMIP\nensemble (reproduced from Haywood et al., 2020, Fig. 7). Here we plot global\ntemperature averages from 60N-60S (because south of 60S, our climate models\nincorporate either a glacial-size, modern, or vastly reduced ice sheet configuration; and\nthese >60S surface temperature differences impact globally averaged temperatures).\nHere we show that our climate model snapshots span the bulk of the wide-ranging\nPLIOMIP SST/SAT space. Basically, by changing insolation and ice sheet topography,\nour warmer/higher-CO Pliocene matrix snapshots (orange) span much of the climate\n2\nvariability represented in the PLIOMIP ensemble of models, although we do not\nsimulate the extreme warmth represented in the warmest three PLIOMIP model\nmembers including CCSM4-Utr. (Note that our cooler/lower-CO Pliocene snapshots\n2\n(blue) expand below the PLIOMIP range of climates; the prescribed climatic inputs for\nPLIOMIP simulations were designed to reflect mid-Pliocene Warm Period conditions,\nrather than a range of glacial/interglacial conditions as we explore here).\nFig. S6. Climate model snapshots compared to PlioMIP2 ensemble members. Climate model\nsnapshots using the GENESIS GCM with 421 ppm CO produce Pliocene climatologies (surface air\n2\ntemperatures and sea surface temperatures) that fall within the PlioMIP2 ensemble range of variability.\nClimate model snapshots with lower CO expand below the PlioMIP2 range; PlioMIP2 simulations were\n2\nconducted under 400 ppm CO concentrations to reconstruct the warm mid-Pliocene, while we explore a\n2\nwider range of glacial/interglacial conditions.\nSo if this model is one of the warmest in the PLIOMIP ensemble, why did we choose to\nuse PLIOMIP CCSM4-Utr model output to establish our SST/subsurface ocean\ntemperature scaling relationship? We made this choice because we use another\nCCSM4 simulation (Trace21k; Liu et al., 2009) to anchor the cooler end of this SAT/SST\nspace in our ocean scaling relationship, and using products from the same modeling\ntool reduces some of the uncertainty associated with different model physics.\nWe have edited the manuscript to better reflect these considerations; we now include\nthe above Figure S6, and add on L428-431: “CCSM4-Utr is one of the warmest\nPlioMIP2 ensemble members (Haywood et al., 2020); we use this simulation as a warm\nend-member to establish the ocean temperature scaling methodology. Our higher-CO\n2\nsnapshot climatologies fall directly in the middle of the PlioMIP2 ensemble spread (Fig.\nS6).”\n(2) We use the ocean scaling relationship in Fig. M1 (now Fig. S7) to apply a uniform\ntemperature anomaly to a modern ocean climatology – but only for warmer-than-present\ntime periods. For cooler-than-present time periods (e.g., Pliocene glacials), we scale\nbetween a modern ocean and a modeled glacial ocean.\nApplying the same approach for warmer time periods (e.g., scaling between a modern\nocean and a CCSM4 Pliocene ocean) would indeed produce much warmer oceans and\nin theory a stronger ocean forcing. Based on the large spread in Pliocene climates\nwithin the PlioMIP2 ensemble highlighted above, though, choosing a CCSM4-Utr\nPliocene ocean as a warmest end-member may not robustly reflect a more “Pliocene-\nlike” climate state; a PlioMIP2 multi-model-mean ocean product is arguably the best\nrepresentation of a “Pliocene-like” climate state, but this MMM product does not include\nsubsurface temperatures that the ice sheet model requires as input.\nWe agree wholeheartedly with Reviewer 1 that applying a different approach for warmer\nversus cooler time periods is methodologically messy and unsatisfying, though we still\nbelieve that it is the best choice for this work given the complexities and limitations in\navailable climate model data discussed here. We outline our logic behind this decision,\nboth here as a response to reviewer concerns, and also in the manuscript to provide\nnecessary information to readers.\n- In our original model setup, during the initial stages of this work, we applied the\nsame methodology to derive ocean input for both cooler and warmer times (using\nthe uniform temperature adjustment for all timesteps). A uniform temperature\nscaling approach applied on top of a modern ocean climatology avoids making\nassumptions about the spatial and vertical heterogeneity of Pliocene ocean\ntemperatures. Given the uncertainty in model products, and the grid resolution\nlimitations of PlioMIP2 coupled global climate models, we feel that the modern\nstructural distribution of ocean temperatures is currently the most appropriate\nassumption to make about the Pliocene ocean.\n- However, we noticed that none of our ice sheet model ensemble members ever\ngrew ice across the Amundsen/Bellinghausen continental shelf. Ice advance in\nthese models is precluded by a zone of enhanced warmth in the Amundsen &\nBellingshausen seas that is a salient feature of the modern ocean climatology\n(World Ocean Atlas; Levitus, 2012; also described by Schmidtko et al 2014).\nEven with a negative uniform temperature anomaly (e.g., during Pliocene glacial\ntimes), a model ice sheet that advances to the Ross and Weddell continental\nshelf edge still cannot grow over the Amundsen Sea region due to these warm\nocean temperatures. Given the geologic evidence that ice did advance across\nthis area in the Pliocene, we concluded that the modern spatial occurrence of\nwarmer ocean temperatures in this region was likely not present during Pliocene\nglacials, and we therefore implemented a different glacial-ocean scaling\nadjustment for cooler-than-present times that modified the spatial structure of\nocean temperatures.\n- Although we opted to implement a glacial ocean weighting approach for cooler-\nthan-present times, we preserve the warmer-than-present uniform-temperature\nscaling because (a) we have no evidence for if the modern spatial temperature\ndistribution is inappropriate for warmer times; and (b) the AIS is mostly terrestrial\nduring the warmest times, where this scaling would be most impactful, and\ntherefore the ocean warmth forcing isn’t expect to significantly impact our results.\n- Previous Pliocene ice sheet modeling work has applied a uniform temperature\ncorrection to a modern ocean climatology (e.g., DeConto & Pollard, 2016;\nDeConto et al., 2021), but no previous work has been able to simulate the\ntransient evolution of the ice sheet (without inverting from the oxygen isotope\nrecord) and so we have no precedent for an ocean temperature input scheme\nthat accounts for changing spatial structure of the oceans across\nglacial/interglacial cycles. This methodological dilemma highlights the difficulties\narising from incomplete proxy based understanding of ocean dynamics in the\ndeep past.\nEdits to L437-446: “A drawback of this ocean scaling approach is that it preserves the\nspatial structure of modern ocean temperatures throughout the model simulation,\ndespite the transient and dynamic nature of ocean structure through time and potential\nfeedbacks with ice growth (Hill et al., 2017). We mitigate this issue by assuming that\nmodern ocean temperatures are representative of warmer-than-modern times\nthroughout the Pliocene (following the scaling approach as described above), but for\ncolder-than-present times, the Liu et al. glacial-state ocean is more representative...\nThis approach avoids extrapolating the significant ocean warmth currently observed\noffshore the Amundsen/Bellingshausen region (Schmidtko et al., 2014) throughout\nPliocene glacial periods.”\nIt needs to be made clear in the main text of the paper exactly what climate forcings are\ngoing into these transient simulations. There are only 2 atmospheric carbon dioxide\nlevels (285 and 421ppm), three orbital configurations (please also state these time\nperiods) and three ice sheet configurations, none of which incorporate any East\nAntarctic Ice Sheet retreat. This means there is very limited ice sheet climate feedbacks\nin East Antarctica and using a slab ocean model means there is no ocean circulation\nchanges in these simulations. The key evidence for the Pliocene being warmer than\ntoday is increases in sea surface temperature and the best explanation for this is\nchanges in ocean circulation, brought about by geological changes in ocean gateways\n(Otto-Bliesner et al., 2016).\nOur smallest-ice topography is characterized by a collapsed West Antarctic Ice Sheet as\nwell as significant East Antarctic Ice Sheet retreat; this has been added to the text, and\nwe have also added a Figure SB showing climate forcings (CO , insolation, and matrix\n2\ntopographies). We believe that this clarification addresses the reviewer concerns above\n(i.e., we are indeed representing ice sheet/climate feedbacks in East Antarctica,\nalthough this was not clearly stated before). We now provide the details of the three\norbital configurations (Fig. S5), and have added the time periods corresponding to these\nmax/min/median configurations in the text (L52 and L409).\nThe slab ocean model does not produce full-depth ocean circulation (using a coupled\nclimate-ocean-ice-sheet model is currently computationally infeasible for so many\nequilibrated matrix snapshots or continuously for millions of years). Given these\nlimitations we think our approach is a reasonable compromise and stress that the\nuncertainties introduced by using a slab ocean are well within the range of uncertainties\nassociated with Pliocene proxy reconstructions of ocean and air temperatures.\nFigure S5: End-member parameters for the climate matrix snapshots. (a) CO 2 end-member values\nwere chosen based on the maximum and minimum CO concentrations reconstructed from an orbital-\n2\nresolution dataset (de la Vega et al., 2020; across the KM5c interval). (b) Time periods characterized by\nmaximum, minimum, and median Antarctic summer insolation values during the G17 interval provided\nobliquity, eccentricity, and precession values to construct ‘maximum-insolation orbit’, ‘minimum-insolation\norbit’, and ‘median-insolation orbit’ matrix snapshots (respective values - obliquity: 23.80°, 22.71°, 22.89°;\neccentricity: 0.0396, 0.0382, 0.0073; precession angle: 68.6°, 245.5°, 155.0°). (c) Surface topography is\neither ‘no marine ice’, with ice loss from marine-based portions of West and parts of East Antarctica\n(DeConto et al., 2021 Pliocene simulation); modern (Bedmap2); or a Pliocene glacial configuration from a\nsimulation forced by a climate model output with 285 ppm CO2, a ‘cold’ astronomical orbit, and eustatic\nsea level at -60m.\nGiven the above and the fact that the Charney sensitivity of this model is on the lowest\nend of modern estimates, I would expect the warming of the Pliocene climate in these\nsimulations to be particularly low. A higher sensitivity climate model with ocean\ncirculation feedbacks and set up with Pliocene feedbacks (palaeogeography, vegetation\nand more ice loss), could potentially see much more warming; close to double from non-\nCO2 and WAIS Pliocene feedbacks (Lunt et al., 2012) and climate sensitivity could be\ntwice as much and still be within IPCC best estimates (IPCC AR6). As such, the\nmodelled EAIS may not require additional physics to collapse and, indeed, in the\nPLISMIP experiments ice sheet models without these additions are able to simulate\nsignificant retreat in East Antarctica (Dolan et al., 2018).\nWe address the two concerns highlighted in the above paragraph separately. First, that\nour Pliocene climate simulations would have a low amount of warming – while the\nclimate sensitivity of our model is on the low end of the PlioMIP2 spectrum, it is within\nthe range of PlioMIP2 ensemble members, and is similar to CCSM4-Utr (see figure\nbelow). Fig. S6 (above) demonstrates that our climate model snapshots fall within the\nPlioMIP2 range of simulations (although we agree that with a higher climate sensitivity,\nour modeled climates would certainly be much warmer). We also note that there is\nsignificant uncertainty associated with geologic proxy reconstructions of paleo-CO , in\n2\naddition to the uncertainty around the ‘real’ Earth system climate sensitivity in the\nPliocene (e.g., the temperature response to CO doubling); this large uncertainty is\n2\nreflected by the large range of climate sensitivity values within the PlioMIP2 ensemble.\nWith vastly more resources and unlimited time, we would perform a more\ncomprehensive model ensemble and include different climate models with different\nclimate sensitivities, impose reconstructed Pliocene paleovegetation, vary more model\nparameters, etc. However, the computational resources to carry out that scope of a\nmodel intercomparison project is not feasible given the multi-million-year simulations we\ncarry out here (the PlioMIP project was focused on snapshot equilibrium simulations\nwith identical boundary conditions) so we can only point out the potential pitfalls of the\nclimate model and parameterizations we employ here. We add in L57-59: “The matrix\nmethod interpolation can account for dynamic ice sheet changes like surface lowering,\nbut it does not include changes to paleogeography or ocean circulation”.\nIn L418-419, we describe our vegetation approach: “The GCM includes a dynamic\nvegetation module (Kaplan et al., 2003) rather than prescribe Pliocene-specific\npaleovegatation”; this is consistent with PlioMIP2 (from Haywood et al., (2020): “Model\ngroups could either prescribe vegetation cover from the PRISM4 dataset… or simulate\nthe vegetation using a dynamic global vegetation model.”)\nComparison of ECS – equilibrium climate sensitivity (degrees per CO doubling) – of our model (green;\nSecond, the reviewer raises the concern that the modelled EAIS may not require\nadditional physics to collapse [if the climate forcing is not warm enough], and that\nprevious work shows that additional physics are not required to simulate EAIS retreat.\nWhile Dolan et al. (2018) do show a suite of ice sheet model experiments characterized\nby East Antarctic mass loss (their Fig. 5), these snapshots were initialized with a pre-\ncollapsed East Antarctic state – and all models actually regrow from that collapsed state\nunder warm Pliocene conditions – so the simulated East Antarctic mass loss in this\nspecific case is simply predetermined rather than triggered by the climate forcing. Thus,\nPLISMIP experiments (Dolan et al., 2018; DeBoer et al., 2015) show that warm Pliocene\nconditions are unable to initiate significant retreat in East Antarctica; ice sheet models\nwith modern initial ice sheet conditions fail to retreat into East Antarctica (Dolan et al.,\n2018, Fig. 6).\nWe note that although PLISMIP model members do not trigger EAIS retreat without\nadditional physics (nor do ISMIP6 models; Seroussi et al., 2020, see response to\nReviewer 3 for their reproduced Fig. 6), the ABUMIP model intercomparison project\n(Sun et al., 2020) does show EAIS retreat under idealized forcing. This experiment\ninvolves non-realistic extreme forcing (“ABUMIP compares ice-sheet model responses\nto decrease in buttressing by investigating the ‘end-member’ scenario of total and\nsustained loss of ice shelves. Although unrealistic, this scenario enables gauging the\nsensitivity of an ensemble of 15 ice-sheet models to a total loss of buttressing… (Sun et\nal., 2020)”.\nGiven that no other models produce EAIS retreat into marine basins (under realistic\nforcing conditions) without additional physics, we argue that the modeled Pliocene EAIS\ndoes require additional physics to collapse, regardless of the climate sensitivity of the\nclimate model. These additional physics, however, do not have to be the MICI\nmechanisms employed here; for example, Golledge et al. (2019) are able to simulate\nsignificant EAIS retreat using a sub-grid ocean melt scheme that applies subglacial melt\nupstream of grounding lines rather than MICI. Our geologic model/data comparison\nprovides a key validation for the overall sensitivity and spatial response of the ice sheet,\nwith the caveat that mechanisms other than (or in addition to) MICI could be at play.\nUsing the direct proximal evidence for the Pliocene Antarctic Ice Sheet seems like a\ngood option for evaluating these simulations. However, the threshold between the\nsimulations that are taken as geological constrained and those dismissed is extremely\nnarrow. The two simulations considered as “best-fit” and “highest fidelity” are\ninconsistent with two of the 10 lines of evidence and a poor fit to two more and you use\nthese to state that the simulations suggest that “the Antarctic Ice Sheet is highly\nsensitive to past and future warm conditions”. However, one of your simulations (OC3-\nHFlow-area), with only 5 msle of AIS variability, is inconsistent with two lines of\nevidence and a poor fit to three. It seems to me that the available evidence is not of\nsufficient quality to rule out scenarios based on one extra unclear misfit to the data,\nespecially with the magnitude of the claims that are produced form this. It needs to be\nmade clear in the manuscript that the available data and these simulations really do not\nrule out a very stable Antarctic Ice Sheet, at least given the assessment made here.\nThe reviewer brings up concerns that the best-fit simulations are not significantly\ndifferent than other ensemble members, including very stable ice sheet configurations\n(e.g., with 5 msle variability from Antarctica). This comment makes it clear that we need\nto better explain the large differences in terms of quality of the geologic data constraints\nin different sectors. With that aim, we first clarify the differences between the ‘most\nconfident’ and least confident’ sectors, and then show quantitatively that the best-fit runs\nscore significantly higher because they are more consistent with higher-confidence\ngeologic constraints. This is reflected in the revised manuscript by significant edits to\nTable 1. (We add to L110-111: “…These geologic criteria, with varied confidence levels\nbased on the robustness of the geologic constraint (Methods), are used to evaluate our\nensemble of model simulations (Table 1).” )\nWe now calculate model scores for each ensemble member by assigning values to our\npreviously qualitative categories of model/data agreement (-1=inconsistent; 0=unclear;\n1=consistent), and appropriate multiplication factors based on our confidence of the\nmodel/data comparison quality for each sector (10=most confident; 5=less confident;\n1=least confident).The large differences in these multipliers reflects the large differences\nin data robustness or appropriateness of model/data comparison. For example, the\nAmundsen Sea sector with a confidence multiplier of ‘10’ has clearly defined and\nrelatively well-dated evidence of the timing and extent of ice sheet advances across the\ncontinental shelf from both seismic and drill-core datasets. The Aurora Basin sector with\na confidence multiplier of ‘1’ has directly conflicting data-based interpretations of ice-\nsheet behavior – we don’t know whether the ice sheet retreated deep into these basins\nin the past, or perhaps didn’t retreat at all – and so we definitely don’t want to choose\nour best-fit simulations based on the constraints there. The Pirrit Hills record of ice\nsheet thickness frequency behavior has a confidence multiplier of ‘1’ because it is\nprobably close to meaningless to compare a mountainside record to model ice\nthicknesses over a 40-km grid cell. (Although these sectors where model/data\ncomparison is designated ‘least confident’ are not particularly useful or meaningful when\nit comes to our goal of constraining model simulations, we feel that it is important to\nkeep these columns in Table 1 because Pliocene geologic records of ice sheet behavior\nare so scarce, and therefore every little bit of information that can be gleaned is\npotentially valuable). This information has now been added to Methods, L480-488:\n“Model simulations are evaluated with respect to each criterion, producing a model-data\nagreement score (‘1’ for simulations deemed consistent with the geologic record, ‘0’ for\na poor fit or unclear model-data comparison, or ‘-1’ for a model that violates the geologic\nrecord). Each criterion is given a confidence weighting that reflects the strength of the\ndata interpretation or the robustness of the model-data comparison; given the wide\nrange of data quality and ambiguity or certainty around proxies and data interpretation,\nthe weighting factor correspondingly varies exponentially, from ‘10’ (most confident) to\n‘5’ (less confident) to ‘1’ (least confident). The model simulation score is calculated from\nthe sum of the model-data agreement scores multiplied by the weighting factor for each\ncriterion (Table 1).”\nNow that we have established the differences between ‘most confident’ and ‘least\nconfident’ model/data comparison between sectors, the top two (tied) best-fit runs can\nbe defended as clearly better-scoring compared to the ‘very stable’ (5 msle) simulations,\nthus supporting our interpretation that “model simulations with the highest fidelity to the\ngeologic record show large-scale [Antarctic ice sheet] fluctuations.”\n(Note: The simulation mentioned above, OC3-HFlow-area, with two ‘inconsistent’ scores\nand poor fit to three, actually has 14 msle of variability, and deglaciates somewhat into\nEast Antarctic basins, so would not be characteristic of a ‘very stable ice sheet’)\nRevised Table 1:\nThere are a couple of significant changes that need to be implemented for this\nmanuscript to be suitable for publication, namely:\n• The close match between the geological record and simulations with very little ice\nsheet fluctuation needs to be recognised much more in the manuscript. Particularly, this\nneeds to be included in the summary statement and the abstract. I would recommend\nfinishing the final sentence of the summary with the clause “without ruling out very\nstable ice sheet parameter sets”. In the abstract, I would rewrite the sentence starting\non line 11 to be “Geologically consistent simulations reveal a dynamic ice sheet that\ncould have contributed up to 25m of sea level change across glacial-interglacial cycles,\nalthough very stable ice sheet configurations cannot be ruled out.”\nWe have now quantified our previously qualitative rankings. With this additional clarity,\nwe note that the three ensemble simulations that are ‘very stable’ (e.g., produce 5 msle\nof sea level fluctuations) systematically violate or produce a poor fit to three out of the\nfive ‘more-confident’ sectors (i.e., the sectors with robust geologic data constraints),\nand\nproduce significantly lower run scores than our best-fit simulations. We therefore add to\nthe abstract (L14): “Geologically consistent simulations reveal a dynamic ice sheet that\ncould have contributed up to 25m of sea level change across glacial-interglacial cycles,\nwhile relatively stable simulations produce worse fits to data constraints.”\n• The climate model simulations that make up the matrix driving the ice sheet model\nneeds to be included in the main text, somewhere around line 48. This description\nshould include the CO2 levels, orbital configurations (including times represented) and\nAntarctic Ice Sheets incorporated (maybe with an image of the latter). It should also\nmention the things that aren’t included, especially Pliocene palaeogeographic changes,\nclimate feedbacks from EAIS loss and changes in ocean circulation.\nThis information has now been added to the main text (L51-54).\n(Note that climate feedbacks from EAIS actually are represented in the matrix method,\nsince the ‘collapsed-WAIS’ climate snapshot also includes ice loss from EAIS marine\nbasins. This is also now stated on L410).\nThere are also a number of changes that would improve the manuscript:\n• The opening sentence of the abstract makes it sound like the Pliocene is rarely as\nwarm as modern, whereas temperature records suggest that it was warmer for very\nlong periods of time, at least 100s of thousands of years, if not millions of years. Please\nrephrase to reflect this.\nRephrased: “Ice sheet dynamics during the Pliocene, when global temperatures\nperiodically exceeded and remained elevated above modern conditions…”\n• Line 5 of the abstract states that you are using coupled ice sheet and climate models.\nThis is only true in the most generalised of senses and not in the normal scientific\nmeaning of this phrase. The following would be a much clearer sentence to include\nhere. “Here we employ an ice sheet model with forcing provided by a matrix of climate\nmodel simulations to simulate transient Antarctic ice sheet cyclicity from 4.5-2.6 Ma.\nThis provides …” Rephrased accordingly.\n• Line 8: From this manuscript I cannot see how the data is used to “constrain the\nmodel” or “ensure geologic compatibility”. Better verbs would be to “evaluate” and “test”\nrespectively. Corrected.\n• Line 10: I do not see any point in the manuscript where the Antarctic geological record\nis extrapolated. Either remove this sentence or make clear what this extrapolation is.\nReplaced ‘extrapolate’ with ‘translate’; e.g., “This model-data comparison approach is\nused to translate the intermittent geologic record into a continuous reconstruction of\nAntarctic contribution to sea level.”\n• Line 26: Do you mean the Pliocene sea level signal? It currently reads as the future\nsea level, since this is the last sea level record you mention, which would make no\nsense. Yes; rephrased accordingly.\n• Line 31: using a constant climate could melt the ice sheet rather than build it up, just\ndepending where in a climate cycle is simulated. Added this consideration.\n• Line 32: Many modelling studies (e.g. PLISMIP) incorporate the uncertainty due to\nhysteresis (and it is relatively straightforward to do so). Maybe add “although modelling\ncan take this into account” at the end of the sentence. Changed to: “the equilibrium\nsnapshot method can introduce additional uncertainty due to hysteresis in initial\nconditions (Dolan et al., 2018).”\n• A better description of how subsurface temperatures are calculated should be included\nin line 49. Changed to: “Ocean temperatures are scaled from a modern climatology\nusing the matrix method weighting scheme to either apply a uniform ocean temperature\nanomaly for warmer-than-present times, or interpolate between a modern and glacial\nocean for colder-than-present times (Methods).”\n• Line 139: I think it is important to point out here that although your model cannot\nretreat into the Wilkes Subglacial Basin without MICI that other models are able to\nreproduce this, albeit with different climate forcings (e.g. PLISMIP)\nTo our knowledge other models cannot simulate retreat into the Wilkes Subglacial Basin\nor other marine East Antarctic basins without MICI or a sub-grid ocean melt scheme\nunder realistic forcing conditions (see response above).\n• Section 2.3: Please could you explain why warming episodes in your model setup do\nnot increase the height of the Antarctic Ice Sheet, through increased precipitation. Is this\nsomething the model is unable to do or does not reproduce? This relationship is\nattested by present day observations, future ice sheet simulations and Pliocene climate\nand ice sheet modelling, so why does this not hold for your modelling setup?\nYes, warming episodes in our model do increase the interior height of the ice sheet (this\nbehavior is later discussed at length in Sect. 3.4). Here in Sect. 2.3, we now clarify:\n“…in these simulations, inland thickening occurs during interglacials due to precipitation,\nwhile coastal thickening occurs during glacials due to marine ice growth.”\n• Line 196: It is fair to say that there are significant sources of uncertainty in oxygen\nisotope based sea level estimates, but I am a little uncomfortable dismissing them as\nsimply unreliable. Maybe something like “though δ18O estimates have significant\nuncertainties, particularly when the multiple factors that can affect the signal are\nconsidered16.” Corrected to “…though δ18O estimates of sea level have significant\nuncertainties.”\n• Line 306: I am confused about your comments on the mPWP, given that the results\npresented here show the highest CO2 concentration (and therefore presumably\ntemperatures) and the longest period of retreated Antarctic Ice Sheets during the mid-\nPliocene warm period. It is very plausible that the early Pliocene was warmer than the\nmid-Pliocene, which is suggested by long climate records, and especially so in the\nSouthern Hemisphere when the Isthmus of Panama was open (although you wouldn’t\ncapture this here). Removed this sentence; also removed Conclusions L349-350\naccordingly.\n• Line 322: I don’t know of anyone who suggested that the mPWP was the warmest\nperiod of the Pliocene, globally or in Antarctica. So please give a reference here, if you\nare trying to refute something in particular. In fact many people have suggested the\nearly Pliocene was significantly warmer than the mid-Pliocene (e.g. ANDRILL or\nOhneiser et al., 2020).\nIn the revised manuscript, we adopt a more nuanced tone in this manuscript regarding\nthe MPWP (see L342-343, and edits to Conclusions as well, L360-361)\n• Line 378: I don’t think it makes any difference, but it is worth noting that CCSM4-Utr is\none of the warmer models in PlioMIP2 and by some estimates has the strongest\nresponse to non-CO2 Pliocene boundary conditions (Haywood et al., 2020). Added\n(L428-429)\n• Line 385: It should be noted that modern melting of Antarctica is driven by upwelling of\ndeeper waters penetrating into glacial troughs on the Antarctic shelf (Walker et al., 2007;\n“Oceanic heat transport onto the Amundsen Sea shelf through a submarine glacial\ntrough”). So I don’t think this formulation would do a good job of simulating observed\nmelting of Antarctica, as it is not closely related to the mixed layer temperatures. Is this\nsomething you have considered?\nYes, this is the reason that we take the extra step of relating mixed layer temperatures\nto 400-m layer temperatures (responsible for melting at grounding lines) using the\ncoupled CMIP model output in Fig S7 (rather than just using ocean temperature\nanomalies from our mixed-layer ocean model).\n• Line 393: You have data from the CCSM4-Utr Pliocene model, so please justify why\nthe same scaling was not applied when temperatures were above modern. Added\nexplanation.\n• Line 559: Please note here that there are marine sediments above the modern ice\nsurface, which suggests significant topographic change that is not included in these\nmodels.\nThe model includes bedrock isostatic response to ice sheet load changes; while we\ndidn’t explicitly require that these locations be underwater (e.g., marine) in the Pliocene\ndue to the coarse grid resolution in our model, this region does experience significant\ntopographic change across glacial/interglacial cycles.\n• Line 599: should say “model-data” rather than “mode-data”. Corrected.\n• Table 1: “Sea level amplitude” is a good metric, but most attempts to produce a\nPliocene sea level estimate aim at maximum sea level. It would be great to have\nanother column with the maximum sea level in metres above present day sea level, so\nthe simulations can be compared to data. Added as a column to the revised Table 1.\n• Figure M3 is an absolutely key figure, in that it shows the forcings that go into the\nmodel, so needs to be in the main paper. It could perhaps be added as panels in Figure\n2. It would also be great to have the three orbital forcing configurations that are actually\nsimulated added to panel c with a symbol. Added to Figure 2 as suggested; additionally,\nthe end-member orbital forcing configurations (and CO concentrations) that we used\n2\nfor our climate model matrix are denoted in Fig. S5 rather than Fig. 2 for simplicity.\n\nReviewer #2 (Remarks to the Author):\nReview Halberstadt et al.\nGEOLOGICALLY CONSTRAINED ANTARCTIC ICE SHEET SIMULATIONS: 2\nMILLION\nYEARS OF GLACIAL RETREAT AND EXPANSION THROUGH THE WARM\nPLIOCENE\nHalberstadt et al. present an ensemble of paleo ice sheet simulations spanning the time\nfrom 4.5 – 2.6 Ma BP. Employing a matrix/climate index offline forcing approach and a\ncoarse resolution model setup (40 km) they asses the evolution of the AIS during\nglacial-interglacial cycles and find extended periods of retreated ice sheet configuration\ndriven by MISI and MICI but also individual periods in which the AIS attains a full glacial\nstate. They calibrate and discuss their model results in light of the available proxy data\nfinding good agreement with i.e. marine sediment cores suggesting sustained ice free\nconditions in different ocean sectors. They conclude that an intermediate\nparameterization of MICI in PSU-ISM yields the best results with respect to the proxy\nrecords, showing that low MICI parameters do not produce sufficient ice retreat\ncompared to proxies and high MICI does not allow for an expanded ice sheet margin.\nBoth cases are not in line with current estimates from paleo-proxy reconstructions. This\nis an important step in calibrating a highly uncertain process such as MICI in an ice\nsheet model. The manuscript is very well written and easy to follow despite the\ncomplexity of the approach and the various proxy records taken into account. The\nmethod section is substantial and very well written shedding light on the more involved\naspects of the model setup employed. The main flaw of the study is the model\nresolution (40 km) which is very coarse for current standards. While I realize that\nsimulations covering such a long time-span can be prohibitive from a computational\nstandpoint it should be possible to at least produce a higher resolution run of the best fit\nensemble-members highlighted in table 1. In his highly cited paper from 2009 Pollard\nand DeConto (both co-authors here) carried out 5 Ma of model simulations on 40 km\nincluding nested 10 km res simulations as well. I am aware that model physics where\nprobably not as complex as the current version of PSU-ISM but say 15-20 km res\nsimulations for a subset (see above) should be manageable and would improve both\nthe quality of the paper as well as the robustness of the results. The second aspect I\nwas missing is a PD or PI calibration of the model setup. As the authors correctly state,\nproxy recs for the Pliocene are sparse and sometimes of qualitative nature while we\nhave a very good coverage of the current AIS state. I am aware that it is almost\ncustomary to reference previous instalments of the model stating that calibration has\nbeen carried out there and the model performs acceptably under well constrained\nreference conditions. It would be nice to see some more information on that in the\nsupplements i.e. equilibrium thickness change and grounding zone positions for e.g. PI-\ncntr conditions.\nThank you for your constructive feedback to improve our work. In response to these\nsuggestions, we have added two supplemental figures. First, we show our PI validation\nof the model setup (Fig. S4, reproduced below), and referenced on L392. For brevity,\nonly a representative subset of the ensemble is plotted in Fig. S4, but all parameter\ncombinations maintain a reasonable modern grounding line after 2000 years under\nconstant preindustrial conditions (e.g., modern orbital configuration, CO , and ice-sheet\n2\ntopography).\nFig. S4. Model simulations under constant preindustrial matrix-interpolated climate. Each\nparameter combination in our ensemble is tested under preindustrial conditions. Simulations are run for\n2000 years with constant modern orbital parameters and CO concentrations set at 280 ppm, and the\n2\nresulting ice thickness and grounding line is compared to Bedmap2 (Fretwell et al., 2013). Only a subset\nof model parameter combinations is shown here; (a) Initial restart file; (b) OC3-HFlo-area; (c) OC3-HFoff-\nvol; (d) OC3-HFmedlo-area; (e) OC4-HFmedlo-vol; (f) OC4-HFmedhi-vol; (g) OC3-HFmax-area; (h) OC3-\nHFmedhi-vol.\nWe also demonstrate that our results are insensitive to model resolution; we have\nperformed higher-resolution (15km) nested model simulations for a best-fit run, and\nshow that the grounding-line patterns across the continental shelf in this 15km-\nsimulation are quite similar to the continent-wide model. Since the patterns of ice\nadvance and retreat across the continental shelf provide the model-based criteria for the\nmajority of our data constraints, we therefore conclude that our results are insensitive to\nthe necessarily coarse resolution of our continental ensemble (Fig S2; and discussed in\nL397-398).\nFig. S2. Insensitivity of results to model resolution. Higher-resolution (15km) nested simulations,\ndriven by boundary conditions from a best-fit continental run (40km resolution), produce very similar\ngrounding-line fluctuation patterns (a,b) and extents (c). (a,b) Grounding-line behavior at each timestep is\nplotted with respect to distance from the continental shelf edge (black star) along the trackline shown in\nthe left-hand-side plot (dashed black line) in the (a) Ross Sea and (b) Amundsen Sea. (c) Nested\nsimulations produce the same grounding line positions through time (right) as the continental simulation\n(left). Grounding line color indicates the ice sheet grounded area with respect to the modern extent.\nThis being said I congratulate the authors on an interesting study and I would suggest\npublication in light of my points above and some minor comments outlined below.\nBest Regards,\nJohannes Sutter\nL70 .. “Best-fit model simulations are used to extrapolate pinpoint geologic\nrecords, disparate in space and time, into a continuous and geologically constrained\nreconstruction of AIS contribution to Pliocene sea level.”\nThis is a nice idea, what about uncertainty assessment due to ISM and boundary\nconditions?\nBest-fit models robustly extrapolate the geologic record into a spatial reconstruction of\nAIS evolution (which translates to a sea level contribution), broadly independent of\nthese uncertainties.\nUncertainties around ISM physics and parameterizations, boundary conditions, etc.\ncome into play more strongly in the next step, when we use geologic constraints to try to\nidentify discrete parameter ranges / specific tipping points for the ice sheet, because\nthose specific values are dependent on the ISM and boundary conditions used. (See\nresponse to Reviewer 3 for more detail).\nL93 “trough mouth fan development”\nmaybe provide a reference as this manuscript is intended to be understandable by a\nwider (non-expert) audience.\nReplaced with: “…the accumulation of glacially-triggered debris flows on the continental\nshelf slope during the Pliocene suggests that the ice sheet periodically advanced to the\nshelf break”\nL188 “GMSL ranges from +18m during an Antarctic interglacial/Greenland glacial period\n(+18m from the Antarctic Ice Sheet plus 0m from an approximately modern-size\nGreenland) to -2m during an Antarctic glacial/Greenland interglacial (-7m from\nAntarctica plus +5m from Greenland), resulting in GMSL amplitude of 20m.”\nI don’t quite understand the reasoning here. During a period of NH glacial advance (I’m\nnot sure whether the term Greenland glacial period is fortunate here as globally your\neither in an interglacial or glacial but not both at the same time) antiphased to SH retreat\nyou can’t exlude a Laurentide or fennoscandian component growing so the assumption\n+0m is a bit unfounded. Likewise, why capping the max Greenland contr. to 5 m, PD\nGIS has about 7m SLE.\n(1) Rephrased glacial/interglacial language here to ‘ice sheet growth/retreat’\n(2) Here we assume that there were no large-scale Northern Hemisphere ice sheets\ngrowing during the Pliocene that could have masked AIS sea level contributions, just\nGreenland; we now make this assumption explicit (L207-209). (However, our ongoing\nand future work investigates this assumption!)\n(3) We cap the ‘Greenland retreat’ sea level contribution at +5m because the\n“Greenland ice sheet contributed up to ~5-7 m sea level equivalent during the Pliocene\n(ref 13: de Boer et al., 2017)” (L201 has now been edited for clarity). However, we feel it\nis equally defensible to use the 7m Present Day ice storage of GrIS, as suggested\nabove; we have redone the calculations accordingly (the resulting GMSL amplitudes are\n18m (antiphased) and 32m (phased)).\nL211 “Specifically, zero or low MICI parameter values cannot generate enough\ngrounding-line\nretreat across Wilkes Subglacial Basin (to produce large pulses of far-traveled iceberg\nrafted debris) or Prydz Bay (to deposit inland open marine sediments),”\nCouldn’t an expanded Wilkes Subglacial Basin ice cover (grounding line adv. + ice\nthickenning) which quickly retreats to PD conditions also produce a large far traveled\niceberg armada signal). At least this would be an alternative explanation which might be\nas conceivable as a large inland retreat of the grounding zone.\nThe goal of this work is to incorporate the original data-based hypotheses as faithfully\nas possible; the primary literature interpreted the IBRD signal as reflective of “large-\nscale destabilization of the East Antarctic Ice Sheet” (Cook et al., 2014). I believe this\ninterpretation is based on both the observation of IBRD pulses occurring during very\nwarm interglacials, as well as the provenance of IBRD (sourced from rocks inland of the\nmodern coastal region).\nL219 If MICI is a process which has driven ice sheet instability in the past, a valid\nparameterization for these past climate states (both on the warmer and colder spectrum\nthan PD) should also be valid for projections such as done e.g. in DeConto et al. 2016.\nOtherwise, one would have to change parameterisations of MICI depending on the\nclimate state which would not be very helpful I guess. Wouldn’t it be possible to at least\ninform the parameter range which DeConto et al. used and therefore reduce the\nassociated uncertainty in hindsight?\nFig S3 places our best-fit MICI parameter values within the context of the DeConto et al.\n(2021) future projections (RCP 8.5 sea level by 2500). We chose not to explicitly\nconstrain the uncertainty in future projections because our parameter variations don’t\nthorougly sample the parameter space (specifically between VCLIFF=3 and 13), but\nFig S1 can be interpreted accordingly. For example, none of our simulations with ‘max’\nMICI parameters (VCLIFF=13 and CALVLIQ=180) were consistent with the Pliocene\ngeologic data; therefore, our work indicates that ~25m of esl from Antarctica by 2500 is\nan overestimate. This is consistent with the best-fit Pliocene simulations producing a\nmaximum of 22m esl above present.\nL226 “Ice sheet fluctuations are driven by variation in climate forcing” as well as internal\ndynamics/feedbacks?\nReplaced with “Ice sheet fluctuations are driven by variation in climate forcing acting in\ntandem with internal feedbacks.”\nL285 “In other words, our simulations indicate that the ice sheet thickening at coastal\nlocations reported in Yamane et al.50 occurred during Pliocene glacials, while\nthickening at interior locations occurred during Pliocene interglacials. During these\ninterglacial warm periods, however, ice sheet contribution to GMSL during warm periods\nwas overwhelmingly dominated by mass loss at marine margins despite increased\nsurface accumulation.”\nI guess this is somewhat to be expected as it basically corresponds to late Qauternery\nglacial interglacial thickening/thinning patterns in the interior EAIS where both ice core\nproxies as well as modelling suggest thickening during interglacials and thinning during\nglacials due to SMB changes.\nReplaced this sentence to better reflect the novel contribution of our work: “… our\nsimulations provide temporal context for the various locations where Yamane et al.\nreconstruct higher-than-present ice elevations: thickening at their coastal sites occurred\nduring Pliocene glacials, while thickening at their interior sites occurred during Pliocene\ninterglacials, similar to ice thickness changes in the late Quaternary. Despite increased\nsurface accumulation during warm periods, however, interglacial ice sheet contribution\nto GMSL was overwhelmingly dominated by mass loss at marine margins.”\nL322 “Antarctica but earlier than the MPWP, suggesting that the globally recognized\nMPWP may not have been the warmest Pliocene period in Antarctica”\nDid you test your model setup on the MPWP, would be interesting to see your model-\nresponse in comparison.\nWe did not perform MPWP equilibrium simulations following the PlioMIP protocols;\nhowever, our model simulations do span the MPWP (Fig. 2). Due to the elevated CO\n2\nforcing during this time, we simulate a time of extended warmth. Here in Sect 3.5 we\ncontextualize this period of MPWP warmth with the rest of the simulated record of\ncontinental shelf ice sheet expansion, and notice an earlier period of extended ice sheet\nretreat from the continental shelf.\nL324/325 You mentioned before that it is not straightforward to assess the DeConto et\nal 2016 parameterization with your approach here, but it would be still nice to\ncontextualize this at least approximately. E.g. the MICI parametertization is on the\nlower/higher/mid range compared to what’s been used the 2016 paper.\nL362: “Our model-data comparison indicates that only intermediate values of modeled\nMICI parameters are consistent with the geologic record…” Here and elsewhere,\n‘intermediate’ describes the MICI parameterization with respect to previous work (e.g.,\nDeConto & Pollard 2016, DeConto et al 2021).\nMethods:\nL336 I suggest to specify this a bit more, acceptably is a very subjective term in this\ncontext. To see caveats it would be good to know whether e.g. grounding line\nreversibility etc. is comparable to the mean model response or belongs to the\nbetter/weaker performing models.\nEdited to: “…with a grounding-line ice-flux formulation that demonstrates grounding-line\nreversibility and reproduces theoretical and full-Stokes modeled grounding-line behavior\nin idealized model intercomparison studies.”\nL340 what parameterization is used for sub-ice basal melt rates and what is the\nreasoning for the scaling parameters? I am aware that this is probably mentioned\nsomewhere in ref 2 but a short sentence here would be helpful.\nEdited to: “…multiplying the sub-ice melt rates (parameterized as a quadratic\ndependence on temperature following Martin et al., 2011).”\nL341 maybe I’m reading this wrong, but does that mean, ice cliffs retreat 100 km, 200\nkm, 300 km or 1200 km per year? I’m sure this is not the case but for readers which are\nnot familiar with the MICI parameterization in PSU-ISM it would be nice to quickly\nspecify what this means in real wastage rates and under which conditions these are\nswitched on.\nVCLIFF wastage rates are km/yr; our parameter variations span 0 – 12km/yr (the typo\nhere has been corrected). We also add: “for context, terminus velocities at Jakobshavn\nIsbrae have been measured up to 12 km/yr (Joughin et al., 2012)”\n∼\nL349 It would be great to see some delta plots (e.g. deltaH) for PD/PI conditions and\ne.g. ice extent for LGM in the supplements. This gives the reader a clearer idea about\npotential biases in the modelling approach.\nAdded Fig. S4.\nL352/353 it would be good to have some sense of resolution-sensitivity. I suggest to run\na shorter timeframe (couple of glacial-interglacial cycles e.g 200 kyr) on e.g. 15-20 km\nto see how the model response will vary. This should be possible computationaly for a\nsubset of parameterizations (e.g. min/max and middle of the road). Alternatively a\nnested approach as in Pollard & DeConto 2009 could be reasonable.\nAdded Fig. S2.\nL354 This means there’s no model spinup? Probably not a larger issue due to multi-\nmillion year period you are investigating but would be good to mention how the\ninitialization shock looks like.\nInitialization shock lasts <5kyr and, for most (but not all) parameter combinations,\nmostly reflecting ice sheet growth due to increased precipitation during the Pliocene.\nL423-426 the geologic constrained compiled by you for this study could be a very\nvaluable resource for the community in general if available from a public repository. If all\nthe data is published somewhere and can be uploaded this would be a great bonus\noutcome of this study.\nWe hope that this resource will have continued utility for the wider community. For each\nregion, we list the datasets and references compiled for each region, along with the\nexplicit model evaluation criteria developed here based on this data compilation.\nL427-430 I think a similar sentence should appear in the main part of the manuscript so\nthe reader can interpret the results with all biases taken into account.\nAdded L111-115: “The computational effort of performing an ensemble of multimillion-\nyear simulations requires a relatively coarse (40 km) model spatial resolution, so our\ninterpretation of the geologic record and model-data comparison efforts are\ncorrespondingly large-scale; however, higher-resolution nested simulations demonstrate\nsimilar patterns of grounding-line fluctuation (Fig. S2).”\nL508 How well can ice retreat be quantified from these sediment records as opposed to\njust saying there was retreat but we don’t really know the baseline state of the ice sheet\nneither how much grl retreated. What I am getting at is whether these records are\nrobust enough to make a statement about MICI parameterization. I assume there is grl\nretreat into the Wilkes Subglacial Basin and Aurora basin also without MICI, what is the\nreasoning behind that this is too small of a change (maybe I have missed this in the\nmain manuscript).\nWhile we have a relatively clear boundary on the maximum extent of collapse, from the\nε data precluding extensive grounding line retreat into the Adelie Craton region, we\nNd\nagree that it is difficult to quantify the grounding line position that could have produced\nthe sedimentologic signal of “large-scale collapse” from the Wilkes Subglacial Basin.\nWe therefore split the Wilkes data into two categories and assess them separately: the\nmore clearly interpretable ε provenance constraint on max grounding line retreat\nNd\n(‘WSB’ under ‘Extent of ice retreat’ category in Table 1); and the sedimentologic/IBRD\nevidence for “large-scale collapse” (‘WSB’ under ‘Ice advance and retreat pattern’\ncategory in Table 1).\nRegarding the sedimentologic/IBRD evidence (i.e., sedimentation patterns indicating a\nsmaller-than-present ice sheet, and evidence of far-traveled iceberg rafted debris\nsourced from this region): because of the lack of a precise spatial extent of retreat,\nwe\ntake a conservative approach and interpret any model simulation with a grounding line\nreceded from the modern position by about 200km (basically, just enough to form a\nslight grounding-line embayment), as consistent with the geologic record. Because this\nis such a small change, most of the model simulations in the ensemble are ranked\n‘consistent’ (green in Table 1) in the WSB ‘Ice advance and retreat’ column. This\nparticular constraint just penalizes a few model members: those with high ocean\nsensitivity (OC4) where the simulations never advance across the continental shelf as\nindicated by the sedimentary record, or those with no or low MICI that do not retreat\nbeyond the modern configuration far enough or frequently enough to match the IBRD\nrecord.\nIs this record “robust enough to make a statement about MICI parameterization”? As\ndescribed above and shown in Table 1, we interpret this record conservatively enough\nthat it doesn’t have a lot of resolving power on its own, but in combination with other\ndatasets in other regions (e.g., the Prydz Bay data constraining min extent of retreat,\nand the ε dataset constraining max extent), we can narrow down the range of MICI\nNd\nparameters that fit the geologic data as a whole. This demonstrates the value of\nconsidering a suite of geologic records from around the continent, rather than focusing\non just one region.\nL522 It is impossible to adequately resolve the transatlantic mountain range at the\nresolution considered here. There will always be too much ice accumulating as small\noutlet glacies are insufficiently resolved and ice drainage is too slow. This is a typical\nproblem with coarser resolution continental scale simulations so I’d be very careful\ninterpreting results from this region.\nThis is a longstanding issue in continental models of Antarctica. To compensate for this,\nwe apply an ad-hoc adjustment that scales the basal sliding coefficient based on the\nsubgrid roughness (the standard deviation of bed topography within a model grid cell,\nbased on Bedmachine). This tuning is described in Pollard & DeConto (2012) and\nallows closer-to-observed flow across the Transantarctic Mountains in modern\nsimulations.\nTable 1 it would be interesting to disentangle in the matrix which element has low\nconfidence/poor fit due to known model biases. As you mention in the manuscript,\ncomparing ice thickness changes at Nunataks at this resolution is not very informative,\nand I’m wondering why one would include it at all in this assessment.\nWe include this equivocal model-data comparison purely for the sake of completeness –\nsince the available data are so limited, including everything (even if the model-data\ncomparison is not very informative) highlights both the need for more data collection as\nwell as the comprehensive nature of our geologic compilation. The nunatak column in\nTable 1 has now been shrunk significantly, which should further highlight the equivocal\nnature of this model-data comparison point (as should the explicitly quantified\nconfidence weights we have now added).\nFigure 1, 2, 5: the 2d AIS plots are cut of at a latitude which masks the Peninsula and\nparts of the EAIS. I suggest to make the cut at a lower latitude.\nOur model domain for the continental runs is quite restricted due to computational\nexpense of these multimillion-year simulations – the established domain is as small as\npossible while still providing robust representation of the ice sheet (see below) – so\nplotting these domains at a relatively high latitude encompasses most of the domain\nand highlights the areas where we performed model-data comparison (we did not use\ndata from the Antarctic Peninsula here).\nThe model domain plotted with a lower-latitude cut.\nFigure 1 and 2: from the time series it seems that some glacial-interglacial-glacial\ntransitions don’t show hysteresis (i.e. fast retreat, slow regrowth). It would be nice to\nhear the authors view on why this is (forcing, model parameterization, combination of\nthe two)?\nWe had originally elaborated on this interesting observation, but decided to remove this\ndiscussion from an earlier manuscript draft both due to length constraints and also\nbecause we plan to pursue this further in future work!\nFigure 3: this is a nice way of visualizing the results. Maybe add cm/yr to the colorbar\nlabel.\nAdded.\nFigure 4: while this is a nice way of illustrating the results I think one could move it to the\nsupplements. It is clear from a physical point of view that deglaciation rates in cases of\nice sheet collapse are much faster than growth rates (the latter being mainly fueled by a\ncessation of high surface and basal melt rates and not by nonlinear effects such as MISI\nand MICI). Seeing this quantified in a figure is nice but does not really show anything\nnew per se.\nThe new information that we intend to convey with this figure is allowing the reader to\nvisually locate the specific tipping point thresholds specified in the main text. We agree\nthat this is our least important figure, and it will be the first to move to the supplement as\nlength limits require.\nFigure 5 I’m a little confused by the figure caption (d) your model results show interior\nEAIS thickening in interglacials and in the manuscript you attribute this to higher precip\nrates which makes sense and is what we see in late quaternary interglacials. However,\nin the caption of (d) you state, that your model produces higher interior precip rates for\ncold glacial orbits (keeping CO2 at the same rate). Is this representative of your glacial\nstates in the matrix method or just a special case? If it is representative than I don’t\nunderstand the thickenning in interglacials even though integrated SMB would be\nsmaller. Please clarify and modify the respective paragraphs if needed.\nThanks for catching this error – this was a ‘typo’ and has been corrected.\nFigure 6 I realize that this is supposed to be a qualitative figure, but it is difficult to\ndisentangle what it illustrates. To begin with, the continental shelf edge is relatively close\nor far away for PD reference conditions. Is the vertical axis calibrated with respect to\nthis? Is distance measured along latitude or nearest point? Does the lower limit of the y-\naxis represent the continental shelf edge and the upper limit is positioned according to\nscale? If so the lower end of the y-axis would be zero right? It would help to add y-axis\nticks as well. From an aesthetic point of view, inlet (c) could be improved.\nThe revised figure indicates the modern (PD) conditions along transect (dashed black\nline), and each transect is now labeled on the spatial map and corresponding plot.\n(Revised Figure 6): Modeled interval of prolonged mid-Pliocene ice sheet recession in all major\ncatchments. Model grounding line position (OC2-HFmedhi-vol) in each major Antarctic catchment is\nplotted along a continental shelf transect indicated in the spatial plot. The modern (modeled) grounding\nline in each catchment region is denoted by a black dashed horizontal line. Yellow shading indicates the\n“Pliocene Amundsen Sea Warm Period”, a time when the ice sheet was primarily receded in the\nAmundsen Sea (4.2-3.2 Ma, although the modeled warm period begins a bit later at ~4.1 Ma). A similar\nperiod of prolonged mid-Pliocene ice sheet recession is evident in other catchments around the Antarctic\nIce Sheet during this time. This time period coincides with the time of greatest warmth in the Ross Sea\nAND-1B record.\n\nReviewer #3 (Remarks to the Author):\nHalberstadt et al. present model experiments to simulate the transient behaviour of the\nAntarctic ice sheet through the Pliocene (4.5-2.6 Ma). They also present a compilation\nof geological evidence to document advances and retreats during this time period.\nThese geological evidence are used to constrain their model results. These simulations\nare clearly interesting are they intrinsically represent the transient nature of this time\nperiod while the majority of previous modelling works have focused on selected\nboundary conditions supposedly representative of typical “Pliocene” climate (e.g.\nwarm/cold orbit, high/low co2, etc.). I generally enjoy reading this paper as I think that\nthe chosen methodology is sound and the results provide a better understanding of the\nAntarctic ice sheet dynamics through the Pliocene. However, I have a few comments\nthat might deserved some discussions.\nWe thank Reviewer 3 for their helpful feedback.\n\nMajor comments\n1. Uncertainties and model constraints. I appreciate the authors' effort to compile\navailable data to constrain their model results. That being said, I found the \"constraints\"\noffered by this data relatively weak as they don't offer precise indication on space nor\ntime. Since these constraints are so weak I feel that many model experiments could in\nprinciple agree with them.\nSome of the geologic data constraints are indeed weak; we have now further clarified\nthe differences between these weaker data constraints (labeled as ‘least confident’) and\nstronger datasets (‘most confident’) in our scoring scheme. We use these confidence\ncategories to ensure that our model-data comparison results are primarily influenced by\nthe most robust datasets.\nBy nature, these geologic datasets do not offer precise indications in space and time\n(e.g., the ice sheet grounding line reached 73\"S 176\"E at 4.2011 Ma), but they do\ncontain valuable information that can be assessed in space and time (e.g., there were at\nleast 13 glacial advances across the ANDRILL-1B drill site in the Ross Sea after 3.4\nMa). This is the first work, to our knowledge, that has attempted to explicitly integrate all\ncurrently available proximal geologic data to inform ice-sheet modeling. Some datasets\nhave less constraining power (e.g., for Weddell Sea: at least one grounding line\nadvance to the shelf edge), but even so, not all model members satisfy this criterion.\nOther datasets have more constraining power (e.g., no grounding-line retreat into Adelie\ncraton region), and many model members violate this constraint.\nThe authors here have explored only a very restricted model parameter ensemble\n(sensitivity to oceanic temperature change, MICI and a methodological parameter). I\nfeel that the ice sheet dynamics shown here is affected by some modelling choices: the\nice sheet model used, the climate forcing used to generate the climate matrix, the\nparameter in the regional climate model used,... On top of that the greenhouse gas\nconcentration reconstructions is far from being perfect and is associated with a\nconsiderable uncertainty. I do not expect the authors to quantify explicitly all these\nuncertainties as the simulations presented here are probably numerically expensive.\nHowever, after reading the paper, I have the feeling that one could have simulate\ndifferent sea level amplitudes and different waxing/waning frequency using a slightly\ndifferent setup. Perhaps at least a discussion on model biases and how these biases\ncould affect the results could be added in the manuscript? And the same for the CO2\ncurve?\nDue to computational constraints, we explore only a restricted number of model\nparameters, as noted here. Our goal for this restricted ensemble was to simulate\ndifferent sea level amplitudes and different waxing/waning frequencies within this\nparameter space, in order to see what patterns of esl / ice sheet fluctuation frequency\nbest match the geologic data.\nThe considerable uncertainty associated with our modeling approach impacts various\naspects of our results differentially. In this work, our primary goal is to simulate\ngeologically-consistent ice sheet evolution across Pliocene glacial/interglacial cycles.\nWe take our ‘best-fit’ simulation(s) and make interpretations about the AIS contribution\nto sea level, the spatial extent of inland retreat, etc. These results are less sensitive to\nmodel uncertainties because they are built only on our observations of modeled ice\nsheet behavior in the simulation that best reflects the geologic record – so, regardless of\nthe model climate sensitivity or potential CO errors or incomplete model physics, if a\n2\nsimulation matches geologic records through space and time around the Antarctic\ncontinent, then we can use this simulation to reconstruct time-evolving AIS behavior and\nmake interpretations about sea level contributions.\nModel uncertainties pose a much larger challenge when we use these ‘best-fit’\nsimulations to back out model parameters or discrete climate variables, for example,\nwhen we identify temperature and CO thresholds for ice sheet collapse. When used for\n2\nthis purpose, we agree that our results need more discussion of uncertainty – a climate\nmodel with a different climate sensitivity might result in a higher or lower CO threshold\n2\nfor ice sheet collapse, for example, and these thresholds may also depend on\ninaccurate or poorly parameterized model physics. We have added caveats in key\nlocations where these uncertainties need to be considered:\nL251-252 (regarding the MICI parameter constraints indicated by our model-data\ncomparison): “…our model-data comparison using spatial geologic information provides\nan upper and lower limit on geologically consistent MICI parameter values (Fig. S3),\nalthough the exact values identified here are specific to this ice sheet and climate model\nsetup with associated uncertainties.”\nL264-266 (regarding the identification of temp and CO thresholds for ice sheet\n2\ncollapse): “The specific thresholds of these forcings vary for each simulated collapse\nevent and are generally characterized as >400 ppm CO and >-20 W/m2 in our\n2\nsimulations (Fig. 3a), although the precise values of these thresholds are model-specific\nand depend on climate model sensitivity and model parameterizations.”\nWe elaborate further in Methods, L468-469: “The timing of ice sheet grounding line\nfluctuations is sensitive to this highly uncertain paleo-CO formulation, though the\n2\namplitude is robust (model simulations across a time interval where orbital-scale\nreconstructions are available (3.3-2.6 Ma) produce similar amplitudes of glacial cyclicity\nas models forced by the Δδ13CP-NA/2 CO2 proxy). Although our δ13C-derived CO\n2\nproxy varies widely, model ice sheet behavior does not directly mirror the CO time\n2\nseries forcing (Fig. 2a,b)....”\n2. Sensitivity to atmospheric and oceanic climate change. From the results presented it\nseems that without the MICI (HF off) the model is unable to simulate grounding line\nretreat in East Antarctica. From Fig. 1c it also seems that no ice volume lower than PI is\nsimulated in HF off.\nThis comment has prompted the addition of a supplemental figure showing the time-\nevolving ice sheet area (below; y-axis is ‘percent of modern grounded area’),\nnot ice\nsheet volume (esl, as in Fig. 1). Here we additionally divide the behavior of the WAIS (b)\nand EAIS (c), with the total AIS in (a). Now that we can view our results in terms of\ngrounding-line retreat rather than just sea level contributions, we demonstrate that the\nmodel does simulate some East Antarctic grounding line retreat even with HF off (blue\nlines in c below).\nFig. S1. Model ensemble runs, colored by hydrofracture parameter combinations, are plotted with respect\nto grounded ice area change compared to modern for (a) the entire ice sheet, (b) WAIS only, and (c) EAIS\nonly.\nFor the forcings representative of the Pliocene we can expect some ice sheet retreat in\nEast Antarctica without invoking the MICI. For example in recent intercomparison\nexercise for future warming scenario in ISMIP6 (Seroussi et al., 2020) or for the\nbuttressing intercomparison exercise (Sun et al., 2020), some models produce East\nAntarctica grounding line retreat. The model used here might present a relatively weak\nsensitivity to the oceanic forcing.\nIn the new supplemental figure above, we now show that the models with no MICI\n(HFoff) can simulate grounding line retreat beyond modern in EAIS basins (Fig. S1c). A\nspatial example is provided below, comparing a no-MICI model result for a warm\nPliocene time period (black) to the modern grounding line (red)\nOur no-MICI model actually produces more EAIS grounding-line retreat than ISMIP6\nprojections for 2100 (Seroussi et al., 2020, Fig, 6 reproduced below):\nIndeed, this no-MICI warm-Pliocene grounding line result is more on par with ABUMIP\n(Sun et al., 2020, Fig. 4 reproduced below), despite their non-physical experiment setup\n(designed to highlight the differential response of models under extreme idealized\nforcing). This reflects that Pliocene conditions far exceed the 2015-2100 projected\nrange of ISMIP6.\nIn summary, our model is responding to ocean forcing in a similar manner as previously\npublished model intercomparisons.\nAlso, I was wondering about the atmospheric sensitivity. There is no information on\nsurface mass balance computation in the manuscript (PDD? lapserate?...), only\nprecipitation changes are discussed. With the forcings used, do you simulate surface\nmelt? Does it explains part of the grounding line retreat simulated?\nYes, much of the MICI-driven retreat is triggered by surface melt. We have now added\nthis information to Methods. L393-394: “Surface mass balance is calculated using a\npositive-degree-day scheme, with a lapse rate correction for topographic differences.”\nAnd we now explicitly specify in L385: “…surface meltwater availability…”\n3. Early Pliocene ice sheet volume. I was surprised to see a systematically larger ice\nsheet than today for the earlier part of the record (Fig. 1 and Fig. 2). Is this linked with\nthe CO2 reconstruction based on delta13C? For this part of the simulation, we can\nexpect a role of Northern Hemisphere ice sheets seem the climate seems favourable to\ntheir growth. I do not think that Northern Hemisphere ice sheets for this time period are\nsupported by geological evidence... More generally in the paper it could have been\nuseful to discuss the fact that to sustain a large scale Antarctic glaciations it might\nrequires a cool climate that can allow for Northern Hemisphere ice growth.\nWe have added explicit mention of this larger ice sheet in the earlier part of the record in\nL221-224: “We also note that our best-fit modeled AIS sea level contributions are\nnegative (i.e., larger than today) during part of the early Pliocene. Despite larger ice\nvolumes, we simulate WAIS and even EAIS grounding line retreat (Fig. S1), as\nprecipitation outweighs interglacial marine mass loss under this early Pliocene\ncombination of relatively low CO proxy values and invariant insolation (Fig. 2).”\nmotivating question for our research. Our ultimate aim is to reconstruct both Southern\nand Northern Hemisphere transient ice sheet dynamics through the Pliocene in order to\ninvestigate the sea level evolution and hemispheric phasing of ice sheet contributions to\nsea level. This manuscript highlights the Antarctic portion of this ongoing project. We set\nthe stage for future work by adding to L207-209 “These calculations assume that the\nonly Northern Hemisphere source of ice was on Greenland, although future work will\ninvestigate potential Northern Hemisphere ice sheet growth.”\n4. Erosion and bedrock evolution. It is hard to imagine that the bedrock has remained\nthe same as today over these long timescales. At present there are some deeply incised\ntrough in some areas that might be less marked during the early Pliocene for example.\nThe amount of subglacial sediment may have been different as well. Is the bedrock\nused the one of the pre-industrial? The bedrock topography is a key driver for ice sheet\ndynamics and as such these considerations could have been more thoroughly\ndiscussed. In the model it would have been possible to use alternative and/or evolving\nbedrocks to partially quantify this effect?\nYes, this is one potential uncertainty in our approach. Although paleotopographic\nreconstructions exist, we ultimately chose to use modern topography as initial\nconditions based on sensitivity testing. L399-403: “Initial conditions for our simulations\nare provided by a modern ice sheet configuration (DeConto et al., 2021). Pliocene\npaleotopographic reconstructions deviate only slightly from modern (Paxman et al.,\n2019; Hochmuth et al., 2020); sensitivity tests conducted using paleotopography\nproduce similar modeled ice sheet fluctuations, but with slightly reduced\nglacial/interglacial variability, potentially due to the slightly shallower Pliocene sub-\nglacial bathymetry and coarser resolution and smoother bed topography of the\npaleotopographic reconstruction.” Using a modern topography also allowed us to\ndynamically downscale our continental models to a nested domain (Fig. S2).\nMinor comments / questions\nL. 370 Which topography the regional climate model use? And how the simulated\nclimate by the RCM is downscaled to the ice sheet model? Using a lapserate? Is there\nany precipitation correction?\nWe varied the GCM/RCM topography, now shown in Fig. S5c. We added to this original\nsentence “GCM output was then dynamically downscaled to a 60 km resolution over the\nAntarctic Ice Sheet using a regional climate model (RegCM3) to provide the\ntemperature and precipitation fields passed to the ice sheet model (with temperature\nand precipitation lapse-rate corrections as needed).“\nFig. M1 The GCM matrix lines for deltaTmixed layer are all below the CCSM4 runs. If I\nunderstand correctly it means that the CCSM4 runs are all warmer than the GENESIS\nruns for this particular variable. Does a direct scaling really make sense then?\nThe CCSM4 runs are among the warmest PlioMIP2 model ensemble members (our\nGENESIS snapshots fall within the PlioMIP2 range; Fig. S6). We choose to use CCSM4\nfor our scaling methodology for consistency, because it uses the same model as the\nLiu et al 2009 simulation anchoring the colder end of our climatologies. We discuss this\nissue in detail in our response to Reviewer 1.\nL. 437-441 There is a potential role of Northern Hemisphere ice sheets here that is\ndisregarded in this paper.\nThis paragraph discusses the absence of any M2 event simulated here. We decided not\nto speculate further on the lack of an M2 signal in our simulated AIS evolution because\n(1) the CO reconstruction has a lot of associated uncertainty; and (2) the M2 event has\n2\nbeen notoriously difficult to reproduce with ice sheet and climate model simulations\n(Dolan et al., 2015; ‘Modelling the enigmatic Late Pliocene Glacial Event - Marine\nIsotope Stage M2’). We therefore feel that a statement concluding that a lack of M2\nsignal in our simulations indicates a Northern Hemisphere M2 contribution would be\nunfounded.\nTechnical corrections\nL.45 and L.335 The Schoof (2007) formulation of the flux at the grounding line is not\nreally a parametrization. Corrected.\nL. 341-342 The value given for VCLIFF are in m/yr right? Corrected.\nL. 376 \"glacial cycle\" -> deglaciation Corrected.\nFig 1. What ESL means here? Volume above flotation in time with respect to the volume\nabove flotation for PI?\nAdded to Fig. 1 caption: “Simulated Pliocene Antarctic esl (equivalent GMSL\ncontribution relative to modern, calculated from the total ice amount in the domain\ndivided by global ocean area)”\nFig. 1 shows poor quality. Corrected.\nReference\nSchoof, C. (2007), Ice sheet grounding line dynamics: Steady states, stability, and\nhysteresis, J. Geophys. Res., 112, F03S28, doi:10.1029/2006JF000664.\nSeroussi, H. et al.: ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet\nevolution over the 21st century, The Cryosphere, 14, 3033–\n3070, https://doi.org/10.5194/tc-14-3033-2020, 2020.\nSun S, Pattyn F, Simon EG, et al. Antarctic ice sheet response to sudden and sustained\nice-shelf collapse (ABUMIP). Journal of Glaciology. 2020;66(260):891-904.\ndoi:10.1017/jog.2020.67\n\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have done an excellent job at meeting the review comments on the original\nmanuscript submission. Particularly, the additional supplementary figures really make\nclear the climatological forcing in these simulations. Figure S6 is an excellent addition,\nbut would be improved with a few technical alterations. It would be helpful to have the\naxes crossing at the origin, so it is obvious the model is getting colder (I assume that all\nthese temperatures are relative to pre-industrial, but that also needs clarifying). It should\nalso be pointed out that the GENESIS model is a slab ocean model, unlike the PlioMIP\nmodels that include a fully coupled ocean atmosphere. If this is not made clear then the\nGENESIS relationship between SAT and SST looks anomalous. The caption should\nalso make clear that the PlioMIP simulations are modelling isotope stage KM5c (not\nsome generic mid-Pliocene), when the orbital forcing is very close to modern and the\nPliocene warming is not at a peak (at least in the benthic oxygen isotope curve). With a\nfew technical alterations to this new figure, I have no hesitation in recommending this\nmanuscript for publication in Nature Communications.\n\nReviewer #3 (Remarks to the Author):\nFor this review, the authors have responded clearly and precisely to all my previous comments. I\nam really grateful for that. I think that the paper can be published in Nature Communications\nwhere it will make a valuable contribution to our understanding of past Antarctic ice sheet\ndynamics.\n\nResponse to second-round reviewer comments on:\n“Geologically constrained Antarctic Ice Sheet simulations: 2 million years of\nglacial retreat and expansion through the warm Pliocene”\nby A. R. Halberstadt et al.\nReviewer #1:\nThe authors have done an excellent job at meeting the review comments on the original\nmanuscript submission. Particularly, the additional supplementary figures really make\nclear the climatological forcing in these simulations. Figure S6 is an excellent addition,\nbut would be improved with a few technical alterations. It would be helpful to have the\naxes crossing at the origin, so it is obvious the model is getting colder (I assume that all\nthese temperatures are relative to pre-industrial, but that also needs clarifying). It should\nalso be pointed out that the GENESIS model is a slab ocean model, unlike the PlioMIP\nmodels that include a fully coupled ocean atmosphere. If this is not made clear then the\nGENESIS relationship between SAT and SST looks anomalous. The caption should\nalso make clear that the PlioMIP simulations are modelling isotope stage KM5c (not\nsome generic mid-Pliocene), when the orbital forcing is very close to modern and the\nPliocene warming is not at a peak (at least in the benthic oxygen isotope curve). With a\nfew technical alterations to this new figure, I have no hesitation in recommending this\nmanuscript for publication in Nature Communications.\n\nEdits to Fig. S6 are shown below.\nFigure S6. Climate model snapshots compared to PlioMIP2 ensemble members.\nClimate model snapshots using the GENESIS GCM with 421 ppm CO produce\nwithin the PlioMIP2 ensemble range of variability7. Note that the GENESIS GCM has a\nslab ocean model, while the PlioMIP2 simulations are fully coupled ocean-atmosphere\nmodels. Climate model snapshots with lower CO expand below the PlioMIP2 range;\n2\nthe KM5c interval (when orbital forcing was close to modern, though Pliocene warming\nwas not peak), while we explore a wider range of glacial/interglacial conditions.\nFrom the Author Checklist (‘One final reviewer comment to address in point-by-\npoint response’):\nThe timing of East Antarctic grounding line retreat at 15 km resolution (Fig. S2, right) is\nmasked. It could be important to show this region because the MICI could play a role\nthere.\nThe Fig S2 spatial plot in question (reproduced below), rather than masking, shows the\ntwo nested domains in color (these domains are restricted in extent and are downscaled\nusing boundary conditions from the continent-wide simulation). The two nested domains\nshown in (a) and (b) are run at a higher (15km) resolution across the ~2 million year\ntime interval, spanning the Amundsen Sea and the Ross Sea regions. We did not nest\nany domains over East Antarctica since the goal of this analysis was to investigate any\nresolution dependence in our simulated grounding-line dynamics. We do not expect the\nMICI formulation to be sensitive model grid resolution because our parameterization is\ndriven by only by meltwater availability and ice cliff height.\nFigure S2. Insensitivity of results to model resolution. Higher-resolution (15km)\nnested simulations, driven by boundary conditions from a best-fit continental run (40km\nresolution), produce very similar grounding-line fluctuation patterns (a,b) and extents\n(c). (a,b) Grounding-line behavior at each timestep is plotted with respect to distance\nfrom the continental shelf edge (black star) along the trackline shown in the left-hand-\nside plot (dashed black line) in the (a) Ross Sea and (b) Amundsen Sea. (c) Nested\nsimulations produce the same grounding line positions through time (right) as the\ncontinental simulation (left). Grounding line color indicates the ice sheet grounded area\nwith respect to the modern extent."}]

Decision letter