Tamsin L Edwards

Sophie M J Nowicki

Ben Marzeion

Regine Hock

Heiko Goelzer

Hélène Seroussi

Nicolas C Jourdain

Donald A Slater

Fiona E Turner

Christopher J Smith

Christine M Mckenna

Erika Simon

Ayako Abe-Ouchi

Jonathan M Gregory

Eric Larour

William H Lipscomb

Antony J Payne

Andrew Shepherd

Cécile Agosta

Patrick Alexander

Torsten Albrecht

Brian Anderson

Xylar Asay-Davis

Andy Aschwanden

Alice Barthel

Andrew Bliss

Reinhard Calov

Christopher Chambers

Nicolas Champollion

Youngmin Choi

Richard Cullather

Joshua Cuzzone

Christophe Dumas

Denis Felikson

Xavier Fettweis

Koji Fujita

Benjamin K Galton-Fenzi

Rupert Gladstone

Nicholas R Golledge

Ralf Greve

Tore Hattermann

Matthew J Hoffman

Angelika Humbert

Matthias Huss

Philippe Huybrechts

Walter Immerzeel

Thomas Kleiner

Philip Kraaijenbrink

Sébastien Le Clec'h 36

Victoria Lee

Gunter R Leguy

Christopher M Little

Daniel P Lowry

Jan-Hendrik Malles

Daniel F Martin

Fabien Maussion

Mathieu Morlighem

James F O'neill

Isabel Nias

Frank Pattyn

Tyler Pelle

Stephen Price

Aurélien Quiquet

Valentina Radić

Ronja Reese

David R Rounce

Martin Rückamp

Akiko Sakai

Courtney Shafer

Nicole-Jeanne Schlegel

Sarah Shannon

Robin S Smith

Fiammetta Straneo

Sainan Sun

Lev Tarasov

Luke D Trusel

Jonas Van Breedam

Roderik Van De Wal

Michiel Van Den Broeke

Ricarda Winkelmann

Harry Zekollari

Chen Zhao

Tong Zhang

Thomas Zwinger

Projected land ice contributions to twenty-first-century sea level rise

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

competing processes of increasing ice loss and snowfall accumulation in a warming climate. However, under risk-averse (pessimistic) assumptions, Antarctic ice loss could be five times higher, increasing the median land ice contribution to 42 cm SLE under current policies and pledges, with the upper end (95 th percentile) exceeding half a metre even under 1.5°C warming. This would severely limit the possibility of mitigating future coastal flooding. Given this large range (13 cm main projections under 1.5°C warming; 42 cm risk-averse projections under current pledges), adaptation must plan for a factor of three uncertainty in the land ice contribution to 21 st century sea level rise until climate policies and the Antarctic response are further constrained.

Land ice has contributed around half of all sea level rise since 1993, and this fraction is expected to increase [START_REF] Oppenheimer | IPCC Special Report on the Ocean and Cryosphere in a Changing Climate[END_REF] . The Ice Sheet Model Intercomparison Project (ISMIP6 [START_REF] Nowicki | Ice Sheet Model Intercomparison Project (ISMIP6) contribution to CMIP6[END_REF][START_REF] Nowicki | Experimental protocol for sea level projections from ISMIP6 standalone ice sheet models[END_REF] ) for CMIP6 [START_REF] Eyring | Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization[END_REF] and the Glacier Model Intercomparison Project (GlacierMIP 5 ) provide the Intergovernmental Panel on Climate Change (IPCC) with projections of Earth's ice sheet and glacier contributions to future sea level. Both projects use suites of numerical models [START_REF] Goelzer | Design and results of the ice sheet model initialisation experiments initMIP-Greenland: an ISMIP6 intercomparison[END_REF][START_REF] Seroussi | initMIP-Antarctica: an ice sheet model initialization experiment of ISMIP6[END_REF][START_REF] Marzeion | Partitioning the Uncertainty of Ensemble Projections of Global Glacier Mass Change[END_REF] and greenhouse gas emission scenarios [START_REF] Van Vuuren | The representative concentration pathways: an overview[END_REF] as the basis of their projections, and a variety of treatments are considered for the interaction between the ice sheets and the ocean [START_REF] Slater | Estimating Greenland tidewater glacier retreat driven by submarine melting[END_REF][START_REF] Slater | Twenty-first century ocean forcing of the Greenland ice sheet for modelling of sea level contribution[END_REF][START_REF] Favier | Assessment of sub-shelf melting parameterisations using the oceanice-sheet coupled model NEMO(v3.6)-Elmer/Ice(v8.3)[END_REF][START_REF] Jourdain | A protocol for calculating basal melt rates in the ISMIP6 Antarctic ice sheet projections[END_REF] . In total, the projects provide 256 simulations of the Greenland ice sheet, 344 simulations of the Antarctic ice sheet, and 288 simulations of the global glacier response to climate change [START_REF] Marzeion | Partitioning the Uncertainty of Ensemble Projections of Global Glacier Mass Change[END_REF][START_REF] Goelzer | The future sea-level contribution of the Greenland ice sheet: a multimodel ensemble study of ISMIP6[END_REF][START_REF] Seroussi | ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century[END_REF][START_REF] Nowicki | Contrasting contributions to future sea level under CMIP5 and CMIP6 scenarios from the Greenland and Antarctic ice sheets[END_REF] (see also Extended Data Table 1). Although these simulations represent an unprecedented effort [START_REF] Nowicki | Experimental protocol for sea level projections from ISMIP6 standalone ice sheet models[END_REF][START_REF] Goelzer | Design and results of the ice sheet model initialisation experiments initMIP-Greenland: an ISMIP6 intercomparison[END_REF][START_REF] Seroussi | initMIP-Antarctica: an ice sheet model initialization experiment of ISMIP6[END_REF][START_REF] Marzeion | Partitioning the Uncertainty of Ensemble Projections of Global Glacier Mass Change[END_REF][START_REF] Slater | Estimating Greenland tidewater glacier retreat driven by submarine melting[END_REF][START_REF] Slater | Twenty-first century ocean forcing of the Greenland ice sheet for modelling of sea level contribution[END_REF][START_REF] Favier | Assessment of sub-shelf melting parameterisations using the oceanice-sheet coupled model NEMO(v3.6)-Elmer/Ice(v8.3)[END_REF][START_REF] Jourdain | A protocol for calculating basal melt rates in the ISMIP6 Antarctic ice sheet projections[END_REF][START_REF] Goelzer | The future sea-level contribution of the Greenland ice sheet: a multimodel ensemble study of ISMIP6[END_REF][START_REF] Seroussi | ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century[END_REF][START_REF] Nowicki | Contrasting contributions to future sea level under CMIP5 and CMIP6 scenarios from the Greenland and Antarctic ice sheets[END_REF][START_REF] Goelzer | Remapping of Greenland ice sheet surface mass balance anomalies for large ensemble sea-level change projections[END_REF][START_REF] Barthel | CMIP5 model selection for ISMIP6 ice sheet model forcing: Greenland and Antarctica[END_REF] , their computational expense and complexity has meant that they (i) focus mainly on previous generation emissions scenarios (Representation Concentration Pathways [START_REF] Van Vuuren | The representative concentration pathways: an overview[END_REF] , RCPs) developed for the IPCC's Fifth Assessment Report, not the more diverse and policy-relevant Shared Socioeconomic Pathways (SSPs [START_REF] Riahi | The Shared Socioeconomic Pathways and their energy, land use, and greenhouse gas emissions implications: An overview[END_REF]20 ) that underpin the IPCC's Sixth Assessment Report, (ii) are driven mostly by a relatively small number of older generation global climate models developed before CMIP6 [START_REF] Taylor | An Overview of CMIP5 and the Experiment Design[END_REF] , and (iii) have incomplete and limited ensemble designs.

To address these limitations, we emulate the future sea level contribution of the 23 regions comprising the world's land ice (see Extended Data Table 2) as a function of global mean surface air temperature change and as a consequence of marine-terminating glacier retreat in Greenland and ice-shelf basal melting and collapse in Antarctica. The ensembles of ice sheet and glacier models are emulated all at once for each region, using their simulations as multiple estimates of sea level contribution for a given set of uncertain input values, and we incorporate the ensemble spread through the use of a 'nugget 'term in Gaussian Process emulation [START_REF] Andrianakis | The effect of the nugget on Gaussian process emulators of computer models[END_REF][START_REF] Gramacy | Cases for the nugget in modeling computer experiments[END_REF] . Gaussian Process regression requires minimal assumptions about the functional form, and provides uncertainty estimates for the emulator predictions [START_REF] Edwards | Revisiting Antarctic ice loss due to marine ice cliff instability[END_REF] ; most previous emulator-type approaches for sea level rise use parametric models, where the functional form is assumed [START_REF] Church | Climate Change 2013: The Physical Science Basis[END_REF][START_REF] Levermann | Projecting Antarctica's contribution to future sea level rise from basal ice shelf melt using linear response functions of 16 ice sheet models (LARMIP-2)[END_REF][START_REF] Bulthuis | Uncertainty quantification of the multi-centennial response of the Antarctic ice sheet to climate change[END_REF][START_REF] Nauels | Synthesizing long-term sea level rise projections -the MAGICC sea level model v2.0[END_REF][START_REF] Palmer | Exploring the drivers of global and local sea-level change over the 21st century and beyond[END_REF] . We then use the emulators to make probabilistic projections for the glacier and ice sheet sea level contributions under five SSPs and under an additional scenario reflecting current climate pledges (Nationally Determined Contributions, NDCs) [START_REF] Mckenna | Stringent mitigation substantially reduces risk of unprecedented near-term warming rates[END_REF] made under the Paris Agreement. Most projections presented are for the year 2100, but we also estimate a full timeseries by emulating each year from 2016 to 2100. The details of our emulation approach are described in the Methods.

Response to temperature and parameters

