We present three types of data in the study: the spatial distribution of basal melt rate (see Figs. 2b and 3) , bathymetry inferred from airborne gravity surveys (shown in Fig. 2a), and locations of subglacial discharge (see Fig. 3) from the grounded ice sheet. We inferred the bathymetry of Getz using profile gravity inversions with the Geosoft GM-SYS software. The subglacial hydrological analysis was generated by the two-dimensional GlaDS model (Werder et al., 2013). The observed basal melt rates were computed using a mass conservation approach from Jenkins (1991) and Gourmelen et al. (2017b) and corrected for melting driven by warm ocean waters using ice bottom elevation and nearby ocean temperature profiles (Holland et al., 2008).

2.1 Helicopter gravity data acquisition The gravity data used in this paper were acquired aboard two aircraft types, one fixed-wing aircraft and one helicopter. Figure 1 shows the data coverage. The OIB data (Cochran and Bell, 2010a) were acquired using the Sander Geophysics Limited (SGL) AIRGrav system aboard NASA's DC-8. More details on this airborne geophysical platform can be found in the literature (Cochran and Bell, 2012; Cochran et al., 2014). The helicopter-based data were acquired using a Canadian Micro Gravity GT-1A in a collaboration between UTIG and KOPRI. Figure S1 shows the helicopter gravity data acquisition platform on the icebreaker Araon. Three dedicated aerogeophysical flights were accomplished in 1 d of helicopter operations from the Araon while off the coast of western Getz, acquiring about 1200 line-kilometers of data. The gravity anomaly (Fig. S2) suggests the effectiveness of combining OIB and helicopter data. The observed gravity anomaly ranges from −60 to 30 mGal (Fig. S2). The high anomaly strongly correlates with the ice rises and grounded icebergs. Large positive gravity anomalies of up to 30 mGal are consistently found over Grant Island, Dean Island, Siple Island, and Wright Island. The areas between ice rises correspond to low gravity anomalies. Both survey data sets show similar repeatability statistics with ∼ 1.4 to 1.6 mGal root mean square (RMS) in the differences at crossovers between lines both internally in each set and between sets. The ship-based UTIG–KOPRI gravity set did not have an absolute gravity tie so that entire survey set was level shifted to minimize the difference in the mean of crossovers with the OIB data; no other adjustments were made.

