the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Uncertainty quantification of the multi-centennial response of the Antarctic ice sheet to climate change

### Maarten Arnst

### Sainan Sun

### Frank Pattyn

Ice loss from the Antarctic ice sheet (AIS) is expected to become the major contributor to sea level in the next centuries. Projections of the AIS response to climate change based on numerical ice-sheet models remain challenging due to the complexity of physical processes involved in ice-sheet dynamics, including instability mechanisms that can destabilise marine basins with retrograde slopes. Moreover, uncertainties in ice-sheet models limit the ability to provide accurate sea-level rise projections. Here, we apply probabilistic methods to a hybrid ice-sheet model to investigate the influence of several sources of uncertainty, namely sources of uncertainty in atmospheric forcing, basal sliding, grounding-line flux parameterisation, calving, sub-shelf melting, ice-shelf rheology and bedrock relaxation, on the continental response of the Antarctic ice sheet to climate change over the next millennium. We provide probabilistic projections of sea-level rise and grounding-line retreat, and we carry out stochastic sensitivity analysis to determine the most influential sources of uncertainty. We find that all investigated sources of uncertainty, except bedrock relaxation time, contribute to the uncertainty in the projections. We show that the sensitivity of the projections to uncertainties increases and the contribution of the uncertainty in sub-shelf melting to the uncertainty in the projections becomes more and more dominant as atmospheric and oceanic temperatures rise, with a contribution to the uncertainty in sea-level rise projections that goes from 5 % to 25 % in RCP 2.6 to more than 90 % in RCP 8.5. We show that the significance of the AIS contribution to sea level is controlled by the marine ice-sheet instability (MISI) in marine basins, with the biggest contribution stemming from the more vulnerable West Antarctic ice sheet. We find that, irrespective of parametric uncertainty, the strongly mitigated RCP 2.6 scenario prevents the collapse of the West Antarctic ice sheet, that in both the RCP 4.5 and RCP 6.0 scenarios the occurrence of MISI in marine basins is more sensitive to parametric uncertainty, and that, almost irrespective of parametric uncertainty, RCP 8.5 triggers the collapse of the West Antarctic ice sheet.

- Article
(6203 KB) -
Supplement
(4790 KB) - BibTeX
- EndNote

The Antarctic ice sheet (AIS) is the largest reservoir of freshwater on Earth (∼60 m sea-level equivalent, Fretwell et al., 2013; Vaughan et al., 2013) and has the potential to become one of the largest contributors to sea level in the next centuries. Yet, studies such as the IPCC Fifth Assessment Report (AR5) (Stocker et al., 2013) and Kopp et al. (2014) have identified the magnitude of the AIS response as the largest source of uncertainty in projecting sea-level rise, when compared with other contributors such as thermal expansion, glaciers and the Greenland ice sheet. Recent observations (Shepherd et al., 2018) have shown an acceleration in the rate of ice loss from the Antarctic ice sheet, especially in West Antarctica, where an irreversible retreat may be underway in the Amundsen Sea Embayment as a consequence of a marine ice-sheet instability (MISI) (Favier et al., 2014; Joughin et al., 2014; Rignot et al., 2014).

Assessing the future response of the Antarctic ice sheet requires numerical
ice-sheet models amenable to large-scale and long-term simulations and
quantification of the impact of modelling hypotheses and parametric
uncertainty. So far, there exist only a limited number of projections
(Golledge et al., 2015; Ritz et al., 2015; DeConto and Pollard, 2016; Schlegel et al., 2018) for the whole
Antarctic ice sheet on a (multi-)centennial timescale. These projections show
similar trends for the next centuries, but they differ in the magnitude of
the mass loss they predict, with differences and uncertainty ranges that can
reach several metres for eustatic sea level. Using a numerical ice-sheet
model supplemented with a statistical approach for the probability of MISI
onset, Ritz et al. (2015) have projected a contribution to sea level around
0.1 m by 2100, with a very low probability of exceeding 0.5 m in the A1B
scenario. Running their simulations with and without sub-grid melt
interpolation at the grounding line, Golledge et al. (2015) have projected that
sea-level rise could range from 0.01 to 0.38 m by 2100 and from 0.21 m to
more than 5 m by 2500 considering all RCP scenarios. Including the marine
ice-cliff instability mechanism (Pollard et al., 2015) in a numerical ice-sheet
model, DeConto and Pollard (2016) have projected that sea level could rise much more
significantly, with projections exceeding 1 m by 2100 and 15 m by 2500 in
the RCP 8.5 scenario. Synthesising recent results for the expected response
of the Antarctic ice sheet to climate change, Pattyn et al. (2018) have
highlighted that the projected contribution to sea level of the Antarctic ice
sheet is limited to well below 1 m by 2500 in the RCP 2.6 scenario, while a
key threshold for the stability of the Antarctic ice sheet is expected to lie
between 1.5 and 2 ^{∘}C mean annual air temperature above present and
the activation of larger ice systems, such as the Ross and Ronne–Filchner
drainage basins, could be triggered by global warming between 2 and
2.7 ^{∘}C. Emulating the numerical ice-sheet model by
DeConto and Pollard (2016) with a Gaussian process emulator, Edwards et al. (2019)
have established new probabilistic projections for the AIS contribution to
sea level by 2100, with the probability of exceeding 0.5 m that reaches
4 % in RCP 2.6 and 71 % in RCP 8.5 when considering the marine ice-cliff
instability mechanism. Coupling an ice-sheet model with a climate model,
Golledge et al. (2019) have shown that freshwater released by the Antarctic ice
sheet can trap warm waters below the sea surface, thus leading to higher
projections by 2100, with a contribution to sea level that could reach
0.05 m in RCP 4.5 and 0.14 m in RCP 8.5.

Despite recent progress in the numerical modelling of ice-sheet dynamics (Pattyn et al., 2017), differences in modelling hypotheses between ice-sheet models remain a major source of uncertainty for sea-level rise projections. Intercomparison projects such as the SeaRise project (Bindschadler et al., 2013) and the ongoing ISMIP6 project (Nowicki et al., 2016; Goelzer et al., 2018) aim at quantifying the impact of such so-called structural uncertainty in ice-sheet models by comparing projections from multiple numerical models. Such multimodel ensembles give some insight into the impact of the structural uncertainty, although challenges remain in combining the different results (Knutti et al., 2010).

Uncertainties in the ice-sheet initial state, climate forcing and parameters in numerical ice-sheet models are another major limitation for accurate projections. To date, the impact of such parametric uncertainty is assessed most often by using large ensemble analysis; that is, the model is run for different values of the parameters and the uncertainty in the projections is estimated from the spread in the model runs. For example, Golledge et al. (2015) ran their model with simplified sensitivity experiments to evaluate qualitatively the sensitivity of the Antarctic ice sheet to temperature, precipitation and sea-surface temperature individually, while DeConto and Pollard (2016) assessed the sensitivity of their model to parametric uncertainty by evaluating the model response at a few samples in the parameter space. By contrast, Ritz et al. (2015) adopted a probabilistic approach and estimated the uncertainty in sea-level rise projections by running an ice-sheet model with samples of parameters drawn randomly from probability density functions for the parameters and subsequently deducing probability density functions for sea-level rise projections.

The field of uncertainty quantification (UQ) develops theory and methods to describe quantitatively the origin, propagation and interplay of sources of uncertainty in the analysis and projection of the behaviour of complex systems in science and engineering; see, for instance, Ghanem et al. (2017) for a recent handbook and Arnst and Ponthot (2014) for a recent review paper. Most of this theory and these methods are based on probability theory, in the context of which uncertain parameters and projections are represented as random variables characterised by their probability density function. Theory and methods are under development to characterise sources of uncertainty by probability density functions inferred from observational data and expert assessment (characterisation of uncertainty), to deduce the impact of sources of uncertainty on projections (propagation of uncertainty) and to ascertain the impact of each source of uncertainty on the projection uncertainty and rank them in order of significance (stochastic sensitivity analysis). These developments have led to new theories and new methods that are of interest to be applied to uncertainty quantification of ice-sheet models, beyond the theory and methods that Golledge et al. (2015), Ritz et al. (2015), DeConto and Pollard (2016) and Schlegel et al. (2018) have already applied.

In this paper, we apply probabilistic methods to assess the impact of uncertainties on the continental AIS response over the next millennium. We use the fast Elementary Thermomechanical Ice Sheet (f.ETISh) model (Pattyn, 2017), a hybrid ice-sheet model that captures the essential characteristics of ice flow and allows large-scale and long-term projections at a reasonable computational cost. To reduce the computational cost of assessing the impact of uncertainties on the change in global mean sea level (GMSL), we draw from UQ methods based on the construction of an emulator (Le Maître and Knio, 2010; Ghanem et al., 2017), also known as a surrogate model, that is, a computational model that mimics the response of the original ice-sheet model at a reduced computational cost. To assess the significance of each individual source of uncertainty in inducing uncertainty in the change in GMSL, we draw from UQ methods for stochastic sensitivity analysis (Saltelli et al., 2008). To express the uncertainty in projections of the retreat of grounded ice, we draw from UQ methods for constructing confidence regions for excursion sets (Bolin and Lindgren, 2015; French and Hoeting, 2015), with a confidence region for grounded ice interpreted as a region of Antarctica that remains covered everywhere with grounded ice for a given level of probability.

On the one hand, our study adds to previous studies (Golledge et al., 2015; Ritz et al., 2015; DeConto and Pollard, 2016) that also provided projections for the multi-centennial response of the whole Antarctic ice sheet: whereas these previous studies provided projections by running the ice-sheet model with samples of parameters from the parameter space, our study differs by the adoption of additional methods from UQ and the analysis of a broader set of parameters that includes uncertainty in climate forcing, basal sliding, grounding-line flux parameterisation, calving, sub-shelf melting, ice-shelf rheology and bedrock relaxation. On the other hand, our study also adds to previous studies (Pollard et al., 2016; Schlegel et al., 2018; Edwards et al., 2019) that also applied methods from UQ for uncertainty propagation in ice-sheet models: whereas these previous studies analysed AIS paleoclimatic responses with Gaussian process (kriging) emulators (Pollard et al., 2016) and multi-decadal forecasts with Latin hypercube sampling (Schlegel et al., 2018), our study differs by its focus on the analysis of the multi-centennial response of the whole Antarctic ice sheet with polynomial emulators; and in addition to uncertainty propagation, we complement our uncertainty analysis with stochastic sensitivity analysis and confidence regions.

This paper is organised as follows. First, Sect. 2 describes the f.ETISh model, the uncertain processes and parameters and the UQ methods, including the use of an emulator, stochastic sensitivity analysis and the construction of confidence regions. Subsequently, Sect. 3 shows and interprets the results for the multi-centennial response of the Antarctic ice sheet as well as the advantages of the adopted UQ methods. Finally, Sect. 4 provides an overall discussion of the results.

## 2.1 Ice-sheet model and simulations

We perform simulations of the response of the Antarctic ice sheet (Fig. 1a) to environmental and parametric perturbations over the period 2000–3000 CE with the fast Elementary Thermomechanical Ice Sheet (f.ETISh) model (Pattyn, 2017) version 1.2. The f.ETISh model is a vertically integrated, thermomechanical, hybrid ice-sheet/ice-shelf model that incorporates essential characteristics of ice-sheet thermomechanics and ice-stream flow, such as the mass-balance feedback, bedrock deformation, sub-shelf melting and calving. The ice flow is represented as a combination of the shallow-ice (SIA) (Hutter, 1983; Greve and Blatter, 2009) and shallow-shelf (SSA) (Morland, 1987; MacAyeal, 1989; Weis et al., 1999) approximations for grounded ice (Bueler and Brown, 2009), while only the shallow-shelf approximation is applied for floating ice shelves. Bedrock deformation is represented as a combined time-lagged asthenospheric relaxation and elastic lithospheric response (Greve and Blatter, 2009; Pollard and DeConto, 2012a), in which the lithosphere relaxes towards isostatic equilibrium due to the viscous properties of the underlying asthenosphere. Calving at the ice front is parameterised based on the large-scale stress field, represented by the horizontal divergence of the ice-shelf velocity field (Pollard and DeConto, 2012a). Prescribed input data include the present-day ice-sheet geometry and bedrock topography from the Bedmap2 dataset (Fretwell et al., 2013), the basal sliding coefficient, and the geothermal heat flux by An et al. (2015).

We perform simulations at a spatial resolution of 20 km while accounting for grounding-line migration at coarse resolution with a flux condition derived from a boundary layer theory at steady state based on either a Weertman (or power) friction law (Schoof, 2007b) or a Coulomb friction law (Tsai et al., 2015) at the grounding line. This flux condition at the grounding line is imposed as an internal boundary condition following the implementation by Pollard and DeConto (2009, 2012a). This implementation has been shown to reproduce the migration of the grounding line and its steady-state behaviour (Schoof, 2007a) at coarse resolution for the SSA and hybrid SSA/SIA models, while these models can only reproduce the migration of the grounding line and its steady-state behaviour at sufficiently fine resolution in the absence of a flux condition (Docquier et al., 2011; Pattyn et al., 2012). Numerical simulations (Pollard and DeConto, 2012a; DeConto and Pollard, 2016) of the Antarctic ice sheet using a flux condition have also been able to reproduce MISI in large-scale ice-sheet simulations. In addition, Pattyn (2017) has shown that a flux condition makes the f.ETISh model rather independent of the resolution for a spatial resolution of the order of 20 km. While the implementation of a flux condition at the grounding line has been shown to reproduce qualitatively ice-sheet dynamics as determined with other ice-sheet models, it cannot reproduce quantitatively changes in ice-sheet mass and the contribution to sea level as determined from models with a higher level of complexity, such as the Blatter–Pattyn (Blatter, 1995; Pattyn, 2003) and full Stokes (Greve and Blatter, 2009) models that include vertical shearing at the grounding line, especially for short transients (decadal timescales) as these flux conditions have been derived at steady state (Drouet et al., 2013; Pattyn et al., 2013; Pattyn and Durand, 2013). A proper representation of grounding-line migration without the need for a flux condition would require a very fine grid resolution (possibly less than 200 m). In addition, the flux condition has been derived for unbuttressed ice shelves (Schoof, 2007b) and may fail to appropriately represent grounding-line migration for buttressed ice shelves (Reese et al., 2018b) as found mostly around Antarctica even when the flux condition is corrected for buttressing with a buttressing factor accounting for back stress at the grounding line. In particular, Pattyn et al. (2013) have shown that models that implement the flux condition cannot reproduce compressional stresses that may apply in the presence of buttressing. Moreover, using a spatial resolution of 20 km does not allow us to capture properly certain mechanisms that control grounding-line migration such as bedrock irregularities and ice-shelf pinning points even with sub-grid parameterisations of these mechanisms. Therefore, we may expect discrepancies between our results and results at higher spatial resolutions (<5 km), especially for important small ice streams such as Pine Island and Thwaites glaciers, which represent only a few grid points with a 20 km resolution. Despite the limitations of the flux condition and our coarse grid resolution, we have adopted these modelling assumptions as an efficient way to capture the essential mechanisms of grounding-line migration in large-scale, long-term and large-ensemble ice-sheet simulations while keeping the computational cost tractable. The computational cost of a single forward simulation with a time step of 0.05 years on the CÉCI clusters (F.R.S-FNRS & Walloon Region) on two threads is approximately 8 h. To investigate the impact of the spatial resolution on the results, we performed additional runs at a spatial resolution of 16 km (Fig. S1 in the Supplement). We found that the uncertainty in the projections due to the spatial resolution is (far) less important than the uncertainty due to the uncertainty in the parameters.

