Derivation of bedrock topography measurement requirements for the reduction of uncertainty in ice-sheet model projections of Thwaites Glacier

Determining the future evolution of the Antarctic Ice Sheet is critical for understanding and narrowing the large existing uncertainties in century-scale global mean sea-level-rise (SLR) projections. One of the most significant glaciers and ice streams in Antarctica, Thwaites Glacier, is at risk of destabilization and, if destabilized, has the potential to be the largest regional-scale contributor of SLR on Earth. This is because Thwaites Glacier is vulnerable to the marine ice-sheet instability as its grounding line is significantly influenced by ocean-driven basal melting rates, and its bedrock topography retrogrades into kilometer-deep troughs. In this study, we investigate how bedrock topography features influence the grounding line migration beneath Thwaites Glacier when extreme ocean-driven basal melt rates are applied. Specifically, we design experiments using the Ice-sheet and Sea-level System Model (ISSM) to quantify the SLR projection uncertainty due to reported errors in the current bedrock topography maps that are often used by ice-sheet models. We find that spread in model estimates of sea-level-rise contribution from Thwaites Glacier due to the reported bedrock topography error could be as large as 21.9 cm after 200 years of extreme ocean warming. Next, we perturb the bedrock topography beneath Thwaites Glacier using wavelet decomposition techniques to introduce realistic noise (within error). We explore the model space with multiple realizations of noise to quantify what spatial and vertical resolutions in bedrock topography are required to minimize the uncertainty in our 200-year experiment. We conclude that at least a 2 km spatial and 8 m vertical resolution would independently constrain possible SLR to ±2 cm over 200 years, fulfilling requirements outlined by the 2017 Decadal Survey for Earth Science. Lastly, we perform an ensemble of simulations to determine in which regions our model of Thwaites Glacier is most sensitive to perturbations in bedrock topography. Our results suggest that the retreat of the grounding line is most sensitive to bedrock topography in proximity to the grounding line's initial position. Additionally, we find that the location and amplitude of the bedrock perturbation is more significant than its sharpness and shape. Overall, these findings inform and benchmark observational requirements for future missions that will measure ice-sheet bedrock topography, not only in the case of Thwaites Glacier but for Antarctica on the continental scale.

1 Introduction

The future of the West Antarctic Ice Sheet (WAIS) is known to be one of the largest sources of uncertainties in global mean sea-level rise (SLR) on a century timescale (Schlegel et al., 2018; Yu et al., 2018). Its inherent instability results from much of its bedrock being below sea level (Fretwell et al., 2013) and its proximity to warm ocean temperatures (Schodlok et al., 2016). Continued acceleration of interior retreat of its grounding line and the loss of ice volume above floatation (VAF) could tip WAIS towards irreversible collapse (Rignot et al., 2014; Milillo et al., 2019). For these reasons, it is critical to understand the regional sensitivities in order to project future changes in the glacier. To characterize how the behavior of various ice shelves throughout WAIS may evolve in the future, the cryosphere community relies heavily on numerical modeling and simulations. Such tools allow physically based prediction and quantification of the sensitivity of grounding line evolution. This aids the community in determining which regional factors and characteristics most significantly contribute to the local stability of individual glaciers.

Within WAIS, Thwaites Glacier is a highly dynamic and sensitive region of the ice sheet (Robel et al., 2019; Rignot, 2001). Unstable and vulnerable to damage (Lhermitte et al., 2020), this glacier holds the potential for a ∼0.59 m global mean SLR and may initiate the collapse of WAIS (Holt et al., 2006). The glacier is held by pinning points present in undersea mountain ridges around its grounding line (Seroussi et al., 2017). Within 100 km upstream of the grounding line, Thwaites Glacier's bedrock slope becomes retrograde, splitting into trenches more than a kilometer in depth. As a result, its grounding line is considered marginally stable, presently between states of equilibrium and retreat (Payne et al., 2004; Seroussi et al., 2017).

Past model-based experiments suggest that two of the strongest controlling mechanisms for Thwaites Glacier are ocean-driven basal melting rates and bedrock topography (Larour et al., 2019; Nias et al., 2016; Schlegel et al., 2018; Waibel et al., 2018). Basal melting rates are difficult to accurately observe and model because stochastic evolution is prevalent in ocean circulation, especially under ice shelves. Temporal variability in ocean forcing further amplifies uncertainty in ocean-driven melt rates and predictions of future ice-sheet responses (Khazendar et al., 2019; Robel et al., 2019). Ocean circulation models physically estimate the dynamic evolution of the ocean-induced basal ice shelf melt rates. Though this can estimate melt rates at the higher spatial and temporal resolutions needed to inform ice-sheet models, accuracy can only be ensured up to the present-day grounding line. Once the current configuration evolves, extrapolation is required, meaning more uncertainty is contributed to the system (Marshall and Clarke, 1997; Nakayama et al., 2019; Zhou and Hattermann, 2020). If dynamic model coupling between the ice and the ocean is used, larger sensitivity issues can arise and amplify uncertainties as they are propagated through a model simulation (de Klerk and Voormeeren, 2008). As a result, a number of more simplistic models have been proposed, such as calculating constant basal melting rates based on incoming ocean heat convection (Rignot et al., 2016; Holland and Jenkins, 1999; Bondzio et al., 2018). Thermal forcing can also be used as functions of ice shelf depth due to ocean salinity and temperature gradients (Yu et al., 2018). Fundamentally, simplified parameterizations do not capture oceanic processes well. Indeed, recent model-based studies suggest that century-scale uncertainty in SLR potential under different basal melting rate or ocean thermal forcing scenarios has a significant spread (Yu et al., 2018; Schlegel et al., 2018; Seroussi et al., 2019), and projection sensitivity to these forcings vary substantially amongst ice-sheet models (Seroussi et al., 2020; Levermann et al., 2020; Edwards et al., 2021). Together, this all suggests that uncertainty in estimates of ocean forcing may remain a significant source (and perhaps the largest) of uncertainty in ice-sheet model simulations of Thwaites Glacier into the near future.

Figure 1 (a) Thwaites bedrock topography from BedMachine (Morlighem et al., 2020) is given with (b) its corresponding uncertainty (considered to represent 3σ) in the region. The data are interpolated onto our ISSM domain (later detailed) from its original grid map.

Bedrock topography, however, is more quantifiable as it is fixed in time on decadal timescales (i.e., aside from longer-term processes such as glacial isostatic adjustment and erosion and sedimentation) and is directly measurable using ice-penetrating radar (see Holt et al., 2006; Holschuh et al., 2020). Decades of remote sensing data from past missions have been used to construct estimates of present-day underlying bedrock topography. The measurements, however, are based on radar tracks, and thus physically informed interpolation procedures are needed to fill in the spatial gaps between measurements (e.g., mass conservation; Morlighem et al., 2020). Even these state-of-the-art products are therefore plagued by statistical uncertainties (especially in areas far from the radar tracks) and potential errors in radiogram interpretation. Consequently, even the most informed estimates of bedrock topography are associated with significant uncertainties in regions where measurements are sparse (Fig. 1).

