Amundsen Sea Embayment ice-sheet mass-loss predictions to 2050 calibrated using observations of velocity and elevation change

Abstract

Mass loss from the Amundsen Sea Embayment of the West Antarctic Ice Sheet is a major contributor to global sea-level rise (SLR) and has been increasing over recent decades. Predictions of future SLR are increasingly modelled using ensembles of simulations within which model parameters and external forcings are varied within credible ranges. Accurately reporting the uncertainty associated with these predictions is crucial in enabling effective planning for, and construction of defences against, rising sea levels. Calibrating model simulations against current observations of ice-sheet behaviour enables the uncertainty to be reduced. Here we calibrate an ensemble of BISICLES ice-sheet model simulations of ice loss from the Amundsen Sea Embayment using remotely sensed observations of surface elevation and ice speed. Each calibration type is shown to be capable of reducing the 90% credibility bounds of predicted contributions to SLR by 34 and 43% respectively.

1. Introduction

Current Antarctic Ice Sheet thinning is concentrated in the Amundsen Sea Embayment (ASE) sector of the West Antarctic Ice Sheet (WAIS), and in Wilkes Basin, East Antarctica (McMillan and others, 2014; Shepherd and others, 2019; Smith and others, 2020). Between 2002 and 2013, the ASE sector contributed an estimated 3.2 ± 0.1 mm to total sea-level rise (SLR) (Sutterley and others, 2014), with ice loss from the sector increasing by a factor of 5 from the decade ending in 2002 to the decade ending in 2017 (Shepherd and others, 2019). The loss has been identified as being clearly dynamic in origin rather than being driven by surface mass-balance (SMB) variability (Shepherd and others, 2019). The areas that are thinning and losing mass are spreading inland, and an estimated 59% of the ASE sector was out of balance by 2016 (Shepherd and others, 2019). By this time, the two largest glaciers draining the ASE sector, Thwaites and Pine Island Glaciers, were each losing mass at rates of 55 ± 4 and 76 ± 6 Gt a⁻¹ of ice, respectively. Much of the region's ice is grounded below sea level, with the ice–bed interface sloping inland away from the grounding line, making it potentially vulnerable to marine ice-sheet instability (MISI), whereby loss of buttressing at the periphery leads to rapid and runaway grounding line retreat (Hughes, 1973; Schoof, 2007; Joughin and others, 2014). Early warning indicators show that two out of three basal melt-rate tipping points for MISI may already have been crossed on Pine Island Glacier, with the third tipping point, one which would lead to a total retreat of the glacier, being an increase in basal melt rate corresponding to an increase in ocean temperatures of 1.2°C (Rosier and others, 2021). MISI in the ASE sector and subsequent collapse could raise global sea levels by 2.3 m (Martin and others, 2019).

Various studies have used stand-alone ice-sheet models to estimate contributions to future SLR from the ASE region. The modelled rates of 21st century SLR include: a modal rate value of 0.3 mm a⁻¹ in an ensemble experiment calibrated with observations (Nias and others, 2019); total contributions by 2100 of 2.0–4.5 cm by 2100 under various ocean forcing scenarios and initialisation conditions (Alevropoulos-Borrill and others, 2020); 1.5–4 cm by 2100 in an experiment that varied future ice-shelf melt rates and meteoric accumulation (Cornford and others, 2015); a median rise by 2100 relative to 2000 of 2 cm using different ocean forcings (Levermann and others, 2020); and a median 50 year SLR of 18.9 mm (with 90% confidence limits of 13.9–24.8 mm) in a study testing ensemble calibration methods (Wernecke and others, 2020). To put these ASE SLRs into context, current observed 2006–2018 rates of global SLR are 3.69 mm a⁻¹ (3.21–4.18 mm a⁻¹) (Fox-Kemper and others, 2021) meaning that estimates of ASE near-future contributions are of the order of 10% of global rates.

