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

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 subshelf 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 sealevel 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.

tively 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 the 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 20 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 theory :::::: theories and new methods that are of interest to be applied to uncertainty quantification of ice-sheet models, beyond the theory and methods that , ,  and  have 25 already applied.
We perform simulations of the response of the Antarctic ice sheet (Fig. 1a) to environmental and parametric perturbations over the period 2000-3000 . Input data include the present-day ice-sheet geometry and bedrock topography from the Bedmap2 15 dataset (Fretwell et al., 2013), the geothermal heat flux by  and the spatially varying effective lithosphere thickness by ? to compute glacial isostatic adjustment more realistically with an elastic lithosphere/relaxed asthenosphere model. We perform simulations at a spatial resolution of 20 km, which allows to keep computational cost tractable. We found that the results with the f. ETISh model are essentially scale independent, as suggested by .
The scenario plays a significant role in the amplitude and speed of the AIS retreat. Recent studies  DeConto and  have shown no substantial retreat of the grounding-line position in the strongly mitigated RCP 2.6 10 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  have also highlighted that the hydrofacturing and ice cliff ::::::::::::: hydrofracturing ::: and ::::::: ice-cliff failure mechanisms (not included in the f.ETISh model 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 15 basins.

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 . Iceberg calving can have strong feedback effects as it affects ice-shelf buttressing ::::::::::::::::  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 :::::: ice-cliff : instability 30 (Pollard et al., 2015).
We introduce a 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 15 an uncertainty range from 0.2 to 1, where a value of 0.5 means that the ice in the ice shelves is two times more viscous than without shelf tuning.

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(Gomez et al., , 2013Adhikari et al., 2014). The amplitude of the glacial 20 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; ?) :::::::::::::::::::::::::::::::::::: (Van der 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. 25 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 τ e :: τ w has an uncertainty range from 1000 to 3500 years.

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 the UQ community ::: UQ : to estimate these statistical descriptors in a nonintrusive manner 20 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 25 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 expansion ::::::::: expansions : require the parameters-to-projection relationship to be sufficiently smooth (no discontinuity or highly nonlinear behaviour) to allow an efficient approximation as a low-degree 30 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 ::::: Monte ::::: Carlo : method (Robert and Casella, 2013), in which samples of the parameters are drawn ran-domly 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 5 that accelerates UQanalysis; (ii) it provides an explicit view on the relationship between the parameters and the projection, highlighting potential linear or nonlinear dependences and interactions between the parameters; (iii) it allows to interpolate efficiently 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).

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 indicates that the uncertain parameter has no impact on the projection uncertainty.

Confidence regions for grounded-ice retreat
To gain insight into the impact of the uncertainty in determining which subregions :::::: regions of Antarctica are most at risk of 5 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 subregion ::::: 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 10 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  for Gaussian random fields to our glaciological context.

15
We present nominal and probabilistic projections (relative to 2000 CE) for short-term (2100), medium-term (2300) and longterm (3000) time scales under different RCP scenarios and sliding laws. All results to follow have been obtained with the SGL parameterisation for the grounding-line migration, except in a discussion at the end of the section, where we discuss the impact of using the TGL parameterisation. 20 We first present nominal projections obtained using the nominal values of the parameters given in Table1 ::: 1. We first present results under nominal conditions in order to assess subsequently the impact of uncertainties on AIS sea-level rise projections.

Nominal projections
Under nominal conditions, we find (Table 2) in RCP 2.6 a nominal :: an AIS contribution to sea-level rise of 2 cm :: sea :::: level ::: of :::: 0.02 :: m by 2100, 7 cm :::: 0.07 :: m by 2300 and 20 cm :::: 0.20 :: m : by 3000 and in RCP 8.5 a nominal :: an : AIS contribution to sea-level rise of 5 cm :: sea ::::: level :: of ::::: 0.05 :: m : by 2100, 59 cm :::: 0.59 :: m : by 2300 and 3.9 :::: 3.90 m by 3000. In Fig. 3a, we represented the 25 nominal projections as a function of time in all RCP scenarios. We find that the nominal AIS contribution to sea-level rise ::: 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.

Parameters-to-projections relationship
One of the advantages of a polynomial chaos expansion is that it provides an explicit approximation to the parameters-toprojections relationship, which can be visualised to gain insight into the relationship between the parameters and the projections.
In Fig. 4 and 5, we used the PC expansions to visualise how the projections depend on each parameter individually (one-at-   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 rather quadratic dependence can be explained by the influence of E shelf on two competing processes: on the one hand, a higher value of E shelf softens the ice, thus leading to faster ice flow in the ice shelves; on the other hand, :: but : a higher value of E shelf increases buttressing at the grounding line, thus making ::: also ::::: leads :: to ::::::: ice-shelf :::::::: thinning, ::: thus :::::::: reducing ice flux at the grounding lineless important. In addition, we find that ∆GMSL depends only little on the bedrock relaxation times . 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, recent studies (Barletta et al., 2018) have : a ::::: recent ::::: study :: by :::::::::::::::::::: Barletta et al. (2018) has : suggested a bedrock relaxation time scale of the order of We find in RCP 8.5 similar trends. However, whereas F calv , F melt and E shelf influence ∆GMSL equally significantly in RCP 2.6, the melt factor influences ∆GMSL most significantly in RCP 8.5. Whereas Fig. 5d-e show ::::: shows : 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 sub-shelf melting ::::: F melt for long-term projections in RCP 8.5.
∆GMSL projections are the median projections with their 5-95 % probability intervals between parentheses. Figure 9 provides the Sobol sensitivity indices for the change in GMSL on short-term, medium-term and long-term time scales 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 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 :: ice :::: flow : and buttressing in ice shelves, as well as calving, which reduces the extent of ice shelves and their buttressing.

Stochastic sensitivity analysis
By contrast, in warmer RCP scenarios and for longer time scales, 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 the sea-level rise projections. As Finally, we find that, in all RCP scenarios and under all sliding laws, the uncertainty in the bedrock relaxation times for 15 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.
18 3.5 Projections of grounded-ice retreat Figure 10 provides insight into the subregions ::::: regions : of Antarctica that are most at risk of ungrounding in different RCP scenarios and at different time scales under the weakly nonlinear sliding law. Figure 10 was obtained as follows. First, we determined the 100 % confidence region for grounded ice, that is, the subregion ::::: 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 retreats to within the grey 5 region. Then, we determined the 95 % confidence region for grounded ice, that is, the subregion ::::: 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 retreats 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. 10 We find that ice is certain to remain :::::: 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 ::::::::: fast-flowing : ice streams, especially in Siple Coast, in the Ronne basin notably the Ellsworth land and in the Amundsen Sea sector. In warm RCP scenarios and at longer time scales, 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 15 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 (differences between between both regions are 20 only of few tenths of kilometresby 3000 :: the :::: 100 :: % ::::::::: confidence ::::: region ::: for :::::::: grounded ::: ice :: by ::::: 3000 :::: only ::::: differ :::: from ::: the :::::::::: present-day :::::::: grounded :: ice :::::: region ::: by : a :::: few :::: tens :: of ::::::::: kilometres). In RCP 4.5, ice is certain to remain :::::: 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 time scales. 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 25 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 Antarctica :::::::: Antarctic 30 ice sheet is underway by 3000 in RCP 8.5, while in Fig. 10l the West Antarctica :::::::: Antarctic ::: ice :::: sheet : has already collapsed by 3000.
Additionally, we compared the projections under the weakly nonlinear sliding law with projections under the other sliding laws, that is the viscous sliding law ( Fig. S6 :: S7) and the strongly nonlinear sliding law ( Fig. S7 :: 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 nonlinear 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. 5 However, Fig. 10i and Fig. S7i :: S8i : suggest that ungrounding may be less significant in Siple Coast under the strongly nonlinear sliding law than the weakly nonlinear 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 10 in the former region to nonlinearity and a more viscous ::::: plastic response.

Comparison of the sea-level rise projections with previous work
Regarding the short-term AIS contribution to sea-level rise :: sea ::::: level, we projected in the RCP 2.6 scenario a median of 2-3 cm parameterisations, our results suggested that the AIS contribution to sea-level rise ::: sea :::: level : does not exceed 50 cm :: 0.5 ::: m by 2100 with a probability of at least 99 %. These results are in agreement with , who determined an AIS contribution to sea-level rise :: sea ::::: level that lies between 5 cm and 30 cm :::: 0.05 :: m :::: and :::: 0.30 :: m, and , who found an AIS contribution to sea-level rise :: sea :::: level : that reaches at most 38 cm ::: 0.38 :: m : in RCP 8.5 with sub-grid interpolation of basal melting at the grounding line. Yet, projections by DeConto and  and , who applied 5 the more sensitive Budd-type friction law , can exceed 50 cm :: 0.5 :: m : 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 rise ::: sea :::: level, our projections under the SGL parameterisation are similar to other estimates by . In particular, both studies suggest that the AIS contribution to sea-level rise :: 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,

Projections of ice loss and grounding-line retreat under parametric uncertainty
The significance of the contribution of the Antarctic ice sheet to sea-level rise ::: sea :::: level : under climate change is primarily 10 controlled by the contribution of ::::::::: 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 rise ::: sea :::: level by 3000 is rather limited (less than 1 metre) in RCP 2.6, RCP 4.5 and RCP 6.0, while an 15 accelerated ice loss that leads to a contribution of several meters ::::: metres : is triggered in RCP 8.5 (Fig. 3a)and that : . :: 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 20 projections hold irrespectively of parametric uncertainty: the AIS contribution to sea-level rise ::: sea :::: level : by 3000 has a 95 % quantile of up to 91 cm ::: 0.91 :: m : (Table 3) and there is a limited risk that the grounding line retreats beyond our nominal projections ( Fig. 3b and Fig. 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 rise :: 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 25 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 rise ::: sea :::: level : by 3000 has an uncertainty range with a 5 % quantile above 82 cm ::: 0.82 ::: m and a 95 % quantile of up to 7.68 m (Table 3) and there is a high risk of triggering a 30 complete disintegration of the West Antarctic ice sheet (Fig. 10l).
In conclusion, the projections hold irrespectively of parametric uncertainty in the strongly mitigated RCP 2.6 scenario: accommodating parametric uncertainty in the ice-sheet model does not question :::: 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 sensi-tive to parametric uncertainty for intermediate scenarios such as RCP 4.5 and RCP 6.0: accommodating parametric uncertainty in the ice-sheet model questions :::: 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 irrespectively of parametric uncertainty.
Data availability. 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 . 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 K. Bulthuis (kevin.bulthuis@uliege.be).
The MATLAB code used to analyse the results is also available on request from the same author.

Appendix A: Polynomial chaos expansion
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 x = (x 1 , . . . , 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 which ::: 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. 10 In order to fit the PC expansion in Eq. (A1) to the ice-sheet model, we calculate the coefficients using a (weighted) leastsquares approach: where x (i) , 1 i N is a set of N training points in the parameter space, y (i) = g(x (i) ), 1 i N is the set of model responses at the training points, c = (c 0 , c 1 , . . . , c K ) collects the PC coefficients and w (i) , 1 i N is the set of weights. 15 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 sampling :::::: 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 correspond ::::::::: corresponds : to K = 56 for d = 5, as we strive to reproduce the overall 20 trend of the response function without overfitting small local variations.
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 10 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 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.

Appendix B: Sobol sensitivity indices
In this appendix, we concisely provide further details about how we used stochastic sensitivity analysis; we refer the reader 25 to, for instance, Ghanem et al. (2017) and Saltelli et al. (2008) for more comprehensive treatments of the theory and various applications.
Assuming ::: 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 28 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 constantand : , the main effects are defined by ::: and ::: the ::::::::: interaction ::::: effect ::: are ::::: given :: by : :::::::::::::: . , 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 10 among all functions that depend only on x i . As a consequence of this least-squares best approximation property, the functions g 0 , g 1 , . . . , g d , g I are orthonormal with respect to the probability density function of the uncertain parameter :::::::: parameters.
As a consequence of the orthonormality of the functions g 0 , g 1 , . . . , g d , g I , the variance σ 2 of the model response can be decomposed as in which :::: with ::::::::::::::::::::::::::::::: (B6) where p i denotes the probability density function for :: 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 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 5 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 .

Appendix C: Confidence regions for grounded ice
In this appendix, we concisely provide further details about how we defined and computed confidence regions for grounded ice.

10
To distinguish between grounded ice and floating ice at a given time, the f.ETISh model  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 subregion ::::: region : of Antarctica that has a probability of at least α to be included in the grounded ice domain, that is, We compute the confidence regions using an adaptation of a thresholding algorithm by . We 20 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  consider Gaussian random fields, here we work with non-Gaussian random fields, which requires to evaluate the marginal probability density functions and the probability of inclusion with Monte-Carlo ::::: Monte ::::: Carlo simulation.          bar ::: and ::: the ::: unit ::::: value ::::::: represents ::: the :::::::: interaction ::::: index.
Probability to remain grounded  Response to the Interactive comment on "Uncertainty quantification of the multicentennial response of the Antarctic Ice Sheet to climate change" by Kevin Bulthuis et al.
We would like to thank anonymous referee #1 for the time dedicated to this manuscript and his/her constructive comments to improve the general quality and readibility of the manuscript. We will try to give a proper response to his/her comments. The manuscript has been revised to include more references, background material and discussions about the description and the limitations of the model and the methodology. For each referee's comment (written in blue), we included below a response (written in black) and proposed means to improve the manuscript.

Summary statement
The manuscript "Uncertainty quantification of the multi-centennial response of the Antarctic Ice Sheet to climate change" by K. Bulthuis et al. provides an assessment of the response of the Antarctic ice sheet and the associated uncertainty to several climate change scenarios over the next few centuries, by combining numerical ice sheet model simulations and probabilistic methods. Using new probabilistic methods, and especially emulators, starts to be used in glaciology and provides an interesting alternative to costly ice flow models. The results show the relative impact of uncertainties in different ice flow model parameters as well as external parameters. They also suggest that the marine ice sheet instability triggers large contribution to sea level for some combinations of parameters for the RCP 4.5 and 6.0 scenarios, while the instability is almost never triggered for the RCP 2.0 scenario and always triggered for RCP 8.5 scenario.
The paper is usually well written and the figures appropriate, but some references and background material are sometimes missing, especially in the introduction and the model references. I am also surprised by how the problems of the relatively low model resolution and simple parameterization of the grounding line are treated and not discussed in more details.
Though it makes sense to use such configurations for large ensemble simulations of the entire Antarctic ice sheet, it seems a bit surprising that using such a configuration is needed when an emulator is used; I though that the emulator was allowing to run fewer but more accurate simulations used only for the calibration and then run a large number of cheap emulated runs to analyze the uncertainties in detail. This should at least be better acknowledged, especially in the limitations.
We agree with this general comment and we have made changes to the manuscript to account for this comment. See answers below to major comments and specific comments to questions about p.4 l.11-12, p.4 l.24 and p.9 l.32.

Major comments
The manuscript suggests that using a 20 km model resolution and a flux condition at the grounding line are reasonable assumptions and that the impact on the results is limited. I think this is a bit misleading as previous intercomparison experiments  have rather different conclusions. So this should be at least better explained and acknowledged in the limitations. Furthermore, I thought that the goal of using emulators was to limit the number of runs needed, and therefore that it would allow using more computationally expensive runs to calibrate the emulator, which does not seem to be the case here.
The referee is right in pointing out that the model resolution and the flux condition should be discussed in more detail given that they represent two major limitations of the simulations. In our opinion, using a 20-km model resolution and a flux condition at the grounding line are acceptable assumptions in large-scale and long-term ice-sheet simulations and large-ensemble simulations (see discussion below).
We are aware of the limitations of the flux condition (short transients, buttressing) but we expect the impact of this parameterisation to be limited when compared to the impact of uncertainties in model parameters and future scenarios. Furthermore, Pattyn et al. (2012 have shown that the use of a flux condition at the grounding line allows to reproduce qualitatively the marine ice-sheet instability (MISI) and simulations of the Antarctic ice sheet using a flux condition (Pollard and DeConto, 2012; have also been able to reproduce qualitatively MISI in large-scale ice-sheet simulations. In addition, it is true that a higher resolution would allow to capture properly certain mechanisms that control grounding-line migration such as bedrock irregularities and ice-shelf pinning points as well as important small glaciers such as Pine Island and Thwaites glaciers, which represent only a few grid points with a 20-km resolution. Yet,  has shown that using the f.ETISh model with higher spatial resolutions gives similar results for given perturbations. We have modified Sect. 2.1 to discuss in more detail the applicability and limitations of using a flux condition at the grounding line (see also answer about p.4 l.11-12 in the specific comments).
The referee is also right in pointing out that using an emulator allows to limit the number of computationally expensive runs. Yet, in this manuscript, we are interested in performing an uncertainty quantification analysis of the model under different sources of uncertainty. A typical uncertainty quantification analysis with a Monte Carlo approach (or related approach such as Latin hypercube sampling) requires in general tenths/hundreds of thousands to millions of runs to achieve proper estimations for the statistical descriptors (mean, variance, . . . ) of the output. This is intractable even with a 20-km model resolution (with a computational cost of ∼ 8h per forward simulation). In addition, we tried to investigate a sufficiently broad ensemble of uncertainties with an ensemble of five parameters and an ensemble of 20 distinct model configurations (4 RCP scenarios with 4 sliding laws and Tsai's flux condition). Each of these 20 distinct model configurations requires the construction of a different emulator because these model configurations represent discrete/categorical variables. For each of these configurations, we ran a set of 500 forward simulations (training set) to build the PC emulator, that is, we had to perform a total of 10000 forward simulations for the 20 model configurations. This is still a large total computational time (given our computational ressources) even with a 20-km model resolution. The emulator does not replace the forward ice-sheet numerical model, which is necessary to capture all possible nonlinear effects of the ice-sheet response. The choice of a 20-km model resolution is a tradeoff between the total computational time and the number of distinct configurations and uncertain parameters investigated in this paper. We clarify the choice for this resolution in the manuscript and its limitations (notably the inability to represent accurately important ice streams such as Pine Island and Thwaites glaciers).
Similarly, basal melting under ice shelves is known to have a very large impact on ice dynamics over both short and long timescales. What is done for future scenarios regarding this melt is however not clear and should be explained in more details. Also, the ocean warming follows a simple parameterisation that mostly depends on atmospheric warming. We know that the future evolution of the Southern Ocean remains very uncertain and that changes in ocean cir-culation rather than ocean warming are expected to cause the largest changes, so that simple parameterisations are unlikely to accurately represent these changes. This should be more clearly acknowledged and it may be good to provide some perspectives on the uncertainties associated to the future ice shelf melt.
This suggestion is totally relevant given the very large impact of sub-shelf melting on ice dynamics and there is a clear need for more discussion about this process in the manuscript. The proper way to include ocean warming would actually require the coupling of the ice-sheet model with an ocean model in order to account for oceanic processes such as changes in ocean currents. Still, the representation of sub-shelf melting heavily depends on how ocean models represent these processes.
We have developed the discussion about sub-shelf melting and the limitations of its representation in Sect 2. We have included a new description of the link between atmospheric and oceanic warmings in Sect. 2.1 and we have also included a more detailed discussion about the physics of sub-shelf melting and how it is implemented and linked with atmospheric forcing in the f.ETISH model in Sect. 2.2.5 and given more details about the limitations in Sect. 4.5 (see also the answers to the specific comments about p.4. l.17, p.5 l.5 and p.7 l.28-30 for more details about these changes in the manuscript).
This manuscript is submitted to a cryosphere journal, so many readers (myself included) don't have a very strong background in probabilistic methods, so it would be great to detail terms that are not common. Also, it is not clear how the emulator is calibrated (with how many runs, how are these runs chosen, etc.) and why it needs to be calibrated separately for each RCP scenario. I naively thought that the goal was to calibrate over a wide range of climate conditions, and then be able to investigate different scenarios easily, so without the addition of new physical simulations.
We agree that the manuscript should be made accessible to the largest possible audience. We tried to reduce technical jargon at a minimum while being consistent with a coherent probabilistic/uncertainty quantification framework.
The emulator is calibrated separately for each RCP scenario because we modelled the uncertainty in the climate forcing as a categorical variable (with one of the RCP scenarios).
Constructing an (PCE) emulator for a categorical variable is not easy. In addition, we think that it is more coherent to build a separate emulator for each RCP scenario because we do not intent to study the uncertainty to climate forcing for a general forcing.
The emulator is built for each of the model configurations (that is RCP scenario and sliding law) from an ensemble of 500 training points (hence 500 forward simulations) in the parameter space of the parameters F calv , F melt , E shelf , τ e and τ w with a maximin Latin hypercube sampling design.
See also answers below to questions about p.9 l.32 for additional details.
References and background material are sometimes missing, especially in the introduction and in the description of the changes made to the model, so it would be good to add some references and to provide more context.
We agree with this suggestion. We have added additional references and background material in the introduction and the description of the model (see answers to the specific comments for more details about these changes in the manuscript).
3 Specific comments -p.1 l.10: "except perhaps the relaxation times": Does it contribute to the uncertainty or not? It would be good to put clearer conclusions in the abstract, or to remove this part.
Based on our experimental set-up, the relaxation times do no contribute to the uncertainty (at least significantly). We have removed the word "perhaps" to have clearer conclusions in the abstract and left the discussion about the impact of the relaxation times on the uncertainty in the main text.
-p.1 l.12: maybe quantify how more dominant it becomes with the different warming scenarios.
We have added in the abstract that the contribution of the uncertainty in sub-shelf melting to the uncertainty in sea-level rise projections goes from 5-25 % in RCP 2.6 to more than 90 % in RCP 8.5.
-p.3 l.19: It would be good to define "confidence regions" and add some references.
We have added a physical interpretation of a confidence region in glaciology to make it clearer. We also have added references  regarding the definition and construction of confidence regions in uncertainty quantification analysis.
-p.4 l.6: Add references for the SIA and SSA approximations.
We have added references for the SIA (Hutter, 1983;) and SSA  approximations.
-p.4 l.11-12: I think it is little biased to simply say that using a flux-condition at the grounding line is an appropriate solution. There has been some recent research  demonstrating that this flux condition does not accurately represent confined ice shelves, which is the case for most ice shelves around Antarctica.
We agree that using a flux condition at the grounding line may be not appropriate to represent grounding-line migration notably in the presence of (strongly) buttressed ice shelves. We agree that a flux condition is an approximation (parameterisation) for the dynamics of the grounding line at coarse resolution but it has been shown to capture the steady-state behaviour of grounding lines under weak buttressing (Schoof 2007;Pattyn et al, 2012). In addition, 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 (especially the SSA model); yet, it cannot reproduce quantitatively changes in ice-sheet mass and sea-level rise contribution as determined from models with a higher level of complexity, such as the Blatter-Pattyn and the full Stokes models that include vertical shearing at the grounding line, especially for short transients (decadal time scales) as these flux conditions have been derived at steady state (Drouet et al., 2013;. We are aware that 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 metres), which is intractable for large-scale AIS modelling within an ensemble. Using a flux condition may also suffer limitations in the presence of buttressed ice shelves as suggested by . Yet, Reese et al. focus primarily on the physical representation of buttressing and not on the parameterisation of buttressing, which aims at translating a complex physical process into a simplified representation. Reese et al. also predict and discuss the appearance of non-physical negative buttressing, which is not allowed with parameterisations (Pollard et DeConto, 2012;. Furthermore, the analysis of Reese et al. is purely diagnostic, while many of the effects described in their paper vanish with an evolving ice sheet (as it is the case with the f.ETISh model).
We have written a new paragraph and added additional references Pattyn et al., 2012Drouet et al., 2013; to discuss in more detail the applicability (reproduce the migration of the grounding-line at coarse resolution for the SSA model) and limitations (short transients, buttressing) of using a flux condition at the grounding line.
-p.4 l.17: "the implementation of an ocean model": is there an ocean model implemented in f.ETISh or a representation of the ocean conditions? I think this part needs to be rephrased.
We have changed this sentence to specify that sub-shelf melt rates are computed with an oceanmodel coupler based on the Postdam Ice-shelf Cavity mOdel (PICO) ocean-coupler model  rather than a simpler parameterisation of sub-shelf melting Pollard and DeConto, 2012;de Boer et al., 2015;Cornford et al., 2016).
-p.4 l.18: How is the positive degree day model changed? Also add references.
In f.ETISh version 1.0 , the positive degree-day model was based on the implementation by , while it is based on the implementation by Janssens and Huybrechts (2000) in f.ETISh version 1.2 for numerical efficiency. We have changed the sentence to give more details about the description of atmospheric forcing in f.ETISh version 1.2 and added corresponding references. We now explain that the description of atmospheric forcing is based on a parameterisation of the changes in atmospheric temperature and precipitation rate (Huybrechts et al., 1998;Pollard and DeConto, 2012), a parameterisation of surface melt with a positive degree-day model (Janssens and Huybrechts, 2000) and the inclusion of meltwater percolation and refreezing ).
-p.4 l.24: The resolution of the model is 20 km, which can be understood for such computations. However, using the argument that the model is "essentially scale independent" is surprising, especially as Frank Pattyn, who is a co-author on this paper, published quite different conclusions following the MISMIP experiments (Pattynet al., 2012;. Furthermore 20 km represents only a few points for important ice shelves such as the Pine Island or Thwaites ice shelves, so it is not clear how accurately these glaciers can be modeled with such a resolution. We agree that stating that the model is "essentially scale independent" is certainly a very strong statement and such statement should be avoided. Yet,  has suggested that his model was able to produce rather scale-independent results for a resolution of the model around 20 km due to the flux condition that makes the model results rather independent from the resolution (at least for the SSA model Pattyn et al., 2012)) and allows to reproduce MISI in large-scale ice-sheet simulations (the flux condition has been shown to pass the tests on reversibility, which is a major criterion in establishing MISI ). However, such spatial resolution cannot capture properly certain mechanisms that control grounding-line migration such as bedrock irregularities and ice-shelf pinning points even using sub-grid parameterisations of these mechanisms. Hence, we may expect discrepancies between our results and results at higher spatial resolutions (<5 km), especially for important small glaciers such as Pine Island and Thwaites glaciers, which represent only a few grid points with a 20-km resolution. Yet, using a higher spatial resolution (even with an emulator and given our computational ressources) remains rather intractable for large-scale and long-term ice-sheet simulations and large-ensemble simulations. We expect the uncertainty in the results due to the choice of the resolution to be limited when compared to the uncertainty in the results due to the uncertainty in the model parameters and forcing.
We have changed this paragraph to discuss more thoroughly the choice and the limitations of using a 20-km resolution.
-p.5 l.5: "linked to the atmospheric forcing": How? What is the link used? A few words should be added to briefly describe this or at least mention that it will be detailed later on.
We have changed the sentence to make more explicit the link between the atmospheric forcing and the oceanic forcing. The link is explicitly written as T oc = T obs oc + F melt ∆T , where T oc is the ocean temperature on the continental shelf, T obs oc is the observed present-day ocean temperature on the continental shelf, F melt the ocean melt factor and ∆T the change in background atmospheric temperature. The equation for T oc 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  for changes in atmospheric and ocean temperatures .
-p.5 l.16: What is the impact of this correction? Can we consider the response to be relatively linear to subtract the reference runs? This needs a bit of clarification. This is a good question. Such initial drifts are not necessarily undesirable (given that the present-day Antarctic ice sheet is not at equilibrium) and simply subtracting the results from the initial drift may be a questionable assumption as it assumes that the model response to perturbations is totally independent of this initial drift. Yet, we have decided to ignore these initial drifts for practical experimental considerations as carried out similarly in other experimental set-ups Goelzer et al., 2018. The initial drifts for the reference runs reach a quasi-equilibrium at the end of the simulation and range from 0.1 to 0.2 metre of sea-level rise by 3000 for the different sliding conditions. The impact of the initial drift is more important for projections on a short-term time scale and has a limited impact on medium-term and long-term time scales. To isolate the impact of the sliding conditions on short-term projections and avoid any spurious impact of the initialisation procedure, we had to correct for the initial drifts, with differences between the reference runs for the different sliding laws comparable to the ordre of magnitude of the projections by 2100.
We have added a bit of clarification about this correction in Sect. 2.1 and Sect. 4.5 and its limitations.
-p.5 l.20: How about ice shelf basal melt, which is a dominant driver of the response of the Antarctic ice sheet?
We did not mention the uncertainty in sub-shelf melting in the list as this source of uncertainty is accounted for through the uncertainty in a physical model parameter (F melt ). We agree that this source of uncertainty should be written more explicitly in the sentence given its significance for the response of the Antarctic ice sheet. Hence, we have added the phrase "subshelf melting" in the sentence. We have changed the overall sentence to list the uncertain physical processes that we consider and added a new sentence that explains that we account for uncertainty in these physical processes by considering an ensemble of uncertain parameters in the f.ETISh model.
-p.6 Table 1: Caption: I don't think using "uncertain" is appropriate given that many other parameters are also uncertain. Consider changing. Also why use m=2 for the nominal value? I had the impression that the standard value is usually m=3.
We have replaced the phrase "uncertain parameters" by "Parameters with probabilistic representation" to avoid any ambiguity. We agree that the standard value is usually considered as m = 3 = n  with n the exponent in Glen's flow law and this value has been used in numerous studies (Schoof, 2007;Pattyn et al., 2012). Yet, a linear relationship (m = 1) is also used commonly in numerical simulations . Here, we have adopted an intermediate value (m = 2) between these two usual values (m = 1 and m = 3) as the nominal value of the sliding exponent in our simulations. We have added a new paragraph in Sect. 2.2.2 to clarify this (maybe unconvential) choice.
-p.6 l.2: I don't agree that the atmospheric forcing can cause large changes on the dynamics of the Antarctic ice sheet. It has a large impact on its volume and therefore sea level rise, but the impact on the dynamics is rather limited, unlike what is observed with changes in oceanic forcing.
We agree with this remark. We have changed the sentence to highlight that the atmospheric forcing is expected to cause changes in the mass balance of the Antarctic ice sheet rather than dynamical changes. We have added some references for this sentence .
-p.7 l.18: it would be interesting to reference the buttressing capabilities of ice shelves around Antarctica (Furst et al., 2016).
We have followed the referee's suggestion and added this reference.
-p.7 l.28-30: It would be good to explain briefly the link between atmospheric and ocean warming (warming of the ocean surface, at depth, ...), especially as changes in ocean circulation rather than ocean warming are expected to happen around Antarctica.
We have added a new paragraph that explains how ocean circulation results in sub-shelf melting and how changes in atmospheric and ocean warmings can affect ocean circulation. We have also added a new paragraph to explain how changes in atmospheric and ocean warmings are linked to sub-shelf melting in the PICO model ).
-p.9 l.32: Why is it necessary to build a different PC expansion for each RCP scenario? I would have imagined that the climate forcing was just a parameter that could be varied, similar to the other ones.
In our experimental set-up, the climate forcing (represented by the change in atmospheric temperature ∆T ) is represented as a discrete/categorical parameter with four different trajectories (the 4 RCP scenarios). This is similar to the sliding exponent, which can only take discrete values (m = 1, 2, 3 and 5) in our experimental set-up. Actually, we consider 20 distinct model configurations given by each combination of RCP scenario with a sliding law (m = 1, 2, 3 or 5) and each combination of RCP scenario with the TGL parameterisation. Building an emulator based on a PC expansion for discrete/categorical parameters/configurations is not straightforward and we prefer keeping them distinct for clarity.
-p.9 l.32: I don't really understand how the emulator is calibrated. Are the 500 forward simulations the training set or the results of the emulator? How many runs are used to calibrate the emulator and how are they chosen? It would be good to add some details in this section.
In our experimental setting, we consider 20 distinct model configurations given by each combination of RCP scenario with a sliding law (m = 1, 2, 3 or 5) and each combination of RCP scenario with the TGL parameterisation. An emulator is built for each of these model configurations from an ensemble of 500 training points (hence 500 forward simulations) in the parameter space of the parameters F calv , F melt , E shelf , τ e and τ w with a maximin Latin hypercube sampling design. In total, we carried out 10000 forward simulations of the f.ETISh model for the 20 model configurations. More details about the construction of the PC expansion are given in Appendix A. We have changed the paragraph and given more details about the construction of the PC expansion.
-p.10 l.14: "subregions" → "regions" or "areas" The word "subregions" has been replaced by "regions" -p.11 l.25-26: "making ice flux at the grounding line less important" ? "reducing ice flux at the grounding line" The sentence has been changed to follow the referee's suggestion.
-p.12 l.4-8: How different is it in this case? How significant are the changes for these additional simulations? It would be good to add numbers here, as the results are not showed.
It is difficult to draw general conclusions as the results are dependent of the characterisation of the uncertainty in the bedrock relation times. Yet, we found in RCP 2.6 that the uncertainty in the bedrock relation time for West Antarctica can account for 10 % of the uncertainty in ∆GMSL when τ w is allowed to vary widely between 50 years and 3500 years. We wrote this number in the text to give the reader an idea about the influence of τ w .
-p.12 l.9: remove "significantly" The word "significantly" has been removed from the manuscript.
-p.12 l.12: How about the order of magnitude difference in Fig.5 for example for several parameters? Where does this difference at different times come from? It would be good to discuss this.
As explained in Sec. 3.2. 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. This means in particular that F melt = 0.3 in Fig. 5(a-c) and Fig. 5(g-o) (this represents a background oceanic forcing). So the response of the Antarctic ice sheet is still driven strongly by ocean warming in RPC 8.5 (even if we assume no uncertainty in F melt ). This explains the difference of magnitude between the difference time scales as ocean warming triggers MISI over time.
-p.13 l.5: Why would it be cooler given the values showed on Fig.2? Beyond 2000 (start of our simulations), the RCP scenarios get warmer (RCP 8.5 warmer than RCP 6.0, RCP 6.0 warmer than RCP 4.5 and RCP 4.5 warmer than RCP 2.6). It is more likely to observe a decrease in sea-level rise in RCP 2.6 than for instance RCP 6.0 or RCP 8.5, that is, under cooler atmospheric conditions.
-p.13 l.7: Is it possible to separate the impact of the ocean and the atmosphere?
In our own opinion, it is not clear how to separate properly the impact of the ocean and the atmosphere because both are linked in our model. The atmospheric forcing (through the RCP scenario) and the ocean melt factor F melt both contribute to the uncertainty in the ocean temperature and therefore to the uncertainty in sub-shelf melting. Hence, the ocean forcing and the atmospheric forcing are dependent variables and their individual impact cannot be directly isolated. On the other hand, the five parameters F calv , F melt , E shelf , τ e and τ w are (assumed) independent and their individual impact on the projections can be isolated. It would be possible to separate the impact of the ocean and the atmosphere using another representation of these processes. However, the direct impact of the atmospheric forcing is quite limited as sub-shelf melting (indirect impact of the atmospheric forcing) is for instance an order of magnitude larger than surface melt. So basically, we witness response to ocean warming.
-p.13 l.16: Maybe add that this is because the forcing is so much larger for the high emission scenarios that the uncertainty caused by model parameters (considered similar in all the cases in your simulations) is therefore automatically reduced.
We have added that this result is due to a much larger increase in the values of the ∆GSLM projections compared to the increase in their dispersion when the RCP scenario gets warmer. Hence, the relative uncertainty in ∆GSLM projections is automatically reduced.
-p.13 l.26: "sliding conditions in RCP 2.6" → "sliding conditions. In RCP 2.6" This sentence has been corrected based on the referee's comment. The phrase "both regions" refer to the 100 % confidence region for grounded ice and the present-day grounded ice region. We have modified the sentence to make it clearer. "differences between both regions are only of few tenths of kilometres by 3000" has been replaced by "the 100 % confidence region for grounded ice by 3000 only differ from the present-day grounded ice region by a few tenths of kilometres".
-p.16 l.16: What is the distribution in these intervals? Is it always a uniform distribution?.
The distribution is still uniform. In the other sections, we considered a uniform distribution with support [F calv,min , F calv,max ]. In Sect. 3.6, the distribution is still uniform but with support [F calv,nom + α(F calv,min − F calv,nom ), F calv,nom + α(F calv,max − F calv,nom )] and the size of the support (thus the uncertainty) is controlled by the parameter α ∈ [0, 1]. Idem for the other parameters. We have changed the first paragraph in Sect. 3.6 to make it clearer.
-p.16 l.21: I don't understand why it would be different.
The emulator is an approximation to the parameters-to-projection relationship determined by the ice-sheet model. The difference between the response of the emulator and the ice-sheet model should of course be as small as possible.
- Table 4 and Table 5 captions: maybe re-write what are the probability intervals.
We have rewritten what are the probability intervals in the captions of Table 4 and Table 5.
-p.19 l.11-12: I don't think the results presented show that the atmospheric forcing has a large impact on the ice sheet dynamics, at least for 2100.
We agree that this sentence is a bit too general and does not reflect the results shown in the manuscript. Hence, we have changed it to point out that the impact of atmospheric forcing becomes important for medium-term and long-term time scales.
-p.19 l.21-22: I think this demonstrated mostly that the results are significantly impacted by the calibration of the basal melt forcing.
The referee is right in his/her interpretation. We have changed the sentence as "  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" to emphasise more clearly the impact of the calibration of the basal melt forcing.
-p.20 l.3: "the contribution of marine drainage basins": I don't understand what you mean here.
We mean that the significance of the contribution of the Antarctic ice sheet to sea-level rise under climate change is primarily controlled by the sensitivity, the response time and the vulnerability of its marine drainage basins. The phrase "the contribution of marine drainage basins" has been replaced by the more explicit phrase "the sensitivity, the response time and the vulnerability of its marine drainage basins" -p.21 l.4-5: Add some references here, this is an important point.
We have followed the referee's suggestion and added some references  -p.21 l.14: What is the impact of this choice? How different would the results be with a different choice for the construction of the emulator? This is a good point. We have checked that considering a PC expansion of order 4 instead of order 3 gave similar results to those in the manuscript. We may expect that for a PC expansion with a higher polynomial degree the results would be degraded due to the fitting of numerical noise (it has not been checked as the number of training points is insufficient to build a PC expansion with a too high polynomial degree) .
A discussion about a different choice for the construction of the emulator (maybe a Gaussian emulator) is beyond the scope of this paper but it would be worthwhile to gain insight into the efficiency of different choices of emulator in glaciology.
Based on our experimental set-up, the relaxation times do no contribute to the uncertainty (at least significantly). We have removed the word "perhaps" to have clearer conclusions.
The term g I is just the remaining term of the decomposition of g(x) in terms of the mean value g 0 and the main effects g i 's. It is simply given as (1) We did not bring attention to this term as it is not required to compute the Sobol indices (and the interaction index S I can be computed as S I = 1 − d i=1 S i ). We have added Eq. (1) in appendix B to clarify what is g I .
- Fig.1 caption: "neglected": Is it neglected because it is negligible or because the melt parameterization used does not allow to compute refreezing?
We agree that this sentence is a bit ambiguous. Here, we do not allow for refreezing under ice shelves. Still, refreezing rates are in general small compare to melt rates. We have changed the caption to clarify that refreezing is not taken into account in the model.
- Fig.3: It would be good to add the ice shelves on panels b to e.
We have added the ice shelves as suggested by the referee.
- Fig.4 and Fig.5: Is it possible to use the same axis for the three columns so that the comparison is easier?
We prefer keeping different axes for the three columns as our goal for these figures is to show the relationship between each individual input parameter and ∆GMSL rather than comparing Response to the Interactive comment on "Uncertainty quantification of the multicentennial response of the Antarctic Ice Sheet to climate change" by Kevin Bulthuis et al.
We would like to thank Andy Aschwanden for the time dedicated to this manuscript and his constructive comments to improve the general quality and readibility of the manuscript. We will try to give a proper response to his comments and revise the manuscript accordingly. For each referee's comment (written in blue), we included below a response (written in black) and proposed means to improve the manuscript.

Summary statement
This manuscript uses a framework comprising an ice sheet model, an emulator, and uncertainty quantification methods to assess the response of the Antarctic Ice Sheet to climate change. The paper is generally well written and boasts beautiful figures. Its strength lies in the comprehensive probabilistic approach, with a thoughtful use various uncertainty quantification methods. While I remain suspicious of the applicability of emulators, as they could miss discontinuities and strong non-linearities, I must also admit that I am not sufficiently familiar with the topic to have an informed opinion.
Best wishes,

Andy Aschwanden
We would like to thank Andy Aschwanden for this kind summary statement. The referee is right in pointing out that discontinuities and nonlinearities could pose challenges to the construction of emulators. However, we believe that our approach is valid when considering the global response of the Antarctic ice sheet (AIS) (e.g. ∆GMSL) as we may expect possible local and regional strong non-linearities and discontinuities to be smoothed out at the continental scale. Our belief is further backed by other studies of the continental response of the Antarctic ice sheet  that also suggest that this response does not exhibit strong non-linearities and discontinuities. Conversely, we did not compute the confidence regions with an emulator as the applicability of an emulator to the local/regional behaviour of the AIS is more questionable and we used directly the training points and forward simulations to estimate them. We leave for a further work the investigation of emulators for estimating confidence regions (Bulthuis et al., in preparation). We have added a sentence in the manuscript to clarify the use of an emulator to represent the continental response of the Antarctic ice sheet, as measured through ∆GMSL.

Major comments
I am a bit puzzled by the model setup. With reported computational costs of 16 CPU hours per simulation, one wonders if an emulator is indeed need as on the order of thousands simulations are computationally tractable? As discussed in the manuscript, the horizontal resolution of the ice sheet model is coarse, and topographic details will be missed. A demonstration of convergence under grid refinement would provide some confidence that the ice sheet model simulations are indeed robust. Nonetheless, one should not hold this too much against the manuscript, because, as I wrote above, the strength is clearly in the uncertainty quantification.
Main-effect Sobel indices could become a standard tool in ice sheet modeling. Aschwanden et al (under review) use a similar approach for Greenland.
We agree with the referee that a resolution of 20 km may be not able to capture properly certain mechanisms that control grounding-line migration such as bedrock irregularities and ice-shelf pinning points as well as important small glaciers such as Pine Island and Thwaites glaciers. We refer the referee to our response to anonymous referee # 1 regarding the limitations of this approach and the changes we brought at the manuscript to discuss the applicability and limitations of this approach (https://doi.org/10.5194/tc-2018-220-AC1).
The referee is also right in pointing out that our ice-sheet numerical model would allow to perform thousands of simulations, thus questioning the (full) need for an emulator. Yet, this approach would be feasible only for a few configurations of the model. In our manuscript, we investigated a set of 20 distinct configurations given by each combination of RCP scenario with a sliding law (m = 1, 2, 3 or 5) and each combination of RCP scenario with the TGL parameterisation in an attempt to investigate a sufficiently broad ensemble of model configurations. For each of the configurations, we ran a set of 500 forward simulations (training set) to build a PC emulator as performing thousands of simulations for each of the configurations would have been less computationally tractable. In addition, we had two definite goals in writing down this manuscript such as (1) providing the reader with a comprehensive probabilistic approach and various uncertainty quantification methods for uncertainty analysis in ice-sheet models (as highlighted by the referee) and (2) investigating a sufficiently broad ensemble of uncertainties to give a large insight into the impact of uncertainties in ice-sheet models.
Following the referee's suggestion, we have added a new figure (Fig. S1 in the supplementary material or Fig. 1 below) in the supplementary material that provides a comparison for the AIS contribution to sea level as a function of spatial resolution (20 km vs 16 km) to give an idea about the impact of the model resolution on our results. This figure suggests that the uncertainty in the projections due to the model resolution is (far) less important than the uncertainty due to the uncertainty in the parameters.
Why this confusing setup? p.10, l.26-26: "All results to follow have been obtained with the SGL parameterisation for the grounding-line migration, except in a discussion at the end of the section, where we discuss the impact of using the TGL parameterisation." Wouldn't it be easier to discuss the experimental design and results if the choice of grounding line parameterization were an additional parameters in the uncertainty quantification (i.e. by prescribing a uniform probability density function).
We agree with the referee that this sentence sounds a bit awkward and does not bring any additional information to the reader, so we have removed it from the manuscript. In addition, we make it clear in Sect. 2.2.3 that most results were obtained with Schoof's grounding-line parameterisation (reference parameterisation) and Tsai's grounding-line parameterisation is only discussed in Sect. 3.8. The choice of the grounding-line parameterisation can be seen as an additional parameter, but as this parameter is categorical (either Schoof's grounding-line parameterisation or Tsai's grounding-line parameterisation) it cannot be prescribed with a probability density function. We could have assigned a weight to the choice of the groundingline parameterisation, but we have no evidence for a suitable choice of weights.
3 Minor comments I find the introduction, discussion and use of "Uncertainty Community" somewhat awkward. While I think I understand what the authors try to express (and I like this approach), some clarification is warranted. What is the "recently formed uncertainty quantification community"? Maybe there is a white paper or similar that could be cited? UQ has been an integral part in numerical modeling outside glaciology but is only now making its way in ice sheet modeling (well, better late than never).
We agree that our terminology might sound a bit awkward, especially for most readers of the Cryosphere journal who are not familiar with uncertainty analysis. Uncertainty analysis has been a long standing part of numerical/experimental studies in science and engineering, but the coherent formalisation of uncertainty analysis in a probabilistic framework is rather recent and is still an ongoing work. The phrase "Uncertainty Community" refers to this community that aims at formalising uncertainty quantification analysis. In this manuscript, we have tried to bridge the gap between UQ outside of glaciology and glaciology. To avoid any ambiguity, we have decided to remove the phrase "Uncertainty Community" from the manuscript and talked more generally about UQ analysis rather than referring to the "Uncertainty Community".

Detailed comments
-Everywhere: say "the contribution to sea-level is..." instead of "the contribution to sea-level rise". "Rise" is not needed here.
We have followed the referee's suggestion and changed the phrase "the contribution to sea-level rise" by "the contribution to sea level" all along the manuscript.
We have changed the sentence in the manuscript based on the referee's suggestion.
We agree with the referee that this phrase sounds a bit awkward. We have replaced it by "as atmospheric and oceanic temperatures rise" following the referee's suggestion.
We have changed the sentence in the manuscript based on the referee's suggestion.
-p.2 l.24: . . . initial state, climate forcing, and parameters . . . We thank the referee for this suggestion. Yet, we tried not to use the Oxford (serial) comma all along the manuscript unless necessary to remain consistent in the manuscript (as suggested in the manuscript preparation guidelines for authors in the Cryosphere jornal).
We have rephrased "Uncertain parameters" as "List of parameters and parameter ranges used in the uncertainty analysis" following the referee's suggestion (The inappropriateness of the phrase "uncertain parameters" was also pointed out by anonymous referee #1 (https: //doi.org/10.5194/tc-2018-220-RC1)). We use the term "nominal value" to refer to the accepted/reference value of a parameter that we would have used in the absence of an uncertainty analysis, without any reference to an underlying probability distribution for the parameters.
In addition, we did not specify the probability distribution for the parameters in Table 1 as we only intended to discuss the sources of uncertainty in Sect. 2.2 and discuss the uncertainty quantification framework we used in Sect. 2.3 (notably the probabilistic characterisation of the parameters is given in Sec. 3.2.1).
We agree that assuming uniform probability distributions for the parameters may not be an appropriate assumption but as discussed in Sec. 3.2.1, a proper characterisation of the probability distributions of the parameters with expert assessment, data and statistical methods is beyond the scope on this paper. While a more refined uncertainty characterisation could be carried out, we used a uniform probability distribution for the parameters as this distribution is known to satisfy the principle of indifference (a.k.a the principle of maximum entropy) in the absence of any prior information about the distribution of the parameters except for their minimum and maximum values (Kapur, 1989), which we assume here for the sake of simplicity. Here, we determined the minimum and maximum values for the parameters based on expert assessment and literature.
-p.7 l.15: Please explain what F calv is. It appears to be a scalar multiplier of something (a calving rate, a stress condition)?
The referee is right in his interpretation of F calv . It is a scalar multiplier factor of the calving rate. Based on the referee's comment, we have changed Sect. 2.2.4 to give more details about the parameterisation of the calving rate in our model and the actual meaning of F calv (scalar multiplier factor of the calving rate).
-p.8 l.29: "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." OK, that's fine. For the exponent of the sliding law, one could perform a calibration and compare simulated and observed surface speeds to assess how well a certain exponent is able to explain the observations. See . This then be used as a prior for describing a PDF (done in Aschwanden et al, under review).
-p.12 l.18: you use "median" but provide a range. Please clarify. Do the ranges represent the 16/84th percentiles, for example? This range is simply the range of values taken by the median values for the three sliding exponents m = 1, 2, 3 as shown in Table 3. A formulation such that "a median AIS contribution to sea-level rise of 0.07-0.13 m in RCP 2.6 by 2300" takes into account the median value 0.07 m for m = 1, 0.13 m for m = 2 and 0.09 m for m = 3. We have changed our sentence following the referee's comment to avoid any ambiguity.
We agree with this suggestion. We have added the suggested reference at the end of the sentence.
-p.12 l.31: ". . . which contributes 3-3.5 metres to sea-level followed by a slower retreat of the East Antarctic ice sheet." We have changed the phrase following the referee's suggestion.
-p.12 l.33: "by the year 3000" We have changed the phrase following the referee's suggestion.
-p.13 l.4: ". . . an increase in sea-level, though a decrease cannot be ruled out for a viscous sliding law and cooler atmospheric conditions." We have simplified the sentence based on the referee's suggestion.
-p.14 l.6: It is interesting that for RCP 2.6 rheology contributes a similar amount 40-60 We agree with this comment. Maybe this suggests the need to properly calibrate E shelf alongside the sliding coefficient in the initialisation procedure.
We agree that using strong words like "certain" is probably not relevant and could lead to misinterpretations. Hence, we have removed the word "certain" as suggested by the referee.
-p.17 l.4-10: Rather than writing down numbers (which are listed in Table 4 anyway), I suggest to tell the reader how the plastic sliding law compares to the intermediate case since this is what one cares about.
We agree with the referee that all these numbers do not bring additional information to the manuscript (as they are listed in Table 4). We have removed all these numbers and written down the paragraph as a discussion about the difference with the other sliding laws.
-p.17 l.18: ". . . , which could explain the smaller ice loss in our results under m = 5 than m = 3." That's interesting, I would not have guessed.
We thank the referee for this comment.
-p.19 l.11: "the pivotal role played by atmospheric forcing". I think you mean the role of the emission scenario.
The referee is right in his interpretation. We have changed our phrase according to his suggestion to make it clearer.
-p.19 l.26: Moreover, the lower sensitivity of the Amundsen Sea sector in our simulations may arise. . .
We have changed the phrase "in our findings" by "in our simulations" as suggested by the referee.
-p.20 l.6-10: Very long sentence, maybe split into two. Figures We agree with this suggestion. The sentence has been split accordingly.
-p.20 l.23-29: I understand what you are trying to say but I'm not sure the formulation "does (not) question the nominal projections". Maybe "in agreement with" or similar? The probabilistic framework expands upon the nominal projections, and your results provide good evidence that a thorough risk assessment must include UQ.
The referee is right in his interpretation of this sentence. We agree with his suggestion to change the formulation of the sentence. We have changed the sentence as ". . . accommodating parametric uncertainty in the ice-sheet model leads to projections in agreement with the nominal projections of limited ice loss . . . " -p.21 l.20: We found that all investigated sources . . . We agree with this suggestion. The sentence has been changed accordingly.
- Fig. 3 b-e: Show outline of present-day grounded area for better visual comparison.
We have added the present-day grounding-line position as suggested by the referee.
- Fig. 9: What does the m in the lower left corners of the plot mean?
The parameter m is the sliding exponent in Weertman's sliding law. We have replaced the phrase "under different sliding laws" by "for different values of the sliding exponent m in Weertman's sliding law" to remind the reader about the meaning of m.