Multidecadal Basal Melt Rates and Structure of the Ross Ice Shelf, Antarctica, Using Airborne Ice Penetrating Radar

Basal melting of ice shelves is a major source of mass loss from the Antarctic Ice Sheet. In situ measurements of ice shelf basal melt rates are sparse, while the more extensive estimates from satellite altimetry require precise information about firn density and characteristics of near‐surface layers. We describe a novel method for estimating multidecadal basal melt rates using airborne ice penetrating radar data acquired during a 3‐year survey of the Ross Ice Shelf. These data revealed an ice column with distinct upper and lower units whose thicknesses change as ice flows from the grounding line toward the ice front. We interpret the lower unit as continental meteoric ice that has flowed across the grounding line and the upper unit as ice formed from snowfall onto the relatively flat ice shelf. We used the ice thickness difference and strain‐induced thickness change of the lower unit between the survey lines, combined with ice velocities, to derive basal melt rates averaged over one to six decades. Our results are similar to satellite laser altimetry estimates for the period 2003–2009, suggesting that the Ross Ice Shelf melt rates have been fairly stable for several decades. We identify five sites of elevated basal melt rates, in the range 0.5–2 m a−1, near the ice shelf front. These hot spots indicate pathways into the sub‐ice‐shelf ocean cavity for warm seawater, likely a combination of summer‐warmed Antarctic Surface Water and modified Circumpolar Deep Water, and are potential areas of ice shelf weakening if the ocean warms.


Introduction
Antarctica's ice shelves, the floating extensions of the Antarctic Ice Sheet, lose mass through basal melting and calving. Integrated around all of Antarctica, each of these processes contributes about half of the total ice shelf loss (Depoorter et al., 2013;Rignot et al., 2013Rignot et al., , 2019. Ice shelves interact with both the atmosphere and the ocean, making them particularly vulnerable to climate variability. Ice shelves play a vital role in controlling the loss of grounded ice from Antarctica through a process called "buttressing," in which they impart a back stress on the upstream grounded ice. Any change in ice shelf geometry that reduces back stress on the upstream grounded ice can lead to increased mass flux from the buttressed grounded-ice basin (Joughin et al., 2014;Minchew et al., 2018;Pritchard et al., 2012;Scambos et al., 2004).
In turn, the accelerated flux of grounded ice into the ocean will increase the ice sheet's contribution to global sea level rise.
Various methods have been used to measure ice shelf basal melt rates. In situ measurements have been obtained using ice penetrating radar (Catania et al., 2010), autonomous phase-sensitive radio-echo sounder radars (e.g., Marsh et al., 2016;Stewart et al., 2019), and moored turbulence sensors (Stanton et al., 2013) and fiber-optic cables (Kobs et al., 2014) deployed through boreholes in the ice shelf. These measurements are, however, spatially sparse and provide estimates only over short periods, generally less than 1 year. In contrast, surface elevation data from satellite altimeters provide ice-shelf-wide estimates of basal melt rates that are averaged over several years (Depoorter et al., 2013;Pritchard et al., 2012;Rignot et al., 2013;Shepherd et al., 2004) to more than two decades Paolo et al., 2015Paolo et al., , 2018. However, the accuracy of these melt rates is limited by uncertainties in models of surface snow and firn density (Davis & Moore, 1993;Ridley & Partington, 1988), incorrect assumption of hydrostatic equilibrium at small scales , and errors in ice advection and strain rate estimates from satellitederived ice velocity (Alley et al., 2018).
For rapidly melting, small ice shelves such as those in the Amundsen Sea (Paolo et al., 2015Pritchard et al., 2012), the area-averaged basal melt rates estimated from satellite data far exceed potential errors, and  (Rignot et al., 2017), with velocity streamlines (black lines): streamline 167 following Byrd Glacier is used in Figure 3. RI: Roosevelt Island, CIR: Crary Ice Rise, J9: ice core that showed marine ice (Zotikov et al., 1980). (c) REMA DEM (Howat et al., 2019) showing surface morphology, streaklines (dark lines), and selected velocity streamlines (yellow). TI and T2 are the locations where the streaklines differ most from the streamlines due to transient events. (d) Modeled change in ice thickness over 1,000 years due to transient events (Hulbe & Fahnestock, 2007). the uncertainty in freshwater production integrated over each ice shelf is small. For example, a 0.05 m a −1 error in correcting for firn-air content for the Pine Island Glacier ice shelf (area~6,000 km 2 ) would result in an integrated mass loss error of~0.3 Gt a −1 after the elevation signal has been corrected to mass by the hydrostatic assumption. This error is negligible compared with the ice shelf's estimated loss of about 100 Gt a −1 by basal melting (Rignot et al., 2013). However, the same error applied over the Ross Ice Shelf (~480,000 km 2 , Figure 1) results in an integrated error of about 25 Gt a −1 . This value is about half of the expected mass loss of~50 Gt a −1 from the Ross Ice Shelf by basal melting, estimated as the residual of the mass balance equation (Rignot et al., 2013) and by integrating satellite-based estimates of basal melt rates (Moholdt et al., 2014). This example demonstrates the value of developing alternative methods for largescale mapping of basal melt rates that are less sensitive to uncertainties in surface firn density.
In this paper, we describe a new method for estimating basal melt rates by imaging changes in an ice shelf's internal structure along ice flowlines with airborne ice penetrating radar. We apply this method to obtain multidecadal melt rate estimates for the Ross Ice Shelf, using ice-penetrating radar data from the ROSETTA-Ice airborne survey (Tinto et al., 2019). Our method takes advantage of the ability of the radar to resolve spatially coherent internal reflectors below the firn layer, so that we can evaluate the contribution of basal mass balance to ice shelf thinning along flowlines independently of ice-shelf surface characteristics. We interpret our derived distribution of melt rates in terms of the oceanographic processes driving mass loss from this ice shelf.

