Statistical emulation of a perturbed basal melt ensemble of an ice sheet model to better quantify Antarctic sea level rise uncertainties

Antarctic ice shelves are vulnerable to warming ocean temperatures, and have already begun thinning in response to increased basal melt rates. Sea level is therefore expected to rise due to Antarctic contributions, but uncertainties in its amount and timing remain largely unquantified. In particular, there is substantial uncertainty in future basal melt rates arising from multi-model differences in thermal forcing and how melt rates depend on that thermal forcing. To facilitate uncertainty quantification in sea level rise projections, we build, validate, and demonstrate projections from a computationally efficient statistical 5 emulator of a high resolution (4 km) Antarctic ice sheet model, the Community Ice Sheet Model version 2.1. The emulator is trained to a large (500-member) ensemble of 200-year-long 4-km resolution transient ice sheet simulations, whereby regional basal melt rates are perturbed by idealized (yet physically informed) trajectories. The main advantage of our emulation approach is that by sampling a wide range of possible basal melt trajectories, the emulator can be used to (1) produce probabilistic sea level rise projections over much larger Monte Carlo ensembles than are possible by direct numerical simulation 10 alone, thereby providing better statistical characterization of uncertainties, and (2) predict the simulated ice sheet response under differing assumptions about basal melt characteristics as new oceanographic studies are published, without having to run additional numerical ice sheet simulations. As a proof-of-concept, we propagate uncertainties about future basal melt rate trajectories, derived from regional ocean models, to generate probabilistic sea level rise estimates for 100 and 200 years into the future. 15