In this study, we design a suite of ice-sheet model experiments to investigate how known uncertainties in bedrock topography affect 200-year simulations of the response of Thwaites Glacier to an extreme increase in ocean-driven basal melt rates. To accomplish this, we first perturb the domain bedrock within its known uncertainty to its minimum (bed minus error) and maximum (bed plus error) possible bedrock configurations. Using these new boundary conditions for simple sensitivity experiments, we characterize the spread in modeled SLR contribution resulting from error in bedrock topography (Sect. 3). Next, we further investigate model sensitivity to bedrock error using wavelet techniques to perturb the bedrock. In this case, the experiments are designed to derive the spatial and vertical resolutions needed to reduce uncertainty in modeled estimates of the glacier's future evolution within the requirements outlined by the 2017 Decadal Survey for Earth Science (Sect. 4). Finally, we perform an uncertainty quantification (UQ) sampling experiment to simulate probabilistic model outcomes followed by a sensitivity test to locate regions where the bedrock topography of Thwaites plays the most significant role in determining grounding line evolution (Sect. 6).

2 Model setup

2.1 Model description

We use the Ice-sheet and Sea-level System Model (ISSM; Larour et al., 2012b), a thermomechanical high-resolution ice-sheet model, to simulate the forward transient evolution of Thwaites Glacier. Our regional model is originally derived by extracting the Thwaites basin from an ISSM model of the Antarctic Ice Sheet (Schlegel et al., 2018; Seroussi et al., 2020). Here, the model uses initial surface and bedrock topography ($\stackrel{\mathrm{̃}}{x}$) from BedMachine v2 digital elevation model (DEM; Morlighem et al., 2020). Our sample space of possible bedrock realizations (±3σ) is derived from the absolute error provided with this product, which we consider to be equivalent to a 99.7 % confidence bound (or 3σ). The initialization and relaxation of our model follow the procedures described for the JPL_ISSM (Jet Propulsion Laboratory ISSM) model by Seroussi et al. (2020), including the determination of the basal friction coefficient over grounded ice and the ice viscosity of the floating ice using data assimilation techniques (Morlighem et al., 2010) to best match observed velocities (Rignot et al., 2011, 2017). For computational efficiency, all simulations are run with a two-dimensional shallow shelf approximation (SSA; MacAyeal, 1989), based on Schlegel et al. (2018; see Sect. 2.3 below), for stress balance approximations. For other detailed information about the model parameterizations and setup, including the treatment of basal friction and the rheology law, we refer the reader to Schlegel et al. (2018).

2.2 Mesh and boundary condition initialization

The upstream ice boundaries of the regional Thwaites Glacier domain are initially determined by the continental-model ice divides and are then modified following Yu et al. (2018). In regions where the boundary is too complex to distinguish at the resolution of the continental model, boundaries were modified to follow those of Rignot et al. (2019). At these boundaries, all thickness and velocity values are held constant as single point constraints during all simulations (Schlegel et al., 2013) and are specified by the thickness and velocity values of the continental model of Schlegel et al. (2018). A minimum ice thickness of 1 m is imposed to ensure a non-zero ice thickness. At the calving front, a free-flux condition is imposed, and the ice front positions are held constant, based on the mask from Morlighem et al. (2020). The grounding line evolves according to a sub-element evolution scheme and assumes hydrostatic equilibrium (Seroussi and Morlighem, 2018).

To create new high-resolution meshes adequate for our sensitivity studies to capture spatial variability in the BedMachine product, we use ISSM's static anisotropic mesh adaptation, informed by the gradient in initial surface velocities. Here, we define two different meshes, one for our vertical test (VT) experiments, with horizontal resolution defined between 200 and 1000 m, and the other for our spatial test (ST) experiments, with horizontal resolution defined between 50 and 200 m (Fig. S1). Due to the large amount of finely resolved elements in the ST model, we further reduce the domain in the upstream interior regions where ice velocities and thickness changes are found to not be significantly perturbed during our most extensive Thwaites grounding line retreat scenarios. This results in large reductions in computational costs. The VT model contains about 400 000 elements, while the ST model contains about 1.3 million elements. Because the mesh resolutions vary between spatial and vertical models, SLR results slightly differ between their control simulations. This result agrees with Seroussi and Morlighem (2018), who conclude that higher-resolution models result in more SLRs as compared to lower-resolution models of the same setup (especially resolution drops between 500 m and 2 km). Note that we observe this divergence primarily in cases where the grounding line migration rate is the largest, specifically when the glacier has already committed to full collapse at simulation year 150.

The friction and ice viscosity are taken from the continental model JPL_ISSM described by Seroussi et al. (2020), and the surface mass balance and ocean basal forcing are based on Schlegel et al. (2018). We interpolate these values onto our new meshes using bicubic interpolation, and based on Schlegel et al. (2018), they are held constant through all simulations.

2.3 Stress balance approximation

All forward simulations conducted here make use of the SSA stress balance approximation, chosen for its computational efficiency, which advantageously decreases the computational costs of running the large number of simulations required by our analysis. In addition, sensitivity experiments conducted on our VT model suggest that the use of a higher-order approximation (Blatter, 1995; Pattyn, 2003) does not significantly affect model results. These results agree with Schlegel et al. (2018), who state that a two-layer thin-film stress balance approximation (Schoof and Hindmarsh, 2010; Hindmarsh, 2004) has a statistically insignificant effect on grounding line sensitivity to perturbations in model boundary conditions with respect to SSA. Furthermore, we assume the effects of thermal variations are slow relative to the grounding line processes modeled here. After running our own sensitivity tests, we have confirmed that they do not significantly affect our results on the 200-year timescale investigated. Therefore, in the absence of a thermal model, a three-dimensional representation of the mesh is not required. As a result, for the purposes of this investigation, we take SSA as an acceptable approximation for stress balance. We set the stress balance mechanical equilibrium residual convergence criterion to 1 % and the stress balance velocity relative convergence criterion to 0.1 %, both tested under decreased convergence criteria to ensure repeatability, stability, and robustness.

2.4 Grounding line migration and forcing

Basal melting rates are calculated using the Massachusetts Institute of Technology General Circulation Model (MITgcm; Marshall and Clarke, 1997) using techniques described in Schodlok et al. (2016) and implemented based on Schlegel et al. (2018). That is, basal melt rates near the grounding line are extended inland using a nearest neighbor extrapolation method. This results in an aggressive basal warming from oceanic heat flux into newly formed ice shelf cavities as the grounding line retreats. Furthermore, we include a multiplier of ×1.8 on the melt rates in order to account for the maximum melting rate realistically achievable within the basin, in the case that Antarctic Bottom Water were to intrude underneath the ice shelf (Schlegel et al., 2018). This choice allows our experiments to explore the domain's full SLR contribution potential over the simulation period of 200 years.

The grounding line migration and grounding line friction interpolation are set using a sub-element on a partially floating elements scheme. The melt interpolation is set to have no melting on partially floating elements, which results in a conservative estimate of sea-level contribution compared to using a sub-element melting scheme. As a result of this setting, all fully ungrounded elements are subject to the melt rates that have been extrapolated into the ice-sheet interior. We choose this option because past studies have shown that the use of no melt on partially floating elements produces realistic results when spatial resolution needs to be compromised for computational efficiency. That is, using the no melt scheme, the modeled grounding line behavior more closely matches that of simulations using much more highly resolved spatial mesh in proximity to the grounding line (Seroussi and Morlighem, 2018).

