Antarctic climate and ice-sheet configuration during the early Pliocene interglacial at 4.23 Ma

Abstract

The geometry of Antarctic ice sheets during warm periods of the geological past is difficult to determine from geological evidence, but is important to know because such reconstructions enable a more complete understanding of how the ice-sheet system responds to changes in climate. Here we investigate how Antarctica evolved under orbital and greenhouse gas conditions representative of an interglacial in the early Pliocene at 4.23 Ma, when Southern Hemisphere insolation reached a maximum. Using offline-coupled climate and ice-sheet models, together with a new synthesis of high-latitude palaeoenvironmental proxy data to define a likely climate envelope, we simulate a range of ice-sheet geometries and calculate their likely contribution to sea level. In addition, we use these simulations to investigate the processes by which the West and East Antarctic ice sheets respond to environmental forcings and the timescales over which these behaviours manifest. We conclude that the Antarctic ice sheet contributed 8.6 $\pm$ 2.8 m to global sea level at this time, under an atmospheric CO${}_{\mathrm{2}}$ concentration identical to present (400 ppm). Warmer-than-present ocean temperatures led to the collapse of West Antarctica over centuries, whereas higher air temperatures initiated surface melting in parts of East Antarctica that over one to two millennia led to lowering of the ice-sheet surface, flotation of grounded margins in some areas, and retreat of the ice sheet into the Wilkes Subglacial Basin. The results show that regional variations in climate, ice-sheet geometry, and topography produce long-term sea-level contributions that are non-linear with respect to the applied forcings, and which under certain conditions exhibit threshold behaviour associated with behavioural tipping points.

Introduction

The response of the Antarctic ice sheets (AIS) to predicted future oceanic and atmospheric warming will dictate the magnitude of global sea level changes for millennia, yet the sensitivity of the AIS system is unclear, leading to a wide range of sea-level predictions for coming centuries . Simulations of future scenarios such as these are most credible when constrained by observations of past changes. Here we simulate the AIS during a warmer period of the geological past for which geological and paleoenvironmental data exist with which to verify simulations. Warm periods of the Pliocene (2.58–5.33 Ma) are considered some of the most appropriate analogues for future environmental conditions , especially the early to mid-Pliocene (4–5 Ma), when atmospheric CO${}_{\mathrm{2}}$ concentrations were in the range 365–415 ppm , similar to present day. Globally averaged surface temperatures during Pliocene interglacials were 2 to $>$ 3 ${}^{\circ }$C warmer than present , comparable to warmings anticipated by 2100 under mid-range emissions scenarios of the Intergovernmental Panel on Climate Change Fifth Assessment Report . Far-field sea-level records imply a collapse of the West Antarctic Ice Sheet (WAIS) together with partial loss of marine-based ice from the East Antarctic Ice Sheet (EAIS) at this time . Geological records from close to the East Antarctic coast suggest that a portion of this sea-level contribution may have originated from the Wilkes Subglacial Basin (WSB), when the ice sheet margin is thought to have migrated hundreds of kilometres inland .

Given the large scale of these inferred ice-sheet reconfigurations during warm periods of the Pliocene, it may be inferred that Antarctica's marine-based ice sheets may be vulnerable to thresholds beyond which an abrupt change in state occurs see. If such thresholds exist, and were crossed during Pliocene glacial–interglacial transitions, simulations of ice-sheet evolution during this period may be able to provide insights into processes that could be relevant under future warmer climate scenarios. To investigate the sensitivity of Antarctica to climate–ice-sheet thresholds under warmer climates, we first run global and regional climate models to simulate changes in Antarctic oceanic and atmospheric temperatures under an orbital configuration representative of an interglacial in the early Pliocene at 4.23 Ma, when austral insolation reached a maximum, and atmospheric CO${}_{\mathrm{2}}$ concentration was around 400 ppm (Figs.  and ). We focus on this period, rather than the commonly investigated “PRISM” interval of the mid-Pliocene , as the higher insolation and atmospheric CO${}_{\mathrm{2}}$ concentration at this earlier time may have resulted in a warmer climate and a smaller ice sheet than the later period, which may help reconcile far-field sea-level records for the period. Outputs from these climate simulations are then used as inputs to ice-sheet model experiments that use the Parallel Ice Sheet Model, a fully coupled ice-sheet–ice-shelf model . These simulations allow us to quantify the continental-scale changes that take place under the prescribed climatology, and the consequent contribution to sea level from the AIS. Although the peak climatic conditions we simulate with the regional-scale climate model (RCM) most likely only persisted for 1–2 kyr, we run our simulations for 10 kyr in order to quantify how the perturbed ice sheet might evolve if such warm conditions were maintained. Although this does not realistically reflect climatic forcings associated with orbital cyclicities, and implies that our sea-level estimates may be upper bounds for the prescribed climates, it is nonetheless a useful approach for the study of “unforced” ice-sheet oscillations, or “tipping points” . In the final part of our study, therefore, we use our long-term ice-sheet simulations to investigate whether evidence for tipping points can be seen, and if so, the mechanisms that give rise to them.

Orbital components and net January insolation at 80${}^{\circ }$ S for the period 3.5–5 Ma. Top panel shows insolation for the interglacial in which peak Pliocene values occur (grey shading). The insolation peak occurs at 4.23 Ma. Data from .

[Figure omitted. See PDF]

Air temperature, precipitation, and sea surface temperature used for our simulations. Top row: the “control” (CTRL) scenario, employing gridded data from observations and modelling of present-day conditions . Middle row: present-day conditions as simulated by the RCM. Bottom row: Pliocene conditions at 4.23 Ma as simulated by the RCM.

[Figure omitted. See PDF]

Methods

Climate and ocean inputs