2.2 Bathymetry inversion approach The gravity data are inverted for depth of targets using GM-SYS Profile Modelling, a 2D gravity modeling and inversion module in Geosoft. In the forward-modeling mode, the module computes the gravity response from a polygon-approximated irregular target model (Talwani et al., 1959). In the inversion mode, the polygon-approximated model is adjusted iteratively to fit the observed gravity data best. Getz is pinned on an array of islands and peninsulas, so our bathymetry inversion is well constrained by the location of the ice rises and the peninsulas. The bathymetry model is updated iteratively until the difference between modeled gravity and observed gravity values is minimized (convergence limit = 0.1 mGal; 0.1 mGal is the standard error of observed gravity data). To better condition the inversion process, we fix the top and bottom of the ice layers, whose depths and topography are obtained from radar data. Similar approaches to infer bathymetry from airborne gravity data have been applied in many regions of Antarctica (Tinto and Bell, 2011; Cochran and Bell, 2012; Muto et al., 2016; Millan et al., 2017; Greenbaum et al., 2015). We first use the gravity data from grounded ice lines to invert for bedrock densities. For those areas covered by the grounded ice lines, we assume a three-layer model: a solid ice layer with density of 917 kg m−3 of known thickness over a bedrock layer, whose density is our free parameter; the third layer is the upper mantle with a density of 3300 kg m−3 at a depth of 20 km. The top, bottom, and thickness of the ice layer is obtained from OIB measurement and thus fixed throughout the inversion. We start the inversion with an initial guess of granitic rock density value 2.75 kg m−3 since west of the survey area has granite outcrop (Mukasa and Dalziel, 2000). The gravity data from floating ice lines are used to invert for the bathymetry under Getz. We use a four-layer model: the first layer is an ice layer with density of 0.917 kg m−3 with a known depth; the second layer is a seawater layer with a density of 1030 kg m−3 – the depth of this layer is our free parameter; the third layer is a bedrock layer with density inferred from grounded-ice-line gravity data; the fourth layer is the upper mantle with a density of 3300 kg m−3 at a depth of 20 km. The lines are processed one by one starting from those that are closer to grounded ice lines. The bathymetry model is updated iteratively until the difference between modeled gravity and observed gravity values is minimized (convergence limit = 0.1 mGal; 0.1 mGal is the standard error of observed gravity data). The polygon densities applied in this region are in Fig. S3. The constructed 2D models can be found in Fig. S4. The different 2D bathymetric profiles are merged through the minimum curvature gridding method, provided by the Grid and Image module from Geosoft Oasis montaj. The details of the minimum curvature method can be found in Briggs (1974). The mesh size of the interpolated grid is 1792 m. To prevent aliasing and high-frequency signals, we increase the low-pass desampling factor (i.e., the number of grid cells that are averaged). This factor (set to 3) removes high-frequency signals since it acts as a low-pass filter by averaging all points in the nearest cell. The distance between grid cells and a valid point greater than the blanking distance (set to 2000 m) are blanked out in the final grid. However, gridding could still introduce artifacts at the intersection points. We calculate the offsets at the profile intersections, and the average offset is about 20 m. The final derived bathymetry (Fig. 2a) includes the Getz Ice Shelf bathymetry from gravity inversion and offshore bathymetry from IBCSO (International Bathymetric Chart of the Southern Ocean; Arndt et al., 2013). We follow the uncertainty estimation approach from Greenbaum et al. (2015). We compare the inversion with the geometry of the grounded ice as a measure of the uncertainty beneath the floating ice assuming that the bed roughness under grounded ice and floating ice is similar. Our estimated root-mean-square error (RMSE) between the ice bottom measured by radar and sampled from the bathymetry model is about 246 m, and the mean offset between the two is about 44 m (see Fig. S5). We also compare the overlapping points where the gravity lines intersect with the ship track (Nitsche et al., 2007). The RMSE between the ship-measured bathymetry and that sampled from the bathymetry model is about 121 m; the mean offset between the two is about 32 m (see Fig. S5). We do not include significant geological or sedimentary signatures in our model since we have insufficient magnetic analysis over Getz. But our methods do account for local geological heterogeneity (see Fig. S3). Published interpretations within the ASE (Amundsen Sea Embayment) imply that we should not expect a significant crustal thickness gradient or sedimentary basin beneath the Getz Ice Shelf (Gohl et al., 2013). Therefore, we expect sediments near the grounding line to be scoured away as seen in other ice shelves of ASE (Gohl et al., 2013; Cochran et al., 2014). However, if the sediment is present, it will cause the gravity-derived bathymetry to be deeper than the actual seafloor. If we are correct and no significant geological structure underlies Getz, then the existence of sediments will shift the bathymetry to be deeper but will not change the shape of the bathymetry and thus will not affect our conclusions.

2.3 Subglacial hydrological model The subglacial hydrological analysis is generated by the two-dimensional GlaDS model (Werder et al., 2013). Distributed flow occurs through linked cavities that are represented as a continuous water sheet of variable thickness. Channels grow along finite element edges and exchange water with the adjacent distributed system, as part of a fully coupled 2D drainage network. The model is run at the steady state over 3000 d with primary outputs being channel discharge over the domain and the grounding line into the Getz cavity. Topography inputs are from airborne radar data; basal velocity is estimated as 90 % of MEaSUREs surface velocity data (Rignot et al., 2017); basal conductivity is assumed constant following other applications of GlaDS in Antarctica (Dow et al., 2016, 2018). The water input rate is set as constant (both spatially and temporally) at 10 mm yr−1 following subglacial melt rate calculations.

2.4 Observed basal melt rate The observed ice shelf basal melt rates are computed using a mass conservation approach from surface elevation, surface mass balance, ice velocity, and ice shelf thickness (Jenkins, 1991; Gourmelen et al., 2017b), using the relation (Jenkins, 1991; Gourmelen et al., 2017b) (1) - 1 - ρ ice ρ ocean m ˙ + SMB = ∂ S ∂ t + S ∇ ⋅ u , where ρ ice is ice density of 917 kg m−3, ρ ocean is the ocean density of 1028 kg m−3, m ˙ is basal melt rate, SMB is surface mass balance, S is surface elevation, and u is ice velocity. SMB is obtained from output of the regional atmospheric climate model RACMO2 (van Wessem et al., 2016). We derive the rates of surface elevation change from a new elevation data set, which is generated by the CryoSat-2 interferometric-swath radar altimetry from 2010 to 2016. Ice velocity is acquired from radar observation of the European Space Agency Sentinel-1A satellite. A detailed discussion of the methodology can be found in Gourmelen et al. (2017a). The observed melt rate of Getz is shown in Fig. 2b.