We obtain an initial grounding line position by calculating where the BedMachine ice thickness mask has buoyant forces exceeding gravitational forces. Then, before applying the ×1.8 basal melt multiplier, we let the model run with its control bedrock (bedrock topography realization given by BedMachine) for 10 years to relax it and stabilize the grounding line position (Schlegel et al., 2018; Seroussi et al., 2011; Gillet-Chaulet et al., 2012). Comparing our grounding line position to the initial grounding line position in Yu et al. (2018), we observe that at worst we are overshooting the initial grounding line by approximately 10 km at the largest offsets (Fig. S2).

Before running a model simulation forward in time, we follow an algorithm on the glacier geometry within ISSM to ensure that the simulation does not experience a shock to the mass transport, stress balance, and grounding line solutions after bedrock perturbation is performed. These details are documented in the Appendix (Sect. A1).

Models are run on two Broadwell nodes with 28 cores each (56 processors total) on the Pleiades supercomputer cluster using ISSM version 4.16. Runs are set for 200 years using a time step of about ∼6.1 d or 60 time steps per year.

3 Experiment 1: minimum and maximum bedrock resulting spreads

To begin our investigation, we run a simple experiment on our 200-year forward simulation to test model SLR contribution sensitivity within our defined bedrock sampling space. That is, we run two simulations, one using the minimum and one using the maximum possible bedrock topography ($\stackrel{\mathrm{̃}}{x}-\mathrm{3}\mathit{\sigma }$ and $\stackrel{\mathrm{̃}}{x}+\mathrm{3}\mathit{\sigma }$), where 3σ represents the map of bedrock error reported by BedMachine. Together, the results of this experiment bound the possible SLR contributions (calculated from the change in VAF from grounding line migration) and final grounding line positions that our forward simulation could yield when forced with variations in bedrock topography within error. The results of these experiments, presented in Fig. 2, are a simple first-order illustration of the magnitude of decadal-scale sea-level projection uncertainty sourced in present day Thwaites bedrock error.

Figure 2 Results of forward model experiments using the control (orange), maximum (yellow), and minimum (red) bedrock configurations are shown. Two timestamps are depicted: 150 years (dashed) and 200 years (solid). At 200 years, the control, maximum, and minimum bedrock configuration models yield 21.1, 4.8, and 26.7 cm of global mean SLR respectively. The resulting SLR difference between the models of minimum and maximum grounding line retreat is 21.9 cm. See the Appendix (Sect. A3) for the 150- and 200-year tabulated SLR values.

The results presented in Fig. 2 suggest that under the maximum bedrock scenario, at least temporary stabilization can occur within the first 50 km of the current grounding line position during the 200-year period investigated in Experiment 1 (Fig. S3, ridge set a). In contrast, the ridge farther upstream (Fig. S3, noted as ridge set b) is not capable of stabilizing Thwaites under our extreme warming conditions; instead it merely slows the widespread retreat. These results support conclusions posed by Morlighem et al. (2020) that continuous retreat in Thwaites occurs (given our modeling assumptions) if the glacier retreats past its initial ridges (Fig. 2a in Morlighem et al., 2020; Fig. S3, noted as ridge set a), especially under an extreme ocean warming condition. Since the majority of the bedrock upstream of these initial ridges is more than a kilometer deep, the presence of significant pinning point features is indeed improbable, even within the sampling space of large errors that exist in the region. Overall, we conclude that the current error in bedrock topography is responsible for a global mean SLR difference of 21.9 cm (maximum SLR − minimum SLR) within 200 years under a forcing of extreme ocean warming.

We also find that the control bed topography simulation itself (without perturbation) forced with extreme basal melt rates results in significant retreat in a 200-year period (Fig. 2); the bedrock elevation is too low for any pinning points to create resistance. The maximum bedrock scenario (as compared to the other models), on the other hand, illustrates that our model results are highly sensitive to the bedrock error. We see here that the grounding line migrates minimally, only about 50 km upstream over our 200-year forcing period. Indeed, we find that the grounding line evolution and consequential sea-level contribution can be significantly altered (between 21.1 and 4.8 cm) simply through an increase in the control bed topography within error. Therefore, within possible realizations of bed topography considered here, unobserved features or realistic noise not captured in our control bedrock may constitute distinct pinning points that could stabilize or delay retreat. This means that it is possible for not-yet-observed bedrock features that exist within the current bedrock error bounds to play critical roles in dictating the future retreat rate of Thwaites Glacier.

It is also important to note that results suggest that there is little difference in retreat and SLR contribution between the minimum and control bedrock topography scenarios. This is due to model sensitivity to the extreme (×1.8) basal melt rates, which leads to the ungrounding of most of the ice in the model domain for both simulations; in other words, the minimum scenario has retreated to higher ridges deep within the domain such that there are no more opportunities for further retreat and SLR. The grounding lines for both runs migrate through Thwaites Glacier's deep trenches over the simulation period, resulting in significant loss of ice VAF and a SLR contribution close to the domain's maximum potential.

The results show that the two extreme bedrock scenarios (maximum vs. minimum) diverge significantly in a 200-year period, suggesting that Thwaites Glacier is highly sensitive to bedrock perturbations given the simulations' aggressive basal melt forcing. Therefore, we find that, within error, perturbations in bedrock topography are capable of slowing grounding line retreat in response to extreme ice shelf melt rates, and we cannot assume that the current vertical error in bedrock topography is negligible within the Thwaites basin. Consequently, our results suggest that in order to obtain high confidence in robust ice-sheet model projections under an aggressive ocean warming scenario, the existing bedrock error bounds within Thwaites Glacier must be reduced.

4 Investigating resolution requirements for uncertainty minimization

Our initial results show that a large projected SLR uncertainty exists with consideration to the currently reported bedrock error. In other words, our forward simulations of Thwaites Glacier evolution are sensitive to perturbations in bedrock topography. We now ask, if the current error ranges are not sufficient, what bedrock uncertainty range would be required to accurately model Thwaites Glacier evolution and at what spatial resolution.

Since an elevation increase or decrease in the entire bedrock topography is not a realistic representation of the effects of error on the reported BedMachine topographical representation, we use discrete wavelet decomposition (Sect. 4.1) in order to create noise amplifications in our bedrock topography array. Two-dimensional discrete wavelet decomposition involves using discrete signals of differing shapes and sizes in order to decompose an array into four subarrays: the A_n array (low-frequency approximation), H_n array (horizontal high-frequency), V_n array (vertical high-frequency), and D_n array (diagonal high-frequency).

We amplify the high-frequency filters (H_n, V_n, and D_n) so that, upon recomposition of the subarrays, we introduce realistic physical noise into our experiments. Because high-frequency filters have large coefficients in locations with high frequencies, many of the pinning points (ridges, mountains, etc.; Fig. S3, noted as ridge set a) that are capable of slowing grounding line retreat in Thwaites Glacier are, by design, targeted by our wavelet image processing. By focusing on these present pinning points, our analysis is more likely to result in a wide range of resultant model SLR contribution and grounding line retreat. An example of how a pinning point may be altered by wavelet amplification techniques may be observed in Fig. S4.

Here, we design two sets of experiments to derive spatial resolution and vertical error requirements for ice-sheet model bedrock topography. The goal of these experiments is to determine what these independent requirements would need to be in order to lower the SLR uncertainty of our forward simulations to ±2 cm. We choose ±2 cm of SLR based on the 2017–2027 Decadal Survey for Earth Science and Applications from Space (National Academies of Sciences, 2018) which states that it is one of the most important science applications to “quantify the rates of sea-level change and its driving processes at global, regional, and local scales, with uncertainty <0.1 mm/yr” (Chapter 3, S-3a). If we assume a maximum ±0.1 mm yr⁻¹ of SLR uncertainty, such a requirement accumulated over 200 years is equivalent to ±2 cm of SLR uncertainty.

