Inverting ice surface elevation and velocity for bed topography and slipperiness beneath Thwaites Glacier

. There is signiﬁcant uncertainty over how ice sheets and glaciers will respond to rising global temperatures. Limited knowledge of the topography and rheology of ice-bed interface is a key cause of this uncertainty, as models show that small changes in the bed can have a large inﬂuence 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 ﬁll in the gaps. Here 5 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 ﬂow 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 10 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. This method could therefore be useful as either an independent test of other interpolation methods such as mass conservation and kriging, or as a complementary technique in regions where those techniques fail. 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 15 Thwaites Glacier. The methods presented here show signiﬁcant promise for using high-resolution satellite datasets, calibrated by the sparser ﬁeld 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. can be resolved is errors in the ice-surface data, which are simulated by adding noise to the synthetic surface generated from the synthetic bed. Noise is added with an amplitude of 2 m for the ice surface elevation (Howat et al., 2019), and an amplitude of 15 ms − 1 for the ice surface velocity (Gardner et al., 2018). Bedforms with an amplitude of less than 10 m, and slipperiness of less than around 1x10 − 4 m yr − 1 Pa − 1 are not well resolved.

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, b is the ice bed elevation, and α is the mean ice surface slope in the x direction.
We consider the response to a small perturbation in basal topography, b, linearising around a reference model (h,s,b,ū,v,c) with h =h + ∆h, s =s + ∆s, b =b + ∆b, u =ū + ∆u, v = ∆v, w = ∆w and c =c.
From depth integration of the Fourier-transformed incompressibility condition w z + u x + v y = 0 we have which, along with the steady-state boundary conditions, yields 85 ih(kû + lv) = −ikūŝ + ikūb.
Equations A7, A8 and A12 form a linear system of equations inŝ,û,v andb which can be solved algebraically (see Appendix A), leading to the steady-state transfer functions which represent the ratio of variability in the Fourier components of the surface to variability in the Fourier components of the bed. The following abbreviations are used for simplicity in the derivation: ξ = γ + 4hj 2 η, γ = τ 1−m d mc , j 2 = k 2 + l 2 , t r = j 2 τ dh cotα ξ , and ν = γ +hj 2 η.

Response of flow to basal slipperiness perturbation 95
Starting once again with the shallow-ice-stream equations (Equations A1 and A2;MacAyeal, 1989), this time we consider the response to a small perturbation in basal slipperiness, c, linearising with h =h + ∆s, s =s + ∆s, b =b, u =ū + ∆u, v = ∆v, w = ∆w and c =c(1 + ∆c) where ∆c is the fractional slipperiness.
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 T uc and T vc 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 115 reproduce the figures in that paper.

Non-dimensionalisation
The form of these transfer functions can be simplified by considering them in a non-dimensional form. For this purpose the same scalings as used in Gudmundsson (2003) and Gudmundsson (2008) are employed. All spatial scales are in units of mean ice thickness (h), and stress components are in units of driving stress (τ d ). Non-dimensional velocity components are in units 120 of mean deformational velocity (ū d ) wherē The scale for slipperiness is given byc/C. From Gudmundsson (2008) we know thatc/C =ū d /τ m d and we also haveC = u b /ū d . The ice surface velocity is the sum of the deformational velocity and the basal velocity, such thatū s =ū d +ū b . With some simple algebra, we can therefore express the scale for slipperiness in terms of the surface velocityū s , which is a known 125 quantity, The non-dimensional form of the equations is obtained using the substitutionsc →C, η → 1/2,h → 1,ū →C, γ → (mC) −1 , and τ d = ρghsinα → 1. The non-dimensional transfer functions can be found in the supplementary information, and are also shown graphically there. Non-dimensionalised parameters are represented by capital letters (B, C, S, U, V, T SB , T U B , T V B ,

130
T SC , T U C and T U V ).