backforce on upstream glaciers which causes grounding line retreat (Konrad et al., 2018;Rignot et al., 2014), increases flow rate (Pattyn, 2018), depresses surface heights of grounded ice (Konrad et al., 2017), and ultimately impacts sea level. Antarctic ice loss is particularly susceptible to a positive feedback due to the so-called marine ice shelf instability (MISI) (Weertman, 1974;Schoof, 2007). Much of West Antarctica's ice is grounded below sea level, with a retrograde bed sloping downward 25 toward the interior of the continent. MISI theory suggests that increased basal melt rates beneath some key WAIS ice shelves (e.g., Pine Island and Thwaites) could result in an unstable grounding line retreat causing runaway ice loss for the entire region.
In fact, there is some evidence through observations and modeling that this process may have already been triggered (Rignot et al., 2014;Joughin et al., 2014;Favier et al., 2014). Forcing due to basal melt is therefore likely to become an increasingly dominant contributor to Antarctic SLR (Bulthuis et al., 2019). Despite its potential to contribute to sea level rise (SLR) vastly 30 more than any other single source (∼ 5 m West Antarctica, ∼ 60 m all Antarctica), and documented ice shelf thinning (e.g. Schroeder et al., 2019;Reese et al., 2018), Antarctica's contribution to future sea level remains highly uncertain (Heal and Millner, 2014).
The primary unknown in how basal melting will affect sea level rise is the uncertainty in sub-shelf melt rates themselves.
Future basal melting is uncertain because it depends on unresolved, coupled ice-ocean processes which are in turn driven by 35 a range of global and regional ocean and atmospheric conditions. Sub-shelf melting can be decomposed into several factors including: changes to the large scale circulation in the Southern Ocean and cross-slope exchange of warm water onto the continental shelf, changes to regional circulation on the shelf, and to the local circulation within ice cavities themselves.
Accurately modeling these changes requires high-resolution ocean models in order to link large-scale ocean circulation to subshelf melt (Pattyn, 2018;Asay-Davis et al., 2017). Development of these modeling capabilities is still considered the forefront 40 of current research. Most coarse-resolution global models do not have ice shelf cavities, and those that do have large uncertainty regarding changes to the influx of relatively warm Circumpolar Deep Water (CDW) mass into the ice shelf cavities. Ice loss resulting from CDW intrusions has already been observed in the Amundsen Sea region in West Antarctica (Hellmer et al., 2017;Pritchard et al., 2012). Mass loss in the Totten region in East Antarctica has also been linked to ocean circulation changes (Greenbaum et al., 2015;Wouters et al., 2015). There is further uncertainty in how ocean eddies -which are unresolved by 45 current standard-resolution climate models -transport heat to the Antarctic coast (Paolo et al., 2015;Stewart and Thompson, 2015). Finally, outside of basal melt uncertainty, there is deep uncertainty in glaciological dynamics (ie. no consensus on what processes to include in an uncertainty analysis, nor how). An example of this is whether ice fracture mechanics, such as the marine ice cliff instability (MICI), could dramatically accelerate ice loss (DeConto and Pollard, 2016;Edwards et al., 2019).
Efforts by the scientific community have surged in hopes of constraining the uncertainty bounds on future SLR from Antarc-50 tica (e.g. initMIP-Antarctica (Seroussi et al., 2019) and ISMIP6 (Seroussi et al., 2020)). The typical approach is to run large ensembles of ice sheet model simulations, perturbing different parameters for each run, and then estimating uncertainty based on the model spread (e.g. Golledge et al., 2015;DeConto and Pollard, 2016). Less conventional techniques have recently been applied, including the use of reduced statistical models (Kopp et al., 2016;Mengel et al., 2016;Fuller et al., 2017;Le Bars et al., 2017) or structured elicitation studies (Kopp et al., 2014;Little et al., 2013;Bamber and Aspinall, 2013). Others have 55 used semi-empirical dynamical models relating GMSL change to global temperature (e.g. Grinsted et al., 2010;Mengel et al., 2 https://doi.org/10.5194/tc-2020-178 Preprint. Discussion started: 17 August 2020 c Author(s) 2020. CC BY 4.0 License. 2016; Kopp et al., 2016). Others still have used very simple reduced-form mechanistic models such as the BRICK (Building blocks for Relevant Ice and Climate Knowledge) model  to simulate changes in global mean surface temperature and sea level as a function of perturbations to radiative forcing. Despite this assortment of methods, there is still deep uncertainty in how the ice sheet itself will respond to forcing in the future . 60 1.2 Benefits of statistical emulation of ice sheet models Kopp et al. (2017) notes that "Ideally, the integration of process models into probabilistic frameworks...would involve the development and use of fast models -or fast statistical emulators of more complex models -in a mode that allows Monte Carlo sampling of key uncertainties and the conditioning of uncertain parameters on multiple observational lines of evidence.
The development of such fast models or model emulators is an involved task." Statistical emulators, sometimes referred to as 65 surrogate models, can be used to fully explore parameter space that would otherwise be too computationally intensive for a process-based model. A typical statistical emulation approach is a response-surface formulation -discussed further in Section 2.3 -such as a Gaussian process or neural network, which interpolates the outputs of a perturbed-parameter ensemble of model runs across its input space (Kennedy and O'Hagan, 2001). Because high-fidelity numerical models are computationally expensive, only a limited number of simulations are typically available for the purposes of uncertainty characterization. The 70 goal of an emulator is to inexpensively predict, from a small training ensemble of an expensive computer model, the output that the model would produce if it were run at a new input setting that was too expensive to simulate. The emulator can then be run in much larger ensembles than the original numerical model in order to explore uncertainties.
Statistical emulators have been used in climate science for some time. For example, Hauser et al. (2012) trained and built a Bayesian artificial neural network with GCM output, using their emulators to calibrate climate models against seasonal 75 climatologies of temperature, pressure and humidity. This generated statistically rigorous probabilistic forecasts for future climate states. Not long after, groups began to apply such statistical methods to ice sheet models as a step toward building ice mass loss projections which included uncertainties. For example, Chang et al. (2014) used spatially-resolved synthetic observations (with data-model fusion) to create a probabilistic calibration of a Greenland Ice Sheet model. Recently, emulation has gained even more traction in the Antarctic ice sheet modeling community. Pollard et al. (2016) used a Bayesian technique 80 involving Gaussian process-based emulations and calibration to provide SLR envelopes based on a 3D hybrid ice-sheet model applied to the last deglaciation of WAIS (∼ 20,000 yr ago). Bulthuis et al. (2019) built an emulator of the continental ice sheet response to a comprehensive set of uncertainties over the next millennium using the Elementary Thermo-mechanical Ice Sheet (f.ETISh) model (Pattyn, 2017) at 20 km resolution. Edwards et al. (2019), built an emulator based on the ice sheet modeling (at 10 km resolution) of DeConto and Pollard (2016) in order to generate probabilistic projections for the Antarctic contribution 85 to sea level rise.
In this paper, we build on these methods by constructing, validating and testing an ice sheet emulator based on the state-ofthe-art Community Ice Sheet Model (CISM) version 2.1 (Lipscomb et al., 2019). Our work is novel primarily in its focus on ocean forcing uncertainty in combination with high-resolution glaciological modeling. Specifically: -We focus on uncertainty in ocean forcing, which is considered the most powerful driver of Antarctic sea level rise in the 90 coming centuries.
-We use realistic transient trajectories of basal melt rates that can be mapped back to ocean models, as opposed to more stylized forcings.
-Statistical emulation is applied to the input uncertainty of the (ocean) climate forcing to the ice sheet model, as opposed to the ice model's parameters.

95
-We use the Community Ice Sheet Model (CISM) at higher resolution (4km) than previously used to construct ice sheet emulators. Benefits of higher resolution ice sheet modeling include, but are not limited to, improved representation of grounding line locations and complex bedrock topography.
-CISM is spun-up in such a way that it is in steady state with the current climate conditions. As such, any forward runs are divorced from issues of drift, so the response can be attributed to forcing and not to internal ice sheet variability.

100
There are several advantages to our approach of statistically emulating an ice sheet model's response to ocean forcing uncertainty. Conventional approaches to providing basal melt boundary conditions to a standalone ice sheet model include: (1) parameterizing melt rates from the thermal forcing of a global climate model (Naughten et al., 2018;Golledge et al., 2019;Seroussi et al., 2020;Jourdain et al., 2019), (2) using basal melt rates calculated from regional ocean models with ice shelf cavities (Cornford et al., 2015;Timmermann and Hellmer, 2013), or (3) providing stylized forcings such as instantaneous or 105 linear ramp melt trajectories, such as those found in the MISMIP+ and SeaRISE community experiments (Asay-Davis et al., 2016;Bindschadler et al., 2013). Each choice is tied to a specific set of physical/modeling assumptions, and typically does not probe the deep uncertainties in these assumptions.
Our statistical emulation method, by contrast, is intended to be agnostic with respect to these assumptions: we design an ice model ensemble that highly samples a wide range of possible basal melt trajectories (a "space-filling" sampling strategy, to 110 be discussed), initially without consideration of where these trajectories come from or which trajectories are most physically plausible. After constructing a statistical emulator of this ensemble, we can then provide the emulator with basal melt assumptions derived from any scientific source. This has two main advantages: (1) If we change our assumptions about future basal melt rates -due to expert disagreement, new scientific discoveries, or simple sensitivity analysis -we can interrogate the emulator to obtain new Antarctic discharge and sea level rise projections, without having to construct a new set of ice sheet 115 model simulations with new ocean forcings. A corollary is that we can use the emulator as a component in a modeling pipeline that can predict new sea level rise distributions (and their downstream impacts such as coastal flooding) with respect to uncertainties in large-scale climate forcing. (2) We can sample the emulator as many times as we wish, allowing a more complete uncertainty characterization of the distribution of future sea level rise, including extensive sampling into the societally relevant low-probability, high-risk tails of the distribution. 120 4 https://doi.org/10.5194/tc-2020-178 Preprint. Discussion started: 17 August 2020 c Author(s) 2020. CC BY 4.0 License.

Outline
The remainder of this paper is structured as follows. In the Methods section, we describe the ice sheet model configuration, ensemble design and the methods to build and validate our Gaussian process emulator. We present results for the CISM ensemble, and then show a simple example of using the emulator to generate a probability distribution function (PDF) of sea level rise by propagating prior distributions of basal melt rate parameters based on fits to two ocean models under the A1B 125 emissions scenario. Figure 1 shows a schematic of the tasks necessary to design, build, validate, and test the Gaussian process emulator.

Ice sheet model configuration
We use CISM, a state-of-the-art 3D, parallel, thermo-mechanical model that runs on a regular mesh grid using a mixture of 130 finite-difference, finite-volume and finite-element methods. As a participant in various ice sheet model intercomparisons (e.g. MISMIP+ (Cornford et al., 2020), LARMIP (Levermann et al., 2019), ABUMIP , and ISMIP6 (Nowicki et al., 2020; Seroussi et al., 2020)), CISM has been shown to perform well compared to other higher-order ice sheet models, some that use resolutions of 1 km or higher in the region containing the grounding line. For the experiments described here, the model was run with the following options:
-A basal sliding law based on Schoof (2005), in the power-law limit where the effective pressure is equal to the ice overburden pressure.
-Application of basal melting to partially floating cells in proportion to the floating fraction of the cell, which is diagnosed 140 from the thickness and basal topography as part of the grounding-line parameterization.
-A no-advance calving criterion that holds the calving front close to its observed location.
-Surface mass balance from late 20 th century simulations with the RACMO2 regional climate model (van Wessem et al., 2018).

145
The model is spun up over 40,000 years, with the modeled ice thickness nudged toward observed thickness by adjusting a 2D basal friction parameter field beneath grounded ice and basal melt rates beneath floating ice. This inversion scheme is similar to that of Pollard and DeConto (2012)) and was used for the CISM contribution to initMIP-Antarctica (Seroussi et al., Task 1: Design Parameterize future basinspecific basal melt rate trajectories into a family of sigmoidal functions, and widely sample the parameters of the trajectories to generate a space-filling statistical design of ocean forcing scenarios. Task 2: Run CISM ensemble Run large Antarctic ice sheet model ensemble over the statistical design, applying basal melt rate anomalies generated from parameterized (500 members, 200 years, 4km resolution).
Task 4: Generate priors Use a small set of regional ocean model projections of basal melting as informative probability distributions over the space of ocean forcing trajectories.
Task 5: Emulation Retrain the emulator on the full 500-member ensemble at key years. Perform Monte Carlo sampling of the emulator predicted future sea level rise over the probability distributions of ocean forcing.   2019) and the ISMIP6 projections (Nowicki et al., 2020;Seroussi et al., 2020). We note that there is no hydrology in the basal friction field, and the basal melt field is noisy, compensating for other errors in the model or observations. However, this 150 work is a proof-of-concept, and the emulator techniques used here would apply equally well to simulations with more realistic physics. Nudging is strong in the first half of the spin-up, and then tapered off for the second half. The ice thickness gradually approaches a quasi-steady state as basal friction parameters and internal temperatures evolve. The model is run on a uniform 4 km grid, resulting in a spun-up state with excellent agreement between observed and modeled surface velocity (Fig. 2), ice shelf extent, and ice thickness (Fig. 3).