We use spatially distributed air temperature, precipitation, and sea-surface temperatures from an established RCM coupled with the GENESIS version 3.0 global climate model (GCM) . The GCM parameterizes fluxes between the land surface boundary and the free atmosphere and includes detailed representations of snow and land ice. Outputs from the GENESIS model have been validated against observed polar climates and ice-sheet mass balance , and are therefore considered reliable for studies such as ours. Our RCM is RegCM3, which takes time-dependent lateral boundary conditions from the GCM. It simulates atmospheric dynamics, radiative transfer, and precipitation, and includes a representation of the open ocean . RCM simulations at 40 km resolution are initialized from coarser-resolution (T31) GCM outputs. We simulate the climate during a peak-insolation interglacial at 4.23 Ma based on an orbital configuration representative of that period (Fig. ), together with a greenhouse gas concentration of 400 ppm. Although previous work has suggested that the relative role of insolation and CO${}_{\mathrm{2}}$ in driving warmer climates varies through time, and that high insolation in one hemisphere may not necessarily lead to elevated temperatures at the same latitudes , palaeoenvironmental proxies for this period are consistent with this being an interval of peak warmth (see “Palaeoenvironmental proxies” below) and fundamentally the ice-sheet responds only to the integrated warming anomaly, not insolation or radiative forcing in isolation.

Temperature and precipitation fields from the Pliocene RCM simulation are used as inputs to the ice-sheet model by calculating the anomalies from a present-day RCM simulation and then adding these to the present-day temperature and precipitation fields used in model initialization and spinup (Figs. –). This approach is preferable to one that uses the RCM fields directly, since the RCM does not simulate the present-day state precisely, when compared to observational data (Fig. ). Since the climate model also uses an ice-sheet topography in which WAIS is already removed, we use a standard lapse rate of 8 K km${}^{-\mathrm{1}}$ to adjust the simulated temperatures based on the elevation difference from our present-day ice geometry. Despite this adjustment, however, it is possible that the air temperatures in these areas may be anomalously warm as a result of ocean–atmosphere heat exchange that would not occur if the removed ice were present in the RCM simulation. This may have an influence of the retreat rate of West Antarctica in our simulations, but since our focus lies mainly on the behaviour of the East Antarctic ice sheet we do not expect this issue to affect our findings. We do not apply any lapse rate adjustment to precipitation values on the basis that precipitation rates are likely to be controlled largely by synoptic patterns (which are captured by the RCM) rather than by a simple elevation-dependent cooling of the near-surface atmosphere. Since the RCM uses an ice-sheet geometry in which WAIS is absent it gives us a sea-surface temperature field for areas of our domain where grounded ice currently persists. However, these ocean temperatures are only used to calculate basal melt beneath ice shelves, not at the bed of grounded ice, and so do not affect the ice sheet unless the modelled grounding line retreats. For areas of grounded ice where we do not have RCM-derived ocean temperatures we prescribe a uniform value of 271.2 K, essentially the sea-water freezing point. This avoids potential errors that could be introduced by interpolating ocean fields landward, but we recognize that using a single, low, value rather than an interpolation may lead to underestimated basal melt in subglacial basins during ice-sheet retreat. Our offline, one-way climate–ice-sheet model coupling is less satisfactory than a two-way coupling, but is computationally easier and allows first-order ice-sheet responses to climate perturbations to be investigated. However, we acknowledge that a fully coupled setup might give more accurate ice-sheet and climate simulations.

Palaeoenvironmental proxies

Local records of palaeo-temperature around the Antarctic continent are relatively sparse, yet those that exist, and can be chronologically constrained with some certainty, offer a means by which climate model simulations can be evaluated. We compiled marine paleoenvironmental data from several circum-Antarctic locations that include drill cores from the western Ross Sea (CIROS-2, DVDP-10, AND-2A), Prydz Bay (ODP 1165), Kerguelen Plateau (ODP 751A), Wilkes Land (IODP U1361), and northwestern Antarctic Peninsula (ODP 1096) and geological outcrop at the Prydz Bay coast (Fig. ). Data were obtained from stratigraphic intervals that span 4.5 to 4.0 Ma (C3n.1r to lowermost C2Ar) which is characterized by relatively depleted values ($<$ 3.2 ‰) in the benthic ${\mathit{\delta }}^{\mathrm{18}}$O stack (i.e. relatively warm glacials and interglacials; MIS CN1 and CN5 are two of the warmest interglacials in the early Pliocene). Sea water temperature estimates through this interval are based on a range of proxies including biological analogues and geochemical data that are described in detail below and which are summarized in Table .

Inferred high southern latitude winter, summer, and mean annual sea-surface temperatures from a range of proxy techniques compared to model predictions. Winter temperatures assume sea ice close to the coast but not in open water. Summer temperatures (and thus $T_{\mathrm{mean}}$ also) are minima, and could be higher. Absolute sea-surface temperature predictions from the RCM shown for each site in italics. Although the considerable uncertainties associated with temperature inferences from palaeoecological proxies make precise comparisons impossible, the data generally indicate an underestimate of mean annual sea-surface temperatures in the RCM that averages approximately 1 ${}^{\circ }$C.

Site	Long	Lat	Winter	Summer	$T_{\mathrm{mean}}$	Reference
			(likely)	(min.)	(min.)
AND-1B	167.083333	$-$77.888889	$-$1.8	5	1.6	,
			–2.22	3.03	–0.22
CIROS-2	163.533333	$-$77.683333	$-$1.8	5	1.6	,
			–1.96	3.29	0.07
DVDP-10	163.511667	$-$77.578612	$-$1.8	5	1.6	,
			–1.96	3.29	0.08
ODP-1165	67.218733	$-$64.379583	3	5	4	,
			–0.02	7.19	3.54
IODP-U1361	143.886653	$-$64.409547	3	5	4
			–0.43	6.24	2.87
ODP-1096	$-$76.96376	$-$67.56681	$-$1.8	5	1.6
			–1.87	4.80	1.26