In considering global and local impacts of and adaptations to rising sea levels, Aschwanden and others (2021), and Durand and others (2022) emphasise the need for a realistic quantification of the uncertainty associated with SLR predictions. Identifying the sources of uncertainty will, further, allow efforts to reduce them to be directed more effectively. By using a probability density function (pdf) to represent the outcome, for example SLR, of ensembles of ice-sheet model runs, the percentage probability bounds can be identified within which predicted outcomes may be expected to lie. The ensembles may represent different initial states, different parameters, different forcing or different boundary conditions. The resulting model pdfs, or priors, should ideally be conditioned or calibrated with observations to produce posterior pdfs (Nias and others, 2019; Aschwanden and others, 2021). Examples of using calibrated pdfs of SLR based on large ensembles of ice-sheet models calibrated with observations to quantify uncertainty include contributions to 2100 SLR from the ASE calibrated using observed rates of surface elevation change (Nias and others, 2019; Wernecke and others, 2020), WAIS contributions to SLR based on a large ensemble calibrated with scores based on modern and reconstructed geological observations (Pollard and others, 2016), a Greenland ice-sheet model calibrated using profiles of surface elevation (Chang and others, 2014) and the use of little ice age constraints to condition an emulated response of the Antarctic Ice Sheet to future high emission scenarios (Gilford and others, 2020). Felikson and others (2022) used observations of thickness change, mass change and velocity change to compare the impact of calibration type on the posterior distribution of an ensemble of model simulations for the Greenland ice sheet, and Aschwanden and Brinkerhoff (2022) explored the use of a two-stage calibration method for ensemble simulations, also of the Greenland ice sheet.

In this report, we consider the uncertainties in predicted sea-level equivalent (SLE) of ice loss from the ASE region associated with parameter specifications in the BISICLES ice-sheet model. In particular, we consider different parameter values in a regularised Coulomb friction law, in the specification of ice-shelf thinning rates, and in coefficients of viscosity and basal drag. We then use observations to calibrate a large ensemble of simulations that run to 2050. The work builds on the study by Nias and others (2019) in that we now have available an improved bedrock and thickness dataset, implement a new sliding law in the model and as well as elevation change observations we now use velocity observations to calibrate the ensemble.

2. Methods

2.1. BISICLES model

We used the BISICLES ice-flow model (Cornford and others, 2013) to simulate the evolution of the ASE sector of the Antarctic ice sheet to the year 2050. BISICLES is a finite-volume model that incorporates longitudinal and lateral stresses and a simplified treatment of vertical shear based on Schoof and Hindmarsh (2010). Adaptive mesh refinement allows the spatial resolution to be fine in the locality of grounding lines and coarser elsewhere, and to evolve with grounding line migration. Here, meshes are constructed from rectangular patches, each of which is a grid of square cells 4 km, 2 km, 1 km or 500 m across. (Fig. 1).

Figure 1. Model domain and nested mesh regions.

Initial-state datasets included bedrock and surface elevation, and ice thickness from MEaSUREs BedMachine Antarctica Version 2 (Morlighem and others, 2020), and a three-dimensional temperature field based on the Jet Propulsion Laboratory Ice Sheet System Model submission to initMIP-Antarctica (Seroussi and others, 2019). Surface accumulation rates were held constant throughout, at the mean of the monthly values from 1981 to 2018 from the Modéle Atmosphérique Régional (MAR) version 3.10 regional climate model (Agosta and others, 2019) (Fig. 2).

Figure 2. (a) Mean and (b) standard deviation of the monthly 1981–2018 Modéle Atmosphérique Régional (MAR) version 3.10 surface mass-balance datasets.

We denote the part of the domain containing floating ice at time t by Ω_f(t). Ice-shelf melt rates were imposed indirectly by setting an ice-shelf thickness change rate ∂h/∂t (Ω_f) < 0, allowing direct comparison with observations, and avoiding the need to invent depth-dependent melt-rate parameterisations. The imposition of ∂h/∂t(Ω_f) < 0 could represent either ocean melt in excess of that needed to counter ice flow, or a weakening as a result of ice-shelf damage.

In addition, we specified ice fronts based on observations (Fig. 3). The complete ice front corresponds to the MEaSUREs V2 coastline (Mouginot and others, 2017c), Pine Island Glacier fronts were adjusted to the observed locations in September 2017, October 2018 and February 2020 coinciding with major calving events; and Thwaites Glacier fronts were defined at monthly intervals based on a convolutional neural network applied to Sentinel-1 SAR backscatter images (Surawy-Stepney and others, 2023).

Figure 3. Ice front positions based on satellite observations and the MEaSUREs V2 coastline.