155
A control run starting from the end of spin-up and going forward 1000 yr (not shown) indicates that there is very little drift in the ice sheet moving forward and that most changes in ice thickness will be a result of forcing as opposed to internal variability or model drift. This is not fully realistic, since the real ice sheet is never truly in equilibrium with the climate, particularly if current observations are used to tune the model. As a result, we consider this model spin-up to be closer to a pre-industrial state, rather than a modern ice sheet. Therefore, henceforth, we do not explicitly state the year corresponding to SLR projections. 160 Rather, we refer to our SLR projections as relative to the number of years run forward in the model from the end of spin-up.
As a result, the sea level rise projections are not tied to a particular year in the future. Rather, they are meant to show that the emulator is a powerful and useful tool, and SLR predictions are considered a proof-of-concept.

Ensemble Design
We use ocean model data from Timmermann and Hellmer (2013) and Cornford et al. (2015) to inform the types of possible 165 basal melt rate trajectory shapes for 200 years forward. The forcing data for the ocean models are generated with the global climate models HadCM3 (Gordon et al., 2000;Collins et al., 2001) and ECHAM5 (Roeckner et al., 2003)) under the A1B and E1 emissions scenarios. (A1B is a moderate scenario similar to Representative Concentration Pathway 6 (RCP6), and E1 is an optimistic, aggressive mitigation scenario.) These are then dynamically downscaled by two high-resolution atmosphere models (RACMO2 and LMDZ4) and two ocean models: the medium resolution BRIOS (Bremerhaven Regional Ice-Ocean 170 Simulation) (Timmermann et al., 2002) and the higher resolution FESOM (Finite-element Sea ice-ocean model) (Wang et al., 2014). We note that the availability of ocean data, particularly melt rate data beyond 2100, was very limited when this study was carried out.
For each CISM ensemble member, we apply a unique basal melt rate anomaly to each of five basins, following the LARMIP region delineation (Levermann et al., 2019). The regions are the Antarctic Peninsula, Weddell Region, East Antarctica, Ross 175 Region, and Amundsen Region (Fig 4). We find that we can accurately capture the behavior of all modeled melt rate shapes 8 https://doi.org/10.5194/tc-2020-178 Preprint. Discussion started: 17 August 2020 c Author(s) 2020. CC BY 4.0 License. with a sigmoidal function. Figure 5 shows that the sigmoids (colored curves) do an excellent job of characterizing the variety of shapes seen in the ocean model data in each basin. The equation describing the sigmoid is: where Equation 1 is a function of three independent parameters: t 0 (inflection point of turnover), τ (timescale of turnover), and 185 M max (melt rate to which the function asymptotes). However, because we are only able to constrain the melt rate 200 years into the simulation (we do not have simulations that go out to infinity), we must invert Eq.1 to be in terms of the melt rate at year 200 of the simulation, M 200 . Using Eqs. 1-5 we can derive coefficients A and B as a function of M 200 . In doing so, we are able to describe the basal melt trajectories in terms of three parameters: t 0 , τ and M 200 . The emulator is built on these three parameters.