Vestfold Hills	78.132	$-$68.631111	$-$1.8	5	1.6
			–0.99	6.00	2.37

(a) Air temperature anomaly for 4.23 Ma from regional climate modelling, adjusted for surface elevation differences between the RCM input orography and the present-day elevations used to initialize the ice-sheet model. Contours show 2.5 ${}^{\circ }$C increments. Biases used in certain simulations are additional to the warming values shown here. (b) As (a) but showing sea-surface temperature anomaly. Contours show 0.25 ${}^{\circ }$C increments. Locations of profile shown in Fig. also shown (red line), as well as the locations (red squares) of the palaeo-environmental proxy data in Table .

[Figure omitted. See PDF]

Surface water temperatures from biological proxies are derived from diatom, silicoflagellate, and calcareous nannofossil data. Diatom assemblages within the target “warm” interval at all of the drill sites contain, and are often dominated by, open ocean species Thalassionema nitzschioides and Shionodiscus tetraoestrupii . Modern descendants of T. nitzschiodes and S. tetraoestrupii are restricted to areas north of the Polar Front where surface water temperatures are greater than 3.5 and 5.5 ${}^{\circ }$C respectively . Numbers of Azpeitia spp. also increase within the early Pliocene section at IODP Site U1361 and reach a maximum abundance of 13 %. Today Azpeitia tabularis is considered cold-tolerant in Southern Ocean waters but species abundance increases as water temperatures increase and generally exceeds 10 % near the modern subtropical front , which suggests that early Pliocene summer sea surface temperatures at Site U1361 were as high as 10 ${}^{\circ }$C. Importantly, sea-ice-associated diatom taxa are either rare to absent at most of the drill sites or indicate much reduced sea ice coverage . Overall, the diatom assemblages indicate that surface water temperatures at locations south of 60${}^{\circ }$ in the early Pliocene were similar to those in modern sub-Antarctic regions and that the circum-Antarctic ocean was sea-ice-free throughout most or all of the annual cycle.

Relatively warm early Pliocene surface water temperatures implied by the diatom assemblages are supported by silicoflagellate data. In the modern marine environment the silicoflagellate genus Distephanus is dominant south of the Antarctic polar front and the genus Dictyocha is rare or absent. The ratio of these genera in ancient sediment samples can be used to infer past ocean temperature. Dictyocha abundance at Ocean Drilling Program Site 1165 increases between 4.3 and 4.4 Ma, which suggests that mean annual SST increased to approximately 4 ${}^{\circ }$C during this interval . A correlative interval in ODP Sites 748 and 751 contains a similar warming signal .

Although temperature reconstructions based on geochemical proxies are rare, TEX${}_{\mathrm{86}}^L$-derived SST estimates of 5 $\pm$ 4 ${}^{\circ }$C were obtained from diatomite that was deposited between 4.5 and 4.2 Ma at the AND-1B drill site . These data support the temperature reconstructions based on biological proxies and indicate that summer surface water temperatures in the western Ross Sea were 5 to 6 ${}^{\circ }$C warmer than today. Furthermore, sedimentary facies indicate increased sediment input during interglacials Motif 2b of. This is interpreted as representing a warmer glacial regime (compared to the late Pliocene) with more sediment-laden meltwater emanating from the margins of the EAIS and discharging into the Ross Sea.

Other evidence for warmth in the early Pliocene comes from fossil invertebrates and a unique assemblage of vertebrates preserved in the Sørsdal Formation, which crops out at Marine Plain in the Vestfold Hills (Fig. ). The formation comprises up to 7.2 m of friable diatomaceous siltstone and sandstone with dark limestone lenses . The age of the deposit is not well constrained, but diatom data suggest it was deposited between 4.4 and 2.74 Ma, and most likely at the earlier end of this range, possible 4.2 to 4.1 Ma . Cetacean fossils include a species of dolphin, a beaked whale, and baleen whale . These noncryophilic species indicate that waters were ice-free during summer . The Sørsdal Formation also contains the scallop Austrochlamys anderssoni. This and similar species of thick-shelled costate scallop thrived in Antarctic waters during interglacial episodes throughout the Neogene but disappeared during the late Pliocene . Modern descendants of these thick-walled species do not live further south than sub-Antarctic regions where surface water temperatures are warmer than 5.5 ${}^{\circ }$C. Diatoms in the formation also indicate an annual range in SST from $-$1.8 to 5.0 ${}^{\circ }$C () and the absence of coccoliths indicates that summer SSTs were likely no higher than 5 ${}^{\circ }$C ().

In summary, the evidence listed above indicates that summertime surface water temperatures in the early Pliocene were at least as high as 4–5 ${}^{\circ }$C during the warmest interglacial periods, and were likely warmer at 64${}^{\circ }$ S (IODP site U1361). Evidence for minimal sea ice suggests that while winter water temperatures may have approached $-$1.8 ${}^{\circ }$C, summers were sea-ice-free. Based on these data we infer that early Pliocene surface water temperatures in the oceans adjacent to Antarctica's continental ice sheets must have been higher than $-$1.8 ${}^{\circ }$C in the winter and at least 5 ${}^{\circ }$C in the summer – a mean annual temperature that was approximately 3 ${}^{\circ }$C warmer than today. Table  compares these proxy-based interpretations with values simulated by our climate model experiments.

The ice-sheet model