The formulation of the model requires that a stiffening coefficient (φ, equivalent to an enhancement cofficient E = φ⁻³) and a basal traction coefficient (β ²) are defined by solving an inverse problem (Cornford and others, 2015). The stiffening coefficient is applied to the vertically integrated effective viscosity to compensate for uncertainties in ice temperature, uncertainties in the rate factor in Glen's flow law and large-scale damage. The goal of the inverse problem is to create smooth continuous fields of φ and effective drag (τ ^b = β ²|u|), that produce a good match between modelled and observed ice speeds (Cornford and others, 2015).

The optimal combination of β ² and φ, subject to 0 < φ < 1, is given by minimising the function

(1)$$J = \Vert ( u_{\rm mi} - u_{\rm ob}) \Vert ^2 + {\sigma^2\over \lambda_{\beta^2}^2}\Vert \nabla \beta^2\Vert ^2 + {\sigma^2\over \lambda_\varphi^2}\Vert \nabla \varphi\Vert ^2,\; $$

where u _mi and u _o are the modelled and observed ice speeds, $\Vert .\Vert$² means $\int \vert.\vert ^2{\rm d}\Omega$, and σ ² represents the combined error of modelled and observed velocity. $\lambda _{\beta ^2}^2$ and $\lambda _\varphi ^2$ are regularisation parameters. Their values are chosen through L-curve analysis (Hansen, 2007). Note that simultaneous estimation of β ² and φ over grounded ice is underdetermined, but is necessary if φ is to be continuous across the grounding line. The regularisation terms, together with the constraint 0 < φ < 1, mitigate against this to some extent (e.g. by prohibiting regions of φ > >1 that would arise in regions of near-zero rate-of-strain) but ultimately mean that an initial guess for both is required and important. Our initial guess assumes that viscous stresses can be neglected over the majority of the ice sheet, so that basal stress τ ^b equals gravitational driving stress τ ^d and that the relationship between temperature and rate factor (Cuffey and Paterson, 2010) holds where viscous stresses are unimportant.

The observed ice speeds used for the inversion were the MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2 (Rignot and others, 2011; Mouginot and others, 2012, 2017a; Rignot aand others, 2017); those for the ASE region are based mostly on data acquired in 2007, hence we assign 2007 as the start date of our simulations.

The complete model initial state comprises ice thickness h, β ² and φ, which should be consistent in the sense that the initial u and ∂h/∂t fields are consistent with observations. Simply taking h to be the BedMachine V2 ice thickness results in (subtle) changes to the ice geometry that lead over a few years to substantial reduction in ice flow. To address this, the grounded ice thickness was relaxed over 10 years (see Fig. 4) with the ice-shelf thickness held constant and β ² and φ periodically recomputed.

Figure 4. (a) Surface elevation from BedMachine V2. (b) Difference in surface elevation from BedMachine V2 after 10 years of model relaxation.

Using the optimal φ and β ² coefficients and the relaxed geometry, we ran the model to the year 2050, defining the start year to be 2007 and imposing the calving-front retreats described above between 2016 and 2021. For these simulations, we switched to a regularised Coulomb friction law (Joughin and others, 2019):

(2)$$\tau^b = -C_{u_0}\left({\left\vert \bf\,{u}\right\vert \over \left\vert \bf\,{u}\right\vert + u_0}\right)^{1/3},\; $$

where

(3)$$C_{u_0} = \beta^2 u_{\rm mi}^{2/3}\left(u_{\rm mi} + u_0 \right)^{1/3},\; $$

and u _mi is the initial model speed, i.e the speed computed when solving the inverse problem.

This formulation of the friction law allows friction to transition towards a plastic Coulomb law at sliding velocities  ≫ u ₀, i.e. for fast-flowing ice streams and glaciers with soft deformable sediments at the bed, where a plastic bed is required to realistically model recent observed increases in flow, for example on Pine Island Glacier (De Rydt and others, 2021).

2.2. Ensemble design

An ensemble of model simulations was designed to explore the range in contributions to SLR originating from uncertainties in the basal traction and ice viscosity fields, and rate of ice shelf thinning.

