Quantifying the potential future contribution to global mean sea level from the Filchner–Ronne basin, Antarctica

Abstract. The future of the Antarctic Ice Sheet in response to climate warming is one of the largest sources of uncertainty in estimates of future changes in
global mean sea level (ΔGMSL). Mass loss is currently concentrated in regions of warm circumpolar deep water, but it is unclear how ice shelves currently surrounded by relatively cold ocean waters will respond to climatic changes in the future. Studies suggest that warm water
could flush the Filchner–Ronne (FR) ice shelf cavity during the 21st century, but the inland ice sheet response to a drastic increase in ice shelf melt rates is poorly known. Here, we use an ice flow model and uncertainty quantification approach to project the GMSL contribution of the FR basin under RCP emissions scenarios, and we assess the forward propagation and proportional contribution of uncertainties in model parameters (related to ice dynamics and atmospheric/oceanic forcing) on these projections. Our probabilistic projections, derived from an extensive sample of the parameter space using a surrogate model, reveal that the FR basin is unlikely to contribute positively to sea level rise by the 23rd century. This is primarily due to the mitigating effect of increased accumulation with warming, which is capable of suppressing ice loss associated with
ocean-driven increases in sub-shelf melt. Mass gain (negative ΔGMSL) from the FR basin increases with warming, but uncertainties in these projections also become larger. In the highest emission scenario RCP8.5, ΔGMSL is likely to range from −103 to 26 mm, and this large spread can be apportioned predominantly to uncertainties in parameters driving increases in precipitation (30 %) and sub-shelf melting (44 %). There is potential, within the bounds of our input parameter space, for major collapse and retreat of ice streams feeding the FR ice shelf, and a substantial positive contribution to GMSL (up to approx. 300 mm), but we consider such a scenario to be very unlikely. Adopting uncertainty quantification techniques in future studies will help to provide robust estimates of potential sea level rise and further identify target areas for constraining projections.



Introduction
Ice loss from the Antarctic Ice Sheet has accelerated in recent decades Shepherd et al., 2018) and the 20 future of the ice sheet under climate warming is one of the largest sources of uncertainty for global mean sea level rise. Current projections suggest that the ice sheet may contribute anywhere between −7.8 and 30 cm to sea level rise by 2100 under Representative Concentration Pathway (RCP) 8.5 scenario forcing . This large spread of potential sea initialised to observations using a model inversion with m = 3 and n = 3, over the grounded portion of the catchment. Blue to yellow shading shows sub-shelf melt rates across the Filchner and Ronne ice shelves, using the ocean box melt parameterisation with point estimates for parameters (see Appendix B). Light to dark blue shading shows sea floor depth from the IBCSO dataset (Arndt et al., 2013). Inset map shows the full extent of our model domain (red).

Figure 2.
Workflow diagram summarising the uncertainty quantification approach used in this study. We first identify uncertain input parameters and represent them in probabilistic framework. A training sample of 500 points is taken from this input parameter space and used as input to an ensemble of simulations in our ice flow model. Using this training sample, and the surrogate modelling capabilities in UQLAB we create a polynomial chaos expansion (PCE) that mimics the behaviour of our ice flow model. This allows us to evaluate a much larger sample from our parameter space and these surrogate models are used to derive predictions and probability density functions for changes in global mean sea level (∆GMSL). Finally, we use sensitivity analysis to identify the proportional contribution of each input parameter on projection uncertainty. warming, on a quantity of interest. We make use of the MATLAB based toolbox, UQLab, and the uncertainty quantification framework of Sudret (2007), on which the MATLAB based toolbox is based (Marelli and Sudret, 2014). UQLab includes an extensive suite of tools encompassing all necessary aspects of uncertainty quantification. Here, we summarise the approach and tools used in this study ( Figure 2) and we refer the reader to the UQLAB documentation (uqlab.com Marelli and Sudret,95 2014) for further details.
We can think of a physical model (M) as a map from an input parameter space to an output quantity of interest, as where our uncertain input parameters are specified as a probabilistic input model (X) with a joint probability distribution function X ∼ f X (x), and Y is a list of model responses. Using this approach we are able to propagate the uncertainties in the 100 inputs X to the outputs Y . We can think of our ice-flow model in the same way, ∆GMSL =Úa(X), where ∆GMSL is our model response or quantity of interest. In the following sections we outline eight uncertain input parameters that are represented in X. These relate to basal sliding and ice rheology (Section 3.2), surface accumulation (Section 3.3) and sub-shelf melting (Section 3.4). Uncertainties in these input parameters are defined in a probabilistic way based on the available information ( Figure 3). For parameters used to force sub-shelf melt rates, we conducted a separate Bayesian analysis to determine their input parameter probability distributions (see Appendix B).
Quantifying the uncertainty in model outputs due to uncertainty in input parameters or forcings, may require a computationally unfeasibly large number of model evaluations. However if, for example, the model response varies slowly as the values of some input parameters are changed, the relationship between model inputs and model outputs may be approximated using a much simpler and computationally faster surrogate model. The uncertainty estimation can then be done in a much more where Ψ α (X) are multivariate polynomials that are orthonormal with respect to the join input probability density function f X , A ⊂ N M is a set of multi-indices of the multivariate polynomials Ψ α , and y α are the coefficients. Here, our PCEs are calculated using the least angle regression (LAR) algorithm in UQLab (Blatman and Sudret, 2011;Marelli and Sudret, 2019) that solves a least-square minimisation problem. This algorithm iteratively moves regressors from a candidate set to an active 120 set and at each iteration a leave-one-out cross-validation error is calculated. After all iterations are complete, the best sparse candidate basis are those with the lowest leave-one-out error. This is designed to reduce the potential for over-fitting, and reduced accuracy when making predictions outside of the training set. This sparse PCE calculation in UQLab also uses the LOO error for: 1) adaptive calculation of the best polynomial degree based on the experimental design and 2) adaptive q-norm setup for truncation scheme. For further details on the PCE algorithm see Marelli and Sudret (2019). We also outline details on 125 how input uncertainties were propagated through our model to create our PCE in Section 3.5.
Once the surrogate model has been created, the moments of the PCE are encoded in its coefficients where the mean (µ P C ) and variance (σ P C ) 2 are as follows Our existing PCE surrogate models can additionally be used in a sensitivity analysis to quantify the proportional contribution of parametric uncertainty on projections of ∆GMSL. This allows us to identify input parameters where improved understanding is needed to constrain future projections. Here, we are using Sobol indices which are a variance-based method where the model can be expanded into summands of increasing dimension, and total variance in model output can be described as the sum of 135 the variances of these summands.
First order indices (S i ), often also referred to as "main-effect", are the individual effect of each input parameter (X i ) on the variability in the model response (Y ), defined as: Total Sobol indices (S T i ) are then the sum of all Sobol indices for each input parameter and encompass the effects of parameter 140 interactions. Values for Sobol Indices are between 0 and 1, where large values of S i indicate parameters that strongly influence the projections of global mean sea level. If S i ≈ S T i then it can be assumed that the effect of parameter interactions is negligible. These Sobol indices can be calculated analytically from our existing PCEs, by expanding portions of the polynomial that depend on each input variable to directly calculate parameter variance using the PCE coefficients. Each of the summands of the PCE can be expressed as Due to the orthonoamlity of the basis, the variance of our truncated PCE reads as The first order Sobol indices in Equation 5 are then calculated as the ratio between the two above terms.

