A high-end sea level rise probabilistic projection including rapid Antarctic ice sheet mass loss

The potential for break-up of Antarctic ice shelves by hydrofracturing and following ice cliff instability might be important for future ice dynamics. One recent study suggests that the Antarctic ice sheet could lose a lot more mass during the 21st century than previously thought. This increased mass-loss is found to strongly depend on the emission scenario and thereby on global temperature change. We investigate the impact of this new information on high-end global sea level rise projections by developing a probabilistic process-based method. It is shown that uncertainties in the projections increase when including the temperature dependence of Antarctic mass loss and the uncertainty in the Coupled Model Intercomparison Project Phase 5 (CMIP5) model ensemble. Including these new uncertainties we provide probability density functions for the high-end distribution of total global mean sea level in 2100 conditional on emission scenario. These projections provide a probabilistic context to previous extreme sea level scenarios developed for adaptation purposes.


Introduction
Sea level rise is one of the main consequences of global warming (Arnell et al 2016). Knowing how fast it can develop given a scenario of future greenhouse gas emissions is crucial for both mitigation and adaptation choices (Wong et al 2014). Unfortunately there remains considerable uncertainty. The Antarctic ice sheet is potentially the largest contributor to future sea level rise but also the most uncertain (Meehl et al 2007, Church et al 2013, Levermann et al 2014, Ritz et al 2015, Deconto and Pollard 2016. In 2013 the fifth assessment report (AR5) from the Intergovernmental Panel on Climate Change (IPCC) assessed the likelihood for an extensive grounding line retreat of the Antarctic ice sheet, that would contribute significantly to sea level rise, to less than 34% (Church et al 2013). In such a case, there was a medium confidence that the magnitude would be 'several tenths of a meter' . However, for long term projects that have a high risk aversion, low probability events also need to be taken into account (Veerman 2008, Ranger et al 2013, Hinkel et al 2015. This can be done in different ways: convening an expert committee to develop extreme scenarios (Katsman et al 2011, Ranger et al 2013, conducting a large expert assessment survey  or combining expert assessment of ice sheet contribution (Bamber andAspinall 2013, de Vries andvan de Wal 2015) with climate models projections (Jevrejeva et al 2014, Grinsted et al 2015. It is difficult with these approaches to capture the correlation between ice sheet mass loss and all the other processes. Also important subjective choices are involved in each of these methods (de Vries and van de Wal 2015). Until recently such choices were unavoidable as climate projections with an ice sheet model were either not available at all, or carried out by models that did not include processes that become important when identifying the high end of the distribution. However, as ice sheet models continue to improve and include new processes it has now become timely to carry out a probabilistic assessment of the high end of the distribution. Therefore we propose here an alternative method where extreme mass loss from numerical ice sheet simulations is used in a 'process-based' method (Church et al 2013).
Recently, Deconto and Pollard (2016) included the potential for break-up of Antarctic ice shelves in a dynamical ice model showing that Antarctica could contribute to global mean sea level by up to 114 cm with a 1s of 36 cm in 2100 relative to 2000 under the Representative Concentration Pathway RCP8.5 scenario. This estimate is much higher than the IPCC AR5 upper bound of the likely range, which was set at 14 cm. Another important conclusion from Deconto and Pollard (2016) is that future mass losses strongly depend on the emission pathway, with an order of magnitude less melt occurring in the RCP2.6 scenario. In particular, Pollard et al (2015) and Deconto and Pollard (2016) argued that ice fracturing of the ice shelves by surface melt (hydrofracturing) and ice cliff failure when cliffs become too high could increase Antarctic mass loss dramatically after they included these effects in their dynamical ice sheet model. These findings need to be taken with caution because they are not confirmed by other studies (Clark et al 2015). Nevertheless the Deconto and Pollard (2016) projections haven't been demonstrated to be unphysical. Therefore it is timely to explore the impact of rapid mass loss from the Antarctic ice sheet on global mean sea level rise and to update the techniques to project future total sea level rise in a high-end scenario.
To be able to construct a probability density function (PDF) instead of a likely range the main two issues are (1) the quantification of the uncertainties in ice-sheet contribution and (2) the uncertainty in climate models projections (Meehl et al 2007, Church et al 2013. Here, we will discuss the modifications of the process-based method that are necessary to quantify the probability of extreme sea level rise. The objective of such a probabilistic scenario is to update previous high end estimates based on expert judgment (Katsman et al 2011, Ranger et al 2013, Sweet et al 2017 and to provide a probability density function that feeds into the policy requirements to assess the likelihood of high risk, low probability events.

Method
The 'process-based' method used to compute total sea level rise projections is based on IPCC AR5 (see supplementary material of Church et al (2013)). We take the same processes into account: ocean dynamics and steric expansion, glaciers and ice caps, Antarctic and Greenland ice sheets and land water storage. The estimates for all processes are also the same as in Church et al (2013) but we successively introduce three key differences to change the focus to high-end projections. Firstly, we look at the impact of the new Antarctic mass loss estimates from Deconto and Pollard (2016) on the global mean total sea level PDF. Secondly, we develop a method to make the ice sheet mass loss dependent on temperature for a given emission scenario. This is important because it broadens the PDF of total sea level rise and therefore increases the values of the high quantiles. Thirdly we explore the implication of the uncertainty related to confidence in climate model projections of global mean temperature. This uncertainty also leads to increase the high quantiles of the PDF. For simplicity we only focus on the sea level rise in 2100 compared to the reference period 1986-2005 and not on the pathway leading there. A detailed description of each process contributing to sea level and a mathematical description of the method to compute the PDFs is given in appendix A.

Antarctic ice sheet mass loss
Ice sheet mass loss may be decomposed into surface mass balance, the difference between snow fall and melt/sublimation, and ice-dynamics (iceberg calving, and basal melt of ice sheet and ice shelves). With climate change the snow fall may be considered to increase with temperature, Church et al (2013) assumed 5.1% per C. Changes in dynamics were based on the probabilistic framework of Little et al (2013), represented with a uniform distribution with a minimum of À2 cm and a maximum of 18.5 cm global sea level equivalent in 2100 relative to 2005.
To estimate the total Antarctic ice sheet mass loss, we use the results from a recent numerical model simulation (Deconto and Pollard 2016) (referred later as DP16) that includes both surface mass balance and ice-dynamics. In their high-end estimate, the expected value for the Antarctic contribution in 2100 relative to 2000 are 58 cm and 114 cm, with standard deviations 28 cm and 36 cm respectively for the RCP4.5 and RCP8.5 scenarios. These numbers represent the most extreme outcome from this study, other hypothesis were also tested and lead to differences (see appendix B). The standard deviations stated above represent the uncertainty in three important parameters of an ice flow model: ocean melting under the ice shelves, hydrofracturing due to surface melting and horizontal wastage due to ice-cliff structural failure. A PDF can then be constructed for the Antarctic contribution to sea level rise at year t : with X DP16,M the expected value, X DP16,SD the standard deviation, r randomly chosen from a normal distribution with zero mean and standard deviation 1 ðN ð0; 1ÞÞ and Sce that can either be the RCP4.5 or the RCP8.5 scenario. We note that the reference period from Deconto and Pollard (2016) is slightly different from the one used by the IPCC AR5, but mass loss in Antarctica from 1996 to 2005 is very small (0.25 cm, Church et al (2013)) so we do not take into account this difference in our computation. These projections of Antarctic mass loss are a lot higher and more uncertain (wider PDF) than the ones considered by Church et al (2013) (figure 1). This is mostly because Church et al (2013) did not try to determine an upper end of Antarctic mass loss since the simulations to do it were not available and their estimates are meant to be a best guess instead of a high-end scenario. As a result for RCP8.5, the median of total sea level rise was 73 cm, but in this high-end estimate it becomes 184 cm and the impact is even higher for the extremes with for example an increase of the 95% quantile from 98 cm to 247 cm (table 1). It should be stresses that by focusing on the high-end distribution we do not undertake the same exercise as in Church et al (2013), as we are not considering other sources of information for Antarctic ice sheet mass loss than DP16.

Temperature dependence of antarctic mass loss
When using the process-based method, the way processes are added is important, especially when computing the extremes. If two or more processes are independent the PDF of their sum has a smaller spread than if they are correlated. For example, when two independent random variables are added, an extreme in one of them, defined here as above the 95% quantile, will only happen together with an extreme of the other variable 5% Â 5% of the time. That is a 0.25% likelihood. In contrast, if these two variables are perfectly correlated the same extreme has a likelihood of 5%.
For each process contributing to sea level rise the uncertainty can be separated into two parts, a methodological part and a part depending on temperature. Our approach concerning the correlation between processes is based on the nature of their uncertainties. The time series of global mean surface air temperature (T ) comes from the CMIP5 climate models and is expressed as an anomaly with respect to its time-mean for 1986-2005. Using the ensemble mean T M and standard deviation T SD we write: with r 2 randomly sampled from N ð0; 1Þ. This is the scalar form of equation (5) discussed in appendix A.
The uncertainty that depends on temperature is fully correlated between all processes while the methodological parts are independent of each other. Processes that have some temperature dependence are land glaciers and ice caps and Greenland surface mass balance. Also the ocean thermal expansion is considered fully correlated with temperature. Greenland ice 1 4 4 3 6 9 6 7 1 4 1  5  51  46  121  104  81  10  56  51  135  121  103  20  62  58  152  143  131  50  73  73  184  184  184  80  85  90  216  225  238  90  92  99  233  247  268  95  98  108  247  265  292  99  111  127  273  299  339 Environ. Res. Lett. 12 (2017) 044013 sheet dynamics and land water storage are considered fully uncorrelated with temperature. In IPCC AR5 the Antarctic surface mass balance is considered temperature dependent but not the Antarctic ice sheet dynamics. Here is where we deviate from AR5, based on Deconto and Pollard (2016) we consider the Antarctic ice sheet mass loss linearly related to global mean surface temperature for a given emission scenario, assuming the following relationship: : For complete consistency with Deconto and Pollard (2016) the mean temperature used to compute d should come from their climate model, GENESIS v3 Global Climate Model. However this model it not part of CMIP5, therefore we have taken the ensemble mean of the CMIP5 models instead. Making the mass loss from Antarctica temperature dependent also makes it partially correlated to all the other processes that depend on temperature. These correlations broaden the PDF of total global sea level rise (figure 1) with largest increases in the high quantiles, for example the 95% quantile increases from 247 cm to 265 cm for RCP8.5 (table 1).

Uncertainty related to confidence in model projections
One assumption that we made is that the uncertainties associated with the variables available from the climate model ensemble are normally distributed with a mean and variance computed from the ensemble (see equation (2) and appendix equations (5) and (6). The mean and variance are then used to build the PDFs. This approach is simple as it requires no further expert judgment. However, it should be noted that in IPCC AR5 the likelihood associated with projections computed from the climate model ensemble is reassessed when formulating the likelihood range for the real world. This is done for at least three reasons. Firstly, other sources of information than the multi-model ensemble are available to constrain climate sensitivity, like paleoclimate proxies, suggesting that the range displayed by the model ensemble might be too narrow. Secondly, climate models share some common limitations like the parameterisation of clouds in the atmosphere, and mesoscale eddies and vertical mixing in the ocean, implying possible common biases. Thirdly, the climate model ensemble does 'not represent a systematically sampled family of models but rely on self-selection by the modelling groups ' Collins et al (2013). All these elements imply that the multi-model ensemble may present an incomplete picture of the range of possible futures and may provide only a limited sampling of the uncertainties (Ranger et al 2013). As a result the projected 5 to 95% range of global mean surface temperature from the climate model ensemble is considered a likely range (Collins et al 2013). In the language of the AR5, this means that it has a likelihood between 66 and 100% and not exactly 90% if one would only consider the model ensemble uncertainty. And Church et al (2013) uses the same approach for other variables originating from the CMIP5 model ensemble like thermal expansion and ocean dynamics.
To include this additional uncertainty in the results we use a parameter g to modify the standard deviation of the normal distribution computed from the climate model ensemble (Temperature and ocean dynamics and expansion). Mathematically this means that parameter r 2 from equation (2) is sampled from N ð0; g 2 Þ instead of N ð0; 1Þ. For g equal to 1 the standard deviation does not change but for g larger than 1 the uncertainty of the ensemble increases leading to a broader PDF (figure 1). We use g equal to 1.64 because this changes the likelihood of the 5%À 95% model range from 90% to 66% which is the minimum likelihood that the AR5 gives to this range. As can be expected this additional source of uncertainty has almost no impact on the median of total sea level rise (table 1). However, it is very important for the extremes. In particular it leads to an increased estimate of the 95% quantile from the AR5 PDF from 98 cm to 108 cm. And in case of a temperature dependent Antarctic mass loss the sensitivity to g is further increased as well. The 95% quantile then changes from 265 cm to 292 cm.

Conclusion and discussion
We have constructed a new high-end projection for global sea level rise in 2100 by modifying and extending the AR5 process-based method in three ways. First, we replaced the AR5 Antarctic contribution by new estimates from Deconto and Pollard (2016). Secondly, we included correlation between temperature and Antarctic mass loss. And finally, we accounted for additional uncertainty in the CMIP5 model ensemble projection of global mean surface temperature. For the RCP8.5 scenario, the PDF obtained has a median of 184 cm and a 95% quantile of 292 cm. The large shift of the median is entirely caused by replacing the AR5 median values by new high-end Antarctic estimate but a considerable part of change of the higher quantiles comes from the other two extensions. The method described here provides a systematic way to build a probabilistic projection of extreme sea level. It depends only on one modelling study of the future Antarctic ice-sheet mass loss but the results will be updated as more projections of the possible impact of hydrofracturing and ice cliff failure on future Antarctic mass loss become available. Other Environ. Res. Lett. 12 (2017) 044013 processes that are not yet parameterised in dynamical ice sheet models, like the developments of rifts from below the ice shelves (Jeong et al 2016, Alley et al 2016, might also be included in the future. We chose to keep the contribution from Greenland and from glaciers and ice caps the same as Church et al (2013) even though our focus is on extreme events. This is because uncertainties from these processes are a lot smaller than those for Antarctica and new projections are mostly in line with Church et al (2013) (Clark et al 2015, Slangen et al 2016. The choice to take a normal distribution for the Antarctic ice sheet contribution is based on the fact that Deconto and Pollard (2016) do not mention important skewness in their projections. This seems to be at odds with other projections. For instance Levermann et al (2014) obtained skewed distributions that seem to arise from the uncertainty in the subsurface ocean temperature forcing basal melt. Deconto and Pollard (2016) is exploring the high end of the distribution by including feedback processes that exclude the more probable low end of the distribution. By just focusing on the small probability high-end sea level rise it seems reasonable to assume a normal distribution as in Deconto and Pollard (2016).
As can be expected from our approach, the highend projection that we present is higher than previously published PDFs of best guess estimates (Jevrejeva et al 2014, Grinsted et al 2015, Jackson and Jevrejeva 2016. The main difference is that we use the highest published projection of mass loss from Antarctica (Deconto and Pollard 2016) instead of an expert judgment assessment (Bamber and Aspinall 2013). Other so called 'extreme' or 'worst' scenarios that are not probabilistic can be compared with our PDF. The Dutch Delta Committee (Veerman 2008) projects 110 cm that is between the 10% and 20% quantiles of our PDF. The UK H++ scenario (Ranger et al 2013) of 270 cm and the NOAA scenario of 250 cm (Kopp et al 2014, Sweet et al 2017 fall between the 80% and 95% quantiles of our PDF, so these values are not radically different from our high end estimate. The main difference is the approach taken to arrive at the high-end estimates. We suggest that our new framing for extreme sea-level scenarios, by rigorously using the process-based method, eliminates the largest sources of subjectivity. As such, it can stimulate further discussion between sea level specialists and policy makers to improve the robustness of long-term planning of adaptation measures, managing environmental risks and help making choices concerning mitigation measures.

Appendix A. Method
In this section we present each process that is expected to contribute to sea level rise in the coming century and the uncertainties associated with them. The processes are evaluated in the same way as Church et al (2013) except for the Antarctic ice sheet dynamics and its surface mass balance. The following method description builds on Vries et al (2014) and is more detailed than in Church et al (2013).
A.1. Global mean surface temperature The temperature fields are derived from 21 climate models that are part of the Coupled Model Intercomparison Project Phase 5 (CMIP5). More than 21 models participated in CMIP5 but only these models provided all the necessary variables for making the sea level projections. No other selection was performed. These 21 models are forced by two different scenarios of greenhouse gas emissions: RCP4.5 for which some mitigation measures are implemented and RCP8.5 which is business as usual.
The number of models is not large enough to determine the shape of the underlying distribution of the time varying global mean surface temperature. Therefore, we assume that this distribution is normal. We represent the global annual mean surface temperature information from all models by a matrix T, whose first dimension is time (t), and second dimension are the member of the model ensemble. A vector (N) of size n (n ¼ 7 Â 10 5 ) is built from randomly sampling the normal distribution of mean 0 and standard deviation 1 ðN ð0; 1ÞÞ. Then for each time t the temperature distribution ðT Þ is computed from the mean temperature ðT Þ and a standard deviation ðsðTÞÞ over the climate model ensemble, as: The temperature is generally used as an anomaly compared to a reference period. In this case the mean temperature during the reference period has to be removed from each model time series before computing the distribution T . This is important because the term sðT ðt; :ÞÞ also depends on the reference period. In the following a reference temperature distribution computed with the reference period 1986-2005 will be written T 1986À2005 .

A.2. Global steric expansion
Many climate models conserve volume and not mass because of the so called 'Boussinesq approximation' . Therefore, in these models an increase in temperature does not lead to a global expansion of the water. This effect is computed off-line from the density fields.
Because climate models have a drift in steric expansion it is necessary to diagnose this drift from each model Environ. Res. Lett. 12 (2017) 044013 using a control experiment that uses a constant forcing. The drift is then removed by subtracting a polynomial fit as a function of time to the control steric expansion time series. Global mean steric expansion from each model and at all time t is stored in a matrix X st . The distribution is computed in the same way as for the global mean temperature: The vector N here is the same as in equation (5). This means that the temperature and steric expansion are assumed to be completely correlated.

A.3. Land glaciers and ice caps
The contribution from land glaciers and ice caps excludes Antarctic glaciers that are included directly in the Antarctic contribution but includes Greenland glaciers. This contribution is derived from four global glacier models (Giesen and Oerlemans 2013, Marzeion et al 2012, Radićet al 2014, Slangen and Van De Wal 2011) that take into account local climate change and its effect on the surface mass balance and the hypsometry of individual glaciers. Each of these models computes the glacier contribution to sea level depending on a temperature pathway. Since these models where originally forced with different temperature pathways we first need to fit the time series of cumulated contribution to fI(t) p , with I(t) the time integral of global mean surface temperature from year 2006 to t. The integrated temperature needs to be used here because the cumulated sea level contribution depend on past temperatures. The fitting parameters f and p obtained for each model are shown in table 2. This method allows to apply these four models for any temperature pathway. In particular for the RCP scenarios: where X gic is the sea level change in cm and i is an index looping over the four sets of parameters from the glacier models. The factor 10 is used to convert from mm to cm. The sum in the second term of the right hand side of equation (8) shows that the average over the four glacier models is taken. The spread of the four models estimates around the mean is about 20%. This uncertainty is included with the vector N 2 of size n whose components are drawn from N ð0; 0:2 2 Þ. N 2 is independent from N which means that glacier modelling uncertainties are not correlated with temperature. The final distribution X gic is still partially correlated with temperature because T 1986À2005 is used to compute ℐ . An additional constant ðX 0 gic ¼ 0:95 cmÞ is added to include the change from 1996 to 2005.

A.4. Greenland ice sheet surface mass balance
The following parameterization is used for the surface mass balance tendency ð _ X Gsmb Þ in terms of global temperature change (Fettweis et al 2013): where the factor 10 À10 is used to convert GT to kg and m to cm, r w ¼ 1 Â 10 3 kg m À3 is the water density and A oc ¼ 3:6704 Â 10 4 m 2 is the ocean surface area. This equation is then integrated in time: where X 0 Gsmb is the observed contribution between 1996 and 2005. Different studies give different estimates. This uncertainty is implemented as L a vector of size n whose components are sampled from the log-normal distribution e N ð0;0:4 2 Þ . The lognormal distribution is used because the estimates of the various Greenland surface mass balance (SMB) models are positively skewed. A positive feedback between SMB and surface topography is also added. As the ice sheet looses mass its altitude decreases and the temperature at its surface increases, leading to increased melt. This is included with the vector U of size n whose components are sampled from the uniform probability distribution between 1 and 1.15.
A.5. Greenland ice sheet dynamics As in Church et al (2013) the range of the Greenland ice sheet dynamical processes contribution for 2100 is 1.4 to 6.3 cm for all scenarios, except RCP8.5 for which it is 2 to 8.5 cm. These ranges are based on an expert assessment of the literature. The mass loss rate at the beginning of the projection is taken as half of the observed rate from 2005 to 2010 (half of 0.46À0.80 mm yr À1 ), the other half being accounted for in the surface mass balance. A maximum (minimum) time series is then built starting in 2006 from the maximum (minimum) estimate of recent mass loss and ending in 2100 at the maximum (minimum) of the range for 2100 and assuming second order in time. These maximum and minimum time series are called X max Gdyn and X min Gdyn respectively. An additional 0.15 cm is added for the contribution before 2006 ðX 0 Gdyn Þ. The Table 2. Parameters for the fits to the global glacier models.

A.6. Antarctic ice sheet
This contribution is taken from Deconto and Pollard (2016). See the main text.

A.7. Groundwater changes
This term is based on projections of future dam constructions and depletion of ground water from human activities. The 5 to 95% quantiles for 2100 are À1 and 9 cm (Wada et al 2012). The time evolution is done with a second order polynomial starting from present observed rate estimates of (0.26, 0.49) [mm yr À1 ] (5%-95% range). A lower (upper) time series is constructed that start at the lower (upper) initial rate and end at the lower (upper) final estimate. These time series are called X lower grw and X upper grw . A central estimate ðX cen grw Þ is obtained as the mean of the two. The final distribution is then computed as: X grw ðtÞ ¼ X cen grw ðtÞ þ N 3 ðtÞ ð 12Þ where N 3 is sampled from N ð0; s 2 grw ðtÞÞ, with s grw ðtÞ ¼ X upper grw ðtÞ À X lower grw ðtÞ a 95 À a 05 ! ð13Þ and a q is the quantile function for a normal distribution. The groundwater contribution is taken as independent of temperature and emission scenario. Table 3. Global mean sea level PDFs in 2100 for RCP8.5 for the four scenarios from Deconto and Pollard (2016). Scenario 2 is discussed in the main text.
A.8. Final combination of processes Once all the contributions have been computed the total is obtained as: A probability density function can then be constructed from X total for each time t.