The first set combined the optimal fields of $C_{u_0}$ and φ derived from the inversion, five values of ∂h/∂t(Ω_f) (0, −2, −5, −10 and −15  m a⁻¹), and five values of u ₀ in the friction law (20, 75, 150, 300 and 600 m a⁻¹) – 25 runs. The set of u ₀ values are centred on and expand upon the value of 300 m a⁻¹ found to perform well for Pine Island Glacier (Joughin and others, 2019). We then ran model simulations in which we scaled each of the optimal $C_{u_0}$ and φ coefficients for each combination of u ₀ and ∂h/∂t(Ω_f) by factors of 0.9 and 1.1, generating a further 100 parameter sets and outcomes. Using a maximin Latin hypercube sampling, a further set of 100 simulations was drawn, subsampling the four-dimensional parameter space defined by all of the above values of ∂h/∂t(Ω_f) = (0, − 2, − 5, − 10 and −15 m a⁻¹), u ₀ = (20,  75,  150,  300 and 600 m a⁻¹), and $C_{u_0}$ and φ coefficient factors between 0.9 and 1.1. The complete ensemble size consisted of 225 simulations, minus 12 that did not complete.

Finally, we ran a set of simulations with a lower overall SMB than MAR in order to reveal the effect of SMB forcing on ice-sheet dynamics. For these simulations, the SMB was set to 0.3 m a⁻¹ everywhere, we used the optimum $C_{u_0}$ and φ fields, five values of ∂h/∂t(Ω_f)(0, − 2, − 5, − 10 and −15 m a⁻¹), and u ₀ values of 20, 75, 150, 300, 600, 1200 and 2400 m a⁻¹.

This final set of simulations with lower SMB was not included in the ensemble used in model calibration, but is presented simply to investigate the impact of varying SMB. The 0.3 m a⁻¹ SMB, integrated over the model domain, is about 25% lower than the MAR SMB. If accumulating this difference over the full model run time and removing it from the final SLE generated by the MAR SMB simulations produced the same SLE as the SMB = 0.3 m a⁻¹ runs, then we could conclude that an increased SMB has no effect on ice-sheet dynamics.

2.3. Model calibration

Our ensemble of simulations were used to calculate the SLE of ice loss by 2050 relative to 2007, from all ice above floatation within the Thwaites and Pine Island Glacier drainage basins. The ensemble SLEs were used to generate a smoothed pdf for SLE using a Gaussian kernel density estimator with a Silverman rule-of-thumb bandwidth (Silverman, 1986). This pdf, or prior, is denoted P(SLE).

According to Bayes’ theorem a posterior probability distribution P(SLE|O) can be found by weighting or conditioning the prior with a likelihood function P(O|SLE).

(4)$$P\left(SLE \vert O \right)\propto P\left(SLE\right)P\left(O \vert SLE \right).$$

By using a kernel density estimator, we are able to estimate the modes and percentiles of the distributions of prior and posterior SLEs.

Following Nias and others (2019), we substituted for P(O|SLE) a set of weights based on the difference between the model output and observations (see also Edwards and others, 2014; Pollard and others, 2016). First we scored the model simulations against an observed spatial distribution of surface elevation change rate, ∂h/∂t, averaged over 2007–2017 (section 2.4). For each model output $m^{j}_{i}$ at location i from simulation j, we used the corresponding observation o _i to calculate a simulation misfit M _j:

(5)$$M_j = {1\over n}\sum_{i}{\left(m^{\,j}_{i} - o_{i}\right)^2\over \sigma^{2}_{i}},\; $$

where n is the total number of locations or pixels involved in the summation.

A likelihood score for each simulation, S _j, was calculated by:

(6)$$S_j = \exp\left[-{1\over 2} M_j \right].$$

We also scored the simulations against observed surface velocity magnitude for 2007 and for 2020 (section 2.4). In this case, the summation in Eqn (5) is over the two years as well as over space. The area used for both calibrations was restricted to the Thwaites and Pine Island Glacier drainage basins.

Equation (6) assumes that the discrepancies are identically distributed in space, and independent (Edwards and others, 2014). In order to maximise the independence, we down-scaled the spatial resolution of model–observation discrepancies to reduce autocorrelation. Using the R (R Core Team, 2022) function variogram from the gstat package, we plotted the spatial autocorrelation of model–observation discrepancies in ∂h/∂t and u (Fig. 5).

