Inverting ice surface elevation and velocity for bed topography and slipperiness beneath Thwaites Glacier
- School of GeoSciences, University of Edinburgh, Drummond St., Edinburgh, EH8 9XP, UK
Correspondence: Helen Ockenden (firstname.lastname@example.org)
There is significant uncertainty over how ice sheets and glaciers will respond to rising global temperatures. Limited knowledge of the topography and rheology of the ice–bed interface is a key cause of this uncertainty as models show that small changes in the bed can have a large influence on predicted rates of ice loss. Most of our detailed knowledge of bed topography comes from airborne and ground-penetrating radar observations. However, these direct observations are not spaced closely enough to meet the requirements of ice-sheet models, so interpolation and inversion methods are used to fill in the gaps. Here we present the results of a new inversion of surface elevation and velocity data over Thwaites Glacier, West Antarctica, for bed topography and slipperiness (i.e. the degree of basal slip for a given level of drag). The inversion is based on a steady-state linear perturbation analysis of the shallow-ice-stream equations. The method works by identifying disturbances to surface flow which are caused by obstacles or sticky patches in the bed and can therefore be applied wherever the shallow-ice-stream equations hold and where surface data are available, even where the ice thickness is not well known. We assess the performance of the inversion for topography with the available radar data. Although the topographic output from the inversion is less successful where the bed slopes steeply, it compares well with radar data from the central trunk of the glacier for medium-wavelength features (5–50 km). This method could therefore be useful as an independent test of other interpolation methods such as mass conservation and kriging. We do not have data to allow us to assess the success of the slipperiness results from our inversions, but we provide maps that may guide future seismic data collection across Thwaites Glacier. The methods presented here show significant promise for using high-resolution satellite datasets, calibrated by sparser field datasets, to generate high-resolution bed topography products across the ice sheets and therefore contribute to reduced uncertainty in predictions of future sea-level rise.
Predicting the rate at which marine sectors of the West Antarctic Ice Sheet will retreat and contribute to globally rising sea levels is of increasing importance due to persistent climate forcing across the region over the last few decades (Scambos et al., 2017; Turner et al., 2017). Ice-sheet modelling studies emphasise the role of bed topography and rheology in understanding future ice loss (Durand et al., 2011; Parizek et al., 2013; Sun et al., 2014; Nias et al., 2016, 2018; Kyrke-Smith et al., 2018; Yu et al., 2018; Koellner et al., 2019). Bed topography is particularly important for marine-terminating glaciers, such as Thwaites Glacier in West Antarctica, which are vulnerable to marine ice-sheet instability (Weertman, 1974; Hughes, 1981; Schoof, 2007; Goldberg et al., 2009; Gudmundsson, 2013). However, bed topography constrained by geophysical surveying at the resolutions required for ice-sheet modelling (Durand et al., 2011; McCormack et al., 2018) is rarely available, so projections of future ice-sheet behaviour have to rely on bed topographies interpolated in a variety of ways between the direct measurements (Vaughan et al., 2006; Fretwell et al., 2013; Rignot et al., 2014; Millan et al., 2017; Morlighem et al., 2020). Over Thwaites Glacier, these interpolations have typically infilled areas of 15 km by 15 km between aerogeophysical flight lines, but 15 km between observations is much coarser than the resolution which Durand et al. (2011) suggest is desirable.
Where ice-penetrating radar surveys have been undertaken with sub-ice-thickness line spacing (e.g. Rutford Ice Stream – King et al., 2016; Pine Island Glacier – Bingham et al., 2017; Thwaites Glacier – Holschuh et al., 2020), they clearly identify details which are important for studying future ice-sheet behaviour that are not present in the interpolated bed topography products. In particular, imaged signatures in the bed often show some similarity to the much subtler topography of the ice surface above them. Theoretical studies based on linear perturbation theory (Gudmundsson, 2003; Raymond, 2005; Gudmundsson, 2008; Gudmundsson and Raymond, 2008; Raymond and Gudmundsson, 2009) have explored the relationship between the bed and the surface. The resulting relations can be used to infer bed properties from those of the surface but have only been applied twice to realistic settings: on 2D surface data from MacAyeal Ice Stream (Thorsteinsson et al., 2003) and on a 1D flow line from Rutford Ice Stream (Pralong and Gudmundsson, 2011). Both studies were undertaken in an era when surface elevation observations over Antarctica were of much lower quality and resolution than they are today.
Bed conditions such as geology, hydrology and sediment distribution also play a role in controlling ice flow and behaviour (Durand et al., 2011; Koellner et al., 2019) and are often poorly constrained. In many ice-sheet models, these bed conditions are combined into one parameter known as slipperiness, which is a measure of how easily the ice can slide over the topography (Rignot et al., 2011). Some seismic lines have been collected on Thwaites Glacier (Muto et al., 2019a, b), allowing a brief glimpse into the sediment distribution. Over the whole glacier, however, there are very few direct measurements of bed conditions which can be combined into slipperiness.
In this paper we exploit the relatively new availability of high-resolution surface elevation (REMA, ∼ 8 m; Howat et al., 2019) and velocity (NASA ITS_LIVE, ∼ 120 m resolution; Gardner et al., 2018) datasets. We apply linear perturbation theories to explore bed topography and slipperiness across the Thwaites Glacier catchment. We use a steady-state version of the shallow-ice-stream equations presented by Gudmundsson (2008) and compare the topography output from the inversion to radar grids and flight lines to assess its performance.
2.1 Derivation of the steady-state shallow-ice-stream transfer functions
Gudmundsson (2008) derived a set of transfer functions which describe the relationship between the time-variant Fourier transforms of bed topography (), bed slipperiness (), surface topography () and horizontal components of surface velocity (). For the purposes considered here, this derivation can be simplified by considering the steady state from the beginning, removing the need to do a Laplace transform. Other than this, the derivation largely follows that of Gudmundsson (2008), but for clarity we state key assumptions and results here.
2.1.1 Response of flow to basal topography perturbations
Following Gudmundsson (2008) and working in a coordinate system tilted forward in the x direction by the mean surface slope, α, we start with the shallow-ice-stream equations of motion (MacAyeal, 1989),
where u and v are the depth-independent velocity components in the x and y directions, respectively; h is the ice thickness; η is the effective ice viscosity; c is the basal slipperiness; m is a sliding-law parameter; ρ is the ice density; g is the acceleration due to gravity; s is the ice surface elevation; and α is the mean ice surface slope in the x direction.
Assuming that ice is a linear viscous medium (n=1) and that there is a non-linear sliding law (m>0), then the shallow-ice-stream equations can be linearised and solved analytically. We consider the spatial response to a small perturbation in basal topography, b, linearising around a reference model () with , , , , v=Δv, w=Δw and . The zero-order solutions are spatially constant, representing uniform flow down an inclined plane.
We, however, are interested in the first-order momentum balance equations:
Also to the first order and importantly in the steady state, we have the following upper and lower kinematic boundary conditions:
Various points about the validity of the steady-state assumption for Thwaites Glacier are raised in the discussion (Sect. 4).
All variables are then Fourier transformed with respect to the spatial variables x and y. Fourier-transformed variables are denoted with a circumflex ( ). In the forward Fourier transform the wavenumbers in the x and y directions are denoted by k and l, respectively. This Fourier transform gives
From depth integration of the Fourier-transformed incompressibility condition , we have
which, along with the steady-state boundary conditions, yields
which represent the ratio of the Fourier components of the surface to the Fourier components of the bed as a function of wavenumber. The following abbreviations are used for simplicity in the derivation: , , , , , , and .
2.1.2 Response of flow to basal slipperiness perturbation
Starting once again with the shallow-ice-stream equations (Eqs. 1 and 2; MacAyeal, 1989), this time we consider the response to a small perturbation in basal slipperiness, c, linearising with , , , , v=Δv, w=Δw and where Δc is the fractional slipperiness.
This gives the following first-order momentum balance equations:
Fourier transforming with respect to the spatial variables x and y gives
As there is no bed topography perturbation, the steady-state boundary conditions become
which represent the ratio of variability in the Fourier components of the surface to variability in the Fourier components of the slipperiness.
Note that the transfer functions Tuc and Tvc are not the same as the steady-state versions of the transfer functions published in Gudmundsson (2008) as there is a typographic error in their paper. However, when plotted graphically, they can be used to reproduce the figures in that paper.
These transfer functions can also be considered in a non-dimensional form, allowing us to make more general statements about the behaviour of the system in terms of key variables, such as ice thickness as the characteristic length scale. For this purpose the same scalings as used in Gudmundsson (2003, 2008) are employed. All spatial scales are in units of mean ice thickness (), and stress components are in units of driving stress (τd). Non-dimensional velocity components are in units of mean deformational velocity () where
The scale for slipperiness is given by , where is the mean dimensional slipperiness and is the mean non-dimensional slipperiness. depends on not only slipperiness but also viscosity and thickness (through its dependence on us and τd). From Gudmundsson (2008) we know that and we also have . The ice surface velocity is the sum of the deformational velocity and the basal velocity, such that . With some simple algebra, we can therefore express the scale for slipperiness in terms of the surface velocity , which is a known quantity,
The non-dimensional form of the equations is obtained using the substitutions , , , , and . The non-dimensional transfer functions can be found in the Supplement and are also shown graphically there. Non-dimensionalised parameters are represented by capital letters (B, C, S, U, V, TSB, TUB, TVB, TSC, TUC and TUV).
2.2 The inverse problem
The non-dimensional transfer functions (TSB, TUB, TVB, TSC, TUC and TVC) describe the relationship between the Fourier transforms of non-dimensionalised bed topography (), bed slipperiness (), surface topography () and surface velocity (), as functions of the wavenumbers k and l. If the bed topography and slipperiness are known, then surface topography and velocity components are given by the forward model:
For each set of wavenumbers in Fourier space (k and l), we have three known variables (, and ) and two unknowns ( and ), so the system is over-determined. We can therefore solve these equations independently for each wavenumber component of non-dimensional bed topography and slipperiness using a weighted least-squares inversion of Eqs. (26)–(28). Short-wavelength features and features aligned with ice flow are problematic because they cause flow disturbances in the ice which do not reach the surface in a measurable way, and so they can not be inverted from the surface data. We include a filtering parameter (pfilt) to remove these problematic wavelength components such that increasing the filtering parameter (pfilt) removes progressively longer wavelength features. This method was first used by Thorsteinsson et al. (2003) in their study of MacAyeal Ice Stream (formerly Ice Stream E). The equations which solve Eqs. (26)–(28) are therefore not repeated here in the main text but are given in notation consistent with this paper in Appendix C.
2.3 Synthetic tests
Synthetic tests allow us to explore which bed features can and can not be resolved using this inversion method. First we create a synthetic bed topography (b) and slipperiness (c) on a 50 km by 50 km grid, with 120 m data spacing, purposely matching the data spacing of the ITS_LIVE velocity product (Gardner et al., 2018). We subtract the mean bed elevation, slope and slipperiness so that the bed varies about 0, and we taper the outer 5 km of the grid linearly to 0 at the edges. This reduces edge effects in the inversion because the Fourier transform requires a periodic domain. We tested a few other sensible tapering functions, including semi-sinusoids but observed negligible differences in the inversion output when compared to the linear function. We then non-dimensionalise the tapered bed using the length scales given in Sect. 2.1.3 and Fourier transform to obtain and . The non-dimensional surface elevation () and the velocity components (, ) are calculated using the forward model and dimensionalised using the length scales given in Sect. 2.1.3. To simulate measurement errors in the real surface data, we add random noise to the generated surface (s, u, v). This noise is white noise with a Gaussian low-pass filter applied in Fourier space to give it a non-random frequency distribution. We then taper, non-dimensionalise and Fourier transform the noise-added surface data. Finally, we invert for non-dimensional bed topography () and slipperiness () using the inversion procedure described in Sect. 2.2 and the Supplement. After dimensionalisation, the inverted bed can be compared to the synthetic bed to study the behaviour of the inversion.
2.3.1 Parameter value choices
When running synthetic tests, several model parameters can be varied, in addition to the synthetic bed topography and slipperiness. Following Gudmundsson (2008) and Thorsteinsson et al. (2003), the sliding-law constant was set to m=1 and the filtering parameter (Eq. C4) to . The mean ice thickness , mean surface slope α and mean ice surface velocity depend on the region studied, but for these synthetic tests, they are set at = 2000 m, α = 0.002 and = 100 m s−1, values thought to be appropriate for the Thwaites Glacier region (Gudmundsson, 2008; Howat et al., 2019; Morlighem et al., 2020; Gardner et al., 2018). After studying the results of the Thwaites Glacier inversion for a variety of values of , the non-dimensional mean slipperiness, we set = 100. When applied to Eq. (25), these values give a dimensional mean slipperiness = 2.7 × 10−3 . We set the weighting factors (Eq. C2) to be Σs=0.001, Σu=1 and Σv=1. This accounts for the mismatch in the relative magnitude of the non-dimensionalised surface elevation and velocity and means the least-squares inversion solves for all three factors equally.
2.3.2 Resolution of bed forms
A two-dimensional Fourier transform decomposes an image into a weighted sum of two-dimensional sinusoidal basis functions. For this reason, all of our synthetic tests used sinusoidal bed topographies and slipperiness as these are the most illustrative of the capabilities of the inversion. Sinusoidal basis functions vary depending on three parameters: the wavenumbers in the x and y directions (k,l) and the weighting or amplitude of the sinusoid. However, rather than considering wavenumbers, it is more intuitive to consider the horizontal wavelength, λ, and angle, θ, to the direction of the flow, where , , k=jcos(θ) and l=jsin(θ).
Figures 1–3 show how well bed topography and slipperiness can be resolved by the inversion when the angle to flow, θ; wavelength, λ; and amplitude are varied. The amplitude of the noise added to the synthetically generated surface is ± 2 m for surface elevation and ± 15 m s−1 for velocity components as these are at the upper limit of the errors for REMA (Howat et al., 2019) and ITS_LIVE (Gardner et al., 2018) in the Thwaites Glacier region.
These synthetic tests show that in this simple least-squares inversion, the bed can be well resolved if the angle to the flow is greater than 15∘, the wavelength is more than 2000 m, and the amplitude is greater than 10 m for topography or 1.34 × 10−4 for variability in slipperiness. It is worth noting, however, that the resolution of wavelengths varies depending on the ice thickness, which is the non-dimensional scale factor for lengths. This means that ice thickness is directly proportional to the wavelength at which variations in bed topography and slipperiness should be resolvable.
2.4 Applying the inversion to real data
We now turn our attention to the methodology used to apply the synthetically tested inversion to real data, using the Thwaites Glacier catchment as our example.
Our base data for surface elevation and velocity were the REMA digital elevation model with 8 m resolution (Howat et al., 2019) and output from the NASA MEaSUREs ITS_LIVE project with 120 m resolution (Gardner et al., 2018), respectively. Based on the latter, we therefore inverted for bed topography and slipperiness at 120 m resolution. An estimate of the ice thickness in each 50 km by 50 km region is obtained from a 50 km averaged version of BedMachine Antarctica ice thickness (Morlighem et al., 2020), and this is the only prior information about ice thickness used in the inversion.
In the synthetic tests discussed above, the inversion was run over a single 50 km by 50 km grid, with the outer 10 % of the grid discarded to reduce edge effects introduced during the Fourier transform. To look at the whole Thwaites catchment, we could use a set of adjacent 50 km by 50 km grids, but instead we chose to use more densely distributed grids which overlap. For each grid point, we calculate nine different ((but overlapping) inverted beds and then the mean bed topography and standard deviation. The standard deviation is not a measure of the error, but since the main approximation in the physics is the linearisation, we interpret the standard deviation to be a measure of non-linearity. In each of the overlapping grids, we use a set of zero-order parameters (such as average ice thickness), and because these zero-order parameters vary between the grids, the linearisation is different and the resulting beds are also different. However, inappropriate application of the shallow-ice-stream approximation or edge effects could also be influencing this. Figure 4a summarises this methodology. The bed topography and slipperiness results presented here are the grid-point by grid-point means of nine overlapping grids where each overlapping region is 1.67 km by 1.67 km. These values were chosen following tests on a small region of the real data (Fig. 4b).
When applied to real surface elevation and velocity data, this method generates four products: the mean and standard deviation of bed topography and the mean and standard deviation of bed slipperiness. The standard deviation is not a measure of the error in the bed topography or bed slipperiness and should not be interpreted as such.
The bed topography product from the inversion is shown in Fig. 5b. On the basin scale, the main basal topographic features identified by the inversion are several sets of parallel ridges which are oriented perpendicularly to the direction of ice flow and smooth basins in between these ridges. The location of these ridges matches well with the BedMachine Antarctica bed (Fig. 5a), particularly around the subglacial lakes, which appear to be between successive sets of ridges. The smoother topography in the basins between ridges is reflected in the inversion, particularly in the basin to the east of the Upper Thwaites region (Basin Y, Fig. 5). Many smaller hills also match BedMachine Antarctica, such as those at the upstream (south) end of the Upper Thwaites radar grid. However, the inversion also generates some notable features that are not present in the BedMachine Antarctica bed, such as the north-eastern extent of the central ridge next to the most upstream subglacial lake (Ridge Z, Fig. 5b).
The standard deviation of the bed topography (Fig. 5c) represents the range in model outputs from overlapping grid squares which use different regions of the ice surface. As might be expected, this standard deviation is lower in the central trunk of the glacier where the topographic gradients are smaller. The standard deviation is higher at the edges of the glacier trunk where the gradient of the topography changes and the shallow-ice-stream approximation breaks down. Standard deviation is also high in the north-west part of the inversion where some of the input surface is the surface of the floating ice shelf, which is subject to different physical processes than the grounded ice. The inversion does not produce results for the north-west corner as there are gaps in the input data where the surface sampled is not ice but open ocean.
The bed slipperiness product from the inversion is shown in Fig. 5d. The pattern in the slipperiness output is similar to the topography, with a dominant east-to-west lineation, although it is slightly difficult to make out due to the strong underlying slipperiness variation across the Thwaites Glacier region. This directional trend in slipperiness is also observed in the slipperiness output of other inversions in the Thwaites Glacier region (Barnes et al., 2021). There is less variability in slipperiness in the basins where the topography is smoother. In contrast the ridges have much more variable slipperiness. Like the standard deviation of bed topography, the standard deviation of bed slipperiness (Fig. 5e) is highest around the edges of the central trunk where there are higher topographic gradients.
3.1 Comparison to radar data
In order to assess the success of the inversion, the bed topography output can be compared to existing radar data. Figure 6 shows a comparison between the inverted bed topography and bed topography sounded by swath radar at sub-ice-thickness resolution across two 20 km by 40 km regions (Holschuh et al., 2020). The inverted bed shows a good match to the swath-radar-imaged bed at larger scales, picking out the locations of all the main hills and valleys. There is a better match for the Upper Thwaites region than the Lower Thwaites region, and the fact that the inversion detects the channel between the ridges in the downstream (left) part of Upper Thwaites is particularly encouraging.
We further compare the inverted bed topography with the bed topography sounded by swath radar along ice flow (Fig. 7) and by airborne radar across ice flow (Fig. 8) over Thwaites Glacier in the 2019/20 field season (Jordan and Robinson, 2021). Further comparisons to more radar flight lines collected in the same field season can be seen in the Supplement. These figures also demonstrate that the inversion performs well in detecting the main hills and valleys and also highlight that the amplitude of topographic features are not always resolved correctly. This is likely due to variability in the local mean slipperiness away from the imposed global value of non-dimensional slipperiness . If a non-dimensional slipperiness is imposed (Figs. 6, 7c and 8c), then the amplitudes of the inverted topography are reduced. Sometimes there is also an offset between the inverted and radar-sounded beds, caused by using a 50 km averaged version of the BedMachine ice thickness as the inversion ice thickness, rather than more detailed prior information.
Our results demonstrate significant promise for being able to invert for bed topography across parts of Antarctica and other polar regions from surface elevation and velocity datasets. Comparisons with existing radar data available from Thwaites Glacier suggest that within the central trunk of the glacier, the bed features identified by the inversion are normally in the correct locations but are not always centred around the correct depth. These average depth differences are primarily due to the mean ice thickness used in the inversion, which is a 50 km averaged version of the BedMachine Antarctica ice thickness (Morlighem et al., 2020). For regions where there are radar flight lines and grids, these radar observations could be used instead of the averaged ice thickness to ensure that the bed depth is correct. For regions where there are very few existing radar data, this method has the potential in the future to identify obstacles to flow which are significant enough to affect surface ice dynamics, even if there is uncertainty about the local or regional ice thickness.
The regions where the inverted bed deviates significantly from the topography picked from radar surveys are of particular interest in assessing the potential of our inversion. Differences between inverted topography and radar lines are likely to be due to physical processes which are not encapsulated by the shallow-ice-stream approximation. One place in which the shallow-ice-stream approximation is known to break down is where the mean slope of the bedrock becomes too steep (Gudmundsson, 2003; Le Meur et al., 2004). The effect of this can be observed around the edges of the central trunk of the glacier, where the topographic slope is steep and the match between the inverted bed and the radar lines is poor (Fig. 5). Gudmundsson (2003) derived a full-system non-hydrostatic momentum balance version of the transfer function used in this work. The full-system approach does not rely on the shallow-ice-stream approximation and should therefore perform better. Future work using this method is likely to incorporate these more complex equations.
A further consideration in comparing the inverted bed to real data is the steady-state assumption made when deriving the transfer functions. Without repeat radar measurements for Thwaites Glacier we can not be sure of the stability of the bed. If the bed beneath Thwaites Glacier is changing rapidly, as observed at Rutford Ice Stream (Smith et al., 2007), then the surface may not represent the current bed but some long-term average. However, observations at Pine Island Glacier (Brisbourne et al., 2017; Davies et al., 2018) suggest that the bed is not changing rapidly there, and it is possible that neighbouring Thwaites Glacier might be behaving similarly. Additionally, the erosion observed at Rutford Ice Stream does not significantly change the shape of the topography on the wavelengths resolved by this inversion. If drumlins or mega-scale glaciation lineations (MSGLs) were forming, we would not be able to detect them with this method, as landforms aligned to flow fall in the null space of the inversion.
The steady-state assumption applies not only to the bed but also to the ice surface. Ice surface lowering due to glacier thinning would also affect the steady-state assumption, but since generally the ice surface lowers in a relatively uniform way, this would not have a significant effect on the first-order variations in the ice surface or the results of the model. More significant would be changes in the ice surface due to the filling and draining of subglacial lakes, but these changes are normally fairly localised and would not propagate to the higher-wavelength Fourier components. For Thwaites Glacier, the location of subglacial lakes is relatively well known (Smith et al., 2017; Hoffman et al., 2020; Malczyk et al., 2020), and we predict troughs in these locations as expected. The ice surface also becomes more unstable closer to the grounding line, with increased crevassing which would affect the surface profile. However, since results in the region immediately adjacent to the grounding line are compromised by the different physics of the ice shelf anyway, this is not a significant concern. With these caveats, we therefore consider the steady-state assumption to be suitable for the purposes of this inversion.
If the steady-state assumption is valid, then the age of the datasets used in the inversion is not important. However, input data from different years or decades could also affect the steady-state assumption. The main surface expressions of known bed features appear to be fairly similar between REMA (Howat et al., 2019) (2008–2018) and the earlier Bamber DEM (Bamber et al., 2009) (2003–2008), supporting the validity of the steady-state assumption for Thwaites Glacier. However, we also note that non-steady-state changes in the ice surface may be the reason for some of the features we observe (such as Ridge Z, Fig. 5) in the inversion output which are not seen in the airborne flight lines.
As in any modelling study, it is important to explore the behaviour of the inversion when the parameters chosen are varied. In this inversion there are just four fixed parameters which are not derived from the input data: the sliding-law exponent, m; the filtering parameter, pfilt; the weighting factor, Σs; and the mean slipperiness, . The filtering parameter pfilt controls which frequencies are suppressed in the inversion to avoid introducing singularities. Higher values of pfilt (closer to 0) will filter out higher frequencies (lower wavelengths), and so a value of is chosen to filter out noisy short wavelengths while maintaining realistic bed features. The weighting factor Σs controls the balance in the inversion between the surface elevation and surface velocity data, with smaller values of Σs weighting the inversion towards the surface data. Varying pfilt and Σs for the inversion of the real surface data confirms the choice of values from the synthetic tests (, Σs=0.001) as sensible values which return the best match with real bed data (Figs. S3 and S4 in the Supplement).
There is less certainty over what is the most suitable value for the non-dimensional mean slipperiness . Although we have some measurements of the bed properties of Thwaites Glacier from seismic lines, gravity and magnetic inversions (Diehl, 2008; Jordan et al., 2010; Muto et al., 2019a), these are spatially limited and it is not currently clear how these properties combine into slipperiness at the bed (Kyrke-Smith et al., 2017). If is higher, then the amplitude of bed variability in the inversion output falls. Given the geological variability likely to be associated with multiple rifted tectonic blocks (Dunham et al., 2020) and the sediments deposited in those rift basins (Muto et al., 2019a, b), it is unlikely that the mean slipperiness, , is the same across the whole region modelled here. Modelling studies which compare the results of different inversion procedures show that slipperiness may be quite variable across the Thwaites Glacier catchment (Barnes et al., 2021). In addition we note that the trend is quite different from features observed in the inverted topography, showing that the slipperiness map is not a result of linear trade-offs with topography in the inversion solution. The three-dimensional radar grids (Holschuh et al., 2020) are both located within regions with more topographic variability, likely unlifted rift blocks. This may explain why a lower value of gives the best match in these regions, whereas a higher value of (more slippery) gives a better match for the radar lines which cover both lithologies. It may also be that the 3D grids, which contain both along- and across-flow variability, are more representative of the bed than the 2D radar lines. For this reason, the catchment-scale bed topography presented in Fig. 5 is from the inversion with .
This uncertainty in means that some prior radar information is useful in order to calibrate the inversion method to give the best results. However, changing only alters the amplitude, so even if there is no prior information, the inversion will still identify hills and troughs. More detailed analysis of seismic and gravity data alongside the results of the inversion could also reveal trends that could be useful when applying this technique elsewhere.
The high-resolution swath radar grids (Holschuh et al., 2020) presented here have already been included in BedMachine Antarctica. Using the radar grids which are currently available, we can not therefore explore how well this inversion method performs compared to BedMachine Antarctica (Morlighem et al., 2020) over a dense radar grid. Both techniques use ice surface elevation and velocity datasets. BedMachine Antarctica uses these datasets and the principles of mass conservation (or streamline diffusion in slow-moving areas) to interpolate between detailed prior information on ice thickness from existing radar measurements. In contrast, the inversion method only requires an estimate of average ice thickness for each 50 km by 50 km grid and uses the linear perturbation theory described. This 50 km averaged ice thickness is subtracted from the ice surface to provide a reference bed to which the inversion adds perturbations. However, even if a single ice-thickness () value is used for the entire catchment, then the inversion method will still identify the location of hills and troughs in the bed topography, although the amplitudes and absolute depths of these features may be affected. Since the mass conservation method used in BedMachine assesses the ice flow through a series of flux gates ideally constrained by topography, it is much more reliant on good-quality, closely spaced, ice-thickness measurements from radar systems.
To explore the role that the 50 km averaged ice thickness plays in the results of the inversion, we computed a 50 km gridded version of the BedMachine Antarctica ice thickness (Fig. 9) and then carried out a new inversion over the Lower Thwaites region (where there is existing swath radar). In this alternate ice-thickness input, each 50 km by 50 km region contains only one ice-thickness value, which is the average over that region. The results of this re-run (Fig. 9) with reduced prior ice-thickness information show the same short-wavelength features as the results using the full ice-thickness input, illustrating that the inversion method presented here is not unduly influenced by the ice thickness derived from BedMachine Antarctica.
In Figs. 7 and 8, it appears that the method estimates shorter-wavelength topography more accurately than longer wavelengths. We demonstrate this in Fig. 10, which shows results after wavelengths greater than 50 km have been removed from all profiles. It is clear that the inversion identifies peaks and troughs in the bed, although the amplitude of these features is not always correct. Fourier components with a wavelength above 50 km mainly represent the prior ice-thickness information supplied to the inversion as this is the large-scale zero-order topography to which the first-order perturbations from the inversion are added. The greater match between the results of the inversion and the PASIN data (Jordan and Robinson, 2021) after this bandpass filter therefore provides further evidence that the ice thickness derived from BedMachine Antarctica does not influence the key results from the inversion.
A comparison of the two methods over an area where radar data have not yet been incorporated into BedMachine would allow an assessment of the reliability of the two techniques and identification of any artificial bed features introduced by each. Since the two radar grids presented (Holschuh et al., 2020) were included in the derivation of BedMachine Antarctica, no independent test is possible until more radar grids are collected.
We present the method and results of an inversion of ice surface elevation and velocity for bed topography and slipperiness in the Thwaites Glacier region. Our method builds on the method used by Thorsteinsson et al. (2003) in their study of MacAyeal Ice Stream but is based on a steady-state linear perturbation analysis of the shallow-ice-stream equations (MacAyeal, 1989; Gudmundsson, 2003). Synthetic tests show that this method can resolve variability in bed topography and slipperiness on wavelengths greater than 1 ice thickness and at amplitudes of more than 10 m for topography or 1 × 10−4 for slipperiness, as long as the variability is not aligned with the ice flow direction. Comparison of the results of the inversion with radar grids and flight lines suggests that the inversion correctly identifies most short-wavelength (< 50 km horizontal) features in the bed, with the correlation coefficient of a linear regression between the inverted bed and the radar bed as high as r=0.93 along some flight lines. This method works best in the central trunk of the glacier, where the gradient of the long-wavelength topography is low and relatively constant. Mismatches between the inverted bed topography and radar measurements are probably due to one of three factors: an incorrect ice thickness for that region, an unusually sticky or slippery bed, or physical processes not accounted for by the steady-state linearised shallow-ice-stream approximation. Future work, including incorporating more local prior ice-thickness data from radar measurements and the non-hydrostatic transfer functions from Gudmundsson (2003), may help to reduce these mismatches. Overall, the inversion provides an additional tool for studying landforms in the bed beneath glaciers which have a significant impact on ice flow. It will be particularly useful in ice streams where radar flight lines are sparse and standard interpolation techniques struggle, potentially reducing uncertainties in modelling the future behaviour of those regions and their contributions to global sea-level rise.
Starting from the shallow-ice-stream equations, linearising around a small perturbation in bed topography and taking the Fourier transform of the first-order equations, we have
The determinant of the left-hand side of these equations is
where the following abbreviations have been used to make the algebra clearer to follow: , and .
Therefore we have
which simplifies to
We then have
In the steady state, the kinematic boundary condition is
Substituting the expression from above gives
which can be rearranged as
Expanding the expression for gives
remembering that .
Expanding the expression for gives
remembering that .
Starting once again with the shallow-ice-stream equations (MacAyeal, 1989), this time we consider the response to a small perturbation in basal slipperiness, linearising and taking the Fourier transform to give
The determinant of the left-hand side is
where the following abbreviations have been used to make the algebra easier to follow: , and .
Therefore we have
which simplifies to
We then have
At a steady state, , so we have the following boundary condition:
Expanding the expression for (Eq. B1) gives
remembering that , and leads to the following steady-state transfer function:
Expanding the expression for (Eq. B2) gives
remembering that , which leads to the following steady-state transfer function:
The transfer functions (Tsb, Tub, Tvb, Tsc, Tuc and Tvc) describe the relationship between the Fourier transforms of bed topography (), bed slipperiness (), surface topography () and surface velocity (). If the bed topography and slipperiness are known, then surface topography and velocity components are given by the forward model:
Non-dimensionalised this gives
Since the system is over-determined, we can use a weighted least-squares inversion of Eqs. (26)–(28) to find the bed topography and slipperiness which are the most consistent with the ice surface. This is the same method used by Thorsteinsson et al. (2003) in their study of MacAyeal Ice Stream (formerly Ice Stream E) but is reproduced here in notation consistent with the rest of the equations we present.
In matrix form we have the following forward model:
where , , and E = .
A least-squares inversion gives
where H is the Hermitian transpose and ∗ is the complex conjugate.
For compactness we define the following:
The least-squares solution is then
This inversion is problematic where is small or zero, which is the case for small wavelengths (when k and l are large) or for topography or slipperiness perturbations which are aligned in the direction of ice flow (k=0). Short-wavelength bed features and features aligned with ice flow are problematic because they cause flow disturbances in the ice which do not reach the surface in a measurable way. They can therefore not be inverted from the surface data.
and pfilt≤0. This filter allows through all wavelengths where is larger than P but gradually filters out other smaller wavelengths. Smaller values of P give more detail but may over-fit the surface data due to errors. Larger values of P under-fit the data and may leave out features actually represented by the data.
The filtered least-squares solution is then
The output data from the inversion are available on Zenodo at https://doi.org/10.5281/zenodo.5105687. The code for the inversion and plotting the figures is available on https://doi.org/10.5281/zenodo.5494600. The surface elevation data from REMA (Howat et al., 2019) and velocity data from ITS_LIVE (Gardner et al., 2018) used as inputs in the inversion are available freely online, as are the swath radar (Holschuh et al., 2020), airborne radar (Jordan and Robinson, 2021) and BedMachine Antarctica (Morlighem et al., 2020) bed datasets to which the results of the inversion are compared.
The supplement related to this article is available online at: https://doi.org/10.5194/tc-16-3867-2022-supplement.
RB, DG and AC developed the concept of the paper. DG advised on the use of the Gudmundsson (2008) equations; HO and DG derived the steady-state transfer functions. AC advised on the inversion methodology. HO wrote the code, adapted the input data to the most useful form and carried out the inversion. RB advised on the comparison to existing radar data. HO wrote the main body of the text, but all authors contributed to the development of the final paper and associated figures.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work is an output of the NERC E4 Doctoral Training Partnership and the Geophysical Habitat of Subglacial Thwaites (GHOST) project, a component of the International Thwaites Glacier Collaboration (ITGC).
This research has been supported by the Natural Environment Research Council (grant nos. NERC NE/S007407, NERC NE/S006672, NERC NE/S006621/1, NERC NE/S006613/1, NERC NE/S006796/1) and the National Science Foundation (grant no. NSF PLR 1738934).
This paper was edited by Ginny Catania and reviewed by three anonymous referees.
Bamber, J. L., Gomez-Dans, J. L., and Griggs, J. A.: A new 1 km digital elevation model of the Antarctic derived from combined satellite radar and laser data – Part 1: Data and methods, The Cryosphere, 3, 101–111, https://doi.org/10.5194/tc-3-101-2009, 2009. a
Barnes, J. M., Dias dos Santos, T., Goldberg, D., Gudmundsson, G. H., Morlighem, M., and De Rydt, J.: The transferability of adjoint inversion products between different ice flow models, The Cryosphere, 15, 1975–2000, https://doi.org/10.5194/tc-15-1975-2021, 2021. a, b
Bingham, R. G., Vaughan, D. G., King, E. C., Davies, D., Cornford, S. L., Smith, A. M., Arthern, R. J., Brisbourne, A. M., De Rydt, J., Graham, A. G., Spagnolo, M., Marsh, O. J., and Shean, D. E.: Diverse landscapes beneath Pine Island Glacier influence ice flow, Nat. Commun., 8, 1618, https://doi.org/10.1038/s41467-017-01597-y, 2017. a
Brisbourne, A. M., Smith, A. M., Vaughan, D. G., King, E. C., Davies, D., Bingham, R. G., Smith, E. C., Nias, I. J., and Rosier, S. H. R.: Bed conditions of Pine Island Glacier, West Antarctica, J. Geophys. Res.-Earth, 122, 419–433, https://doi.org/10.1002/2016JF004033, 2017. a
Davies, D., Bingham, R. G., King, E. C., Smith, A. M., Brisbourne, A. M., Spagnolo, M., Graham, A. G. C., Hogg, A. E., and Vaughan, D. G.: How dynamic are ice-stream beds?, The Cryosphere, 12, 1615–1628, https://doi.org/10.5194/tc-12-1615-2018, 2018. a
Diehl, T. M.: Gravity analyses for the crustal structure and subglacial geology of West Antarctica, particularly beneath Thwaites Glacier, The University of Texas at Austin ProQuest Dissertations Publishing, 2008. a
Dunham, C. K., O'Donnell, J. P., Stuart, G. W., Brisbourne, A. M., Rost, S., Jordan, T. A., Nyblade, A. A., Wiens, D. A., and Aster, R. C.: A joint inversion of receiver function and Rayleigh wave phase velocity dispersion data to estimate crustal structure in West Antarctica, Geophys. J. Int., 223, 1644–1657, https://doi.org/10.1093/gji/ggaa398, 2020. a
Durand, G., Gagliardini, O., Favier, L., Zwinger, T., and le Meur, E.: Impact of bedrock description on modeling ice sheet dynamics, Geophys. Res. Lett., 38, L20501, https://doi.org/10.1029/2011GL048892, 2011. a, b, c, d
Fretwell, P., Pritchard, H. D., Vaughan, D. G., Bamber, J. L., Barrand, N. E., Bell, R., Bianchi, C., Bingham, R. G., Blankenship, D. D., Casassa, G., Catania, G., Callens, D., Conway, H., Cook, A. J., Corr, H. F. J., Damaske, D., Damm, V., Ferraccioli, F., Forsberg, R., Fujita, S., Gim, Y., Gogineni, P., Griggs, J. A., Hindmarsh, R. C. A., Holmlund, P., Holt, J. W., Jacobel, R. W., Jenkins, A., Jokat, W., Jordan, T., King, E. C., Kohler, J., Krabill, W., Riger-Kusk, M., Langley, K. A., Leitchenkov, G., Leuschen, C., Luyendyk, B. P., Matsuoka, K., Mouginot, J., Nitsche, F. O., Nogi, Y., Nost, O. A., Popov, S. V., Rignot, E., Rippin, D. M., Rivera, A., Roberts, J., Ross, N., Siegert, M. J., Smith, A. M., Steinhage, D., Studinger, M., Sun, B., Tinto, B. K., Welch, B. C., Wilson, D., Young, D. A., Xiangbin, C., and Zirizzotti, A.: Bedmap2: improved ice bed, surface and thickness datasets for Antarctica, The Cryosphere, 7, 375–393, https://doi.org/10.5194/tc-7-375-2013, 2013. a
Gardner, A. S., Moholdt, G., Scambos, T., Fahnstock, M., Ligtenberg, S., van den Broeke, M., and Nilsson, J.: Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years, The Cryosphere, 12, 521–547, https://doi.org/10.5194/tc-12-521-2018, 2018. a, b, c, d, e, f, g, h, i
Goldberg, D., Holland, D. M., and Schoof, C.: Grounding line movement and ice shelf buttressing in marine ice sheets, J. Geophys. Res.-Earth, 114, F04026, https://doi.org/10.1029/2008JF001227, 2009. a
Gudmundsson, G. H.: Analytical solutions for the surface response to small amplitude perturbations in boundary data in the shallow-ice-stream approximation, The Cryosphere, 2, 77–93, https://doi.org/10.5194/tc-2-77-2008, 2008. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q
Gudmundsson, G. H. and Raymond, M.: On the limit to resolution and information on basal properties obtainable from surface data on ice streams, The Cryosphere, 2, 167–178, https://doi.org/10.5194/tc-2-167-2008, 2008. a
Hoffman, A. O., Christianson, K., Shapero, D., Smith, B. E., and Joughin, I.: Brief communication: Heterogenous thinning and subglacial lake activity on Thwaites Glacier, West Antarctica, The Cryosphere, 14, 4603–4609, https://doi.org/10.5194/tc-14-4603-2020, 2020. a
Holschuh, N., Christianson, K., Paden, J., Alley, R., and Anandakrishnan, S.: Linking postglacial landscapes to glacier dynamics using swath radar at Thwaites Glacier, Antarctica, Geology, 48, 268–272, https://doi.org/10.1130/G46772.1, 2020. a, b, c, d, e, f, g, h
Howat, I. M., Porter, C., Smith, B. E., Noh, M.-J., and Morin, P.: The Reference Elevation Model of Antarctica, The Cryosphere, 13, 665–674, https://doi.org/10.5194/tc-13-665-2019, 2019. a, b, c, d, e, f, g, h, i
Jordan, T. and Robinson, C.: Bed, surface elevation and ice thickness measurements derived from radar data acquired during the Thwaites Glacier airborne survey (2019/2020) (Version 1.0) [Data set], NERC EDS UK Polar Data Centre, https://doi.org/10.5285/7C12898D-7E55-458C-BA7D-ECEC8252F3B5, 2021. a, b, c, d
Jordan, T. A., Ferraccioli, F., Vaughan, D. G., Holt, J. W., Corr, H., Blankenship, D. D., and Diehl, T. M.: Aerogravity evidence for major crustal thinning under the Pine Island Glacier region (West Antarctica), Bull. Geol. Soc. Am., 122, 714–726, https://doi.org/10.1130/B26417.1, 2010. a
King, E. C., Pritchard, H. D., and Smith, A. M.: Subglacial landforms beneath Rutford Ice Stream, Antarctica: detailed bed topography from ice-penetrating radar, Earth Syst. Sci. Data, 8, 151–158, https://doi.org/10.5194/essd-8-151-2016, 2016. a
Koellner, S., Parizek, B. R., Alley, R. B., Muto, A., and Holschuh, N.: The impact of spatially-variable basal properties on outlet glacier flow, Earth Planet. Sc. Lett., 515, 200–208, https://doi.org/10.1016/J.EPSL.2019.03.026, 2019. a, b
Kyrke-Smith, T., Gudmundsson, G., and Farrell, P.: Relevance of detail in basal topography for basal slipperiness inversions: a case study on Pine Island Glacier, Antarctica, Front. Earth Sci., 6, 33, https://doi.org/10.3389/feart.2018.00033, 2018. a
Kyrke-Smith, T. M., Gudmundsson, G. H., and Farrell, P. E.: Can Seismic Observations of Bed Conditions on Ice Streams Help Constrain Parameters in Ice Flow Models?, J. Geophys. Res.-Earth, 122, 2269–2282, https://doi.org/10.1002/2017JF004373, 2017. a
Le Meur, E., Gagliardini, O., Zwinger, T., and Ruokolainen, J.: Glacier flow modelling: A comparison of the Shallow Ice Approximation and the full-Stokes solution, C. R. Phys., 5, 709–722, https://doi.org/10.1016/j.crhy.2004.10.001, 2004. a
MacAyeal, D. R.: Large-scale ice flow over a viscous basal sediment: theory and application to ice stream B, Antarctica, J. Geophys. Res., 94, 4071–4087, https://doi.org/10.1029/jb094ib04p04071, 1989. a, b, c, d
Malczyk, G., Gourmelen, N., Goldberg, D., Wuite, J., and Nagler, T.: Repeat Subglacial Lake Drainage and Filling Beneath Thwaites Glacier, Geophys. Res. Lett., 47, e2020GL089658, https://doi.org/10.1029/2020GL089658, 2020. a
McCormack, F., Galton-Fenzi, B., Seroussi, H., and Roberts, J.: The impact of bed elevation resolution on Thwaites Glacier ice dynamics, Twenty-Fifth Annual WAIS Workshop, 16–19 September 2018, Stony Point, New York, USA, 2018. a
Millan, R., Rignot, E., Bernier, V., Morlighem, M., and Dutrieux, P.: Bathymetry of the Amundsen Sea Embayment sector of West Antarctica from Operation IceBridge gravity and other data, Geophys. Res. Lett., 44, 1360–1368, https://doi.org/10.1002/2016GL072071, 2017. a
Morlighem, M., Rignot, E., Binder, T., Blankenship, D., Drews, R., Eagles, G., Eisen, O., Ferraccioli, F., Forsberg, R., Fretwell, P., Goel, V., Greenbaum, J. S., Gudmundsson, H., Guo, J., Helm, V., Hofstede, C., Howat, I., Humbert, A., Jokat, W., Karlsson, N. B., Lee, W. S., Matsuoka, K., Millan, R., Mouginot, J., Paden, J., Pattyn, F., Roberts, J., Rosier, S., Ruppel, A., Seroussi, H., Smith, E. C., Steinhage, D., Sun, B., van den Broeke, M. R., van Ommen, T. D., van Wessem, M., and Young, D. A.: Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet, Nat. Geosci., 13, 132–137, https://doi.org/10.1038/s41561-019-0510-8, 2020. a, b, c, d, e, f, g, h, i
Muto, A., Alley, R. B., Parizek, B. R., and Anandakrishnan, S.: Bed-type variability and till (dis)continuity beneath Thwaites Glacier, West Antarctica, Ann. Glaciol., 60, 82–90, https://doi.org/10.1017/aog.2019.32, 2019a. a, b, c
Muto, A., Anandakrishnan, S., Alley, R. B., Horgan, H. J., Parizek, B. R., Koellner, S., Christianson, K., and Holschuh, N.: Relating bed character and subglacial morphology using seismic data from Thwaites Glacier, West Antarctica, Earth Planet. Sc. Lett., 507, 199–206, https://doi.org/10.1016/j.epsl.2018.12.008, 2019b. a, b
Nias, I. J., Cornford, S. L., and Payne, A. J.: New Mass-Conserving Bedrock Topography for Pine Island Glacier Impacts Simulated Decadal Rates of Mass Loss, Geophys. Res. Lett., 45, 3173–3181, https://doi.org/10.1002/2017GL076493, 2018. a
Parizek, B. R., Christianson, K., Anandakrishnan, S., Alley, R. B., Walker, R. T., Edwards, R. A., Wolfe, D. S., Bertini, G. T., Rinehart, S. K., Bindschadler, R. A., and Nowicki, S. M.: Dynamic (in)stability of Thwaites Glacier, West Antarctica, J. Geophys. Res.-Earth, 118, 638–655, https://doi.org/10.1002/jgrf.20044, 2013. a
Pralong, M. and Gudmundsson, G. H.: Bayesian estimation of basal conditions on Rutford Ice Stream, West Antarctica, from surface data, J. Glaciol., 57, 315–324, https://doi.org/10.3189/002214311796406004, 2011. a
Raymond, M. J. and Gudmundsson, G. H.: Estimating basal properties of ice streams from surface measurements: a non-linear Bayesian inverse approach applied to synthetic data, The Cryosphere, 3, 265–278, https://doi.org/10.5194/tc-3-265-2009, 2009. a
Rignot, E., Mouginot, J., Morlighem, M., Seroussi, H., and Scheuchl, B.: Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, and Kohler glaciers, West Antarctica, from 1992 to 2011, Geophys. Res. Lett., 41, 3502–3509, https://doi.org/10.1002/2014GL060140, 2014. a
Scambos, T., Bell, R., Alley, R., Anandakrishnan, S., Bromwich, D., Brunt, K., Christianson, K., Creyts, T., Das, S., DeConto, R., Dutrieux, P., Fricker, H., Holland, D., MacGregor, J., Medley, B., Nicolas, J., Pollard, D., Siegfried, M., Smith, A., Steig, E., Trusel, L., Vaughan, D., and Yager, P.: How much, how fast?: A science review and outlook for research on the instability of Antarctica's Thwaites Glacier in the 21st century, Global Planet. Change, 153, 16–34, https://doi.org/10.1016/J.GLOPLACHA.2017.04.008, 2017. a
Smith, A. M., Murray, T., Nicholls, K. W., Makinson, K., Adalgeirsdóttir, G., Behar, A. E., and Vaughan, D. G.: Rapid erosion, drumlin formation, and changing hydrology beneath an Antarctic ice stream, Geology, 35, 127–130, https://doi.org/10.1130/G23036A.1, 2007. a
Smith, B. E., Gourmelen, N., Huth, A., and Joughin, I.: Connected subglacial lake drainage beneath Thwaites Glacier, West Antarctica, The Cryosphere, 11, 451–467, https://doi.org/10.5194/tc-11-451-2017, 2017. a, b
Thorsteinsson, T., Raymond, C., Gudmundsson, G., Bindschadler, R., Vornberger, P., and Joughin, I.: Bed topography and lubrication inferred from surface measurements on fast-flowing ice streams, J. Glaciol., 49, 481–490, https://doi.org/10.3189/172756503781830502, 2003. a, b, c, d, e, f
Turner, J., Orr, A., Gudmundsson, G. H., Jenkins, A., Bingham, R. G., Hillenbrand, C.-D., and Bracegirdle, T. J.: Atmosphere-ocean-ice interactions in the Amundsen Sea Embayment, West Antarctica, Revi. Geophys., 55, 235–276, https://doi.org/10.1002/2016RG000532, 2017. a
Vaughan, D. G., Corr, H. F., Ferraccioli, F., Frearson, N., O'Hare, A., Mach, D., Holt, J. W., Blankenship, D. D., Morse, D. L., and Young, D. A.: New boundary conditions for the West Antarctic ice sheet: Subglacial topography beneath Pine Island Glacier, Geophys. Res. Lett., 33, L09501, https://doi.org/10.1029/2005GL025588, 2006. a
Yu, H., Rignot, E., Seroussi, H., and Morlighem, M.: Retreat of Thwaites Glacier, West Antarctica, over the next 100 years using various ice flow models, ice shelf melt scenarios and basal friction laws, The Cryosphere, 12, 3861–3876, https://doi.org/10.5194/tc-12-3861-2018, 2018. a