190
We use a quasi-random Sobol' sequence as a space-filling design method for each parameter (Sobol', 1967). Compared to pseudo-random Monte Carlo sampling, space-filling designs reduce the likelihood of clustering leading to uneven sampling.
Sobol' sequences also have an advantage compared to some other common space-filling designs, such as Latin hypercubes (McKay et al., 1979), in that they are sequences designed to fully cover the parameter space at each point, recursively filling the space more densely as points continue to be added to the sequence. This sequential feature of quasi-random sampling allows 195 us to extend the sequence if we desire more ensemble members, while maintaining a space-filling design. A Latin hypercube, by contrast, is not a sequential design and its size must be specified in advance: more points cannot be later added to an existing design without violating the properties of a Latin hypercube.
The ranges for t 0 and τ are determined by maximizing the space-filling properties of the parameters whilst capturing all of the sigmoidal characteristics seen in the modeled ocean melt rate projections in Timmermann and Hellmer (2013) and Cornford 200 et al. (2015). Therefore, using a Sobol' sequence, we sample uniformly between the following lower and upper bounds for the An example of some basin-specific Sobol' generated sigmoids that would be fed to CISM are shown in Figure 6. For each CISM run, a random basal melt rate curve generated with the Sobol' sequence is chosen for each basin. We have chosen not to assume any correlation between basins. Of course, general ocean warming occurs with global warming, but there may be unique regional circulation patterns that could cause very different basal melt rates in one region compared to another. 215 We run an ocean-forcing-perturbed CISM ensemble (with 500 members) over the entire Antarctic domain, where CISM surface mass balance is a climatology, and unique ocean melt rate anomalies are added to the background basal melt rates from the end of the CISM spin-up. Since the spun-up basal melt rates are resolved along the ice draft, the original state of basal melt rates has a realistic character, in that it varies along a 2-D surface that is a function of depth so has spatial variability with increasing basal melt rates near the grounding line (Fig 7). We note that by imposing a uniform basal melt rate anomaly 220 for each basin, we are neglecting to account for complex patterns in sub-shelf ocean circulation changes or depth-dependence in the anomaly itself. However, most available melt rate projection data, particularly any estimates that go out 200 years, are regional averages. Furthermore, we must limit the number of parameters we vary in order to be able to run a large enough ensemble to appropriately sample the parameter space and subsequently build an emulator.