Figure 5. Semi-variograms for model–observation discrepancy in directions 0, 45, 90 and 135 degrees. (a) Mean ∂h/∂t between 2007 and 2017, (b) year 2020 u. The plots are based on the simulation with u ₀ = 20 m a⁻¹, and ∂h/∂t(Ω_f) = −2 m a⁻¹.

For ∂h/∂t the range, the point beyond which the semi-variance no longer increases was ~100 km. The sill was less well determined for u but the semi-variance increases most steeply over the first 25 km (Fig. 5b). Therefore, we resampled the ∂h/∂t and u discrepancies to 100 and 25 km, respectively (Fig. 6).

Figure 6. Example simulation-minus-observation discrepancy maps used for scoring. (a) Mean ∂h/∂t discrepancy at 100 km resolution, (b) 2020 u discrepancy at 25 km resolution. The plots are based on the simulation with u ₀ = 20 m a⁻¹, and ∂h/∂t(Ω_f) = −5 m a⁻¹.

As described in Nias and others (2019), $\sigma ^2_{i}$ in Eqn (5) is the discrepancy variance and it represents how well we expect the model to agree with an observation. $\sigma ^2_{i}$ needs to include both model and observational uncertainty. To set the uncertainty in observations of the 11-year mean ∂h/∂t we relied on the estimates of uncertainty accompanying the supplied data, these estimates are space and time dependent and were summed in quadrature for the period 2007–17. For surface speed, we also used the errors supplied with the data which are dependent on the sensor and the degree of data stacking within the measurement (Mouginot and others, 2017a).

Model uncertainty is harder to put a number on, and it is common to choose a model error that is some multiple of observation error (e.g Nias and others, 2019; Wernecke and others, 2020; Felikson and others, 2022). Here, we decided to base the model error on the magnitude of the observation on the basis that observation error is mostly a function of observing technique and available observations, and not necessarily relevant to the magnitude of error we can tolerate in the simulations. However, we can expect and tolerate a larger model error where velocities are high and ∂h/∂t is large. The goal is to avoid concentrating the weights on a small number of models, and yet to provide some discrimination between simulations: a very large discrepancy variance would result in a posterior distribution no narrower than the prior. For each calibration type, we tested model errors equivalent to 1, 5, 10 and 50% of the observed values.

Having calculated the set of likelihood scores, the weight for each ensemble member j is given by the normalised likelihood score:

(7)$$W_j = {S_j\over \sum_{\,j}{S_j}},\; $$

which is P(O|SLE) in Eqn (4).

2.4. Observations of surface elevation change and velocity

Observed surface elevation change rates used in the calibration consist of annual ∂h/∂t rates computed over 3-year periods using data from ERS-1/2, ENVISAT and CryoSat-2. We used data from 2007 to 2017. The method used for the derivation of these elevation change rates is described in Shepherd and others (2019). The 2007–17 mean detected ∂h/∂t is shown in Figure 7a.

Figure 7. (a) Observed 2007–2017 mean surface elevation change rates. (b) Example simulated 2007–2017 mean elevation change rates, for u ₀ = 20 m a⁻¹ and ∂h/∂t(Ω_f) = −5 m a⁻¹. Values are masked by the drainage basin boundaries for Pine Island and Thwaites Glaciers.

As well as the MEaSUREs Version 2 ice velocity, we also used the 2019–2020 velocity from the MEaSUREs Version 1, 1 km annual velocities dataset (Mouginot and others, 2017b) (Fig. 8a), to calibrate the model response. The MEaSUREs Version 2 velocities were used to infer the initial optimum $C_{u_0}$ and φ fields; we would therefore expect these simulations to have a low mismatch with these early velocities meaning we are mainly testing the non-optimal $C_{u_0}$ and φ simulations. By including a second, later in time velocity, we are also testing the ability of all simulations to match the acceleration of the ice.

Figure 8. (a) MEaSUREs V1 2019–2020 surface velocities. (b) Example of 2020 simulated velocities, for u ₀ = 20 ma⁻¹ and ∂h/∂t(Ω_f) = −5 m a⁻¹.