Ice-flow model
Here we use the vertically integrated ice-flow model Úa  to solve the ice dynamics equations using the shallow-ice stream approximation (SSTREAM), also commonly referred to as the shallow-shelf approximation (SSA) and the 'shelfy-stream' approximation. (MacAyeal, 1989). Úa has been used in previous studies on grounding line migration and ice 155 shelf buttressing and collapse (De Rydt et al., 2015;Reese et al., 2018b;Gudmundsson et al., 2012;Gudmundsson, 2013;Hill et al., 2018) and model results have been submitted to a number of intercomparison experiments (Pattyn et al., 2008(Pattyn et al., , 2012Levermann et al., 2020;Cornford et al., 2020).
Our model domain extends across the two major drainage basins that feed into the FRIS (Figure 1). Within this domain, we generated a finite-element mesh with ∼ 92, 000 nodes and ∼ 185, 000 linear elements using the Mesh2D Delaunay-based 160 unstructured mesh generator (Engwirda, 2015). Element sizes were refined based on effective strain rates and distance of the grounding line, ranging from 900 m close to the grounding line, and in the shear margins, to a maximum of 50 km further inland. Outside of our uncertainty analysis, we tested the sensitivity of our results to mesh resolution by repeating our median and maximum ∆GMSL simulations under RCP 8.5 forcing, and dividing or multiplying the aforementioned element sizes by two. Our results are largely insensitive to the mesh resolution, with a percentage deviation of only 3%. Finally, we linearly 165 interpolated ice surface, thickness and bed topography from BedMachine Antarctica v1 (Morlighem et al., 2020) onto our model mesh. We initialise our model to match observed velocities using an inverse approach (see Section 3.2 and Appendix A).
During forward transient simulations, Úa allows for fully implicit time integration, and the non-linear system is solved using the Newton-Raphson method. Úa includes automated time-dependent mesh refinement, allowing for high mesh resolution 170 around the grounding line as it migrates inland. We also impose a minimum thickness constraint using the active-set method to ensure that ice thicknesses remain positive. Throughout all simulations our calving front remains fixed in its originally prescribed position. At the end of each forward simulation we calculate the final change in global mean sea level (∆GMSL) as the ice volume above flotation that will contribute to sea level change based on the area of the ocean .

175
There are two components of surface glacier velocities; internal deformation and basal sliding. Úa uses inverse methods to optimise these velocities components based on observations by estimating the ice rate factor (A) in Glen's flow law and basal slipperiness parameter (C) in the sliding law. This section introduces uncertainties related to the exponents of the flow law and basal sliding law, whereas details of the inverse methodology are included in Appendix A.
Glen's flow law (Glen, 1955) is used to relate strain rates and stresses as a simple power relation where˙ ij are the elements of the strain rate tensor, τ e is effective stress (second invariant of the deviatoric stress tensor), τ ij are the elements of the deviatoric stress tensor, A is the temperature dependent rate factor, and n is the stress exponent.
This stress exponent (n) controls the degree of non-linearity of the flow law and most ice flow modelling studies adopt n = 3, as it is considered applicable to a number of regimes (See review in: Cuffey and Paterson, 2010). However, experiments 185 reaching high-stresses (Kirby et al., 1987;Goldsby and Kohlstedt, 2001;Treverrow et al., 2012), or analysing borehole measruements, and ice velocities (e.g. Gillet-Chaulet et al., 2011;Cuffey and Kavanaugh, 2011;Bons et al., 2018) have suggested that n > 3. It is also possible that at low stresses, the creep regime may become more linear n < 3 (Jacka, 1984;Pettit and Waddington, 2003;Pettit et al., 2011), which is supported by ice shelf spreading rates n = 2 − 3 (Jezek et al., 1985;Thomas, 1973). While it can be considered that n = 3 is appropriate in most dynamical studies, the exact numerical value is not known and it appears plausible that it can range between 2 and 4. To capture the uncertainty in the stress exponent we take n ∈ [2, 4] and sample continuously from a uniform distribution within this range ( Figure 3).
Basal sliding is considered the dominant component of surface velocities in fast flowing ice streams. The Weertman sliding law is defined as where C is a basal slipperiness coefficient and v b the basal sliding velocity. The Weertman sliding law typically captures hardbed sliding, in which case m = n and is normally set equal to 3 (Cuffey and Paterson, 2010). However, using different values for m alters the non-linearity of the sliding law, and can thus be used to capture different sliding processes, i.e. viscous flow for m = 1 and plastic deformation for m = ∞. There are limited in situ observations of basal conditions, and the value of m relies on numerical estimates of basal sliding based on model fitting to observations.

200
A number of studies have tested different values of m to fit observations of grounding line retreat or speedup at Pine Island Glacier (Gillet-Chaulet et al., 2016;Joughin et al., 2010;De Rydt et al., 2021). These studies show that m = 3 can underestimate observations, and more plastic like sliding (m > 3) is needed in at least some parts of the catchment to replicate observations (Joughin et al., 2010;De Rydt et al., 2021). This uncertainty in the value of m can ultimately affect projections of sea level rise (Ritz et al., 2015;Bulthuis et al., 2019;Alevropoulos-Borrill et al., 2020) by altering the length and time taken 205 for perturbations (e.g. ice shelf thinning or grounding line retreat) to propagate inland.
While additional sliding laws have been proposed and are now implemented within a number of existing ice flow models, in this study we use the Weertman sliding law, as it remains the most common. This narrows the parameter space, allowing us to fully integrate the influence of m on projections of future sea level rise into our uncertainty assessment (by performing an inverse model run prior to each perturbed run, see Section 3.5). This is an advancement over previous Antarctic wide 210 studies, that given domain size have no choice but to invert the model for a handful of different m values prior to uncertainty propagation (e.g. Bulthuis et al., 2019;Ritz et al., 2015). To capture uncertainty in m and to sample from a range of possible methods of basal slip, we take m ∈ [2, 9] and sample from a uniform distribution ( Figure 3).

Surface accumulation
To capture uncertainties in future climate forcing, we use projections from four Representative Concentration Pathways (RCPs) 215 presented in the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. These pathways capture plausible changes in anthropogenic greenhouse gas emissions for the 21st and 22nd centuries. RCP 2.6 is a strongly mitigated scenario with a global temperature increase of less than 2 • C above pre-industrial levels by 2100, and is the goal of the 2016 Paris Agreement. Two intermediate scenarios (RCP 4.5 and RCP 6.0) represent global temperature increases of ∼ 2.5 • C and ∼ 3 • C with reductions in emissions after 2040 and 2080 respectively. Finally, RCP 8.5 projects a global temperature increase of 220 ∼ 5 • C by 2100 and is now often referred to as an "extreme" or "worst-case" climate change scenario. Global mean temperature changes (∆T g ) from 1900 to 2300 relative to pre-industrial (1860-99) were obtained from the atmosphere-ocean general circulation model emulator MAGICC6.0 (live.magicc.org: Meinshausen et al. (2011)). For each RCP scenario we obtain 600 (historically-constrained) model simulations between 2000 and 2100 (see Meinshausen et al. (2009) for details on the probabilistic set-up). We then use the ensemble median and uncertainty bounds within a "very likely" 225 range between the 25th and 75th percentiles. To extend the record to 2300, we use a single model realisation, using the default climate parameter settings used to produce the RCP greenhouse gas concentrations for each RCP scenario (Meinshausen et al., 2009) , and keep the upper and lower bounds constant from 2100 to 2300 ( Figure 4). Uncertainty in projections from 2100 to 2300 may well be larger, but we choose not to make an assumption on how errors will propagate up to 2300. Global temperatures from MAGICC 6.0 were also used in the Antarctic linear response model inter-comparison (LARMIP-2) experiment 230  and are consistent with projections used in other Antarctic wide simulations (Bulthuis et al., 2019;Golledge et al., 2015).
Following the work of a number of previous studies (e.g. Pattyn, 2017;Bulthuis et al., 2019;DeConto and Pollard, 2016;Garbe et al., 2020), global temperature changes (∆T g ) are used to force annual changes in surface mass balance through our forward-in-time simulations, by prescribing changes in surface temperature (T air ) and precipitation (P ) as follows:
where T air obs and A obs are surface temperatures and accumulation rates from RACMO2.3 respectively (Van Wessem et al., 2014). Temperature changes through time are corrected for changes in surface elevation (s) from initial observations (s obs ), 240 using a lapse rate of γ = 0.008 • C m −1 (Pattyn, 2017;DeConto and Pollard, 2016), and subsequently used to force changes in precipitation using an expected percentage increase in precipitation (p) per degree of warming (Aschwanden et al., 2019). This captures the rise in snowfall expected with the increased moisture content of warmer air, suggested by climate models (e.g. Palerme et al., 2017;Frieler et al., 2015).
To capture further uncertainties associated with atmospheric forcing, we introduce two uncertain parameters into our anal-245 ysis, 1) a scaling factor to select a temperature realisation between the 25th and 75th percentiles of the ensemble median temperature for each RCP scenario ( Figure 4) and 2) uncertainty in the expected changes in precipitation across the Antarctic ice sheet with increased air temperatures.
First, instead of only using the ensemble median change in temperature for each RCP scenario, we capture the spread within each forcing scenario by incorporating a temperature scaling parameter (T ) as follows: ∆T g (n) = ∆T median g +T ·T err g ,