Gaussian Process Emulator of the CISM ice sheet model 225
Our knowledge (through observations and modeling) of Antarctic ocean melt rates is changing rapidly, giving rise to frequently revised estimates of sub-shelf melt rates (e.g. through efforts such as ISMIP6 (Nowicki et al., 2020;Seroussi et al., 2020)).
Model emulators, or fast surrogates, are useful tools as they allow quick propagation of prior parameter uncertainties through the dynamics of the glaciological system without the need to re-run large ensembles of ice sheet models every time new information about those uncertainties becomes available. Instead, very large Monte Carlo ensembles of surrogate predictions 230 from the emulator, which approximate the ice sheet simulator, can be run in order to highly sample the input uncertainties.  The emulator constructed here predicts a single output (Antarctic SLR anomaly in a specified year of simulation), as a function of a 15-dimensional input vector representing the perturbed regional basal melt trajectories (5 regions x 3 sigmoid parameters, (t 0 , τ , M 200 )). The training set is the full 500-member CISM ensemble.
Statistical emulation in this paper is of the response surface type (Box and Wilson, 1951). Our training ensemble -con- is computed via a statistical regression procedure. Usually the smoothness of data being fit is estimated as part of the interpolating procedure. A Gaussian process's interpolating surface is 'optimal' in a technical sense, going back to the geostatistics 245 literature where it is known as 'kriging' (Matheron, 1962;Krige, 1951). The modern Bayesian interpretation of Gaussian processes provides a second useful statistical feature, namely they can provide the uncertainty in their predictions, i.e., a built-in estimate of interpolation error.