The inverse problem
The non-dimensional transfer functions (T SB , T U B , T V B , T SC , T U C and T V C ) describe the relationship between the Fourier transforms of non-dimensionalised bed topography (B), bed slipperiness (Ĉ), surface topography (Ŝ) and surface velocity Since the system is over-determined, we can solve for non-dimensional bed topography and slipperiness using a weighted 140 least-squares inversion of equations 26, 27 and 28, with a filter applied to remove problematic wavelength components. Shortwavelength 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. Increasing the filtering parameter (p f ilt ) 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 Equation 26, 27, 28 are therefore 145 not repeated here in the main text, but are given in notation consistent with this paper in Appendix C.

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 a 120 m data spacing, purposely matching the data spacing of the ITSLIVE velocity product (Gardner et al., 2018). We subtract the mean bed elevation, slope 150 and slipperiness so that the bed varies about 0, and 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 Section 2.1.3, and Fourier transform to getB andĈ. The non-dimensional surface elevation (Ŝ) and the velocity components (Û ,V ) are calculated using the 155 forward model, and dimensionalised using the length scales given in Section 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 (B) and slipperiness (Ĉ) using the inversion procedure described in Section 2.2 and the supplementary information. After dimensionalisation, the inverted bed 160 can be compared to the synthetic bed to study the behaviour of the inversion.

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 (Equation C4) to p f ilt = −2. The mean ice thicknessh, mean surface slope α and mean ice surface velocity ū depend on the region studied, but for these synthetic tests are set ath = 2000 m, α = 0.002,ū = 100ms −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 ofC, the non-dimensional mean slipperiness, we set to beC = 100. When applied to Equation 25, these values give a dimensional mean slipperiness c = 2.7 x10 −3 myr −1 Pa −1 . We set the weighting factors (Equation C2) to be Σ s = 0.001, Σ u = 1 and Σ v = 1. This accounts 170 for the mismatch in the relative magnitude of the non-dimensionalised surface elevation and velocity, and means the leastsquares inversion solves for all three factors equally.

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 175 of the capabilities of the inversion. Sinusoidal basis functions vary depending on three parameter: 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 j 2 = l 2 + k 2 , λ/h = 2π/j, k = jcos(θ), and l = jsin(θ). 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 ITSLIVE (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 degrees, the wavelength is more than 2000 m and the amplitude is greater than 10 m for topography or 185 1.34 x10 −4 myr −1 Pa −1 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.

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 190 Glacier catchment as our example.
Our base data for surface elevation and velocity were, respectively, 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).
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 195 (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   Figure 3. The effect of amplitude on how well landforms and slipperiness (rows 1 and 3, respectively; created synthetically) can be resolved by the inversion. These tests are presented on a 50 km by 50 km grid, where in the inversion results (rows 2 and 4, respectively) the outer 5 km is greyed out to hide edge effects will should be neglected. In these simulations mean ice thicknessh = 2000 m, mean slipperiness C = 100, surface slope α = 0.02, wavelength λ = 20 km and angle to flow θ= 60 • . The limiting factor on how well small amplitudes can be resolved is errors in the ice-surface data, which are simulated by adding noise to the synthetic surface generated from the synthetic bed.
Noise is added with an amplitude of 2 m for the ice surface elevation (Howat et al., 2019), and an amplitude of 15 ms −1 for the ice surface velocity (Gardner et al., 2018). Bedforms with an amplitude of less than 10 m, and slipperiness of less than around 1x10 −4 m yr −1 Pa −1 are not well resolved.
could use a set of adjacent 50 km by 50 km grids, but instead we chose to use more densely distributed grids which overlap. This is because in the region of overlap, the variability in the inverted results is a measure of how well the physics that we 200 have assumed fits reality. Since the dominant approximation made in the physics is linearisation, we assume that variability in regions of overlap is mainly a measure of non-linearity in the physics of ice flow, although there could also be a contribution from the inappropriate use of the shallow-ice approximation. Variability is calculated as the standard deviation of overlapping solutions at each grid point. Discarding more of the outer part of each overlapping region further reduces edge effects. Figure   4a summarises this methodology. The bed topography and slipperiness results presented here are the grid-point by grid-point 205 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 ( Figure 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.
3 Results for Thwaites Glacier 210 Figure 5 shows the bed conditions we inverted from REMA (Howat et al., 2019) and ITSLIVE (Gardner et al., 2018) over a 280 km by 160 km region of the main trunk of Thwaites Glacier.
The bed topography product from the inversion is shown in Figure 5b. On the basin scale, the main basal topographic features identified by the inversion are several sets of parallel ridges which are oriented perpendicular 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 ( Figure   215 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, Figure 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, Figure   220 5b).
The standard deviation of the bed topography (Figure 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 225 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 Figure 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

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 swathradar-imaged bed at larger scales, picking out the locations of all the main hills and valleys. There is a better match for the 240 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 (Figure 7) and across ice flow (Figure 8) by airborne radar over Thwaites Glacier in the 2019/2020 field season (Jordan and Robinson, 2021). Further comparisons to more radar flight lines collected in the same field season can be seen in the supplementary 245 information. These figures also demonstrate that the inversion performs well in detecting the main hills and valleys, but also highlights that their amplitudes 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 slipperinessC = 100. If a non-dimensional slipperinessC = 150 is imposed (Figures 6, 7c, 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 250 the inversion ice thickness, rather than more detailed prior information.  Inverted bed for C = 150, r = 0.86, slope = 0.80 PASIN radar bed (data) PASIN radar bed (filtered) Bedmachine bed, r = 0.93, slope = 1.08 Inverted bed ±1 s.d. Figure 7. a) Comparative plot of inverted bed topography with mean non-dimensional slipperinessC = 100 with along-flow radar-sounded bed topography. Bed topography is given in three forms: as unfiltered bed picks from the 2019/20 airborne surveys (Jordan and Robinson, 2021); a version of the same filtered to 2 km wavelengths to be more representative of the detail we might expect to image in our inversion; and the bed profile extracted from BedMachine Antarctica (Morlighem et al., 2020). The envelope around the inverted bed topography shows plus or minus one standard deviation. The correlation coefficients (r) and slopes given are the results of a linear regression between the inverted bed or the BedMachine Antarctica bed and the filtered radar bed. b) Profile location within the inverted grid. c) As for panel a) but with mean non-dimensional slipperinessC = 150.  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 data sets. 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 255 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 260 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,265 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 ( Figure 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.

270
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 Davies et al., 2018) suggest that the bed is not changing rapidly there, and it is possible that neighbouring Thwaites 275 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 (MSGL) 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 does not only apply to the bed but also to the ice surface, which becomes more unstable closer to the grounding line. However, since results in the region immediately adjacent to the grounding line 280 are compromised by the different physics of the ice shelf anyway, this is not a significant concern. We therefore consider the steady-state assumption to be suitable for the purposes of this inversion.
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, p f ilt , the weighting factor Σ s and the mean slipperiness,C. The filtering parameter p f ilt (Equation There is less certainty over what is the most suitable value for the non-dimensional mean slipperinessC. 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). IfC is higher, then the amplitude of bed variability in the inversion 295 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,C, 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 additional we note that the trend is quite different from features observed in the inverted topography, showing that the slipperiness map is not a result 300 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 C = 100 gives the best match in these regions, whereas a higher value ofC = 150 (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 305 Figure 5 is from the inversion withC = 100.
This uncertainty inC means that some prior radar information is useful in order to calibrate the inversion method to give the best results. However, changingC 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.

310
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 315 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, as shown in Figure 9. However, even if a single ice-thickness (h) 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 Figure 9. a) The inverted bed topography, b) The reference bed topography (calculated from the ice surface with the 50-km-averaged ice surface subtracted) to which the inversion adds perturbations, and c) the Bedmachine Antarctica bed topography (Morlighem et al., 2020) for the region containing the Upper Thwaites radar grid (Holschuh et al., 2020). Panel d) shows in red the outline of this region within the inverted grid shown in Figure 5, and the locations of four subglacial lakes  and the two radar grids (Holschuh et al., 2020) . 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.
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 325 independent test is possible until more radar grids are collected.