The main changes in f.ETISh version 1.2 as compared to version 1.0 (Pattyn, 2017) consist of the computation of sub-shelf melt rates with an ocean-model coupler based on the Potsdam Ice-shelf Cavity mOdel (PICO) ocean-model coupler (Reese et al., 2018a) rather than simpler parameterisations of sub-shelf melt rates (Beckmann and Goosse, 2003; Holland et al., 2008; Pollard and DeConto, 2012a; de Boer et al., 2015; Cornford et al., 2016), the inclusion of a laterally varying flexural rigidity for the lithosphere (Chen et al., 2018) to compute glacial isostatic adjustment more realistically with an elastic lithosphere/relaxed asthenosphere model, an improvement of the SSA numerical scheme based on the implementation by Rommelaere and Ritz (1996) and a description of atmospheric forcing based on a parameterisation of the changes in atmospheric temperature and precipitation rate (Huybrechts et al., 1998; Pollard and DeConto, 2012a), a parameterisation of surface melt with a positive degree-day model (Janssens and Huybrechts, 2000) and the inclusion of meltwater percolation and refreezing (Huybrechts and de Wolde, 1999).

We drive our simulations with both atmospheric and oceanic forcings.
Present-day mean surface air temperature and precipitation are obtained from
Van Wessem et al. (2014), based on the regional atmospheric climate model RACMO2.
Changes in atmospheric temperature and precipitation rate induced by a
forcing temperature change Δ*T* are applied in a parameterised way that
accounts for elevation changes (Huybrechts et al., 1998; Pattyn, 2017). We use a
positive degree-day (PDD) model to calculate surface melt
(Janssens and Huybrechts, 2000), assuming a PDD factor of 3 and
8 mm ^{∘}C^{−1} d^{−1} for snow and ice, respectively. The PDD
model also includes meltwater percolation and refreezing
(Huybrechts and de Wolde, 1999).

Basal melting underneath ice shelves is determined from the PICO ocean-model
coupler (Reese et al., 2018a), which evaluates sub-shelf melting from the ocean
temperature *T*_{oc} and salinity fields on the continental shelf
via an ocean box model that captures the basic overturning circulation within
ice-shelf cavities. We employ data by Schmidtko et al. (2014) for present-day
ocean temperature ${T}_{\mathrm{oc}}^{\mathrm{obs}}$ and salinity on the
continental shelf. For a change in background atmospheric temperature Δ*T*, we assume that the ocean temperature *T*_{oc} on the continental
shelf changes as

where *F*_{melt} is an ocean melt factor that represents the ratio
between oceanic and atmospheric temperature changes
(Maris et al., 2014; Golledge et al., 2015). Equation (1) with
*F*_{melt}=0.25 has been shown to reproduce trends in ocean
temperatures following an analysis of the Climate Model Intercomparison
Project phase 5 (CMIP5) data set (Taylor et al., 2012) for changes in
atmospheric and ocean temperatures (Golledge et al., 2015).
Figure 1b shows the initial sub-shelf melt rate as computed with
the PICO model. The PICO model is able to reproduce the general pattern of
sub-shelf melting with higher melt rates near the grounding line and lower
melt rates near the calving front. Sub-shelf melt rates are also higher in
the Amundsen Sea sector and lower underneath the large Ross and
Filchner–Ronne ice shelves.

The calibration of the basal sliding coefficient follows the data assimilation method of ice-sheet geometry by Pollard and DeConto (2012b). This approach is based on a fixed-point iteration scheme that adjusts the basal sliding coefficient iteratively so as to match the present-day ice-sheet configuration while assuming that the present-day configuration is in steady state. After applying this approach, we further adjust the basal sliding coefficient in ice streams with a sliding multiplier factor similar to Bindschadler et al. (2013) to reduce the initial drift. We carry out the calibration step independently for the different sliding and grounding-line flux conditions investigated in this paper. For the nominal values of the model parameters and without forcing anomaly, the initial drifts (reference runs) range from 0.1 to 0.2 m of sea-level rise by 3000. To make our analysis insensitive to model initialisation, we correct all our results for this initial bias by subtracting the reference runs from the results. Therefore, after corrections for present-day dynamical changes, our results reflect the model response to perturbations. These corrections have in general little impact on medium-term and long-term projections as the initial drift is small compared to the projections but leads to a more significant bias for short-term projections.

## 2.2 Sources of uncertainty

We aim at quantifying the response of the Antarctic ice sheet to climate change while accounting for uncertainty in atmospheric forcing, basal sliding, grounding-line flux parameterisation, calving, sub-shelf melting, ice-shelf rheology and bedrock relaxation. We account for uncertainty in these physical processes by introducing uncertainty in parameters in the f.ETISh model (see Table 1). Here, we choose the uncertainty ranges for the parameters sufficiently large so as to encompass extreme conditions. The effect of choosing narrower uncertainty ranges for the parameters is discussed in Sect. 3.6. We assume the parameters to take the same value everywhere over the whole Antarctic ice sheet. This approach is likely to be conservative as the parameters may vary at least regionally. We refer to Schlegel et al. (2018) for a discussion about uncertainty quantification based on a partitioning of the ice-sheet domain. In the remainder of this section, we list the uncertain parameters and briefly discuss their influence on the AIS response.

### 2.2.1 Atmospheric forcing

Alongside oceanic forcing, atmospheric forcing is generally considered to be
the primary driver of future changes in the AIS mass balance
(Lenaerts et al., 2016; Pattyn et al., 2018). Simulating atmospheric forcing over several
centuries with regional climate models is computationally prohibitive. In
addition, the accuracy of the climate models is limited by uncertainties such
as the possible trajectories of anthropogenic greenhouse gas emissions.
Therefore, we adopt here the four schematic extended representative
concentration pathway (RCP) scenarios RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5
introduced by Golledge et al. (2015) for temperature changes Δ*T* in the
atmosphere (Fig. 2) as a means of representing different
atmospheric forcings relevant for policymakers.

The scenario plays a significant role in the amplitude and speed of the AIS retreat. Recent studies (Golledge et al., 2015; DeConto and Pollard, 2016) have shown no substantial retreat of the grounding-line position in the strongly mitigated RCP 2.6 scenario. The other scenarios lead to a reduction in the extent of the major ice shelves (Ross, Filchner–Ronne and Amery ice shelves) within 100–300 years, leading to accelerated grounding-line recession due to reduced buttressing. DeConto and Pollard (2016) have also highlighted that the hydrofracturing and ice-cliff failure mechanisms (not included in f.ETISh version 1.2) driven by increased surface melt and sub-shelf melting could potentially lead to an accelerated collapse of the West Antarctic ice sheet and a deeper grounding-line retreat in the East Antarctic subglacial marine basins.

### 2.2.2 Basal sliding

Basal sliding controls the motion of fast-flowing ice streams, which drain about 90 % of the total Antarctic ice flux (Bennett, 2003). Several studies have shown the importance of basal sliding for the behaviour of ice streams and stressed the need for understanding of physical processes at play at the ice–bedrock interface (Joughin et al., 2009, 2010; Ritz et al., 2015; Brondex et al., 2017, 2019). In particular, Ritz et al. (2015) have shown that, under a power-law rheology for basal sliding, the contribution to future sea level is an increasing function of the sliding exponent.

We introduce basal sliding as a Weertman sliding law, that is,

where *τ*_{b} is the basal shear stress,
*v*_{b} the basal
velocity, *A*_{b} the basal sliding coefficient and *m* a sliding
exponent. The sliding exponent is often related to Glen's flow law exponent
*n* as *m*=*n* for sliding over hard bedrock (Weertman, 1957). The value
*m*=3, related to the usual value *n*=3 for Glen's flow law exponent, has
been applied in a number of studies
(Schoof, 2007a; Pattyn et al., 2012, 2013; Brondex et al., 2017; Gladstone et al., 2017), but the
value *m*=1 has also been commonly used in ice-flow models
(Larour et al., 2012; Schäfer et al., 2012; Gladstone et al., 2014; Yu et al., 2018).

In addition to the usual exponents *m*=1 and *m*=3, we consider for the
sliding exponent a nominal value of 2 as an intermediate sliding condition
between the two usual exponents used in ice-flow models. Hereafter, we refer
to the exponents $m=\mathrm{1},\mathrm{2}$ and 3 as the viscous (or linear), the weakly
non-linear and the strongly non-linear sliding law, respectively. We consider
the discrete values $m=\mathrm{1},\mathrm{2}$ and 3 as representative of the most common
values in large-scale ice-sheet modelling and discuss in
Sect. 3.7 the impact of a more plastic sliding law
(*m*=5) to represent a quasi-plastic deformation of the till in ice streams
(Gillet-Chaulet et al., 2016).

### 2.2.3 Grounding-line flux parameterisation

The f.ETISh model employs a parameterisation of the grounding-line flux based on a boundary layer theory at steady state by either Schoof (2007b) (SGL) or Tsai et al. (2015) (TGL). For the TGL parameterisation, Pattyn (2017) has shown an increased AIS contribution to sea level and a more significant retreat of the grounding line. We consider SGL parameterisation as the reference parameterisation for most of the simulations and discuss the impact of TGL parameterisation in Sect. 3.8.

We applied the TGL parameterisation under the weakly non-linear sliding law. There is no consensus on the compatibility between the TGL parameterisation and the Weertman sliding law, as the TGL parameterisation was derived from the Coulomb friction law near the grounding line (Tsai et al., 2015). Yet, the Coulomb friction law is applicable in a narrow transition region in the vicinity of the grounding line not resolved on our coarse mesh, which lends validity to the combination of the TGL parameterisation and Weertman's sliding law.

### 2.2.4 Calving

Ice loss due to ice calving at the edges of ice shelves is responsible for almost half of the present-day ice mass loss of the Antarctic ice sheet (Depoorter et al., 2013; Rignot et al., 2013). Iceberg calving can have strong feedback effects as it affects ice-shelf buttressing (Fürst et al., 2016) and therefore ice flux at the grounding line and the stability of marine ice sheets (Schoof et al., 2017). It can also lead to a total disintegration of ice shelves followed by a potential marine ice-cliff instability (Pollard et al., 2015).

The nominal calving rate *C*_{r} (in m yr^{−1}) in the f.ETISh
model is evaluated with the following parameterisation (Pollard and DeConto, 2012a; Pattyn, 2017):

where ** v** is the vertical mean of the horizontal velocity,

*h*

_{e}is the sub-grid ice thickness within a fraction of the ice-edge grid cell that is occupied by ice (Pollard and DeConto, 2012a), Δ

*x*is the spatial resolution and ${w}_{\mathrm{c}}=min\left(\mathrm{1},{h}_{\mathrm{e}}/\mathrm{200}\right)$ is a weight factor.

We introduce uncertainty in calving by controlling the magnitude of
the calving rate with a scalar multiplier factor *F*_{calv}. This
approach is similar to Briggs et al. (2013) and Pollard et al. (2016). Here, we
consider for *F*_{calv} a nominal value of 1.0 and an uncertainty
range from 0.5 to 1.5; that is, we consider the calving rate to vary
between 50 % and 150 % of the nominal calving rate *C*_{r}.

### 2.2.5 Oceanic forcing

Sub-shelf melting is mainly controlled by sub-shelf ocean circulation, which can be affected by atmospheric changes. Ice-shelf thinning caused by increased sub-shelf melting leads to a reduction in ice-shelf buttressing. West Antarctica, where the bedrock lies mainly below sea level, is particularly vulnerable, as suggested by observational (Rignot et al., 2014) and modelling (Favier et al., 2014; Joughin et al., 2014) studies.

High melt rates at the base of ice shelves result from the inflow of relatively warm Circumpolar Deep Water in ice-shelf cavities (Hellmer et al., 2012; Schmidtko et al., 2014). Changes in ocean circulation resulting in stronger sub-ice-shelf circulation are expected to increase basal melt rates at the base of ice shelves (Jacobs et al., 2011; Hellmer et al., 2017). Increase in atmospheric temperature leading to the presence of warmer deep water on the continental shelf is expected to strengthen sub-ice-shelf circulation, thus leading to an increase in sub-shelf melting. Yet, the link between global climate change and changes in sub-shelf melting is not always clear: it has been suggested that future climate change could lead to an increase in sub-shelf melting (Hellmer et al., 2012; Timmermann and Hellmer, 2013), a positive meltwater feedback that enhances sub-ice-shelf circulation and can trigger a climate tipping point (Hellmer et al., 2017) and even to a decrease in sub-shelf melting (Dinniman et al., 2012). Golledge et al. (2019) have also highlighted the need to couple ice-sheet models to climate models as the increased discharge of freshwater from the Antarctic ice sheet could trap warm waters of the Southern Ocean below the sea surface.

Here, we capture the basic overturning circulation in ice-shelf cavities with
the PICO box model. In the PICO model, the strength of the overturning flux
is represented by a single parameter that depends on the density difference,
or equivalently on both the salinity and temperature differences, between the
incoming water masses on the continental shelf and the water masses near the
deep grounding line of the ice shelf. An increase in the ocean temperature on
the continental shelf leads to a stronger overturning flux and higher melt
rates at the base of ice shelves. The ocean temperature on the continental
shelf is determined from the present-day ocean temperature
${T}_{\mathrm{oc}}^{\mathrm{obs}}$ and the change in background atmospheric
temperature Δ*T* via the linear relationship in
Eq. (1).

We introduce uncertainty in sub-shelf melting by controlling the strength of
the overturning flux through uncertainty in the ocean temperature on the
continental shelf. Hence, we consider the ocean melt factor
*F*_{melt} as an uncertain parameter with a nominal value of 0.3
(Maris et al., 2014; Golledge et al., 2015) and an uncertainty range from 0.1 to 0.8.

### 2.2.6 Ice-shelf rheology

Ice rheology in large-scale ice-sheet models is usually described as an isotropic material obeying Glen's flow law (Greve and Blatter, 2009). However, ice is known to be an anisotropic material whose fabric is dependent on the temperature field, the strain-rate history and the stress history (Ma et al., 2010; Calonne et al., 2017). For a given fabric, the anisotropic response depends on the stress regime, which explains the ice stiffening when moving from a shear-dominated stress regime for grounded ice to an extension-dominated stress regime for floating ice.

We introduce an ice-shelf tune parameter *E*_{shelf} that accounts
for anisotropy between grounded and floating ice. A lower value makes ice
shelves more viscous (Briggs et al., 2013; Maris et al., 2014). We consider for
*E*_{shelf} a nominal value of 0.5 and an uncertainty range from 0.2
to 1, where a value of 0.5 means that the ice in the ice shelves is 2 times
more viscous than without shelf tuning.

### 2.2.7 Bedrock relaxation

Bedrock relaxation due to deglaciation may induce a negative feedback that promotes stability in marine portions and mitigates the effect of a marine ice-sheet instability (Gomez et al., 2010, 2013; Adhikari et al., 2014). The amplitude of the glacial isostatic uplift is determined by the flexural rigidity of the lithosphere and the viscous relaxation time of the asthenosphere. Recent studies (Van der Wal et al., 2015; Chen et al., 2018) have shown significant differences in the properties of the lithosphere and the asthenosphere between West and East Antarctica. Van der Wal et al. (2015) have found a lower viscosity and therefore a lower relaxation time for Earth's mantle underneath West Antarctica, making the glacial isostatic uplift in this region more sensitive to changes in ice thickness.