We acknowledge that this is a strict constraint for an ice-sheet model projection system and that our results will represent stringent ideal high-end requirements. Nonetheless, we consider it a meaningful benchmark for the quantification of uncertainty in regional projections of glacial contribution to sea level. Following such, spatial and vertical resolution requirements are quantified throughout the proceeding sections (Sect. 4.2 and 4.3) as they relate to restricting SLR uncertainty within the ±2 cm range.

4.1 Wavelet decomposition setup

We use discrete wavelet transform (DWT) for a two-dimensional decomposition and recomposition of bedrock topography. Wavelets form basis functions for projection that reveal the waveform (or frequency content) of the signal around a given location. The method relies on first multiplying the signal by an oscillating function (the wavelet) that vanishes away from the location at which it is applied and then integrating that product over space. The result is a coefficient that measures the similarity between the signal and wavelet around that point. The wavelet can be spatially rescaled to analyze the signal at different wavelengths; hence it is widely used as a tool to study the space–frequency content of a signal. Here, we also use the filtering capability of wavelets to generate new realizations of the bedrock topography used as input for our ice model runs.

DWT specifically operates using a pair of wavelets labeled the low-pass and high-pass filters which are combined in two dimensions and used to convolve the signal into four sets of coefficients: the approximation coefficients (A_n; dual low-pass), which represent a low-frequency (smoothed) version of the signal, and the three detail coefficients according to their spatial orientation, i.e., x-axis coefficients (H_n; high-pass x axis, low-pass y axis), y-axis coefficients (V_n; low-pass x axis, high-pass y axis), and diagonal coefficients (D_n; dual high-pass). The wavelet decomposition can in turn be reapplied using rescaled wavelets on the resulting A_n matrix to provide a new set of detail coefficients H_n+1, V_n+1, and D_n+1 of lower frequency than H_n, V_n, and D_n, as well as a new approximation matrix A_n+1. This iterative process is used to characterize the signal at different decomposition and resolution levels here indexed with n (Daubechies, 1992; Mallat, 1989; Meyer, 1995).

We use the second Daubechies wavelet (hereafter referred to as db2; Daubechies, 1992) and the discretized Meyer wavelet (hereafter referred to as dmey; Meyer, 1990; Abry, 1997) for our noise study. We decide to use these wavelet bases due to their opposing strengths and weaknesses (Abry, 1997; Daubechies, 1992). The db2 wavelet can localize signals well in the space domain but poorly in the frequency domain, while the dmey wavelet localizes signals poorly in the space domain but well in the frequency domain. Therefore, they jointly analyze the full effect of perturbation on our DEM, helping us to observe the SLR sensitivity of Thwaites in both space and frequency domains.

Using these wavelets to obtain a post-decomposition state, we introduce a multiplier to the H_n, V_n, and D_n coefficient matrices in order to amplify existing bedrock topography noise. Therefore, upon recomposition of the original image (A_n−1) using the A_n approximated coefficients, we obtain a similar bedrock topography that we previously had with the addition of introduced noise on areas with existing high-frequency features.

In some cases, we further decompose the A_n coefficient matrix using the same wavelet in order to achieve a multiresolution analysis. Upon higher-level decompositions (n≥2), the high-frequency features targeted change to be of larger spatial scale due to single pixels now representing larger areas. Therefore, noise amplification of solely the n layer at these large orders produces geological landforms such as entire ridges and trenches (Daubechies, 1992; Mallat, 1989; Meyer, 1995). We take advantage of this mathematical notion to develop a function to decompose our DEM to an n level, amplify the H_n, V_n, and D_n coefficient matrices, and then recompose accordingly.

We choose not to conserve ice volume when creating the bedrock realizations in order to simplify the wavelet amplification process. We also maintain the ice surface DEM in places above floatation (see Sect. A1 for situations in which the ice surface DEM may be altered). As a result, more divergent results are expected between different bedrock realizations as some models may have more or less bedrock added at critical pinning points. Sensitivity tests suggest that this divergence, however, is minimal and does not affect our results.

4.2 Experiment 2: spatial resolution degradation results

To understand the spatial resolution required for modeling Thwaites within our given regional domain, we use a mesh that captures all features present in the BedMachine DEM; that is, the ST mesh has a higher spatial resolution (200 m; Sect. 2.2) than BedMachine (500 m). To determine the bedrock topography and the bedrock error for the ST mesh, we use bicubic interpolation. Using such a finely resolved mesh affords us a bedrock topography that we can manipulate at a higher resolution in order to test the impact of noise below the current 500 m spatial resolution.

There is observational evidence that high-resolution features, not fully captured by the current regional-scale bedrock topography maps used by ice-sheet models, may exist in Thwaites Glacier (e.g., Boon, 2011; Schroeder et al., 2014; Chu et al., 2021) and that high-resolution basal features within glaciers neighboring Thwaites Glacier, such as in Pine Island Glacier, also exist (e.g., Rippin et al., 2011; Bingham et al., 2017). In Pine Island Glacier, it has been observed that the basal drag inversion would improve if the bedrock topography is considered noisy, suggesting that more accurate modeling of glacier stress balance may depend on whether or not a simulation can capture high-frequency basal terrain (Kyrke-Smith et al., 2018).

In the case that uncaptured high-frequency bumps are present in the bedrock topography, the DEM needs to be interpolated to a higher resolution before wavelet techniques are performed. We begin by upgrading the bedrock realization to a 200 m resolution, amplifying the high-frequency filter noise (through techniques and variations described in Sect. 4.1), and then degrading the resolution to 400 m. We use this 400 m resolution as the control resolution for each specific wavelet–resolution combination. A 400 m spatial resolution is chosen as our control resolution since, at this resolution, a single feature (around 400 m) will be described by multiple mesh elements. At higher resolutions than this, we would only be representing features close to a 200 m spatial resolution with a single element which shocks the system upon degradation. Next, image degradation of the DEM resolution to various spatial resolutions ranging from 400 m to 30 km follows through bicubic interpolation. Results of SLR contribution difference from the control run at year 200 for the 258 forward simulations for this experiment are shown in Fig. 3.

Figure 3 Final SLR is plotted against the current downgraded resolution for various noisy bedrock topography (a, results for all spatial scales; b, results for spatial scales of 6 km and finer). The model results (dots) are spline interpolated with trends depicted with a mean resultant line (dark blue line) and a percentile band (shaded region) of ±1σ (68.3 %). See Appendix (Sect. A2) for further details on bedrock perturbation codes from the legend.

In Fig. 3b, we find that results converge under 2 km for all of the constructed maps. Within the first 2 km, many of the positive high-frequency noise perturbations about the realization (previous wavelet noise amplifications performed at a 200 m level) are either amplified or reduced from elements due to the degradation of resolution (resultant from spatial interpolation). Therefore, relatively minor changes in SLR are observed as compared to the respective control runs due to the compensation of bedrock about these elements. However, after passing the ∼2 km threshold, much of the noise is completely lost by the wide grids used in the interpolation approximation, and after this 2 km threshold, we observe a large sudden gain in SLR difference. Since we observe a clear compensation before 2 km, while beyond 2 km they are no longer compensated for, we conclude that the spatial degradation threshold for proper accuracy lies around this point.