Most land ice regions show a fairly linear relationship of increasing mass loss with global mean surface air temperature. Figure 1 shows the temperature-dependence of the sea level contribution at 2100 for the ice sheets and peripheral glaciers (Fig. 1 a-f) and eleven other glacier regions: four with large maximum contributions (Alaska, Arctic Canada North and South, Russian Arctic: Fig. 1g-j), two with non-linear temperature-dependence, giving near or total disappearance at high temperatures (Central Europe and Caucasus: Fig. 1k,l), and the three regions comprising High Mountain Asia (Fig. 1m-o), which are important for local water supply [START_REF] Biemans | Importance of snow and glacier meltwater for agriculture on the Indo-Gangetic Plain[END_REF] . Values of ice sheet parameters are fixed at two possible values for Greenland glacier retreat and Antarctic basal melting, with no Antarctic ice shelf collapse; only simulations using these values are shown. The ensemble designs are not complete -for example, many fewer ice sheet simulations were performed under RCP2.6 than RCP8.5 -so some of the apparent patterns in the simulation data are artefacts of the gaps, which the emulator is intended to account for. Greenland and the glaciers, which are dominated by surface melting [START_REF] Marzeion | Partitioning the Uncertainty of Ensemble Projections of Global Glacier Mass Change[END_REF][START_REF] Goelzer | The future sea-level contribution of the Greenland ice sheet: a multimodel ensemble study of ISMIP6[END_REF][START_REF] Nowicki | Contrasting contributions to future sea level under CMIP5 and CMIP6 scenarios from the Greenland and Antarctic ice sheets[END_REF] , show clear dependence on temperature. Fourteen of the nineteen glacier regions show approximately linear relationships, and five are nonlinear (Fig. 1f,k,l; also Western Canada & U.S. and North Asia, which have weaker nonlinearity: not shown). In contrast, East Antarctica (Fig. 1c) shows a slight decrease in sea level contribution with temperature: snowfall increases, because warmer air can hold more water vapour, and this dominates over the increase in mass loss due to melting [START_REF] Seroussi | ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century[END_REF][START_REF] Nowicki | Contrasting contributions to future sea level under CMIP5 and CMIP6 scenarios from the Greenland and Antarctic ice sheets[END_REF] . Finally, West Antarctica and the Peninsula (b, e) show little detectable temperature-dependence, due to an approximate cancellation across varying climate and ice sheet model predictions of snowfall accumulation and ice loss. Antarctic ice sheet results are discussed in detail later (see 'Antarctic focus').

The ice sheet contributions depend strongly on the Greenland glacier retreat and Antarctic sub-shelf basal melting parameters, which determine the sensitivity of the marine-terminating glaciers to ocean temperatures (and surface meltwater runoff for Greenland). Figure 2 shows these relationships; the Greenland parameter is defined such that more negative values correspond to further retreat inland.

Land ice contributions in 2100

We use probability distributions for global mean surface air temperature (Fig. 3a: FaIR simple climate model [START_REF] Mckenna | Stringent mitigation substantially reduces risk of unprecedented near-term warming rates[END_REF] ) and ice-ocean parameters (Figs. 3b and 3c show κ and g, which are derived from the original parameterisation studies; ice shelf collapse is assigned equal probability off/on) as inputs to the emulators. Time series projections for the land ice contribution under all scenarios are shown in Fig. 3d, and probability density functions at 2100 for the Greenland ice sheet, Arctic Canada North, the glacier total, and West and East Antarctica in Fig. 3e-i. The Antarctic ice sheet total under the NDCs is shown in (j). ('Riskaverse' projections in (d) and (j) are discussed later.) Density estimates are less smooth for the glacier and Antarctica totals than individual regions, because sums of regions are estimated by random sampling rather than deterministic integration; these samples are shown for Antarctica (j).

Our projections show that reducing greenhouse gas emissions from current and projected pledges under the Paris Agreement (NDCs) enough to limit warming to 1.5 °C (SSP1- 19) would nearly halve the land ice contribution to sea level at 2100 (Table 1: median decreases from 25 cm to 14 cm SLE). This halving is not evenly distributed across the three ice sources: Greenland ice sheet mass losses would reduce by 70%, glacier mass losses by about half, and Antarctica shows no significant difference between scenarios; this is not due to a lack of change in the Antarctica simulations themselves, but rather to the cancellation of mass gains and losses mentioned above.

Average rates of mass loss for each ice sheet and the glacier total are within 1-2 cm/century, of those of the 2013 IPCC Fifth Assessment Report [START_REF] Church | Climate Change 2013: The Physical Science Basis[END_REF] (see Methods: Comparison with IPCC assessments), and the updated assessment for RCP2.6 in the 2019 IPCC Special Report on the Oceans and Cryosphere in a Changing Climate (SROCC) [START_REF] Oppenheimer | IPCC Special Report on the Ocean and Cryosphere in a Changing Climate[END_REF] . However, SROCC revised the projection for Antarctica under RCP8.5 up to 11 cm/century, close to the upper end of our 66% interval for SSP5-85 (though our projections may omit a commitment contribution of up to about 2 cm/century; see Methods). Our results are therefore closer to the 2013 than 2019 IPCC assessment regarding the magnitude and unclear scenario-dependence for Antarctica.

Our 66% uncertainty intervals are narrower than the IPCC 66% (SROCC) and ³ 66% (AR5) uncertainty intervals, as would be expected from the latter being open-ended, except those for Greenland under SSP1-26: too few Greenland simulations were performed under low scenarios (RCP2.6, SSP1-26) to constrain the emulator variance (see Fig. 1a; Methods: 'Parameter interactions').

Emulation allows us to additionally assess the sensitivity of projections to uncertainties in their inputs as well as their robustness. If we use CMIP6 global climate models for the projections (Extended Data Figure 3), instead of FaIR, we find a slight increase in sea level contributions due to the larger proportion of models with high climate sensitivity to carbon dioxide [START_REF] Forster | Latest climate models confirm need for urgent mitigation[END_REF][START_REF] Meehl | Context for interpreting equilibrium climate sensitivity and transient climate response from the CMIP6 Earth system models[END_REF] : the 95 th percentile increases by 7 cm under SSP5-85. We estimate the potential impact of reducing uncertainty with future knowledge by using fixed values for temperature, or for the ice sheet retreat and basal melt parameters: the width of the 5-95% ranges reduce by up to 13% and 17% respectively (tests 2-4 in Methods: Sensitivity tests; Extended Data Table 3 and Extended Data Figure 4). In other words, the ice-ocean interface is a similar magnitude contributor to, or larger, uncertainty for these projections as global warming under a particular emissions scenario. When we assess the robustness of the projections to different selections and treatments of the ice sheet simulations, we find this makes very little difference (tests 2-4 in Methods: Robustness checks; Extended Data Table 4; Extended Data Figure 5).

Antarctic focus

No clear dependence on emissions scenario emerges for Antarctica. This is partly due to the opposite scenario-dependencies of West and East Antarctica regions (Fig. 3f andg). But the average response to emissions scenario for each region is also small. A key reason is the wide variety of changes in the atmosphere and ocean in the global climate models. Figure 4 shows ice sheet model simulations where both the high and low emissions scenario were run (two climate models for Greenland, three for Antarctica). For the Greenland ice sheet, all simulations predict increased mass loss under higher emissions (Fig. 4a: red shaded region).

For Antarctica, the picture is more complex, and mostly clustered according to the climate model. Many West Antarctica simulations show the same straightforward response as Greenland (Fig. 4b), particularly those that do not use the ISMIP6 basal melting parameterisation (see Methods). However, the West Antarctica simulations driven by CNRM-CM6-1 show the reverse, where mass gain through snowfall accumulation increases more under high emissions than mass loss (which is predominantly ocean-induced). (Note fewer simulations were driven by IPSL-CM5A-MR and CNRM-CM6-1 than by NorESM1-M, so their spread is necessarily smaller). East Antarctica and the Peninsula mostly also show this latter response, though some simulations show other combinations: more mass loss under low emissions than high, or mass loss under low emissions and mass gain under high.

It is challenging to evaluate which of these three climate models, or others used by ISMIP6, are most reliable for Antarctic climate change. Ocean conditions and accumulation show large spatio-temporal variability and are sparsely observed; models imperfectly represent important processes, and it is unclear whether the newer CMIP6 models have improved relative to CMIP5 [START_REF] Jourdain | A protocol for calculating basal melt rates in the ISMIP6 Antarctic ice sheet projections[END_REF][START_REF] Meredith | IPCC Special Report on the Ocean and Cryosphere in a Changing Climate[END_REF][START_REF] Naughten | Future Projections of Antarctic Ice Shelf Melting Based on CMIP5 Scenarios[END_REF][START_REF] Mottram | What is the Surface Mass Balance of Antarctica? An Intercomparison of Regional Climate Model Estimates[END_REF][START_REF] Roussel | Brief communication: Evaluating Antarctic precipitation in ERA5 and CMIP6 against CloudSat observations[END_REF] . Most of the climate models were from CMIP5, including NorESM1-M and IPSL-CM5A-MR, and were selected by their success at reproducing southern climatological observations (while also sampling a range of future climate responses) [START_REF] Barthel | CMIP5 model selection for ISMIP6 ice sheet model forcing: Greenland and Antarctica[END_REF] . NorESM-1M has a lower than average atmospheric warming, hence less snowfall, while IPSL-CM5A-MR is higher than average (particularly for East Antarctica) [START_REF] Barthel | CMIP5 model selection for ISMIP6 ice sheet model forcing: Greenland and Antarctica[END_REF] .

The newer CMIP6 models, including CNRM-CM6-1, were selected only by their availability.

Changing the selection or treatment of Antarctica simulations -e.g. using subsets of climate models, or rejecting simulations with net mass gain early in the projections -do not result in any substantial scenario-dependence (see tests 7-10 in Methods: Robustness checks; Extended Data Table 4; Extended Data Figure 5).