250
where for each RCP scenario ∆T median g is the median, ∆T err g is the distance either side of the median within the 25th and 75th percrentiles, and ∆T g (n) is the resultant temperature realisation used to force both surface accumulation (P ) and ocean temperature (see Section 3.4). We assume that there is decreasing likelihood of temperature profiles further away from the median, and so prescribe a Gaussian distribution for T between −1 (25th percentile) and 1 (75th percentile) and centered around 0 (median: Figure 3).

255
Secondly, we capture uncertainty associated with precipitation by varying the amount by which precipitation increases per degree of warming (p). While it is generally accepted that accumulation will increase with future warming, the value of p remains uncertain. Snow accumulation could prevent runaway ice discharge from the Antarctic ice sheet, which means that parameterisations of precipitation increase with warming have implications for accurate projections of mass change across the ice sheet. Previous studies using ice core records, historical global climate model (GCM) simulations, and future GCM 260 simulations as part of the CMIP5 ensemble, have estimated anywhere between 3.7-9% increase in Antarctic accumulation per degree of warming (Krinner et al., 2007(Krinner et al., , 2014Gregory and Huybrechts, 2006;Bengtsson et al., 2011;Ligtenberg et al., 2013;Frieler et al., 2015;Palerme et al., 2017;Monaghan et al., 2008). To capture this range of possible values for (p) we sample from p ∈ [4, 8] and make no assumption of the distribution (likelihood) of the value of p within this range by sampling from a uniform distribution (Figure 3).

Sub-shelf melt
Ice shelf thinning due to ocean-induced melt can reduce buttressing forces on grounded ice and accelerate ice discharge.
Such feedbacks may already be taking place in parts of West Antarctica. However, future changes in ocean conditions remain uncertain, owing to poor understanding and the challenges of modelling interactions between global atmospheric warming and ocean circulation/temperature changes (Nakayama et al., 2019;Thoma et al., 2008). In particular, the likelihood that the 270 Filchner-Ronne ice shelf cavity will be flushed with modified warm deep water in the future is unclear (Hellmer et al., 2012).
To parameterise basal melting beneath the ice shelf, we use an implementation of the PICO ocean box model (Reese et al., 2018a) for use in Úa, which we hereafter refer to as the ocean box model. This provides a computationally feasible alternative to fully coupled ice-ocean simulations for large ensemble analysis, but more physically based than simple depth dependent parameterisations (e.g. Favier et al., 2014). The basic overturning circulation in ice shelf cavities is captured using a series of 275 ocean boxes, calculated based on their distance from the grounding line. The overturning flux q is then calculated as the density difference between the far-field (p 0 ) and grounding line (p 1 ) water masses using a constant overturning strength parameter (c).
The melt parameterisation also includes a turbulent heat exchange coefficient γ * T that controls the strength of melt rates by varying the heat flux across the ice-ocean boundary. For a detailed description of the physics of the PICO box model, see Reese et al. (2018a). To calculate sub-shelf melt rates, the box model requires inputs of temperature (T ocean ) and salinity (S) on the 280 continental shelf to drive the ocean cavity circulation. We use S = 34.65 psu and the initial observed ocean temperature for the Weddell Sea T ocean obs = −1.76 • C from Schmidtko et al. (2014) which was proposed for use in PICO (Reese et al., 2018a). For the FR basin we use five ocean boxes, and only apply sub-shelf melting to nodes that are fully afloat (no connecting grounded nodes) to avoid overestimating grounding line retreat .
with an ocean temperature anomaly: It is often assumed that atmospheric temperature changes ∆T g can be translated to ocean temperature changes ∆T o using some scaling factor (α) (Maris et al., 2014;Golledge et al., 2015;Levermann et al., 2014Levermann et al., , 2020. Here, we use the linear scaling proposed in Levermann et al. (2020) which additionally includes a time delay τ to capture the assumed time lag between 290 atmospheric and subsurface ocean warming.
To obtain suitable values for α and τ , Levermann et al. (2020)  We identify a further two uncertain parameters in the ocean box model that additionally control the strength of sub-shelf 300 melt. These are the turbulent heat exchange coefficient γ * T and the strength of the overturning circulation c. Values for these parameters presented in Reese et al. (2018a) have been optimised to present day ocean temperatures and observations of melt rates for a circum-Antarctic set-up. While upper and lower bounds for these parameters have also been proposed (Reese et al., 2018a;Olbers and Hellmer, 2010), little information exists on the likelihood of parameter values within these ranges, particularly for different regions of Antarctica, with varying ocean conditions.

305
For the four parameters that control the sub-shelf melt rates (α , τ ,γ * T and c) we decided to constrain their uncertainty (probability distributions) using a Bayesian approach. Using the a priori information on the distributions for α and τ

Propagating uncertainty
In this section we explain how uncertainties in the input parameters introduced in the previous sections are propagated through our model to obtain projections of global mean sea level ( Figure 2). We began by generating an experimental design (training set) for the surrogate model. An input parameter sample of 500 points was extracted from the parameter space using Latin of model drift at the start of the simulation, characterised by a slowdown and thickening of many of the ice streams in our domain. We found model drift to be similar between parameter sets, but was affected by the basal boundary conditions from our inversion. Hence, rather than specify a single baseline for the entire experimental design, we perform a control simulation for each set of basal boundary conditions. This control run uses selected values of m and n and inverted fields of C and A but holds all other input parameters fixed to their point estimates, as well as using a constant temperature forcing. Each control run 325 is followed by four forward runs, one for each RCP forcing scenario, in which surface accumulation and sub-shelf melt rates are updated at annual intervals based on global temperature changes ( Figure 4). The final calculated change in global mean sea level for each RCP scenario is with respect to the preceding control run (∆GMSL rcp − ∆GMSL ctrl ). For our 500 member ensemble we perform 500 model inversions, and 2500 forward simulations.
Model responses (∆GMSL) and input parameter sample X = x i , ..., x N are used to train four surrogate models (one for 330 each RCP scenario). This is done using the Polynomial Chaos module in UQLab  using the LAR algorithm previously described in Section 2. We allow the LAR algorithm to choose a PCE with a degree anywhere between 3 and 15 and q-norm between 0.1 and 1. Predictions (mean and variance) of ∆GMSL are then directly extracted from each surrogate model. To estimate the accuracy of our PCE predictions and calculate quantiles we use bootstrap replications. We use 1000 replications (B) each with the same number of sample points as the original experimental design (500) to create an 335 additional set of B PCEs and associated responses. Quantiles (5th and 95th) were extracted from the bootstrap evaluations (see Supplementary Figure S1). To assess the performance of our surrogate model, we generated an additional and independent validation sample of 20 points in the parameter space, evaluated for each RCP scenario (total of 80 perturbed simulations).
Predictions made by the surrogate models are close to the responses generated by our ice-flow model for all RCP scenarios (Supplementary Figure S1).

Projections of sea level rise from FR basin
We begin by presenting probabilistic projections of global mean sea level change from the Filchner-Ronne basin for four RCP scenarios. These projections were derived from surrogate models that were trained with a 500 member ensemble of forwardin-time ice flow model simulations. We then evaluated these surrogate models with a 1,000,000 point sample from our input 345 parameter space to derive model responses, and calculated probability distributions using kernel density ( Figure 5).
Our projections indicate it is most likely that the FR basin will undergo limited change or contribute negatively to global mean sea level by the year 2300. Under the lowest warming scenario (RCP 2.6: 0.77 − 1.7 • C), ∆GMSL is limited, ranging  Table 1. Contribution to global mean sea level [mm] at years 2100, 2200 and 2300. The first number is the median projection and values in brackets are 5-95% confidence intervals between −19.1 mm (5th) and 9.57 mm (95th percentiles: Table 1). The probability distribution ( Figure 5) takes a near-tonormal shape, with a median projection close to zero (−5.49 mm), but has a weak positive skew of 0.16 (calculated using the 350 moment coefficient of skewness), with a tail extending towards a maximum sea level contribution of ∼ 50 mm. Projections of ∆GMSL under the medium warming scenarios RCP 4.5 and 6.0 range from −36.1 to 9.7 mm and −49.6 to 11.1 mm respectively (Table 1). These distributions are more positively skewed than RCP 2.6 (skewness coefficients of 0.29 and 0.39 respectively), with tails extending towards ∼ 100 mm of global sea level rise ( Figure 5). Extreme warming leads to the greatest uncertainty in projections, which range from −103 to 26 mm under RCP 8.5 (Table 1). The median projection indicates 355 a greater negative contribution to sea level rise under higher warming. However, the probability distribution is asymmetric, with a long-tail (high positive skew = 0.77), that decreases exponentially away from the median and reaches a maximum ∆GMSL of 332 mm ( Figure 5). This long-tail represents the potential for low-probability, but high magnitude contributions to sea level rise under certain parameter combinations.
As our surrogate modelling is based around a single quantity of interest (∆GMSL) it does not allow us to evaluate temporal 360 changes in ice loss directly. Figure 6 instead presents projections through time from our 500 member ensemble alongside the final ∆GMSL from each surrogate model (PCE). We also generate two additional surrogate models for each RCP scenario at years 2100 and 2200 (Table 1) to evaluate projections at these time intervals, and identify the temporal importance of parameters on projection uncertainty (see Section 4.2 and Figure 7).
Both Table 1 and Figure 6 show that the contribution to ∆GMSL and associated uncertainties increase through time. Within 365 the next 100 years (up to 2100) we project little change in ice mass from the FR basin. This constitutes a small negative contribution to sea level rise of < 10 mm in all warming scenarios (Table 1), with a maximum range of −14.1 to 0.11 mm.
By 2200 the spread of ∆GMSL has diverged based on warming scenario, with little change under limited forcing (RCP 2.6 = −5.05), and a greater negative contribution under with higher warming (−30.2 in RCP 8.5). Between 2200 and 2300 uncertainties increase dramatically in all warming scenarios, particularly in RCP 8.5. Box plots on Figure 6 show the projections 370 generated from our surrogate model in 2300 alongside our ice flow model ensemble. We note that the tail for RCP 8.5 extends substantially beyond the maximum ∆GMSL from our ensemble (150 mm). While, this is expected with more extensive sampling of our parameter space, we test the feasibility of ∆GMSL = 332 mm, by taking the parameter values that led to this, and re-evaluating the "true" ice flow model. This gives a slightly lower value of ∆GMSL = 250 mm, but one that is still considerably higher than in our original ensemble despite its relatively large size (N = 500), highlighting the benefits of our Dark shading is the interquartile range (IQR) defined between the 25th and 75th percentiles. Lighter shading shows 5th-95th percentiles.
Box plots show the projections from the surrogate models (∆GMSLPCE) for each RCP scenario at 2300. Extreme values are located at 1.5 times the interquartile range away from the 25th and 75th percentiles. Values outside of these extreme bounds are considered to be outliers. surrogate modelling approach. Recalculating the surrogate model for RCP 8.5 including this "extreme" sample point, reduces the maximum contribution to sea level rise to 288 mm.

Parametric Uncertainty
In this section we present the results of our sensitivity analysis, in which we determine how uncertainties in our input parameters (parametric uncertainty) impact our projections of ∆GMSL. To do this, first order Sobol indices were decomposed from each 380 of our PCE models (four RCP forcing scenarios) and for three timesteps: 2100, 2200 and 2300, which are presented in Figure   7. We additionally assessed the individual parameter to projection relationship, by re-evaluating our surrogate model for each parameter, while all other parameters were held at their point estimates (see Supplementary Figure S2).
By 2300 (dark shaded bars in Figure 7) uncertainties in our four ocean forcing parameters collectively have the greatest fractional contribution to the uncertainty in our projections of global mean sea level contribution. This ranges from 60% in 385 RCP 8.5 to 75% for RCP 2.6. Projection uncertainty in all RCP scenarios is primarily driven by ocean temperature forcing, and the value of α used to scale atmospheric to ocean temperatures. Uncertainties attributed to α appear to increase both with warming scenario and through time. In all RCP scenarios, fractional uncertainty associated with α increases from 2100 to 2300 (light shaded bars in Figure 7), coincident with an increase in the spread of ∆GMSL contribution ( Figure 6). In 2100, α has a greater impact on projection uncertainty in the lower warming scenario. However, by 2300, the fractional uncertainty is 390 greatest in RCP 8.5, accounting for almost half of projection uncertainty (0.44) compared to 0.34 in RCP 2.6. Re-evaluating the surrogate models varying only the value of α reveals a quadratic dependency of ∆GMSL on the value of the scaling coefficient ( Supplementary Fig S2), which is consistent with the quadratic sensitivity of sub-shelf melt rates to ocean temperature forcing observed for the FR ice shelf cavity by Reese et al. (2018a). Under extreme warming (RCP 8.5), this quadratic relation becomes stronger, and variability in α alone can cause ∆GMSL to range between −86 and 73 mm by 2300. Under all RCP warming 395 scenarios the value of α contributes to a greater range of ∆GMSL than any other parameter (Supplementary Fig S2), and encompasses almost all of the 5-95% spread of projections (Table 1).
Of our two ocean box model parameters, overturning strength (c) accounts for more projection uncertainty than the turbulent heat exchange coefficient (γ * T ) in all warming scenarios. This is consistent with the theory that sub-shelf melting at large and cold-cavity ice shelves is predominantly driven by overturning strength (Reese et al., 2018a). The fractional importance of c has the greatest variability between forcing scenarios than any other parameter. Unlike α, uncertainty associated with c decreases with warming, from 0.32 in RCP 2.6 to 0.1 for RCP 8.5 (Figure 7). The importance of the overturning strength, c, also increases with time, which is most pronounced in lower warming scenarios, e.g. RCP 2.6 where α and c are similar by 2300. This suggests that with greater ocean warming (in RCP 8.5) and a transition to warm cavity conditions, uncertainties in temperature (associated with the value of α) outweigh uncertainties in sub-shelf melt rates driven by the overturning strength 405 alone. Whereas in colder conditions (RCP 2.6) variability in c has a greater control on heat supply for sub-shelf melt. A similar trend exists for uncertainties associated with γ * T ; greater importance for lower warming scenarios and in all scenarios, increasing importance with time. However, in contrast to c, there is a greater relative increase in the Sobol index for the highest warming scenario (RCP8.5) from 2100 to 2300; the importance of γ * T doubled from 0.033 to 0.066 versus only a 54% increase in RCP 2.6. This suggests that as the FR ice shelf transitions to warm cavity conditions (∼ 2 • C in RCP 8.5) the heat exchange 410 in the turbulent boundary layer may become a more important driver of sub-shelf melt than under colder conditions. Atmospheric forcing parameters account for the second largest proportion of uncertainty in ∆GMSL by 2300. This is primarily driven by variability in the percentage increase in precipitation per degree of warming (p), and to a lesser extent the temperature scaling parameter (T ). At all time intervals (2100, 2200, 2300), projection uncertainty attributed to p is largest based on warming scenario. Unlike α, Sobol indices for p decrease through time, which is asynchronous to increased uncer-415 tainty in ∆GMSL contribution ( Figure 6). In 2100, p accounts for over half of projection uncertainty in all scenarios (except RCP 2.6) reaching a maximum of 0.62 in RCP 8.5. p remains the dominant parameter in 2200 for higher warming scenarios, but for RCP 2.6 the Sobol index for p decreases to 0.26, less than both α and c. By 2300, the fractional importance of p is < 0.3 and lower than α for all four warming scenarios. Evaluating the surrogate model (at year 2300) for p only, reveals a linear dependency on the value of p, where, as expected, increases in precipitation lead to a decrease in the contribution to 420 GMSL, or in this case a greater negative contribution to GMSL (Supplementary Figure S2). In RCP 8.5, p alone contributes to between −8 and −80 mm of ∆GMSL (Supplementary Figure S2). This suggests that even with a limited (p = 4%) increase in precipitation, and fixed melt rates, the FR basin is unlikely to contribute positively to ∆GMSL. However, in RCP 2.6, increased accumulation with p < 0.05 does not outweigh mass loss associated with sub-shelf melting and could lead to a small positive contribution to sea level rise.

425
Finally, uncertainties in our ice dynamical parameters relating to the non-linearity in the sliding (m) and flow (n) laws used in our model, have a limited contribution to uncertainties in our projections of ∆GMSL (Figure 7). The combined contribution of these parameters by 2300 under all warming scenarios (0.02 in RCP8.5) is an order of magnitude less than uncertainties associated with atmospheric and oceanic forcing (see Supplementary Figure S3 for Sobol Indices for just m and n). Of the two parameters, m accounts for the most uncertainty in ∆GMSL, which is unsurprising given that basal sliding is likely to be the  Figure S3). Increasing the values of m and n in isolation, reduces the negative contribution to ∆GMSL, i.e. less mass gain (Supplementary Fig S2). In both cases, a stronger non-linearity in the ice flow (n), or more plastic like flow (m), allows for faster delivery of the ice to the grounding line in response to a perturbation.

Partitioned Mass Change
Our Sobol indices reveal that the percentage change in precipitation and ocean temperature scaling are the main drivers of uncertainty in changes in global mean sea level. To further examine the relative importance of precipitation and sub-shelf melt parameters on mass change in the FR basin, we partition components of mass balance (accumulation and discharge) using the input-output method. We calculated the integrated input accumulation (P ) across the grounded area and the total integrated 440 discharge (D) output across the grounding line with respect to our control runs. These mass balance components, as well as total mass change (M = P − D), are shown for the low (RCP 2.6) and high (RCP 8.5) warming scenarios in Figure 8) and for intermediary scenarios in Supplementary Figure S4.
Mass change under the lowest warming scenario (RCP 2.6) closely follows the temperature anomaly trend and appears primarily driven by increases (and subsequent decreases) in accumulation with warming. In the first 50 years, mass balance 445 increases to 12.6 (6.56 − 18.6) Gt yr −1 (where values in brackets here and in the remainder of this section are 5-95%). This is primarily due to an increase in accumulation at a rate of 15.2 Gt yr −1 in 2050, which is offset by a limited increase in discharge across the grounding line (2.5 Gt yr −1 ) during this period. Between 2050 and 2100 accumulation remains constant and discharge increases, which consequently reduces the rate of total mass gain. Uncertainties associated with accumulation are greater than those for discharge during this period (Figure 8c), which is consistent with the high contribution of the percentage 450 increase in precipitation on projection uncertainty in 2100 (Figure 7). The rate of mass gain continues to decrease after 2100, alongside a reduction in the temperature perturbation in RCP 2.6, and decelerating accumulation. During this period, discharge across the grounding line stabilises at a median of 10 Gt yr −1 , but the uncertainty range increases dramatically to −8 to 26 Gt yr −1 , which coincides with an increasing importance of parameters relating to sub-shelf melt from 2200 to 2300 ( Figure   7). By 2300 this drives the mass balance towards zero, at which accumulation is approximately balanced by ice discharge.