Our ice-sheet simulations use the Parallel Ice Sheet Model (PISM) version 0.6.3, an open-source, three-dimensional, thermodynamic, coupled ice-sheet–ice-shelf model. Both the model and our implementation of it are described in detail elsewhere , so only a summary is given here. In brief, the model combines equations of the shallow-ice and shallow-shelf approximations (SIA and SSA respectively) for grounded ice, and uses the SSA for floating ice . Superposing the SIA and SSA velocity solutions allows basal sliding to be simulated according to the “dragging shelf” approach , and enables a consistent treatment of stress regime across the grounded to floating ice transition . Ice streams develop naturally as a consequence of plastic failure of saturated basal till , depending on the thermal regime and volume of water at the ice-sheet bed. Because we employ a large number of simulations in our study, we adopt a relatively coarse model grid of 20 km. Mesh dependence of results can be an important issue under certain circumstances , but in previous work we have shown that this is not the case in our simulations . Furthermore, the principal aim of our experiments is to identify differences between scenarios, rather than to define absolute ice geometries, and for this purpose we believe our approach is both suitable and robust.

Mass loss in sea-level equivalent metres from the simulated Pliocene Antarctic ice sheet after 10 kyr, using environmental conditions from the regional climate model and augmented with a range of air and ocean temperature bias corrections. Numbers in parentheses represent duplicate experiments in which a more aggressive grounding-line scheme is employed (see “Methods”).

		Sea-surface temperature bias
		0	1	2
Air	0	4 (4.6)	6 (8.7)	7.9 (12.2)
temp.	1	5.6 (6.2)	7.7 (9.4)	8.8 (13)
bias	2	7.7 (8.4)	9.1 (10.9)	10.3 (14.1)

Grounding-line migration is a key component of Antarctic ice-sheet simulations, and is a much-debated modelling challenge. By default we adopt here a novel grounding-line scheme that uses a sub-grid interpolation method to more accurately track migrations . In this scheme we allow the sub-grid interpolation of driving stress at the bed, and the calculation of one-sided derivatives to better characterize the ice-sheet–shelf junction. We run duplicate experiments both with and without an additional interpolation scheme that allows basal melt rates to be smoothly propagated across the grounding line from the first floating cell upglacier to the last grounded cell . This mechanism tends to accelerate grounding-line retreat, leading to rapid mass loss in the initial centuries of the experiments, and greater long-term (near-equilibrium) sea-level-equivalent mass loss from the simulated ice sheets (Table ). It is conceptually supported by geophysical studies of modern-day grounding lines that infer an “estuarine”-type environment at ice-stream grounding zones , and although less catastrophic than the cliff-collapse mechanism used in other models the approach is consistent with recent grounding-line process analyses . Basal melt beneath ice shelves (and at the grounding line) is calculated from a three-equation model that uses temperature, salinity, and pressure to determine the freezing point in the boundary layer . Depending on whether there is melt, freeze-on, or neither, a different approximation of the temperature at the ice-shelf base is used. Outputs from RegCM3 do not include subsurface ocean temperatures, and furthermore, PISM is currently set up to read in just a single ocean temperature field. Consequently we are not able to investigate how vertical differences in ocean warming may impact the ice sheet. Our model does, however, allow the pressure effects of deeper bathymetry to be accounted for, giving rise to enhanced basal melt near grounding lines and less melt in central and outer sectors of the shelf seeFig. S5.

Surface mass balance depends on monthly climatological data and a positive degree-day model that tracks snow thickness and allows for melting of snow and ice at 3 and 8 mm ${}^{\circ }$C${}^{-\mathrm{1}}$ day${}^{-\mathrm{1}}$ respectively. We incorporate a white noise signal (normally distributed, mean zero random temperature increment) into the calculation of daily temperature variations. The standard deviation of daily temperature variability is set at 2 ${}^{\circ }$C, somewhat lower than the commonly employed value of 5 ${}^{\circ }$C, on the basis that the latter has a tendency to overestimate melt . Surface temperatures are adjusted for elevation according to an altitudinal lapse rate of $-$8 K km${}^{-\mathrm{1}}$, and a refreezing coefficient of 0.6 is used to mimic meltwater capture within the snowpack.

Experimental methods

There are two approaches to dealing with ice model parameter uncertainty in the kind of study we present here. One approach undertakes thousands of low-resolution experiments (a large ensemble) with incremental changes in each of several key parameters, such as flow enhancement factors. The results are then subsequently analysed with respect to observational constraints to establish which ensemble members are consistent with the data. We do not adopt this kind of approach. Instead we follow a more targeted methodology in which model parameter choice is incrementally refined through an iterative procedure in which we constrain our model to fit the present-day ice-sheet geometry and surface velocity field. To achieve a good fit we adjust ice flow parameters based on expert judgement, not in an unguided manner as is done with ensemble approaches. We start with an initial 20-year smoothing run that uses only the shallow-ice approximation to derive velocities. This removes spurious surface irregularities. Secondly, we implement a 150 000-year run in which the ice-sheet geometry is held fixed but where internal thermal fields are allowed to evolve. The third phase uses a 25 000-year simulation in which full model physics are employed and the ice sheet is allowed to evolve on all boundaries – that is, it is entirely unconstrained. Through carefully guided parameter iteration our procedure results in a spun-up, thermally and dynamically equilibrated ice-sheet simulation that is the best fit to observational constraints that is possible by tuning available model parameters . Thus although parameter uncertainty can be a large source of error under certain circumstances see, we argue that our approach significantly reduces this uncertainty prior to our undertaking the prognostic experimentation. All experiments start from the same spun-up present-day ice-sheet simulation.