In Fig. 3a, we find that results lack an overall general trend, with random changes in SLR resulting as the resolution is degraded and no obvious pattern is observed. This results from changes in bedrock feature characteristics (e.g., ridge omission, mountain expansion) due to our degradation of spatial resolution (accomplished through spatial interpolation) which leads to highly unrealistic topography as it crosses the bedrock error range in BedMachine (only when the spatial degradation, resulting from spatial interpolation, occurs at extremely low resolutions).

With consideration to the desired ±2 cm of SLR requirement, we evaluate our minimum spatial resolution limit to be approximately 2 km, particularly at Thwaites' pinning points, as SLR does not approach the 2 cm threshold again for any of the lower spatial resolutions investigated here (Fig. 3a).

4.3 Experiment 3: vertical resolution noise construction results

To derive vertical resolution constraints in Thwaites, we design a set of experiments to systematically perturb our bedrock topography with noise of increasing amplitudes. We take BedMachine's bedrock topography to be our control model and then use wavelet decomposition to amplify the high-frequency noise (similar to Sect. 4.2). Using either the db2 or dmey wavelet in combination with four decomposition levels, we create eight bedrock topography realizations for the model. We then take the absolute value of the noise created for each realization, resulting in perturbations that are all positive. This allows us to apply a noise realization as either a positive (through addition) or a negative (through subtraction) set of perturbations and then to examine the impact of both types of noise on the model simulations independently as opposed to having to deconvolve the impact of the two types of perturbations in combination.

To test the effect of vertical changes in bedrock topography, we establish a set of limitations on the maximum height to which new perturbation features can be built. For each mesh vertex in our model bedrock topography map, the bedrock altitude is equal to the minimum of (a) the perturbed bedrock realization, (b) the maximum feature height change (Fig. 4, x axis), and (c) the bedrock error limit ($\stackrel{\mathrm{̃}}{x}-\mathrm{3}\mathit{\sigma }$ and $\stackrel{\mathrm{̃}}{x}+\mathrm{3}\mathit{\sigma }$ for negative and positive noise respectively) on that vertex. For our experiments, we systematically change the amplitude of our bed topography noise by varying (b) the maximum feature height change between 10 and 100 m with 10 m increments, giving a total of 20 final bedrock topography realizations per wavelet-decomposition-level combination after both positive and negative noise is included, for a total of 160 model simulations.

Figure 4 Final SLR plotted against the maximum bedrock change threshold designated for each model simulation. Points represent the model results connected using a spline interpolation. Eight different bedrock topographic realizations are made with high-frequency amplification for this experiment. See Appendix (Sect. A2) for further details on bedrock perturbation codes from the legend.

In Fig. 4, we present results from all forward simulation experiments for the vertical resolution tests. We find that the results are nearly linear within ±20 m of bedrock addition with asymptotic behavior dominating with increasing max height perturbation for both the positive and negative addition of noise. The results are biased towards the bedrock error due to different bedrock realizations from various perturbations being leveled off by (c) the bedrock error limit of each pixel.

For negative bedrock changes (bedrock subtraction), the asymptotic behavior observed from about −100 to −20 m is due to the simulation reaching the limit of maximum glacial retreat. At this point in the modeled evolution of the grounding line, a complete Thwaites destabilization and retreat is achieved, meaning that there is little ice left in the domain that can unground into deep bedrock troughs. Therefore, SLR potential beyond this point is limited since our SLR calculation is calculated from changes in VAF. Between −20 and 0 m of bedrock subtraction, however, the SLR contribution behaves more linearly as full retreat has not yet occurred by the time these model simulations have been completed.

For positive bedrock changes (bedrock addition), we still observe the destabilization of Thwaites between about 0 and 20 m of maximum bedrock addition. Less retreat has occurred by the time 200 years have passed since the topographically higher bedrock realizations create a delay in the retreat rate of the grounding line (Fig. 2). There are not sufficient pinning points present in these simulations to completely prevent retreat. Between 20 and 100 m, however, the glacier obtains much more effective stabilization due to the topographically higher bedrock levels. We find these results to be the only ones in which the different wavelet-decomposition-level combinations diverge from one another. As a result, there is almost 5 cm of SLR difference between the least and most stable models at 100 m of maximum bedrock addition.

Our results from these experiments strongly suggest that amplitude of bedrock change matters more than the shape of the perturbation itself. There is little variance between the differing runs and significant agreement in terms of general trend despite our application of different wavelet shapes within the ensemble. As the two types of wavelet shapes are used at multiple resolutions, no shape seems to have a significantly larger effect on our overall model results. Moreover, because different volumes of bedrock are added to each run, one would also expect larger divergences from the general trend due to alteration of the available VAF for contribution to SLR. However, we find that this effect is negligible and that the small spread between the various simulations at the same bedrock alteration limit only lends further support that it is the amplitude of change rather than the type of wavelet, which more heavily determines simulation SLR contribution.

With respect to our ±2 cm of SLR requirement, we focus on the control run at 0 m of bedrock addition, where within Fig. 4, the tangential slope is the largest. From this point, we observe that vertical bedrock perturbation of a magnitude within approximately ±8 m at areas where the model is most sensitive to bedrock would hold the simulation SLR uncertainty within our target range. It is important to note that, in comparison to variations in spatial resolution, we find that vertical perturbation in bedrock topography contributes more significantly to uncertainty in ice-sheet model estimates of SLR contribution from Thwaites Glacier, especially since even a meter of perturbed vertical bedrock elevation has the potential to alter model results significantly.

5 ISSM-DAKOTA framework

For further investigation of the sensitivity of model simulations to topographic data accuracy in Thwaites, we take advantage of the Design Analysis Kit for Optimization and Terascale Applications (DAKOTA) software from Sandia National Laboratories (Eldred et al., 2008) that has been embedded within the ISSM framework (Larour et al., 2012a; Larour and Schlegel, 2016). DAKOTA, a tool for UQ analysis and statistical error quantification, has traditionally been used to perturb ISSM model input and boundary conditions within various regional domains in order to perform a variety of sensitivity tests on model output diagnostics, like mass flux (Schlegel et al., 2013, 2015; Schlegel and Larour, 2019) and regional mass balance (Schlegel et al., 2018).

In this study, we first use the DAKOTA software to create a statistical sampling of bed error perturbations in order to isolate bedrock pinning points. Here, we launch an ensemble of forward models, applying the predetermined set of perturbations to the bed using the initialization procedure detailed in the Appendix (Sect. A1). We treat the entire domain as one single partition such that the bedrock everywhere is perturbed by the same percent error for any given ensemble member (Schlegel et al., 2018). The perturbation set consists of a normal sampling distribution of 300 samples within a ±3 % bedrock standard deviation, the mean BedMachine standard error determined for the Thwaites basin based on the error supplied by Morlighem (2020). We use the Latin hypercube sampling (LHS) algorithm (Swiler and Wyss, 2004) to generate the sampling distribution, based on Schlegel et al. (2015).