455
Under RCP 8.5 forcing the spread of mass change in 2300 is driven by anomalies in ice discharge. During the first 150 years, surface accumulation steadily increases at an average rate of 50 Gt yr −1 , which is consistent with increased temperature forcing of 6.4 • C. During this period, increases in discharge lag that of accumulation, averaging only 9.6 Gt yr −1 , which can partly be explained by the prescribed time delay between atmospheric and oceanic warming (τ : Equation 14). Hence, it appears likely that total mass balance will remain positive up to 2150, as no parameter combinations (ensemble members) are 460 able to sufficiently increase ice discharge above that of accumulation. Consistent with other forcing scenarios, uncertainties in accumulation are also greater than those associated with discharge, which corresponds to greater projection uncertainty attributed to p up to 2200 (Figure 7).
After 2150, the rate of temperature increase is reduced, which leads to a reduction in surface accumulation, that plateaus at ∼ 140 (86 − 202) Gt yr −1 between 2200 and 2300. Despite a limited change in temperature (+1.6 • C) between 2150 and 465 2300 (relative to 2000 to 2150), discharge continues to increase linearly from 31 to 92 Gt yr −1 . Simultaneously, uncertainties in discharge increase substantially and span 0.05 − 0.76 mm yr −1 (5-95th percentiles) of sea level equivalent volume in 2300.
This suggests that the atmospheric temperature anomaly itself becomes less important than the amount by which atmospheric temperatures are scaled to ocean warming, i.e. the value of α, where variability in α alone accounts for most of the range of sea  Figure S2). Indeed, the spread of total mass change (M ) closely follows the 470 uncertainties in ice discharge, where it is possible, albeit unlikely, that certain combinations of parameters (within the 5-95th percentiles) allow for mass imbalance (P < D) and a positive contribution to sea level.
Variability in ice discharge alone also reveals the spread of potential positive contribution to sea level rise that would have occurred in our simulations if surface accumulation had remained unchanged. To explore this, we rerun our median simulation under RCP 8.5 forcing using the same parameter values, but keep surface mass balance fixed at it's initial value. This reveals 475 that increases in sub-shelf melt alone would contribute to 84 mm of global mean sea level rise (as opposed to −50 mm), and highlights the important and compensating effect accumulation has on the sign of our sea level projections.