We make an assumption regarding the state of the AIS prior to the interglacial at 4.23 Ma based on geological evidence that indicates ice extent greater than present during early Pliocene glacial episodes , implying that the AIS retreated through a present-day configuration during interglacial transitions. However, since the initial condition of the ice sheet is not known, our assumption of a present-day configuration may be erroneous. However, differences in initial geometry would most likely affect the rate at which new equilibria to the imposed climatologies occurred, rather than the geometries of those steady states. Consequently we consider our equilibrium simulations to be representative of the likely long-term Antarctic response to the prescribed orbital and greenhouse gas configuration. We use outputs from the RCM (described above) to define anomalies to our present-day climate grids . Simulations are run for 10 kyr, which is long enough for most simulations to approach a steady state (see “Results” below). We first model ice-sheet evolution using climate fields simulated by the RCM, using both implementations of the sub-grid grounding-line scheme described above. Then we explore the threshold response of the EAIS by running additional experiments with uniform increments of 1 and 2 ${}^{\circ }$C added to the air and sea-surface temperature fields, which address both a known cold bias in the RCM as well as mismatches between RCM and proxy-based temperature reconstructions. By bracketing a range of temperatures we are also able to investigate more easily the existence and sensitivity of the AIS to climate–ice-sheet thresholds. To isolate the effects of the imposed climatic perturbations most clearly we keep sea level at its modern level, rather than the $+$10 to $+$30 m thought likely for the Pliocene . We note, however, that in separate experiments not shown here the response of the ice sheet appears to be unaffected by sea level changes of these magnitudes. We run duplicates of all experiments using the two variants of the sub-grid grounding-line scheme described above.

Tipping point analysis

The aim of this technique is to analyse time series data (in our case ice mass evolution) to find early warning signs of impending tipping points. Such events are characterized by a non-linear response to an underlying forcing, based on the phenomenon of “critical slowing down”, and relies on the gradual shallowing and widening of the state-space of the system (also known as the “basin of attraction”) . This “slowing down” can be mathematically detected by looking at the pattern of fluctuations in the short-term trends of the data before the threshold is passed . This change in shape of the state-space allows the system to travel further from its point of equilibrium , whereupon the system takes increasingly longer to recover from small perturbations , which can be detected as an increase in the lag-1 autocorrelation and variance of the time series . A change in skewness may also be interpreted as a precursor to tipping .

This tipping point analysis technique has been applied to a range of climate and palaeoclimate data using both natural archives and model outputs . Importantly, the early warning signs, or precursors, are only expected in the presence of two or more quasi-stationary states, separated by an unstable equilibrium, and are not expected in the case of tipping induced by stochastic fluctuations only, so their identification tends to imply a real threshold change. The method only analyses data preceding the tipping point. These data are detrended using a Gaussian kernel smoothing filter over a suitable bandwidth, so that long-term trends are removed without overfitting the data. The resulting residuals are then measured for autocorrelation at lag-1 and variance over sliding windows of two different sizes. The Kendall tau rank correlation coefficient is used to provide a quantitative measure of the trend; this metric varies between $+$1 and $-$1, where higher values indicate a greater concordance of pairs and a stronger increasing trend. Importantly, it is the presence of a parallel increasing trend in both autocorrelation and variance, rather than the absolute values of the indicators, that indicates critical slowing down .

For the purposes of this study we make an important distinction between a tipping point and a threshold, which are sometimes used synonymously. We define a tipping point as a transient feature, whereas a threshold is non-temporal. During the evolution of an ice sheet under a constant forcing it may be that a point is reached in which the trajectory of evolution changes, i.e. the system “tips” into a new state of behaviour. By analysing the time series data from a constant forcing experiment, tipping point analysis is able to not only show where genuine system instabilities occur, but also to provide information on the timescale over which this instability evolves. This is what the final part of our study investigates. However, this point of change, under a steady forcing, is not the same as the identification of a single temperature at which the ice sheet may be stable or unstable in a given area and over a discrete period of time. Consequently it is not possible from our tipping point analysis to provide information on a threshold, as the two phenomena are simply different entities. A detailed explanation of tipping points in Earth systems can be found in , whereas threshold temperatures for individual Antarctic ice-sheet catchments have been quantified in .

Results

Histogram representation on the anomaly data shown in Fig. . Biases used in certain simulations are additional to the warming values shown here. Note the very long tail that skews the mean value of modelled land surface warming, indicating that the modal value might be most representative of the continent as a whole. Sea-surface temperature anomalies are also higher in the open ocean than those of a subset representing ocean areas where present-day ice shelves exist.

[Figure omitted. See PDF]

For the first of our experiments, we run our RCM to simulate the early Pliocene climate at 4.23 Ma. The model yields spatially variable annual air and sea-surface temperature changes whose means are approximately 5 and 2 ${}^{\circ }$C above present respectively, for presently ice-covered areas of Antarctica (Figs.  and ). The advantage of using a known period of the past such as this, rather than a forward projection, is that we can use empirical data to evaluate our modelled climate fields and to constrain our ice-sheet simulations. Proxy data suggest that our simulated climate may be 1–2 ${}^{\circ }$C cooler than actually occurred (Table ); thus, we run duplicate simulations allowing for additional warming of 1 and 2 ${}^{\circ }$C in both the atmosphere and ocean. The sea-level equivalent ice volume loss from Antarctica in these simulations ranges from 4 to 14 m (Table ; Fig. ). Assuming a 5 to 7 m contribution from the Greenland ice sheet and 0.5 to 1 m from thermal expansion of the oceans (assuming 1 to 3 ${}^{\circ }$C ocean warming) , all of our simulations are consistent with a Pliocene sea level highstand range of approximately 10 to $>$ 20 m reconstructed from proxy records . However, by considering likely air–ocean temperature relationships , as well as empirical records of an absent WAIS and a retreated EAIS at this time , we are able to identify two scenarios as most plausible. Figure a and b illustrate the modelled Antarctic ice-sheet geometry that arises under constant forcing with the Pliocene climate, modified to incorporate a $+$2 ${}^{\circ }$C air temperature bias and a 0 or 1 ${}^{\circ }$C bias in the ocean (see “Methods”). In both cases we reproduce the smaller-than-present ice geometry, as well as the pattern of basal ice velocities that would have controlled bedrock erosion and sediment transport, that are necessary to be consistent with geological interpretations .

