Geological sketch map and implications for ice flow of Thwaites Glacier, West Antarctica, from integrated aerogeophysical observations

The geology beneath Thwaites Glacier, the Antarctic glacial catchment most vulnerable to climate change, is unknown. Thwaites Glacier lies within the West Antarctic Rift System, but details of the subglacial geology relevant to glacial flow, including sediment availability, underlying lithology, and heat flux, are lacking. We present the first sketch map of the subglacial geology of Thwaites Glacier, interpreted from maps of airborne gravity, magnetic and radar data, supported by 2D models and 3D inversion of subsurface properties, and the regional geological context. A zone of Cretaceous mafic magmatism extending ~200 km inland from the coast is interpreted, while sedimentary basins are restricted to a region 150 to 200 km inboard of the coast, underlying just 20% of the catchment. Several granitic subglacial highlands are identified, forming long-lived topographic highs. Our geological interpretation places constraints on the basal properties of Thwaites Glacier, laying the foundation for both improved predictions of ice sheet change and studies of West Antarctic tectonics.


S2. Aeromagnetic compilation
The aeromagnetic compilation was created from line data made available through the ADMAP2 project (61) and incorporating more recent ITGC data (73,74). All data was merged into a single data-base and levelled using the ADMAP line data as the initial reference. Combination of all line data into a single database and re-gridding ensures anomalies at survey boundaries are accurately reproduced. The ITGC data underwent standard processing to account for the global reference field, errors introduced by operational procedures and diurnal variations, before it was merged and statistically levelled to the older datasets. Visual inspection of the resulting grid revealed some areas where additional levelling, or removal of inconsistent data was required. This most importantly included removal of all data collected in the 1970s by SPRI, which in this region is superseded by more modern and better geo-located data. Data were upward/downward continued to 500 m above the ice/sea surface to place all observations on a common surface (Fig. S3a). This continuation level was chosen as it preserves high frequency anomalies, while minimizing the extent that data is downward continued. An alternative continuation to 2500 m above the underlying rock was also carried out to give a view with a consistent distance to the closest possible source (Fig. S3b). The final continued line data were interpolated onto a 1 km mesh and reduced-to-the-pole (Fig. 2C), a procedure which theoretically aligns anomalies with their sources.

S3 Tilt depth magnetic depth to source method
For the tilt depth method the tilt angle was first calculated from the horizontal and vertical gradients of the magnetic field (77). The distance from every node on the zero tilt angle contour to the nearest +30° and -30° tilt angle contour was then extracted using ArcGIS 'near' function. The distance to the up and down contour were transformed to estimates of source depth and averaged (64,78). Two tests were then applied to isolate robust depth solutions, assuming the distance to the positive and negative tilt contour are independent estimates of source depth. First, the vectors to the closest positive and negative contours were within ±30° of being diametrically opposite. This prevents local kinks and areas with complex contours biasing the result. Second, points where depth estimates from the positive and negative contours differed by more than 50% of the average depth were excluded.
S4. Airborne gravity compilation and corrections.
Line gravity data from the surveys noted above was used to construct the free air gravity grid for this study. Detailed data editing and levelling was carried out prior to gridding to ensure a consistent and high quality product free of apparent artefacts. OIB data collected >1500 m above the ice surface was discarded, as the lack of resolution of short wavelength features on these transit flights distorts the resulting integrated gravity field. The line data from ITGC, AGASEA and BBAS were initially statistically levelled using the OIB dataset as a reference. Subsequent additional statistical levelling was applied to AGASEA, ITGC, and BBAS datasets to improve the gridded product. This included in places editing the AGASEA line data where erroneous anomalies became apparent on gridding. These regions typically corresponded to areas of aircraft elevation change where the L&R gravity system is known to perform poorly. Prior to gridding the line free air data was upward, or in some cases downward continued to a uniform level of 2500 m. This elevation was chosen as it is above the flight altitude of 95% of the dataset, minimizing the need for application of downward continuation, which is inherently less robust (Fig. S5a).
The full 3D gravity effect of the topography, bathymetry and ice was accounted for by calculating the Bouguer correction. For this the 2km mesh raster of topography was used as an input, including data up to 150 km from our study area to mitigate any lateral terrain effects. Ice surface elevations were taken from the BEDMAP2 data compilation (79). An observation altitude of 2500 m was used, as it is above the surface elevation across most of the study area. This is not true for the peak of Mt Takahea and some other of the volcanic peaks which lie outside our study region, but within the region of calculation. These discrepancies have limited impact on the calculated correction grid, but anomalies over Mt Takahea cannot be interpreted. Use of an observation level above all topography (~4000 m) would require additional continuation, degrading the observed gravity signal in the areas of most interest. We also note that the Bouguer correction is not robust over the floating ice shelves, as the bathymetry is initially derived from gravity data.
The Bouguer correction was calculated using a 3-D Gauss-Legendre quadrature (GLQ) method (80), which calculates the gravity effect of prisms of material with a specific density bounded by an upper and lower surface at each point on an observational mesh. The gravity effect of four layers was considered. Ice (915 kgm -3 ), water (1028 kgm -3 ), topography above sea level (2670 kgm -3 ) and topography below sea level (-2670 kgm -3 ). Summation of the gravity effect of these four layers provides the Bouguer correction. This correction was sampled onto the flight lines and subtracted from the free air anomaly. The resulting Bouguer anomaly was interpolated onto a 2 km mesh raster. The Bouguer gravity anomaly grid was filtered with 3 passes of a 9x9 Hanning filter to minimize residual line to line noise (Fig. S5b).
The final correction applied to the gravity data aims to isolate gravity anomalies associated with crustal geology from the gravity effect of isostatic compensation of the surface topographic, ice and water loads. Isostatic compensation is typically associated with a buoyant low-density crustal root beneath elevated topography. A simple way to calculate the magnitude of the compensating crustal root is to assume that the surface loads are compensated locally, i.e. by a buoyant root directly below. Such Airy compensation is consistent with the low elastic thickness calculated by previous studies in the adjacent Pine Island Glacier region (23). In the Thwaites Glacier region, especially towards Mary Byrd Land, the topography may be supported by the buoyancy of warm mantle, rather than a crustal root, i.e. so called Pratt type isostatic support (27).
However, we argue that an Airy model will adequately remove the first order gravity effect of isostatic support of the topography. The long wavelength gravity effect of the low density buoyant mantle required to support the topography will be approximately equivalent to the effect of a low density crustal root, as the cumulative mass deficiency required to support the topography is the same. The size and gravity effect of the compensating root was calculated using a fast Fourier transform routine within the GMT software package, assuming a uniform crustal density of 2670 kgm -3 and a mantle density of 3330 kgm -3 . This calculation requires assumption of a reference depth for the mantle where the surface topography has zero elevation. We chose a reference depth of 29 km, as this gives a mean predicted Moho depth matching the mean result from passive seismic estimates in this area (81). The Airy isostatic gravity model was subtracted from the Bouguer anomaly to give the final residual Airy isostatic anomaly (Fig. 2D).

