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, Reference McMillan2014; Shepherd and others, Reference Shepherd2019; Smith and others, Reference Smith2020). Between 2002 and 2013, the ASE sector contributed an estimated 3.2 ± 0.1 mm to total sea-level rise (SLR) (Sutterley and others, Reference Sutterley2014), 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, Reference Shepherd2019). The loss has been identified as being clearly dynamic in origin rather than being driven by surface mass-balance (SMB) variability (Shepherd and others, Reference Shepherd2019). 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, Reference Shepherd2019). 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−1 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, Reference Hughes1973; Schoof, Reference Schoof2007; Joughin and others, Reference Joughin, Smith and Medley2014). 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, Reference Rosier2021). MISI in the ASE sector and subsequent collapse could raise global sea levels by 2.3 m (Martin and others, Reference Martin, Cornford and Payne2019).
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−1 in an ensemble experiment calibrated with observations (Nias and others, Reference Nias, Cornford, Edwards, Gourmelen and Payne2019); total contributions by 2100 of 2.0–4.5 cm by 2100 under various ocean forcing scenarios and initialisation conditions (Alevropoulos-Borrill and others, Reference Alevropoulos-Borrill, Nias, Payne, Golledge and Bingham2020); 1.5–4 cm by 2100 in an experiment that varied future ice-shelf melt rates and meteoric accumulation (Cornford and others, Reference Cornford2015); a median rise by 2100 relative to 2000 of 2 cm using different ocean forcings (Levermann and others, Reference Levermann2020); 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, Reference Wernecke, Edwards, Nias, Holden and Edwards2020). To put these ASE SLRs into context, current observed 2006–2018 rates of global SLR are 3.69 mm a−1 (3.21–4.18 mm a−1) (Fox-Kemper and others, Reference Fox-Kemper2021) 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 (Reference Aschwanden, Bartholomaus, Brinkerhoff and Truffer2021), and Durand and others (Reference Durand2022) 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, Reference Nias, Cornford, Edwards, Gourmelen and Payne2019; Aschwanden and others, Reference Aschwanden, Bartholomaus, Brinkerhoff and Truffer2021). 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, Reference Nias, Cornford, Edwards, Gourmelen and Payne2019; Wernecke and others, Reference Wernecke, Edwards, Nias, Holden and Edwards2020), WAIS contributions to SLR based on a large ensemble calibrated with scores based on modern and reconstructed geological observations (Pollard and others, Reference Pollard, Chang, Haran, Applegate and DeConto2016), a Greenland ice-sheet model calibrated using profiles of surface elevation (Chang and others, Reference Chang, Applegate, Haran and Keller2014) 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, Reference Gilford2020). Felikson and others (Reference Felikson2022) 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 (Reference Aschwanden and Brinkerhoff2022) 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 (Reference Nias, Cornford, Edwards, Gourmelen and Payne2019) 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, Reference Cornford2013) 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 (Reference Schoof and Hindmarsh2010). 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, Reference Morlighem2020), and a three-dimensional temperature field based on the Jet Propulsion Laboratory Ice Sheet System Model submission to initMIP-Antarctica (Seroussi and others, Reference Seroussi2019). 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, Reference Agosta2019) (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, Reference Mouginot, Scheuchl and Rignot2017c), 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, Reference Surawy-Stepney, Hogg, Cornford and Davison2023).
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 = φ−3) and a basal traction coefficient (β 2) are defined by solving an inverse problem (Cornford and others, Reference Cornford2015). 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 = β 2|u|), that produce a good match between modelled and observed ice speeds (Cornford and others, Reference Cornford2015).
The optimal combination of β 2 and φ, subject to 0 < φ < 1, is given by minimising the function
where u mi and u o are the modelled and observed ice speeds, $\Vert .\Vert$
2 means $\int \vert.\vert ^2{\rm d}\Omega$
, and σ 2 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, Reference Hansen2007). Note that simultaneous estimation of β 2 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, Reference Cuffey and Paterson2010) 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, Reference Rignot, Mouginot and Scheuchl2011; Mouginot and others, Reference Mouginot, Scheuchl and Rignot2012, Reference Mouginot, Rignot, Scheuchl and Millan2017a; Rignot aand others, Reference Rignot, Mouginot and Scheuchl2017); 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, β 2 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 β 2 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 β 2 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, Reference Joughin, Smith and Schoof2019):
where
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 0, 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, Reference De Rydt, Reese, Paolo and Gudmundsson2021).
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−1), and five values of u 0 in the friction law (20, 75, 150, 300 and 600 m a−1) – 25 runs. The set of u 0 values are centred on and expand upon the value of 300 m a−1 found to perform well for Pine Island Glacier (Joughin and others, Reference Joughin, Smith and Schoof2019). We then ran model simulations in which we scaled each of the optimal $C_{u_0}$
and φ coefficients for each combination of u 0 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−1), u 0 = (20, 75, 150, 300 and 600 m a−1), 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−1 everywhere, we used the optimum $C_{u_0}$
and φ fields, five values of ∂h/∂t(Ωf)(0, − 2, − 5, − 10 and −15 m a−1), and u 0 values of 20, 75, 150, 300, 600, 1200 and 2400 m a−1.
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−1 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−1 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, Reference Silverman1986). 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).
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 (Reference Nias, Cornford, Edwards, Gourmelen and Payne2019), 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, Reference Edwards2014; Pollard and others, Reference Pollard, Chang, Haran, Applegate and DeConto2016). 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:
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:
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, Reference Edwards2014). 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 0 = 20 m a−1, and ∂h/∂t(Ωf) = −2 m a−1.
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 0 = 20 m a−1, and ∂h/∂t(Ωf) = −5 m a−1.
As described in Nias and others (Reference Nias, Cornford, Edwards, Gourmelen and Payne2019), $\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, Reference Mouginot, Rignot, Scheuchl and Millan2017a).
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, Reference Nias, Cornford, Edwards, Gourmelen and Payne2019; Wernecke and others, Reference Wernecke, Edwards, Nias, Holden and Edwards2020; Felikson and others, Reference Felikson2022). 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:
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 (Reference Shepherd2019). 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 0 = 20 m a−1 and ∂h/∂t(Ωf) = −5 m a−1. 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, Reference Mouginot, Scheuchl and Rignot2017b) (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 0 = 20 ma−1 and ∂h/∂t(Ωf) = −5 m a−1.
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 (Reference Zwally, Giovinetto, Beckley and Saba2012), e.g. Fig. 7). Both observed and modelled volume losses were converted to SLE by equating 103 km3 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−1, the other five a rate of 15 m a−1. 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−1.
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 0, ∂h/∂t(Ωf) and SMB
Figure 12 shows the 2050 SLE of ice loss according to u 0 and ∂h/∂t(Ωf) for the MAR SMB and SMB =0.3 m a−1 fields. Removing the excess SMB that MAR provides compared with 0.3 m a−1 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 0 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−1, and the black sections show the additional SLE after the excess (MAR - 0.3 m a−1) accumulation has been removed from the MAR SMB simulations. (Note that the simulations for ∂h/∂t(Ωf) = −15 m a−1 with MAR SMB and ∂h/∂t(Ωf) = −10 m a−1 with SMB = 0.3 m a−1, failed to complete.)
Figure 12 also reveals that for ∂h/∂t(Ωf) ≥ −5 m a−1, 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−1 (Fig. 12). This shelf thinning rate is greater than observed thickness changes of −19.4 m/decade for the region (Paolo and others, Reference Paolo, Fricker and Padman2015). 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, Reference Alley, Scambos, Alley and Holschuh2019; Lhermitte and others, Reference Lhermitte2020) 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 0 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, Reference Jourdain, Mathiot, Burgard, Caillet and Kittel2022). A doubling of shelf-thinning rates from 5 to 10 m a−1 in these simulations increases 2050 SLE ice loss by about 30%.
The correlation between the magnitude of u 0 in the regularised Coulomb friction law and the total SLE switches sign between ∂h/∂t(Ωf) = −5 m a−1 and ∂h/∂t(Ωf) = −10 m a−1 (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 0 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 0). 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 0, 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−1 everywhere, corresponding to a reduction in total accumulated mass of 61 Gt a−1 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, Reference Donat-Magnin2021) 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, Reference Felikson2022). 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−1 highlighted that ice loss by 2050 was more sensitive to ice-shelf thinning or damage, than to the chosen value of u 0 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.
Data
The surface elevation change observational datasets are available at https://zenodo.org/record/8117577 (doi: 10.5281/zenodo.8117577). A summary file of the ensemble simulated annual SLE values, and annual NetCDF files of land-ice thickness, and the u and v components of velocity at 1 km spatial resolution are available at https://zenodo.org/record/8120460 (doi: 10.5281/zenodo.8120460).
Acknowledgements
This publication was supported by PROTECT. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 869304, PROTECT contribution number 71. Support for this work was provided through the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy (DOE), Office of Science, Biological and Environmental Research and Advanced Scientific Computing Research programs, as a part of the ProSPect SciDAC Partnership. Work at Berkeley Lab was supported by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award ASCR-ERCAPm1041.
Author contributions
Suzanne Bevan designed and ran the model simulations, analysed the results and wrote most of the text. Stephen Cornford contributed to data analysis and writing of the text. Lin Gilbert, Inés Otosaka and Trystan Surawy-Stepney contributed data. Stephen Cornford and Daniel Martin supported Suzanne Bevan in use of the BISICLES model and Berkeley Lab supercomputing facilities.
Open access