Simulated Antarctic ice-sheet contributions to global mean sea level under a range of RCM-based climatologies, using two different grounding-line parameterizations. Values from previous studies are also shown (grey bars), as well as global sea-level highstand values inferred from far-field sites.

[Figure omitted. See PDF]

Basal ice velocities for the simulated early Pliocene ice sheet under environmental conditions from regional climate model simulations using two different bias adjustments. Fastest-flowing outlets occur in the subglacial basins and troughs of East Antarctica. These zones correspond to inferred areas of subglacial erosion during past warm climates . Red line identifies catchment transect shown in Fig. .

[Figure omitted. See PDF]

In both of the scenarios shown in Fig. , substantial grounding-line retreat is evident in the WSB, but the Aurora and Recovery basins are less affected. To better understand the timescales and rates involved in retreat in the WSB, we take timeslice ice and bed geometries along the WSB flowline from the two simulations considered representative of Pliocene conditions at 4.23 Ma (Fig. a, b). In both scenarios, we find that margin retreat under constant climate forcing is punctuated by periods of stability. Figure a illustrates a gradual surface lowering at the margin of the ice sheet that continues for 2000 years before flotation and grounding-line retreat ensue. At this point, rapid retreat takes place across the deepest part of the WSB, and grounded ice is replaced by a floating ice shelf. Peak retreat rates are sustained for two to three centuries during this period (Fig. c) before declining. Continued retreat across the inner portion of the WSB proceeds more slowly, but is similarly punctuated by alternating episodes of acceleration and relative stability that continue for 7 to 8 kyr (Fig. c). Under warmer oceanic conditions (Fig. b) retreat proceeds more quickly, but initial destabilization from the seaward pinning point still requires thinning and lowering of the ice-sheet margin, which in this case takes place over the first 1000 years.

Mechanism and timescale of ice-sheet retreat across Wilkes Subglacial Basin. (a) Ice-sheet margin initially occupies stable location pinned on bedrock high, unaffected by marine influence. Gradual surface lowering from rising air temperatures promotes thinning, flotation, and subsequent rapid inland retreat of margin. (b) Identical simulation to (a) but employing a warmer ocean. Coloured lines denote ice geometries at 200-year intervals. Note also the associated bedrock uplift following ice retreat. (c) Rate of retreat across the WSB is non-uniform, despite time-invariant forcings, and is governed primarily by the location of bedrock pinning points.

[Figure omitted. See PDF]

Figure  illustrates that thresholds clearly exist in areas like the Wilkes subglacial basin where the topography beneath the ice sheet allows for rapid retreat after prolonged thinning under a warmer climate destabilizes the grounded margin. However, in order to robustly assess whether genuine “tipping points” exist in the Antarctic ice-sheet system it is necessary to use a statistical treatment of the time series data in which statistical signatures of instability are sought. To do this we first plot sea-level equivalent whole-continent ice volume trends (Fig. a, h), using air and ocean temperature bias corrections as described above. Viewing ice volume evolution of these six simulations in semi-log space it is evident that the modelled ice sheets evolve through three distinct phases. First, in all experiments there is a period of rapid mass loss in which the marine-based sectors of the WAIS collapse. This phase takes one to three centuries, depending on the magnitude of the applied forcing. The second phase is characterized by a slower, but sustained, sea-level contribution from the EAIS, arising from both dynamic adjustment to the loss of fringing ice shelves, and from surface lowering resulting from negative mass balance in coastal areas. It is during this phase that trajectories begin to diverge, depending on the air temperature forcing applied. Using the unadjusted RCM air temperatures, the EAIS stabilizes and actually accumulates sufficient mass from the associated increase in precipitation to offset domain-integrated losses. In this scenario, the ice sheet effectively removes water from global oceans and leads to a slight lowering of sea level (Fig. a, blue line). Perturbing the climate with an additional one or two degrees of atmospheric warming, however, means that surface melt is no longer completely offset by increased precipitation, and the sea-level contribution from the ice sheet continues to increase. The third and final phase of ice-sheet evolution in our simulations takes place once ice in major subglacial basins begins to retreat. Interestingly, in terms of ice-sheet thresholds the key tipping point identified in the Wilkes Subglacial Basin appears to be triggered in very similar ways by both atmospheric and oceanic forcing (Fig. ), suggesting that several climate scenarios may produce equilibrium ice-sheet geometries that are consistent with the empirical constraints.

(a) Sea-level equivalent ice mass loss for ice-sheet simulations forced by regional climate model outputs allowing for 0–2 ${}^{\circ }$C air temperature bias. Inset box identifies region shown in (b, c). Precursor signatures of tipping points: (b) increasing autocorrelation (AR1), and (c) increasing variance (Var). Grey lines show precursors over 1500-year sliding window, black lines show precursors during the 250-year sliding window immediately preceding the tipping point. Additional detail shown in Fig. . Mean annual thinning rate over precursor (d) and response (e) periods. (f, g) Thinning rate (dHdt) and surface velocity (Vel.) at locations in the Recovery and Wilkes subglacial basins (white boxes in d, e). Note the lagged response of peak ice velocity compared to thinning maxima. Panels (h–n) as (a–g) but for scenarios with 1 ${}^{\circ }$C sea-surface temperature bias.

[Figure omitted. See PDF]