Conclusions
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;330 Gudmundsson, 2003). Synthetic tests show that this method can resolve variability in bed topography and slipperiness on wavelengths greater than one ice thickness, and at amplitudes of more than 10 m for topography or 1 x10 −4 myr −1 Pa −1 for slipperiness, as long as the variability is not aligned to the ice flow direction. Comparison of the results of the inversion with radar grids and flight lines suggests that the inversion does a good job of picking up 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 335 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 340 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.
Code and data availability. The output data from the inversion is available on Zenodo at https://doi.org/10.5281/zenodo.5105687. The code for the inversion and plotting the figures is available on Zenodo at 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.
Appendix A: Derivation of transfer functions from a topography perturbation 350 Following Gudmundsson (2008), we start with the shallow-ice-stream equations (MacAyeal, 1989), where u, v, w are the velocity components in the x, y and z 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 355 is the ice surface elevation, b is the ice bed elevation, and α is the mean ice surface slope in the x direction.
This gives the first order equations 4ηh∂ 2 xx ∆u + 3ηh∂ 2 xy ∆v + ηh∂ 2 yy ∆u − γ∆u = ρghcosα∂ x ∆s − ρgsinα∆h (A3) 360 4ηh∂ 2 yy ∆v + 3ηh∂ 2 xy ∆u + ηh∂ 2 xx ∆v − γ∆v = ρghcosα∂ y ∆s (A4) To the first order, and importantly in the steady state, we also have the upper and lower boundary conditions All variables are then Fourier transformed with respect to the spatial variables x and y. In the forward Fourier transform 365 the two space variables both carry a positive sign and 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
which, along with the steady-state boundary conditions, yields Equations A7, A8 and A12 form a linear system of equations inŝ,û,v andb which can be solved algebraically using standard techniques: The determinant of the left-hand side of these equations: where the following abbreviations have been used to make the algebra clearer to follow: ν =hηj 2 + γ, j 2 = l 2 + k 2 , and which simplifies to: We then have: In the steady state, the kinematic boundary condition is Substituting the expression from above: In agreement with Gudmundsson (2008), this leads to the steady-state transfer function Expanding the expression forû: remembering that ξ = 4ηhj 2 + ν.
In agreement with Gudmundsson (2008), this leads to the steady-state transfer function: Expanding the expression forv: remembering that ξ = 4ηhj 2 + ν.
In agreement with Gudmundsson (2008), this leads to the steady-state transfer function: Appendix B: Derivation of transfer functions from a slipperiness perturbation Starting once again with the shallow-ice-stream equations (Equations A1 and A2;MacAyeal, 1989), this time we consider the response to a small perturbation in basal slipperiness, c, linearising with h =h + ∆s, s =s + ∆s, b =b, u =ū + ∆u, v = ∆v, w = ∆w and c =c(1 + ∆c) where ∆c is the fractional slipperiness.