Uncertainty about the scenario-dependence of Antarctic projections is not new. The IPCC Fifth Assessment Report (2013) stated 'the current state of knowledge does not permit an assessment' of the dependence of rapid dynamical change on scenario. Some studies that show strong scenario-dependence neglect the compensating accumulation part [START_REF] Levermann | Projecting Antarctica's contribution to future sea level rise from basal ice shelf melt using linear response functions of 16 ice sheet models (LARMIP-2)[END_REF][START_REF] Reese | The role of history and strength of the oceanic forcing in sea level projections from Antarctica with the Parallel Ice Sheet Model[END_REF] , use extreme 1 ice shelf collapse scenarios [START_REF] Edwards | Revisiting Antarctic ice loss due to marine ice cliff instability[END_REF] , or the basal melt parameterisation uncertainty is the same order as, or larger than, the scenario-dependence [START_REF] Bulthuis | Uncertainty quantification of the multi-centennial response of the Antarctic ice sheet to climate change[END_REF][START_REF] Golledge | The multi-millennial Antarctic commitment to future sea-level rise[END_REF][START_REF] Golledge | Global environmental consequences of twenty-first-century icesheet melt[END_REF] . To be clear, we do not assert that Antarctica's future does not depend on future greenhouse emissions or global warming: only that the relationship between global and Antarctic climate change, and the ice sheet's response, are complex, only partially understood, and involve compensating factors of increasing mass loss and gain which result in a balance we are not yet confident about.

We test the sensitivity of the Antarctica projections to the basal melting parameter. The main projections combine two distributions [START_REF] Jourdain | A protocol for calculating basal melt rates in the ISMIP6 Antarctic ice sheet projections[END_REF] for g derived from observations of mean Antarctic basal melt rates or the ten highest melt rates for Pine Island Glacier (see Methods). Using the mean distribution decreases the median to ~0 cm SLE and the 95 th percentile to ~8 cm SLE for all scenarios; using the high distribution has less effect, increasing the median to 6 cm SLE and the 95 th percentile to ~16 cm SLE (Extended Data Table 3 and Extended Data Figure 4: tests 5 and 6). We also try and reproduce the higher projections of ref. [26] using a similar approach to sampling basal melt (see Methods), and find we only obtain similar projections when using extreme values of our parameter range (Extended Data Table 3 and Extended Data Figure 4: tests 7 and 8). This suggests ref. [26] could be interpreted as more pessimistic projections: they use values of basal melt sensitivity to ocean temperature consistent with those estimated for the Amundsen Sea region [START_REF] Reese | The role of history and strength of the oceanic forcing in sea level projections from Antarctica with the Parallel Ice Sheet Model[END_REF] , which is currently undergoing most change.

However, other factors can lead to similarly high projections. In particular, the sensitivity of an individual ice sheet model to the basal melt parameter can have a large effect. This differs widely across ice sheet models, and also depends on the climate model (Extended Data Figure 6). Emulator projections based on a single model with high or low sensitivity are shown in Extended Data Figure 5 (tests 4 and 5; Extended Data Table 4). These also do not show strong scenario-dependence -just a 2-3 cm decrease under high emissions for the low sensitivity model, because the snowfall effect is more apparent -but instead predict a high or low sea level contribution, respectively, regardless of scenario (95 th percentiles: 29-30 cm and 7-9 cm, respectively). The high sensitivity of the first model (SICOPOLIS) is probably due to the way that sub-shelf melting is applied: over entire grid cells along the grounding line, rather than just the parts detected as floating [START_REF] Levermann | Projecting Antarctica's contribution to future sea level rise from basal ice shelf melt using linear response functions of 16 ice sheet models (LARMIP-2)[END_REF] . We also show results from the four most sensitive models, which are similarly high (Extended Data Table 4 and Extended Data Figure 5: test 6). We do not have sufficient observations to evaluate which ice sheet models have the most realistic response, nor sufficient understanding to confidently predict how basal melt sensitivity might change in future [START_REF] Jourdain | A protocol for calculating basal melt rates in the ISMIP6 Antarctic ice sheet projections[END_REF][START_REF] Naughten | Future Projections of Antarctic Ice Shelf Melting Based on CMIP5 Scenarios[END_REF] , and therefore use all models in the main projections (see also 'Risk-averse projections' below).

The ice shelf collapse scenario has little effect on our projections. Switching it on increases the Antarctic Peninsula and East Antarctic median contributions by 1 cm and 0-1 cm SLE from 2015-2100, with no change for West Antarctica (Extended Data Table 3 and Extended Data Figure 4: test 9-10). This is similar, within uncertainties, to the ice sheet simulations (Extended Data Figure 7). The effect is small because surface meltwater is not projected to be enough to cause collapses until the second half of the century, and even then only for small number of shelves, mostly around the Peninsula [START_REF] Seroussi | ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century[END_REF] . Some combinations of climate and ice sheet models do project larger sea level contributions -in particular, 5 cm for East Antarctica from the SICOPOLIS ice sheet model driven by HadGEM2-ES. The HadGEM2-ES climate model projects extreme ocean warming in the Ross Sea [START_REF] Barthel | CMIP5 model selection for ISMIP6 ice sheet model forcing: Greenland and Antarctica[END_REF] , while SICOPOLIS has one of the largest responses among the ice sheet models (as described above). If these two were found to be the most realistic models, then the ISMIP6 ensemble and emulator may underestimate the effect of ice shelf collapse by a few centimetres. Further results are in the Methods ('Parameter interactions').

Risk-averse projections

Given the wide range and cancellations of responses across models and parameters, we present alternative 'pessimistic but physically plausible' Antarctica projections for risk-averse stakeholders, by combining a set of assumptions that lead to high sea level contributions.

These are: the four ice sheet models most sensitive to basal melting; the four climate models that lead to highest Antarctic sea level contributions, and the one used to drive most of the ice shelf collapse simulations; the high basal melt (Pine Island Glacier) distribution; and with ice shelf collapse 'on' (i.e. combining robustness tests 6 and 7 and sensitivity tests 6 and 10). This storyline would come about if the high basal melt sensitivities currently observed at Pine Island Glacier soon become widespread around the continent; the ice sheet responds to these with extensive retreat and rapid ice flow; and atmospheric warming is sufficient to disintegrate ice shelves, but does not substantially increase snowfall. The risk-averse projections are more than five times the main estimates: median 21 cm (95 th percentile range 7 to 43 cm) under the NDCs (Fig. 3j), and essentially the same under SSP5-85 (Table 1; regions shown in Extended Data Figure 4: test 11), with the 95 th percentiles emerging above the main projections after 2040 (Fig. 3d). This is very similar to projections [START_REF] Edwards | Revisiting Antarctic ice loss due to marine ice cliff instability[END_REF] under an extreme scenario of widespread ice shelf collapses for RCP8.5 (median 21 cm; 95 th percentile range 9 to 39 cm). The median is higher than ref. [26] for RCP8.5, though the 95 th percentile is smaller. No models that include a representation of rapid ice cliff collapse through the proposed 'Marine Ice Cliff Instability' [START_REF] Deconto | Contribution of Antarctica to past and future sea-level rise[END_REF] mechanism participated in ISMIP6. This hypothesis is the process with the largest estimated systematic impact on projections: it could increase projections by tens of centimetres, if both the mechanism and projections of extreme ice shelf collapse are found to be robust [START_REF] Edwards | Revisiting Antarctic ice loss due to marine ice cliff instability[END_REF][START_REF] Clerc | Marine Ice Cliff Instability Mitigated by Slow Removal of Ice Shelves[END_REF] .

Our risk-averse Antarctica projections increase the total land ice sea level contribution to 42 cm (95 th percentile 25 to 67 cm) SLE under current policies and pledges (NDCs), and to 30 cm (95 th percentile 12 to 56 cm) SLE even under SSP1-19. This means that plausible modelling choices for Antarctica could change the median land ice contribution by more (17 cm SLE) than the difference between these emissions scenarios (12 cm SLE). This ambiguity limits confidence in assessing the effectiveness of mitigation on the response of global land ice to climate change. When combined, the effects of uncertain emissions and Antarctic response lead to a threefold spread in median projections of the land ice contribution to sea level rise, ranging from 13 to 42 cm SLE over 2015-2100, implying that flexible adaptation under substantial uncertainty will be essential until either can be further constrained. Not all modelling uncertainties could be systematically assessed here. Aside from the ice cliff instability hypothesis, these include ice sheet basal hydrology and sliding; glacier model parameters, ice-water interactions, and meltwater routing; model initialisation; and the use of coarse resolution global climate models (and a single high-resolution regional model for the Greenland ice sheet). The probabilities we present are therefore specific to our ensembles, and adding new climate and ice sheet models, or exploration of new parameters, could shift or broaden their distributions [START_REF] Williamson | How are emergent constraints quantifying uncertainty and what do they leave behind?[END_REF] . However, our projections demonstrate the importance of systematic design to assess as many uncertainties as feasible, and represent the current stateof-the art in estimating the land ice contribution to global mean sea level rise. All authors contributed to the manuscript. West Antarctica East Antarctica

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • SSP1-26 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • SSP2-45 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • NDCs • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • SSP3-70 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • SSP5-85 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Land ice
• MIROC5 • CNRM-CM6-1 • Not using parameterisation a • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -20 -10 0
• NorESM1-M • CNRM-CM6-1 • IPSL-CM5A-MR • Not using parameterisation b • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0
• NorESM1-M • CNRM-CM6-1 • IPSL-CM5A-MR • Not using parameterisation c • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
Sea level contribution at 2100 for RCP2. The glacier regions are listed in Extended Data Table 2 and all simulations are described in ref [8]. Greenland ice sheet projections have the peripheral glaciers (region 5) masked out, so there is no double-counting. The Antarctic periphery glaciers (region 19) are located only on the surrounding islands, not on the mainland ice sheet; ice sheet models include some of the larger islands, so there is some overlap in area, but the effect of this is estimated to be small (see 'Comparison with IPCC assessments' for an estimate of this and other limitations).

NorESM1-M • CNRM-CM6-1 • IPSL-CM5A-MR • Not using parameterisation d • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Methods

Simulations

All projections are calculated as annual global mean sea level contributions since 2015, converting mass (for the glaciers) or mass above flotation (for the ice sheets) to sea level contribution using 362.5 Gt per mm SLE.

Global climate model simulations

We use projections of annual global mean surface air temperature change since 2015 from the CMIP5 and CMIP6 global climate models used to drive the ice sheet and glacier models to build the emulator. If multiple realisations (different initial conditions) for a model were available, we use the mean of these. Data from 1850-2100 were downloaded from the JASMIN/CEDA archive and ESGF on the 7 th November 2019 and and 4 th December 2019; the CMIP6 snapshot was updated 28 th -29 th July 2020.