Using a sliding window to carry out tipping point analyses of these data we can establish whether or not statistical signatures of instability are evident during any of these three phases. We analysed six ice-sheet trajectories that encompass the range of environmental conditions considered most consistent with RCM simulations and empirical constraints, using air temperature biases of 0, 1, and 2 ${}^{\circ }$C and ocean temperature biases of 0 and 1 ${}^{\circ }$C added to the RCM-simulated climatologies. In these six experiments we found evidence of critical slowing down in three. The loss of WAIS occurs too quickly in our simulations to determine whether tipping points for that ice sheet exist or not. However, in the coolest of the six scenarios, using the unadjusted RCM values, we found evidence for a tipping point that leads to accelerated growth of the East Antarctic ice sheet (Fig. a, blue line). The next warmest scenario yields a stable EAIS (Fig. a, yellow line) with no clear tipping point, and the two warmest scenarios (Fig. h, yellow and red lines) lead to such rapid mass loss that no tipping points are detected. In the two intermediate scenarios (Fig. a, red line; h, blue line), however, two obvious inflections in mass loss are visible. Focusing on these two scenarios, we employ sliding analysis windows of two different widths, 250 and 1500 years, to see if precursor signals indicative of the build-up to a tipping point can be seen in advance of the actual threshold. Figure b, c, i, j illustrate the trends in autocorrelation at lag-1 (AR1) and variance (Var) for the two scenarios, overlain on the corresponding mass loss trajectories. In all cases, analyses using the 250-year window (black lines) reveal abruptly increasing AR1 and Var values (early warning signals of critical slowing down) immediately preceding the tipping point. A change in skewness is also observed (Fig. ), which suggests that the system may have reached a region of asymmetry in its basin of attraction. What is surprising, however, is that even the analyses using the 1500-year sliding window (Fig. b, c, i, j, grey lines) show evidence of increasing trends (Kendall tau values shown in Fig. ), suggesting that the system may be showing signs of instability far in advance of the observed tipping point, and on a timeframe comparable to the gradual thinning observed in the Wilkes subglacial basin analysis described above (Fig. ).

Tipping point analyses for two scenarios. Upper row shows results of analyses of dT $+$2 ${}^{\circ }$C, dSST $+$0 ${}^{\circ }$C, lower row shows dT $+$0 ${}^{\circ }$C, dSST $+$1 ${}^{\circ }$C. In each panel, top row shows modelled mass loss curve with inset identifying period used in the tipping point analysis. The four lower rows each show one key indicator metric (residuals, autocorrelation, variance, skewness) as described in “Methods”. Left and right columns show results based on 1500-year (from the solid line; dashed line shows analysis prior to this window) and 250-year sliding windows respectively. Note the steeper trends and higher tau values in the shorter analyses, indicative of a stronger signal over the shorter timeframe. Inset contour plots display the range of Kendall tau values for different sliding window lengths (25–50 %) and smoothing bandwidths (5–15 %).

[Figure omitted. See PDF]

To determine whether the results are sensitive to parameter choices such as the length of the smoothing bandwidth or sliding window size, we ran repeats of the analysis with a range of 25 smoothing bandwidths (ranging from 5 to 15 % of the time-series length) and 25 sliding window sizes (ranging from 25 to 50 % of the time-series length). The results are visualized using contour plots of the Kendall tau values of these repeats; a more homogeneous colour indicates increased robustness (Fig. ). In order to test the significance of these results we created a surrogate dataset by randomizing the original data over 5000 permutations based on the null hypothesis that the data are generated by a stationary Gaussian linear stochastic process. This method guarantees the same amplitude distribution as the original time series, but removes any ordered structure or linear correlation and makes no assumption about the statistical distribution of the data. The Kendall tau values for autocorrelation at lag-1 and variance were computed for each of the surrogate time series, and the probability of making a type I statistical error (false positive) for the original data was computed by comparing to the probability distribution of the surrogate data (Table  and Fig. ).

Calculated significance statistics for two precursor metrics of tipping points. Kendall tau and $p$ values are shown for autocorrelation (“AR1”) and variance (“Var”) trends over 250 (“short”) and 1500 (“long”) sliding windows, for two scenarios in which distinct inflections occur in their mass loss time series. See Fig.  for more information.

		Kendall tau	$p$ value
Scenario 3 short	AR1	0.805	$<$ 0.05
	Var	0.626	$<$ 0.12
Scenario 3 long	AR1	0.588	$<$ 0.15
	Var	0.584	$<$ 0.15
Scenario 4 short	AR1	0.771	$<$ 0.05
	Var	0.911	$<$ 0.01
Scenario 4 long	AR1	0.656	$<$ 0.1
	Var	0.721	$<$ 0.06

Histograms showing Kendall tau correlation coefficients for surrogate time series, generated from the two scenarios for which the existence of tipping points was investigated. Values relate to autocorrelation (blue lines) and variance (red lines) scores shown in Table . Dashed lines mark 90 and 95 % confidence intervals.

[Figure omitted. See PDF]

To understand the glaciological changes taking place as these thresholds are crossed, we divided the analysis time periods into “precursor” and “response” episodes, and calculated mean thinning rates across the ice sheet (Fig. d, e, k, l). There are clear indicators of ice-sheet thinning in the major subglacial basins of East Antarctica, but when velocity and thinning rate data are extracted from two of the most dynamic areas there appears to be no clear distinction between the two episodes – in one scenario thinning and acceleration seem to occur during both the precursor and the response phases, whereas in the other scenario dynamic thinning takes place only during the response phase (Fig. f, g, m, n). What is clear, however, is that wherever ice-sheet thinning is observed, it is associated with an acceleration of ice flow at the same location. In our examples the lag between peak thinning rates and maximum velocities ranges from 200 to 1200 years.

Discussion

Far-field records of Pliocene sea-level changes are sparse, and their interpretation is complicated by glacio-isostatic effects and dynamic topography . Various techniques used to define an envelope of the likely Pliocene sea-level highstand imply a GMSL contribution from the EAIS, but the uncertainties involved are sufficiently great that the sea-level contribution from East Antarctica remains hard to quantify, even though some loss from this ice sheet seems likely . Novel isotope-enabled climate and ice-sheet modelling has allowed the likely GMSL contribution from Antarctica to be constrained to 3–12 m during the warmest parts of the mid-Pliocene, with an absolute maximum of 13 m , which requires at least some loss of ice from East Antarctica. The 4.23 Ma interglacial in our study was likely warmer than mid-Pliocene interglacials, however, suggesting mid-Pliocene sea-level estimates should be treated as minima in our comparison. Warmer air temperatures during all Pliocene interglacials most likely led to increased precipitation in East Antarctica, thickening the ice-sheet interior . Consequently, any sea-level-equivalent mass loss from the EAIS must have been over and above these mass gains, perhaps supporting geologically based inferences of large-scale ice-sheet margin retreat in areas such as the Wilkes Subglacial Basin .