We account for the differences between East and West Antarctica by
introducing two characteristic relaxation times *τ*_{e} and
*τ*_{w} that we consider to both have a nominal value of
3000 years. We consider that *τ*_{e} has an uncertainty range
from 2500 to 5000 years and that *τ*_{w} has an uncertainty range
from 1000 to 3500 years.

## 2.3 Uncertainty quantification methods

### 2.3.1 Characterisation of uncertainty

To quantify the impact of uncertainties on the AIS response, we adopt a
probabilistic framework. Here, we assume, in the absence of any prior
information other than the aforementioned nominal values and minimal and
maximal values of the uncertainty ranges, that the parameters
*F*_{calv}, *F*_{melt}, *E*_{shelf},
*τ*_{e} and *τ*_{w} are uniform independent random
variables with bounds given by the minimal and maximal values of their
uncertainty ranges. We explore the uncertainty in the sliding exponent *m* by
considering its nominal and extremal values separately, thus allowing us to
reduce the dimension of the parameter space while being consistent with other
studies (Ritz et al., 2015). We explore the uncertainty in the atmospheric
forcing Δ*T* by considering the four RCP scenarios. We limit the
probabilistic characterisation to assuming uniform probability density
functions, and we do not address how this probabilistic characterisation
could be refined by using expert assessment, data and statistical methods
such as Bayesian inference. Yet, we provide results that give some insight
into the impact of the choice of this probabilistic characterisation later in
Sect. 3.6. We refer the reader to Petra et al. (2014),
Isaac et al. (2015), Ruckert et al. (2017), Gopalan et al. (2018) and
Conrad et al. (2018) for applications of Bayesian inference in glaciology and
to Aschwanden et al. (2016) and Gillet-Chaulet et al. (2016) for examples of a
calibration of the sliding exponent based on a comparison between simulated
and observed surface velocities that can be used to prescribe a probabilistic
characterisation of the sliding exponent.

### 2.3.2 Propagation of uncertainty

Given the probabilistic characterisation of the uncertainty in the parameters, the propagation of uncertainty serves to assess the impact of the uncertainty on the global mean sea-level change. In particular, its intent is to estimate the probability density functions for the change in GMSL as well as some of its statistical descriptors such as its mean, variance and quantiles. Various methods have been developed in UQ to estimate these statistical descriptors in a non-intrusive manner treating the ice-sheet model as a black box. Here, we use emulation methods based on a polynomial chaos (PC) expansion.

An emulator, also known as a surrogate model, is a computational model that mimics the ice-sheet model at low computational cost. Although emulators can also be obtained by Gaussian process regression (Rasmussen and Williams, 2006), we use polynomial chaos (PC) expansions (Ghanem and Spanos, 2003; Le Maître and Knio, 2010), which involve approximating the parameters-to-projection relationship as a polynomial in the parameters. We write this polynomial as a linear combination of polynomial basis functions and use least-squares regression (Appendix A) to evaluate the coefficients from a limited number of ice-sheet model runs at an ensemble of training points in the parameter space. The training points must be adequately chosen in the parameter space (Hadigol and Doostan, 2018), and the convergence of the PC expansion must be assessed to ensure accuracy. PC expansion may suffer from limitations: PC expansions require the parameters-to-projection relationship to be sufficiently smooth (no discontinuity or highly non-linear behaviour) to allow an efficient approximation as a low-degree polynomial, and PC expansions may be inefficient in high-dimensional problems.

The emulator of the relationship between the parameters and the projection is then used as a substitute for the ice-sheet model in a Monte Carlo method (Robert and Casella, 2013), in which samples of the parameters are drawn randomly from their probability density function and mapped through the emulator into corresponding samples of the projections. Approximations to the probability density function of the projection and its statistical descriptors are then obtained from these samples of the projections by using statistical estimation methods: for instance, the mean and quantiles are approximated with the sample mean and quantiles.

The use of an emulator has the following advantages: (i) it provides an inexpensive approximation of the ice-sheet model that accelerates UQ; (ii) it provides an explicit view of the relationship between the parameters and the projection, highlighting potential linear or non-linear dependences and interactions between the parameters; (iii) it allows efficient interpolation of the projections in the parameter space; (iv) it can be used to carry out stochastic sensitivity analysis to assess the influence of each parameter on the projections; (v) under certain conditions, the same emulator can be reused between UQ analyses with different probability density functions for the parameters; and (vi) it can be used for Bayesian calibration (Ruckert et al., 2017).

We consider an ensemble of 20 distinct model configurations given by each
combination of RCP scenario with a sliding law ($m=\mathrm{1},\mathrm{2},\mathrm{3}$ or 5) and each
combination of RCP scenario with the TGL parameterisation under the weakly
non-linear sliding law (*m*=2). For each model configuration, we built a
separate PC expansion to investigate the impact of uncertainty in the five
parameters *F*_{calv}, *F*_{melt}, *E*_{shelf},
*τ*_{e} and *τ*_{w} on the uncertainty in
ΔGMSL. Here, we assume the continental response of the Antarctic ice
sheet, as measured through ΔGMSL, to be sufficiently smooth to be
represented with a PC expansion. For each model configuration, we generated
an ensemble of 500 training points in the parameter space with a
maximin Latin hypercube
sampling design (Stein, 1987) and performed 500 forward simulations at
these training points. In total, we carried out 10 000 forward simulations
of the f.ETISh model. For each model configuration and time instant of
analysis, we fitted to the corresponding ensemble of 500 training points and
forward simulations a polynomial chaos expansion of degree 3, which we then
used as an emulator to evaluate the statistical descriptors, either directly
from the coefficients of the PC expansion or by running the emulator with an
ensemble of 10^{6} independent and identically distributed samples from the
parameter space. We present a convergence study and cross-validation for the
PC expansion in Appendix A. Finally, the confidence
regions for each model configuration are determined directly from the
corresponding ensemble of 500 training points and forward simulations without
the need for an emulator.

### 2.3.3 Stochastic sensitivity analysis

Stochastic sensitivity analysis serves to identify which uncertain parameters and their associated physical phenomenon are most influential in inducing uncertainty in the ice-sheet response. Here, we adopt the variance-based sensitivity indices (Saltelli et al., 2008), also called Sobol indices, described in more detail in Appendix B. Variance-based sensitivity indices rely on the decomposition of the variance of the projections as a sum of contributions from each uncertain parameter taken individually and an interaction term. Then, the Sobol index of a given uncertain parameter represents the fraction of the variance of the projections explained as stemming from this sole uncertain parameter. A Sobol index takes values between 0 and 1, whereby a value of 1 indicates that the entire variance of the projections is explained by this sole uncertain parameter and a value of 0 indicates that the uncertain parameter has no impact on the projection uncertainty.

We compute the Sobol indices directly from the PC coefficients (Crestaux et al., 2009; Le Gratiet et al., 2017).

### 2.3.4 Confidence regions for grounded-ice retreat

To gain insight into the impact of the uncertainty in determining which regions of Antarctica are most at risk of ungrounding, we construct confidence regions for grounded ice for several probability levels. We define a confidence region for grounded ice for a given probability level as a region of Antarctica that remains covered everywhere with grounded ice with a probability of at least the given probability level under the uncertainty introduced in the ice-sheet model (see Appendix C for the mathematical definition). The differences between these confidence regions for grounded ice for different probability levels provide insight into the risk of ungrounding (see Sect. 3.5). Such confidence regions are useful because confidence regions with boundaries far from the initial grounding line may indicate an important MISI and large differences between confidence regions for different probability levels indicate a significant impact of the uncertainty on the ice-sheet ungrounding. We construct these confidence regions for grounded ice based on an extension of previous work by Bolin and Lindgren (2015) for Gaussian random fields to our glaciological context.

We present nominal and probabilistic projections (relative to 2000 CE) for short-term (2100), medium-term (2300) and long-term (3000) timescales under different RCP scenarios and sliding laws.

## 3.1 Nominal projections

We first present nominal projections obtained using the nominal values of the parameters given in Table 1. We first present results under nominal conditions in order to assess subsequently the impact of uncertainties on AIS sea-level rise projections. Under nominal conditions, we find (Table 2) in RCP 2.6 an AIS contribution to sea level of 0.02 m by 2100, 0.07 m by 2300 and 0.20 m by 3000 and in RCP 8.5 an AIS contribution to sea level of 0.05 m by 2100, 0.59 m by 2300 and 3.90 m by 3000. In Fig. 3a, we represented the nominal projections as a function of time in all RCP scenarios. We find that the nominal AIS contribution to sea level is rather small in the first decades and starts to increase more significantly around 2100 with a rather constant growth rate. In Fig. 3b–e, we represented the nominal grounded ice region by 3000 for all RCP scenarios. We find that there is little ungrounding by 3000 in RCP 2.6 and RCP 4.5, while we observe a more significant ungrounding in Siple Coast and the Ronne basin in RCP 6.0 and a much more significant ungrounding in Siple Coast and the Amundsen Sea sector in RCP 8.5.

## 3.2 Parameters-to-projection relationship

One of the advantages of a polynomial chaos expansion is that it provides an explicit approximation to the parameters-to-projection relationship, which can be visualised to gain insight into the relationship between the parameters and the projections.

In Figs. 4 and 5, we used the PC expansions to visualise how the projections depend on each parameter individually (one-at-a-time) while keeping the other parameters fixed at their nominal value under the weakly non-linear sliding law in RCP 2.6 and RCP 8.5, respectively.

Figure 4a–c show that in RCP 2.6 ΔGMSL
increases with an increase in the calving factor non-linearly. The slope is
steeper for small values of this parameter, thus suggesting that ΔGMSL
is more sensitive to small changes about small values than about higher
values. Figure 4d–f indicate that ΔGMSL
increases rather linearly with an increase in the melt factor.
Figure 4g–i indicate a rather quadratic dependence
on the shelf anisotropy factor for short-term and medium-term projections,
with small and large values of this factor leading to more significant ice
loss than the nominal value. This quadratic dependence can be explained by
the influence of *E*_{shelf} on two competing processes: a higher
value of *E*_{shelf} softens the ice, thus leading to faster ice
flow in the ice shelves, but a higher value of *E*_{shelf} also
leads to ice-shelf thinning, thus reducing ice flux at the grounding line. In
addition, we find that ΔGMSL depends only little on the bedrock
relaxation times (Fig. 4j–o). In fact, lower bedrock
relaxation times do tend to stabilise the ice sheet and lower sea-level rise,
but this impact is weak compared to the influence of the other parameters.
This result may be explained by our orders of magnitude of the bedrock
relaxation times being of the order of a few millennia, thus preventing any
significant uplift in the next few centuries. Yet, a recent study by
Barletta et al. (2018) has suggested a bedrock relaxation timescale of the
order of decades to a century in the Amundsen Sea sector, thus making glacial
isostatic adjustment significant in the next decades and centuries in this
region. To assess the influence of shorter relaxation times, we performed
additional numerical experiments, albeit not reported in this paper, with a
relaxation time for the whole West Antarctic ice sheet that varies widely
from a few decades to a few millennia. For this range of values, we found
that bedrock relaxation has a more significant influence on the AIS response,
with a contribution to the uncertainty in ΔGMSL that can reach 10 %
in RCP 2.6 when *τ*_{w} is allowed to vary widely between 50 and
3500 years.

We find in RCP 8.5 similar trends. However, whereas *F*_{calv},
*F*_{melt} and *E*_{shelf} influence ΔGMSL equally
in RCP 2.6, the melt factor influences ΔGMSL most significantly in
RCP 8.5. Whereas Fig. 5d–e show that ΔGMSL
increases rather linearly with an increase in the melt factor,
Fig. 5f shows that ΔGMSL rather levels off at
a plateau for large values of *F*_{melt} for long-term projections
in RCP 8.5.

Additionally, Fig. S2 shows the emulators for several pairs of parameters with the other parameters fixed at their nominal values in RCP 2.6 and RCP 8.5. These figures show essentially the same trends as those identified in Figs. 4 and 5, in addition to interaction effects between the parameters. In particular, we find that smaller values of the calving and melt factors lead to mass gain (Fig. S2c), while larger values of the shelf anisotropy and melt factors lead to an important mass loss (Fig. S2d).

## 3.3 Sea-level rise projections

Under parametric uncertainties, we find (Table 3, Fig. 6) in RCP 2.6 a median AIS contribution to sea level that ranges from 0.02 to 0.03 m by 2100, from 0.07 to 0.13 m by 2300 and from 0.18 to 0.30 m by 3000 for the three sliding laws, as compared with the nominal projections of 0.02 m by 2100, 0.07 m by 2300 and 0.20 m by 3000; and we find in RCP 8.5 a median AIS contribution to sea level that ranges from 0.09 to 0.11 m by 2100, from 0.78 to 1.15 m by 2300 and from 3.18 to 6.12 m by 3000, as compared with the nominal projection of 0.05 m by 2100, 0.59 m by 2300 and 3.90 m by 3000. The median AIS sea-level rise projections are higher than the nominal projections, except for certain cases under the viscous sliding law. The nominal projections are not equal to the median projections because the ice-sheet model exhibits non-linearities, as illustrated in Figs. 4 and 5, and the probability density functions for the parameters are not symmetric about their nominal value. As in the nominal projections, we find that in all RCP scenarios and under all sliding laws, the median AIS contribution to sea level is rather small in the first decades and starts to increase around 2100. As in the nominal projections, the median AIS contribution to sea level shows a rather constant growth rate in RCP 2.6, RCP 4.5 and RCP 6.0; contrary to the nominal projection, the median AIS contribution to sea level exhibits initially an acceleration before subsequently exhibiting a deceleration in RCP 8.5 under both non-linear sliding laws (Fig. 6). This behaviour in the median projections in RCP 8.5 is a consequence of the initial rapid collapse of the West Antarctic ice sheet, which contributes 3–3.5 m to sea level, followed by a slower retreat in the East Antarctic ice sheet. By contrast, the nominal projections, which involve a melt factor of 0.3, indicate that the West Antarctic ice sheet does not completely collapse by the year 3000.

In RCP 2.6, we find (Table 3) 5 %–95 % probability
intervals for sea-level rise projections that range from −0.06 to 0.10 m
by 2100, from −0.14 to 0.31 m by 2300 and from −0.36 to 0.91 m by 3000.
In RCP 8.5, these probability intervals range from −0.03 to 0.23 m by
2100, 0.17 to 2.01 m by 2300 and 0.82 to 7.68 m by 3000. The nominal
projections are inside the 5 %–95 % probability intervals. We see that
the 5 %–95 % probability intervals indicate an increase in sea level,
though a decrease cannot be ruled out for the viscous sliding law and cooler
atmospheric conditions. Figure 6 highlights that the
33 %–66 % probability intervals become wider under warmer atmospheric
conditions and, to a lesser extent, more non-linear sliding conditions. The
uncertainty in the projections due to parametric uncertainty is rather
significant, with possible overlaps between the 5 %–95 % probability
intervals for different RCP scenarios. For instance,
Fig. 4i shows that a contribution to sea level of
about 0.7 m can be reached by 3000 in RCP 2.6 for the extreme value
*E*_{shelf}=1, while Fig. 5f shows that this
value is reached in RCP 8.5 for very limited sub-shelf melting
(*F*_{melt} of about 0.12). This suggests that projections with
similar contributions to sea level can arise in different RCP scenarios with
different combinations of parameter values. Measuring the relative dispersion
in ΔGMSL via the coefficients of variation, that is, the ratio between
the standard deviation and the mean value, we find a coefficient of variation
of 0.74 in RCP 2.6, 0.67 in RCP 4.5, 0.59 in RCP 6.0 and 0.33 in RCP 8.5 by
3000 under the weakly non-linear sliding law. This demonstrates that the
relative uncertainty in ΔGSLM projections is higher in cooler RCP
scenarios. This results from a much larger increase in the values of
ΔGSLM projections compared to the increase in their dispersion when
the RCP scenario gets warmer.