440
This gives the first order equations 4ηh∂ 2 xx ∆u + 3ηh∂ 2 xy ∆v + ηh∂ 2 yy ∆u − γ∆u = ρghcosα∂ x ∆s − ρgsinα∆s − γū∆c (B1) 4ηh∂ 2 yy ∆v + 3ηh∂ 2 xy ∆u + ηh∂ 2 xx ∆v − γ∆v = ρghcosα∂ y ∆s (B2) 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 Equations B3, B4 and B5 form a linear system of equations which can be solved using standard algebraic techniques: The determinant of the left-hand side: where the following abbreviations have been used to make the algebra easier to follow: j 2 = l 2 + k 2 , ξ = 4hηj 2 + γ and Which simplifies to: We then have: At steady stateŝ t = 0, so we have the boundary condition: In agreement with Gudmundsson (2008) this leads to the steady-state transfer function: Expanding the expression forû (Eq B6): remembering that ξ = 4ηhj 2 + ν, leads to the steady-state transfer function: which is not as stated by Gudmundsson (2008).
remembering that ξ = 4ηhj 2 + ν, leads to the steady-state transfer function: which is also not as stated by Gudmundsson (2008).

505
Appendix C: The inverse problem The transfer functions (T sb , T ub , T vb , T sc , T uc and T vc ) describe the relationship between the Fourier transforms of bed topography (b), bed slipperiness (ĉ), surface topography (ŝ) and surface velocity (û,v). If the bed topography and slipperiness are known then surface topography and velocity components are given by the forward model: Non-dimensionalised this gives: S = T SBB + T SCĈ (Eq. 26 non-dimensionalised) U = T U BB + T U CĈ (Eq. 27 non-dimensionalised)

515V
= T V BB + T V CĈ (Eq. 28 non-dimensionalised) Since the system is over-determined, we can use a weighted least squares inversion of equations 26, 27 and 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.

520
In matrix form we have the forward model: 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 LM − KK * is small or zero, which is 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 wavelengths bed features and features aligned with ice flow are problematic because they cause flow disturbances in the ice which do not reach 540 the surface in a measurable way. They can therefore not be inverted from the surface data. To avoid this problem with an ill-conditioned inverse, Thorsteinsson et al. (2003) use a truncated version of L M − KK * as a filter, F to remove problematic wavelengths. Their filter is and p f ilt ≤ 0. This filter allows through all wavelengths where LM − KK * 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: