the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Uncertainty quantification of the multicentennial response of the Antarctic ice sheet to climate change
Kevin Bulthuis
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 icesheet models remain challenging due to the complexity of physical processes involved in icesheet dynamics, including instability mechanisms that can destabilise marine basins with retrograde slopes. Moreover, uncertainties in icesheet models limit the ability to provide accurate sealevel rise projections. Here, we apply probabilistic methods to a hybrid icesheet model to investigate the influence of several sources of uncertainty, namely sources of uncertainty in atmospheric forcing, basal sliding, groundingline flux parameterisation, calving, subshelf melting, iceshelf 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 sealevel rise and groundingline 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 icesheet 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)  Fulltext XML

Supplement
(4790 KB)  BibTeX
 EndNote
The Antarctic ice sheet (AIS) is the largest reservoir of freshwater on Earth (∼60 m sealevel 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 sealevel 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 icesheet 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 icesheet models amenable to largescale and longterm 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 icesheet 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 subgrid melt interpolation at the grounding line, Golledge et al. (2015) have projected that sealevel 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 icecliff instability mechanism (Pollard et al., 2015) in a numerical icesheet 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 icesheet 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 icecliff instability mechanism. Coupling an icesheet 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 icesheet dynamics (Pattyn et al., 2017), differences in modelling hypotheses between icesheet models remain a major source of uncertainty for sealevel 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 socalled structural uncertainty in icesheet 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 icesheet initial state, climate forcing and parameters in numerical icesheet 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 seasurface 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 sealevel rise projections by running an icesheet model with samples of parameters drawn randomly from probability density functions for the parameters and subsequently deducing probability density functions for sealevel 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 icesheet 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 icesheet model that captures the essential characteristics of ice flow and allows largescale and longterm 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 icesheet 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 multicentennial response of the whole Antarctic ice sheet: whereas these previous studies provided projections by running the icesheet 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, groundingline flux parameterisation, calving, subshelf melting, iceshelf 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 icesheet models: whereas these previous studies analysed AIS paleoclimatic responses with Gaussian process (kriging) emulators (Pollard et al., 2016) and multidecadal forecasts with Latin hypercube sampling (Schlegel et al., 2018), our study differs by its focus on the analysis of the multicentennial 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 multicentennial 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 Icesheet 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 icesheet/iceshelf model that incorporates essential characteristics of icesheet thermomechanics and icestream flow, such as the massbalance feedback, bedrock deformation, subshelf melting and calving. The ice flow is represented as a combination of the shallowice (SIA) (Hutter, 1983; Greve and Blatter, 2009) and shallowshelf (SSA) (Morland, 1987; MacAyeal, 1989; Weis et al., 1999) approximations for grounded ice (Bueler and Brown, 2009), while only the shallowshelf approximation is applied for floating ice shelves. Bedrock deformation is represented as a combined timelagged 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 largescale stress field, represented by the horizontal divergence of the iceshelf velocity field (Pollard and DeConto, 2012a). Prescribed input data include the presentday icesheet 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 groundingline 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 steadystate 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 steadystate 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 largescale icesheet 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 icesheet dynamics as determined with other icesheet models, it cannot reproduce quantitatively changes in icesheet 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 groundingline 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 groundingline 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 groundingline migration such as bedrock irregularities and iceshelf pinning points even with subgrid 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 groundingline migration in largescale, longterm and largeensemble icesheet 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.SFNRS & 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 subshelf melt rates with an oceanmodel coupler based on the Potsdam Iceshelf Cavity mOdel (PICO) oceanmodel coupler (Reese et al., 2018a) rather than simpler parameterisations of subshelf 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 degreeday 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. Presentday 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 degreeday (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 oceanmodel coupler (Reese et al., 2018a), which evaluates subshelf 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 iceshelf cavities. We employ data by Schmidtko et al. (2014) for presentday 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 subshelf melt rate as computed with the PICO model. The PICO model is able to reproduce the general pattern of subshelf melting with higher melt rates near the grounding line and lower melt rates near the calving front. Subshelf 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 icesheet geometry by Pollard and DeConto (2012b). This approach is based on a fixedpoint iteration scheme that adjusts the basal sliding coefficient iteratively so as to match the presentday icesheet configuration while assuming that the presentday 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 groundingline 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 sealevel 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 presentday dynamical changes, our results reflect the model response to perturbations. These corrections have in general little impact on mediumterm and longterm projections as the initial drift is small compared to the projections but leads to a more significant bias for shortterm 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, groundingline flux parameterisation, calving, subshelf melting, iceshelf 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 icesheet 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 groundingline 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 groundingline recession due to reduced buttressing. DeConto and Pollard (2016) have also highlighted that the hydrofracturing and icecliff failure mechanisms (not included in f.ETISh version 1.2) driven by increased surface melt and subshelf melting could potentially lead to an accelerated collapse of the West Antarctic ice sheet and a deeper groundingline retreat in the East Antarctic subglacial marine basins.
2.2.2 Basal sliding
Basal sliding controls the motion of fastflowing 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 powerlaw 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 iceflow 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 iceflow models. Hereafter, we refer to the exponents $m=\mathrm{1},\mathrm{2}$ and 3 as the viscous (or linear), the weakly nonlinear and the strongly nonlinear sliding law, respectively. We consider the discrete values $m=\mathrm{1},\mathrm{2}$ and 3 as representative of the most common values in largescale icesheet modelling and discuss in Sect. 3.7 the impact of a more plastic sliding law (m=5) to represent a quasiplastic deformation of the till in ice streams (GilletChaulet et al., 2016).
2.2.3 Groundingline flux parameterisation
The f.ETISh model employs a parameterisation of the groundingline 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 nonlinear 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 presentday 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 iceshelf 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 icecliff 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 subgrid ice thickness within a fraction of the iceedge 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
Subshelf melting is mainly controlled by subshelf ocean circulation, which can be affected by atmospheric changes. Iceshelf thinning caused by increased subshelf melting leads to a reduction in iceshelf 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 iceshelf cavities (Hellmer et al., 2012; Schmidtko et al., 2014). Changes in ocean circulation resulting in stronger subiceshelf 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 subiceshelf circulation, thus leading to an increase in subshelf melting. Yet, the link between global climate change and changes in subshelf melting is not always clear: it has been suggested that future climate change could lead to an increase in subshelf melting (Hellmer et al., 2012; Timmermann and Hellmer, 2013), a positive meltwater feedback that enhances subiceshelf circulation and can trigger a climate tipping point (Hellmer et al., 2017) and even to a decrease in subshelf melting (Dinniman et al., 2012). Golledge et al. (2019) have also highlighted the need to couple icesheet 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 iceshelf 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 presentday 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 subshelf 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 Iceshelf rheology
Ice rheology in largescale icesheet 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 strainrate 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 sheardominated stress regime for grounded ice to an extensiondominated stress regime for floating ice.
We introduce an iceshelf 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 icesheet 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 GilletChaulet 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 sealevel 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 nonintrusive manner treating the icesheet 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 icesheet 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 parameterstoprojection relationship as a polynomial in the parameters. We write this polynomial as a linear combination of polynomial basis functions and use leastsquares regression (Appendix A) to evaluate the coefficients from a limited number of icesheet 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 parameterstoprojection relationship to be sufficiently smooth (no discontinuity or highly nonlinear behaviour) to allow an efficient approximation as a lowdegree polynomial, and PC expansions may be inefficient in highdimensional problems.
The emulator of the relationship between the parameters and the projection is then used as a substitute for the icesheet 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 icesheet model that accelerates UQ; (ii) it provides an explicit view of the relationship between the parameters and the projection, highlighting potential linear or nonlinear 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 nonlinear 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 crossvalidation 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 icesheet response. Here, we adopt the variancebased sensitivity indices (Saltelli et al., 2008), also called Sobol indices, described in more detail in Appendix B. Variancebased 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 groundedice 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 icesheet 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 icesheet 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 shortterm (2100), mediumterm (2300) and longterm (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 sealevel 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 Parameterstoprojection relationship
One of the advantages of a polynomial chaos expansion is that it provides an explicit approximation to the parameterstoprojection 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 (oneatatime) while keeping the other parameters fixed at their nominal value under the weakly nonlinear 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 nonlinearly. 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 shortterm and mediumterm 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 iceshelf 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 sealevel 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 longterm 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 Sealevel 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 sealevel 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 icesheet model exhibits nonlinearities, 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 nonlinear 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 sealevel 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 nonlinear 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 subshelf 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 nonlinear 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 sealevel 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 sealevel 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 nonlinear 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 shortterm, mediumterm and longterm 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 iceshelf rheology (Sobol indices ranging from 40 % to 60 %) followed by the uncertainty in the calving rate (Sobol indices ranging from 20 % to 40 %) and subshelf melting (Sobol indices ranging from 5 % to 25 %). Indeed, in RCP 2.6, subshelf 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 iceshelf 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 subshelf melting, which accounts in RCP 8.5 for more than 90 % of the uncertainty in sealevel 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 groundedice 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 nonlinear 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 fastflowing 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 groundingline 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 presentday position. The risk of retreat in Wilkes basin may partially explain the acceleration in sealevel rise that we observed in Fig. 6 in RCP 8.5. The risk of groundingline 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 presentday grounded ice region (the 100 % confidence region for grounded ice by 3000 only differs from the presentday 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 nonlinear sliding law with projections under the other sliding laws, that is, the viscous sliding law (Fig. S7) and the strongly nonlinear 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 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. However, Figs. 10i and 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 in the former region to nonlinearity 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 icesheet model without parametric uncertainty, that is, the icesheet 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 sealevel 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 icesheet model.
For the longterm 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 nonlinear 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 nonlinearly 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 sealevel rise projections on shortterm and mediumterm timescales, thus suggesting a more significant and faster response to perturbations under more plastic sliding conditions. On a longterm timescale, the icesheet 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 GilletChaulet 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 nonlinearities 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 nonlinear 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 sealevel rise projections compared to our results under the SGL parameterisation. This result was expected as the TGL parameterisation has been shown to increase groundingline 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 sealevel 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 groundingline retreat in the marine sectors and an additional mass loss from East Antarctica, especially in the Aurora basin.
4.1 Comparison of the sealevel rise projections with previous work
Regarding the shortterm 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 subgrid 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 sealevel 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 sealevel targets between 10 and 20 m. In all RCP scenarios, under all sliding laws and under both groundingline 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 subgrid 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 Buddtype 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 longterm 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 longterm projections are generally below projections by DeConto and Pollard (2016) with hydrofracturing and icecliff 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 sealevel 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 subshelf melting as critical drivers of the future changes in the AIS mass balance on mediumterm and longterm timescales through iceshelf 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 groundingline 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 groundedice retreat with previous work
Consistent with our results, Golledge et al. (2015), with a 10 km resolution model, projected that groundingline 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 subkilometre resolution around the grounding line, projected that groundingline retreat is most significant in the Amundsen Sea sector. Schlegel et al. (2018) found that groundingline retreat is most significant in the Amundsen Sea sector under generalised ocean warming experiments for the Antarctic ice sheet, but, after calibrating subshelf 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 subshelf 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 groundingline 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 groundingline retreat is determined by bedrock topography and iceshelf buttressing which depends on the importance of iceshelf 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 icesheet model leads to projections in agreement with the nominal projections of limited ice loss and limited groundingline 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 icesheet model leads to projections in disagreement with the nominal projections and indicates instead some risk of triggering accelerated ice loss and significant groundingline 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 icesheet model. The f.ETISh model is an icesheet model that focuses on essential marine icesheet mechanisms, similarly to the icesheet model by Pollard and DeConto (2012a). Certain physical processes, especially smallscale 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, groundingline migration and thresholds for instabilities may not be captured properly even with a parameterisation of the groundingline flux, especially in lowforcing scenarios and shortterm 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 groundingline remains an acceptable assumption in largescale and longterm icesheet simulations and largeensemble 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, subshelf melting may not be captured properly despite the use of the PICO oceanmodel 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 subshelf melting through the ocean melt factor. However, as discussed in Sect. 2.2.5, atmospheric forcing affects the Circumpolar Deep Water circulation in iceshelf cavities and subsequently modifies subshelf melting (Pritchard et al., 2012). Still, the fate of the Southern Ocean and the evolution of subshelf 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 subshelf melting in driving the future response of the Antarctic ice sheet, there is a clear need to better constrain future subshelf melt and incorporate it properly into icesheet 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 shortterm 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 parameterstoprojection relationship, but not local variations in the parameter space that may stem from numerical noise and errors.
We studied the multicentennial response of the Antarctic ice sheet under uncertainty using methods from the field of UQ. We investigated uncertainties in atmospheric forcing, basal sliding, groundingline parameterisation, subshelf melting, calving, iceshelf rheology and bedrock relaxation. We used emulationbased methods to represent the parameterstoprojection 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 groundingline 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 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 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 icesheet 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 icesheet 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, $\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 illconditioning of the leastsquares 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 spacefilling 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 leastsquares approach is its good ability to handle numerical noise in lowdegree 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 icesheet 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 crossvalidation and convergence tests. We generate a new set of samples in the parameter space and we compare the exact response of the icesheet model with the approximate response of the PC expansion. Figure A1a–c show results of such a crossvalidation 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 meansquared error between the exact response of the icesheet 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 socalled main effect associated with the corresponding parameter x_{i} and g_{I} the socalled 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 leastsquares best approximation of g among all functions that depend only on x_{i}. As a consequence of this leastsquares 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 ith 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 ith 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 icesheet 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 nonconstant 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 socalled 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 nonGaussian 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/tc1313492019supplement.
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 coauthors.
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.SFNRS) 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/se55692014, 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, lithosphereasthenosphere 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., RoviraNavarro, 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 icesheet 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/s14635003(02)000197, 2003. a
Bennett, M. R.: Ice streams as the arteries of an ice sheet: their mechanics, stability and significance, EarthSci. Rev., 61, 309–339, https://doi.org/10.1016/s00128252(02)001307, 2003. a
Bindschadler, R. A., Nowicki, S., AbeOuchi, 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.: Icesheet 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/tc719492013, 2013. a, b
Brondex, J., Gagliardini, O., GilletChaulet, 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., GilletChaulet, 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/tc131772019, 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.: Centuryscale simulations of the response of the West Antarctic Ice Sheet to a warming climate, The Cryosphere, 9, 1579–1600, https://doi.org/10.5194/tc915792015, 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., AbeOuchi, A., Saito, F., and van de Wal, R. S. W.: Simulating the Antarctic ice sheet in the latePliocene warm period: PLISMIPANT, an icesheet model intercomparison project, The Cryosphere, 9, 881–903, https://doi.org/10.5194/tc98812015, 2015. a
DeConto, R. M. and Pollard, D.: Contribution of Antarctica to past and future sealevel 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/jclid1100307.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/s1071201191333, 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/tc73952013, 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 icecliff instability, Nature, 566, 58–64, https://doi.org/10.1038/s4158601909014, 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., GilletChaulet, F., Zwinger, T., Payne, A. J., and Le Brocq, A. M.: Retreat of Pine Island Glacier controlled by marine icesheet 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., RigerKusk, 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/tc73752013, 2013. a, b, c, d
Fürst, J. J., Durand, G., GilletChaulet, 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/9783319123851, 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
GilletChaulet, 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/tc813932014, 2014. a
Gladstone, R. M., Warner, R. C., GaltonFenzi, B. K., Gagliardini, O., Zwinger, T., and Greve, R.: Marine ice sheet model performance depends on basal sliding physics and subshelf melting, The Cryosphere, 11, 319–329, https://doi.org/10.5194/tc113192017, 2017. a
Goelzer, H., Nowicki, S., Edwards, T., Beckley, M., AbeOuchi, A., Aschwanden, A., Calov, R., Gagliardini, O., GilletChaulet, 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 initMIPGreenland: an ISMIP6 intercomparison, The Cryosphere, 12, 1433–1460, https://doi.org/10.5194/tc1214332018, 2018. a
Golledge, N. R., Kowalewski, D. E., Naish, T. R., Levy, R. H., Fogwill, C. J., and Gasson, E. G. W.: The multimillennial Antarctic commitment to future sealevel 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 twentyfirstcentury icesheet melt, Nature, 566, 65–72, https://doi.org/10.1038/s4158601908899, 2019. a, b
Gomez, N., Mitrovica, J. X., Huybers, P., and Clark, P. U.: Sea level as a stabilizing factor for marineicesheet grounding lines, Nat. Geosci., 3, 850–853, https://doi.org/10.1038/ngeo1012, 2010. a
Gomez, N., Pollard, D., and Mitrovica, J. X.: A 3D 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/tc1222292018, 2018. a
Greve, R. and Blatter, H.: Dynamics of Ice Sheets and Glaciers, Advances in Geophysical and Environmental Mechanics and Mathematics, SpringerVerlag Berlin Heidelberg, https://doi.org/10.1007/9783642034152, 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.: Twentyfirstcentury warming of a large Antarctic iceshelf 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/jclid160420.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/9789401511674, 1983. a
Huybrechts, P. and de Wolde, J.: The Dynamic Response of the Greenland and Antarctic Ice Sheets to MultipleCentury Climatic Warming, J. Climate, 12, 2169–2188, https://doi.org/10.1175/15200442(1999)012<2169:tdrotg>2.0.co;2, 1999. a, b
Huybrechts, P., AbeOuchi, 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 largescale 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 oilspill 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 massbalance 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 oceaninduced 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 sealevel projections at a global network of tidegauge 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.: MetamodelBased 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/9789048135202, 2010. a, b, c, d
Lenaerts, J. T. M., Vizcaino, M., Fyke, J., van Kampenhout, L., and van den Broeke, M. R.: Presentday 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/s0038201529074, 2016. a
Ma, Y., Gagliardini, O., Ritz, C., GilletChaulet, F., Durand, G., and Montagnat, M.: Enhancement factors for grounded ice and ice shelves inferred from an anisotropic iceflow model, J. Glaciol., 56, 805–812, https://doi.org/10.3189/002214310794457209, 2010. a
MacAyeal, D. R.: LargeScale 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/tc813472014, 2014. a, b, c
Moholdt, G., Padman, L., and Fricker, H. A.: Basal mass budget of Ross and FilchnerRonne 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 IceShelf 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/9789400937451_6, 1987. a
Nowicki, S. M. J., Payne, A., Larour, E., Seroussi, H., Goelzer, H., Lipscomb, W., Gregory, J., AbeOuchi, A., and Shepherd, A.: Ice Sheet Model Intercomparison Project (ISMIP6) contribution to CMIP6, Geosci. Model Dev., 9, 4521–4545, https://doi.org/10.5194/gmd945212016, 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.14679868.2004.05304.x, 2004. a
Pattyn, F.: A new threedimensional higherorder 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.: Sealevel response to melting of Antarctic ice shelves on multicentennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0), The Cryosphere, 11, 1851–1878, https://doi.org/10.5194/tc1118512017, 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/tc65732012, 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.: Groundingline migration in planview marine icesheet 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 IceSheet Dynamics, Curr. Clim. Change Rep., 3, 174–184, https://doi.org/10.1007/s4064101700697, 2017. a
Pattyn, F., Ritz, C., Hanna, E., AsayDavis, 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/s4155801803058, 2018. a, b
Petra, N., Martin, J., Stadler, G., and Ghattas, O.: A Computational Framework for InfiniteDimensional 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 sheetshelf model, and application to Antarctica, Geosci. Model Dev., 5, 1273–1295, https://doi.org/10.5194/gmd512732012, 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/tc69532012, 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/gmd916972016, 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 icesheet 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., AsayDavis, X., and Winkelmann, R.: Antarctic subshelf melt rates via PICO, The Cryosphere, 12, 1969–1985, https://doi.org/10.5194/tc1219692018, 2018a. a, b, c, d
Reese, R., Winkelmann, R., and Gudmundsson, G. H.: Groundingline 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/tc1232292018, 2018b. a
Rignot, E., Jacobs, S., Mouginot, J., and Scheuchl, B.: IceShelf 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 sealevel rise from Antarctic icesheet 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/9781475741452, 2013. a
Rommelaere, V. and Ritz, C.: A thermomechanical model of iceshelf 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 100year contribution to sea level rise and associated model uncertainties using the ISSM framework, The Cryosphere, 12, 3511–3534, https://doi.org/10.5194/tc1235112018, 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 icesheet 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/tc1122832017, 201 a
Schäfer, M., Zwinger, T., Christoffersen, P., GilletChaulet, 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/tc67712012, 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., SandbergSø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/s415860180179y, 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/bamsd1100094.1, 2012. a, b
Timmermann, R. and Hellmer, H. H.: Southern Ocean warming and increased ice shelf basal melting in the twentyfirst and twentysecond centuries based on coupled iceocean finiteelement modelling, Ocean Dynam., 63, 1011–1026, https://doi.org/10.1007/s1023601306420, 2013. a, b
Timmermann, R., Wang, Q., and Hellmer, H.: Iceshelf basal melting in a global finiteelement seaice/iceshelf/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 icesheet 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/tc1238612018, 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