Figure 7 shows the probability density functions for the change
in GMSL at different timescales. The results display essentially unimodal
probability density functions with wider tails for warmer scenarios and
longer timescales. In RCP 2.6, the probability density functions resemble
Gaussian probability density functions. In RCP 8.5 and at the short and
medium timescales, the probability density functions are rather flat, which
can be explained by the dependence of ΔGMSL on *F*_{melt}
being rather linear (Fig. 5d and e). In RCP 8.5 and
at the long timescale, the probability density functions exhibit a more
localised mode at higher values of ΔGMSL, which can be explained by
the collapse of the West Antarctic ice sheet and thus the presence of the
plateau for higher values of the melt factor
(Fig. 5f).

Figure 8 gives the probability of exceeding the threshold values of 0.5, 1.0 and 1.5 m as a function of time. We find that the probability of exceeding 0.5 m by 2100 is negligible (probability of less than 1 %) in all RCP scenarios and under all sliding conditions. In RCP 2.6, the AIS contribution to sea level in the next centuries is strongly limited, with a probability of exceeding 0.5 m by 3000 reaching at most 30 %. In RCP 4.5 and RCP 6.0, we found nominal sea-level rise projections well below 1.5 m, while in RCP 8.5 this threshold is exceeded. In the presence of uncertainties, we find that the probability of exceeding 1.5 m of sea-level rise by 3000 can reach about 35 % in RCP 4.5, about 70 % in RCP 6.0 and about 95 % in RCP 8.5. The last result can be seen from Fig. 5, which indicates that ΔGMSL is below 1.5 m only in a small region of the parameter space associated with small values of the melt factor. Furthermore, the shape of the exceedance curves in Fig. 8 provides some insight into the uncertainty in the time when a certain threshold value is exceeded. The time when ΔGMSL exceeds 0.5 m with a probability of 33 % under the weakly non-linear sliding law is 2415 in RCP 4.5, 2270 in RCP 6.0 and 2185 in RCP 8.5, while this value is exceeded with a probability of 66 % at 2790 in RCP 4.5, 2430 in RCP 6.0 and 2245 in RCP 8.5. We find that nominal projections overestimate the time when ΔGMSL exceeds 0.5 m, with the exceedance time being beyond 3000 in RCP 4.5, around 2620 in RCP 6.0 and around 2280 in RCP 8.5.

## 3.4 Stochastic sensitivity analysis

Figure 9 provides the Sobol sensitivity indices for the change in GMSL on short-term, medium-term and long-term timescales in different RCP scenarios and under different sliding laws. We find that in RCP 2.6, the largest contribution to the uncertainty in ΔGMSL stems from the uncertainty in the ice-shelf rheology (Sobol indices ranging from 40 % to 60 %) followed by the uncertainty in the calving rate (Sobol indices ranging from 20 % to 40 %) and sub-shelf melting (Sobol indices ranging from 5 % to 25 %). Indeed, in RCP 2.6, sub-shelf melting plays only a limited role because ocean conditions remain essentially unchanged. Therefore, in RCP 2.6, the dispersion in ΔGMSL is mainly controlled by the ice-shelf rheology, which controls ice flow and buttressing in ice shelves, as well as calving, which reduces the extent of ice shelves and their buttressing.

By contrast, in warmer RCP scenarios and for longer timescales, the dominant
source of uncertainty becomes the uncertainty in sub-shelf melting, which
accounts in RCP 8.5 for more than 90 % of the uncertainty in sea-level rise
projections. As shown in Fig. 5, in RCP 8.5 at 3000
ΔGMSL varies by several metres over the range of values of
*F*_{melt}, while ΔGMSL varies only by a few tens of
centimetres over the range of values of the other parameters. Hence, the
dominant influence of the uncertainty in the melting factor is also a
consequence of the rather wide uncertainty range that we chose for this
parameter.

Finally, we find that, in all RCP scenarios and under all sliding laws, the uncertainty in the bedrock relaxation times for West and East Antarctica has a limited impact (Sobol index smaller than 1 %), which is a direct consequence of the very limited dependence of the projections on the bedrock relaxation times. Moreover, the interactions between the parameters have a negligible impact as the sums of the individual Sobol indices account almost entirely for the variances of the projections.

## 3.5 Projections of grounded-ice retreat

Figure 10 provides insight into the regions of Antarctica that are most at risk of ungrounding in different RCP scenarios and at different timescales under the weakly non-linear sliding law. Figure 10 was obtained as follows. First, we determined the 100 % confidence region for grounded ice, that is, the region of Antarctica where ice is certain to remain grounded, and we coloured it in grey. Thus, there is no risk that the grounding line will retreat to within the grey region. Then, we determined the 95 % confidence region for grounded ice, that is, the region of Antarctica that remains covered everywhere with grounded ice with a probability of more than 95 %, and we coloured the portion of the 95 % confidence region that extends beyond the 100 % confidence region in dark blue. Thus, there is a low risk of (less than) 5% that the grounding line will retreat to within the dark blue region. We continued this procedure for decreasing values of the confidence level and using different colours as indicated in the legend in Fig. 10.

We find that ice remains grounded in regions above sea level. By contrast, in all RCP scenarios, the risk of ungrounding is highest in marine sectors of West Antarctica with fast-flowing ice streams, especially in Siple Coast, in the Ronne basin, notably Ellsworth Land, and in the Amundsen Sea sector. In warm RCP scenarios and at longer timescales, we also observe a risk of grounding-line retreat in the Wilkes marine basin in East Antarctica, where the grounding line could retreat between 100 km (with a risk of 95 %) and 500 km (with a risk of 5 %) from its present-day position. The risk of retreat in Wilkes basin may partially explain the acceleration in sea-level rise that we observed in Fig. 6 in RCP 8.5. The risk of grounding-line retreat in the Antarctic Peninsula is very limited due to the bedrock topography being above sea level, the marine glaciers being small and a high increase in precipitation in this region.

In RCP 2.6, we observe that the grounding line is quite stable over the next millennium, with the 100 % confidence region for grounded ice being almost unchanged from the present-day grounded ice region (the 100 % confidence region for grounded ice by 3000 only differs from the present-day grounded ice region by a few tens of kilometres). In RCP 4.5, ice remains grounded in most of the West Antarctic ice sheet over the next centuries, but our results also suggest a risk of retreat of the grounding line in some sectors of West Antarctica on longer timescales. In RCP 6.0, we find that the West Antarctic ice sheet belongs by 3000 to the 66 % confidence region for grounded ice, while in RCP 8.5, it belongs by 3000 only to the 5 % confidence region. This suggests a risk of 33 % that a major collapse of the West Antarctic ice sheet might occur in RCP 6.0 by 3000 and a risk of 95 % that a major collapse of the West Antarctic ice sheet might occur in RCP 8.5 by 3000.

As compared with the nominal projections of a limited retreat of the grounding line in West Antarctica in RCP 6.0 by 3000, we find that the impact of the parametric uncertainty is that a complete collapse of the West Antarctic ice sheet may occur with a risk of 33 % in RCP 6.0 by 3000. Moreover, Fig. 3e suggests that a complete disintegration of the West Antarctic ice sheet is underway by 3000 in RCP 8.5, while in Fig. 10l the West Antarctic ice sheet has already collapsed by 3000.

Additionally, we compared the projections under the weakly non-linear sliding law with projections under the other sliding laws, that is, the viscous sliding law (Fig. S7) and the strongly non-linear sliding law (Fig. S8). We find a lower risk of retreat of the grounding line under the viscous sliding law with a slower disintegration of the West Antarctic ice sheet compared to the other sliding laws, especially in the drainage basins of Thwaites and Pine Island glaciers, which belong to the 50 % confidence region for grounded ice by 3000 in RCP 8.5, while they belong to the 5 % confidence region by 3000 in RCP 8.5 for the other sliding laws. The strongly non-linear sliding law seems to favour a faster and deeper retreat of the grounding line, especially in the marine sectors of East Antarctica and the drainage basins of Thwaites and Pine Island glaciers. However, Figs. 10i and S8i suggest that ungrounding may be less significant in Siple Coast under the strongly non-linear sliding law than the weakly non-linear sliding law. Actually, MISI may occur when driving stresses overcome resistive stresses (Waibel et al., 2018). Driving stresses are primarily determined by the surface slopes, while resistive stresses depend on the basal sliding coefficient and the ice velocity through the sliding exponent. Driving stresses at the grounding line are higher in the Amundsen Sea sector than in Siple Coast due to steeper surface slopes, leading to a greater sensitivity of the grounding line in the former region to non-linearity and a more plastic response.

## 3.6 Influence of the parameter probability density function

So far, we have represented the uncertain parameters *F*_{calv},
*F*_{melt}, *E*_{shelf}, *τ*_{e} and
*τ*_{w} by a uniform distribution on a fixed support. For
instance, the uncertain parameter *F*_{calv} is represented by a
uniform distribution with support
$[{F}_{\mathrm{calv},\mathrm{min}},{F}_{\mathrm{calv},\mathrm{max}}]$, where
*F*_{calv,min} and *F*_{calv,max} are the minimum and
maximum values in Table 1. We now address the influence of
the probabilistic characterisation of the parametric uncertainty on our
probabilistic projections by controlling the size of these supports. Hence,
we represent the uncertain parameters *F*_{calv},
*F*_{melt}, *E*_{shelf}, *τ*_{e} and
*τ*_{w} by a family of uniform probability density functions
indexed by a unique scaling factor $\mathit{\alpha}\in [\mathrm{0},\mathrm{1}]$ that controls the
supports of the uniform probability density functions. For instance, the
uncertain parameter *F*_{calv} is represented for a given value of
*α* by a uniform distribution with support
$[{F}_{\mathrm{calv},\mathrm{nom}}+\mathit{\alpha}({F}_{\mathrm{calv},\mathrm{min}}-{F}_{\mathrm{calv},\mathrm{nom}}),{F}_{\mathrm{calv},\mathrm{nom}}+\mathit{\alpha}({F}_{\mathrm{calv},\mathrm{max}}-{F}_{\mathrm{calv},\mathrm{nom}}\left)\right]$, where
*F*_{calv,nom}, *F*_{calv,min} and *F*_{calv,max}
are the nominal, minimum and maximum values in Table 1. For
the other parameters, the supports are defined similarly. The value
*α*=0 represents the ice-sheet model without parametric uncertainty,
that is, the ice-sheet model for the nominal values of the parameters, and
the value *α*=1 represents the full uncertainty ranges considered so
far. We propagated the uncertainty from the parameters to the sea-level rise
projections for different values of the scaling factor, reusing the emulator
that we had built for *α*=1 over the whole parameter space. Hence, the
projections to follow for *α*=0 based on this emulator are not exactly
equal to the nominal projections determined directly from the ice-sheet
model.

For the long-term projections, Fig. 11a–d show the median and
the 33 %–66 % and 5 %–95 % probability intervals as a function of
the scaling factor in the different RCP scenarios under the weakly non-linear
sliding law and Fig. 11e–h show the probability density
functions for the values $\mathit{\alpha}=\mathrm{0.2},\mathrm{0.6}$ and 1.0. We find that the median
projections increase with an increase in the scaling factor and range in
RCP 2.6 from 0.16 to 0.30 m, in RCP 4.5 from 0.37 to 0.93 m, in RCP 6.0
from 1.11 to 2.16 m and in RCP 8.5 from 3.80 to 5.08 m. In addition, the
width of the probability intervals increases with increasing uncertainty in
the parameters. While the width of the 33 %–66 % probability interval
increases rather linearly with the scaling factor and is rather symmetric
about the median, the width of the 5 %–95 % probability interval
increases more non-linearly with an increase in the scaling factor, as
illustrated in Fig. 11a in RCP 2.6 and Fig. 11d
in RCP 8.5. The probability density functions attribute higher weight to
larger values of ΔGMSL under increased parametric uncertainty, with
increased weight given to larger values of in particular *E*_{shelf}
in RCP 2.6 (Fig. 4i) and *F*_{melt} in
RCP 8.5 (Fig. 5f). Figure 11a–d show
that the sensitivity of the amount of uncertainty in the projections
(probability intervals) to the amount of uncertainty in the parameters
(scaling factor) is higher for warmer scenarios, with an upper bound between
RCP 6.0 and RCP 8.5 as a consequence of the collapse of the West Antarctic
ice sheet.

## 3.7 Projections under a more plastic sliding law

We ran the same ensemble of simulations under a more plastic sliding law
(*m*=5). Table 4 and Fig. 12 give for
the AIS contribution to sea level the median and the 33 %–66 % and
5 %–95 % probability intervals. As compared with less plastic sliding
laws ($m=\mathrm{1},\mathrm{2},\mathrm{3}$), we find an increase in sea-level rise projections on
short-term and medium-term timescales, thus suggesting a more significant and
faster response to perturbations under more plastic sliding conditions. On a
long-term timescale, the ice-sheet mass loss can be less important under
*m*=5 than *m*=3, as observed in the median projections by 3000 in RCP 4.5
and RCP 8.5.

We find that overall the AIS contribution to sea level is an increasing
function of the sliding exponent, with the differences between successive
exponents getting smaller as *m* increases, as already pointed out by
Gillet-Chaulet et al. (2016); for instance, we observe a greater difference in
the projections between *m*=1 and *m*=2 than between *m*=3 and *m*=5. On
longer timescales, tipping points and non-linearities associated with MISI
may trigger a slightly different response of the ice sheet depending on the
initial conditions, which could explain the smaller ice loss in our results
under *m*=5 than *m*=3.

## 3.8 TGL parameterisation

We ran the same ensemble of simulations under the weakly non-linear sliding law using, this time, the TGL parameterisation instead of the SGL parameterisation. Under the TGL parameterisation, Table 5 and Fig. 13 give for the AIS contribution to sea level the median and the 33 %–66 % and 5 %–95 % probability intervals. We find an overall increase in sea-level rise projections compared to our results under the SGL parameterisation. This result was expected as the TGL parameterisation has been shown to increase grounding-line sensitivity to environmental changes (Tsai et al., 2015; Pattyn, 2017). The probability of exceeding 0.5 m by 2100 is still negligible (probability of less than 1 %) in all RCP scenarios. However, the probability of exceeding 0.5 and 1 m by 3000 in RCP 2.6 can reach more than 40 % and 10 %, respectively. For other scenarios, the retreat of the Antarctic ice sheet is much more pronounced and faster than under the SGL parameterisation, and the probability of exceeding 1.5 m of sea-level rise by 3000 can reach more than 50 % in RCP 4.5, 80 % in RCP 6.0 and 99 % in RCP 8.5.

