ARTICLE doi:10.1038/nature17145 Contribution of Antarctica to past and future sea-level rise Robert M. DeConto1 & David Pollard2 Polar temperatures over the last several million years have, at times, been slightly warmer than today, yet global mean sea level has been 6–9 metres higher as recently as the Last Interglacial (130,000 to 115,000 years ago) and possibly higher during the Pliocene epoch (about three million years ago). In both cases the Antarctic ice sheet has been implicated as the primary contributor, hinting at its future vulnerability. Here we use a model coupling ice sheet and climate dynamics— including previously underappreciated processes linking atmospheric warming with hydrofracturing of buttressing ice shelves and structural collapse of marine-terminating ice cliffs—that is calibrated against Pliocene and Last Interglacial sea-level estimates and applied to future greenhouse gas emission scenarios. Antarctica has the potential to contribute more than a metre of sea-level rise by 2100 and more than 15 metres by 2500, if emissions continue unabated. In this case atmospheric warming will soon become the dominant driver of ice loss, but prolonged ocean warming will delay its recovery for thousands of years. Reconstructions of the global mean sea level (GMSL) during past warm climate intervals including the Pliocene (about three million years ago)1 and late Pleistocene interglacials2–5 imply that the Antarctic ice sheet has considerable sensitivity. Pliocene atmospheric CO2 concentrations were comparable to today’s (∼400 parts per million by volume, p.p.m.v.)6, but some sea-level reconstructions are 10–30 m higher1,7. In addition to the loss of the Greenland Ice Sheet and the West Antarctic Ice Sheet (WAIS)2, these high sea levels require the partial retreat of the East Antarctic Ice Sheet (EAIS), which is further supported by sedimentary evidence from the Antarctic margin8. During the more recent Last Interglacial (LIG, 130,000 to 115,000 years ago), GMSL was 6–9.3 m higher than it is today2–4, at a time when atmospheric CO2 concentrations were below 280 p.p.m.v. (ref. 9) and global mean temperatures were only about 0–2 °C warmer10. This requires a substantial sea-level contribution from Antarctica of 3.6–7.4 m in addition to an estimated 1.5–2 m from Greenland11,12 and around 0.4 m from ocean steric effects10. For both the Pliocene and the LIG, it is difficult to obtain the inferred sea-level values from ice-sheet models used in future projections. Marine ice sheet and ice cliff instabilities Much of the WAIS sits on bedrock hundreds to thousands of metres below sea level (Fig. 1a)13. Today, extensive floating ice shelves in the Ross and Weddell Seas, and smaller ice shelves and ice tongues in the Amundsen and Bellingshausen seas (Fig. 1b) provide buttressing that impedes the seaward flow of ice and stabilizes marine grounding zones (Fig. 2a). Despite their thickness (typically about 1 km near the grounding line to a few hundred metres at the calving front), a warming ocean has the potential to quickly erode ice shelves from below, at rates exceeding 10 m yr−1 °C−1 (ref. 14). Ice-shelf thinning and reduced backstress enhance seaward ice flow, grounding-zone thinning, and retreat (Fig. 2b). Because the flux of ice across the grounding line increases strongly as a function of its thickness15, initial retreat onto a reverse-sloping bed (where the bed deepens and the ice thickens upstream) can trigger a runaway Marine Ice Sheet Instability (MISI; Fig. 2c)15–17. Many WAIS grounding zones sit precariously on the edge of such reverse-sloped beds, but the EAIS also contains deep subglacial basins with reverse-sloping, marine-terminating outlet troughs up to 1,500 m deep (Fig. 1). The ice above floatation in these East Antarctic basins is much thicker than in West Antarctica, with the potential to raise GMSL by around 20 m if the ice in those basins is lost13. Importantly, previous ice-sheet simulations accounting for migrating grounding lines and MISI dynamics have shown the potential for repeated WAIS retreats and readvances over the past few million years18, but could only account for GMSL rises of about 1 m during the LIG and 7 m in the warm Pliocene, which are substantially smaller than geological estimates. So far, the potential for MISI to cause ice-sheet retreat has focused on the role of ocean-driven melting of buttressing ice shelves from below16,18–20. However, it is often overlooked that the major ice shelves in the Ross and Weddell seas and the many smaller shelves and ice tongues buttressing outlet glaciers are also vulnerable to atmospheric warming. Today, summer temperatures approach or just exceed 0 °C on many shelves21, and their flat surfaces near sea level mean that little atmospheric warming would be needed to dramatically increase the areal extent of surface melting and summer rainfall. Meltwater on ice-shelf surfaces causes thinning if it percolates through the shelf to the ocean. If refreezing occurs, the ice is warmed, reducing its viscosity and speeding its flow22. The presence of rain and meltwater can also influence crevassing and calving rates23 (hydrofracturing) as witnessed on the Antarctic Peninsula’s Larson B ice shelf during its sudden break-up in 200224. Similar dynamics could have affected the ice sheet during ancient warm intervals25, and given enough future warming, could eventually affect many ice shelves and ice tongues, including the major buttressing shelves in the Ross and Weddell seas. Another physical mechanism previously underappreciated at the icesheet scale involves the mechanical collapse of ice cliffs in places where marine-terminating ice margins approach 1 km in thickness, with >90 m of vertical exposure above sea level26. Today, most Antarctic outlet glaciers with deep beds approaching a water depth of 1 km are protected by buttressing ice shelves, with gently sloping surfaces at the grounding line (Fig. 2d). However, given enough atmospheric warming above or ocean warming below (Fig. 2e), ice-shelf retreat can outpace its dynamically accelerated seaward flow as buttressing is lost and 1 Department of Geosciences, University of Massachusetts, Amherst, Massachusetts 01003, USA. 2Earth and Environmental Systems Institute, Pennsylvania State University, University Park, Pennsylvania 16802, USA. 3 1 M A RC H 2 0 1 6 VO L 5 3 1 NAT U R E 5 9 1 © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE 2,500 a 2,000 1,500 1,000 500 BS Co as t 0 Si pl e AS Ross Sea –500 Aurora Basin Bedrock elevation (m) Weddell Sea –1,000 Wilkes Basin –1,500 –2,000 –2,500 1,000 b The Antarctic Ice Sheet in the Pliocene 900 800 Larson C ice shelf Ronne ice shelf 700 WAIS Pine Island Thwaites 500 EAIS WDIC Law Dome Totten Ross ice shelf 600 400 Ice speed (m yr–1) Amery ice shelf Recovery 300 200 Ninnis Mertz in the few places where such conditions exist today (the Helheim and Jakobsavn glaciers on Greenland and the Crane Glacier on the Antarctic Peninsula), hinting that a Marine Ice Cliff Instability (MICI) in addition to MISI could be an important contributor to past and future ice-sheet retreat. Our three-dimensional ice sheet–ice shelf model25,27 (Methods) predicts the evolution of continental ice thickness and temperature as a function of ice flow (deformation and sliding) and changes in mass balance via precipitation, runoff, basal melt, oceanic melt under ice shelves and on vertical ice faces, calving, and tidewater ice-cliff failure. The model captures MISI (Fig. 2a–c) by accounting for migrating grounding lines and the buttressing effects of ice shelves with pinning points and side-shear. To capture the dynamics of MICI (Fig. 2d–f), new physical treatments of surface-melt and rainwater-enhanced calving (hydrofracturing) and grounding-line ice-cliff dynamics have been added25. Including these processes was found to increase the model’s contribution to Pliocene GMSL from +7 m (ref. 18) to +17 m (ref. 25). The model formulation used here is similar to that described in ref. 25, but with improvements in the treatment of calving, thermodynamics, and climate–ice–ocean coupling (Methods). 100 0 Figure 1 Antarctic sub-glacial topography and ice sheet features. a, Bedrock elevations13 interpolated onto the 10-km polar stereographic ice-sheet model grid and used in Pliocene, LIG, and future ice-sheet simulations. b, Model surface ice speeds and grounding lines (black lines) show the location of major ice streams, outlet glaciers, and buttressing ice shelves (seaward of grounding lines) relative to the underlying topography in a. Features and place names mentioned in the text are also shown. AS, Amundsen Sea; BS, Bellingshausen Sea; WDIC, WAIS Divide Ice Core. The locations of the Pine Island, Thwaites, Ninnis, Mertz, Totten, and Recovery glaciers are shown. Model ice speeds (b) are shown after equilibration with a modern atmospheric and ocean climatology (see Methods). retreating grounding lines thicken15. In places where marine-terminating grounding lines are thicker than 800 m or so, this would produce >90 m subaerial cliff faces that would collapse (Fig. 2f) simply because longitudinal stresses at the cliff face would exceed the yield strength (about 1 MPa) of the ice26. More heavily crevassed and damaged ice would reduce the maximum supported cliff heights. If a thick, marine-terminating grounding line began to undergo such mechanical failure, its retreat would continue unabated until temperatures cooled enough to reform a buttressing ice shelf, or the ice margin retreated onto bed elevations too shallow to support the tall, unstable cliffs25. If protective ice shelves were suddenly lost in the vast areas around the Antarctic margin where reverse-sloping bedrock is more than 1,000 m deep (Fig. 1a), exposed grounding-line ice cliffs would quickly succumb to structural failure, as is happening The warm mid-Pliocene and LIG provide complementary targets for model performance, via the ability to produce ∼5–20 m and ∼3.5–7.5 m GMSL from Antarctica, respectively. These two time periods highlight model sensitivities to different processes, because Pliocene summer air temperatures were capable of producing substantial surface meltwater, especially during warm austral summer orbits28. Conversely, LIG temperatures were cooler29, with limited potential for surface meltwater production. Instead, ocean temperatures30 could have been the determining factor in LIG ice retreat31. To simulate Pliocene and LIG ice sheets, we couple the ice model to a high-resolution, atmospheric regional climate model (RCM) adapted to Antarctica and nested within a global climate model (GCM; see Methods). The RCM captures the orographic details of ice shelves and adjacent ice-sheet margins, which is critical here because the new calving and grounding line processes are mechanistically linked to the atmosphere. High-resolution ocean modelling beneath time-evolving ice shelves on palaeoclimate timescales exceeds existing capabilities. Instead, we use a modern ocean climatology32 interpolated to our ice-sheet grid, with uniformly imposed sub-surface ocean warming providing melt rates on sub-ice-shelf and calving-front surfaces exposed to sea water. The RCM climatologies and imposed ocean warming are applied to quasi-equilibrated initial ice-sheet states, with atmospheric temperatures and the precipitation lapse-rate corrected as the ice sheet evolves. As in ref. 25, the Pliocene simulation uses a RCM climatology with 400 p.p.m.v. CO2, a warm austral summer orbit28, and 2 °C imposed ocean warming to represent maximum mid-Pliocene warmth (Extended Data Fig. 1). The model produces an 11.3-m contribution to GMSL rise, reflecting a reduction in its sensitivity of about 6 m relative to the formulation in ref. 25, but within the range of plausible sea-level estimates1,7. Pliocene retreat is triggered by meltwaterinduced hydrofracturing of ice shelves, which relieves backstress and initiates both MISI and MICI retreat into the deepest sectors of WAIS and EAIS marine basins. The Antarctic Ice Sheet during the LIG Summer air temperatures in the RCM are slightly warmer at 116 kyr ago than 128 kyr ago, but remain below freezing in both cases, with little to no surface melt (Extended Data Fig. 2). As a result, substantial oceanic warming >4 °C is required to initiate WAIS retreat at 128 kyr ago, which occurs once an ocean-melt threshold is reached in the stability of the Thwaites grounding line (Extended Data Fig. 3a and d). Allowing two-way coupling between the RCM and the ice-sheet model 5 9 2 NAT U R E VO L 5 3 1 3 1 M A RC H 2 0 1 6 © 2016 Macmillan Publishers Limited. All rights reserved ARTICLE RESEARCH MISI a Ice flux Reve d Ice shelf Grounded ice lope rse-s MICI Buttressing lt Me h bed Ocean Melt Warm deep water Surface meltwater/rain Crevassing, hydrofracturing Buttressing Ice flux elt M Calving Melt Warm deep water Shelf cavity b e Melt Ice flux Grounding-line retreat lt Me Melt c Grounding-line retreat f Melt Calving Melt Unstable cliff-face >90 m >800 m Melt Ice margin retreat Melt M el t Ice flux Grounding-line retreat Ice flux Figure 2 Schematic representation of MISI and MICI and processes included in the ice model. Top-to-bottom sequences (a–c and d–f) show progressive ice retreat into a subglacial basin, triggered by oceanic and atmospheric warming. The pink arrow represents the advection of warm circumpolar deep water (CDW) into the shelf cavity. a, Stable, marineterminating ice-sheet margin, with a buttressing ice shelf. Seaward ice flux is strongly dependent on grounding-line thickness h. Sub-ice melt rates increase with open-ocean warming and warm-water incursions into the ice-shelf cavity. b, Thinning shelves and reduced buttressing increase seaward ice flux, backing the grounding line onto reverse-sloping bedrock. c, Increasing h with landward grounding-line retreat leads to an ongoing increase in ice flow across the grounding line in a positive runaway feedback until the bed slope changes. d, In addition to MISI (a–c), the model physics used here account for surface-meltwater-enhanced calving via hydrofracturing of floating ice (e), providing an additional mechanism for ice-shelf loss and initial grounding-line retreat into deep basins. f, Where oceanic melt and enhanced calving eliminate shelves completely, subaerial cliff faces at the ice margin become structurally unstable where h exceeds 800 m, triggering rapid, unabated MICI retreat into deep basins. (Methods) captures dynamical atmospheric feedbacks as the ice margin retreats. This enhances retreat (Extended Data Fig. 3b, e), but still requires >4 °C of ocean warming to produce a >3.5 m increase in GMSL. We find that by accounting for the additional influence of circum-Antarctic ocean warming on the RCM atmosphere (Methods), the GMSL contribution increases to >6.5 m with just 3 °C sub-surface ocean warming (Extended Data Fig. 3c and f), despite the cooler orbit of the Earth 128 kyr ago. The ocean-driven continental warming at 128 kyr ago agrees with ice core records29 and supports a Southern Ocean control on the timing of ice-sheet retreat30,31, possibly through Northern Hemisphere influences on the ocean meridional overturning circulation33. Alternative simulations (Fig. 3) use time-evolving atmospheric and oceanic climatologies (Methods) based on marine and ice-core proxy reconstructions29. These time-continuous simulations produce GMSL contributions of 6–7.5 m early in the interglacial, followed by a prolonged plateau and rapid recovery of the ice sheet beginning around 115 kyr ago. This result matches the magnitude, temporal pattern, and rate of LIG sea-level change in ref. 3. (Fig. 3a), and the simulated recovery of the WAIS satisfies the presence of ice >70 kyr ago at the bottom of the WAIS Divide Ice Core34. Combined with estimates of Greenland ice loss11,12,35 and ocean thermal effects10, the simulated, Antarctic contributions to Pliocene and LIG sea level are in much better agreement with geological estimates2–4 than previous versions of our model18,27, which lacked these new treatments of meltwater-enhanced calving and ice-margin dynamics, suggesting that the new model is better suited to simulations of future ice response. following three extended Representative Carbon Pathway (RCP) scenarios (RCP2.6, RCP4.5 and RCP8.5)36. Future circum-Antarctic ocean temperatures used in our time-evolving sub-ice melt-rate calculations come from matching, high-resolution (1°) National Center for Atmospheric Research (NCAR) CCSM4 simulations (ref. 37, Extended Data Fig. 5). The simulations begin in 1950 to provide some hindcast spinup, and are run for 550 years to 2500. The RCP scenarios (Fig. 4) produce a wide range of future Antarctic contributions to sea level, with RCP2.6 producing almost no net change by 2100, and only 20 cm by 2500. Conversely, RCP4.5 causes almost complete WAIS collapse within the next five hundred years, primarily owing to the retreat of Thwaites Glacier into the deep WAIS interior. The Siple Coast grounding zone remains stable until late in the simulation, thanks to the persistence of the buttressing Ross Ice Shelf (see Supplementary Video 2). In RCP4.5, GMSL rise is 32 cm by 2100, but subsequent retreat of the WAIS interior, followed by the fringes of the Wilkes Basin and the Totten Glacier/Law Dome sector of the Aurora Basin produces 5 m of GMSL rise by 2500. In RCP8.5, increased precipitation causes an initial, minor gain in total ice mass (Fig. 4d), but rapidly warming summer air temperatures trigger extensive surface meltwater production38 and hydrofracturing of ice shelves by the middle of this century (Extended Data Fig. 4). The Larsen C is one of the first shelves to be lost, about 2055. Around the same time, major thinning and retreat of outlet glaciers commences in the Amundsen Sea Embayment, beginning with Pine Island Glacier (Fig. 4h), and along the Bellingshausen margin. Massive meltwater production on shelf surfaces, and eventually on the flanks of the ice sheet, would quickly overcome the buffering capacity of firn39. In the model, the meltwater accelerates WAIS retreat via its thermomechanical influence on ice rheology (Methods) and the influence of hydrofacturing on crevassing and structural failure of the retreating margin. Antarctica contributes 77 cm of GMSL rise by 2100, and continued loss of the Ross and Weddell Sea ice shelves drives WAIS retreat from three sides simultaneously (the Amundsen, Ross, and Weddell seas), all with Future simulations Using the same model physics and parameter values as used in the Pliocene and LIG simulations, we apply the ice-sheet model to longterm future simulations (Methods). Here, atmospheric forcing is provided by high-resolution RCM simulations (Extended Data Fig. 4) 3 1 M A RC H 2 0 1 6 VO L 5 3 1 NAT U R E 5 9 3 © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE a 10 Time-continuous LIG simulations forced by proxy climatology Δ(GMSL) (m) c 5 b Model (modern initial conditions) Model (glacial initial conditions) Observations 0 −5 130 125 b Maximum retreat 120 Year (kyr ago) c Maximum retreat (modern initial conditions) 115 (glacial initial conditions) +7.54 m +6.09 m 110 d Maximum retreat (old model physics) +1.07 m 4,500 4,000 3,000 2,500 2,000 1,500 Ice thickness (m) 3,500 1,000 500 0 Figure 3 Ice-sheet simulations and Antarctic contributions to GMSL through the LIG driven by a time-evolving, proxy-based atmosphere– ocean climatology. a, Change in GMSL in LIG simulations starting at 130 kyr ago and initialized with a modern ice sheet (blue) or a bigger LGM ice sheet representing glacial conditions at the onset of the LIG (red). A probabilistic reconstruction of Antarctica’s contribution to GMSL is shown in black3 with uncertainties (16th and 84th percentiles) as dashed lines. b, c, Ice-sheet thickness at the time of maximum retreat using modern initial conditions (b) and using glacial initial conditions (c). Ice-free land surfaces are brown. The bigger sea-level response when initialized with the ‘glacial’ ice sheet is caused by deeper bed elevations and the ∼3,000-yr lagged bedrock response to ice retreat50, which enhances bathymetrically sensitive MISI dynamics. d, The same simulation as b without the new model physics accounting for meltwater-enhanced calving or ice-cliff failure27. GMSL contributions are shown at top left. reverse-sloping beds into the deep ice-sheet interior. As a result, WAIS collapses within 250 years. At the same time, steady retreat into the Wilkes and Aurora basins, where the ice above floatation is >2,000 m thick, adds substantially to the rate of sea-level rise, exceeding 4 cm yr−1 (Fig. 4c) in the next century, which is comparable to maximum rates of sea-level rise during the last deglaciation40. At 2500, GMSL rise for the RCP8.5 scenario is 12.3 m. As in our LIG simulations, atmosphere–ice sheet coupling accounting for the warming feedback associated with the retreating ice sheet adds an additional 1.3 m of GMSL to the RCP8.5 scenario (Fig. 4b). The CCSM4 simulations providing the model’s sub-ice-shelf melt rates (Extended Data Fig. 5) underestimate the penetration of warm Circum-Antarctic Deep Water into the Amundsen and Bellingshausen seas observed in recent decades41. As a result, the model fails to capture recent, 21st-century thinning and grounding-line retreat along the southern Antarctic Peninsula42 and the Amundsen Sea Embayment43. Correcting for the ocean-model cool bias along this sector of coastline improves the position of Pine Island and Thwaites grounding lines relative to observations42,43 (Fig. 4h) and increases GMSL rise by 9 cm at 2100 (mainly due to the accelerated retreat of Pine Island Glacier), but the correction has little effect on longer timescales (Extended Data Table 1). Ocean warming is important to the behaviour of individual outlet glaciers early in the simulations, but we find that most of the long-term sea-level rise in RCP4.5 and RCP8.5 scenarios is caused by atmospheric warming and the onset of extensive surface meltwater production, rather than ocean warming as implied by other recent studies44–46. Without atmospheric warming, the magnitude of RCP8.5 ocean warming in CCSM4 is insufficient to cause the major retreat of the WAIS or East Antarctic basins; and even with >3 °C additional warming in the Amundsen and Bellingshausen seas it takes several thousand years for WAIS to retreat via ocean-driven MISI dynamics alone (Extended Data Fig. 6). We note that despite the 10-km grid resolution, the model simulates major ice streams well (Fig. 1), including their internal variability18. However, during drastic subglacial-basin retreat the internal variability is quickly overtaken as grounding lines recede into deep interior catchments (see Supplementary Video 10). Large Ensemble analysis To better utilize Pliocene and LIG geological constraints on model performance, we perform a Large Ensemble analysis (Methods) to explore the uncertainty associated with the primary parameter values controlling (1) relationships between ocean temperature and subice-shelf melt rates, (2) hydrofracturing (crevasse penetration in relation to surface liquid water supply), and (3) maximum rates of marineterminating ice-cliff failure. The combination of Pliocene and LIG sea level targets is ideal, because Pliocene retreat is dominated by processes associated with (2) and (3), while the LIG is dominated by process (1). Both Pliocene and LIG ensembles are run with combinations of widely ranging parameter values associated with the three processes, and the combinations are scored by their ability to simulate target ranges of Pliocene and LIG Antarctic sea-level contributions (Methods). The filtered subsets of parameter values capable of reproducing both targets are then used in ensembles of future RCP scenarios (Extended Data Table 2), providing both an envelope of possible outcomes and an estimate of the model’s parametric uncertainty (Fig. 5). Importantly, the ensemble analysis supports our choice of ‘default’ model parameters used in the nominal Pliocene, LIG, and future simulations (Fig. 4, Extended Data Table 2). The lack of substantial ice-sheet retreat in the optimistic RCP2.6 scenario remains unchanged, but the Large Ensemble analysis substantially increases our RCP4.5 and RCP8.5 2100 sea-level projections to 49 ± 20 cm and 105 ± 30 cm, if higher (>10 m instead of >5 m) Pliocene sea-level targets are used. Adding the ocean temperature correction in the Amundsen and Bellingshausen seas (Fig. 4d and h) further increases the 2100 projections in RCP2.6, RCP4.5 and RCP8.5 to 16 ± 16 cm, 58 ± 28 cm and 114 ± 36 cm, respectively (see Methods and Extended Data Tables 1 and 2). 5 9 4 NAT U R E VO L 5 3 1 3 1 M A RC H 2 0 1 6 © 2016 Macmillan Publishers Limited. All rights reserved a e 10 8 6 4,000 4 3,000 2 2,500 2100 2200 Year 15 Δ(GMSL) (m) 4,500 3,500 2000 b RCP2.6 2500 RCP2.6 RCP4.5 RCP8.5 2300 2400 2500 2,000 1,500 RCP2.6 RCP4.5 RCP8.5 RCP8.5 with ice–climate feedback 10 Ice thickness (m) CO2 equivalent (PAL) ARTICLE RESEARCH 1,000 500 f 5 RCP4.5 2500 0 4,500 4,000 3,500 0 2300 2400 3,000 2500 2,500 Rate of GMSL rise (cm yr–1) TG Interior WAIS 6 WS SC T 4 Wilkes RCP2.6 RCP4.5 RCP8.5 RCP8.5 with ice–climate feedback 2,000 1,500 1,000 Recovery Aurora ASE AP 2 500 g RCP8.5 2500 0 2100 2200 Year 2300 2400 0 4,500 4,000 2500 3,500 3,000 RCP8.5 RCP8.5 with 3 ºC ocean adjustment in AS–BS sector 2,500 0.5 2,000 1,500 0 2000 Year 2050 Ice thickness (m) Δ(GMSL) (m) 1.0 1950 h 2200 Year 8 2000 d 2100 Ice thickness (m) 2000 c 1,000 2100 500 RCP8.5 Amundsen Sea sector of the WAIS 2015 2050 0 2100 Figure 4 Future ice-sheet simulations and Antarctic contributions to GMSL from 1950 to 2500 driven by a high-resolution atmospheric model and 1° NCAR CCSM4 ocean temperatures. a, Equivalent CO2 forcing applied to the simulations, following the RCP emission scenarios in ref. 36, except limited to 8 × PAL (preindustrial atmospheric level, where 1 PAL = 280 p.p.m.v.). b, Antarctic contribution to GMSL. c, Rate of sea-level rise and approximate timing of major retreat and thinning in the Antarctic Peninsula (AP), Amundsen Sea Embayment (ASE) outlet glaciers, AS–BS, Amundsen Sea–Bellingshausen Sea; the Totten (T), Siple 2150 2200 2300 Coast (SC) and Weddell Sea (WS) grounding zones, the deep Thwaites Glacier basin (TG), interior WAIS, the Recovery Glacier, and the deep EAIS basins (Wilkes and Aurora). d, Antarctic contribution to GMSL over the next 100 years for RCP8.5 with and without a +3 °C adjustment in ocean model temperatures in the Amundsen and Bellingshausen seas as shown in Extended Data Fig. 5d. e–g, Ice-sheet snapshots at 2500 in the RCP2.6 (e), RCP4.5 (f) and RCP8.5 (g) scenarios. Ice-free land surfaces are shown in brown. h, Close-ups of the Amundsen Sea sector of WAIS in RCP8.5 with bias-corrected ocean model temperatures. 3 1 M A RC H 2 0 1 6 VO L 5 3 1 NAT U R E 5 9 5 © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE Δ(GMSL) (m) 20 a 15 15.65 ± 2.00 m c RCP2.6 RCP4.5 RCP8.5 15 13.11 ± 3.04 m 10 10 5.69 ± 1.00 m 5 5 0 2000 1.5 Δ(GMSL) (m) 20 RCP2.6 RCP4.5 RCP8.5 2100 2200 2300 2400 0 2500 2000 b 1.0 3.97 ± 1.97 m 0.25 ± 0.23 m 1.5 RCP2.6 RCP4.5 RCP8.5 1.05 ± 0.30 m 0.49 ± 0.20 m 0.5 0.19 ± 0.42 m 2100 2200 2300 2400 2500 d 1.0 RCP2.6 RCP4.5 RCP8.5 0.64 ± 0.49 m 0.5 0.26 ± 0.28 m 0 2000 0.11 ± 0.11 m 2020 2040 2060 2080 0 2100 2000 2020 2040 Year 2060 0.02 ± 0.13 m 2080 2100 Year Figure 5 Large Ensemble model analyses of future Antarctic contributions to GMSL. a, RCP ensembles to 2500. b, RCP ensembles to 2100. Changes in GMSL are shown relative to 2000, although the simulations begin in 1950. Ensemble members use combinations of model parameters (Methods) filtered according to their ability to satisfy two geologic criteria: a Pliocene target of 10–20 m GMSL and a LIG target of 3.6–7.4 m. c and d are the same as a and b, but use a lower Pliocene GMSL target of 5–15 m. Solid lines are ensemble means, and the shaded areas show the standard deviation (1σ) of the ensemble members. The 1σ ranges represent the model’s parametric uncertainty, while the alternate Pliocene targets (a and b versus c and d) illustrate the uncertainty related to poorly constrained Pliocene sea-level targets. Mean values and 1σ uncertainties at 2500 and 2100 are shown. Long-term commitment to elevated sea level because of their strong dependency on bathymetry. Future simulations should include coupling with Earth models that account for these processes. Improved ancient sea-level estimates are also needed to further constrain model physics and to reduce uncertainties in future RCP scenarios (Fig. 5). Despite these limitations, our new model physics are shown to be capable of simulating two very different ancient sea-level events: the LIG, driven primarily by ocean warming and MISI dynamics, and the warmer Pliocene, in which surface meltwater and MICI dynamics are also important. When applied to future scenarios with high greenhouse gas emissions, our palaeo-filtered model ensembles show the potential for Antarctica to contribute >1 m of GMSL rise by the end of this century, and >15 m metres of GMSL rise in the next 500 years. In RCP8.5, the projected onset of major ice-sheet retreat occurs sooner (about 2050), and is substantially faster (>4 cm yr−1 after 2100) and higher (Figs 4 and 5) than implied by other recent studies44,45,49. These differences are mainly due to our addition of model physics linking surface meltwater and ice dynamics via hydrofracturing of buttressing ice shelves and structural failure of marine-terminating ice cliffs. In addition, we use (1) freely evolving grounding-line dynamics that preclude the need for empirically calibrated retreat rates49, (2) highly resolved atmosphere and ocean model components rather than intermediate-complexity climate models45 or simplified climate forcing44, and (3) calibration based on major retreat during warm palaeoclimates rather than recent minor retreat driven by localized ocean forcing. As in these prior studies, we also find that ocean-driven melt is an important driver of grounding-line retreat where warm water is in contact with ice shelves, but in scenarios with high greenhouse gas emissions we find that atmospheric warming soon overtakes the ocean as the dominant driver of Antarctic ice loss. Surface meltwater may lead to the ultimate demise of the major buttressing ice shelves (Supplementary Videos 8 and 9) and extensive grounding-line retreat, but it is the long thermal memory of the ocean that will inhibit the recovery of marinebased ice for thousands of years after greenhouse gas emissions are curtailed. Ocean warming alone may be limited in its potential to trigger massive, widespread ice loss, but the multi-millennial thermal response time of the ocean47 will have a profound influence on the ice sheet’s recovery. In simulations run 5,000 years into the future, we conservatively assume no ocean warming beyond 2300 and simply maintain those ocean temperatures while the atmosphere cools assuming different scenarios of CO2 drawdown beginning in 2500 (Methods). For RCP8.5 and natural CO2 drawdown, GMSL continues to rise until 3500 with a peak of about 20 m, after which the warm ocean inhibits the re-advance of grounding lines into deep marine basins for thousands of years (Extended Data Fig. 7). Even in the moderate RCP4.5 scenario with rapidly declining CO2 after 2500, WAIS is unable to recover until the global ocean cools, implying a multi-millennial commitment to several metres of sea-level rise despite human-engineered CO2 drawdown. Given uncertainties in model initial conditions, simplified hybrid ice dynamics, parameterized sub-ice melt, calving, structural ice-margin failure, and the ancient sea-level estimates used in our Large Ensemble analysis, the rates of ice loss simulated here should not be viewed as actual predictions, but rather as possible envelopes of behaviour (Fig. 5) that include processes not previously considered at the continental scale. These are among the first continental-scale simulations with model physics constrained by ancient sea-level estimates, simultaneously accounting for high-resolution atmosphere–ice sheet coupling and ocean model temperatures. However, several important processes are lacking and should be included in future work. In particular, the model lacks two-way coupling between the ice sheet and the ocean. This is especially relevant for RCP8.5, in which >1 Sv of freshwater and icebergs would be supplied to the Southern Ocean during peak retreat (Extended Data Fig. 8). Rapid calving and ice-margin collapse also implies ice mélange in restricted embayments that could provide buttressing and a negative feedback on retreat. The loss of ice mass would also have a strong effect on relative sea level at the margin owing to gravitational and solid-earth deformation effects48, which could affect MISI and MICI dynamics 5 9 6 NAT U R E VO L 5 3 1 3 1 M A RC H 2 0 1 6 © 2016 Macmillan Publishers Limited. All rights reserved ARTICLE RESEARCH Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper. Received 27 May 2015; accepted 12 January 2016. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. Rovere, A. et al. The Mid-Pliocene sea-level conundrum: glacial isostasy, eustasy and dynamic topography. Earth Planet. Sci. Lett. 387, 27–33 (2014). Dutton, A. et al. Sea-level rise due to polar ice-sheet mass loss during past warm periods. Science 3491, http://dx.doi.org/10.1126/science.aaa4019 (2015). Kopp, R. E., Simons, F. J., Mitrovica, J. X., Maloof, A. C. & Oppenheimer, M. Probabilistic assessment of sea level during the last interglacial stage. Nature 462, 863–868 (2009). O'Leary, M. J., Hearty, P. J., Thompson, W. G., Raymo, M. E. & Mitrovica, J. X. Ice sheet collapse following a prolonged period of stable sea level during the last interglacial. Nature Geosci. 6, 796–800 (2013). Raymo, M. E. & Mitrovica, J. X. Collapse of polar ice sheets during the stage 11 interglacial. Nature 483, 453–456 (2012). Seki, O. et al. Alkenone and boron-based Pliocene pCO2 records. Earth Planet. Sci. Lett. 292, 201–211 (2010). Miller, K. G. et al. High tide of the warm Pliocene: implications of global sea level for Antarctic deglaciation. Geology 40, 407–410 (2012). Cook, C. P. et al. Dynamic behaviour of the East Antarctic ice sheet during Pliocene warmth. Nature Geosci. 6, 765–769 (2013). Lüthi, D. et al. High-resolution carbon dioxide concentration record 650,000– 800,000 years before present. Nature 453, 379–382 (2008). McKay, N. P., Overpeck, J. & Otto-Bliesner, B. The role of ocean thermal expansion in Last Interglacial sea level rise. Geophys. Res. Lett. 38, L14605 (2011). NEEM community members. Eemian interglacial reconstructed from a Greenland folded ice core. Nature 493, 489–494 (2013). Stone, E. J., Lunt, D. J., Annan, J. D. & Hargreaves, J. C. Quantification of the Greenland ice sheet contribution to Last Interglacial sea level rise. Clim. Past 9, 621–639 (2013). Fretwell, P. et al. Bedmap2: improved ice bed, surface and thickness datasets for Antarctica. Cryosphere 7, 375–393 (2013). Shepherd, A., Wingham, D. & Rignot, E. Warm ocean is eroding West Antarctic Ice Sheet. Geophys. Res. Lett. 31, L23402 (2004). Schoof, C. Ice sheet grounding line dynamics: steady states, stability, and hysteresis. J. Geophys. Res. 112, F03S28 (2007). Favier, L. et al. Retreat of Pine Island Glacier controlled by marine ice-sheet instability. Nature Geosci. 7, 874–878 (2014). Mercer, J. H. West Antarctic Ice Sheet and CO2 greenhouse effect—threat of disaster. Nature 271, 321–325 (1978). Pollard, D. & DeConto, R. M. Modeling West Antarctic Ice Sheet growth and collapse through the last 5 million years. Nature 458, 329–332 (2009). Pritchard, H. D. et al. Antarctic ice-sheet loss driven by basal melting of ice shelves. Nature 484, 502–505 (2012). Joughin, I., Smith, B. E. & Medley, B. Marine ice sheet collapse potentially under way for the Thwaites Glacier basin, West Antarctica. Science 344, 735–738 (2014). Tedesco, M. & Monaghan, A. J. An updated Antarctic melt record through 2009 and its linkages to high-latitude and tropical climate variability. Geophys. Res. Lett. 36, L18502 (2009). Phillips, T., Rajraram, H. & Steffen, K. Cryo-hydrologic warming: a potential mechanism for rapid thermal response of ice sheets. Geophys. Res. Lett. 37, L20503 (2010). Nick, F. M., Van der Veen, C. J., Vieli, A. & Benn, D. I. A physically based calving model applied to marine outlet glaciers and implications for the glacier dynamics. J. Glaciol. 56, 781–794 (2010). Banwell, A. F., MacAyeal, D. R. & Sergienko, O. V. Breakup of the Larsen B Ice Shelf triggered by chain reaction drainage of supraglacial lakes. Geophys. Res. Lett. 40, 5872–5876 (2013). Pollard, D., DeConto, R. M. & Alley, R. B. Potential Antarctic Ice Sheet retreat driven by hydrofracturing and ice cliff failure. Earth Planet. Sci. Lett. 412, 112–121 (2015). Bassis, J. N. & Walker, C. C. Upper and lower limits on the stability of calving glaciers from the yield strength envelope of ice. Proc. R. Soc. Lond. A 468, 913–931 (2012). Pollard, D. & DeConto, R. Description of a hybrid ice sheet-shelf model, and application to Antarctica. Geosci. Model Dev. 5, 1273–1295 (2012). 28. DeConto, R. M., Pollard, D. & Kowalewski, D. Modeling Antarctic ice sheet and climate variations during Marine Isotope Stage 31. Glob. Planet. Change 96–97, 181–188 (2012). 29. Capron, E. et al. Temporal and spatial structure of multi-millennial temperature changes at high latitudes during the Last Interglacial. Quat. Sci. Rev. 103, 116–133 (2014). 30. Duplessy, J. C., Roche, D. M. & Kageyama, M. The deep ocean during the Last Interglacial period. Science 316, 89–91 (2007). 31. Overpeck, J. T. et al. Paleoclimatic evidence for future ice-sheet instability and rapid sea-level rise. Science 311, 1747–1750 (2006). 32. Levitus, S. et al. World ocean heat content and thermosteric sea level change (0–2000 m), 1955–2010. Geophys. Res. Lett. 39, http://dx.doi.org/10.1029/ 2012GL051106 (2012). 33. Galaasen, E. V. et al. Rapid reductions in North Atlantic Deep Water during the peak of the Last Interglacial period. Science 343, 1129–1132 (2014). 34. Fudge, T. J. et al. Onset of deglacial warming in West Antarctica driven by local orbital forcing. Nature 500, 440–444 (2013). 35. Koenig, S. J., DeConto, R. M. & Pollard, D. Impact of reduced Arctic sea ice on Greenland ice sheet variability in warmer than present climate. Geophys. Res. Lett. 41, 3933–3942 (2014). 36. Meinshausen, N. et al. The RCP greenhouse gas concentrations and their extensions from 1765 to 2300. Clim. Change 109, 213–241 (2011). 37. Gent, P. R. et al. The Community Climate System Model Version 4. J. Clim. 24, 4973–4991 (2011). 38. Trusel, L. D. et al. Divergent trajectories of Antarctic surface melt under two twenty-first-century climate scenarios. Nature Geosci. 8, 927–932 (2015). 39. Munneke, P. K., Ligtenberg, S. R. M., van den Broeke, M. R. & Vaughan, D. G. Firn air depletion as a precursor of Antarctic ice-shelf collapse. J. Glaciol. 60, 205–214 (2014). 40. Carlson, A. E. & Clark, P. U. Ice sheet sources of sea level rise and freshwater discharge during the last deglaciation. Rev. Geophys. 50, RG4007 (2012). 41. Schmidtko, S., Heywood, K. J., Thompson, A. F. & Aoki, S. Multidecadal warming of Antarctic waters. Science 346, 1227–1231 (2014). 42. Wouters, B. et al. Dynamic thinning of glaciers on the Southern Antarctic Peninsula. Science 348, 899–903 (2015). 43. Rignot, E., Mouginot, J., Morlighem, M., Seroussi, H. & Scheuchl, B. Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, and Kohler glaciers, West Antarctica, from 1992 to 2011. Geophys. Res. Lett. 41, 3502–3509 (2014). 44. Golledge, N. R. et al. The multi-millennial Antarctic commitment to future sea-level rise. Nature 526, 421–425 (2015). 45. Winkelmann, R., Levermann, A., Ridgwell, A. & Caldeira, K. Combustion of available fossil fuel resources sufficient to eliminate the Antarctic Ice Sheet. Sci. Adv. 1, e1500589 (2015). 46. Feldmann, J. & Levermann, A. Collapse of the West Antarctic Ice Sheet after local destabilization of the Amundsen Basin. Proc. Natl Acad. Sci. 112, 14191–14196 (2015). 47. Li, C., von Storch, J.-S. & Marotzke, J. Deep-ocean heat uptake and equilibrium climate response. Clim. Dyn. 40, 1071–1086 (2013). 48. Gomez, N., Pollard, D. & Holland, D. Sea-level feedback lowers projections of future Antarctic Ice-Sheet mass loss. Nature Commun. 6, 8798, http://dx.doi. org/10.1038/ncomms9798 (2015). 49. Ritz, C. et al. Potential sea-level rise from Antarctic ice-sheet instability constrained by observations. Nature 528, 115–118 (2015). 50. Huybrechts, P. The Antarctic Ice Sheet and Environmental Change: a Threedimensional Modelling Study PhD thesis, http://epic.awi.de/1463/, Vrije Univ. Brussel (1992). Supplementary Information is available in the online version of the paper. Acknowledgements We thank C. Shields at NCAR for providing CCSM4 ocean model data. NCAR is sponsored by the NSF. We also thank R. Kopp for providing LIG sea-level data, and R. Alley, A. Dutton, and M. Raymo for discussions. This research was supported by the NSF under awards OCE 1202632 PLIOMAX project and AGS 1203910/1203792. Author Contributions R.M.D. and D.P. conceived the model experiments, developed the models, and wrote the manuscript. Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper. Correspondence and requests for materials should be addressed to R.M.D. (deconto@geo.umass.edu). 3 1 M A RC H 2 0 1 6 VO L 5 3 1 NAT U R E 5 9 7 © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE METHODS Ice sheet–ice shelf model. We use an established ice-sheet model, with hybrid ice dynamics following the formulation described in ref. 27, and an internal condition on ice velocity at the grounding line15 that captures MISI (Fig. 2a–c) by accounting for migrating grounding lines and the buttressing effects of ice shelves with pinning points and side shear. Bedrock deformation under changing ice loads is modelled as an elastic lithospheric plate above local isostatic relaxation. A grid resolution of 10 km is used for all simulations, the finest resolution computationally feasible for long-term continental simulations. The model includes newly added treatments of hydrofracturing and ice cliff failure (Fig. 2d–f) described in ref. 25 and extended here. Basal sliding coefficients are determined by an inverse method51, iteratively matching ice-surface elevations to observations until a quasi-equilibrium is reached. In this case, inverted sliding coefficients are derived from a modern (preindustrial) surface climatology, using the same RCM used in our Pliocene, LIG, and future simulations. In addition to the Pliocene and LIG targets highlighted here, the ice sheet– ice shelf model has been shown capable of simulating: (1) the modern ice sheet, including grounding-line positions, ice thicknesses, velocities, ice streams, and ice shelves (Fig. 1b), (2) the Last Glacial Maximum (LGM) extent27, (3) the timing of post-LGM retreat18, and (4) the ability of the ice sheet to regrow to its modern extent following retreat25. Calving and hydrofracturing. Calving depends on the combined penetration depths of surface and basal crevasses, relative to total ice thickness23,26,52,53. Crevasse depths are parameterized according to the divergence of the ice velocity field52, with an additional contribution depending on the logarithm of ice speed that crudely represents the accumulated strain history (ice damage) along a flow path25. Rapid calving is imposed as ice thickness falls below 200 m for unconfined embayments. The 200-m criterion is decreased in confined embayments according to 200 × max[0, min[1, (α − 40)/20]], where α is the ‘arc to open ocean’ (in degrees), crudely representing the effects of ice mélange in narrow seaways. The unconfined onset thickness of 200 m was increased from its value of 150 m in ref. 25 in order to improve modern Ross and Weddell Sea calving-front locations. A similar dependence on α is imposed for oceanic sub-ice-shelf melt rates, as described below. Surface crevasses are additionally deepened (hydrofractured) as they fill with liquid water, which is assumed to depend on the grid-scale runoff of surface melt and rainfall available after refreezing23,53. The crevasse-depth dependence on surface runoff plus rainfall rate R (in metres per year) has been modified slightly for low R values. The R used in equation (B.6) of ref. 25 is changed to: 0 for R < 1.5 m yr−1 4 × 1.5 × (R − 1.5) for 1.5 m yr−1 < R < 3 m yr−1 R2 for R > 3 m yr−1 (as before) This supposes that minimal hydrofracturing occurs for relatively small R values. The linear segment between 1.5 m yr−1and 3 m yr−1 intersects the R2 parabola as a tangent at R = 3. This modification prevents small amounts of recession in some East Antarctic basins for modern conditions, where small amounts of summer melt and rainfall occur. Structural failure of ice cliffs. To account for structural ice-cliff failure26,54 (MICI in Fig. 2), a wastage rate of ice W is applied locally to the grid cell adjacent to tidewater grounding lines with no floating ice, if the required stresses at the exposed cliff face exceed the yield strength of ice. This condition depends on the subaerial cliff height at the interpolated grounding line relative to the maximum ice thickness that can be supported, modified locally to account for any meltwater-enhanced crevasse penetration (hydrofracturing), and any reductions in crevassing caused by backstress. For dry crevassing at an ice margin with no hydrofracturing and no buttressing (backstress), the maximum exposed cliff height is 100 m, assuming an ice yield strength of 1 MPa25,26. The formulation of W results in a steep ramp in wastage rates of 0–3 km yr−1, where exposed ice cliffs ramp from 80 m to 100 m. The maximum wastage rate of 3 km yr−1 used as our default is conservatively chosen, based on recent observations of the Jakobshavn Isbrae Glacier (up to ∼12 km yr−1) and the Crane Glacier (∼5 km yr−1) following the loss of their ice-buttressing shelves55–57. Other modifications to ice-sheet model physics. The model is modified from ref. 25 to include a more physically based parameterization of the vertical flow of surface mobile liquid water (runoff and rainfall) through moulins and other fracture systems towards the base22,58, which affects the vertical temperature profiles within the ice sheet. Vertical sub-grid-scale columns of liquid water are assumed to exist, through which the water freely drains while exchanging heat by conduction with the surrounding ambient ice that cools and can freeze some or all of the liquid water within the ice interior. We use uniform parameter values everywhere: we set the fractional area of sub-grid columns to overall area to be 0.1, and the horizontal scale of drainage elements to be 10 m (R in ref. 22, used in the calculations of conductive heat exchange with ambient ice). The fractional area includes both large moulins and any downward movement of liquid water in crevasses or cracks of all scales, which would be prevalent in the future melting scenarios investigated here. Offline sensitivity tests show low sensitivity of our model behaviour to these values, but further investigation is warranted. For reasonable numerical behaviour, the horizontal heat exchange needs to be part of the time-implicit vertical diffusive heat solution for ambient ice temperature in the main model. To avoid an iterative procedure in cases where all liquid water is frozen before reaching the bed, a time-explicit calculation of the water penetration is made first, and one of the following measures is applied in the time-implicit ice-temperature step: (1) the conductive heat exchange coefficient at all levels is reduced by a constant factor for the column, so that the liquid penetrates to the lowest layer but no further; and (2) the conductive coefficient is set to zero below the depth of furthest penetration. Both methods give very similar results in idealized single-column tests; method (1) was used for all runs here. In cases with greater surface liquid flux, there is no reduction of coefficients and some water reaches the base. A minor bug fix is corrected in the calculation of vertical velocities within the ice (w′ in ref. 27), which previously did not account for the removal of ice at the base due to oceanic melting. This only affects advection of temperature in ice shelves, and has negligible effects on results. Ice-sheet initial conditions. Ice-sheet initial conditions and basal sliding coefficients are provided by a 100-kyr inverse simulation following the methodology in ref. 51, using mass-balance forcing provided by a bias-corrected RCM climatology and modern observed ocean temperatures (described below). In the inverse procedure, basal sliding coefficients under modern grounded ice are adjusted iteratively to reduce the misfit with observed ice thickness, with grounding-line positions fixed to observed locations. The LIG simulation using ‘glacial’ initial conditions (Fig. 3) uses the same basal sliding coefficients (along with a relatively slippery value for modern ocean beds), but initialized from a previous simulation of the LGM with a prescribed, cold glacial climate representing conditions at ∼20 kyr ago. The total ice volume in the modern and glacial ice sheets is 26.55 × 106 km3 and 32.30 × 106 km3, respectively, equivalent to bedrock-compensated GMSL values of 56.80 m and 62.28 m. Atmospheric coupling. Atmospheric climatologies providing surface massbalance inputs to the ice model are provided by decadal averages of meteorological fields from the RegCM3 RCM59, adapted to Antarctica with a polar stereographic grid and small modifications of model physics for polar regions. The RCM uses a 40-km grid, over a generous domain spanning Antarctica and surrounding oceans, nested within the GENESIS v3 Global Climate Model60,61. The GCM and RCM share the same radiation code62 and orbitally dependent calculations of shortwave insolation, important for the Pliocene and LIG palaeoclimate simulations. Anomaly methods are used to correct a small <2 °C Antarctic cold bias in the RCM: T = Texp + Tobs − Tctl P = Pexp × Pobs / Pctl where T is monthly surface air temperature and P is monthly precipitation. Subscripts ‘exp’, ‘obs’ and ‘ctl’ refer to model experiment, observed modern climatology, and model modern control, respectively. A modern (1950) RCM simulation is used for the model modern control, and the ALBMAP data set63 is used for observed modern climatology. In the climatic correction for the difference between the ice-model surface elevation and the interpolated elevation in the climate model or observational data set27, precipitation is now corrected as well as temperature. As before, air temperature T (in degrees Celsius) is shifted by ∆T = γ∆z, where γ = −0.008 °C m−1 is the lapse rate (that is, the decrease in atmospheric temperature with respect to altitude) and ∆z is the elevation difference. Now, precipitation P is multiplied by a Clausius–Clapeyron-like factor: P × 2∆T / 10 Rates of surface snowfall and rainfall are now consistently multiplied by a factor ρw/ρi ≈ 1.1, where ρw and ρi are the densities of liquid water and of ice respectively. This consistently converts between the units of most climate models and climatological databases (metres of liquid water equivalent per year) and the ice-model surface budget terms (metres of ice equivalent per year). Oceanic sub-ice shelf and calving-face melt rates. Direct coupling of highresolution ocean models and ice sheets remains challenging. For present-day simulations we use a parameterization of sub-ice shelf melt rates, similar to that © 2016 Macmillan Publishers Limited. All rights reserved ARTICLE RESEARCH used by other model groups64. The parameterization27 links oceanic melt rates to the nearest observed (or modelled) ocean temperatures: OM = KTρwC w To − Tf (To − Tf ) ρiL f where To is ocean temperature interpolated from the nearest point in an observational (or ocean model) gridded data set, Tf is the local freezing-point temperature at the depth of the ice base, and Cw is the specific heat of ocean water. The transfer factor KT = 15.77 m yr−1 °C−1 results in a combined coefficient (KTρwCw/ρiLf ) of 0.224 m yr−1 °C−2. The depth dependence on Tf produces higher melt rates at the grounding line, as observed, and the dependence on T0 − Tf is quadratic65. Although spatially coarse observational data sets and standard GCM ocean models fail to capture detailed ocean current systems below ice-shelf cavities, this approach (Extended Data Fig. 6e and f) is preferable to the ad hoc prescription of single temperatures and transfer coefficients along individual sectors of the Antarctic margin as in ref. 27. The effects of confined geography on ocean currents are represented by reducing basal melting depending on the total arc to open ocean α, representing the concavity of the coastline25. The melt rate computed from ocean temperatures as above is multiplied by the factor: max[0, min[1, (α − 20)/20]] This effect, combined with the reduction of thin-ice calving with a similar dependence on α described above, allows ice to expand into interior basins during cool-climate recovery after major retreats of marine-based ice, as presumably occurred many times in West Antarctica over the last several million years66. Melting of vertical ice surfaces in direct contact with ocean water is derived from the oceanic melt rate (OM) of surrounding grid cells, but is increased by a scaling factor of 10, producing more realistic calving front positions and in better agreement with hydrographic melt rate observations and detailed modelling67. Present-day sub-ice shelf and calving-face melt rates described here use the 1° resolution World Ocean Atlas32,68 temperatures at 400-m depth, interpolated to the time-evolving ice model grid and propagated under ice-shelf surfaces using contiguous neighbour iteration to provide To. The depth of 400 m represents typical observed levels of Circum-Antarctic Deep Water, a main source of warm-water incursions into the Amundsen Sea Embayment today69. Pliocene simulation. Our default Pliocene simulation uses the same nested GCM–RCM climatology used in a prior study25, with 400 p.p.m.v. CO2 and a generic warm austral summer orbit28 (Extended Data Fig. 1). Ocean temperatures are increased uniformly by 2 °C everywhere in the Southern Ocean. The resulting Antarctic contribution of 11.3 m GMSL implies >15 m GMSL rise if an additional ∼5 m contribution from Greenland70 and the steric effects of a warm Pliocene ocean are also considered. This result is ∼6 m less than in ref. 25, reflecting a reduction in the sensitivity of the model with the changes described above. LIG simulations. The LIG spans a ∼20-kyr interval with greenhouse-gas atmospheric mixing ratios comparable to the pre-industrial Holocene9. Opportunities for Antarctic ice-sheet retreat within this interval include a peak in the duration of Antarctic summers coeval with a boreal summer insolation maximum at 128 kyr ago, and an Antarctic summer insolation maxima one half-precession cycle later at 116 kyr ago (Extended Data Fig. 2). We target these two orbital time slices because they contrast radiatively long and weak (128 kyr ago) versus short and intense Antarctic summers (116 kyr ago), both of which have been postulated to be important drivers of ice volume on glacial–interglacial timescales71. LIG simulations that include climate–ice sheet feedback asynchronously couple the GCM–RCM and the ice-sheet model. In this case, the nested RCM land (ice) surface boundary conditions are updated at the end of the initial retreat at ice model-year 5000 and the ice-sheet model is rerun using the updated climatology. This improves the representation of ice-climate feedbacks via albedo, ocean surface conditions (sea surface temperatures and sea ice), and dynamical effects of the changing topography on the atmosphere. We find that explicitly including climate–ice feedbacks improves model performance, relative to simple lapse-rate adjustments. LIG simulations (Extended Data Table 1; Extended Data Fig. 3d, e) apply anomaly-corrected RCM mass-balance forcing at each LIG time slice, using the appropriate greenhouse gas9,72 and orbital values73 in the nested GCM–RCM. Ocean temperatures are provided by the World Ocean Atlas data set32, with incremental warming of 1–5 °C applied uniformly over the Southern Ocean grid domain. To allow the RCM atmosphere to respond to a warmer Southern Ocean in addition to applying elevated ocean temperatures to the ice model, we increase the southward ocean-heat convergence in the nested GCM–RCM using the methodology described in ref. 28, effectively warming the Southern Ocean sea surface temperatures by ∼2 °C and reducing sea-ice extent. Accounting for the effect of a warmer Southern Ocean on the overlying atmosphere produces more LIG icesheet retreat for a given ocean warming, improving our model–data fit. With this technique, only 3 °C of assumed sub-surface ocean warming is required to produce >6 m GMSL rise from Antarctica at either LIG orbital time slice, reinforcing the notion of a dominant oceanic control on LIG ice-sheet retreat. The two time-continuous LIG simulations using prescribed climatologies (Fig. 3) use bias-corrected, present-day RCM climatologies with a uniform, time-evolving perturbation derived from the average of Antarctic ice-core climatologies compiled in ref. 29. Southern Ocean temperatures are treated similarly, with World Ocean Atlas temperatures32 increased according to the average of circum-Antarctic LIG anomalies29. Only records from marine drill-cores poleward of 45° S are used in the averages, but we note that there is considerable uncertainty in the proxy sea surface temperature estimates (>2 °C)29. This approach also assumes that the proxy sea surface temperatures reflect changes at sub-surface depths (∼400 m), which is uncertain. The resulting anomalies are applied to the ice sheet model at 130 kyr ago, 125 kyr ago, 120 kyr ago, and 115 kyr ago and the ice-sheet model is run continuously from 130 kyr ago to 115 kyr ago. The pairs of air and ocean temperature perturbations applied at each 5-kyr LIG timestep are 1.97° and 1.70°, 1.41° and 1.51°, 0.83° and 1.09°, −1.57° and 0.31°, respectively. The time-continuous LIG simulations are initialized from either a present-day initial ice state (Fig. 1b), or from a prior Last Glacial Maximum simulation with 5.76 × 106 km3 more ice than today. The latter initial condition may better represent the ice sheet at the onset of the LIG and leads to a greater potential sea-level rise owing to the deeper bed conditions early in the deglaciation, which enhances the bathymetrically sensitive MISI dynamics. The proxy-forced LIG simulation clearly supports a maximum Antarctic contribution to GMSL early in the interglacial period (Fig. 3). However, we note that owing to the demonstrated influence of Southern Ocean temperature on the timing of retreat and the uncertain magnitude and chronology of our imposed forcing29, these results cannot definitively rule out maximum Antarctic retreat at the end of the LIG, as has also been proposed4,74 Future simulations. Because of the new ice-model physics that directly involve the atmosphere via meltwater enhancement of crevassing and calving, highly resolved atmospheric climatologies are needed at spatial resolutions beyond those of most GCMs. However, multi-century RCM simulations are computationally infeasible. To accommodate the need for long but high-resolution climatologies, the nested GCM–RCM is run to equilibrium with 1 × PAL, 2 × PAL, 4 × PAL and 8 × PAL CO2. In the ice-sheet simulations, CO2 follows the extended RCP greenhouse gas emissions36 to the year 2500, and the climate at any time is the average of the two appropriate surrounding RCM solutions, weighted according to the logarithm of the concentration of CO2. The RCM climatologies follow total equivalent CO2, which accounts for all radiatively active trace gases in the RCP timeseries. In RCP8.5, equivalent CO2 forcing exceeds 8 × PAL after 2175, but it is conservatively limited here to a maximum of 8 × PAL (Fig. 4a). A 10-yr lag is imposed in the RCM climatologies to reflect the average offset between sea surface temperatures and surface air temperatures in the equilibrated RCM (with equilibrated sea surface temperatures from the parent GCM) and the transient response of the real ocean’s mixed layer. Ocean temperatures in the RCP scenarios are provided by high-resolution (0.5° atmosphere and 1° ocean) NCAR CCSM437 ocean model output, following the RCP2.6, RCP4.5, and RCP8.5 greenhouse gas emissions scenarios run to 2300. Ocean temperatures beyond the limit of the CCSM4 simulations at 2300 are conservatively maintained at their 2300 values. As with the World Ocean Atlas, water temperatures at 400-m depth (between ocean model z-levels 30 and 31) are used in the parameterization of oceanic sub-ice melt (oceanic melt rate) described above. The CCSM4 underestimates the wind-driven warming of Antarctic Shelf Bottom Water41 in the Amundsen and Bellingshausen seas associated with recent increases in melt rates and grounding-line retreats20,42,43. To account for this, additional warming is added to the Amundsen and Bellingshausen sectors of the continental margin. We find the addition of 3 °C to the CCSM4 ocean temperatures increases melt rates to 25–30 m yr−1 (Extended Data Fig. 5f). While still less than observed, this substantially improves grounding-line positions in the Amundsen Sea (Pine Island Glacier in particular) from 1950 to 2015. When applied to RCP4.5 and RCP8.5, the ocean-bias correction accelerates twenty-first-century WAIS retreat (Fig. 4d, g, h) but is found to have little effect beyond 2100 (Extended Data Table 1). Extended RCP greenhouse gas scenarios36 are available up to 2500, beyond which we assume two different scenarios: (1) natural decay of CO275,76 and no further anthropogenic emissions, or (2) engineered, fast drawdown towards pre-industrial levels with an e-folding time of 100 years. These choices are not intended to be definitive, but serve to illustrate the ice-sheet response to a wide range of possible long-term future forcings. © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE Future high-resolution ocean-model output is not available on multi-millennial timescales. In our long (5,000-year) future simulations (Extended Data Fig. 7), CCSM4 ocean temperatures at 400 m depth are assumed to remain at their 2300 values for thousands of years beyond 2300 (until 7000). This assumption is based on the thermal inertia of the deep ocean (thousands of years)47, its longwave radiative feedback on atmospheric temperatures77, and its relative isolation from surface variations. The response of the intermediate and deep ocean to atmospheric and surface-ocean warming before 2300 is heavily lagged in time, and consequently deep-ocean temperatures would continue to rise long after CO2 levels and surface temperatures began to decline after 250077. However, at some point several thousand years later, intermediate- and deep-ocean waters would start to cool if CO2 levels decay as in Extended Data Fig. 7. The trajectory of these temperatures would vary spatially and depend on details of the ocean circulation. To our knowledge, the state of the ocean as it recovers from a greenhouse gas perturbation over these timescales is largely unknown, as relevant coupled atmosphere–ocean global climate model simulations at the resolution and duration appropriate to our ice model have not been run. Consequently, our assumption of constant 400-m ocean temperatures after 2300, although likely to be conservative beyond 2500, may be questionable for the latter parts of the simulations assuming fast, engineered CO2 drawdown. However, assuming the slow, natural pace of CO2 recovery76, atmospheric concentrations would remain above twice the current level of carbon dioxide (2 × CO2) for thousands of years in the RCP8.5 scenario (Extended Data Fig. 7). Assuming a global temperature sensitivity of ∼3 °C per doubling of CO2, our ocean temperatures applied to the long RCP8.5 scenario are probably conservative over the duration of the simulation. Geologically constrained Large Ensemble analysis of future ice-sheet retreat. To quantify model uncertainty due to poorly known parameter values, ensembles of future RCP scenarios are performed with varying model parameters affecting sub-ice oceanic melt rates, meltwater-enhanced calving (hydrofracturing) and marine-terminating ice cliff failure. Ensemble members use the high-resolution atmospheric and ocean forcing described in the main text and above. Alternative ensembles are run both with and without the bias correction of CCSM4 ocean temperatures in the Amundsen and Bellingshausen Seas. The three parameters and four values used for each are as follows. OCFAC is the coefficient in the parameterization of sub-ice-shelf oceanic melt, which is proportional to the square of the difference between nearby ocean water temperature at 400-m depth, and the pressure-melting point of ice. It corresponds to K in equation (17) of ref. 27. The relationship between proximal ocean conditions and melting at the base of floating ice shelves remains a challenging topic of ongoing research78, and a simple parameterization64 is used here. Ensemble values of OCFAC are 0.1, 1, 3 and 10 times the default value of 0.224 m yr−2 °C−2. CREVLIQ is the coefficient in the parameterization of hydrofracturing due to surface liquid. It replaces the constant 100 in equation (B.6) of ref. 25, and is the additional crevasse depth due to surface melt plus rainfall rate, with a quadratic dependence. This crudely represents the complex relationship between surface water and crevasse propagation, and basic model sensitivity is shown in supplementary figure 7b of ref. 25. Values of CREVLIQ are 0 m, 50 m, 100 m and 150 m per (m yr−1)−2. VCLIF is the maximum rate of horizontal wastage due to ice-cliff structural failure. It replaces the default value of 3,000 (3 km yr−1) in equation (A.4) of ref. 25. Its magnitude is based on observed retreat rates of modern large ice cliffs, and basic model sensitivity is shown in supplementary figure 7a of ref. 25. Values of VCLIF are 0 km yr−1, 1 km yr−1, 3 km yr−1 and 5 km yr−1. Medium-range, default values of OCFAC, CREVLIQ, and VCLIF used in our nominal Pliocene (Extended Data Fig. 1), LIG (Fig. 3), and Future (Fig. 4) simulations are OCFAC = 1 (corresponding to 0.224 m yr−2 °C−2), CREVLIQ = 100 m per (m yr−1)−2, and VCLIF = 3 km yr−1, respectively. Simulations for the Pliocene and LIG scenarios are run with all possible combinations of these parameter values, that is, 64 (=43) runs (Extended Data Table 2). Each run is subject to a pass/fail test that its equivalent GMSL rise falls within the observed ranges for the LIG (3.6–7.4 m) and the Pliocene (10–20 m). The filtered subset of parameter combinations that pass (15 out of 64) are then used in an ensemble of future RCP scenarios. An additional ensemble calculation is performed using the same LIG criteria, but a lower accepted range for Pliocene sea-level rise (5–15 m), to reflect the large uncertainty in Pliocene sea-level reconstructions1 (29 out of 64 passed this test). The mean and 1σ range of each ensemble are shown for the three RCP scenarios in Fig. 5, providing both an envelope of possible outcomes and an estimate of the model’s parametric uncertainty. Two alternative sets of future RCP ensembles are run with the ocean-temperature bias correction in the Amundsen and Bellingshausen seas shown in Extended Data Fig. 5. This increases Antarctica’s GMSL contribution by ∼9 cm over the next century in both RCP8.5 and RCP8.5, but has almost no effect on longer timescales (Extended Data Tables 1, 2). In the RCP2.6 ensemble calibrated against the higher >10 m Pliocene sea-level targets, the ocean-bias correction increases both the ensemble-mean and 1σ standard deviation to 16 ± 16 cm in 2100 and 62 ± 76 cm in 2500 (Extended Data Table 1). The increased variance is caused by three simulations in the RCP2.6 ensemble set, in which the stability of the Thwaites Glacier grounding line is exceeded and the WAIS retreats into the deep interior. Although the ensemble members with bias-corrected ocean temperatures are generally more consistent with observations of recent retreat in the Amundsen–Bellingshausen sector, the validity of the bias correction in the long-term future is unknown. Code availability. Ice sheet and climate model codes, results from Pliocene, LIG, and future simulations, and tabulated ensemble results are freely available from the corresponding author. 51. Pollard, D. & DeConto, R. M. A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica. Cryosphere 6, 953–971 (2012). 52. Benn, D. I., Warren, C. R. & Mottram, R. H. Calving processes and the dynamics of calving glaciers. Earth Sci. Rev. 82, 143–179 (2007). 53. Nick, F. M. et al. Future sea-level rise from Greenland’s main outlet glaciers in a warming climate. Nature 497, 235–238 (2013). 54. Bassis, J. M. & Jacobs, S. Diverse calving patterns linked to glacier geometry. Nature Geosci. 6, 833–836 (2013). 55. Joughin, I. et al. Continued evolution of Jakobshavn Isbrae following its rapid speedup. J. Geophys. Res. 113, F04006 (2008). 56. Joughin, I. et al. Seasonal to decadal scale variations in the surface velocity of Jakobshavn Isbrae, Greenland: observation and model-based analysis. J. Geophys. Res. 117, F02030 (2012). 57. Scambos, T. A., Bohlander, J. A., Shuman, C. A. & Skvarca, P. Glacier acceleration and thinning after ice shelf collapse. Geophys. Res. Lett. 31, L18402 (2004). 58. Phillips, T., Rajraram, H., Colgan, W., Steffen, K. & Abdalati, W. Evaluation of cryo-hydrologic warming as an explanation for increased ice velocities in the wet snow zone, Sermeq Avannarleq, West Greenland. J. Geophys. Res. 118, 1241–1256 (2013). 59. Pal, J. S. et al. Regional climate modeling for the developing world: the ICTP RegCM3 and RegCNET. Bull. Am. Meteorol. Soc. 88, 1395–1409 (2007). 60. Alder, J. R., Hostetler, S. W., Pollard, D. & Schmittner, A. Evaluation of a present-day climate simulation with a new coupled atmosphere-ocean model GENMOM. Geosci. Model Dev. 4, 69–83 (2011). 61. Thompson, S. L. & Pollard, D. Greenland and Antarctic mass balances for present and doubled atmospheric CO2 from the GENESIS Version-2 Global Climate Model. J. Clim. 10, 871–900 (1997). 62. Kiehl, J. T. et al. The National Center for Atmospheric Research Community Climate Model: CCM3. J. Clim. 11, 1131–1149 (1998). 63. Le Brocq, A., Payne, A. J. & Vieli, A. An improved Antarctic dataset for high resolution numerical ice sheet models (ALBMAP v1). Earth Syst. Sci. Data 2, 247–260 (2010). 64. Martin, M. A. et al. The Potsdam Parallel Ice Sheet Model (PISM-PIK)—part 2: dynamic equilibrium simulation of the Antarctic ice sheet. Cryosphere 5, 727–740 (2011). 65. Holland, P. R., Jenkins, A. & Holland, D. The response of ice shelf basal melting to variations in ocean temperature. J. Clim. 21, 2558–2572 (2008). 66. Naish, T. et al. Obliquity-paced, Pliocene West Antarctic Ice Sheet oscillations. Nature 458, 322–328 (2009). 67. Slater, D. A., Nienow, P. W., Cowton, T. R., Goldberg, D. N. & Sole, A. J. Effect of near-terminus subglacial hydrology on tidewater glacier submarine melt rates. Geophys. Res. Lett. 42, 2861–2868 (2015). 68. Locarnini, R. A. et al. World Ocean Atlas 2013 Vol. 1 Temperature 1–40 http://data.nodc.noaa.gov/woa/WOA13/DOC/woa13_vol1.pdf (2013). 69. Dutrieux, P. et al. Strong sensitivity of Pine Island ice-shelf melting to climatic variability. Science 343, 174–178 (2014). 70. Koenig, S. J. et al. Ice sheet model dependency of the simulated Greenland Ice Sheet in the mid-Pliocene. Clim. Past 11, 369–381 (2015). 71. Raymo, M. E. & Huybers, P. Unlocking the mysteries of the ice ages. Nature 451, 284–285 (2008). 72. Spahni, R. et al. Atmospheric methane and nitrous oxide of the late Pleistocene from Antarctic ice cores. Science 310, 1317–1321 (2005). 73. Berger, A. Long-term variations of daily insolation and Quaternary climatic changes. J. Atmos. Sci. 35, 2362–2367 (1978). 74. Hearty, P. J., Hollin, J. T., Neumann, A. C., O'Leary, M. J. & McCulloch, M. Global sea-level fluctuations during the Last Interglaciation (MIS 5e). Quat. Sci. Rev. 26, 2090–2112 (2007). 75. Archer, D. et al. Atmospheric lifetime of fossil fuel carbon dioxide. Annu. Rev. Earth Planet. Sci. 37, 117–134 (2009). 76. Archer, D. & Maier-Reimer, E. Effect of deep-sea sedimentary calcite preservation on atmospheric CO2 concentration. Nature 367, 260–263 (1994). 77. Matthews, H. D. & Caldeira, K. Stabilizing climate requires near-zero emissions. Geophys. Res. Lett. 35, L04705 (2008). 78. Hellmer, H. H., Kauker, F., Timmermann, R., Determann, J. & Rae, J. Twentyfirst-century warming of a large Antarctic ice-shelf cavity by a redirected coastal current. Nature 485, 225–228 (2012). 79. Huybers, P. Early Pleistocene glacial cycles and the integrated summer insolation forcing. Science 313, 508–511 (2006). 80. Holden, P. B. et al. Interhemispheric coupling, the West Antarctic Ice Sheet and warm Antarctic interglacials. Clim. Past 6, 431–443 (2010). © 2016 Macmillan Publishers Limited. All rights reserved ARTICLE RESEARCH Extended Data Figure 1 Warm mid-Pliocene climate and ice-sheet simulation. a, January (warmest monthly mean) difference in 2-m (surface) air temperature simulated by the RCM relative to a preindustrial control simulation with 280 p.p.m.v. CO2 and present-day orbit. The temperature difference is lapse-rate-corrected to account for the change in ice-sheet geometry and surface elevations. The Pliocene simulation uses 400 p.p.m.v. CO2, a warm austral summer orbit, and assumes a retreated WAIS to represent maximum Pliocene warm conditions. b, The Pliocene ice-sheet is shown after 5,000 model years, driven by the RCM climate in a, and assuming 2 °C ocean warming relative to a modern ocean climatology32. In the model formulation used here, maximum Pliocene ice-sheet retreat with default model parameters is equivalent to 11.26 m GMSL, about 6 m less than in ref. 25. © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE Extended Data Figure 2 LIG greenhouse gases, orbital parameters, and RCM climates. a, Greenhouse gas concentrations9,72 converted to radiative forcing shows the LIG interval (light red bar) and the best opportunity for ice-sheet retreat. b, Summer insulation at 70° latitude in both hemispheres73 (red, south; blue, north) and summer duration at 70° S (black)79 shown over the last 150 kyr, and the two orbital time slices (vertical dashed black lines at 128 kyr ago and 116 kyr ago). c, Table showing the greenhouse gas atmospheric mixing ratios (CO2 in parts per million by volume; CH4 and N2O in parts per billion by volume) and orbital parameters (eccentricity, obliquity, precession) used in the GCM–RCM at the LIG time slices (dashed lines 1 and 2 in a and b), respectively. d–f, January (warmest monthly mean) differences in 2-m surface air temperature relative to a preindustrial control simulation at 128 kyr ago (d), 116 kyr ago (e), and the present-day (2015) (f). Simulated austral summer temperatures at 116 kyr ago (e) with relatively highintensity summer insolation is warmer than the long-duration summer orbit at 128 kyr ago (d), but unlike the Pliocene (Extended Data Fig. 1a), neither LIG climatology is as warm as the present day, producing little to no rain or surface melt on ice-shelf surfaces. © 2016 Macmillan Publishers Limited. All rights reserved ARTICLE RESEARCH Extended Data Figure 3 Effect of Southern Ocean warming on Antarctic surface air temperatures and the ice sheet at 128 kyr ago. a–c, January (warmest monthly mean) differences in 2-m surface air temperature at 128 kyr ago, relative to a preindustrial control simulation (top row). GHG, greenhouse gas; SST, sea surface temperature. d, e, Ice-sheet thickness (m) after 5,000 model years, driven by the corresponding climate in a–c. a and d, Without climate–ice sheet coupling (present-day ice extent and surface ocean temperatures in the RCM), and prescribed 5 °C sub-surface ocean warming felt only by the ice sheet. b and e, With asynchronous coupling between the RCM atmosphere and ice sheet, and prescribed 5 °C sub-surface ocean warming felt only by the ice sheet. c and f, With asynchronous coupling between the RCM atmosphere and ice sheet, prescribed 3 °C sub-surface ocean warming felt by the ice sheet, and ∼2 °C surface ocean warming felt by the RCM atmosphere. c shows the locations of East Antarctic ice cores (EDC, EPICA Dome C; V, Vostock; DF, Dome F; EDML, EPICA Dronning Maud Land) indicating warming early in the interglacial29 and previously attributed to WAIS retreat80; this warming is similar to that simulated in c from a combination of ice-sheet retreat and warmer Southern Ocean temperatures, supporting the notion that the timing of LIG retreat was largely driven by far-field ocean influences, rather than local astronomical forcing. © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE Extended Data Figure 4 RCM climates used in future, timecontinuous RCP scenarios and evolving ice-surface melt rates linked to hydrofracturing model physics. a–d, January surface (2-m) air temperatures simulated by the RCM at the present-day (2015) (a), twice the present level of carbon dioxide, 2 × CO2 (b), 4 × CO2 (c), and 8 × CO2 (d) with the retreating ice sheet. The colour scale is the same in all panels. Yellow to red colours indicate temperatures above freezing with the potential for summer rain, and surface meltwater production. e–h, Evolving ice-surface meltwater production (in metres per year) in the time-evolving RCP8.5 ice-sheet simulations, driven by a timecontinuous RCM climatology (Methods) following the RCP8.5 greenhouse gas time series (Fig. 4a). Black lines show the positions of grounding lines and ice-shelf calving fronts at discrete time intervals—e, 2050; f, 2100; g, 2150; and h, 2500—with superposed meltwater production rates. © 2016 Macmillan Publishers Limited. All rights reserved ARTICLE RESEARCH Extended Data Figure 5 NCAR CCSM4 ocean temperatures and oceanic sub-ice-shelf melt rates. a, RCP2.6 ocean warming at 400-m depth, shown as the difference of decadal averages from 1950–1960 to 2290–2300. b, Same as a but for RCP4.5. c, Same as a but for RCP8.5. d, CCSM4 RCP8.5 ocean warming from 1950–1960 to 2010–2020 showing little to no warming in the Amundsen and Bellingshausen seas. The red line shows the area of imposed, additional ocean warming. e, f, Oceanic melt rates at 2015 calculated by the ice-sheet model from interpolated CCSM4 temperatures (e), and with +3 °C adjustment in the Amundsen and Bellingshausen seas (f), corresponding to the area within the red line in d. © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE Extended Data Figure 6 Effect of future ocean warming only. a, Antarctic contribution to future GMSL rise in long, 5,000-yr ice-sheet simulations driven by sub-surface ocean warming simulated by CCSM4, following RCP8.5 (black line), with a +3 °C adjustment in the Amundsen and Bellingshausen seas (blue line; see Extended Data Fig. 5) and a warmer +5 °C adjustment (red line). Atmospheric temperatures and precipitation are maintained at their present values. b–d, Ice-sheet thickness at model-year 5,000, driven by sub-surface ocean forcing from CCSM4 (b) and from CCSM4 with a +3 °C (c) or +5 °C (d) adjustment in the Amundsen and Bellingshausen seas. Note the near-complete loss of ice shelves, but modest grounding-line retreat in b, the retreat of Pine Island Glacier in c, and the near-complete collapse of WAIS once a stability threshold in the Thwaites Glacier grounding line is reached in d. © 2016 Macmillan Publishers Limited. All rights reserved ARTICLE RESEARCH Extended Data Figure 7 The long-term future of the ice sheet and GMSL over the next 5,000 years following RCP8.5 and RCP4.5. a, Equivalent CO2 forcing following RCP8.5 until the year 2500, and then assuming zero emissions after 2500 and allowing a natural relaxation of greenhouse gas levels (red) or assuming a fast, engineered drawdown (blue) with an e-folding timescale of 100 years. b. Antarctic contribution to GMSL over the next 5,000 years, following the greenhouse gas scenarios in a. c, The same as a, except showing long-term RCP4.5 greenhouse gas forcing. d, The same as b, except following the RCP4.5 scenarios in c. The insets in b and d show the ice sheet (and remaining sea-level rise) after 5,000 model years in RCP8.5 and RCP4.5, respectively, assuming fast CO2 drawdown (blue lines), highlighting the multi-millennial commitment to a loss of marine-based Antarctic ice, even in the moderate RCP4.5 scenario. Note the different y-axes in RCP8.5 versus RCP4.5. © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE Extended Data Figure 8 Freshwater input to the Southern Ocean. Total freshwater and iceberg flux from 1950 to 2500, following the future RCP scenarios shown in Fig. 4b. Freshwater input calculations include contributions from ice loss above and below sea level and exceed 1 Sv in RCP8.5. © 2016 Macmillan Publishers Limited. All rights reserved ARTICLE RESEARCH Extended Data Table 1 Summary of Antarctic contributions to GMSL during the Pliocene, LIG, future centuries, and future millennia Antarctic contributions to GMSL for the Pliocene and LIG simulations (rows 1–9) with +2 °C ocean warming in the Pliocene and incrementally imposed ocean warming in the LIG simulations. Values shown represent ice retreat at the end of quasi-equilibrated 5000-yr simulations. Time-continuous LIG simulations forced by proxy-based atmosphere and ocean climatologies (rows 10–12) list maximum GMSL contributions occurring early in the LIG (Fig. 3a). The remaining rows list Antarctic contributions to GMSL at specific times (years as shown) in time-dependent future simulations. Ensemble means and standard deviations (1σ) of the RCP ensemble members listed in Extended Data Table 2 are also shown. Future GMSL contributions are shown relative to 2000. © 2016 Macmillan Publishers Limited. All rights reserved RESEARCH ARTICLE Extended Data Table 2 Ensemble simulations of Pliocene, LIG, and future Antarctic contributions to GMSL Varying combinations of three model parameters (first column) correspond to OCFAC, CREVLIQ and VCLIF, respectively (see Methods). The resulting GMSL contributions of each ensemble member driven by Pliocene and LIG climatologies are shown in the second and third columns. Combinations of model parameters satisfying Pliocene and LIG sea-level targets are assigned a Large Ensemble number (LE#) in the fourth column. Default model parameter values (LE# 12) and resulting Pliocene and LIG GMSL rise are in bold type. Four future ensembles using alternative sets of the palaeofiltered Large Ensemble members and following RCP2.5, RCP4.5 and RCP8.5 emissions scenarios are shown at right. The top two ensembles use 29 combinations of parameter values that satisfy LIG sea-level targets and a range of Pliocene sea-level targets between 5 m and 15 m. The bottom two ensembles use a more restricted set of 15 parameter combinations that satisfy a higher Pliocene target range >10 m. Future RCP ensembles at left correspond to the GMSL time series in Fig. 5. The two ensemble sets at far right include the ocean-temperature bias correction described in the text, Fig. 4 and Extended Data Fig. 5. Antarctic GMSL contributions for each ensemble member are shown at 2100 and 2500. Ensemble means and 1σ standard deviations are also shown. GMSL contributions in future ensembles are relative to 2000. © 2016 Macmillan Publishers Limited. All rights reserved