This first sampling experiment informs us at which locations, under various bedrock configurations within the reported BedMachine standard error (i.e., 3σ), the grounding line prefers to reside. Next, we use these results to isolate the most influential region of bed topography within our domain. We then design an additional spatial sensitivity sampling experiment to determine which specific features within that region have local errors that would influence Thwaites grounding line stability the most. To do so, we partition the targeted region into 400 equal-area subsections using the Chaco software (Hendrickson and Leland, 1995), based on Schlegel et al. (2013). We then perturb the bedrock in each partition one at a time by adding the maximum local BedMachine standard error for that partition (Morlighem et al., 2019) and finally run a forward transient, resulting in a total of 400 simulations. Comparison of the grounding line behavior between these runs allows us to identify pinning points by quantifying which individual bedrock regions (partitions) have the largest effect on modeled sea-level contribution when perturbed within error.

5.1 Experiment 4: UQ probabilistic distribution results

As previously discussed in Sect. 5, we take advantage of the DAKOTA software in order to understand the probability of divergence from our control model. More specifically, we use DAKOTA functionality to derive a normally distributed perturbation set, with a 3σ sampling range of about ±9 %, equivalent to the mean uncertainty of BedMachine error within the Thwaites basin. The perturbation set consists of 300 values, by which we multiply our control bedrock topography over the entire model domain. Each of these new bedrock realizations are run with the basal melt rate multiplier of ×1.8 and a multiplier of ×1.4, resulting in a total of 600 forward simulations.

In order to compare our results to other SLR drivers, we also use DAKOTA to derive a uniform sampling of ice shelf basal melting rate multipliers ranging between ×1 and ×1.8 in order to represent a full range of possible future ocean-warming scenarios, based on Schlegel et al. (2018). We then repeat a third series of sampling experiments with 300 of these melting rate multiplier perturbations and a consistent control bedrock. In this way, we are able to compare the impact of errors in bedrock topography against that of uncertainty in ice shelf basal melting rates, which have been found to contribute a significant amount of uncertainty to decadal-scale simulations of the Antarctic Ice Sheet (Schlegel et al., 2018).

Figure 5 Distributions of event probability as a function of SLR contribution after 200 years of simulation for three different boundary condition sampling experiments, each experiment consisting of 300 model simulations. Experiments include sampling of bedrock topography under two different ocean warming scenarios: one model scenario with a basal melt multiplier of ×1.4 (pink) and one model scenario with a basal melt multiplier of ×1.8 (blue). For comparison, we also present results from sampling of basal melt rates between ×1 and ×1.8, with no changes in bedrock topography (yellow).

We find that under the control bedrock and variable basal melt rates, results exhibit a bimodal distribution absolutely ranged between 4 and 21 cm of SLR contribution. This suggests that if no error is present in bedrock topography, the second set of ridges inland (Fig. S3, noted as ridge set b) can still play an important role in stopping, or delaying, Thwaites' grounding line retreat on the 200-year timescale investigated here. The uncertainty range due to basal melting rates is comparable to the simulation uncertainty in SLR resulting from bedrock variation under our most extreme ocean warming scenario (Fig. 5, ×1.8 Melt). We find that under this extreme forcing and without alteration to bedrock topography full grounding line retreat is almost achieved, and under lower-end basal melt forcing, i.e., present-day basal melt rates, stability is maintained. The resulting distribution curve, however, is not normal, and its bimodal nature suggests that melt rates must exceed a specific threshold to ensure a SLR contribution of about 15 cm or greater. The bimodal response also implies that the melt multiplier threshold is dependent on the location of bedrock pinning points beyond Fig. S3, ridge set (a), that must be surpassed in order to achieve full retreat (i.e., Fig. S5).

Under the maximum basal melting rate (×1.8) and a variable bedrock configuration, Thwaites achieves a 21 cm distribution between an absolute range of 6 to 27 cm of SLR. This large spread in SLR contribution after 200 years of simulation suggests that our model simulations are highly sensitive to the currently reported bedrock error. Furthermore, this distribution's SLR peak probability has a larger magnitude than that of our control SLR contribution. This implies that the grounding line retreat under our most extreme ocean forcing scenario is inherently unstable. Previous results (Fig. 4) suggest that the negative skew of the curve and tail with lower SLR contributions is likely a response to positive bedrock perturbations. Specifically, the addition of bedrock at influential pinning points promotes glacier stability. Since the sensitivity of modeled grounding line retreat rates is also highly dependent on ocean forcing (i.e., basal melt rates), as discussed above, it is clear that model results are strongly driven by a combination of bedrock topography configuration and ocean-driven floating ice melt rates. This is illustrated in Fig. 5 by bedrock sampling performed with the lower melt multiplier of ×1.4. The use of this multiplier results in a range of SLR contributions that span the area between the two peaks of the melt sampling's distribution and results in a smaller simulation uncertainty than that exhibited by the other two curves. This ranges between an absolute 4 and 14.5 cm of SLR contribution.

5.2 Experiment 5: pinning point sensitivity test

To isolate the pinning points and understand their weight on retreat, we design a spatial sensitivity test to determine which areas of bedrock topography are most significant in determining the resultant SLR. Due to the large number of elements within our domain, it is not suitable to perform the sensitivity test over the entire domain as there is a large computational cost. Therefore, we use the distribution results from the bedrock sampling with ×1.8 melt (Fig. 5) to determine where the most influential grounding line pinning points may be, or in other words, at what model mesh vertices the grounding line prefers to rest during the 300 bedrock sampling simulations.

Taking grounding line positions at 10-year intervals for all 300 simulations, we calculate the percent time that the grounding line falls within 100 m of each mesh vertex (Fig. S5). We take a threshold of 20 % or above to denote locations of strong topographical influence on the grounding line such as pinning points and depressions. We appropriate this percentage by examining the sensitivity of expected SLR contribution to perturbations in bedrock at locations above various values of percent thresholds. Specifically, we conduct an experiment, in which we add the maximum bedrock error (3σ) to each mesh vertex above a minimum percent threshold and then run a set of simulations varying the threshold value. Results suggest that a threshold of 20 % would capture the majority of the influential bedrock pinning points (Fig. S6).

Because all resulting points of interest are in proximity to the initial grounding line, we isolate a single subdomain in which to perturb the domain's bedrock topography for our sensitivity test. We decide on 25 km² for the spatial resolution, yielding 400 partitions total and requiring a total of 400 forward simulations. Ideally, we would conduct this study using a spatial resolution of about 4 km² (∼2 km spatial resolution; Sect. 4.2). A spatial resolution of 25 km² is chosen instead, however, because it drastically decreases the number of partitions, the required number of simulations (one per partition), and therefore the computational price of the sensitivity test, lowering the computation cost by over 75 %. Because we expect our sensitivity experiment to consider three independent variables – bedrock error, mean feature altitude, and spatial location – we decide to perturb the bed by the maximum bedrock error. Choosing a constant perturbation change or a ratio instead would not have fully taken all three of these variables into account.

Figure 6 (a) A plot of the SLR change created by a partition's independent perturbation. (b) A plot of the SLR change created by a partition's independent perturbation is divided by the magnitude of the bedrock topography change for each mesh element. The quantity dSLR/dZ is unitless.

Figure 6 suggests that model results are most sensitive to changes in bedrock topography closest to the present-day grounding line and local perturbations to bedrock in proximity to the grounding line result in the greatest divergence in modeled SLR contribution. We also note that some locations too are more sensitive to ice thickness losses rather than bedrock volume additions, particularly in the areas where SLR rises in response to a bedrock increase.

