Assessing Uncertainty in the Dynamical Ice Response to Ocean Warming in the Amundsen Sea Embayment, West Antarctica

Ice mass loss from the Amundsen Sea Embayment ice streams in West Antarctica is a major source of uncertainty in projections of future sea level rise. Physically based ice flow models rely on a number of parameters that represent unobservable quantities and processes, and accounting for uncertainty in these parameters can lead to a wide range of dynamic responses. Here we perform a Bayesian calibration of a perturbed parameter ensemble, in which we score each ensemble member on its ability to match the magnitude and broad spatial pattern of present‐day observations of ice sheet surface elevation change. We apply an idealized melt rate forcing to extend the most likely simulations forward to 2200. We find that diverging grounding line response between ensemble members drives an exaggeration in the upper tail of the distribution of sea level rise by 2200, demonstrating that extreme future outcomes cannot be excluded.


Introduction
Despite considerable advances in physically based models of ice dynamics over the last decade (Pattyn et al., 2017), there are still large uncertainties in the projections of future sea level rise from the Antarctic ice sheets. One major focus of uncertainty is the dynamic response of fast-flowing ice streams in regions that are grounded well below sea level, in particular the Amundsen Sea Embayment (ASE). Specifically, there is uncertainty regarding the onset, speed, and extent of large-scale grounding line retreat given the marine ice sheet instability theory (Favier et al., 2014;Joughin et al., 2014;Seroussi et al., 2014). Quantifying likely sea level rise over the coming centuries is critical to the adequate provision of coastal defenses.
The aim of this work is to demonstrate how spatial data of present-day observations can be used to calibrate an ensemble of ice flow model simulations, in order to construct a probability distribution of future sea level rise from the ASE. Quantification of uncertainty has been an integral part of global climate model projection for a number of years and features heavily in the Intergovernmental Panel on Climate Change assessment reports (Collins et al., 2013). However, only in the last few years has the ice sheet modeling community begun to formally consider uncertainty when estimating future contributions to sea level (Applegate et al., 2012;Chang et al., 2014;Edwards et al., 2014aEdwards et al., , 2014bEdwards et al., , 2019Gladstone et al., 2012;Levermann et al., 2014;Little et al., 2013;Ritz et al., 2015;Ruckert et al., 2017;Schlegel et al., 2018;Tsai et al., 2017). This delay is due, in part, to computational issues making it difficult to produce sufficiently large ensembles of simulations to investigate parameter uncertainty with available computational resources (Chang et al., 2014).
Several of the previous studies that do consider uncertainty focus on Greenland, where the fate of the ice sheet tends to be dependent on the modeled relationship between surface mass balance and surface elevation (Applegate et al., 2012;Edwards et al., 2014aEdwards et al., , 2014b; ocean-driven dynamics, while under-resolved, play a less important role in the ice sheet's behavior (Fürst et al., 2015;Goelzer et al., 2018), compared with the ASE. Studies that focus on Antarctica vary in model resolution, complexity, and spatial extent. Many Antarctic-wide studies use low-resolution models, which has consequences for the treatment of grounding line migration, often relying on parameterization (DeConto & Pollard, 2016;Edwards et al., 2019;Levermann et al., 2014;Ritz et al., 2015).
Depending on the magnitude of the step change in basal sliding and melt rate at the grounding line (Gladstone et al., 2017), explicitly simulating fine-scale grounding line dynamics (i.e. without relying on parameterization) requires sub-kilometer grid resolution at the grounding line , which is computationally expensive. We use the BISICLES ice flow model, which relies on adaptive mesh refinement, where the vicinity of the grounding line is modeled at a considerably higher resolution (250 m) than the interior of the ice sheet (4 km). We focus on a smaller area -the ASE, rather than the whole of Antarctica -as this region is likely to dominate the Antarctic mass loss signal in the next one to two centuries (DeConto & Pollard, 2016;Levermann et al., 2014;Ritz et al., 2015). These decisions allow us to perform a sufficient number of simulations using a sophisticated, high-resolution model to explore the likely range of dynamic response of the ASE to an idealized increase in sub-ice shelf melting.
Firstly, we perform a Bayesian calibration of a perturbed parameter ensemble of ice sheet model simulations, by comparing the model results with observations of surface elevation change. We then extend the calibrated ensemble to 2200 using an idealized melt rate forcing and explore the uncertainty in the ice sheet response given this forcing scenario.

Model and Observation Data
The perturbed parameter ensemble is described in Nias et al. (2016); here we will outline the relevant details. The BISICLES ice flow model is initialized to present-day conditions by performing an iterative procedure to find unknown quantities such as the basal traction coefficient (C), the ice viscosity stiffening factor ( ), and the sub-ice shelf melt rate (M b ) (Nias et al 2016). C and are found by solving an inverse problem given observations of velocity . M b is determined from the flux divergence over floating ice, but it is parameterized to be spatially smooth and to have the highest rates close to the grounding line. Two alternative initial states are created; one that uses the Bedmap2 geometry (Fretwell et al., 2013) and one that modifies the bed topography and grounded ice thickness to smooth spurious thickening signals in the flux divergence, which have been attributed to incorrect thickness measurements (Morlighem et al., 2011;Nias et al., 2018).
In total there are four optimal parameter sets, for all combinations of the two bedrock geometries and two Weertman sliding laws (m = 1 and m = 1∕3), which form the basis of the ensemble. Nias et al. (2016) use Latin hypercube sampling to create 64 distinct parameter vectors in which C, , and M b vary between a halving and a doubling of the optimized values, which, with the addition of the optimal member and six end members, produces a 284-member ensemble. The 50-year simulations are run under present conditions, that is, there is no time-dependent climate forcing.
Observed rates of surface elevation change (dh∕dt) over grounded ice are obtained from swath processing of CryoSat-2 radar altimetry measurements from 2010 to 2015 (Gourmelen et al., 2017). This processing technique is able to capture thinning rates in the swath rather than just the point of closest approach (POCA). In doing so it provides a greater spatial coverage of dh∕dt measurements compared to traditional POCA technique (Foresta et al., 2016). The spatial resolution of the data is approximately 500 m.

Bayesian Calibration
Bayes' theorem states that the posterior probability distribution (P[ |Y] -the probability of given Y ) is proportional to the prior probability distribution (P[ ] -the probability of ) multiplied by a likelihood function (P[Y| ] -probability of Y given ): (1) In other words, we are trying to find the probability of an event (e.g. a magnitude of sea level rise) produced by a particular parameter vector (e.g. the scaling factors used to vary C, , and M b ), given observations Y (e.g. dh∕dt). We are not trying to find a single estimate of , rather a distribution.
Each ensemble member is assigned a likelihood score based on discrepancies between the model output and observed dh∕dt, assuming Gaussian, independent errors (Edwards et al., 2014b). The likelihood score s j for the j th simulation in the ensemble is where f is the modeled dh∕dt, and z is the observed dh∕dt, and i is a spatial index.
2 is the discrepancy variance, which is a combination of observational error and structural error and represents the mismatch between the model, given the optimum parameter set and the real world (Edwards et al., 2019;Murphy et al., 2009). Observational error is found from the covariance matrix of the parameters used to derive the swath-processed dh∕dt (Foresta et al., 2016). Structural error has numerous sources related to the structural properties of the model; for example missed physical processes, spatial resolution of the grid, and the numerical representation. Structural error is poorly constrained, and so we conservatively assign it a value of double the observational error ( Figure S4).
Often in model-data evaluation, spatial comparisons are made at every available location. While this is appealing in terms of maximizing the number of data points, the spatial correlation inherent in most environmental variables means that this tends to overly penalize models in regions of coherent spatial patterns. The model error is "double counted" for each neighboring grid cell, even though they arise from a common source. One approach is to model this spatial correlation explicitly, but this is challenging and requires assumptions about the precise features of grid cell-to-cell correlations everywhere. A more common approach is to remove the spatial correlation by averaging or sub-sampling the data at a spatial scale at which they are reduced so the model-data discrepancies are sufficiently independent.
Using semi-variograms, we empirically investigate the length scales at which the covariance is reduced to an acceptable value and use this to decide upon an appropriate spatial scale on which to sample the discrepancies. We find this to be approximately 100 km, both in the x-and y-directions (see the supporting information).
The score for each simulation is normalized to create a weight w, where the simulation with the lowest discrepancy with observations has the largest weight. These weights, which are akin to P(Y| ) in equation 1, are used to weight the prior distribution, P( ), to produce a calibrated (posterior) distribution of sea level contribution, P( |Y) (Figure 1a). The most likely (modal) sea level rise estimate according to the prior distribution is 0.26 mm/year (50-year mean), which shifts to 0.30 mm/year in the posterior distribution (Table 1). The similarity between these estimates and observed rates of mass loss from the ASE (Figure 1a) indicates that these independent methods for quantifying present-day mass loss are in good agreement; whether it is BISICLES tuned using velocity observations (the prior), the observed spatial field of surface elevation change (the posterior), or the methods used to estimate total mass loss from the ice sheet (vertical lines, Figure 1a). The spread of the posterior distribution is reduced from the prior distribution, indicating that the calibration process is useful for reducing uncertainty in sea level rise projections. Future work could test the impact of using different types of observational data in the calibration process. For example, maps of observed velocity change could be a good candidate, as the dynamic signal is not influenced by changes in surface mass balance.

Extended Simulations
Of the 284-member ensemble, 187 members were within the 5-95% probability interval of the 50-year mean rate of sea level contribution (0.02-0.72 mm/year, Table 1). We exclude the extremes because their implied rates of elevation change perform poorly in the comparison with present-day observations (at the 10% probability level). Satellite observations have consistently shown that the ASE has been losing mass (McMillan et al., 2014;Mouginot et al., 2014), and so it is particularly appropriate to discard those members with mass gain. In addition, previous regional modeling has suggested that a linear viscous law is not suitable for describing sliding over bedrock and is prone to underestimate the sensitivity to changes in basal traction at the grounding line (Joughin et al., 2009(Joughin et al., , 2010. Therefore, given limited computational resources, we chose to extend to 2,200 only the 71 ensemble members that fall within the 5-95% probability interval of sea level contribution and use the non-linear (m = 1∕3) Weertman sliding law.
The sub-ice shelf melt rate is perturbed within the perturbed parameter ensemble described above. In addition, we apply a simplified melt rate forcing anomaly to ensure the direction of change is consistent with the expected behavior. In the Amundsen Sea, we expect there to be an increase in sub-ice shelf melt rates over the coming centuries due to increasing Circumpolar Deep Water (CDW) intrusions onto the continental shelf, as well as potentially direct warming (Holland et al., 2019;Timmermann & Hellmer, 2013). Therefore, we apply an idealized melt rate forcing, based loosely on regional ocean modeling, given a "business-as-usual" emissions scenario (Timmermann et al., 2002). The ice shelf averaged melt rate anomaly increases linearly to 15 m/year by the end of the 21st century-as ice shelf contact with the CDW increases-and remains constantly elevated in the 22nd century-representing continued CDW intrusion. The mean melt rate anomaly is in addition to the melt rates of the perturbed parameter ensemble (∼5-20 m/year mean) and is distributed to be highest near to the grounding line. Further details about the melt rate forcing can be found in the supporting information.
All the extended simulations continued to lose mass from the ASE, and by the end of the two centuries the modal contribution is 12 cm sea level equivalent (Table 1). The probability distributions of cumulative sea level contribution (Figure 1b) broaden over time, particularly in the upper tail of the distribution where the contribution in the second century is larger than in the first. This super-linear response persists, even when the melt rate remains constant in the second century. However, other simulations maintain an approximately linear response at a lower rate throughout the 200-year simulations. These two response types can be seen in Figure 2c; approximately 40% of the simulations exhibit a super-linear trend (blue lines).
By the end of the 21st century, all ensemble members have experienced a reduction in the total ASE grounded area (Figure 2a). However, retreat is not ubiquitous across all ice streams: Some members result in grounding

Discussion
During the 200-year simulations, the high-end ensemble members diverge from the modal behavior, creating a skewed distribution toward higher values of sea level contribution (Figure 1b). This is despite the high-end tail of the original ensemble being down weighted (Figure 1a), resulting in the most extreme simulations being removed from the longer century-scale simulations. Non-normal distributions, characterized by a long tail at the high end, are also found in other studies (DeConto & Pollard, 2016;Edwards et al., 2019;Kopp et al., 2017;Levermann et al., 2014;Robel et al., 2019).
The super-linear response of the high-end members means that while the mode of the distribution increases linearly, the total sea level contribution after 200 years is approximately double the total after 100 years-the 95th percentile increases disproportionately (Table 1). Given our idealized melt rate forcing, there is 5% probability that the ASE will contribute more than 12 cm of sea level rise by ∼2100 and 42 cm by ∼2200. This is in contrast to a study by Ritz et al. (2015), in which the response at the 95th percentile is quasi-linear, with 25 cm of sea level rise from the ASE in the 21st century and 48 cm by 2200. In their model, the representation of ice dynamics is simpler and at a lower resolution (15 km) than the model used here. In particular, the grounding line retreat is imposed rather than computed, which may dampen non-linear behavior.
We find that the grounding line behavior regulates the linearity of the sea level response. Indeed, the long tail at the high end of the sea level rise distribution is mirrored in the grounding line retreat: Some simulations achieve extreme retreat, whereas many simulations experience more modest retreat, and advance is limited (Figure 2). As the ice stream grounding lines retreat further into the deep basins they inhabit, the flux across them increases with ice thickness and with the lengthening of the fluxgate. The relationship between grounding line retreat and the non-linearity of the sea level response over time is illustrated by the colored lines in Figure 2. As alluded to above, the non-linearity in the rate of grounding line retreat is related to the bedrock topography, as demonstrated for Pine Island Glacier in Figure 3. Approximately 18% of our simulations maintain their initial grounding line position (in the case of the Bedmap2 simulations), or close to it (in the case of the modified bedrock, which has a topographic rise ∼15 km upstream) for much, if not all, of the simulations. Others experience only limited retreat, with a number of topographic rises, most notably at ∼25 km upstream of the initial grounding line position, producing a step-like pattern in retreat (Figure 3). These modest responses are similar to the first mode of retreat described by Gladstone et al. (2012), in which retreat is gradual on the order of 0.1 km/year. Gladstone et al. (2012) find a second mode of grounding line behavior characterized by rapid accelerating retreat. In their flow-line model simulations the initial grounding line retreat off the bedrock high occurs quickly and, once it reaches an uninterrupted retrograde slope, the rate of retreat can reach up to 10 km/year and be sustained for up to ∼10 years, which is similar to our most extreme results ( Figure 3). The similarities between our results and those shown in Figure 3 of Gladstone et al. (2012), despite the significant differences in model physics, indicates that topography, as well as the forcing, exerts a strong control on the temporal form of Pine Island Glacier grounding line retreat. Robel et al. (2019) demonstrate that an ensemble becomes progressively more skewed toward greater retreat when the grounding line is located on a predominantly retrograde bedrock slope, because the rate of retreat in the extreme ensemble members diverges further away from the more moderate members, which is seen here in Figure 3. The skewness in the distribution toward the high end of sea level rise is fundamentally linked to the non-linearity in the rate of grounding line retreat (Robel et al., 2019).
Missed processes in the model contribute to its structural error. For example, we do not include calving in our model-the ice front is fixed, and we impose a minimum ice thickness of 10 m. This could have implications for stability as ice shelves can provide a buttressing effect on the grounded ice sheet (Gudmundsson, 2013); although in many cases here the simulated ice shelf, across the vast majority of the area, is close to or at the minimum thickness constraint of 10 m, the buttressing effect of which is negligible. The lack of calving and ice shelf collapse, precludes any potential loss through marine ice cliff instability (DeConto & Pollard, 2016), although on the timescales of this study, it is unlikely that sufficient surface melt will occur to cause ice shelf collapse in the ASE (Kuipers Munneke et al., 2014;Trusel et al., 2015). Another source of model uncertainty is the sliding law used to determine basal shear stress-the choice of which can lead to different ice sheet responses (Brondex et al., 2017(Brondex et al., , 2019Nias et al., 2018). Given the amount of grounding line retreat experienced by the extreme simulations, we would expect interactions with neighboring West Antarctic drainage basins Feldmann & Levermann, 2015). However our model configuration has a fixed boundary so these dynamics have not been explored here.
Ocean-driven melt is likely to be a major source of uncertainty in future projections of sea level rise (Nowicki & Seroussi, 2018;Schlegel et al., 2018), for example there is uncertainty in the future ocean temperature projection, its relationship to melt rate, and how this is parameterized in BISICLES. Here, we have tested the impact of uncertainty in the melt rate obtained during the initialization of BISICLES: halving and doubling the optimal melt rate field results in a 4-cm difference in sea level contribution by 2200, when all other parameters are held at their optimal values. This is an order of magnitude less than the total spread of the distribution given in Table 1. However, we have not investigated the impact of uncertainty in the melt rate forcing anomaly added during the extended simulations. Accounting for different ocean temperature projections is likely to add considerable spread to the distribution of future sea level contribution (Holland et al., 2019), whereas the precise form of the ocean melt parameterization is likely to be less influential (Favier et al., 2019).
For the PSK group of ice streams, retreat occurs in all ensemble members, and there is less ambiguity in the future outlook, compared with Pine Island and Thwaites glaciers (Figure 2). Smith Glacier in particular has seen rapid retreat over the last two decades, and although there has been a recent slowdown in the retreat, Scheuchl et al. (2016) predict that it will continue unabated in the coming years, and they have attributed the recent stabilization to a locally prograde slope. The uncertainty in future sea level rise from the ASE lies in the vast range of responses exhibited in Pine Island and Thwaites glaciers, not the PSK group, which are relatively consistent with one another and where the potential for retreat is much more limited by the topographic constraints.
Simulations that behave similarly at the beginning of the 200 year simulations do not necessarily follow a similar trajectory ( Figure S5). This demonstrates that it is essential to use process-based models, which can predict the changing evolution of the ice sheet, instead of extrapolation methods when making projections; as also found by others (Kopp et al., 2017;Ritz et al., 2015).

Conclusion
Here we have attempted to constrain and quantify the uncertainty in sea level rise from the ASE using BISI-CLES, a high-resolution ice-flow model capable of capturing grounding line dynamics. Using present-day (2010-2015) observations we calibrated the perturbed parameter ensemble of Nias et al. (2016) by scoring each member on its ability to match the magnitude and spatial pattern of surface elevation change. Based on the resulting posterior distribution of sea level change rates, the extreme ensemble members, which matched poorly with observations, were discarded. Simulations that start out in agreement with the present day can end up contributing more than 42 cm (5% probability) by 2200, although the modal estimate is 12 cm. The long high-end tail of the sea level distribution becomes more exaggerated over time due to extreme members exhibiting a super-linear sea level response and is mirrored in the divergence in grounding line response between ensemble members. Our results only reflect uncertainty in the ice dynamics, and the range of potential outcomes would be greater if the uncertainty in the projected forcing is also included. Overall the uncertainty in the response of this particularly dynamic region of Antarctica is a major challenge in provisions for mitigating the impact of sea level rise, and at this time extreme future outcomes cannot be excluded.