S5. Detail of 3D inversion inputs and parameters.
The 3D inversion used the VOXI 3D inversion module of the Geosoft software suite (68). The inversion predicts a three-dimensional distribution of susceptibility or density from input magnetic or gravity observations. It does this using a Tikhonov minimum gradient regularisation (82) to produce a model conforming to a set of reference Earth properties, where the modelled response closely matches the input data. The inversion included a proprietary iterative reweighting inversion (IRI) focusing function, which acts to sharpen the boundaries of source bodies. By limiting the user-intervention to global parameters such as the minimum/maximum values for rock properties and depth weighting, use of this inversion produces an output for interpretation relatively unbiased by user intervention. However, we acknowledge that the inversion output, like other minimally constrained potential field models, is non-unique and other geometries could provide an equivalent fit to the data. The input magnetic field for the 3D inversion was the magnetic anomaly continued to 500 m above the ice surface, with a mean and best fit linear trend removed. Susceptibility values were constrained to be >0 and <75x10 -3 SI consistent with values for typical igneous rocks (38). The recovered source properties implicitly include the impact of any magnetic remanence and are therefore apparent susceptibility values. An active magnetic model volume extending from 0 km to -10 km depth, with a mesh resolution of 1.25/1.25/0.5 km was used, with an expansion ratio of 1.5 with depth. Five padding cells were included below the active model. As many of the depth to source solutions indicate shallow sources, and to counter the tendency of the inversion to assign higher susceptibility values to deeper levels in the model, we imposed an additional linear weighting scheme favoring higher susceptibility values at shallower depths.
When inverting for density the residual Airy isostatic gravity anomaly, with a mean and best fit linear trend removed was considered as the input data. Density was constrained to lie between +/-1000 kgm -3 , and recovered values reflect a density contrast, rather than absolute density values. An active model volume extending from ~2.5 km to -12 km depth, with a mesh resolution of 2.5/2.5/1.25 km was used, with an expansion ratio of 1.5 with depth. Five padding cells below the active model were included. Finer mesh resolution, or weighting the results to favor shallow results failed to give valid results within the inversion scheme. For the magnetic and gravity inversion cells within the ice column were set to either zero susceptibility or zero density contrast, ensuring recovered geophysical sources were below the ice-bed interface.  Table S1).     Additional 2D forward models and slices of 3D inversion volume over down-stream magnetic and gravity anomalies. Map shows Airy isostatic anomaly and locates profiles L1 to L5, and the key profile K-K' shown in the main text Fig. 4. For each profile upper panels show observed and forward calculated magnetic and gravity anomalies while lower panels show 2D forward models (black bodies) and inverted magnetic susceptibility or density structure.

Fig. S7.
Input and output data, together with difference maps for the 3D inversion. a) Input magnetic anomaly grid. b) Input Airy isostatic gravity residual grid. c) Predicted magnetic anomaly map. d) Predicted gravity anomaly map. e) Magnetic anomaly error after inversion (input -output). Overall standard deviation is 8.5 nT, but a negative of up to -120 nT is present indicating magnetic remanence may be an issue in this area. f) Gravity anomaly error after inversion (inputoutput). Overall standard deviation is 0.35 mGal, suggesting the model is over-fit compared to the quality of the input data (2.94 mGal standard deviation).