Our study set out to shed further light on this period of uncertainty by investigating the sensitivity of the Antarctic ice sheet to climate-related thresholds, with a particular focus on a peak-warmth interglacial of the early Pliocene period at 4.23 Ma. To establish the likely climatic conditions of this time we use global and regional climate models, but lean heavily on empirical biogeochemical proxy data to validate these simulations. Because there are considerable uncertainties in both the modelled and proxy-inferred temperatures, we defined an envelope of air and ocean temperatures and ran a small ensemble of ice-sheet simulations to capture the range of possible responses. We find that, under the applied climate conditions, the Antarctic ice sheet contributed 8.6 $\pm$ 2.8 m to global mean sea level in the early Pliocene, consistent with some recent studies of the mid-Pliocene but higher than , and lower than , .

The reduction in AIS volume that we simulate arose from a loss of marine-based portions of WAIS, and an inland migration of the EAIS grounding line into the Wilkes subglacial basin. Lesser retreat occurred in the eastern Weddell Sea, in the Amery Basin, and at the margins of the Aurora subglacial basin. In our simulations it appears that the significant retreat in the WSB arises as a consequence of prolonged surface lowering due to atmospheric warming and increasing surface melt, which after 1000 to 2000 years allows flotation of the present-day ice margin and retreat from the primary pinning point that currently maintains stability of this sector. An increase in surface melting at this time is consistent with facies of this age in the Ross Sea that record increased deposition of terrigenous sediment . Subsequent retreat proceeds rapidly into the inland-deepening basin, but is punctuated by periods of slower retreat where bedrock highs allow temporary stabilization of the ice margin.

The delayed atmospheric warming control on ice-sheet behaviour that occurs in the WSB is also mirrored at the continental scale. Our whole-continent ice volume trajectories show that, over multi-millennial timescales, the magnitude of atmospheric warming is critical in determining whether or not ice-sheet mass loss exhibits tipping points (i.e. non-linear behaviour). We find the clearest evidence of tipping under climate scenarios that are moderate, rather than extreme. This is because of the sensitive interplay between the climate forcing that drives retreat and the topographic restraint afforded by pinning points that aid stability. This balance is overwhelmed under extremely warm climates and ice-sheet retreat proceeds rapidly and linearly.

Conclusions

Recent studies have alluded to the idea that the Antarctic ice sheets may be susceptible to “tipping points” in their stability that could give rise to non-linear behaviour and rapid contributions to sea level . We have critically examined this by simulating the Antarctic climate and ice-sheet system under warmer-than-present conditions of a period of the past when geological evidence indicates a much reduced ice sheet, despite CO${}_{\mathrm{2}}$ levels similar to present. Our results indicate an Antarctic contribution to global sea level of approximately 8.5 m during the 4.23 Ma early Pliocene interglacial, sourced primarily from West Antarctica and the WSB of East Antarctica. These findings are consistent with modern studies that implicate ocean thermal forcing as the dominant driver of marine-based ice-sheet retreat in West Antarctica , but we find in addition that over centennial to millennial timescales, atmospheric warming plays a key role in destabilizing sectors of East Antarctica . Whereas the rate of ice-sheet retreat within EAIS subglacial basins is controlled by the interplay between ocean temperature, local bedrock elevation, and ice thickness, initial destabilization appears to be governed by the prolonged and gradual surface lowering in response to warmer-than-present air temperatures. Under certain circumstances this long period of top-down melting constitutes the “precursor” to the eventual tipping point, beyond which ice loss accelerates and produces a non-linear contribution to sea level. Whether such “early warning” signals could be identified in the short time series of satellite-era records of the modern ice sheet remains an research area yet to be explored.

The Parallel Ice Sheet Model is freely available as open-source code from the PISM github repository (git://github.com/pism/pism.git). RegCM3 is available from https://users.ictp.it/~regcm/model.html. Bedrock topography and ice thickness data are from the BEDMAP2 compilation, available at http://www.antarctica.ac.uk//bas_research/our_research/az/bedmap2/. Information on surface mass balance data is available at http://www.projects.science.uu.nl/iceclimate/models/antarctica.php#racmo21. Air temperature and geothermal heat flux inputs were taken from the ALBMAP v1 compilation and can be downloaded from 10.1594/PANGAEA.734145.

The datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request.

NRG devised and carried out the ice-sheet modelling experiments and led the research. ZAT undertook the tipping point analysis. RHL, TRN, and RMM compiled and interpreted paleoenvironmental proxy data. EGWG and DEK undertook climate modelling experiments. All authors contributed to the development of ideas and writing of the paper.

The authors declare that they have no conflict of interest.

Acknowledgements

We are grateful to the two anonymous reviewers for their comments on a previous version of this manuscript. This work was funded by contract VUW1203 of the Royal Society of New Zealand's Marsden Fund, with support from the Antarctic Research Centre, Victoria University of Wellington, ANDRILL, GNS Science, NSF Polar Programs grants ANT-1043712 and PLR-1245899, and the Australian Research Council (ARC) including a Laureate Fellowship (FL100100195) supporting Zoë A. Thomas. We are very grateful to PISM developers for their continued support. PISM is supported by NASA grants NNX13AM16G and NNX13AK27G.Edited by: Zhengtang Guo Reviewed by: two anonymous referees