Grounded Ice Loss
In this section we explore the changes in grounded area throughout our simulations to see how our projections of mass change correspond to the retreat of the grounding line. Figure 9 presents changes in grounded area with respect to the control runs (a) Figure 9. Changes in grounded area and grounding line position for RCP 8.5. Top panel shows change in grounded area (∆GA calculated as GArcp − GA ctrl ) in ×10 4 km 2 . Coloured lines represent the 5th, 50th and 95th percentiles of the projections of ∆GMSL rather than the percentiles of the change in grounded area itself. However, they are generally close to the grounded area results. Lower panel shows the FR basin and bed elevation in metres above (green to brown) and below (light to dark blue) sea level. Coloured lines show grounding line positions from our ensemble simulations that lie closest to our percentiles (5, 50 and 95%) from our surrogate model projections, with respect to control runs (dashed grey lines). Two additional grounding line positions are shown; the maximum ∆GMSL from our ensemble (orange) and the maximum ∆GMSL from our surrogate model (magenta), which we evaluated separately to our initial ensemble of simulations. and grounding line positions from ensemble members closest to our 5, 50 and 95% percentile projections of ∆GMSL from our surrogate models (b) for RCP 8.5, while additional RCP scenarios are shown in Supplementary Figures S5. S7.
Despite the negative contribution to ∆GMSL likely under all warming scenarios (50th percentiles: Table 1), our results show that these median projections correspond to simulations that all experience a reduction in grounded area by 2300 ( Figure   9 and Supplementary Figures S5-S7). In all scenarios there is limited grounding line retreat in the next 100 years (up to 2100), amounting to only −716 km 2 (−4800 to 2050 km 2 : 5 − 95%) change in grounded area in RCP 8.5. After 2150, grounded area decreases more rapidly, and the spread of change in grounded area increases within each forcing scenario. This coincides with the timing of a reduction in the rate of accumulation and increases in the rate of grounded ice discharge (Figure 8). This is particularly the case in RCP 8.5 (Figure 9), where more substantial ungrounding coincides with a sharp increase in uncertainty in ice discharge in 2200 and the increasing importance of sub-shelf parameters on projection uncertainty (Figure 7). This  Figures S5-S7). This suggests that these ice streams are prone to some ungrounding even with limited increases in sub-shelf melting. The likelihood of retreat and complete collapse of these ice streams increases dramatically with warming (RCPs 6.0 and 8.5). In RCP 8.5 perturbations become large enough for possible (but less likely) retreat of additional regions, in particular the Rutford and Evans Ice Streams draining into western part of the Ronne Ice Shelf and some retreat of ice streams feeding the Filchner Ice Shelf (Figure 9). Our separate additional simulation that validated the upper end of 500 sea level projections from our surrogate model (∆GMSL = 250 mm) shows that there is potential within our parameter space (beyond our initial ensemble) for the grounding line to retreat much further inland (Figure 9). This is characterised by runaway retreat of Möller and Institute ice streams, likely due to being topographically unconfined and rested upon a retrograde bed below sea level (Figure 9). This suggests that increases in sub-shelf melt with climate warming, have the potential to reduce ice shelf buttressing and force substantial grounded ice loss.