For the purposes of calculating modelled and observed losses of volume of ice above floatation, the region was restricted to the combined Thwaites and Pine Island Glacier drainage basins (basins 21 and 22 in Zwally and others (2012), e.g. Fig. 7). Both observed and modelled volume losses were converted to SLE by equating 10³ km³ of ice to 2.5 mm of sea level.

3. Results

Examples of modelled ∂h/∂t (mean 2007–17 rate) and u (2020) are shown in Figures 7b and 8b, respectively, where they can be compared with observed values.

The ensemble of simulations described in section 2.2 resulted in losses/gains of ice ranging from 63.7 to −4.4 mm SLE by 2050 (Fig. 9). The histogram of SLEs and the derived pdf is shown in all panels of Figure 11. The median SLE of the ensemble of runs was 16.2 mm. The 5th and 95th percentiles of the pdf are −0.2 and 39.1 mm SLE, respectively, and the modal value is 10.7 mm. The observation-based SLE of 14.5 mm by 2050 is based on a linear extrapolation of the mean rate of observed volume change between 2007 and 2017 (black dashed line in Fig. 11).

Figure 9. Time evolution of SLE ice loss for the large ensemble described in section 2.2. The observations are based on observed dh/dt.

3.1. Calibration using ∂h/∂t and u

We used pair plots of SLE, and S _j for ∂h/∂t and u to trial the different scaling factors for model error. Using a 10% scaling produced reasonable results in that scores are non-zero over a large number of simulations (Figs 10a,b) and the two scoring methods are strongly correlated with a Pearson's correlation coefficient of 0.90, significant at the 99% level (Figs 10c). Figure 10 also highlights that the simulations with optimum $C_{u_0}$ and φ score highly. The latter being true indicates that these simulations are capturing well both the rate of mass loss and the acceleration of ice during the calibration period.

Figure 10. Pair plots of the full ensemble, axis labels are given by the text in the diagonal panels. (a) S _j (vel) distribution as a function of SLE. (b) S _j (dh/dt) distribution as a function of SLE. (c) Relationship between S _j (vel) and S _j (dh/dt). SLE is the 2050 SLE of mass loss. S _j (vel) is the score for each simulation calculated using velocity and a model error equivalent to 10% of the observed value of velocity. S _j (dh/dt) is the score using ∂h/∂t and a model error also equivalent to 10% of the observed value of ∂h/∂t. CC = 0.9 is the Pearson's correlation coefficient between S _j (dh/dt) and S _j (vel).

All ten of the highest scoring simulations scored using ∂h/∂t are with optimum $C_{u_0}$ and φ values, of these five have a shelf thinning rate of 10 m a⁻¹, the other five a rate of 15 m a⁻¹. Of the top ten highest scoring simulations scored by u, six have optimum $C_{u_0}$ and φ values but have shelf thinning rates ranging from 0 to 10 m a⁻¹.

Weighting the prior pdfs using likelihood scores based on observations of ∂h/∂t and u and with model errors set equal to 10% of the observed values narrows the 90% credibility bands of likely 2050 SLE of ice loss to 1.7–27.5 and 2.8–25.3 mm (Fig. 11a); with median (mode) values of 13.2 and 13.0 mm (12.1 and 12.1 mm) (Table 1). Increasing (decreasing) the allowance for model error broadens (narrows) the credibility range (Figs 11b,c). Note that in Figure 11b, the calibrated curve with model error equal to 1% of observed value is barely distinguishable from the 10% curve. The calibration process has the effect of suppressing many of the simulations that resulted in extreme high SLEs.

Figure 11. Normalised histogram and prior probability density functions of the full ensemble of simulations. The vertical black dashed line shows the mean 2007–17 rate of observed SLE volume loss extrapolated to 2050. (a) With ∂h/∂t and u calibrations for model error equal to 10% of the observational error. (b) Calibrated SLE pdfs for the ∂h/∂t calibration method with varying model error. (c) As for (b) but for the u calibration.

Table 1. Ensemble statistics for prior and posterior (calibrated) SLE of ice loss shown in Figure 11a

3.2. Sensitivity of total SLE to the combined effect of varying u ₀, ∂h/∂t(Ω_f) and SMB