Emulator Validation
To validate the emulator, we randomly withhold 20% of the ensemble members, and build an emulator based on the remaining respectively, against the withheld CISM output (Fig 8). For year 100, 8% of hold-out validation points lie outside the 2-sigma predictive intervals, which is plausible (one would expect 5%). Based on these highly correlated results, we are confident in the emulator's ability to approximate CISM's SLR output 100 and 200 years into the simulation.
In this section we present results from the CISM ensemble, showing the SLR contributions at year 100 and 200. We note that the SLR results from the ensemble itself should not be considered physically realistic as the prior melt parameters are uniformly sampled (as described in the Methods Section).
We also present results from a simple example illustrating the propagation of parametric uncertainties (M 200 , t 0 and τ ) de-260 scribing melt rate trajectories (by basin) derived from ocean model projections. These priors are then run through the emulator, resulting in a PDF of SLR at year 100 and 200.

CISM ensemble results
The CISM ensemble consists of 500 members where each member is forced by five melt rate trajectories, one in each basin. Figure 9 shows the sea level rise time series (top) and corresponding mass loss (Gt/yr) (bottom), resulting from the full 500-265 member ice sheet ensemble (blue shading), with the ensemble mean shown in red. The distributions of SLR at year 100 are more constrained (ranging from 0.5 -96 mm) than those in year 200 (ranging from 33 -543 mm) (Figure 9 top inset). Again, we note that these SLR projections are not physically meaningful since the parameter sampling over which the ensemble is created is uniform. The ensemble is designed purely to be used for the creation of an emulator.