Discussion
Here, we have used an uncertainty quantification approach to assess the spread of future changes in global mean sea level (∆GMSL) contribution from the Filchner-Ronne region of Antarctica under different RCP emissions scenarios. We have taken a large and extensive sample of parameter space using a novel surrogate modelling approach, and our results show it is highly likely that, within the bounds of our input parameter space, ∆GMSL from the FR basin will be negative. Under RCP 2.6 510 forcing, during which atmospheric temperatures increase by 2 • C (Figure 4: in line with the targets of the Paris agreement), this region is likely to remain close to balance (accumulation ≈ discharge) and contribute to −19.1 to 9.57 mm (5-95%) of sea level rise by 2300. Under higher warming scenarios, projections of sea level rise become increasingly negative, but the uncertainties in these projections also increase dramatically. In the highest warming scenario (RCP 8.5), the FR basin could contribute anywhere between −103 to 26 mm (5-95%) to global mean sea level. Our projections are predominantly negative 515 due the the mitigating effect of increased accumulation with warming on sub-shelf melt driven increases in ice discharge.
Increases in precipitation across the Antarctic ice sheet, are known to have an important mitigating effect on the contribution to global sea level rise (Medley and Thomas, 2019;Winkelmann et al., 2012). Unlike other parts of West Antarctica, the in accumulation during the 21st century, driven by warming, can outweigh slow increases in ice discharge associated with sub-shelf melting (Figure 8). This may be enough to stabilise some of the major ice streams, in particular the Institute and Möller (see 50% line in Figure 9). In most cases, mass gain continues through to 2300, despite increases in the rate of ice discharge and slow-down in the accumulation trend, both of which are not enough to switch to mass loss/positive ∆GMSL by 2300. We continued our median RCP 8.5 simulation for a further 200 years with forcing held at 2300, and found that total 525 mass balance (accumulation − discharge) remained positive (mass gain) by 2500. These findings are consistent with previous regional modelling studies which showed that when applying both accumulation and sub-shelf melt anomalies, it is possible for accumulation to suppress the effects of increased sub-shelf melt and keep the FR basin close to balance (Wright et al., 2014;Cornford et al., 2015). This is in contrast to Antarctic wide simulations that do not impose changes in surface accumulation through time (Ritz et al., 2015;Schlegel et al., 2018) and hence project greater ice loss and sea level rise contribution from this 530 region. Overall, our results demonstrate that the FR basin is particularly sensitive to future accumulation changes, which are capable of stabilising this region in response to climate-warming.
Despite the important role that precipitation plays in suppressing ice loss from the FR basin, our sensitivity analysis reveals that the percentage increase in precipitation per degree of warming p, is the second largest contributor (approx 30% in RCP 8.5) to uncertainties in our projections. Hence, we can identify the representation of accumulation changes with warming in 535 ice-sheet models as a target area for further research, in order to better constrain projections of sea level contribution. However, modelling future precipitation trends is challenging and CMIP5 models themselves show large temporal and spatial variability in projected precipitation trends with warming (Tang et al., 2018;Palerme et al., 2017;Rodehacke et al., 2020). Despite this, it has been assumed that there is a general correlation between precipitation and temperature anomalies. As a result, temperature scaling of precipitation (as used in this study) is a common approach in ice sheet modelling, most of which use 540 a Clausius-Clapeyron relation, equivalent to a 5% spatially uniform increase in precipitation (Golledge et al., 2015;Gregory and Huybrechts, 2006;Garbe et al., 2020;DeConto and Pollard, 2016). But, as we have shown here, small variations in p can have a large effect on the spread of ∆GMSL projections, and better constraints on the value of p are needed. Recent work by Rodehacke et al. (2020) sought to improve Antarctic estimates for p, and found that the assumed correlation/or lack thereof between temperature and precipitation anomalies, has strong regional differences, which may invalidate the use of a simple, 545 uniform scaling. Instead they propose spatially variable values of p across the Antarctic Ice Sheet. Over the FRIS, they showed p = 4 − 6%, which is consistent with values used in this study, but may reach up to 10% in inland regions of the FR basin (Rodehacke et al., 2020). Thus, using a spatially uniform value for p, may lead to under or over estimates of precipitation across the catchment. This is even more important when conducting Antarctic wide simulations, and future studies should move towards using a spatially variable value for p, or ultimately conduct coupled atmosphere-ocean-ice-sheet simulations.

550
Our probabilistic projections have show that a negative ∆GMSL is most likely from the FR basin. Nonetheless, uncertainties associated with our input parameters reveal that it is also possible within our parameter space for the FR basin to contribute positively to sea level rise. This occurs predominantly under RCP 8.5 forcing, where the the long tail of the projections in Figure 5 reflects these more "extreme" (maximum of 332 mm by 2300) yet unlikely contributions to ∆GMSL. Hence, it is possible for sub-shelf melt to increase enough to outweigh 21st century accumulation, by forcing substantial increases in ice 555 discharge (Figure 8), and un-grounding of the ice streams feeding the FRIS (Figure 9). These high-magnitude contributions to sea level rise are characterised by the rapid retreat of the lightly-grounded Möller and Institute ice streams, which, once initiated, continues unabated across the Robin subglacial basin throughout our simulations (maximum GMSL simulation: Figure 9). Additional retreat occurs predominantly in the ice streams feeding the Ronne ice shelf (Evans and Rutford Ice streams and Carlson Inlet), suggesting that these regions are likely to be the dominant contributors to future sea level rise.

560
There is also some potential for un-grounding in the ice streams feeding the Filchner ice shelf, but this is less than elsewhere in the region (Figure 9). Retreat of the grounding line in our simulations is consistent with the magnitude and spatial pattern of retreat simulated in other studies in response to increased ocean-driven melt rates, with comparable projections of sea level rise contribution (approx 150-160 mm) to the high-end of our results (up to 300 mm) (Schlegel et al., 2018;Wright et al., 2014;Cornford et al., 2015).