Emulation

An emulator is a fast statistical approximation of a computationally expensive simulator. This can be used to predict the simulator response at untried input values -to explore the uncertain input space far more thoroughly -for sensitivity analysis, to adjust the chosen inputs, and to estimate probability distributions. We construct statistical models of the simulated ice sheet and glacier sea level contribution as a function of the global mean surface air temperature of the driving climate models -and also different representations of the ice sheet-ocean interface -to make predictions under new emissions scenarios that incorporate these uncertainties, as well as those arising from the different structures of the climate and ice sheet models (and the emulators themselves).

Typically emulation is performed for one model at a time [START_REF] Edwards | Revisiting Antarctic ice loss due to marine ice cliff instability[END_REF] , but here we emulate each multimodel ensemble all at once. This is made possible by the systematic design of the ISMIP6 and GlacierMIP projects, which explore uncertainties in global climate change and three iceocean parameters simultaneously, and by our approach of applying emulation to multiple models rather than (as is usual) one. The three ice-ocean parameters control: (1) how much Greenland marine-terminating glaciers retreat (κ) with increasing local ocean temperatures and meltwater runoff; (2) how much Antarctic ice-shelf basal melting (g) increases with increasing local ocean temperature; and (3) an on/off scenario of Antarctic ice shelf collapse (C), which can increase glacier flow into the ocean when atmospheric temperatures rise [START_REF] Trusel | Divergent trajectories of Antarctic surface melt under two twentyfirst-century climate scenarios[END_REF] .

We predict the 23 land ice regions separatelythe Greenland ice sheet, the West and East Antarctic ice sheets and Antarctic Peninsula, and 19 glacier regionsso the spatial distribution of meltwater can be used in regional sea level projections.

We choose and evaluate emulator structures using the year 2100 (Extended Data Table 2;

Extended Data Figures 1 and2). Global mean surface air temperature projections are taken from the FaIR simple climate model [START_REF] Mckenna | Stringent mitigation substantially reduces risk of unprecedented near-term warming rates[END_REF] , because it can explore uncertainties more thoroughly than the relatively small CMIP6 ensemble of (computationally expensive) general circulation models. We use the same global mean temperature value across all land ice sources for each individual estimate: in other words, we include any co-dependence arising from global temperature. Full details are described in the following sections.

Global mean surface air temperature

Previous sea level emulation studies [START_REF] Church | Climate Change 2013: The Physical Science Basis[END_REF][START_REF] Levermann | Projecting Antarctica's contribution to future sea level rise from basal ice shelf melt using linear response functions of 16 ice sheet models (LARMIP-2)[END_REF][START_REF] Nauels | Synthesizing long-term sea level rise projections -the MAGICC sea level model v2.0[END_REF][START_REF] Palmer | Exploring the drivers of global and local sea-level change over the 21st century and beyond[END_REF] have typically used global mean temperature as the main input, rather than regional climate variables. We follow this approach for several reasons: to include correlation of land ice regions induced by global climate change (i.e. no need to assume/estimate their correlations, or to treat them as independent), and to have a larger sample of climate change projections. Using regional climate variables would improve the signal to noise for the emulator, but would restrict us to using computationally expensive general circulation models from CMIP5/6, for which there only a few tens of models. The simple climate model FaIR can be used to explore uncertainties in each scenario thoroughly, using the latest assessments of equilibrium climate sensitivity.

Global mean temperature is the only regressor for the glacier regions. For the ice sheets, there are additional terms derived from the ISMIP6 parameterisations of ice-ocean interactions.

Ice sheet model parameters

The Greenland glacier retreat parameter κ (Fig. 3a; units km (m 3 s -1 ) -0.4 °C-1 ) is a scaling coefficient relating marine-terminating glacier retreat to ocean temperatures and meltwater runoff [START_REF] Slater | Estimating Greenland tidewater glacier retreat driven by submarine melting[END_REF][START_REF] Slater | Twenty-first century ocean forcing of the Greenland ice sheet for modelling of sea level contribution[END_REF] , where larger negative values indicate greater retreat of the glacier terminus in response to warming. This is a continuous variable, but most simulations use one of three values: the default, which is the median of the distribution in the parameterisation [START_REF] Slater | Twenty-first century ocean forcing of the Greenland ice sheet for modelling of sea level contribution[END_REF] , κ50 = -0.17, and the quartiles κ25 = -0.37 and κ75 = -0.06. One model uses 5 th and 95 th percentile values, κ5 = -0.9705 and κ95 = 0.0079. For ice sheet models that did not use this parameterisation (N = 29 simulations) [START_REF] Goelzer | The future sea-level contribution of the Greenland ice sheet: a multimodel ensemble study of ISMIP6[END_REF] , we assign the mean value from the other simulations robust parameter estimation [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] . There are 23 emulators for the 2100 projections (Greenland ice sheet, three Antarctic ice sheet regions, and 19 glacier regions) and 1955 emulators for the full land ice time series (23 regions for each year from 2016 to 2100). An alternative to predicting each year separately would be to model the temporal correlation explicitly, but we prefer to use the simpler method, with fewer judgments, and allow temporal correlation to emerge.

Nugget

We use a 'nugget 'term to incorporate simulations from each multi-model ensemble. The nugget is usually zero for deterministic models -the emulator predicts each simulation in the ensemble exactly, i.e. the regression curve goes through all points -or a very small value, to improve numerical stability or other properties [START_REF] Andrianakis | The effect of the nugget on Gaussian process emulators of computer models[END_REF][START_REF] Gramacy | Cases for the nugget in modeling computer experiments[END_REF] . Here we allow the emulator to estimate the nugget, and treat each multi-model ensemble as a set of outputs from a single stochastic simulator or set of noisy observations. This approach has previously been used for emulating stochastic simulators [START_REF] Van Beers | Kriging for interpolation in random simulation[END_REF] and for emulating climate models accounting for internal variability, other inert inputs (uncertainties not explicitly modelled in the emulator), and approximations of the model outputs [START_REF] Salter | A comparison of statistical emulation methodologies for multi-wave calibration of environmental models[END_REF][START_REF] Williamson | Evolving Bayesian Emulators for Structured Chaotic Time Series, with Application to Large Climate Models[END_REF][START_REF] Williamson | Identifying and removing structural biases in climate models with history matching[END_REF][START_REF] Araya-Melo | Global sensitivity analysis of the Indian monsoon during the Pleistocene[END_REF][START_REF] Bounceur | Global sensitivity analysis of the climate-vegetation system to astronomical forcing: an emulator-based approach[END_REF][START_REF] Lord | Emulation of long-term changes in global climate: application to the late Pliocene and future[END_REF] . Our method is similar to the use of 'emergent constraints' for climate models [START_REF] Clerc | Marine Ice Cliff Instability Mitigated by Slow Removal of Ice Shelves[END_REF][START_REF] Bowman | A hierarchical statistical framework for emergent constraints: Application to snow-albedo feedback[END_REF] , seeking relationships between past and future simulations across multimodel ensembles to constrain them with observations, but here the predictors are inputs to the models rather than their outputs for the past.

This approach does not require the simulations to be normally distributed but does assume they are independent, which has been a long-standing difficulty of interpreting multi-model climate ensembles. But with ice sheet models, although model names may be the same across groups, each one has a very different set up, including physics approximations, parameterisations, tuning, grid resolution, and -in particular -initialisation methods, which have been shown to produce very different results even for simulations produced by the same group [START_REF] Goelzer | Design and results of the ice sheet model initialisation experiments initMIP-Greenland: an ISMIP6 intercomparison[END_REF][START_REF] Seroussi | initMIP-Antarctica: an ice sheet model initialization experiment of ISMIP6[END_REF][START_REF] Goelzer | The future sea-level contribution of the Greenland ice sheet: a multimodel ensemble study of ISMIP6[END_REF][START_REF] Seroussi | ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century[END_REF][START_REF] Nowicki | Insights into spatial sensitivities of ice mass response to environmental change from the SeaRISE ice sheet modeling project I: Antarctica[END_REF][START_REF] Nowicki | Insights into spatial sensitivities of ice mass response to environmental change from the SeaRISE ice sheet modeling project II: Greenland[END_REF][START_REF] Saito | SeaRISE experiments revisited: potential sources of spread in multi-model projections of the Greenland ice sheet[END_REF] . For glacier models, their structures are also vastly different, ranging from simple scaling parameterisations to dynamic physical models [START_REF] Marzeion | Partitioning the Uncertainty of Ensemble Projections of Global Glacier Mass Change[END_REF] . We test two approaches to account for any model dependence: a dummy variable (see below) and random effects ('Antarctic cross-check model').

Statistical model

Let y denote the simulated global mean sea level contribution for given region and year (in cm SLE), and x the simulator inputs (see below). Following ref. [22], we write the simulator as a function y = f(x), for which the Gaussian Process emulator is described by a mean function:

E[f(x)] = h(x) T β,
where h(x) is a vector of regression functions and β the corresponding regression coefficients, and a covariance function, with variance σ 2 and correlation function c(x, x′),

Cov[ f(x), f(x′) ] = σ 2 (c(x, x′) + νI),
where ν is the nugget term and I the identity matrix. So the prior for f(x) is:

p(f(x) | β, σ 2 , δ, ν)) ∼ N(h(x) T β, σ 2 (c(x,x′) + νI)),
where x are whichever model inputs are used for a given region, δ are the correlation lengths of the covariance function, and σ 2 ν is the variability not explained by the inputs. Parameters (β, σ 2 , δ, ν) are estimated from the simulation data.

The inputs x used in the regression functions are global mean temperature change, T, and, for the ice sheets, the ice-ocean parameter values (κ for Greenland; γ, C for Antarctica), plus a dummy variable denoting whether Greenland models used the retreat parameterisation. These are discussed in the next section. All inputs are rescaled to have zero mean and unit variance.

Mean functions

The Gaussian Process mean function describes the large-scale response of the simulator to its inputs, usually specified as a linear trend with the remainder described by a zero-mean Gaussian process.