Horizontal Structure
The Ross Ice Shelf (Figure 1), the largest ice shelf in Antarctica, receives ice from the grounded catchments of both the West Antarctic and East Antarctic ice sheets (WAIS and EAIS, respectively). These catchments represent a total of about 12 m of potential global sea level rise (Tinto et al., 2019). The ice shelf is fed by several broad ice streams draining WAIS, and by EAIS outlet glaciers flowing through narrow valleys in the Transantarctic Mountains. Based on modern ice velocity maps derived from satellite remote sensing (Figure 1b), the typical time taken for an ice parcel to flow from the grounding line to the calving front is about 1,000 years (Tinto et al., 2019, their Figure ED3c).
On the ice shelf, there is distinct topographic structure normal to ice flow, forming along-flow "streaklines" that are visible in satellite products such as the Moderate Resolution Imaging Spectroradiometer Mosaic of Antarctica (Hulbe & Fahnestock, 2007), the Landsat Image Mosaic of Antarctica (Bindschadler et al., 2008;Fahnestock et al., 2016), and the Reference Elevation Model of Antarctica (REMA: Howat et al., 2019). These streaklines delineate paths of ice parcels from the grounding line to the ice front ( Figure 1c). Differences between streaklines and modern streamlines calculated from satellite-based measurements of ice velocity indicate changes in ice flow across the grounding line associated with decadal to millennial variability of grounded ice streams and glaciers (e.g., Campbell et al., 2018;Hulbe & Fahnestock, 2007).
On the WAIS side of the Ross Ice Shelf, Kamb Ice Stream and its tributaries have stagnated (Conway, 2002;Retzlaff & Bentley, 1993), Whillans Ice Stream is slowing down (Beem et al., 2014), and MacAyeal and Bindschadler ice streams have changed speeds (Hulbe & Fahnestock, 2007). The response of the ice shelf to these changing ice inputs is reflected in the ice thickness: ice is thinner downstream of sites where the incoming ice has slowed, and thicker where incoming ice has sped up (Campbell et al., 2018;Hulbe & Fahnestock, 2007;Pritchard et al., 2012) (Figure 1d).
On the EAIS side of the Ross Ice Shelf, streaklines and modern streamlines diverge near the outlet of Nimrod Glacier ( Figure 1c); however, this feature may represent the effect of changing flow of WAIS ice on the EAIS portion of the ice shelf in the past rather than changing inputs from the grounded portion of EAIS (Reese et al., 2018). Short-term changes in mass flux have been observed for Byrd Glacier, which delivers about 22 Gt a −1 of ice to the Ross Ice Shelf from EAIS (Stearns, 2011); this glacier experienced a 14-month increase in surface velocity and mass flux in 2005-2007 that was associated with a subglacial lake discharge event (Stearns et al., 2008). However, we do not expect such short-term fluctuations to lead to substantial offsets between streaklines and modern streamlines.

Vertical Structure
An ice shelf's vertical structure can consist of up to three distinct ice units based on where and how they formed, and each unit may consist of multiple layers resolved in a radar profile (Craven et al., 2009). The continental meteoric ice (CMI) unit is the package of ice that has flowed from the grounded ice sheet to form the ice shelf. The stratigraphy of this unit may be substantially deformed as it encounters variable stress regimes while still grounded (Wright et al., 1990). As CMI advects downstream on the ice shelf, it is progressively buried by local snow accumulation, which forms an upper unit of local meteoric ice (LMI). The LMI unit generally has horizontal layers because it forms on the almost flat CMI unit, which experiences little basal drag after going afloat. In some regions, marine ice accretes to the bottom of the ice shelf when basal conditions favor direct freezing or formation of buoyant frazil ice within the water column (Craven et al., 2009;Fricker et al., 2001;Joughin & Vaughan, 2004;Lewis & Perkin, 1983;Zotikov et al., 1980).

Mass Balance
The Ross Ice Shelf gains mass from ice flow across the grounding line and snowfall on the ice shelf surface and loses mass by iceberg calving and basal melting; overall, the ice shelf is approximately in balance (Depoorter et al., 2013;Moholdt et al., 2014;Rignot et al., 2013). We summarize each term of mass balance below; values in parentheses are ice-shelf-integrated mass flux per year from Rignot et al. (2013), their Table 1.

Ice Flux Across the Grounding Line (+129 ± 8 Gt a −1 )
This process accounts for about two thirds of the total mass input for the Ross Ice Shelf. About 73.0 ± 1 Gt a −1 of this ice flux is via WAIS ice streams, while 56.1 ± 4 Gt a −1 is via EAIS glaciers (Rignot et al., 2013).

Precipitation Onto the Ice Shelf (+65 ± 12 Gt a −1 )
Averaged over the whole ice shelf, the surface accumulation is about 0.14 ± 0.025 m of water equivalent per annum (m w.e. a −1 ). Models suggest that rates vary across the ice shelf by about a factor of two (Moholdt et al., 2014, their Figure 7a). A 2,700-year record of annual accumulation at the Roosevelt Island Climate Evolution ice core (Winstrup et al., 2019) revealed a range of annual precipitation rates from about 0.12 to 0.40 m w.e. a −1 between 1900 and 2010, but with the decadal running-mean values remaining within the range 0.20-0.25 m w.e. a −1 .

Iceberg Calving (−146 ± 12 Gt a −1 )
This process accounts for about 70-75% of total mass loss and is dominated by intermittent production of large tabular icebergs, with length scales of tens to hundreds of kilometers (e.g., Lazzara et al., 1999). Note that Rignot et al. (2013) evaluated mass loss by calving from the ice flux through a fixed "flux gate" near the ice front; basal melting for any region of ice shelf downstream of this flux gate was counted as calving rather than melting.