Figure 12 shows the 2050 SLE of ice loss according to u ₀ and ∂h/∂t(Ω_f) for the MAR SMB and SMB =0.3 m a⁻¹ fields. Removing the excess SMB that MAR provides compared with 0.3 m a⁻¹ everywhere, an excess equivalent to 6.7 mm SLE by 2050, shows that the dynamic forcing (shown in black in Fig. 12) provided by this excess SMB is of the order of 10% of total discharge.

Figure 12. Total sea-level equivalent of loss of ice above floatation for varying u ₀ values, grouped by the imposed rate of thinning of floating ice (∂h/∂t(Ω_f)). The grey bars show the SLE using the mean MAR SMB, the grey-plus-coloured bars show the SLE for SMB = 0.3 m a⁻¹, and the black sections show the additional SLE after the excess (MAR - 0.3 m a⁻¹) accumulation has been removed from the MAR SMB simulations. (Note that the simulations for ∂h/∂t(Ω_f) = −15 m a⁻¹ with MAR SMB and ∂h/∂t(Ω_f) = −10 m a⁻¹ with SMB = 0.3 m a⁻¹, failed to complete.)

Figure 12 also reveals that for ∂h/∂t(Ω_f) ≥ −5 m a⁻¹, increasing u _o results in an increased total SLE by 2050. At higher thinning rates the reverse is true.

4. Discussion

4.1. Sliding law parameters, shelf thinning rates and SMB

Whilst the aim of these experiments was to demonstrate the utility of using observations to reduce uncertainty in modelling results, we can also consider which model parameter sets lead to simulations that agree best with extrapolated observations of current ASE ice loss. With the optimum $C_{u_0}$ and φ fields and the MAR SMB, the simulations that best agree with the extrapolated observed SLE of 14.5 mm by 2050 are those with ∂h/∂t(Ω_f), of between −5 and −10 m a⁻¹ (Fig. 12). This shelf thinning rate is greater than observed thickness changes of −19.4 m/decade for the region (Paolo and others, 2015). However, as indicated in section 2.1, the imposed thinning rates also act as a proxy for ice-shelf damage, weakening the back stress buttressing the inland grounded ice. Damage such as fractures and crevasses in shear zones is observed to be increasing on both Thwaites and Pine Island Glacier ice shelves (Alley and others, 2019; Lhermitte and others, 2020) and may explain why we need to impose quite high thinning rates to keep pace with extrapolated observed ice loss rates. A model that evolves damage could circumvent this requirement for artificially high melt rates.

Within the range of values tested, the specification of shelf-thinning rate has a much greater effect on ice loss than the u ₀ value in the sliding law (Fig. 12). Recent projections for 2100 under RCP8.5 indicate ice-shelf melt rates increasing by factors of 1.4–2.2 compared with present day (Jourdain and others, 2022). A doubling of shelf-thinning rates from 5 to 10 m a⁻¹ in these simulations increases 2050 SLE ice loss by about 30%.

The correlation between the magnitude of u ₀ in the regularised Coulomb friction law and the total SLE switches sign between ∂h/∂t(Ω_f) = −5 m a⁻¹ and ∂h/∂t(Ω_f) = −10 m a⁻¹ (Fig. 12). This relationship indicates that greater bed plasticity leads to greater sensitivity to loss of buttressing. At higher shelf thinning rates, greater overall volume loss is seen when u ₀ is lower; that is, when a larger portion of the ice experiences plastic, Coulomb-like drag. At lower thinning rates, the converse is true. In the higher melt-rate simulations, the ice tends to flow faster, with the speed increase limited more by the increasing basal traction in Weertman-like cases (large u ₀). At lower thinning rates, the ice tends to slow down, and in these cases the speed decrease is limited more by the decreasing basal traction in Weertman-like models. In everyday terms, the higher u ₀, the more the bed will both increase its grip on the ice as it accelerates, and decrease its grip when the ice slows.

Simulations run with an SMB of 0.3 m a⁻¹ everywhere, corresponding to a reduction in total accumulated mass of 61 Gt a⁻¹ by 2050, or about 30% compared with the MAR SMB, resulted in around a 10% decrease in the dynamic discharge of ice to the ocean. These experiments with reduced SMB suggest that the increases in SMB predicted for the sector by 2080–2100 of 30–40% (Donat-Magnin and others, 2021) may be slightly offset by increases in discharge forced by the additional accumulation.

4.2. Large ensemble

The results of the large ensemble (Fig. 9) exhibit a wide range of 2050 SLE ice loss and encompass current observations, indicating that we have adequately sampled the parameter space; a conclusion confirmed by the calibrated curves (Fig. 11). Confidence in the calibration method is strengthened by the fact that the two independently calibrated posterior distributions are in good agreement with each other in terms of the maximum likelihood SLE. This agreement between calibration types is in contrast to an experiment using an ensemble of model simulations for the Greenland ice sheet (Felikson and others, 2022). In that experiment using velocity change to calibrate the ensemble resulted in a distinctly higher SLE by 2100 compared with using thickness change, which actually indicated that the ice sheet was gaining mass. There are many differences between that study and ours, but the authors concluded that using a firn densification model to convert observed thickness change to dynamic thickness may have had an effect on the calibration. In our methodology, as the majority of the region is in dynamic imbalance, we make the assumption that observed elevation changes are dominated by dynamics and thus avoid the need for any firn modelling.

The credibility bands of the calibrated curves are affected by the chosen σ values. We have little knowledge of the model structural error but using spatially variable estimates of observational uncertainty concentrates the contribution of grid cells to the likelihood values on regions where observational uncertainties are lowest. Increasing model error to be a greater proportion of observed values broadens the confidence intervals, and vice versa. Parameter pair plots (Fig. 10) confirm that for both ∂h/∂t and u, including a model σ value equal to 10% of the observed ∂h/∂t and u results in scores that avoid heavily weighting only a small number of simulations and yet do provide a useful constraint on the ensemble. Assuming a very large permissible model error would effectively make the calibration process without value.

The full 213-member ensemble included varying $C_{u_0}$ and φ by factors ranging between 0.9 and 1.1 to explore underdetermination in the initialisation. The calibration process gives low credibility to the simulations that result in extreme SLEs of ice loss. At the low end, six of the seven simulations resulting in a negative SLE were a result of setting both $C_{u_0}$ and φ to 1.1 times the optimum value obtained by the inverse problem. Similarly, at the high end, six of the seven highest SLE simulations were a result of setting both $C_{u_0}$ and φ to 0.9 times the optimum value. Thus, together with the result that the highest scoring simulations have optimum C and φ, we conclude that the calibration process adds confidence to the solution obtained to the inverse problem.

The calibrated curves assign a higher likelihood to 2050 SLE values lower than those indicated by an extrapolation of observations. This discrepancy could be a result of imposing calving front retreats for years that overlap with the calibration periods: during these periods, the simulations can match the observed rates of loss through a combination of shelf thinning and loss of shelf area. Beyond 2021, we impose no calving-front retreat and back-stress reduction is then achieved only by shelf thinning, resulting in a lower 2050 SLE compared with extrapolating observations. Alternatively, a simple extrapolation of observed mass loss cannot take into account any changes in grounding line retreat rates owing to upstream variations in bed geometry. In addition, the model is relaxing towards a stable equilibrium in which losses decay over time, whereas extrapolating observations assumes that the current imbalance is preserved.

5. Conclusions

An ensemble of decadal-scale simulations of ice-sheet evolution in the ASE region of WAIS was created using the BISICLES model, with a regularised Coulomb friction law, and with sliding and rheology parameters defined via the solution of the inverse problem.

Reproducing current observed mass loss rates requires the inclusion of some form of reduction of ice-shelf buttressing. In the experiments presented in this study, that loss of buttressing was supplied by a combination of imposed ice-shelf thinning and, in the early years, a short period of calving front retreats based on observations. Simulations with shelf thinning rates varying from 0 to 15 m a⁻¹ highlighted that ice loss by 2050 was more sensitive to ice-shelf thinning or damage, than to the chosen value of u ₀ in the sliding law.

Simulations also revealed that SMB changes have an impact on ice flow dynamics that may partially offset any predicted future increases in SMB.

An expanded ensemble including up to 10% variation in sliding and viscosity parameters was then calibrated using spatial fields of misfits between model output and remotely sensed observations of changes in surface elevation and ice speed. Calibration by u and by ∂h/∂t produced consistent posterior likelihood distributions, and reduced uncertainty compared with the prior distribution by 43 and 34%, respectively.