For the glaciers, the linear regressor is simply global mean temperature in the same year (T).

For the ice sheets, the additional ice sheet model parameters are κ for Greenland, and g and C

for Antarctica. We also try two types of dummy variable. The first is for the ice sheet and glacier model names, so these can be treated distinctly in the emulator, but this leads to clear overfitting (i.e. the model is too flexible in Figs. 1 and2).

Covariance functions

The covariance function describes the smoothness of the Gaussian Process. As in any statistical modelling, there is a trade-off between improving accuracy and over-fitting. We assess this using the usual leave-one-out procedure [START_REF] Rougier | Analyzing the Climate Sensitivity of the HadSM3 Climate Model Using Ensembles from Different but Related Experiments[END_REF][START_REF] Bastos | Diagnostics for Gaussian Process Emulators[END_REF] . We fit the emulator to all ensemble members but one, then predict the sea level contribution from this simulation; we repeat this for every combination, noting the emulator error (residual) and uncertainty for each prediction. We perform this for each of the 23 regional emulators for the year 2100 with five covariance functions of varying smoothness -Matérn(5/2), which is the default in RobustGaSP, Matérn (3/2), and three members of the power exponential family with high, medium and low exponent values (a = 1.9, i.e. close to a squared exponential, the default value; a = 1.0, exponential, and a = 0.1, for which the covariance function has a small effect so the emulator approaches linear regression).

For 18 of the 19 glacier regions, we use the covariance function with the smallest standardised Euclidean distance between the emulator predictions and simulations (standardised because, unlike simpler metrics such as root mean square error or mean absolute error, it does not penalise larger errors if the emulator uncertainty intervals are sufficiently large), as in ref [24]. For the Southern Andes (region 17), all covariance functions give identical distances, so we use the default for RobustGaSP. For the ice sheets, we use the covariance function that gives close to linear regression (power exponential, a = 0.1), rather than the one with the minimum Euclidean distance, for various reasons. For Greenland, West Antarctica, and the Antarctic Peninsula, the minimum distance covariance functions (power exponential a = 1.0 for Greenland; Matérn(3/2) for the Antarctic regions) result in overfitting for temperature (i.e. too much flexibility in Fig. 1). For East Antarctica, the minimum distance covariance functions (Matérn(5/2)) result in an incorrect sign prediction under the ice shelf collapse switch. Using the alternative covariance function solves all of these issues and does not increase the standardised Euclidean distance by much: 4% for the Peninsula, and 0.4-1% for the other three regions. The resulting covariance functions are given in Extended Data Table 2.

Evaluating the emulators

After selecting the covariance functions for each regional emulator at 2100, we evaluate the emulators further by plotting the emulator predictions against the simulations from the leaveone-out procedure, and the standardised residuals (the difference between the emulator prediction and the simulator, divided by the emulator standard deviation), and calculating the percentage of simulations falling within ± 2 s.d. (Extended Data Table 2 and Extended Data Figures 1 and2). We would not expect exactly 95% of the simulations to fall within 2 s.d., in part because the predictions are not independent, but very low or high values would suggest emulator over-or under-confidence. The region with the lowest percentage of predictions within the uncertainty intervals is North Asia (region 10) with 89%, indicating slightly too small emulator uncertainty estimates, and the highest is 98% (Scandinavia: region 8), indicating the reverse.

Mean absolute errors for each emulator are given in Extended Data Table 2 and Extended Data Figures 1 and2: for the ice sheet regions they are 0.28 cm (Peninsula), 1.4 cm (Greenland) and 1.5 cm (East Antarctica) and 2.0 cm (West Antarctica), and for the For the mean function (also linear), we use the logarithm of g as a regressor, so it is always positive. Consequently we use the geometric mean as the missing value, rather than the arithmetic mean. We use a dummy variable to denote these models, as for Greenland in the

Sea level projections

We use probability distributions for global temperature and the ice sheet model parameters as inputs to each emulator to make the projections.

Global mean temperature projections

We use projections of global annual mean surface air temperature change since 2015 from the FaIR (Finite amplitude Impulse Response) simple climate model for the main projections.

We take the 500-member ensemble from reference [30]: SSP1-19, SSP1-26, SSP3-70, SSP5-85 and a scenario estimated for the 2019 Nationally Determined Contributions. We also use projections for SSP-245 generated with the same ensemble.

Ice sheet model parameter distributions

For Greenland, we sample from a kernel density estimate of the original k distribution (N = 191) with the same bandwith used in deriving the parameterisation [START_REF] Slater | Estimating Greenland tidewater glacier retreat driven by submarine melting[END_REF][START_REF] Slater | Twenty-first century ocean forcing of the Greenland ice sheet for modelling of sea level contribution[END_REF] (0.0703652) (Fig. 1b).

The dummy variable is always set to represent the standard ISMIP6 parameterisation.