We conclude based on these results that the placement of the pinning points matters more than the shape of the pinning points. As observed in Sect. 4.2 (“Spatial resolution degradation”), there was little SLR divergence between different wavelets and spatial resolutions under 2 km spatial representation of the bedrock topography. In Sect. 4.3 (“Vertical resolution noise construction”), a similar trend was seen as different wavelets and different spatial resolutions resulted in little divergence. Despite using two opposing wavelets, dmey (which focuses on frequency retention) and db2 (which focuses on signal localization), there is little evidence to suggest that either had any role in creating divergent results. The only SLR differences between different wavelet perturbations (arguably negligible) can be attributed to the different amplification multipliers and spatial resolutions.

We also conclude that our results suggest the spatial size of perturbation is negligible compared to the location of perturbation. As stated above, neither of the results of the experiments discussed in Sect. 4.2 (“Spatial resolution degradation”) and Sect. 4.3 (“Vertical resolution noise construction”) suggests that the spatial resolution of perturbations has any major effect on results. Figure 6 similarly does not support any clear conclusion on the effect of spatial perturbation size. Therefore, with limited evidence available to defend such conclusions, we believe further tests will be necessary to examine the model's sensitivity and fully justify this claim.

These sensitivity results suggest the primary bedrock controls of Thwaites exist near the present-day grounding line features. The sharpness and shape of perturbation do not strongly impact model results, but rather the amplitude and location of the chosen perturbation are more influential bedrock characteristics.

6 Discussion

Our experiments aim to quantify the uncertainty in an ice-sheet model forward simulation of Thwaites Glacier over a 200-year period of extreme ocean warming. Specifically, we focus on characterizing the bedrock's potential to stop or delay grounding line retreat and glacial collapse. From these experiments, we also aim to determine how uncertainty in ice-sheet model estimates of Thwaites' SLR contribution might be reduced through more accurate and precise data measurements.

Our results suggest that bedrock is an important source of uncertainty in Thwaites Glacier model projections, and in order to constrain projection uncertainty, it is important to minimize present day bedrock data uncertainty when modeling the behavior of an unstable glacier in response to extreme ocean-driven ice shelf basal melt rates. This is especially the case for our aggressive ocean warming scenario for Thwaites Glacier, which suggests that bedrock topography error could produce an overall spread of 21.9 cm in sea-level contribution (Fig. 2). This methodology assumes that the minimum and maximum bedrock topographies represent the maximum and minimum grounding line retreats respectfully. However, there are cases where additions of bedrock may encourage grounding line retreat and where bedrock removal may discourage grounding line retreat (Fig. 6; Nias et al., 2016; Koellner et al., 2019). We have analyzed the curves from our DAKOTA runs (Sect. 5) and observed that the minimum and maximum grounding line retreats sit at the extremes of these curves. This suggests we capture most of the responses of the model with our assumption.

In order to identify the requirements for spatial and vertical bedrock resolution that would reduce our modeled Thwaites SLR uncertainty to ±2 cm, in accordance with goals outlined by the most recent Decadal Survey for Earth Sciences, we design a unique set of experiments that take advantage of discrete wavelet decomposition techniques. Experiments that test spatial resolution requirements suggest that we need bedrock geometry to be measured at a spatial resolution of 2 km or less in order to accurately characterize the bedrock topography pinning points and to minimize ice-sheet model uncertainty on 200-year timescales. In particular, we find that at a resolution finer than approximately 2 km (at sensitive glacial regions), estimated SLR contribution converges (Fig. 3). For our assessment of vertical accuracy requirements, we conclude that vertical error needs to be known to an accuracy of ±8 m in order to constrain simulation uncertainty to ±2 cm, especially at major glacial pinning points (Fig. 4) where the grounding line is most likely to be affected by perturbations in bedrock topography. In fact, we find that our model simulation is highly sensitive to all perturbations in bedrock topography regardless of the spatial resolutions and shape of the bedrock noise applied.

We also find that, with respect to estimated SLR contribution, perturbation sharpness and shape of bedrock features are negligible in comparison to their amplitude and geographic location. This is because variations in perturbation shape and resolution for both the spatial and vertical resolution experiments result in little divergence in SLR contribution (Figs. 3 and 4). Results also suggest that vertical resolution results are driven by the amplitude of noise, while the spatial resolution results are affected by the mitigation of high-frequency vertical noise (in sensitive glacial regions). These results offer further justification that the amplitude of change in perturbations plays a larger role than perturbation shape. In addition to the amplitude of a bedrock feature, the location of perturbation also strongly influences our model estimates of the SLR contribution from Thwaites Glacier. This is because results of our spatial sensitivity experiment suggest that changing the ice thickness locally within a partition does impact our simulation beyond the consequences of altering just the bedrock topography. We also find that almost all of the most influential partitions are within ∼50 km of the present-day grounding line, suggesting that pinning points are abundant in proximity to Thwaites' rough undersea mountain ridges that persist just upstream of where the glacier is currently grounded.

It is important to restate that our model simulations are forced with highly aggressive basal melting rates. In a perfect world, we would want to reach bed data at the extreme precision we detail in order to cover all possible situations. However, as the basal melting rate used in our model is extreme, a lower melting rate would yield less spread in SLR contribution and would therefore allocate less precision in spatial and vertical bedrock to reach our ±2 cm SLR goal. In that case, DEMs with worse resolution could be acceptable depending on the modeling data sets and configurations used.

Results suggest that our findings are not dependent on melt rates but rather the grounding line sensitivity to feature size. Figure S7 represents an additional trial of Experiment 2 (Sect. 4.2) with a ×3.6 basal melt multiplier (a physically improbable melt rate multiplier as detailed in Schlegel et al., 2018). The resolution threshold observed decreases to 4 km due to larger melt rates impeding the ability for pinning points before the second set of ridges (Fig. S3b) to prevent grounding line retreat. Therefore, all scenarios in the spatial resolution trials with a ×3.6 basal melt multiplier witness Thwaites Glacier's complete collapse. This narrows the possible uncertainty range, in which collapse is highly probable, and decreases the spatial resolution required. We believe bedrock roughness does not have a major effect on grounding line retreat if the melt rates encourage grounding line retreat far beyond the indicated ridges prior to bedrock perturbation. Our previous spatial resolution suggestion of 2 km (Sect. 4.2) is reasonable since the melt rate will likely not exceed ×1.8 within the next century, and care should be taken in improving the certainty of whether or not the ridges of highest influence (Fig. S3b) are capable of slowing down, or stabilizing, grounding line retreat within the coming centuries.

It is also important to note that wavelet decomposition focuses on the amplification of noise that is already present in the bedrock topography map, and therefore it does not create new perturbations. As a result, if pinning points are present in low-frequency zones, it may not be sensed by the wavelet tests. We, however, assume this does not strongly affect our conclusions as we find that the most influential bedrock regions are almost exclusively near the grounding line where high-frequency noise is prevalent. We also recognize some regions may witness a “clipping off” of the top of perturbation shapes in Sect. 4.3 due to changes in amplitudes going above the maximum bedrock error. Therefore, if SLR has more dependence on the perturbation's curvature or shape than previously suggested, we may lose information regarding the overall shape effect. For spatial resolution degradation experiments (Sect. 4.2), we observe gridding effects when building noise on top of the 200 m interpolated DEM. Based on additional sensitivity experiments, we do not believe this creates error in our experiments, but it is important to note this could have some effect on the resulting bedrock realizations.