Figure 10m–o show the confidence regions for grounded ice under the TGL parameterisation in RCP 8.5. As also pointed out by Pattyn (2017), the TGL parameterisation leads to a faster and more significant grounding-line retreat in the marine sectors and an additional mass loss from East Antarctica, especially in the Aurora basin.

## 4.1 Comparison of the sea-level rise projections with previous work

Regarding the short-term AIS contribution to sea level, we projected in the RCP 2.6 scenario a median of 0.02–0.03 m under the SGL parameterisation and 5 %–95 % probability intervals from −0.06 to 0.10 m. These projections are similar to other estimates based on other mechanisms. Golledge et al. (2015) reported an AIS contribution to sea level between −0.01 and 0.10 m with the lower and higher bounds corresponding to the absence and presence of sub-grid interpolation of basal melting at the grounding line. DeConto and Pollard (2016) reported an AIS contribution to sea level between −0.11 and 0.15 m based on a model calibration with a range of Pliocene sea-level targets between 5 and 15 m higher than today. In the same scenario, we found an increased AIS response under the TGL parameterisation with a median AIS contribution to sea level of 0.06 m and a 5 %–95 % probability interval between 0.00 and 0.14 m by 2100. These higher projections are similar to the higher projections (0–0.22 m) reported by DeConto and Pollard (2016) for a higher range of Pliocene sea-level targets between 10 and 20 m. In all RCP scenarios, under all sliding laws and under both grounding-line parameterisations, our results suggested that the AIS contribution to sea level does not exceed 0.5 m by 2100, with a probability of at least 99 %. These results are in agreement with Ritz et al. (2015), who determined an AIS contribution to sea level that lies between 0.05 and 0.30 m, and Golledge et al. (2015), who found an AIS contribution to sea level that reaches at most 0.38 m in RCP 8.5 with sub-grid interpolation of basal melting at the grounding line. Yet, projections by DeConto and Pollard (2016) and Schlegel et al. (2018), who applied the more sensitive Budd-type friction law (Brondex et al., 2017), can exceed 0.5 and even 1 m, but these higher projections are under extreme and maybe unrealistic warming conditions.

Regarding the long-term AIS contribution to sea level, our projections under the SGL parameterisation are similar to other estimates by Golledge et al. (2015). In particular, both studies suggest that the AIS contribution to sea level by 3000 in RCP 2.6 is limited to less than 1 m with a probability estimated to be at least 95 % in our study, while an AIS contribution to sea level above 1.5 m by 3000 may arise in all other RCP scenarios. Yet, our long-term projections are generally below projections by DeConto and Pollard (2016) with hydrofracturing and ice-cliff failure mechanisms especially under warmer RCP scenarios, but the discrepancies between the projections of both models are reduced under the TGL parameterisation. In RCP 2.6, our projections by 2500 under the TGL parameterisation range from 0.04 to 0.73 m, which is similar to projections by DeConto and Pollard (2016), which range respectively from −0.23 to 0.61 m and 0.02 to 0.48 m for the lower and higher ranges of Pliocene of sea-level targets.

## 4.2 Comparison of the impact of parametric uncertainty with previous work

Similarly to Golledge et al. (2015, 2017), our study emphasised the pivotal role played by the emission scenario and sub-shelf melting as critical drivers of the future changes in the AIS mass balance on medium-term and long-term timescales through ice-shelf thinning and subsequent reduced buttressing. As in Ritz et al. (2015), we found that the AIS contribution to sea level is an increasing function of the sliding exponent, thus meaning that more plastic sliding conditions speed up the ice flow and consequently ice loss. Following Pattyn (2017), we highlighted the greater sensitivity of the grounding-line migration under the TGL parameterisation, thus stressing the key role played by physical processes in the vicinity of the grounding line.

## 4.3 Comparison of projections of grounded-ice retreat with previous work

Consistent with our results, Golledge et al. (2015), with a 10 km resolution model, projected that grounding-line retreat is most significant in the Siple Coast region. However, Ritz et al. (2015), based on the probability of MISI onset, as well as Cornford et al. (2015), with a sub-kilometre resolution around the grounding line, projected that grounding-line retreat is most significant in the Amundsen Sea sector. Schlegel et al. (2018) found that grounding-line retreat is most significant in the Amundsen Sea sector under generalised ocean warming experiments for the Antarctic ice sheet, but, after calibrating sub-shelf melt rates with bounds that vary region by region and are assigned values deduced from the literature and model sensitivity studies, they found that the western Ronne basin has the greater sensitivity. These discrepancies between our findings and those of Cornford et al. (2015), Ritz et al. (2015) and Schlegel et al. (2018) may be explained by our ocean model which may overestimate sub-shelf melting underneath the Ross ice shelf and underestimate ocean circulation in the Amundsen Sea and by our initialisation method which may underestimate the basal sliding coefficients for Thwaites and Pine Island glaciers. Moreover, the lower sensitivity of the Amundsen Sea sector may arise in our simulations from shortcomings in the buttressing parameterisation, our low resolution not capturing properly the bedrock topography, the small pinning points and the flow dynamics in the narrower sectors of the ice sheet and our representation of calving, which may increase the instability threshold; see, for instance, Arthern and Williams (2017) and Waibel et al. (2018) for more thorough discussions about the instability threshold in the Amundsen Sea sector.

## 4.4 Projections of ice loss and grounding-line retreat under parametric uncertainty

The significance of the contribution of the Antarctic ice sheet to sea level under climate change is primarily controlled by the sensitivity, the response time and the vulnerability of its marine drainage basins, with the West Antarctic ice sheet more sensitive and vulnerable than the East Antarctic ice sheet. The instability of marine drainage basins and their ability to trigger accelerated ice loss and significant grounding-line retreat is determined by bedrock topography and ice-shelf buttressing which depends on the importance of ice-shelf thinning. Our nominal projections showed that the AIS contribution to sea level by 3000 is rather limited (less than 1 m) in RCP 2.6, RCP 4.5 and RCP 6.0, while an accelerated ice loss that leads to a contribution of several metres is triggered in RCP 8.5 (Fig. 3a). In addition, the nominal retreat of the grounding line by 3000 is rather limited in RCP 2.6, RCP 4.5 and RCP 6.0 (Fig. 3b–d), while a significant retreat of the grounding line is triggered in the Siple Coast, Ronne–Filchner and Amundsen Sea sectors in RCP 8.5 (Fig. 3e).

Our probabilistic results provide insight into the impact of parametric uncertainty on these projections. In RCP 2.6, the projections hold irrespective of parametric uncertainty: the AIS contribution to sea level by 3000 has a 95 % quantile of up to 0.91 m (Table 3) and there is a limited risk that the grounding line will retreat beyond our nominal projections (Figs. 3b and 10c). In RCP 4.5 and RCP 6.0, the projections are more sensitive to parametric uncertainty than in RCP 2.6: the AIS contribution to sea level by 3000 has a 95 % quantile of up to 2.56 m in RCP 4.5 and up to 5.17 m in RCP 6.0 (Table 3), and both scenarios entail a risk of triggering a more significant retreat of the grounding line beyond our nominal projections (Fig. 3c–d). This risk is present especially in the Amundsen Sea sector, and it is less significant in RCP 4.5 than in RCP 6.0, in which a complete disintegration of the West Antarctic ice sheet may be triggered (Fig. 10f and i). Finally, in RCP 8.5, our probabilistic results suggest an accelerated ice loss and a significant retreat of the grounding line in West Antarctica, as in our nominal projections: the AIS contribution to sea level by 3000 has an uncertainty range with a 5 % quantile above 0.82 m and a 95 % quantile of up to 7.68 m (Table 3), and there is a high risk of triggering a complete disintegration of the West Antarctic ice sheet (Fig. 10l).

In conclusion, the projections hold irrespective of parametric uncertainty in the strongly mitigated RCP 2.6 scenario: accommodating parametric uncertainty in the ice-sheet model leads to projections in agreement with the nominal projections of limited ice loss and limited grounding-line retreat in RCP 2.6. However, the projections are more sensitive to parametric uncertainty for intermediate scenarios such as RCP 4.5 and RCP 6.0: accommodating parametric uncertainty in the ice-sheet model leads to projections in disagreement with the nominal projections and indicates instead some risk of triggering accelerated ice loss and significant grounding-line retreat for intermediate scenarios such as RCP 4.5 and RCP 6.0. Finally, the warm RCP 8.5 scenario triggers the collapse of the West Antarctic ice sheet, almost irrespective of parametric uncertainty.

## 4.5 Structural uncertainty and limitations

A first limitation of our study is associated with the modelling hypotheses inherent to our ice-sheet model. The f.ETISh model is an ice-sheet model that focuses on essential marine ice-sheet mechanisms, similarly to the ice-sheet model by Pollard and DeConto (2012a). Certain physical processes, especially small-scale processes, may be represented imperfectly, especially with the 20 km resolution adopted for our simulations. This may reduce the ability to simulate important ice streams such as Pine Island and Thwaites glaciers (only represented by a few grid points), whose stability is controlled by local bedrock features (Waibel et al., 2018). Hence, grounding-line migration and thresholds for instabilities may not be captured properly even with a parameterisation of the grounding-line flux, especially in low-forcing scenarios and short-term projections. As discussed in Sect. 2.1, the flux condition at the grounding line is also questionable for buttressed ice shelves (Pattyn et al., 2013; Reese et al., 2018a) found around Antarctica and short transients (Drouet et al., 2013; Pattyn et al., 2013; Pattyn and Durand, 2013). Yet, we think that using a 20 km resolution and a flux condition at the grounding-line remains an acceptable assumption in large-scale and long-term ice-sheet simulations and large-ensemble simulations. We expect discrepancies between our results and results at a higher spatial resolution or with a higher level of complexity to be limited when compared to the uncertainty in the results due to the uncertainty in the model parameters and forcing. Besides, sub-shelf melting may not be captured properly despite the use of the PICO ocean-model coupler, especially in the Amundsen Sea sector and underneath the Ross ice shelf (Timmermann et al., 2012; Depoorter et al., 2013; Rignot et al., 2013; Moholdt et al., 2014). A second limitation comes from the hypotheses relevant to our characterisation of uncertainties. We adopted rather large uncertainty ranges for the parameters. As discussed in Sect. 3.6, projections can be strongly affected by extreme conditions. Moreover, we chose the bounds of the uncertainty ranges quite heuristically. A third limitation comes from the fairly simple way in which certain sources of uncertainty were introduced. We assumed a direct influence of the atmospheric forcing on sub-shelf melting through the ocean melt factor. However, as discussed in Sect. 2.2.5, atmospheric forcing affects the Circumpolar Deep Water circulation in ice-shelf cavities and subsequently modifies sub-shelf melting (Pritchard et al., 2012). Still, the fate of the Southern Ocean and the evolution of sub-shelf melting under global climate change remains unclear (Hellmer et al., 2012; Dinniman et al., 2012; Timmermann and Hellmer, 2013; Hellmer et al., 2017). Given the importance of sub-shelf melting in driving the future response of the Antarctic ice sheet, there is a clear need to better constrain future sub-shelf melt and incorporate it properly into ice-sheet models. We also introduced uncertainty in calving with a simple multiplier factor that does not take into account the stress regime in the ice shelves. A fourth limitation concerns the correction for the initial drift that adds a bias to the projections. We adopted this correction to make our results insensitive to model initialisation, but our corrected results are likely to underestimate the AIS response on a short-term timescale, as our approach does not take into account any current transient changes in the Antarctic ice sheet. Another limitation is related to the construction of our emulator. To avoid overfitting the training data, we tuned the emulator to reproduce the overall trend in the parameters-to-projection relationship, but not local variations in the parameter space that may stem from numerical noise and errors.

We studied the multi-centennial response of the Antarctic ice sheet under uncertainty using methods from the field of UQ. We investigated uncertainties in atmospheric forcing, basal sliding, grounding-line parameterisation, sub-shelf melting, calving, ice-shelf rheology and bedrock relaxation. We used emulation-based methods to represent the parameters-to-projection relationship, stochastic sensitivity analysis to assess the significance of each source of uncertainty in inducing uncertainty in the projections and confidence regions for excursion sets to assess the risk of grounding-line retreat. We found that all investigated sources of uncertainty, except bedrock relaxation time, contribute to the uncertainty in the projections. We showed that the sensitivity of the projections to uncertainties increases and the contribution of the uncertainty in sub-shelf melting to the uncertainty in the projections becomes more and more dominant as atmospheric and oceanic temperatures rise, with a contribution to the uncertainty in sea-level rise projections that goes from 5 % to 25 % in RCP 2.6 to more than 90 % in RCP 8.5. We showed that the significance of the AIS contribution to sea level is controlled by MISI in marine basins, with the biggest contribution stemming from the more vulnerable West Antarctic ice sheet. We found that, irrespective of parametric uncertainty, the strongly mitigated RCP 2.6 scenario prevents the collapse of the West Antarctic ice sheet, that in both the RCP 4.5 and RCP 6.0 scenarios the occurrence of MISI in marine basins is more sensitive to parametric uncertainty, and that, almost irrespective of parametric uncertainty, RCP 8.5 triggers the collapse of the West Antarctic ice sheet.

All datasets used in this paper are publicly available, including Bedmap2 (Fretwell et al., 2013), ocean temperature and salinity (Schmidtko et al., 2014) and geothermal heat flux (An et al., 2015). Results of the RACMO2 model were kindly provided by Melchior Van Wessen. Results from the f.ETISh model for this study are available on request from Kevin Bulthuis (kevin.bulthuis@uliege.be). The MATLAB code used to analyse the results is also available on request from the same author.

In this Appendix, we concisely provide further details about how we used PC expansions in our study; we refer the reader to, for instance, Ghanem et al. (2017) and Le Maître and Knio (2010) for more comprehensive treatments of the theory and various applications of PC expansions.

Let us represent the ice-sheet model as an abstract model
*y*=*g*(** x**)
where $\mathit{x}=({x}_{\mathrm{1}},\mathrm{\dots},{x}_{d})$ is a vector of

*d*parameters,

*y*the model response and

*g*the response function. In our study,

*d*=5,

*x*

_{1}=

*F*

_{calv},

*x*

_{2}=

*F*

_{melt},

*x*

_{3}=

*E*

_{shelf},

*x*

_{4}=

*τ*

_{e},

*x*

_{5}=

*τ*

_{w}and

*y*=ΔGMSL at a given time. In our study, the parameters are uncertain and have a probability density function that we denote by

*p*.

A polynomial chaos expansion is an approximation of the response function *g*
with a polynomial *g*^{K} as

where the *ψ*_{k} are a basis of predefined polynomials of increasing degree
and orthonormal with respect to the probability density function of the
parameters, by which we understand that

and *K*+1 is the number of predefined polynomials in the expansion.

In order to fit the PC expansion in Eq. (A1) to the ice-sheet model, we calculate the coefficients using a (weighted) least-squares approach:

where {*x*^{(i)},1*⩽**i**⩽**N*} is a set of *N*
training points in the parameter space,
$\left\{{y}^{\left(i\right)}=g\left({\mathit{x}}^{\left(i\right)}\right),\mathrm{1}\mathit{\u2a7d}i\mathit{\u2a7d}N\right\}$ is the set
of model responses at the training points,
$\mathit{c}=({c}_{\mathrm{0}},{c}_{\mathrm{1}},\mathrm{\dots},{c}_{K})$ collects the PC coefficients and
{*w*^{(i)},1*⩽**i**⩽**N*} is the set of weights. We
normalise the parameters to accommodate the different orders of magnitude of
the parameters and reduce the potential ill-conditioning of the least-squares
problem. We generate the set of training points with a maximin Latin
hypercube sampling design (Stein, 1987; Fajraoui et al., 2017; Hadigol and Doostan, 2018), which
is a space-filling design that aims at maximising the smallest distance
between neighbouring points, thus ensuring a proper coverage of the parameter
space. We consider a PC expansion of degree 3, which corresponds to *K*=56
for *d*=5, as we strive to reproduce the overall trend of the response
function without overfitting small local variations.

We solve Eq. (A3) by solving the normal equations

where $\mathit{y}=({y}^{\left(\mathrm{1}\right)},\mathrm{\dots},{y}^{\left(N\right)})$ collects the model responses at the
training points, **G** is the measurement matrix whose entries are
given by **G**_{ki}=*ψ*_{k}(*x*^{(i)}) and **W** is a
diagonal weight matrix whose diagonal entries are the weights
(**W**_{ii}=*w*^{(i)}). As we construct the training points by using a
Latin hypercube sampling design, the training points *x*^{(i)} have
equal weights of ${\mathbf{W}}_{ii}={w}^{\left(i\right)}=\mathrm{1}/N$. We do note that there exist
other methods to determine the coefficients in
Eq. (A1), including methods involving deterministic
quadrature rules (Le Maître and Knio, 2010). Here, one of our motivations for
choosing a least-squares approach is its good ability to handle numerical
noise in low-degree expansions (Iskandarani et al., 2016).

We estimate statistical descriptors of the uncertain model response with
Monte Carlo simulation in which the PC expansion is used as a computationally
efficient substitute for the ice-sheet model. For instance, we estimate the
probability density function of the response through kernel density
estimation (Scott, 2015), while we approximate the mean *μ* and the
variance *σ*^{2} of the model response as

where {*x*^{(i)},1*⩽**i**⩽**ν*} now denotes
an ensemble of *ν* independent and identically distributed samples from the
probability density function of the parameters. Because the PC expansion is
based on orthonormal polynomials with respect to the probability density
function of the uncertain parameters and of increasing degree with
*ψ*_{0}=1, we can evaluate some of the statistical descriptors of the
uncertain model response directly from the PC coefficients. For example, the
mean *μ* and the variance *σ*^{2} of the model response can also be
approximated as follows:

The accuracy of the PC expansion has to be assessed with respect to the degree of the PC expansion and the number of training points. We validate the accuracy of the PC expansion using cross-validation and convergence tests. We generate a new set of samples in the parameter space and we compare the exact response of the ice-sheet model with the approximate response of the PC expansion. Figure A1a–c show results of such a cross-validation for a PC expansion of degree 3. These figures suggest that the PC expansion represents the overall model response with sufficient accuracy. Lower accuracy is achieved near the boundaries of the parameter space. We also carry out convergence tests, as illustrated in Fig. A1d in RCP 8.5 and at time 3000. This figure represents the maximum absolute error and the mean-squared error between the exact response of the ice-sheet model and the approximate response of the PC expansion at the training points as a function of the number of training points. We see that for the 500 training points considered in this paper, reasonable convergence of the PC expansion is achieved.

In this Appendix, we concisely provide further details about how we used stochastic sensitivity analysis; we refer the reader to, for instance, Ghanem et al. (2017) and Saltelli et al. (2008) for more comprehensive treatments of the theory and various applications.

We first assume that the uncertain parameters are statistically independent.
Sobol indices are based on the decomposition of the response function *g* in
terms of the main effects associated with the parameters individually and an
interaction effect associated with all parameters together:

where *g*_{0} is constant, each *g*_{i} the so-called main effect associated with
the corresponding parameter *x*_{i} and *g*_{I} the so-called
interaction effect. The constant, the main effects and the interaction effect
are given by

where ${\mathit{x}}_{-i}=({x}_{\mathrm{1}},\mathrm{\dots},{x}_{i-\mathrm{1}},{x}_{i+\mathrm{1}},\mathrm{\dots},{x}_{d})$ denotes the
subset of parameters including all the parameters except *x*_{i} and *p*_{−i}
denotes the probability density function for this subset of parameters. The
constant *g*_{0} is the mean value of *g*. Using the calculus of variations, it
can be shown that the main effect *g*_{i} is such that the function *g*_{0}+*g*_{i}
is the least-squares best approximation of *g* among all functions that
depend only on *x*_{i}. As a consequence of this least-squares best
approximation property, the functions ${g}_{\mathrm{0}},{g}_{\mathrm{1}},\mathrm{\dots},{g}_{\mathrm{d}},{g}_{\mathrm{I}}$ are
orthonormal with respect to the probability density function of the uncertain
parameters.

As a consequence of the orthonormality of the functions
${g}_{\mathrm{0}},{g}_{\mathrm{1}},\mathrm{\dots},{g}_{\mathrm{d}},{g}_{\mathrm{I}}$, the variance *σ*^{2} of the model
response can be decomposed as

with

where *p*_{i} denotes the probability density function of *x*_{i}. Here, *S*_{i},
which takes a value between 0 and 1, is the Sobol index for the *i*th
parameter and *S*_{I} is the interaction index. The Sobol index *S*_{i}
can be interpreted either as the relative contribution of the uncertainty in
the sole *i*th parameter to the variance of the model response or as the
reduction in the variance of the model response that we may expect by
learning the exact value of this parameter (Oakley and O'Hagan, 2004). Sobol indices
allow us to rank the uncertain parameters in terms of their contribution to
the variance of the model response, thus indicating which parameters are most
influential in inducing the uncertainty in the model response.

Here, we substituted the ice-sheet model by a PC expansion based on orthonormal polynomials and estimated the Sobol indices directly from the PC coefficients (Crestaux et al., 2009), that is,

where *A*_{i} is the set of indices associated with the non-constant
polynomials that only depend on *x*_{i}.

In this Appendix, we concisely provide further details about how we defined and computed confidence regions for grounded ice. To distinguish between grounded ice and floating ice at a given time, the f.ETISh model (Pattyn, 2017) evaluates the so-called buoyancy imbalance:

where *b* is the bedrock elevation, *h* the ice thickness, *ρ*_{w}
the water density and *ρ*_{i} the ice density. The buoyancy
imbalance is negative for floating ice, positive for grounded ice and null at
the grounding line. Therefore, the grounded ice domain *D*_{g} is

In the presence of uncertainties in the model, we define an *α* %
confidence region *D*_{g}(*α*) for the grounded ice domain as a
region of Antarctica that has a probability of at least *α* of being
included in the grounded ice domain, that is,

We compute the confidence regions using an adaptation of a thresholding algorithm by Bolin and Lindgren (2015). We seek the confidence regions in a parametric family indexed by a threshold parameter and based on the marginal probability density functions, and we determine the threshold parameter so as to achieve the required level of confidence. While Bolin and Lindgren (2015) consider Gaussian random fields, here we work with non-Gaussian random fields, which requires us to evaluate the marginal probability density functions and the probability of inclusion with Monte Carlo simulation.

The supplement related to this article is available online at: https://doi.org/10.5194/tc-13-1349-2019-supplement.

All the authors discussed the results presented in this paper. KB conducted the design, execution and UQ analysis of the experiments, with relevant inputs from MA for the UQ methodology and SS and FP for the physical interpretation of the results. The paper was written by KB and MA with relevant comments from all the co-authors.

The authors declare that they have no conflict of interest.

We would like to thank Andy Aschwanden and one anonymous referee for their very helpful comments that helped improve the overall quality and readability of the manuscript. Kevin Bulthuis would like to acknowledge Andreas Wernecke for personal communication about the manuscript. Kevin Bulthuis would like to acknowledge the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) for its financial support (F.R.S.-FNRS Research Fellowship). Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S-FNRS) under grant no. 2.5020.11.

This paper was edited by Olivier Gagliardini and reviewed by Andy Aschwanden and one anonymous referee.

Adhikari, S., Ivins, E. R., Larour, E., Seroussi, H., Morlighem, M., and Nowicki, S.: Future Antarctic bed topography and its implications for ice sheet dynamics, Solid Earth, 5, 569–584, https://doi.org/10.5194/se-5-569-2014, 2014. a

An, M., Wiens, D. A., Zhao, Y., Feng, M., Nyblade, A., Kanao, M., Li, Y., Maggi, A., and Lévêque, J.-J.: Temperature, lithosphere-asthenosphere boundary, and heat flux beneath the Antarctic Plate inferred from seismic velocities, J. Geophys. Res.-Sol. Ea., 120, 8720–8742, https://doi.org/10.1002/2015jb011917, 2015. a, b

Arnst, M. and Ponthot, J.-P.: An overview of nonintrusive characterization, propagation, and sensitivity analysis of uncertainties in computational mechanics, Int. J. Uncertain. Quan., 4, 387–421, https://doi.org/10.1615/int.j.uncertaintyquantification.2014006990, 2014. a

Arthern, R. J. and Williams, C. R.: The sensitivity of West Antarctica to the submarine melting feedback, Geophys. Res. Lett., 44, 2352–2359, https://doi.org/10.1002/2017gl072514, 2017. a

Aschwanden, A., Fahnestock, M. A., and Truffer, M.: Complex Greenland outlet glacier flow captured, Nat. Commun., 7, 10524, https://doi.org/10.1038/ncomms10524, 2016. a

Barletta, V. R., Bevis, M., Smith, B. E., Wilson, T., Brown, A., Bordoni, A., Willis, M., Khan, S. A., Rovira-Navarro, M., Dalziel, I., Smalley Jr., R., Kendrick, E., Konfal, S., Caccamise II, D. J., Aster, R. C., Nyblade, A., and Wiens, D. A.: Observed rapid bedrock uplift in Amundsen Sea Embayment promotes ice-sheet stability, Science, 360, 1335–1339, https://doi.org/10.1126/science.aao1447, 2018. a

Beckmann, A. and Goosse, H.: A parameterization of ice shelf–ocean interaction for climate models, Ocean Model., 5, 157–170, https://doi.org/10.1016/s1463-5003(02)00019-7, 2003. a

Bennett, M. R.: Ice streams as the arteries of an ice sheet: their mechanics, stability and significance, Earth-Sci. Rev., 61, 309–339, https://doi.org/10.1016/s0012-8252(02)00130-7, 2003. a

Bindschadler, R. A., Nowicki, S., Abe-Ouchi, A., Aschwanden, A., Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Khroulev, C., Levermann, A., Lipscomb, W. H., Martin, M. A., Morlighem, M., Parizek, B. R., Pollard, D., Price, S. F., Ren, D., Saito, F., Sato, T., Seddik, H., Seroussi, H., Takahashi, K., Walker, R., and Wang, W. L.: Ice-sheet model sensitivities to environmental forcing and their use in projecting future sea level (the SeaRISE project), J. Glaciol., 59, 195–224, https://doi.org/10.3189/2013jog12j125, 2013. a, b

Blatter, H.: Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients, J. Glaciol., 41, 333–344, https://doi.org/10.3189/s002214300001621x, 1995. a

Bolin, D. and Lindgren, F.: Excursion and contour uncertainty regions for latent Gaussian models, J. R. Stat. Soc. B, 77, 85–106, https://doi.org/10.1111/rssb.12055, 2015. a, b, c, d

Briggs, R., Pollard, D., and Tarasov, L.: A glacial systems model configured for large ensemble analysis of Antarctic deglaciation, The Cryosphere, 7, 1949–1970, https://doi.org/10.5194/tc-7-1949-2013, 2013. a, b

Brondex, J., Gagliardini, O., Gillet-Chaulet, F., and Durand, G.: Sensitivity of grounding line dynamics to the choice of the friction law, J. Glaciol., 63, 854–866, https://doi.org/10.1017/jog.2017.51, 2017. a, b, c

Brondex, J., Gillet-Chaulet, F., and Gagliardini, O.: Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law, The Cryosphere, 13, 177–195, https://doi.org/10.5194/tc-13-177-2019, 2019. a

Bueler, E. and Brown, J.: Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model, J. Geophys. Res., 114, F03008, https://doi.org/10.1029/2008jf001179, 2009. a

Calonne, N., Montagnat, M., Matzl, M., and Schneebeli, M.: The layered evolution of fabric and microstructure of snow at Point Barnola, Central East Antarctica, Earth Planet. Sc. Lett., 460, 293–301, https://doi.org/10.1016/j.epsl.2016.11.041, 2017. a

Chen, B., Haeger, C., Kaban, M. K., and Petrunin, A. G.: Variations of the effective elastic thickness reveal tectonic fragmentation of the Antarctic lithosphere, Tectonophysics, 746, 412–424, https://doi.org/10.1016/j.tecto.2017.06.012, 2018. a, b

Conrad, P. R., Davis, A. D., Marzouk, Y. M., Pillai, N. S., and Smith, A.: Parallel Local Approximation MCMC for Expensive Models, SIAM/ASA J. Uncertainty Quantification, 6, 339–373, https://doi.org/10.1137/16m1084080, 2018. a

Cornford, S. L., Martin, D. F., Payne, A. J., Ng, E. G., Le Brocq, A. M., Gladstone, R. M., Edwards, T. L., Shannon, S. R., Agosta, C., van den Broeke, M. R., Hellmer, H. H., Krinner, G., Ligtenberg, S. R. M., Timmermann, R., and Vaughan, D. G.: Century-scale simulations of the response of the West Antarctic Ice Sheet to a warming climate, The Cryosphere, 9, 1579–1600, https://doi.org/10.5194/tc-9-1579-2015, 2015. a, b

Cornford, S. L., Martin, D. F., Lee, V., Payne, A. J., and Ng, E. G.: Adaptive mesh refinement versus subgrid friction interpolation in simulations of Antarctic ice dynamics, Ann. Glaciol., 57, 1–9, https://doi.org/10.1017/aog.2016.13, 2016. a

Crestaux, T., Le Maître, O., and Martinez, J.-M.: Polynomial chaos expansion for sensitivity analysis, Reliab. Eng. Syst. Safe., 94, 1161–1172, https://doi.org/10.1016/j.ress.2008.10.008, 2009. a, b

de Boer, B., Dolan, A. M., Bernales, J., Gasson, E., Goelzer, H., Golledge, N. R., Sutter, J., Huybrechts, P., Lohmann, G., Rogozhina, I., Abe-Ouchi, A., Saito, F., and van de Wal, R. S. W.: Simulating the Antarctic ice sheet in the late-Pliocene warm period: PLISMIP-ANT, an ice-sheet model intercomparison project, The Cryosphere, 9, 881–903, https://doi.org/10.5194/tc-9-881-2015, 2015. a

DeConto, R. M. and Pollard, D.: Contribution of Antarctica to past and future sea-level rise, Nature, 531, 591–597, https://doi.org/10.1038/nature17145, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Depoorter, M. A., Bamber, J. L., Griggs, J. A., Lenaerts, J. T. M., Ligtenberg, S. R. M., van den Broeke, M. R., and Moholdt, G.: Calving fluxes and basal melt rates of Antarctic ice shelves, Nature, 502, 89–92, https://doi.org/10.1038/nature12567, 2013. a, b