565
Greater variability in ice discharge from the 22nd century onwards (Figure 8), coincides with an increase in the spread of our projections, suggesting that sub-shelf melt could strongly influence the regions potential contribution to sea level rise.
Indeed, our sensitivity analysis clearly reveals that ocean forcing parameters are the dominant component of uncertainties in our projections of sea level contribution from the FR basin (Figure 7). This sensitivity corroborates the well established theory that ocean forcing and the impact on sub-shelf melt rates, is a dominant, yet uncertain, driver of Antarctic Ice Sheet mass 570 loss. (Seroussi et al., 2019Cornford et al., 2015;Bulthuis et al., 2019). Of our ocean forcing parameters, the magnitude by which global atmospheric temperature anomalies are scaled to ocean temperature changes α, appears highly uncertain, amounting to 44% of projection uncertainty in RCP 8.5 (Figure 7). This is consistent with the high sensitivity to the value of α across the entire Antarctic Ice Sheet (Bulthuis et al., 2019). Going forward, given the impact of a linear scaling and the value of α on projection uncertainty, it may be more suitable to instead force ocean temperature changes with the results of CMIP ocean 575 models directly (e.g. the approach used in ISMIP6 proposed by Jourdain et al., 2020). In addition to uncertainties in ocean temperature forcing, it remains challenging to accurately represent ice shelf melt rates, and their sensitivity to temperature changes, in ice-sheet models. While melt parameterisations such as the PICO box model (Reese et al., 2018a) are a substantial advancement in our ability to efficiently apply sub-shelf melting in a physically plausible way, they remain a simplification of observed melt rate patterns, and those simulated by ocean models. Hence, it is possible we do not capture the same spatial 580 distribution or magnitude of melt in highly buttressed regions of the ice shelf as shown in observations, which could ultimately impact the (in)stability of the grounding line.
Alongside the mitigating effect of accumulation, using a linear scaling of ocean temperatures, and a simple melt parameterisation, may both be responsible for our simulations not projecting a substantial increase in sub-shelf melt or contribution to global mean sea level rise. Crucially, it appears that we are not capturing the regime shift from 'cold' to 'warm' cavity 585 conditions as seen in ocean model results (Hellmer et al., 2012;Hazel and Stewart, 2020;Hellmer et al., 2017). Simulations by Hellmer et al. (2012) showed relatively warm ocean waters could flush the ice shelf cavity, and increase the area-integrated (fixed ice shelf extent) basal mass loss from 80 Gt yr −1 to 1600 Gt yr −1 by the year 2100. By comparison, basal mass loss during our RCP 8.5 ensemble of simulations (integrated over the initial ice shelf area), reached a maximum of only 1200 Gt yr −1 some 200 years later (2300). Beyond our initial ensemble of 500 simulations, we may start to capture the same magnitude 590 of basal mass loss, which is reflected in the long-tail of positive contribution to sea level rise in RCP 8.5 projections. However, the key difference is that melt rates increase at a slow and steady rate over the 150 year period, and do not impose a rapid switch from cold to warm conditions which may be possible with a sudden flushing of warm water into the ice shelf cavity.
To fully capture such a regime shift, and the effects this would have on ice shelf thinning, loss of buttressing, and increases in ice sheet discharge, it is necessary to run fully coupled ice-ocean model simulations. Recent work performing coupled simula-595 tions on the FR region (Naughten et al., 2021) found that ice-shelf melt rates are unlikely to increase over the next century, and thus the region will have a limited contribution to sea level rise until ocean temperatures increase substantially (+7 • C). While more coupled modelling studies emerge, they are currently computationally feasible for regional configurations, and have yet to be accomplished on an ice sheet scale. In the meantime, melt parameterisations will remain important for future ice-sheet simulations, and so work should still focus on improving their ability to capture the physical behaviour of ocean models, as 600 well as the choice of ocean temperature forcing used to perturb those melt rates.
In contrast to atmospheric and oceanic forcing parameters, those related to ice flow dynamics in our model appear to play a less important role on uncertainties in projections of ∆GMSL. Consistent with other studies (Ritz et al., 2015;Bulthuis et al., 2019;Gillet-Chaulet et al., 2011;Alevropoulos-Borrill et al., 2020) we have shown that stronger non-linearity in our basal sliding and ice flow laws (increasing values of m and n) reduces the response time to a temperature perturbation, allowing 605 for faster delivery of ice to the grounding line, and a greater contribution to sea level rise. By varying the value of m in the Weertman sliding law we have captured a large range of amounts of basal sliding, and this has been fully integrated into our uncertainty analysis. However, this may not have captured the full spread of basal sliding possible under different sliding laws and/or spatially variable fields of m. Different sliding laws, (e.g. Budd sliding) may allow for even faster delivery of ice to the grounding line and thus greater contributions to sea level rise (Schlegel et al., 2018;Brondex et al., 2019). We are also 610 not accounting for any transient variability in our basal slipperiness and ice rheology fields, which is not yet captured in most ice flow models, but may additionally increase sea level rise. Progress is being made towards assessing the sensitivity of sea level projections to the choice of sliding law (Brondex et al., 2017(Brondex et al., , 2019Cornford et al., 2020). Future work will benefit from choosing the form of the sliding law (Ritz et al., 2015;Gillet-Chaulet et al., 2016), and/or determining spatially variable values for m Joughin et al., 2010), that best replicate regional observations of ice loss. This will help to constrain 615 uncertainties associated with the prescription of basal sliding, but this remains an active area of research.
Overall, our simulations have shown that regional increases in accumulation assumed with warming are likely to have an important stabilising effect on the ice loss from the FR basin under scenarios of future climate change. There is still some potential for a positive contribution to global mean sea level rise under high sub-shelf melt scenarios. This means that the sign of ∆GMSL projections from the FR basin cannot be fully constrained. Parameters driving both accumulation and sub-shelf melt-620 ing are highly uncertain and we identify them as priority areas for research, where more accurate parameterisations will help to constrain future projections, not only from the FR basin, but the entire Antarctic Ice Sheet. Future coupled atmospheric-oceanice-sheet simulations will help to more accurately capture feed-backs between the atmosphere and ocean on the evolution of the ice-sheet, but remain computationally challenging on a Pan-Antarctic scale. In addition to coupled simulations, model calibrations with observations can substantially reduce uncertainties and constrain future projections by narrowing the parameter 625 space for future simulations based on their fit to observations (Wernecke et al., 2020;Ritz et al., 2015;DeConto and Pollard, 2016;Ritz et al., 2015;Reese et al., 2020). As the number and timespan of observations increases, we will be able to better inistalise our ice sheet models to present-day conditions prior to future simulations. Overall, employing uncertainty quantification techniques in future studies will help to provide more robust estimates of potential sea level rise, and identify priority areas for future better constraining these projections.

Conclusions
This study set out to implement an uncertainty quantification framework (UQLAB) for use alongside the ice flow model Úa and use this to quantify uncertainties in projections of mass loss from the Filchner-Ronne region of Antarctica. We used a novel surrogate modelling approach to extensively sample an input parameter space to determine the forward propagation of uncertainties. Our probabilistic projections indicate that this region may not undergo dramatic ice loss under climate warming 635 scenarios, and instead have a negative contribution to global mean sea level rise. This is primarily due to the effects of increased accumulation assumed with greater moisture content in a warmer climate, that is capable of suppressing mass loss attributed to ocean-driven increases in sub-shelf melt rates. Despite this, we find that there is the potential, albeit highly unlikely, within the bounds of our input parameter space, for a substantial positive contribution to global mean sea level. In these high mass loss scenarios, sub-shelf melting increases enough to outweigh accumulation and force major retreat of some of the ice streams 640 flowing into the FRIS. Uncertainties associated with parameters driving accumulation and sub-shelf melt account for most of the spread of future changes in global mean sea level, and we highlight these as priority areas for constraining projections of ice loss. Future work would benefit from employing uncertainty quantification techniques similar to those used in this study, to fully assess the spread of future projections of sea level rise, not only from the FR basin, but across the entire ice sheet.
To estimate the rate factor (A) and basal slipperiness coefficient (C) for each of our randomly sampled combination of m and n we use the inverse capabilities of Úa to minimise the misfit between observed (u obs ) and modelled (u mod ) velocities. Observed velocities are MEaSUREs InSAR-Based Antarctica ice velocities (Version 2) from 1996 to 2016 and with a spatial resolution of 450 m Mouginot et al., 2012). Ice velocities were linearly interpolated onto our model mesh. Úa uses a standard inverse methodology in which a cost function J, which is the sum of a misfit (I) and regularisation (R) term, is 650 minimized. The gradients of J with respect to A and C are determined in a computationally efficient way using the adjoint method and Tikhonov-type regularisation. The misfit (I) and regularisation (R) terms are defined as: Prior to our uncertainty quantification routine (see Section 3.5) we generated a 75 member 'library' of inversions for every half integer between 2 and 9 for m and between 2 and 4 for n. For these 75 inversions we defined prior values as follows:Â = /τ n with = 10 −4 yr −1 and τ e = 80 kPa which for n = 3 givesÂ ≈ 2 × 10 −9 kPa −3 yr −1 equivalent to an ice temperature 665 of approx. −25 • C using an Arrhenius temperature relation (Cuffey and Paterson, 2010).Ĉ = u b /τ m b with u b = 10 m yr −1 and τ b = 80 kPa. This library of inversions was designed to make it computationally feasible to incorporate a model inversion before to every forward model run into our uncertainty analysis. For each of our randomly sampled values of m and n we select the closest inversion from our library as the a priori values for A and C. These priors provide a good initial estimate of the spatial fields of A and C which means the subsequent inversions need far less iterations to converge. After each inversion 670 we advected C beneath the ice shelf to avoid a sharp gradient in C downstream of the grounding line in the case of glacier advance. We note that the model calculated velocities for each model inversion will vary slightly based on the value of m and n used and the resultant fields of A and C. However, we find all ensemble members (N = 500) to provide an optimal fit to observations and that the misfit between observed and modelled velocities varies by only 1 m yr −1 , which is small with respect to measurement errors. based on observations. This can be done using Bayes theorem: where the posterior probability distribution of θ (a hyperparameter) given Y observations is equal to the likelihood ( ) of θ given Y multiplied by the prior probability distribution π(θ). We conduct this analysis on the four 'hyperparameters' used in the box model. These are: the time delay (τ ) and scaling coefficient (α) used to force changes in ocean temperature  (2010)). While some information exists on all these parameters, their bounds and distributions are not well known. The primary aim is not to find single point estimates for these parameters, but obtain an optimal range of parameter values that fit model predicted melt rates to observations. These posterior distributions are then used as input to our uncertainty analysis (see Figure 3). We 690 conduct this Bayesian optimisation using the tools in UQLab, including specifying a prior input model, surrogate models, and the Bayesian inversion itself (Marelli and Sudret, 2014).