Basal melting (−48 ± 34 Gt a −1 )
This process accounts for about 25-30% of the total mass loss from the Ross Ice Shelf south of the flux gate typically used for estimating calving flux (e.g., Rignot et al., 2013). The area average of the basal melt rate (w b ) is about 0.10 ± 0.07 m w.e. a −1 , similar to the surface mass balance. However, w b varies spatially over a wide range, from close to 0 over much of the central portion of the ice shelf to over 10 m a −1 near the deep grounding line of Byrd Glacier (Kenneally & Hughes, 2004) and in a channel near Whillans Ice Plain (Marsh et al., 2016). Relatively high basal melt rates of about 1-2 m a −1 have been estimated near the ice shelf front from Lagrangian analyses of ICESat repeat-track laser altimetry (Horgan et al., 2011;Moholdt et al., 2014) and from an autonomous phase-sensitive radio-echo sounder survey in a small region near Ross Island (Stewart et al., 2019). Satellite-based estimates of melt rate identify extensive regions where ice appears to be accreting to the base of the ice shelf (i.e., w b < 0); however, the rates are usually smaller than the associated uncertainties (Moholdt et al., 2014). A 6 m layer of marine ice was observed under the Ross Ice Shelf at the J9 drill site north of Crary Ice Rise (Zotikov et al., 1980; see Figure 1b for location). In addition, thickening of up to 0.3 ± 0.15 m a −1 is observed in some flowbands south of Crary Ice Rise (Thomas & Bentley, 1978).

Overview of Method: Deriving Basal Melt Rate From Airborne Radar Data and Satellite Ice Velocities
Our primary goal is to estimate the spatial distribution of basal melt rate (w b : positive for melting) from measurements of changing CMI thickness (H CMI ) along ice flowlines. The local (Eulerian) rate of thickness change of the CMI layer is where ∇ · V is the divergence of horizontal ice velocity V(x). We assume that the three dimensional ice divergence is zero so that ∇ · V is equal to the negative of the vertical strain rate. Values of ∇ · V greater than 0 represent divergent flow (ice shelf thinning); values less than zero identify convergent flow (thickening). No surface processes such as accumulation or firn densification are included in equation (1) since most of our estimates of depth of the CMI unit are greater than the expected depth of the firn layer (~50-70 m; Cuffey & Paterson, 2010;Ligtenberg, 2014). That is, we assume that the CMI unit is always composed of fully compacted solid ice.
Assuming steady state (H CMI (x) and V(x) are constant with time), and converting to Lagrangian coordinates (i.e., following a column of ice), the thickness change (ΔH CMI ) between one survey line and a second line where t 1 and t 2 are the times at which the ice column passes through the first and second lines of each pair of the ROSETTA-Ice survey lines ( Figure 1a).
The full expected change in H CMI for a specific column of ice can be evaluated over an arbitrary distance along-flow, allowing for time-dependent variability in V and w b , by integrating equation (2) between two arbitrary points on the same flowline. Here, however, we restrict our estimates to time differences of a few decades or less, for which we expect the time-dependent variability of most of the Ross Ice Shelf to be small  compared with the signals from steady-state ice divergence and basal melting.
With this assumption, we approximate equation (2) as where Δt = t 2 − t 1 is the time taken for the ice column to travel from the first to the second survey lines, and the overbars denote averages over this interval of the flowline. We reformulate equation (2) to estimate the time-averaged basal melt rate as We refer to the first term on the right-hand side as the "strain-induced thickness change rate." Estimating basal melt rates from equation (4) requires measurements of the thickness of the CMI unit, velocity of ice flow along flowlines (to provide Δt = Δx/|V|, where V is the mean velocity between the survey linesÞ; and horizontal divergence of ice flow that contributes to the net Lagrangian thinning rate.

Data Sources
Our primary data sources are ice penetrating radar measurements from an airborne survey of the Ross Ice Shelf and satellite observations of ice velocity.

Total Ice Thickness and Ice Shelf Structure
We used airborne ice penetrating radar data acquired during the 3 years (2015, 2016, and 2017) of the ROSETTA-Ice survey of the Ross Ice Shelf (Tinto et al., 2019). Our survey was flown with 10-30 km spacing between survey lines (Figure 1a), for a total flight distance of about 61,000 km. Two radar systems were deployed in the IcePod instrument package that was mounted on a New York Air National Guard LC-130 aircraft. The shallow ice radar (SIR) operates at 2 GHz center frequency and 600 MHz bandwidth to image internal structures closer to the ice surface and up to a depth of about 500 m (Figures 2a and 2b). The deep ice radar (DICE) has a 188 MHz center frequency and 60 MHz bandwidth designed to identify the base of the ice shelf ( Figure 2c). The range resolution of the DICE in ice is 1.4 m and that of SIR is 0.14 m.
For both sets of radar data, we converted the two-way travel time to depth using the velocity of electromagnetic waves in solid ice, 1.68 × 10 8 m s −1 . We did not apply any correction to account for increased velocities in firn because firn depths and density profiles are not precisely known throughout the ice shelf, and because our primary interest is in thickness changes of the solid ice of the CMI unit for which the radar velocity is well known.
We used the DICE data to identify the base of the ice shelf and to derive total ice shelf thickness, H T ( Figure 2c). The base of the ice shelf is usually identified by a bright reflector (Figure 2c), which was picked by a semi-automatic picker also used in earlier studies Das et al., 2013). The SIR data identified distinct internal layers indicating the stratigraphy of the ice column (Figures 2a and 2b). The dominant feature in the SIR profiles was a bright reflector separating two distinct units. The depth of these units vary for different outlet glaciers and ice streams ( Figure 2b). The first documented imaging of this layer may have been close to the grounding line of Whillans Ice Stream by Wright et al. (1990).
We interpret the upper unit, which exhibits generally horizontal stratification, as the LMI unit that forms as snowfall onto the ice shelf accumulates and is compacted (Figures 2a and 2b). The thickness of the LMI unit (H LMI ) is usually greater than the expected firn depth in that location. The boundary between the LMI and CMI units was usually easy to identify, with the exception of highly crevassed regions near shear margins (Figures 2a and 2b). The CMI-LMI interface on the Ross Ice Shelf was also resolved in IceBridge radar, using a depth sounder that has roughly the same frequency (Paden et al., 2014) as the ROSETTA-Ice DICE ( Figure 2d, and see Figure 1a for location). This demonstrates the consistency in identifying the same interface by two independent surveys. We identified and traced the CMI-LMI transition using SIR radargrams using the same semiautomated picker used for identifying the base of the ice shelf in DICE radar. We calculated CMI thickness (H CMI ) by subtracting H LMI measured with SIR from H T measured with DICE.

Ice velocity and Strain Rates
We used ice velocity maps from two sources: satellite measurements, and a regional ice-sheet model tuned to observed velocities and ice thicknesses. The observed ice velocity field (Figure 1b) was obtained from the Making Earth System data records for Use in Research Environments (MEaSUREs) project, Version 2 Rignot et al., 2011Rignot et al., , 2017, which was developed from data acquired during 2011-2016 from RADARSAT-2, Sentinel-1, TerraSAR-X, and TanDEM-X synthetic aperture radar, and Landsat-8 optical imagery. This data set is provided at a 450 m grid spacing. To reduce noise in gradient-based strain rate estimates, we smoothed the velocity field using a 10 km × 10 km kernel and inverse square weighting following Moholdt et al. (2014). We assumed that the velocity through the ice column is uniform since there is no stress at the base of the ice column to support vertical shear (Weertman, 1957).
We also used the vertically integrated shallow-shelf component of the open-source ELMER/Ice model (Gagliardini et al., 2013) to generate a steady-state model of the ice shelf using methods based on Fürst et al. (2016). We inverted satellite ice velocity and thickness for the flow-law rate factor (the parameter governing the relationship between stresses and resulting strain rates) using a control method that minimized the difference between observed and simulated velocity (e.g., Fürst et al., 2016). The model physics results in the velocity field being smoother than the MEaSUReS 450 m product and removes some of the artifacts in that product. We used these two strain rate fields to assess the sensitivity of basal melt rates to our smoothing techniques and strain-induced thickness change using equation (4).

Journal of Geophysical Research: Earth Surface
We derived an independent estimate of vertical strain rates using ground-based measurements collected in 1973-1975 during the RIGGS survey (Bentley, 1990;Thomas et al., 1984). We used the principal components of strain rates provided by Thomas et al. (1984) to compute the divergence of the horizontal velocities. The RIGGS strain rate measurements were made by optical and electronic survey of the orientation and distances of strain rosettes in a region of about 1-3 km about the nominal sites on a grid of roughly of 55 km. The errors in the principal strain rates are about ±5 × 10 −5 a −1 (Thomas et al., 1984).

Implementation of Method
We used the method described in section 3.1 and applied it to the data described in section 3.2 to obtain basal melt rates. Figure 3 shows transects, along a single flowline (streamline 167 on Byrd Glacier, see Figure 1b for location), of the data and derived quantities following equations (1)-(4). The smoothing windows and length scales for extracting parameters along velocity streamlines for deriving basal melt rates are provided in Table 1. Values of H LMI and H T (from SIR and DICE, respectively) were evaluated as averages of the individual thickness estimates within 200 m of the streamline's intersection with the survey lines (i.e., 400 m along track with the center on the velocity streamline; see Table 1). This window size or length scale was chosen so that the flowlines sampled the same glacier packets from one survey line to the next to derive thickness changes.
Strain rate values were averaged within a sampling window of 20 km × 3 km (Table 1). This sampling window size was chosen to accommodate the velocity grid's smoothing of 10 km × 10 km, to get a reasonable spread of values around the survey lines when they are spaced 10 km apart (Figure 1a), and to suppress large values of strain rates when a survey line intersects shear margins and crevasses. Strain rates do not change very rapidly around any given point, unless they are located over shear margins, rifts and crevasses.
Each time a flowline intersects a survey line, the radar thicknesses were averaged over~150 values and the strain rates were averaged over~1,200 values. The process for calculating the basal melt rate is summarized in the flowchart in Figure 4a. The combination of the satellite-derived ice velocity and the spacing of ROSETTA-Ice survey lines (Figure 1a) results in a time scale between along-flow measurement points, Δt = Δx/|V|, that spans multiple decades (Figure 4b). The time scale, Δt, is~10-20 years near the ice front, where Δx is 10 km and |V| is about 1 km a −1 , and longer (up to~60 years) in the central and southern ice shelf where Δx is typically 20-30 km and |V| is about 0.5 km a −1 .

Error Sources
We identified two principal sources of error in our calculation of w b using equation (4) and the data described above: uncertainties in ΔH CMI arising from picking the horizon of the ice-penetrating radar data (section 3.2.1), and errors from temporal variability of the ice velocity and strain fields on the 10-60 year time scales of ice flow from one survey line to the next (sections 3.2.2 and 3.3).
Misidentifying the base of the ice shelf or the CMI-LMI interface results in errors in H T and H LMI , respectively . We used the standard deviation of the two ice thicknesses around each intersection of the flowlines with survey lines as a measure of the variability of ice thicknesses including the picking variability around that point. We tested the internal consistency of thickness estimates from both SIR (H LMI ) and DICE (H T ) by comparing measurements at crossovers between the east-west and the north-south survey lines that were spaced about 55 km apart to fly over RIGGS sites (see Figure 1a). Crossover estimates during the 3 years of data collection were consistent to 2 m for H T , and 10 m for H LMI . Other sources of errors include the velocity uncertainties and the errors in calculating strain-induced thickness change rates. We estimated total uncertainties in basal melt rates, ϵ total , from the combination of error sources as In equation (5), ∈ HCMI is the standard deviation of H CMI around each streamline intersection with a survey line (Table 1) and H CMI is the mean CMI thickness in the same averaging space. ∈ v is the local surface velocity error, assumed to be 30 m a −1 following Moholdt et al. (2014), and V is the mean velocity within a 20 km × 3 km along-track segment (choice of length scale explained in section 3.3). ∈ strain is the standard deviation of satellite-derived and ELMER/Ice-derived strain rates within a 20 km × 3 km segment (choice explained in section 3.3) of a streamline point. The largest source of error in equation (5) is the error in strain-induced thickness change (the last term in equation (5)). We estimated an uncertainty for strain-induced thickness change by dividing ∈ strain with the mean strain rates around the same sampling window (Table 1), and multiplied by the mean w b within a 400 m length scale along a streamline point (third term in equation (5)).
Equation (4) evaluates melt rates using the change of thickness of the CMI unit, ΔH CMI , between two survey lines; therefore, melt rates represent an average value for the time it takes for ice to flow between the two lines (10-60 years; see Figure 4b). Thomas et al. (2013) documented changes in ice-shelf velocity over the four-decade period from the time of RIGGS in situ measurements (1973)(1974)(1975) to the modern epoch of velocity measurements by satellites. These changes were largest near the grounding lines of the southern WAIS ice streams (Mercer and Whillans), where speeds decreased by up to 40% between 1975 and 2009 (Thomas  Figure 2). Over most of the ice shelf, however, the differences between RIGGS and modern satellite-observed velocities were less than 15%. For most cases, we excluded areas with large observed velocity changes (>15%) over the last four decades, since uncertainties associated with both the ice advection time (Δt) between survey lines and the strain-induced thickness change rate of the lower unit, H CMI ∇·V, would be large.
The recent changes in ice shelf velocities are associated with thickness transients due to past variation of ice stream and glacier flow into the ice shelf (Campbell et al., 2018;Catania et al., 2012;Hulbe & Fahnestock, 2007). The area around Crary Ice Rise (near T1 and T2 marked in Figures 1c and 1d) is of particular interest. This area was identified, and the magnitude of the spatially variable, transient dH T /dt was quantified, using the numerical model of ice shelf flow in Hulbe and Fahnestock (2007). That study used a diagnostic/prognostic model of ice shelf flow to simulate past ice discharge events implied by observed distortions of streaklines on the ice shelf. Matching streakline patterns required stagnation and reactivation of tributary ice streams and this, in turn, generated the dynamic thickness transient (unrelated to basal melting) that we recognize here. For the remainder of the ice shelf, we assumed that errors in w b caused by errors in Δt were small compared with uncertainties arising from radar thickness errors and strain, so that we could use the values of Δt from modern velocity fields.
We investigated the variability in strain rates over the past 40 years by comparing our strain rates from modern, satellite-derived velocity fields with the RIGGS strain rates. For some areas downstream of Whillans and Kamb and within our domain, models show effects of the transient thickness changes driven by the stagnation of Kamb Ice Stream~180 years ago, and ongoing flow changes on Whillans Ice Stream (Figure 1d), with additional signals from other, older events (Campbell et al., 2018;Hulbe et al., 2016;Hulbe & Fahnestock, 2007). These are the locations where the velocity fields are not in steady state, and the modern velocity streamlines may not represent the true pathways for thickness change.
To assess the effect of nonsteady-state conditions near T1 and T2 on dH CMI /dt we differenced dH CMI /dt derived using satellite-based strain-induced thickness change and the nonsteady-state model from Hulbe and Fahnestock (2007) described earlier. Figure 1c shows that, near Crary Ice Rise (near T2), the velocity streamlines are significantly different from the streaklines. The difference in dH CMI /dt from the two methods shows that, near T1 and T2, a steady-state assumption while a nonsteady state prevails will likely result in errors in calculating basal melt rates.

Ice Shelf Structure: Total Ice Thickness, and Thicknesses of Local and CMI Ice Units 4.1.1. Total Ice Thickness
The DICE radar resolved total ice thickness, H T (Figure 5a), from the grounding line to the ice shelf front and shows large variability in ice thickness. In general, the ice shelf is thicker on the WAIS side than on the EAIS side.
Measurements of H T were obtained during the RIGGS campaign (1973)(1974)(1975)) also using airborne radar (Robertson & Bentley, 2013;Thomas et al., 1984). Comparison of H T from RIGGS and ROSETTA-Ice radar (Figure 5a) at the locations of the RIGGS sites showed good agreement between the two data sets, although RIGGS measurements included corrections for~7 m of firn air on the Ross Ice Shelf. The mean difference was 3.7 m (ROSETTA-Ice thinner than RIGGS) with a standard deviation of 4.8 m calculated at the RIGGS locations. The low mean ice thickness difference between ROSETTA-Ice and RIGGS radar suggests that the ice shelf has been fairly stable during the last four decades between the two surveys.
Values of H T exceed 500 m near the WAIS grounding lines and major EAIS glacier outlets. The central Ross Ice Shelf has typical values of H T of about 300-500 m. Total ice thickness generally decreases along ice flowlines; that is, thinning due to the combination of basal melting and strain usually exceeds thickening by surface accumulation. We measured values of H T of about 300 m along the ice shelf front, with localized regions of much thinner ice (marked A-E on Figure 5a). Thin ice (<300 m) extends up to 80 km from the front near B and 100 km at E, where it forms a corridor of thin ice along the flow paths of Byrd and Mulock glaciers (locations shown in Figure 1b).

Thickness of the LMI Unit
The map of H LMI (Figure 5b) derived from the SIR shows well-defined bands following flowlines originating at the grounding line of various glaciers such as the Kamb Ice Stream in WAIS, and Beardmore, Nimrod, and Byrd glaciers in EAIS. H LMI generally increases along-flow (Figure 5b) since surface accumulation exceeds strain thinning of the LMI unit for most of the ice shelf. For example, the base of the LMI unit along the Nimrod flowline is close to the surface near the grounding line and deepens to 100 m over the next 300 km, finally reaching a maximum depth of 150 m close to the calving front. The character of the reflector at the base of the LMI unit varies with ice source (Figures 2a and 2b). In the Byrd Glacier band, this reflector is undulating and shallow (Figure 2b). The Mercer Ice Stream band also contains a relatively shallow, but smoother reflector, while that of Kamb Ice Stream is characterized by a deeper, rougher reflector.

Thickness of the CMI Unit
At each location along the survey lines, we calculated H CMI as the difference between H T from DICE, and H LMI from SIR. The map of H CMI (Figure 5c) shows thick CMI units originating at Kamb and Whillans ice streams in the WAIS, and Nimrod, Byrd, and Mulock glaciers in the EAIS. In each case, H CMI generally decreases downflow. For some flowpaths such as for Byrd Glacier, the CMI unit disappears before the ice reaches the ice front. H CMI reaches particularly low values at five locations (A-E) along the ice front ( Figure 5c).

Modern Steady State Strain and Strain-Induced Thickness Change Rates
We estimated the contribution of vertical strain rates (∇·V) to ΔH CMI from the observed velocity ( Figure 1b) and ELMER/Ice modeled velocity fields. Our maps of strain rates (Figures 6a and 6b) (Figures 6a and 6c). The typically small differences between satellite-derived and RIGGS strain rates (Figure 6c) suggest that the fundamental characteristics of the ice shelf's strain field have not changed significantly over the last four decades. We attribute much of the difference between the observed strain rate field and the RIGGS measurements to the different spatial scales of each measurement type: we smoothed the satellite-observed velocities at 10 × 10 km length scales, while RIGGS vertical strain rates were estimated from relative ice motion measured at 1-3 km scales.
An additional source of disagreement between satellite-based strain rate estimates and RIGGS may be associated with significant uncertainties in RIGGS station locations, which were located either by optical sightings of the Sun or transit satellite navigation observations. In regions of large gradients in strain rates, the relatively low accuracy of the RIGGS positions may contribute to the apparent disagreement between the satellite-based and RIGGS strain rate estimates. The mean difference between RIGGS and the modern satellite-derived strain rates is 2.7 × 10 −4 a −1 , although the linear correlation is weak (r = 0.31). Significant differences between the two estimates (Figures 6a and 6c) are observed near Whillans Ice Stream, and closer to the ice shelf front where large crevasses and rifts are visible on the ice surface in REMA images (Figure 1d). We interpret these changes as the combined result of slowing down of Whillans Ice Stream and active rift and calving processes during the four decades between RIGGS and ROSETTA-Ice.
We converted strain rates (units of a −1 ) from observed and modeled velocity fields to strain-induced thickness change (strain thinning or thickening) rates (units of m a −1 ) for the CMI unit by multiplying vertical strain rate by H CMI . These fields (Figures 6d and 6e) show that the strain-induced thinning rates exceed 0.1 m a −1 for most of the Ross Ice Shelf. Smaller strain-induced thinning occurs downflow from Crary and Steershead ice rises, and between the flow paths for Nimrod and Byrd glaciers. Strain-induced thickening was only observed in a few small regions, notably in the compressive flow of ice from Byrd and Mulock glaciers as it passes Minna Bluff and Ross Island (Figure 1b). The convergence zones downstream of Crary and Steershead ice rises (Figures 6a and 6b) do not appear in the maps of strain-induced thickness change (Figures 6d and 6e) because we often cannot trace the CMI-LMI interface there.  (Figures 6b and 6e) are based on an ice-sheet model for which steady state is assumed. As shown by Thomas et al. (2013), however, there are significant differences in the velocity fields from RIGGS compared with those from the satellite era, consistent with the expected ice shelf response to changes in inflowing WAIS ice streams. We investigated the effects of a nonsteady ice flow on basal melting using a numerical simulation of past flow events inferred from streakline distortions for 1,000 years (Hulbe & Fahnestock, 2007). Thickness transients due to the stagnation of Kamb Ice Stream~180 years ago and ongoing flow changes of Whillans Ice Stream (Figure 1d) are of particular interest here due to their magnitudes, but signals from other, older events also persist (Hulbe et al., 2013;Campbell et al., 2018).
We quantified the potential contribution of this nonsteady-state term to our estimated basal melt rates by first evaluating thickness change in a model scenario (Hulbe & Fahnestock, 2007). A steady-state model was perturbed with the velocity changes associated with a stopping of the Kamb Ice Stream and the stopping, starting and eventual slowing of the Whillans Ice Stream in the course of 1,000 model years. The modeled values of ΔH T were then scaled by the fraction H CMI /H T to obtain the thickness change ΔH CMI (Figure 6f) caused by nonsteady-state processes on grounded ice. Given that substantial velocity changes could have occurred over the four decades between RIGGS and ROSETTA-Ice surveys, and that this is a typical averaging time for our basal melt rate estimates (Figure 4b), the nonsteady contribution to strain-induced thickness change rates may be significant near Crary Ice Rise. However, Figure 6f indicates the effect of transients on basal melt rates for most of the ice shelf is small.

Basal Melt Rates
We estimated the basal melt rate, w b (Figure 7a) from equation (4) by applying the corrections for straininduced thickness change rates discussed in section 4.2. The basal melt rates represent an average over the time taken (Δt) for ice to flow between adjacent ROSETTA-Ice survey lines; in the range of about 10-

10.1029/2019JF005241
Journal of Geophysical Research: Earth Surface 60 years (Figure 4b). The melt rate at the ice shelf front is for an average time of~10-20 years and up to~60 years in the interior. Uncertainties of w b over most of the ice shelf are less than 0.2 m a −1 (Figure 8).
The highest values of w b , exceeding 0.5 m a −1 , were found along the ice front, roughly collocated with the regions of relatively small H T (Figures 5a and 7a) and H CMI (Figure 5c) identified on these figures as Locations A-E. A transect of w b along the ice front with 20 km spatial averaging (Figure 7c) shows elevated melt rates from Ross Island to about 170°W. The most prominent region of high melt rates (Site E) includes an extensive region with w b~2 -2.5 m a −1 extending southward from the ice front just east of Ross Island to south of Minna Bluff, along the flowlines coming from Byrd and Mulock glaciers (Figure 1b).
South of the ice front, there are extensive regions of negative values of w b suggesting marine ice accumulation, and whose absolute values often exceed the formal errors of our method (Figure 8a). Some marine ice is known to be present: A 6 m layer of marine ice was observed at the base of a hole drilled through the ice shelf during RIGGS at Site J-9 north of Crary Ice Rise (Zotikov et al., 1980). These authors estimated, however, that the mean Lagrangian basal accumulation rate for ice flow leading to this site was only 0.02 m a −1 . Thomas and Bentley (1978) estimated that, for a region south of Crary Ice Rise, during the RIGGS program the ice shelf should have been thickening at a rate of about 0.35 ± 0.15 m a −1 . However, they attributed this thickening to ice convergence associated with ice shelf grounding on Crary Ice Rise rather than basal ice freezeon, which they assumed was negligible.

Validation of the Radar-Based Method for Evaluating Basal Melt Rates
We compared the basal melt rates from our method and from the Lagrangian analysis of ICESat (2003ICESat ( -2009) laser altimetry (Moholdt et al., 2014). The gridded, ICESat-derived values were averaged over a 2,000 m radius (Table 1) centered around each of our radar-derived basal melt points (Figure 9a). In general, our estimates are consistent with the rates from Moholdt et al. (2014) (Figure 9b) within the typical uncertainty of our radar-based method (±0.1-0.3 m a −1 , Figure 8b). At the ice shelf front, our radar-based basal melt rates were usually lower than the ICESat-based values. However, our radar-derived values show larger areas of basal melting compared with the ICESat grid, partly due to the higher melt rates estimated near Crary Ice Rise (Figure 9a). It is possible that the higher melt rates estimated near Crary Ice Rise are a result of Whillans Ice Stream slowing down such that strain-induced thickness change is nonsteady, causing inferred higher basal melt rates in our estimates. In the Byrd-Mulock corridor, our melt rates are lower than the ICESat-derived values. This may indicate that there was more melting in the ICESat era (2003)(2004)(2005)(2006)(2007)(2008)(2009) compared with our longer (~20 years) average; however, the differences are very close to our noise levels ( Figure 8). Close to the western edge of Roosevelt Island, our melt rate estimates of 0.3-0.4 ± 0.2 m a −1 agree with values estimated by Catania et al. (2010) based on drawdown of internal layers in ground-based radar profiles.

Possible Formation Mechanism for the CMI-LMI Interface
Coherent internal reflectors that can be traced for long distances have often been used for calculating surface mass balance and basal melt rates, as long as they could be dated (e.g., Catania et al., 2010;Das et al., 2015). However, identifying the specific cause of the reflectors provides more confidence in our interpretation of them as age markers. In this work, we have assumed that the set of bright coherent reflectors that form  (Scambos et al., 2007), also showing ice front sites of melting hot spots (A-E), ice shelf front bathymetric contours extending to the continental shelf (gray lines; Tinto et al., 2019), maximum ocean temperature (T max ) at depth (z) >120 m, primarily the mCDW affecting melt rates at B, the 250 m ice draft contour (white lines). (b) Modeled summer ocean temperature along the ice shelf front (following Tinto et al., 2019). (c) Basal melt rates averaged over 20 km along the ice shelf front (dark lines), and standard deviation of the melt rates (gray lines).

10.1029/2019JF005241
Journal of Geophysical Research: Earth Surface the CMI-LMI interface within each outlet glacier and ice stream are caused by surface or englacial processes happening on the grounded ice or close to the grounding line. On the EAIS side, the subsequent burial of wind-scoured ice surface resulting from strong katabatic winds (Das et al., 2013;Scambos et al., 2012) channeled through the outlet glaciers may cause the bright reflectors. Satellite imagery showing regions of blue ice in the Transantarctic Mountains (e.g., Hui et al., 2014) support this hypothesis. Compressive forces on a small glacier (Skelton Glacier) in the Transantarctic Mountains may also contribute to the exposed blue ice there (Crary & Charles, 1961). On the WAIS side, compressed ice and buried crevasses likely produce an internal layer with sufficient dielectric contrast for it to be identified in the radar images (Wright et al., 1990). Variations in the roughness of the CMI-LMI interface between glacier packets (Figure 2b) are likely due to a combination of surface processes and different stress regimes of each glacier.

Causes of Spatial Variability of Melt Rates
The amount of basal melting depends on the heat content of the underlying ocean determined by the ocean circulation under the ice shelf, and the intensity of mixing near the ice base to maintain an upward flux of ocean heat against the stabilizing effect of buoyant meltwater production (Holland & Jenkins, 1999). For the Ross Ice Shelf, the sub-ice-shelf circulation is described in detail by Assmann et al. (2003), Dinniman et al. (2011Dinniman et al. ( , 2016, Jendersie et al. (2018), Tinto et al. (2019), and others. The low basal melt rates, and possible marine ice accumulation, that we estimate under most of the ice shelf ( Figure 7a) is consistent with previous satellite-derived estimates (Moholdt et al., 2014;Rignot et al., 2013). Numerical simulations of ocean circulation, including those of Tinto et al. (2019) based on their updated seabed bathymetry under the Ross Ice Shelf, indicate that low rates are associated with the dominance of cold High Salinity Shelf Water in the ocean flows into the sub-ice-shelf cavity. Production of cold Ice Shelf Water released by melting of the ice base in contact with High Salinity Shelf Water near the deep grounding lines of the major EAIS glaciers further reduces basal melting for most of the ice shelf.
Near Crary Ice Rise, the ocean models do not predict the relatively high melt rates suggested by our analyses of ROSETTA-Ice radar data. We suggest that this discrepancy is caused by transient changes in ice thickness from past dynamic changes in the grounded ice streams (Catania et al., 2012). In these places, it is likely that the ice column is not in steady state, and variations in velocity streamlines over the past decades (e.g., Figure 1d) may drive anomalies in strain rates. The difference between steady-state and nonsteady-state thickness change (Figure 9c) shows larger values downstream of Kamb and Whillans Ice Stream, leading to larger uncertainties in our estimated melt rates (Figure 8a). The ice in front of Roosevelt Island also shows large differences between steady and nonsteady thickness change rates (Figure 9c). We propose that this anomaly is caused by large rifts and crevasses on that section of the ice shelf (see Figure 1c).
The observed heightened melting along most of the ice front (Figure 7) has been previously identified from analyses of satellite laser altimetry (Horgan et al., 2011;Moholdt et al., 2014). Direct measurements of basal melt rate near Ross Island (Stewart et al., 2019) support the hypothesis that long-term melt rates at Site E are dominated by rapid melting during summers when Antarctic Surface Water (AASW) has been warmed locally by insolation after the sea ice near the ice front has disappeared (e.g., Assmann et al., 2003;Stern et al., 2013;Stewart et al., 2019;Tinto et al., 2019). The subsequent advection of warmed AASW under the outer ice shelf may be driven by a combination of numerous processes including tides (Arzeno et al., 2014;MacAyeal, 1984;Makinson & Nicholls, 1999), eddies and other instabilities of the westward flowing The melting hotspot at Site B is associated with the southward subsurface flow of modified Circumpolar Deep Water (mCDW) along Hayes Bank to the ice front; this stream of mCDW also undergoes a seasonal cycle driven by changing wind fields and variations in the rate of heat loss as the upper ocean cools and becomes more saline through sea ice growth in winter (Dinniman et al., 2011).
The melting hotspot at Site D is a region of the ice front where the depth of the warm AASW layer is similar to the ice draft ( Figure 7b). Therefore, there may be a direct path for AASW under the ice shelf, similar to the process occurring at Site E. Hot spots at Sites C and D are west of embayments in the ice front, which may assist in diverting the westward flow of the ice front current under the ice shelf, a process similar to the topographic control of mCDW flows across the shelf break documented by Dinniman et al. (2003).
We have no explanation for the elevated melting at Site A at this time but speculate that it may point to a source of oceanic heat from the Antarctic Coastal Current passing Cape Colbeck. If so, it suggests potential for future changes in melting of the Ross Ice Shelf near the Marie Byrd Land coast, if the coastal current was to undergo significant changes in future climate as modeled for Filchner Ice Shelf in the eastern Weddell Sea (Hellmer et al., 2012). (c) Difference in dynamic thinning rates between steady-state and nonsteady-state transients indicating regions where our steady state assumption may not hold on the WAIS side of the ice shelf. Also note large differences on the ice shelf in front of Roosevelt Island where large crevasses and rifts are located (see Figure 1c).

Conclusions
We have developed a novel method for estimating multidecadal average basal melt rates from airborne ice penetrating radar surveys. Our method leverages the ability to locate, in the radar data, a strong internal reflector for each outlet glacier, that can be traced from the grounding line to the ice front. For the Ross Ice Shelf, which was surveyed by the ROSETTA-Ice project during 2015-2017, we found a strong and spatially coherent reflector between an upper unit of Local Meteoric Ice and a lower unit of Continental Meteoric Ice (CMI). The Local Meteoric Ice unit, which we assume is composed of snow falling onto the relatively flat surface of the ice shelf, has almost horizontal internal layers that were imaged by our shallow ice radar. The CMI unit, which likely originates from the grounded part of the ice sheet, is almost featureless in our deep-ice radar. We used thickness changes of the lower CMI unit, corrected to account for straininduced thickness change, to infer the integrated basal melt rate along ice flow lines between pairs of survey lines.
Over most of the ice shelf, the estimated basal melt rate is close to zero, within the uncertainties associated with radar signal processing and the known, but poorly resolved, temporal variability of the ice shelf's velocity and strain rate fields. Higher values of melt rates are found near the ice front, centered on five hot spots; melt rates at these hotspots often exceed 0.5-2 m a −1 . This distribution of melt rates is similar to previous estimates based on laser altimetry, confirming that the new method is correctly identifying regions of rapid melting. However, in the ROSETTA-Ice survey, the radar-derived melt rates represent averages over one to six decades, a longer time scale than for currently available satellite-derived products. These results suggest that we may be able to use comparisons between satellite and radar-derived estimates to identify whether short-term melt rates are representative of long-term trends. Our new estimates indicate long-term persistence of the oceanographic features driving current melting of the Ross Ice Shelf, including seasonal warming of the upper ocean layer of AASW and intrusions of mCDW along Hayes Bank.
This study highlights the value of airborne surveys with ice penetrating radar systems for directly estimating basal melt rates without the need for models of surface mass balance, ice density profiles, and firn compaction rates. The method generates melt rates averaged over multidecadal timescales set by ice speeds and survey-line separation. Although significant sources of error remain, the new method provides a complementary view of basal melting for large ice shelves and may provide a mechanism for validating developing models of ice-shelf dynamic responses to processes such as changing flows of grounded ice streams and glaciers and large-scale iceberg calving events.