Dinniman, M. S., Klinck, J. M., and Hofmann, E. E.: Sensitivity of Circumpolar Deep Water Transport and Ice Shelf Basal Melt along the West Antarctic Peninsula to Changes in the Winds, J. Climate, 25, 4799–4816, https://doi.org/10.1175/jcli-d-11-00307.1, 2012. a, b

Docquier, D., Perichon, L., and Pattyn, F.: Representing Grounding Line Dynamics in Numerical Ice Sheet Models: Recent Advances and Outlook, Surv. Geophys., 32, 417–435, https://doi.org/10.1007/s10712-011-9133-3, 2011. a

Drouet, A. S., Docquier, D., Durand, G., Hindmarsh, R., Pattyn, F., Gagliardini, O., and Zwinger, T.: Grounding line transient response in marine ice sheet models, The Cryosphere, 7, 395–406, https://doi.org/10.5194/tc-7-395-2013, 2013. a, b

Edwards, T. L., Brandon, M. A., Durand, G., Edwards, N. R., Golledge, N. R., Holden, P. B., Nias, I. J., Payne, A. J., Ritz, C., and Wernecke, A.: Revisiting Antarctic ice loss due to marine ice-cliff instability, Nature, 566, 58–64, https://doi.org/10.1038/s41586-019-0901-4, 2019. a, b

Fajraoui, N., Marelli, S., and Sudret, B.: Sequential Design of Experiment for Sparse Polynomial Chaos Expansions, SIAM/ASA J. Uncertainty Quantification, 5, 1061–1085, https://doi.org/10.1137/16m1103488, 2017. a

Favier, L., Durand, G., Cornford, S. L., Gudmundsson, G. H., Gagliardini, O., Gillet-Chaulet, F., Zwinger, T., Payne, A. J., and Le Brocq, A. M.: Retreat of Pine Island Glacier controlled by marine ice-sheet instability, Nat. Clim. Change, 4, 117–121, https://doi.org/10.1038/nclimate2094, 2014. a, b

French, J. P. and Hoeting, J. A.: Credible regions for exceedance sets of geostatistical data, Environmetrics, 27, 4–14, https://doi.org/10.1002/env.2371, 2015. a

Fretwell, P., Pritchard, H. D., Vaughan, D. G., Bamber, J. L., Barrand, N. E., Bell, R., Bianchi, C., Bingham, R. G., Blankenship, D. D., Casassa, G., Catania, G., Callens, D., Conway, H., Cook, A. J., Corr, H. F. J., Damaske, D., Damm, V., Ferraccioli, F., Forsberg, R., Fujita, S., Gim, Y., Gogineni, P., Griggs, J. A., Hindmarsh, R. C. A., Holmlund, P., Holt, J. W., Jacobel, R. W., Jenkins, A., Jokat, W., Jordan, T., King, E. C., Kohler, J., Krabill, W., Riger-Kusk, M., Langley, K. A., Leitchenkov, G., Leuschen, C., Luyendyk, B. P., Matsuoka, K., Mouginot, J., Nitsche, F. O., Nogi, Y., Nost, O. A., Popov, S. V., Rignot, E., Rippin, D. M., Rivera, A., Roberts, J., Ross, N., Siegert, M. J., Smith, A. M., Steinhage, D., Studinger, M., Sun, B., Tinto, B. K., Welch, B. C., Wilson, D., Young, D. A., Xiangbin, C., and Zirizzotti, A.: Bedmap2: improved ice bed, surface and thickness datasets for Antarctica, The Cryosphere, 7, 375–393, https://doi.org/10.5194/tc-7-375-2013, 2013. a, b, c, d

Fürst, J. J., Durand, G., Gillet-Chaulet, F., Tavard, L., Rankl, M., Braun, M., and Gagliardini, O.: The safety band of Antarctic ice shelves, Nat. Clim. Change, 6, 479–482, https://doi.org/10.1038/nclimate2912, 2016. a

Ghanem, R., Higdon, R., and Owhadi, H. (Eds.): Handbook of Uncertainty Quantification, Springer, https://doi.org/10.1007/978-3-319-12385-1, 2017. a, b, c, d

Ghanem, R. G. and Spanos, P. D.: Stochastic Finite Elements: A Spectral Approach, Dover Publications, Mineola, New York, revised edition, 2003. a

Gillet-Chaulet, F., Durand, G., Gagliardini, O., Mosbeux, C., Mouginot, J., Rémy, F., and Ritz, C.: Assimilation of surface velocities acquired between 1996 and 2010 to constrain the form of the basal friction law under Pine Island Glacier, Geophys. Res. Lett., 43, 10311–10321, https://doi.org/10.1002/2016gl069937, 2016. a, b, c

Gladstone, R., Schäfer, M., Zwinger, T., Gong, Y., Strozzi, T., Mottram, R., Boberg, F., and Moore, J. C.: Importance of basal processes in simulations of a surging Svalbard outlet glacier, The Cryosphere, 8, 1393–1405, https://doi.org/10.5194/tc-8-1393-2014, 2014. a

Gladstone, R. M., Warner, R. C., Galton-Fenzi, B. K., Gagliardini, O., Zwinger, T., and Greve, R.: Marine ice sheet model performance depends on basal sliding physics and sub-shelf melting, The Cryosphere, 11, 319–329, https://doi.org/10.5194/tc-11-319-2017, 2017. a

Goelzer, H., Nowicki, S., Edwards, T., Beckley, M., Abe-Ouchi, A., Aschwanden, A., Calov, R., Gagliardini, O., Gillet-Chaulet, F., Golledge, N. R., Gregory, J., Greve, R., Humbert, A., Huybrechts, P., Kennedy, J. H., Larour, E., Lipscomb, W. H., Le clec'h, S., Lee, V., Morlighem, M., Pattyn, F., Payne, A. J., Rodehacke, C., Rückamp, M., Saito, F., Schlegel, N., Seroussi, H., Shepherd, A., Sun, S., van de Wal, R., and Ziemen, F. A.: Design and results of the ice sheet model initialisation experiments initMIP-Greenland: an ISMIP6 intercomparison, The Cryosphere, 12, 1433–1460, https://doi.org/10.5194/tc-12-1433-2018, 2018. a

Golledge, N. R., Kowalewski, D. E., Naish, T. R., Levy, R. H., Fogwill, C. J., and Gasson, E. G. W.: The multi-millennial Antarctic commitment to future sea-level rise, Nature, 526, 421–425, https://doi.org/10.1038/nature15706, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p

Golledge, N. R., Levy, R. H., McKay, R. M., and Naish, T. R.: East Antarctic ice sheet most vulnerable to Weddell Sea warming, Geophys. Res. Lett., 44, 2343–2351, https://doi.org/10.1002/2016gl072422, 2017. a

Golledge, N. R., Keller, E. D., Gomez, N., Naughten, K. A., Bernales, J., Trusel, L. D., and Edwards, T. L.: Global environmental consequences of twenty-first-century ice-sheet melt, Nature, 566, 65–72, https://doi.org/10.1038/s41586-019-0889-9, 2019. a, b

Gomez, N., Mitrovica, J. X., Huybers, P., and Clark, P. U.: Sea level as a stabilizing factor for marine-ice-sheet grounding lines, Nat. Geosci., 3, 850–853, https://doi.org/10.1038/ngeo1012, 2010. a

Gomez, N., Pollard, D., and Mitrovica, J. X.: A 3-D coupled ice sheet – sea level model applied to Antarctica through the last 40 ky, Earth Planet. Sc. Lett., 384, 88–99, https://doi.org/10.1016/j.epsl.2013.09.042, 2013. a

Gopalan, G., Hrafnkelsson, B., Aðalgeirsdóttir, G., Jarosch, A. H., and Pálsson, F.: A Bayesian hierarchical model for glacial dynamics based on the shallow ice approximation and its evaluation using analytical solutions, The Cryosphere, 12, 2229–2248, https://doi.org/10.5194/tc-12-2229-2018, 2018. a

Greve, R. and Blatter, H.: Dynamics of Ice Sheets and Glaciers, Advances in Geophysical and Environmental Mechanics and Mathematics, Springer-Verlag Berlin Heidelberg, https://doi.org/10.1007/978-3-642-03415-2, 2009. a, b, c, d

Hadigol, M. and Doostan, A.: Least squares polynomial chaos expansion: A review of sampling strategies, Comput. Method. Appl. M., 332, 382–407, https://doi.org/10.1016/j.cma.2017.12.019, 2018. a, b

Hellmer, H. H., Kauker, F., Timmermann, R., Determann, J., and Rae, J.: Twenty-first-century warming of a large Antarctic ice-shelf cavity by a redirected coastal current, Nature, 485, 225–228, https://doi.org/10.1038/nature11064, 2012. a, b, c

Hellmer, H. H., Kauker, F., Timmermann, R., and Hattermann, T.: The Fate of the Southern Weddell Sea Continental Shelf in a Warming Climate, J. Climate, 30, 4337–4350, https://doi.org/10.1175/jcli-d-16-0420.1, 2017. a, b, c

Holland, P. R., Jenkins, A., and Holland, D. M.: The Response of Ice Shelf Basal Melting to Variations in Ocean Temperature, J. Climate, 21, 2558–2572, https://doi.org/10.1175/2007jcli1909.1, 2008. a

Hutter, K.: Theoretical Glaciology: Material Science of Ice and the Mechanics of Glaciers and Ice Sheets, D. Reidel Publishing Company, https://doi.org/10.1007/978-94-015-1167-4, 1983. a

Huybrechts, P. and de Wolde, J.: The Dynamic Response of the Greenland and Antarctic Ice Sheets to Multiple-Century Climatic Warming, J. Climate, 12, 2169–2188, https://doi.org/10.1175/1520-0442(1999)012<2169:tdrotg>2.0.co;2, 1999. a, b

Huybrechts, P., Abe-Ouchi, A., Marsiat, I., Pattyn, F., Payne, T., Ritz, C., and Rommelaere, V.: Report of the Third EISMINT Workshop on Model Intercomparison, European Science Foundation (Strasbourg), 1998. a, b

Isaac, T., Petra, N., Stadler, G., and Ghattas, O.: Scalable and efficient algorithms for the propagation of uncertainty from data through inference to prediction for large-scale problems, with application to flow of the Antarctic ice sheet, J. Comput. Phys., 296, 348–368, https://doi.org/10.1016/j.jcp.2015.04.047, 2015. a

Iskandarani, M., Wang, S., Srinivasan, A., Thacker, W. C., Winokur, J., and Knio, O. M.: An overview of uncertainty quantification techniques with application to oceanic and oil-spill simulations, J. Geophys. Res.-Oceans, 121, 2789–2808, https://doi.org/10.1002/2015jc011366, 2016. a

Jacobs, S. S., Jenkins, A., Giulivi, C. F., and Dutrieux, P.: Stronger ocean circulation and increased melting under Pine Island Glacier ice shelf, Nat. Geosci., 4, 519–523, https://doi.org/10.1038/ngeo1188, 2011. a

Janssens, I. and Huybrechts, P.: The treatment of meltwater retention in mass-balance parameterizations of the Greenland ice sheet, Ann. Glaciol., 31, 133–140, https://doi.org/10.3189/172756400781819941, 2000. a, b

Joughin, I., Tulaczyk, S., Bamber, J. L., Blankenship, D., Holt, J. W., Scambos, T., and Vaughan, D. G.: Basal conditions for Pine Island and Thwaites Glaciers, West Antarctica, determined using satellite and airborne data, J. Glaciol., 55, 245–257, https://doi.org/10.3189/002214309788608705, 2009. a

Joughin, I., Smith, B. E., and Holland, D. M.: Sensitivity of 21st century sea level to ocean-induced thinning of Pine Island Glacier, Antarctica, Geophys. Res. Lett., 37, L20502, https://doi.org/10.1029/2010gl044819, 2010. a

Joughin, I., Smith, B. E., and Medley, B.: Marine Ice Sheet Collapse Potentially Under Way for the Thwaites Glacier Basin, West Antarctica, Science, 344, 735–738, https://doi.org/10.1126/science.1249055, 2014. a, b

Knutti, R., Furrer, R., Tebaldi, C., Cermak, J., and Meehl, G. A.: Challenges in Combining Projections from Multiple Climate Models, J. Climate, 23, 2739–2758, https://doi.org/10.1175/2009jcli3361.1, 2010. a

Kopp, R. E., Horton, R. M., Little, C. M., Mitrovica, J. X., Oppenheimer, M., Rasmussen, D. J., Strauss, B. H., and Tebaldi, C.: Probabilistic 21st and 22nd century sea-level projections at a global network of tide-gauge sites, Earth's Future, 2, 383–406, https://doi.org/10.1002/2014ef000239, 2014. a

Larour, E., Schiermeier, J., Rignot, E., Seroussi, H., Morlighem, M., and Paden, J.: Sensitivity Analysis of Pine Island Glacier ice flow using ISSM and DAKOTA, J. Geophys. Res., 117, F02009, https://doi.org/10.1029/2011jf002146, 2012. a

Le Gratiet, L., Marelli, S., and Sudret, B.: Metamodel-Based Sensitivity Analysis: Polynomial Chaos Expansions and Gaussian Processes, in: Handbook of Uncertainty Quantification, edited by: Ghanem, R., Higdon, R., and Owhadi, H., chap. 38, 1289–1325, Springer, 2017. a

Le Maître, O. P. and Knio, O. M.: Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics, Springer Science & Business Media, https://doi.org/10.1007/978-90-481-3520-2, 2010. a, b, c, d

Lenaerts, J. T. M., Vizcaino, M., Fyke, J., van Kampenhout, L., and van den Broeke, M. R.: Present-day and future Antarctic ice sheet climate and surface mass balance in the Community Earth System Model, Clim. Dynam., 47, 1367–1381, https://doi.org/10.1007/s00382-015-2907-4, 2016. a

Ma, Y., Gagliardini, O., Ritz, C., Gillet-Chaulet, F., Durand, G., and Montagnat, M.: Enhancement factors for grounded ice and ice shelves inferred from an anisotropic ice-flow model, J. Glaciol., 56, 805–812, https://doi.org/10.3189/002214310794457209, 2010. a

MacAyeal, D. R.: Large-Scale Ice Flow over a Viscous Basal Sediment: Theory and Application to Ice Stream B, Antarctica, J. Geophys. Res., 94, 4071–4087, https://doi.org/10.1029/jb094ib04p04071, 1989. a

Maris, M. N. A., de Boer, B., Ligtenberg, S. R. M., Crucifix, M., van de Berg, W. J., and Oerlemans, J.: Modelling the evolution of the Antarctic ice sheet since the last interglacial, The Cryosphere, 8, 1347–1360, https://doi.org/10.5194/tc-8-1347-2014, 2014. a, b, c

Moholdt, G., Padman, L., and Fricker, H. A.: Basal mass budget of Ross and Filchner-Ronne ice shelves, Antarctica, derived from Lagrangian analysis of ICESat altimetry, J. Geophys. Res.-Earth, 119, 2361–2380, https://doi.org/10.1002/2014jf003171, 2014. a

Morland, L. W.: Unconfined Ice-Shelf Flow, in: Dynamics of the West Antarctic Ice Sheet, edited by: van der Veen, C. J. and Oerlemans, J., Kluwer Academic Publishers, 99–116, https://doi.org/10.1007/978-94-009-3745-1_6, 1987. a