Furthermore, our methods may overrepresent the SLR response as Thwaites Glacier is affected by vertical land motion (VLM; Larour et al., 2019). GPS measurements in regions of bedrock uplift have shown trends of VLM in the range of tens of millimeters per year which has the potential to affect results (Barletta et al., 2018). However, as the aggression from our high basal melting rate rapidly recedes from the grounding line, we do not believe the uplift response has a major effect on our model as pinning points likely would not translate quickly enough to stabilize the glacier. Results likely would be mitigated by this contribution in milder circumstances. As modern uplift couplings become online, we hope to further explore their contribution.

Model SLR contribution results can have variable results when other bedrock topography DEMs of Thwaites are used. Different DEMs have significant SLR changes as we find that the model is highly sensitive to even slight changes to the mean bedrock elevation. In this case, our results may be biased towards BedMachine. To test this, we did try to compare our model results against a simulation that uses a DEM of Thwaites made from geostatistical data from paleo-observations (Mackie and Schroeder, 2019) instead. However, all simulations using this other bedrock realization result in full retreat. We do note that this alternative DEM, however, falls outside of BedMachine's bedrock error range (both above its maximum limit and below its minimum limit) and would therefore not be correctly configured for the sensitivity experiments we perform.

Other stress balance models of higher sensitivity to bedrock perturbations may show differing results when compared to simulations that solve the SSA equations. This is due to SSA approximation solutions being less sensitive to small bedrock perturbations than other higher-order stress balance solutions. Therefore, we perhaps may be overestimating our resulting bedrock uncertainty ranges as other stress balance models may produce larger SLR mitigation effects from various bedrock perturbations seen throughout this study.

With respect to biases in our model setup, results could vary based on how the 10-year relaxation period of the control model is performed. Though the relaxation period we run brings the grounding line fairly close to the observed present day Thwaites grounding line, small but potential pinning points may already be passed by the grounding line. The converse too is possible: potential pinning points may not yet be passed by the model grounding line, while the observed grounding line is actually already upstream of these pinning points.

We choose to not conduct a new friction inversion for every bedrock configuration due to the high computational price addition to our model. Parameters for basal slipperiness and basal sliding laws are also held constant and not fully investigated in this study as we seek to investigate the effect of the isolated bedrock bumpiness without adding the complexity of other changes. This ensures that all ice responses are forced only by the perturbations and noisiness of the basal topography and that even small sensitivities can be robustly calculated against our control simulation. It is important to note, however, that conducting an inversion for basal sliding for each topographical realization could alter our results as modeled stress balance (Kyrke-Smith et al., 2018) and, consequently, modeled grounding line stability are sensitive to basal drag. This is especially the case at pinning points, such as the current grounding line position and before the first set of ridges (Fig. S3, ridge set a). As basal slipperiness and basal sliding laws are not fully investigated in this study, we leave this open question for future investigation as to how much of an impact neglecting the transient friction recalculation may bias our model results.

By design, we use a simplified but highly aggressive basal melt rate approximation in order to capture a wide range of retreat scenarios. Alternative forms of basal melt approximations (interpolation techniques, bedrock altitude dependents, etc.) have been observed to change ice-sheet model projections and can affect grounding line migration (Favier et al., 2019). Furthermore, alternative sliding laws and sliding coefficients can also influence results depending on the approximation made (e.g., Bulthuis et al., 2019; Alevropoulos-Borrill et al., 2020).

7 Conclusion

We conclude that within the currently reported bedrock topography error, it is possible for a simulation of Thwaites Glacier to remain stable despite aggressive basal melting rates. In particular, we find that positive bedrock perturbations equivalent to the error itself applied throughout the basin promote stability within the 200-year time period examined (Sect. 3). This suggests that unknown bedrock features within the basin could render Thwaites glacier more stable than would be estimated by a model using BedMachine topography to define its geometry. Finely resolved bedrock features that may not be captured by current observations should not be neglected in forward simulations. Through the wavelet amplification techniques detailed above (Sect. 4.1), we determine that accurately resolving bedrock at 2 km and 8 m respectively for spatial and vertical resolutions, especially in proximity to pinning point locations (Sect. 4), would be required to constrain our model uncertainty in 200-year extreme warming projections to within ±2 cm of SLR contribution. We find that the most influential pinning points are located near the present-day grounding line (Sect. 5). Therefore, in order to accommodate finely resolved observations of bedrock topography observations, model meshes should also be at least as fine as 2 km in areas where the grounding line of the glacier may migrate. Additionally, future bedrock data should aim to limit bedrock uncertainty to ±8 m near the present-day grounding line. Lastly, we conclude that the location and amplitude of a bedrock feature is more important than its sharpness and shape (Sect. 5).

As Thwaites Glacier is unstable, sensitive, and potentially may contribute one of the largest SLRs from the Antarctic Ice Sheet, our resolution requirements (Sect. 4) may represent a lower bound throughout the continent for bedrock data measurement; more studies may be required to confirm this assumption. This is particularly true since Thwaites is highly sensitive to retreat in comparison to the rest of the continent.

Our experiments suggest that Thwaites Glacier requires precise data measurements in order to obtain a full projection of its stability. With the glacier being particularly sensitive to bedrock perturbation amplitude and location, it is critical to constrain how much of the uncertainty in simulated SLR contribution is sourced from bedrock topography. Even small perturbations in highly influential locations can have large implications for ice-sheet model estimates of grounding line evolution and projections of future regional contribution to SLR. Overall, improving uncertainty in measurements of bedrock topography will improve certainty in projections of the future evolution of the entire Antarctic Ice Sheet.

Appendix A

A1 Initialization algorithm

Due to the perturbed bedrock geometry, a series of initialization steps occurs before all runs, and we adjust the bed and ice thickness to ensure model stability and avoid instantaneous shock upon restart of the forward model. These steps are as follows:

1. Locations at which the ice sheet is ungrounded do not sustain bedrock perturbations to prevent a sudden spontaneous ice-sheet advance.

2. The new bedrock elevation is set as the base height for all previously grounded features.

3. Anywhere that the new base goes above the surface has a recalculation to make the surface at least 1 m greater than the new base (i.e., the model minimum ice thickness of 1 m exists everywhere).

4. We force all grounded areas to be above hydrostatic equilibrium by adding an appropriate amount of ice to unstable element surfaces (allowing a momentary stability to occur within the first few time steps).

5. Locations that have an ice surface less than 1 m are set equal to 1 m.

6. Thickness is calculated to be the new surface subtracted from the new base (we implement in this fashion due to ice thickness often being found as a function of known bedrock depth).

A2 Figure code combination legend

Figures 3 and 4 use wavelet amplification techniques to create perturbations of various shapes and amplitudes. The naming techniques in the legends are as follows: “Wvlt: 〈Wavelet name〉| Lvl: 〈resolution level〉| %: 〈additional percent amplification〉”. The wavelet name is either db2 (Daubechies' second-order wavelet) or dmey (discrete Meyer wavelet), the resolution level (varying integer between 1 and 4) is the amount of wavelet decompositions that occur on the low-frequency or original image before high-frequency amplification occurs, and the additional percent amplification is the amplification multiplier added onto a base multiplier of 1.

A3 Tabulated SLR values for Experiment 1 (Fig. 2): minimum and maximum bedrock resulting spreads

Table A1 The SLR contribution spread between the minimum, maximum, and control simulations for Experiment 1 (Fig. 2). SLR contribution for each simulation is calculated by converting the change in the volume of ice above floatation in the domain between simulation years 150 or 200 and time 0 to its equivalent contribution to global mean sea level.