For Antarctica, we combine the Mean Antarctic and Pine Island Glacier g distributions (N = 10,000 each), and sample from a kernel density estimate using three times the automatic bandwidth (Silverman's 'rule of thumb' [START_REF] Silverman | Density Estimation[END_REF] ) to merge and smooth them into a near-unimodal distribution that we truncate at zero (Fig. 1c). For the collapse switch C, we sample randomly from 0 or 1 with equal probability (8% of the ISMIP6 simulations have ice shelf collapse).

The ice shelf collapse scenario does not include the possibility of surface meltwater draining efficiently from some ice shelves under certain conditions, thereby avoiding collapse, so we feel this is a reasonable judgement.

Sampling

For the 2100 projections, we sample from the FaIR ensemble (N=500) with replacement (N = 5000 for main and risk-averse projections; N = 1000 for robustness and sensitivity tests). For the full time series, we use the 500 FaIR projections directly without resampling. We make one set of emulator predictions (23 regions) for each temperature value in a given year, randomly sampling the relevant ice-ocean parameters (k, g0, C) once for each FaIR ensemble member.

for RCP8.5. The difference between scenarios for Antarctica in AR5 arises only from additional accumulation, because the dynamic contributions are assumed to be the same.

Glacier projections could be overestimated because meltwater routing to the ocean is not accounted for (not all volume lost from the glaciers reaches the oceans), or underestimated because only one glacier model includes ice-water interactions (i.e. frontal ablation of marine-and lake-terminating glaciers). For the latter, we compare mean projections for the GloGEM model to the emulator for RCP8.5/SSP5-85 and RCP4.5/SSP2-45 for key regions, and find they are larger by less than 1 cm for Alaska and Russian Arctic (regions 1 and 9), by less than 0.5 cm for Svalbard (7) and Arctic Canada South (4), and smaller than the emulator for Arctic Canada North (3). All are within the emulator 95 th percentile estimates. We may slightly underestimate uncertainty in the global glacier total due to correlated errors across models 8 by emulating the regions independently, though there are compensating advantages (more accurate emulation; spatial pattern of meltwater); a similar argument applies to Antarctica.

Sensitivity tests

We perform a number of checks to test the sensitivity of the ice sheet projections to changes in the chosen inputs, predominantly the input distributions, but also the dataset in the final test (see Extended Data Table 3 and refs [25, 26,30, 34, 39]). All results are shown for the SSP5-85 scenario in Extended Data Figure 4 under the index given (where 1 is the main projection); numerical values in the text refer to changes in the median and [5,95 th ] percentile estimates for the ice sheet under this scenario unless otherwise stated.

Robustness checks

We perform a number of checks to test robustness of the ice sheet projections to changes in the simulation dataset (see Extended Data Table 4 and refs [14, 16, 24, 66]). Results are shown for the NDCs scenario in Extended Data Figure 5 under the test index given (where 1

is the main projection); numerical values in the text refer to changes in the median and [5,95 th ] percentile estimates under this scenario unless otherwise stated. The full datasets are 256 simulations for Greenland and 344 simulations for Antarctica.

Figure 1 .

 1 Figure 1. Ice sheet and glacier mass loss generally increases linearly with global mean temperature. Projected mass changes from 2015-2100 in sea level equivalent (SLE) as a function of global mean surface air temperature change over the same period for (a) Greenland ice sheet, (b, c) West and East Antarctic ice sheets, (d) Greenland peripheral glaciers, (e, f) the Antarctic Peninsula and Antarctic peripheral glaciers, (g-j) four glacier regions with large maximum sea level contributions (Alaska, Arctic Canada North and South, Russian Arctic), (k, l) two regions with nonlinear temperature-dependence and total or near-total disappearance projected at high temperatures (Central Europe and Caucasus); and (m-o) three regions comprising High Mountain Asia. Central solid lines show the emulator mean, and shaded regions the mean ± 2 s.d.. For the ice sheets (a-c, e), darker shaded regions use parameter values fixed at their default values (Greenland glacier retreat: median; Antarctic sub-shelf basal melting: median of Mean Antarctic distribution; Antarctic ice shelf collapse off), and lighter shaded regions use alternative values (Greenland: 75 th percentile; Antarctica: median of Pine Island Glacier distribution). See Methods for details. Points show ice sheet and glacier simulations under RCP2.6/SSP1-26 (blue), RCP4.5 (yellow), RCP6.0 (orange) and RCP8.5/SSP5-85 (red). Solid circles for the ice sheets use the default ice-ocean parameter value and open circles use the alternative value (other simulations are not shown). Glacier

Figure 2 .

 2 Figure 2. Ice sheet mass loss strongly depends on ice-ocean parameters. Projections of sea level contribution from 2015-2100 as a function of (a) Greenland glacier retreat parameter (κ), and basal melt parameter (g) for (b) West Antarctica, (c) East Antarctica, (d) Peninsula. Solid line shows emulator mean estimate using fixed global temperature (projected by the global climate model most used for simulations, under RCP8.5), and shaded regions show the mean ± 2 s.d. Symbols show ice sheet models forced by this climate model for which simulations for at least three (Greenland) or four (Antarctic) melt parameter values were available: circles use the ISMIP6 parameterisation for the iceocean interface; crosses use other representations, and are assigned ensemble mean values of the parameter; triangles show the Greenland ice sheet model for which two additional values of κ were run.

Figure 3 .

 3 Figure 3. Projected land ice contribution to 21 st century sea level rise and for selected regions at 2100. (a) Probability distributions for global mean surface air temperature change from 2015-2100 from the FaIR simple climate model under the five Shared Socioeconomic Pathways (SSPs) and current Nationally Determined Contributions (NDCs) (N = 5000 each). (b) Greenland ice sheet retreat parameter (κ) distribution (N = 10,000): vertical lines show the five values used for simulations: median (solid), 25 th and 75 th percentiles (dashed), and 5 th and 95 th percentiles (dotted). (c) Antarctic basal melt parameter (g) distribution (N = 8200): vertical lines show the six values used for simulations: median (solid), 5 th and 95 th percentiles (dashed) of the Mean Antarctic (black) and Pine Island Glacier (grey) distributions (see Methods). (d) Projected land ice contribution to sea level (cm SLE) from 2015-2100 under the five SSPs and NDCs. Solid lines and shaded regions: median and 5-95 th percentiles (N = 11,500 per year per scenario): 5 year smoothing applied, with original data shown as dots (interannual variation arises from annual sampling of emulator uncertainties). Pale solid lines: 95 th percentiles of risk-averse projections. Box and whiskers show [5, 25, 50, 75, 95] th percentiles at 2100 (N = 115,000 per scenario) for main projections (left) and risk-averse projections for Antarctica (right).(e-j). Probability density functions for 2100 estimated for: (e) Greenland ice sheet, (f) Arctic Canada North, (g) total for glaciers, (h, i) West and East Antarctica for all scenarios, and (j) total for Antarctic ice sheet under main and risk-averse projections for the NDCs. Glacier and Antarctic totals are less smooth because they are estimated from a sum of Monte Carlo samples from each region, rather than deterministic integration (see Methods); these samples are shown for SSP1-19 and NDCs (N = 5000). Ice sheet projections do not include pre-2015 response, which is estimated to add less than 1 cm to the Greenland contribution and up to ~2 cm to the Antarctic (see Methods).

Figure 4 .

 4 Figure 4. Climate and ice sheet projections show a wide range of responses to greenhouse gas emissions scenario. Sea level contribution at 2100 under high greenhouse gas emissions scenarios (RCP8.5 or SSP5-85) versus low scenarios (RCP2.6 or SSP1-26), categorised by climate model forcing (NorESM1-M and IPSL-CM5A-MR use RCPs; CNRM-CM6-1 use SSPs), without ice shelf collapse. a, Greenland. b, West Antarctica. c, East Antarctica. d, Antarctic Peninsula. Filled circles show ice sheet models that use the ISMIP6 parameterisations of the ice-ocean interface, while open circles show models that used their own. Simulations in the red shaded regions have more mass loss under high emissions (RCP8.5/SSP5-85) than low (RCP1-26/SSP1-26); those in the green shaded regions have more mass gain under high emissions scenarios than low. Two regions with other possible combinations are also labelled.

1 Frequency

 1 parameter, κ (km (m 3 s -1 ) -0.4 °C-1 ) Sea level contribution at 2100 (cm SLE) parameter, γ (x10 m a -1 )Sea level contribution at 2100 (cm SLE) parameter, γ (x10 3 m a -1 ) Sea level contribution at 2100 (cm SLE) parameter, γ (x10 m a -1 )Sea level contribution at 2100 (cm SLE)Greenland glacier retreat parameter, κ (km (m 3 s -1 ) -0.4 °C-

  05 4e-05 6e-05 8e-05 1e-04 Sea level contribution at 2100 (cm SLE) at 2100 for RCP2.6/SSP1-26 (cm SLE) Sea level contribution at 2100 for RCP8.5/SSP5-85 (cm SLE) MASS LOSS DOMINATES More mass loss under high emissions Greenland

  at 2100 for RCP2.6/SSP1-26 (cm SLE) Sea level contribution at 2100 for RCP8.5/SSP5-85 (cm SLE) MASS LOSS DOMINATES More mass loss under high emissions ACCUMULATION DOMINATES More mass gain under high emissions More mass loss under low emissions Mass loss under low emissions, mass gain under high emissions

5 - 5 0 5 Sea

 555 level contribution at 2100 for RCP2.6/SSP1-26 (cm SLE) Sea level contribution at 2100 for RCP8.5/SSP5-85 (cm SLE) MASS LOSS DOMINATES More mass loss under high emissions ACCUMULATION DOMINATES More mass gain under high emissions More mass loss under low emissions Mass loss under low emissions, mass gain under high emissions

  GP emulator. The full global mean temperature change trajectories are used instead of only the total change at 2100. To increase the signal-to-noise ratio, the annual means are reduced to decadal means (2015-2029, 2030-2039, . . . , 2090-2100). There are thirteen distinct forcings, each one the product of a global climate model and a scenario, so we represent the forcing variables as twelve bisquare basis functions. These start as thirteen bisquare basis functions, each one centred at one of the thirteen forcings, but one is dropped because otherwise the model matrix becomes rank deficient when a constant is added. The one dropped is the one with the smallest mean Euclidean distance to the other twelve. We use bisquare kernels, where the standard deviation of each kernel is set to one tenth of the maximum Euclidean distance between all pairs of forcings, to cover the forcing space with non-zero values for the forcing regressors. We use the same distributions for temperature, basal melt and collapse as the main projections, and set the dummy variable to represent standard parameterisation models. This emulator predicts 50[5, 95] th percentiles for the West Antarctic sea level contribution at 2100 of 2 [-4, 8] cm SLE for SSP1-26 and 3 [-4, 10] cm SLE for SSP5-85, which are very similar to the GP emulator predictions of 2 [-5, 10] cm SLE and 3[-4, 11] cm SLE. We test the effect of changing the kernel standard deviation to one twelfth or one fourteenth of the maximum Euclidean distance; the largest change is a 2 cm decrease in the 95 th percentile under SSP5-85. For East Antarctica, the emulator with random effects predicts 2 [-3, 6] cm SLE for both scenarios; the GP emulator predicts a small scenario-dependence, 2[-4, 7] cm SLE for the low emissions scenario and 0[-5, 6] cm SLE for the high. For the Antarctic Peninsula, the random effects predictions are 0 [-1, 2] cm SLE for both scenarios, and the GP are the same. These similarities give us confidence that model dependence is not substantially affecting our projections -i.e. that differences in model structure, resolution, calibration and initialisation dominate over the similarities -although it would be worth investigating this in more detail.

Sea level contribution from 2015-2100 (cm SLE)

  Energy (DOE), Office of Science, Advanced Scientific Computing Research and Biological and Environmental Research Programs. Nicholas R. Golledge, Daniel P. Lowry and Brian Anderson were supported by NZ Ministry for Business, Innovation and . contributed to the ISMIP6 experimental design, organisation and analysis as members of its steering committee, and R.S and W.H.L. led the ISMIP6 atmosphere focus group. C.M.L., A.B. and C.A. selected the CMIP5 models for ISMIP6, X.F. and P.A. ran the surface mass balance model for the Greenland and R.Cu. prepared the Antarctic surface mass balance, and L.D.T. and M.v.d.B. provided the ice shelf collapse forcing. For Antarctica: T.K. and A.H.

		Main projections		Risk-averse projections
		50 [5, 95]%	[17, 83]%	50 [5, 95]%	[17, 83]%
		percentiles	percentiles	percentiles	percentiles
	Employment contracts RTUV1705 (`NZSeaRise') and ANTA1801 (`Antarctic Science
	Global glaciers Platform'). Jonathan Gregory and Robin S. Smith were supported by the National Centre for
	SSP119 Atmospheric Science, funded by the UK National Environment Research Council. Reinhard 7 [4, 10] [5, 9] contributed the AWI/PISM simulations; M.H., T.Z. and S.P. contributed the DOE/MALI
	SSP126 Calov was funded by the PalMod project of the Bundesministerium für Bildung und 8 [5, 12] [6, 10] simulations; R.G. and R.Ca. contributed the ILTS_PIK/SICOPOLIS simulations; H.G. and
	SSP245 Forschung (BMBF) with the grants FKZ 01LP1502C and 01LP1504D. Daniel Martin and 11 [7, 15] [9, 13] R.v. d. W. contributed the IMAU/IMAUICE simulations; N.-J.S. and H.S. contributed the
	NDCs	13 [9, 18]	[11, 16]
	SSP370 Courtney Shafer were supported by the Director, Office of Science, Offices of Advanced 14 [10, 19] [12, 17] JPL/ISSM simulations; C.D. and A.Q. contributed the LSCE/GRISLI simulations; G.L. and
	SSP585 Scientific Computing Research (ASCR) and Biological and Environmental Research (BER), 16 [12, 21] [14, 19] W.L. contributed the NCAR/CISM simulations; R.R., T.A. and R.W. contributed the
	Greenland ice sheet
	SSP1-19	2 [-6, 11]
	SSP1-26	3 [-4, 12]
	SSP2-45	5 [-2, 14]
	NDCs	7 [0, 16]
	SSP3-70	8 [0, 17]
	SSP5-85	10 [2, 20]
	Helene Seroussi was supported by grants from NASA Cryospheric Science, Sea Level contributed the GloGEMflow simulations; S.S. contributed the JULES simulations; P.K. and
	Change Team, and Modeling, Analysis, and Predictions Programs. Torsten Albrecht and W.I. contributed the KRA2017 simulations; B.M. and J.M. contributed the MAR2012
	Ricarda Winkelmann are supported by the Deutsche Forschungsgemeinschaft (DFG) in the simulations; F.M. and N.C. contributed the OGGM simulations; D.R. and R.H. contributed
	framework of the priority program "Antarctic Research with comparative investigations in the PyGEM simulations; A.B. and V.R. contributed the RAD2014 simulations; R.v.d.W.
	contributed the WAL2001 simulations; and A. Bl. and J.-H. M. assisted with data handling.

of of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, as a part of the ProSPect SciDAC Partnership. BISICLES simulations used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. Chen Zhao and Ben Galton-Fenzi were supported under the Australian Research Council's Special Research Initiative for Antarctic Gateway Partnership (Project ID SR140300001) and received grant funding from the Australian Government for the Australian Antarctic Program Partnership (Project ID ASCI000002). Work was performed by Eric Larour, Nicole Schlegel, and Helene Seroussi at the California Institute of Technology's Jet Propulsion Laboratory under a contract with the National Aeronautics and Space Administration's Cryosphere, Sea Level Change Team, and Modeling, Analysis and Prediction (MAP) Programs. They acknowledge computational resources and support from the NASA Advanced Supercomputing Division. The CMIP5 and CMIP6 projection data were processed by Christine McKenna with funding from the European Union's CONSTRAIN project as part of the Horizon 2020 Research and Innovation Programme under grant agreement number 820829. Alice Barthel was supported by the DOE Office of Science HiLAT-RASM project and Early Career Research program. Arctic ice areas" by grants WI4556/2-1 and WI4556/4-1, and within the framework of the PalMod project (FKZ: 01LP1925D) supported by the German Federal Ministry of Education and Research (BMBF) as a Research for Sustainability initiative (FONA). Ronja Reese is A.SPIK/PISM simulations;. T.P., M.M. and H.S. contributed the UCIJPL/ISSM simulations; F.P. and S.S. contributed the ULB/fETISh simulations; C.Z., R.G., B.G-F. and T.Z. contributed the UTAS/Elmer/Ice simulations; J.V.B. and P.H. contributed the VUB/AISMPALEO simulations; N.R.G. and D.L. contributed the VUW/PISM simulations; and D.F.M. and C.S. contributed the CPOM/BISICLES simulations. For Greenland: M.R. and A.H. contributed the AWI/ISSM simulations; V.L. and A.J.P. contributed the BGC/BISICLES simulations; I.N., D.F. and S.N. contributed the GSFC/ISSM simulations; R.G., R.Ca. and C.C. contributed the ILTS_PIK/SICOPOLIS simulations; H.G., R.v.d.W. and M.v.d.B. contributed the IMAU/IMAUICE simulations; N.-J.S. and H.S. contributed the JPL/ISSM simulations; J.C. and N.-J.S. contributed the JPL/ISSMPALEO simulations; A.Q. and C.D. contributed the LSCE/GRISLI simulations; L.T. contributed the MUN/GSM simulations; W.H.L. and G.R.L. contributed the NCAR/CISM simulations; A.A contributed the UAF/PISM simulations; Y.C., H.S. and M.M. contributed the UCIJPL/ISSM simulations; S.L.c. and P.H. contributed the VUB/GISM simulations; and D.P.L. and N.R.G. contributed the VUW/PISM simulations. For global glaciers: B.A. contributed the AND2012 simulations; K.F. and A.S. contributed the GLIMB simulations; M.H. contributed the GloGEM simulations; H.Z.

Table 1 . Projected land ice contributions to sea level rise in 2100 under different greenhouse gas scenarios and Antarctic modelling assumptions

 1

. Projected changes to global glaciers, Greenland and Antarctic ice sheets and land ice total from 2015-2100 in sea level equivalent (cm SLE) for five Shared Socioeconomic Pathways (SSPs) and predicted emissions under the 2019 Nationally Determined Contributions (NDCs). Ice sheet projections do not include pre-2015 response, which is estimated to add less than 1 cm to the Greenland contribution and ~2 cm to the Antarctic (see Methods). The glaciers include the Greenland and Antarctic ice sheet peripheral glaciers; the overlap of Antarctic periphery glaciers with the ice sheet contribution is estimated to be less than 1 cm SLE.

  The second represents whether an ice sheet model uses the ISMIP6 retreat or basal melt parameterisation, to absorb any misalignment between the imputed value and the effective value. Bayesian Information Criterion (BIC) from a stepwise model selection (testing up to first-order interactions) suggests this dummy variable is informative for Greenland, so we retain it (o, for open parameterisation), but not for the Antarctic regions. The stepwise model selection suggests we could reasonably include terms for the interaction between temperature and retreat for Greenland, temperature and basal melt for West Antarctica, and temperature and collapse for East Antarctica, but we choose not to, to avoid the risk of overfitting. The selection also shows that collapse strongly dominates the Antarctic Peninsula response, and is may not be needed for West Antarctica, but we retain all terms (i.e. Ti, g0, C) because we otherwise find the covariance matrix is poorly conditioned. The resulting mean functions are hGrIS(x)i ~ (Ti,

k, o) for Greenland, hAIS(x)i ~ (Ti, g0, C) for the Antarctic regions, and hGlaciers(x)i ~ (Ti) for the glaciers, where h ~ (a,b) means h is a linear function of a and b, and i is the index for the year.

Acknowledgements

We thank Jonathan Rougier for generously providing advice and support throughout, and writing the original random effects model. We also thank Baylor Fox-Kemper, Helene Hewitt, Robert Kopp, Sybren Drijfhout and Jeremy Rohmer for useful discussions, suggestions and support. We thank Daniel Williamson, Nicholas Barrand and two anonymous referees for their thorough and constructive comments, which greatly improved the manuscript. We thank the Climate and Cryosphere (CliC) effort, which provided support for ISMIP6 and GlacierMIP through sponsoring of workshops, hosting the websites and ISMIP6 wiki, and promotion. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP5 and CMIP6. We thank the climate modeling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the CMIP data and providing access, the University at Buffalo for ISMIP6 data distribution and upload, and the multiple funding agencies who support CMIP5 and CMIP6 and ESGF. We thank the ISMIP6 steering committee, the ISMIP6 model selection group and the ISMIP6 dataset preparation group for their continuous engagement in defining ISMIP6. This is ISMIP6 contribution No. 13. This publication was supported by PROTECT, which has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 869304. This is PROTECT contribution number XX. Individual author acknowledgements follow. Tamsin Edwards was supported by PROTECT and the UK Natural Environment Research Council grant NE/T007443/1. Fiona Turner was supported by PROTECT. James O'Neill was supported by the UK Natural Environment Research Council London Doctoral Training Partnership. Rupert Gladstone's contribution was supported by Academy of Finland grants 286587 and 322430. William Lipscomb and Gunter Leguy were supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement No. 1852977. Computing and data storage resources for CISM simulations, including the Cheyenne supercomputer (doi:10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. Support for Xylar Asay-Davis, Matthew J. Hoffman, Stephen Price, and Tong Zhang was provided through the Scientific Discovery through Advanced Computing (SciDAC) program funded by the US Department supported by the Deutsche Forschungsgemeinschaft (DFG) by grant WI4556/3-1 and through the TiPACCs project that receives funding from the European Union's Horizon 2020 Research and Innovation program under grant agreement no. 820575. Ralf Greve and Christopher Chambers were supported by Japan Society for the Promotion of Science (JSPS) KAKENHI grant Nos. JP16H02224 and JP17H06323. Ralf Greve was supported by JSPS KAKENHI grant No. JP17H06104, by a Leadership Research Grant of Hokkaido University's Institute of Low Temperature Science (ILTS), and by the Arctic Challenge for Sustainability (ArCS) project of the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) (program grant number JPMXD1300000000). Frank Pattyn and Sainan Sun were supported by the MIMO project within the STEREO III programme of the Belgian Science Policy Office, contract SR/00/336 and the Fonds de la Recherche Scientifique (FNRS) and the Fonds Wetenschappelijk Onderzoek-Vlaanderen (FWO) under the EOS Project number O0100718F. Andrew Shepherd was supported by the UK Natural Environment Research Council in partnership with the Centre for Polar Observation and Modelling and the British Antarctic Survey and by the European Space Agency Climate Change Initiative. Denis Felikson was supported by an appointment to the NASA Postdoctoral Program at the NASA Goddard Space Flight Center, administered by Universities Space Research Association under contract with NASA.

Code availability

37 R code and input data are available at https://github.com/tamsinedwards/emulandice. Each simulation in the sea level projections file has a label in the 'publication' column for the reference (Goelzer2020, Seroussi2020, Nowicki2020 or Marzeion2020), or 'New' if previously unpublished.

Data availability

All global climate, simple climate, ice sheet and glacier model data used as inputs to this study are provided with the code as described above. Main and risk-averse projections from the analysis are provided in the Supplementary Information as annual quantiles for each of the 23 regions, and the Antarctic, glacier and land ice sums.

Author contributions

T.L.E. conceived the idea, carried out all statistical analysis except the random effects model, produced the figures, and wrote the manuscript. S.N. led ISMIP6, including experimental design, organisation and analysis, and provided scientific interpretation. B.M and R.H. co-led GlacierMIP and contributed simulations (below), and provided data and interpretation. H. G. and H.S. led the processing and analysis in ISMIP6 for the Greenland and Antarctic ice sheets, respectively, contributed simulations (below), and provided scientific interpretation and advice. N.J. and. D.S. co-derived with T.L.E. the ice sheet continuous parameter distributions for the emulator, and also derived the corresponding ocean forcing parameterisation studies with X.A.-D. and T.H. for Antarctica and F.S., D.F. and M.M. for Greenland. F.T. performed the random effects model cross-check for Antarctica. C.S. provided the FaIR projections and C.M. provided the CMIP5 and CMIP6 projection data for the emulator. E.S. led the ISMIP6 data processing. A.A.O., J.M.G., E.L., W.H.L., A.J.P., to minimise the impact on the emulator (κ = -0.2073). One of these models (BISICLES) also ran 'high' and 'low' retreat experiments by doubling and halving the ocean thermal forcing, to which we assign the κ25 and κ75 values.

The Antarctic sub-shelf basal melt parameter g (Fig. 3b; units m a -1 ) is the 'ocean heat exchange velocity' scaling coefficient relating sub-shelf basal melting to ocean temperatures [START_REF] Favier | Assessment of sub-shelf melting parameterisations using the oceanice-sheet coupled model NEMO(v3.6)-Elmer/Ice(v8.3)[END_REF][START_REF] Jourdain | A protocol for calculating basal melt rates in the ISMIP6 Antarctic ice sheet projections[END_REF] . Two alternative distributions for g were derived in the parameterisation [START_REF] Jourdain | A protocol for calculating basal melt rates in the ISMIP6 Antarctic ice sheet projections[END_REF] : the first from mean Antarctic melt rates, and the second from the 10 highest observations of melt rate at the grounding line of Pine Island Glacier, where melt rates are currently highest.

The values of g estimated from Pine Island Glacier are an order of magnitude larger, and the two distributions do not overlap. This is a continuous variable, but most simulations use one of three values: the default, which is the median of the Mean Antarctic distribution, MeanAnt50 = 14477, and the 5 th and 95 th percentiles, MeanAnt5 = 9619 and MeanAnt95 = 21005. Further simulations used the same percentiles from the Pine Island Glacier distribution: PIG50 = 159188, PIG5 = 86984 and PIG95 = 471264. Some models [START_REF] Seroussi | ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century[END_REF] used an alternative variant of the parameterisation in which only local ocean temperatures were used, rather than a combination of local and regional, which uses a different tuning for g. However, the values used are also the 50 [5, 95] th percentiles of those distributions, so we consider them equivalent. For ice sheet models that did not use this parameterisation (N = 62 simulations), we again assign the ensemble mean value (g = 59317).

The Antarctic ice shelf collapse parameter C is a switch that indicates whether a scenario of ice shelf collapse was used, which can lead to glacier speed-up. A timeline of collapses was derived according to the presence of surface meltwater on ice shelves above a threshold (725 mm a -1 ) for 10 years, estimated from surface air temperature projections [START_REF] Trusel | Divergent trajectories of Antarctic surface melt under two twentyfirst-century climate scenarios[END_REF] in the global climate model driving the ice sheet model (mostly CCSM4). This method does not predict whether meltwater may be efficiently drained from the surface for a given ice shelf [START_REF] Bell | Antarctic ice shelf potentially stabilized by export of meltwater in surface river[END_REF] , thus avoiding collapse. We use values of 1 or 0 indicating whether the scenario is implemented or not.

Gaussian Process emulation

Gaussian Process emulation [START_REF] O'hagan | Bayesian analysis of computer code outputs: A tutorial[END_REF] is non-parametric, treating the simulator as an unknown mathematical function of its inputs. We use the R package RobustGaSP [START_REF] Gu | RobustGaSP: Robust Gaussian Stochastic Process Emulation in R[END_REF] for its numerically individual glacier regions they range from 0.0020 cm to 0.87 cm (Antarctic periphery: region 19). Mean absolute standardised errors are all less than 0.006.

The emulator underestimates the three to four highest West and East Antarctic contributions by around 10-15 cm (Extended Data Figure 1b and1c). The five highest of these are from the SICOPOLIS model, which has a much greater sensitivity to basal melting than other models (see main text, Robustness checks and Extended Data Figure 6), and use the highest value of this parameter (g = PIG95). These simulations are therefore extreme: 1% of the 344 simulations, and the 97.5 th percentile value of the basal melt parameter. There are processbased reasons to expect that SICOPOLIS is an upper bound or overestimate (see main text).

When the emulator is calibrated with this model alone, it does not underestimate its highest contributions (not shown). The resulting projections under the NDC scenario are shown in Robustness checks (test 4); the difference with the main projections may be interpreted as the maximum possible impact of this emulator underestimate, if SICOPOLIS were the sole realistic ice sheet model. These are lower than the 'risk-averse' projections, which are made with a subset of high sensitivity ice sheet models and other pessimistic assumptions (see main text).

We therefore consider the emulators to be adequate for the predictions of large-scale sea level contribution presented here.

Antarctic cross-check model

We perform a cross-check for the Antarctic ice sheet regions at 2100 using a linear mixed model, with the ice sheet model name included as a random effect to deal with any systematic uncertainty arising from dependence of ensemble members. This attributes some of the uncertainty in the response to the ice sheet model used, and this uncertainty can then be removed from the predicted PDF. We thus model the ensemble members as 'similar but not identical', using a mean function of temperature and ice sheet parameters, plus a structured error term which includes a systematic component according to the ice sheet model and a noise component to capture other sources of variability such as initialisation.

We integrate over the uncertain inputs (temperature in a given year, and ice-ocean parameters) to obtain the final probability density functions (PDFs). Each regional emulator predicts a Student-t distribution for a given set of these input values, defined by a mean and standard deviation; we approximate this with a normal distribution, as in refs [55, 57], which is accurate enough for this application. We use different integration methods for the 23 individual regional PDFs compared with the regional sums (Antarctica, global glaciers, and land ice total). For the individual regional estimates, we use deterministic numerical integration (the midpoint rule: we sum the Gaussian distributions for each emulator prediction, then normalise). For regional sums we must use Monte Carlo sampling, because the three ice sources (Greenland, Antarctica and glaciers) have different parameters, and we also desire traceability of predictions to input values within a given ice source. We sample once from the Gaussian distribution for each emulator prediction, then sum the regional samples for a given temperature to estimate the PDF, smoothing with kernel density estimation for figures (again using Silverman's 'rule of thumb' [START_REF] Silverman | Density Estimation[END_REF] for the bandwidth). Sampling is a more noisy method of integration than deterministic methods, so the PDFs for regional sums are less smooth than those for individual regions.

Glacier maximum cap

We apply a cap to the glacier projections using estimates of their maximum sea level contribution [START_REF] Farinotti | A consensus estimate for the ice thickness distribution of all glaciers on Earth[END_REF] . Glacier model projections often exceed this cap in some regions, if near or total loss is projected under high emissions, either because they report changes in total mass, not mass above flotation, or because of errors in initial mass [START_REF] Marzeion | Partitioning the Uncertainty of Ensemble Projections of Global Glacier Mass Change[END_REF] , or both. We restrict values to the maximum in the emulator mean predictions and then the PDFs (the latter exceeding the cap due to emulator uncertainty).

Time series smoothing

Interannual variability arises in the time series due to sampling the emulator uncertainty for each annual regional prediction. We apply a five year running mean in Fig. 3d to visualise the expected smoothness of sea level contributions; projections provided in the Supplementary Information are unsmoothed.

Comparison with IPCC assessments

The ice sheet projections are made relative to control simulations with a constant recent climate. This control includes both the model drift and, depending on the initialisation method, any background contribution arising from forcing before 2015. This background contribution should be added to the ice sheet projections, but is difficult to quantify. Five year mean rates of sea level contribution since 1992/3 range from 0.1-0.8 mm/yr for the Greenland ice sheet 65 and 0.1-0.6 mm/yr for Antarctica 66 , but they would decrease in the absence of forcing after 2014. Modelling work to quantify the background contribution from Greenland 67 suggests a contribution of 0.6 ± 0.2 cm SLE by 2100. Estimates made for this study range from 0.3-0.8 cm under a range of retreat parameter values, κ75 -κ25 (IMAU/IMAUICE1: 0.3-0.4 cm; CISM variant similar to NCAR/CISM: 0.4-0.8 cm). For Antarctica, the dynamic commitment has been estimated to be 2 cm SLE at 2100 for the Amunden Sea Embayment region of West Antarctica, where most mass loss is currently occurring [START_REF] Alevropoulos-Borrill | Ocean-forced evolution of the Amundsen Sea catchment, West Antarctica, by[END_REF] . Part of these trends may still be due to residual model drift. The committed contribution could therefore add up to ~1 cm/century to our Greenland projections and ~2 cm/century to the Antarctic.

The Antarctic ice sheet models include some of the larger islands that are also included in region 19, potentially leading to double-counting. However, median projections for region 19 range from 1-2 cm under different emissions scenarios, and the ice sheet models are much lower resolution (i.e. the glaciers are likely less responsive), so the effect is expected to be of order 0.5-1 cm SLE or less.

We average our projections over the 86 years and compare them with the average IPCC AR5 [START_REF] Church | Climate Change 2013: The Physical Science Basis[END_REF] and SROCC 1 projections over 95 years (the midpoints of 1986-2005 to 2081-2100) as rates of cm SLE per century. For the glaciers, we project 8 cm/century SLE for SSP1-26 and 16 cm/century for SSP5-85 excluding the Antarctic peripheral glaciers (region 19: 1 cm and 2 cm, respectively), compared with 10 cm for RCP2.6 and 17 cm for RCP8.5 in AR5. For the Greenland ice sheet, we project 4 cm/century SLE for SSP1-26 and 11 cm for SSP5-85, compared with 6 cm for RCP2.6 and 13 cm for RCP8.5 in AR5. For Antarctica, we project 5 cm/century SLE for both scenarios; the AR5 projections are 5 cm/century SLE for RCP2.6 and 4 cm for RCP8.5, while those for SROCC are 4 cm/century SLE for RCP2.6 and 11 cm

Parameter interactions

Retreat and basal melt vs temperature

Ice sheet projection uncertainties are constant across scenarios. However, tests with three ice sheet models show that the range of projections from high to low values of the retreat parameter (κ95 -κ5) and basal melt parameter (PIG95 -MeanAnt50) is consistently smaller under RCP2.6 than RCP8.5, so the emulator uncertainty should be smaller at lower temperatures. The ratios of ranges, RCP2.6/RCP8.5, for each group/model + GCM are:

Greenland

• IMAU/IMAUICE + MIROC5 = 1.4097/8.3069 = 0.17 The emulator does not have sufficient data from lower emissions scenarios to reduce the variance, particularly for Greenland. If other ice sheet models respond the same way as the above, then adding more simulations may reduce the uncertainty for low SSPs.

Ice shelf collapse vs basal melt

The contribution due to ice shelf collapse does not increase with higher values of the basal melt parameter in the models JPL1/ISSM and CPOM/BISICLES (0.1 cm difference for the Peninsula in BISICLES; all other regional differences for both models £ 0.02 cm).

Author information

The authors declare no competing financial or non-financial interests. Correspondence and requests for materials should be addressed to T.L.E. (tamsin.edwards@kcl.ac.uk). Reprints and permissions information is available at www.nature.com/reprints.

Extended Data

Extended Data Table 1. The additional 22 Greenland and 37 Antarctic ice sheet model experiments not previously described elsewhere. Retreat parameter values κ5 and κ95 are the 5 th and 95 th percentile values of the retreat (κ) distribution; basal melt parameter values MeanAnt [5, 50, 95] and PIG [5, 50, 95] 3). Box and whiskers show [5, 25, 50, 75, 95] 4). Box and whiskers show [5, 25, 50, 75, 95]