Nowicki, S. M. J., Payne, A., Larour, E., Seroussi, H., Goelzer, H., Lipscomb, W., Gregory, J., Abe-Ouchi, A., and Shepherd, A.: Ice Sheet Model Intercomparison Project (ISMIP6) contribution to CMIP6, Geosci. Model Dev., 9, 4521–4545, https://doi.org/10.5194/gmd-9-4521-2016, 2016. a

Oakley, J. E. and O'Hagan, A.: Probabilistic sensitivity analysis of complex models: a Bayesian approach, J. R. Stat. Soc. B, 66, 751–769, https://doi.org/10.1111/j.1467-9868.2004.05304.x, 2004. a

Pattyn, F.: A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes, J. Geophys. Res., 108, 2382, https://doi.org/10.1029/2002jb002329, 2003. a

Pattyn, F.: Sea-level response to melting of Antarctic ice shelves on multi-centennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0), The Cryosphere, 11, 1851–1878, https://doi.org/10.5194/tc-11-1851-2017, 2017. a, b, c, d, e, f, g, h, i, j, k

Pattyn, F. and Durand, G.: Why marine ice sheet model predictions may diverge in estimating future sea level rise, Geophys. Res. Lett., 40, 4316–4320, https://doi.org/10.1002/grl.50824, 2013. a, b

Pattyn, F., Schoof, C., Perichon, L., Hindmarsh, R. C. A., Bueler, E., de Fleurian, B., Durand, G., Gagliardini, O., Gladstone, R., Goldberg, D., Gudmundsson, G. H., Huybrechts, P., Lee, V., Nick, F. M., Payne, A. J., Pollard, D., Rybak, O., Saito, F., and Vieli, A.: Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP, The Cryosphere, 6, 573–588, https://doi.org/10.5194/tc-6-573-2012, 2012. a, b

Pattyn, F., Perichon, L., Durand, G., Favier, L., Gagliardini, O., Hindmarsh, R. C. A., Zwinger, T., Albrecht, T., Cornford, S., Docquier, D., Fürst, J. J., Goldberg, D., Gudmundsson, G. H., Humbert, A., Hütten, M., Huybrechts, P., Jouvet, G., Kleiner, T., Larour, E., Martin, D., Morlighem, M., Payne, A. J., Pollard, D., Rückamp, M., Rybak, O., Seroussi, H., Thoma, M., and Wilkens, N.: Grounding-line migration in plan-view marine ice-sheet models: results of the ice2sea MISMIP3d intercomparison, J. Glaciol., 59, 410–422, https://doi.org/10.3189/2013jog12j129, 2013. a, b, c, d, e

Pattyn, F., Favier, L., Sun, S., and Durand, G.: Progress in Numerical Modeling of Antarctic Ice-Sheet Dynamics, Curr. Clim. Change Rep., 3, 174–184, https://doi.org/10.1007/s40641-017-0069-7, 2017. a

Pattyn, F., Ritz, C., Hanna, E., Asay-Davis, X., DeConto, R., Durand, G.,
Favier, L., Fettweis, X., Goelzer, H., Golledge, N. R., Munneke, P. K.,
Lenaerts, J. T. M., Nowicki, S., Payne, A. J., Robinson, A., Seroussi, H.,
Trusel, L. D., and van den Broeke, M.: The Greenland and Antarctic ice
sheets under 1.5 ^{∘}C global warming, Nat. Clim. Change, 8,
1053–1061, https://doi.org/10.1038/s41558-018-0305-8, 2018. a, b

Petra, N., Martin, J., Stadler, G., and Ghattas, O.: A Computational Framework for Infinite-Dimensional Bayesian Inverse Problems, Part II: Stochastic Newton MCMC with Application to Ice Sheet Flow Inverse Problems, SIAM J. Sci. Comput., 36, A1525–A1555, https://doi.org/10.1137/130934805, 2014. a

Pollard, D. and DeConto, R. M.: Modelling West Antarctic ice sheet growth and collapse through the past five million years, Nature, 458, 329–332, https://doi.org/10.1038/nature07809, 2009. a

Pollard, D. and DeConto, R. M.: Description of a hybrid ice sheet-shelf model, and application to Antarctica, Geosci. Model Dev., 5, 1273–1295, https://doi.org/10.5194/gmd-5-1273-2012, 2012a. a, b, c, d, e, f, g, h, i

Pollard, D. and DeConto, R. M.: A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica, The Cryosphere, 6, 953–971, https://doi.org/10.5194/tc-6-953-2012, 2012b. a

Pollard, D., DeConto, R. M., and Alley, R. B.: Potential Antarctic Ice Sheet retreat driven by hydrofracturing and ice cliff failure, Earth Planet. Sc. Lett., 412, 112–121, https://doi.org/10.1016/j.epsl.2014.12.035, 2015. a, b

Pollard, D., Chang, W., Haran, M., Applegate, P., and DeConto, R.: Large ensemble modeling of the last deglacial retreat of the West Antarctic Ice Sheet: comparison of simple and advanced statistical techniques, Geosci. Model Dev., 9, 1697–1723, https://doi.org/10.5194/gmd-9-1697-2016, 2016. a, b, c

Pritchard, H. D., Ligtenberg, S. R. M., Fricker, H. A., Vaughan, D. G., van den Broeke, M. R., and Padman, L.: Antarctic ice-sheet loss driven by basal melting of ice shelves, Nature, 484, 502–505, https://doi.org/10.1038/nature10968, 2012. a

Rasmussen, C. E. and Williams, C. K. I.: Gaussian Processes for Machine Learning, Adaptive Computation and Machine Learning, The MIT Press, Cambridge, Massachussetts, 2006. a

Reese, R., Albrecht, T., Mengel, M., Asay-Davis, X., and Winkelmann, R.: Antarctic sub-shelf melt rates via PICO, The Cryosphere, 12, 1969–1985, https://doi.org/10.5194/tc-12-1969-2018, 2018a. a, b, c, d

Reese, R., Winkelmann, R., and Gudmundsson, G. H.: Grounding-line flux formula applied as a flux condition in numerical simulations fails for buttressed Antarctic ice streams, The Cryosphere, 12, 3229–3242, https://doi.org/10.5194/tc-12-3229-2018, 2018b. a

Rignot, E., Jacobs, S., Mouginot, J., and Scheuchl, B.: Ice-Shelf Melting Around Antarctica, Science, 341, 266–270, https://doi.org/10.1126/science.1235798, 2013. a, b

Rignot, E., Mouginot, J., Morlighem, M., Seroussi, H., and 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, https://doi.org/10.1002/2014gl060140, 2014. a, b

Ritz, C., Edwards, T. L., Durand, G., Payne, A. J., Peyaud, V., and Hindmarsh, R. C. A.: Potential sea-level rise from Antarctic ice-sheet instability constrained by observations, Nature, 528, 115–118, https://doi.org/10.1038/nature16147, 2015. a, b, c, d, e, f, g, h, i, j, k, l

Robert, C. P. and Casella, G.: Monte Carlo Statistical Methods, Springer Texts in Statistics, Springer Science & Business Media, 2nd Edn., https://doi.org/10.1007/978-1-4757-4145-2, 2013. a

Rommelaere, V. and Ritz, C.: A thermomechanical model of ice-shelf flow, Ann. Glaciol., 23, 13–20, https://doi.org/10.3189/s0260305500013203, 1996. a

Ruckert, K. L., Shaffer, G., Pollard, D., Guan, Y., Wong, T. E., Forest, C. E., and Keller, K.: Assessing the Impact of Retreat Mechanisms in a Simple Antarctic Ice Sheet Model Using Bayesian Calibration, PLoS ONE, 12, e0170052, https://doi.org/10.1371/journal.pone.0170052, 2017. a, b

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., and Tarantola, S.: Global Sensitivity Analysis: The Primer, John Wiley & Sons, West Sussex, United Kingdom, 2008. a, b, c

Schlegel, N.-J., Seroussi, H., Schodlok, M. P., Larour, E. Y., Boening, C., Limonadi, D., Watkins, M. M., Morlighem, M., and van den Broeke, M. R.: Exploration of Antarctic Ice Sheet 100-year contribution to sea level rise and associated model uncertainties using the ISSM framework, The Cryosphere, 12, 3511–3534, https://doi.org/10.5194/tc-12-3511-2018, 2018. a, b, c, d, e, f, g, h

Schmidtko, S., Heywood, K. J., Thompson, A. F., and Aoki, S.: Multidecadal warming of Antarctic waters, Science, 346, 1227–1231, https://doi.org/10.1126/science.1256117, 2014. a, b, c

Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability and hysteresis, J. Geophys. Res., 112, F03S28, https://doi.org/10.1029/2006JF000664, 2007a. a, b

Schoof, C.: Marine ice-sheet dynamics. Part 1. The case of rapid sliding, J. Fluid. Mech., 573, 27–55, https://doi.org/10.1017/s0022112006003570, 2007b. a, b, c

Schoof, C., Davis, A. D., and Popa, T. V.: Boundary layer models for calving marine outlet glaciers, The Cryosphere, 11, 2283–2303, https://doi.org/10.5194/tc-11-2283-2017, 201 a

Schäfer, M., Zwinger, T., Christoffersen, P., Gillet-Chaulet, F., Laakso, K., Pettersson, R., Pohjola, V. A., Strozzi, T., and Moore, J. C.: Sensitivity of basal conditions in an inverse model: Vestfonna ice cap, Nordaustlandet/Svalbard, The Cryosphere, 6, 771–783, https://doi.org/10.5194/tc-6-771-2012, 2012. a

Scott, D. W.: Multivariate Density Estimation: : Theory, Practice, and Visualization, John Wiley & Sons, 2nd Edn., https://doi.org/10.1002/9781118575574, 2015. a

Shepherd, A., Ivins, E., Rignot, E., Smith, B., van den Broeke, M., Velicogna, I., Whitehouse, P., Briggs, K., Joughin, I., Krinner, G., Nowicki, S., Payne, T., Scambos, T., Schlegel, N. A. G., Agosta, C., Ahlstrøm, A., Babonis, G., Barletta, V., Blazquez, A., Bonin, J., Csatho, B., Cullather, R., Felikson, D., Fettweis, X., Forsberg, R., Gallee, H., Gardner, A., Gilbert, L., Groh, A., Gunter, B., Hanna, E., Harig, C., Helm, V., Horvath, A., Horwath, M., Khan, S., Kjeldsen, K. K., Konrad, H., Langen, P., Lecavalier, B., Loomis, B., Luthcke, S., McMillan, M., Melini, D., Mernild, S., Mohajerani, Y., Moore, P., Mouginot, J., Moyano, G., Muir, A., Nagler, T., Nield, G., Nilsson, J., Noel, B., Otosaka, I., Pattle, M. E., Peltier, W. R., Pie, N., Rietbroek, R., Rott, H., Sandberg-Sørensen, L., Sasgen, I., Save, H., Scheuchl, B., Schrama, E., Schröder, L., Seo, K.-W., Simonsen, S., Slater, T., Spada, G., Sutterley, T., Talpe, M., Tarasov, L., van de Berg, W. J., van der Wal, W., van Wessem, M., Vishwakarma, B. D., Wiese, D., and Wouters, B.: Mass balance of the Antarctic Ice Sheet from 1992 to 2017, Nature, 558, 219–222, https://doi.org/10.1038/s41586-018-0179-y, 2018. a

Stein, M.: Large Sample Properties of Simulations Using Latin Hypercube Sampling, Technometrics, 29, 143–151, https://doi.org/10.2307/1269769, 1987. a, b

Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boshung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M. (Eds.): Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, United Kingdom, https://doi.org/10.1017/cbo9781107415324, 2013. a

Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An Overview of CMIP5 and the Experiment Design, B. Am. Meteorol. Soc., 93, 485–498, https://doi.org/10.1175/bams-d-11-00094.1, 2012. a, b

Timmermann, R. and Hellmer, H. H.: Southern Ocean warming and increased ice shelf basal melting in the twenty-first and twenty-second centuries based on coupled ice-ocean finite-element modelling, Ocean Dynam., 63, 1011–1026, https://doi.org/10.1007/s10236-013-0642-0, 2013. a, b

Timmermann, R., Wang, Q., and Hellmer, H.: Ice-shelf basal melting in a global finite-element sea-ice/ice-shelf/ocean model, Ann. Glaciol., 53, 303–314, https://doi.org/10.3189/2012aog60a156, 2012. a

Tsai, V. C., Stewart, A. L., and Thompson, A. F.: Marine ice-sheet profiles and stability under Coulomb basal conditions, J. Glaciol., 61, 205–215, https://doi.org/10.3189/2015jog14j221, 2015. a, b, c, d

Van der Wal, W., Whitehouse, P. L., and Schrama, E. J. O.: Effect of GIA models with 3D composite mantle viscosity on GRACE mass balance estimates for Antarctica, Earth Planet. Sc. Lett., 414, 134–143, https://doi.org/10.1016/j.epsl.2015.01.001, 2015. a, b

Van Wessem, J., Reijmer, C. H., Morlighem, M., Mouginot, J., Rignot, E., Medley, B., Joughin, I., Wouters, B., Depoorter, M. A., Bamber, J. L., Lenaerts, J. T. M., van de Berg, W. J., van den Broeke, M. R., and van Meijgaard, E.: Improved representation of East Antarctic surface mass balance in a regional atmospheric climate model, J. Glaciol., 60, 761–770, https://doi.org/10.3189/2014jog14j051, 2014. a

Vaughan, D. G., Comiso, J. C., Allison, I., Carrasco, J., Kaser, G., Kwok, R., Mote, P., Murray, T., Paul, F., Ren, J., Rignot, E., Solomina, O., Steffen, K., and Zhang, T.: Observations: Cryosphere, in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xi, Y., Bex, V., and Midgley, P. M., Cambridge University Press, 317–382, 2013. a

Waibel, M. S., Hulbe, C. L., Jackson, C. S., and Martin, D. F.: Rate of Mass Loss Across the Instability Threshold for Thwaites Glacier Determines Rate of Mass Loss for Entire Basin, Geophys. Res. Lett., 45, 809–816, https://doi.org/10.1002/2017gl076470, 2018. a, b, c

Weertman, J.: On the Sliding of Glaciers, J. Glaciol., 3, 33–38, https://doi.org/10.3189/s0022143000024709, 1957. a

Weis, M., Greve, R., and Hutter, K.: Theory of shallow ice shelves, Continuum Mech. Therm., 11, 15–50, https://doi.org/10.1007/s001610050102, 1999. a

Yu, H., Rignot, E., Seroussi, H., and Morlighem, M.: Retreat of Thwaites Glacier, West Antarctica, over the next 100 years using various ice flow models, ice shelf melt scenarios and basal friction laws, The Cryosphere, 12, 3861–3876, https://doi.org/10.5194/tc-12-3861-2018, 2018. a

- Abstract
- Introduction
- Model description and methods
- Results
- Discussion
- Conclusions
- Data availability
- Appendix A: Polynomial chaos expansion
- Appendix B: Sobol sensitivity indices
- Appendix C: Confidence regions for grounded ice
- Author contributions
- Competing interests
- Acknowledgements
- Review statement
- References
- Supplement

- Abstract
- Introduction
- Model description and methods
- Results
- Discussion
- Conclusions
- Data availability
- Appendix A: Polynomial chaos expansion
- Appendix B: Sobol sensitivity indices
- Appendix C: Confidence regions for grounded ice
- Author contributions
- Competing interests
- Acknowledgements
- Review statement
- References
- Supplement