270
We test the emulator by running a simple example of prior distributions (based on the FESOM and BRIOS ocean model data described in Timmermann and Hellmer (2013) and Cornford et al. (2015)) through the emulator. We produce a probability density function (PDF) of SLR at years 100 and 200. We generate a distribution of prior parameters based on the 'best fits' found by a least-squares optimization. In the cases that the best fit to a sigmoid has one of the parameters (t 0 , τ , or M max ) near the edge of the allowed parameter range (defined in Section 2.2), we find a fit that is within 1-sigma of the best fit which 275 relaxes the parameter away from the edge of the allowed range. The sigmoid parameters that describe the data fits are then used to generate a normal distribution that serves as the prior.
The parameter prior distributions in Fig. 10 can be combined and presented as a distribution of sigmoid anomalies (Fig 11).
By propagating the priors in Fig 10 for the A1B scenario through the emulator, we can predict SLR probability distributions for year 100 and 200 (Fig 12). The likeliest SLR at year 100 is found to be ∼ 7.5 mm. This falls within the range of century-280 scale future projections of sea level rise from Antarctica in the literature, albeit on the low end of most estimates. The SLR uncertainty in year 100 ranges from 1 mm (2.5%) to 31 mm (97.5%). At year 200, the SLR uncertainty ranges from 37 mm (2.5%) to 255 mm (97.5%).
As explained in the Model Configuration section, we do not consider this to be a prediction for the year 2100, simply because our spin-up is more in line with a pre-industrial state rather than a modern ice sheet. In the future, one option to eliminate this 285 problem might be to run a historical simulation before the perturbation simulations in order to bring the model into a more sensitive state. However, preliminary attempts to do this with CISM have been focused on global climate model output, which  spond to the 'best fit' to the ocean model data in Timmermann and Hellmer (2013) and Cornford et al. (2015).
did not capture the recent melting in the Amundsun Sea. Another reason for the small response is that the melt rate perturbation is applied uniformly over the basin. Focusing melt near the grounding line (for a given basin mean) would give greater retreat (Lipscomb et al., in review). Section 4 elaborates on both of these points.