B1 Priors
Prior probability distributions of our four parameters take into account any available information on their values before our Bayesian calculation. For the scaling coefficient (α) and time delay τ we used the values presented in Levermann et al. (2020) 695 (see outline in Section sec:oceanforcing) as the a priori information on the probability distributions of these parameters ( Figure   B1). While some ranges for the heat exchange coefficient (γ * T ) and overturning strength (c) , have been proposed (Reese et al., 2018a;Olbers and Hellmer, 2010), their probability distributions are unknown. Therefore, we use non-informative priors, i.e. we do not prescribe any prior information about these parameters and use a uniform distribution within the bounds c ∈ [0.1, 9] Sv m 3 kg −1 and γ * T ∈ [5 × 10 −6 , 1 × 10 −4 ] m s −1 given by Reese et al. (2018a); Olbers and Hellmer (2010) (see 700 Figure B1). However, we do know that these parameters are related, and certain combinations mean that the physics in the box model no longer hold Reese et al. (2018a). To specify the dependence between values of γ * T and c and ensure that values outside of these bounds are not sampled, we prescribe a Gaussian copula with a correlation of p = 0.9. To test the sensitivity of our posterior distributions to our specification of the priors we repeated the analysis for several prior distributions, with or without copulas and found our results were largely insensitive to our priors ( Figure B1).

B2 Surrogate modelling
Rather than use a single area-integrated value of sub-shelf melt, and to preserve some of the spatial distribution of melt rates across the shelf, we chose to tune our parameters to observations, using average melt rates within each melt box. As Bayesian analysis requires a large number of iterations to settle on a posterior distribution, we first construct five surrogate models to emulate calculated melt rates in each ocean box. For this we use the same method as our surrogate modelling for changes in 710 global mean sea level (see Section 2). We sample 2000 points from our prior probability distributions using Latin hypercube sampling and use these to directly evaluate sub-shelf melt rates in the box-model. Each simulation is run in a diagnostic mode for our nominal start year of 2000, using observations of topography. We assessed the performance of our surrogate models by taking a separate validation sample from the parameter space, and found a good fit between true and surrogate modelled melt rates (Supplementary Figure S8).

B3 Bayesian inversion
To derive posterior distributions for our hyperparameters in Equation B1 we require three things: 1) prior probability distributions, 2) observations of ice shelf melt rates, and 3) a likelihood function that specifies the likelihood of parameter probability distributions given observed and modelled melt rates and associated errors. Prior probability distributions have been outlined above and are shown in Figure B1. Secondly, we take observations of sub-shelf melt rates from Moholdt et al. (2014) and av-720 erage these within each of our five ocean boxes across the Filchner-Ronne ice shelf ( Figure B2). We then assume that average melt rates within each box are independent (uncorrelated) with one another and use a log likelihood function defined in the common format as: where i is the box number, y are observed melt rates, f (x) are modelled melt rates from the surrogate models, and ε is discrep-725 ancy term for the melt of each box. Errors associated with model physics are difficult to quantify, so we instead incorporate errors from both measurements σ 2 obs (Moholdt et al., 2014) and the surrogate models σ 2 pce . We then weight this error term (w) with the normalised box area with respect to the total ice shelf area.
We performed a number of sensitivity experiments in which we varied the specification of the discrepancy term in the log 730 likelihood function and found that our posterior distributions are similar, regardless of the choice of discrepancy (see Figure   B1).
Observations for the ocean box closest to the calving front showed high average melt rates (0.5 m yr −1 ) which is not related to the overturning circulation in the cavity, but instead seasonal warm surface water intrusions, and will not be replicated by the box model. We therefore choose to replace this with the box average melt from an initial run (using default parameters from  Table B1) and reduce the weighting so that it is largely excluded from the analysis. As most ice shelf melt takes place close to the grounding line, far-field melt rates are less important for the total mass balance of the shelf.
We specify our priors, likelihood function, and observed melt rates in the Bayesian Inference module in UQLab (Wagner et al., 2019). The posterior probability distributions (π(θ|Y )) for our hyperparameters (θ) are then estimated using Monte Carlo Markov Chains (MCMC). We use 1000 independent parallel chains (or sometimes referred to as walkers that move 740 randomly around the parameter space), the starting points of which are randomly sampled initial estimates for θ from our prior probability distributions (π(θ)). Then an initial step is made from the current position and the posterior probability distribution at that point (which is the product of the likelihood and prior probability see Equation B1) is accepted or rejected using an adaptive metropolis Hastings algorithm (Haario et al., 2001). This is based on whether they are in the right direction from the last sampled point, using an acceptance probability and then the posterior distribution is updated along the way using the 745 information accumulated so far. This process is repeated for 10000 steps for each Markov Chain, by which time the posterior distributions have converged. For more details on MCMC and the adaptive metropolis Hasting algorithm see the UQLab Bayesian Inference Manual (Wagner et al., 2019). We estimate whether the chains have converged on the same point using multivariate potential scale reduction factor (MPSRF: see Brooks and Gelman (1998)

B4 Posterior distributions
Our posterior probability distributions for all four parameters in the ocean box model are shown in Figure B1 and are input 755 to our uncertainty propagation (Figure 3). These show the distribution of possible values for these parameters that can lead to melt rates closer to observations than non-informative priors. Posterior distributions for τ and α remain close to their priors . There is a decreasing likelihood of the delay between increases in atmospheric and ocean temperatures from τ ≈ 10 to 100. The scaling coefficient is centered around 0.24, which is consistent with the scaling factor found to provide a good fit to CMIP5 model data (Taylor et al., 2012;Golledge et al., 2015). We note that the range for α is In contrast, distributions for γ * T and c have shifted significantly from their prior distributions, and in both cases favour values towards the lower end of prior ranges. Our point estimate for c (1.2 Sv m 3 kg −1 ) is close to the value proposed in Reese et al. (2018a), while γ * T is lower but the probability density function still extends to the value in Reese et al. (2018a) (Figure B1). To examine the performance of our Bayesian inversion, we compare the total area-integrated basal mass balance, i.e. mean 765 specific surface mass balance, (see Table B1) and spatially averaged melt rates within each box ( Figure B2a  To( • C) -1.66 -1.64 BMB (m yr −1 ) -0.087 -0.133 Table B1. Comparison between the forcing temperature and integrated basal mass balance across the total ice shelf for prior and posterior parameter values. In the prior case we use the parameter values for γT * and c proposed in Reese et al. 2018 and use mean point estimates from our prior distributions for τ and α. Our posterior case uses point estimates for all parameters from our posterior distributions.
(−0.13 m yr −1 ), primarily by improving the mean melt rates in boxes 2 and 3 to within the error of observations ( Figure B2a).

770
Additionally, we take our model optimised velocities (using an inversion with m = 3 and n = 3) and calculate sub-shelf melt rates from ice flux divergence assuming steady-state conditions, i.e. negligible SMB and no surface thinning/thickening. In this case, BMB is −0.15 m yr −1 , which suggests that our point parameter estimates are producing near-steady-state melt rates (−0.13 m yr −1 ).
We re-evaluate the surrogate model for the entire prior and posterior sample sets to examine the distribution of melt rates 775 in each ocean box. Figure B2b shows the range of initial melt rates (at start year 2000) close to the grounding line (Box 1) that would have occurred in our uncertainty propagation if we had chosen to use non-informative priors, particularly uniform distributions for γ * T and c. Given the prior information we have on sub-shelf melt rates (Moholdt et al., 2014), this range of melt rates suggests we would have been sampling unlikely regions of the parameter space. Our Bayesian analysis has successfully tightened the posterior distribution of melt rates, where 5-95% fall within the standard deviation of observations. We can now 780 be confident that our initial starting point for melt rates under the ice shelf at the beginning of our forward simulations is reasonable with respect to observations. We note that our approach estimates these distributions at a single snapshot in time and does not take into account variations in these parameters that may occur under future warming e.g. an increase in the strength of the overturning circulation in ice shelf cavities with warmer ocean temperatures.