Weddell
Year Basal melt rate anomaly (m/y) Figure 11. Prior basal melt rate anomaly (m/yr) trajectories as seen by CISM for each region: Amundsen Region, East Antarctic ice sheet (EAIS), Antarctic Peninsula, Ross Region and Weddell Region. Red lines show mean of distribution, blue shaded zones correspond to 50th, 80th and 90th % confidence intervals (from dark blue to light blue respectively), orange lines correspond to best fits (or within 1 sigma of best fit to relax parameters from the edge of their allowable range) to ocean model data.
19 https://doi.org/10.5194/tc-2020-178 Preprint. Discussion started: 17 August 2020 c Author(s) 2020. CC BY 4.0 License. Over a range of future melt-rate trajectories derived from a small collection of high-resolution regional ocean models, the CISM model projects from 1 mm (2.5%) to 31 mm (97.5%) of Antarctic SLR in year 100 and 37 mm (2.5%) to 255 mm (97.5%) in year 200. We caution that these results should be interpreted as a proof-of-concept of a method to quantify SLR 310 uncertainty with respect to uncertainties in ocean forcing, rather than a reliable SLR projection tool. There are several avenues for improvement.
One important limitation of this study is that we apply a uniform melt rate perturbation to an entire basin, neglecting to account for melt rate depth-dependence. This assumption could be relaxed in future work, for example, if we were to add (or swap a sigmoid parameter with) another parameter for ice-ocean physics parameters that control the extent to which melting 315 is focused near the grounding line. Further improvements would be to increase the ice sheet model resolution, use more realistic melt parameterizations, and include novel physical mechanics such as hydrofracture and cliff collapse (DeConto and Pollard, 2016). Eventually, the ice sheet emulator could be included in a larger system linking different sources of uncertainty with multiple emulators, such as a SLR/coastal flooding integrated assessment. As noted previously, the spin-up we use is advantageous because it is in equilibrium with the modern forcing and therefore allows us to isolate the effects of forcing as 320 opposed to model drift in forward runs. However, the drawback of such a spin-up is that one expects the present-day ice sheet to be out of equilibrium with present-day ocean forcing. Therefore, the model could have a lag in the response to basal melt rate anomalies, and therefore underestimate future sea level rise.
One way to address this shortfall in future work might be to insert a 'historical forcing' to the spin-up in order to ramp up the ice sheet to a more sensitive state before applying further basal melt anomalies.  (Newman et al., 2019), sub-shelf circulation and melt processes are still largely unknown in the existing glaciological literature 330 that focuses on land ice dynamics (Ritz et al., 2015;Pollard et al., 2016).
Another issue we encountered was the lack of data available to use as probability distributions over future basal melt rate tajectories. This was particularly difficult given our attempt at projecting out 200 years instead of the standard 100 years where most CMIP5/CMIP6 simulations terminate. One possibility to collect more priors is in the case of moving to a thermal forcing spin-up, such as in Lipscomb et al. (in review). In such a scenario, the ocean forcing consists of thermal forcing anomalies 335 (derived from ocean temperature and salinity) instead of basal melt rates. The ice sheet model has a parameterization to convert thermal forcing to melt rates instead of being given melt rates directly. This approach was used for ISMIP6 Antarctic projections. With such a method, one could use much more data from the CMIP6 dataset, particularly since the ISMIP6 effort has already generated thermal forcing files for certain projections/models (Jourdain et al., 2019). Again, these data cease at year 2100 but would be a convenient starting point.

340
Our method for deriving continuous probability distributions over melt rate parameters from a small collection of ocean models is somewhat ad-hoc. One limitation is that there are many possible, equally good, parameter fits of a sigmoidal curve to ocean model melt rate trajectory data. Selecting the 'best' fit found by least-squares optimization leads, in some cases, to parameter estimates on the boundary of the plausible range of values, requiring manual adjustment of the fits. A more principled statistical approach would calculate the full distribution of possible parameter fits to each ocean model, and combine 345 the distributions found for each ocean model into a single multi-model prior distribution using hierarchical Bayesian model fusion as in Jonko et al. (2018).
The general issues described above, particularly how to move from coarse (global climate model scale) knowledge of ocean temperatures to higher resolution sub-shelf melt rates, is a confounding issue in the current state of Antarctic ice sheet modeling.
We need global and regional ocean models to help address how ocean circulation will change, as well as how eddies transport 350 water from the open ocean to the continental shelf. Not only do we not know how general increases in ocean temperatures will translate to sub-shelf melt rates, but also changes in ocean circulation could impact the transport of relatively warm water to the continental shelf, thereby increasing sub-shelf melt rates as well (e.g. CDW). Progress in this direction will require larger ensembles of high-resolution regional and global ocean models that sample a wide range of climate scenarios driving Southern Ocean circulation change and variability. Indeed, regional ocean models would be a good target for emulation since these 355 models are particularly expensive, especially when run beyond a few decades.
Code and data availability. Code and Data are available at: https://github.com/mberdahl-uw/EmulatorPaper and MB ran them. BL developed and provided an ice sheet spin-up. NU wrote the statistical model that was run by MB. MB prepared the manuscript with contributions from all co-authors. 23 https://doi.org/10.5194/tc-2020-178 Preprint. Discussion started: 17 August 2020 c Author(s) 2020. CC BY 4.0 License.