Detecting high spatial variability of ice shelf basal mass balance, Roi Baudouin Ice Shelf, Antarctica

Abstract. Ice shelves control the dynamic mass loss of ice sheets through buttressing and their integrity depends on the spatial variability of their basal mass balance (BMB), i.e. the difference between refreezing and melting. Here, we present an improved technique – based on satellite observations – to capture the small-scale variability in the BMB of ice shelves. As a case study, we apply the methodology to the Roi Baudouin Ice Shelf, Dronning Maud Land, East Antarctica, and derive its yearly averaged BMB at 10 m horizontal gridding. We use mass conservation in a Lagrangian framework based on high-resolution surface velocities, atmospheric-model surface mass balance and hydrostatic ice-thickness fields (derived from TanDEM-X surface elevation). Spatial derivatives are implemented using the total-variation differentiation, which preserves abrupt changes in flow velocities and their spatial gradients. Such changes may reflect a dynamic response to localized basal melting and should be included in the mass budget. Our BMB field exhibits much spatial detail and ranges from −14.7 to 8.6 m a−1 ice equivalent. Highest melt rates are found close to the grounding line where the pressure melting point is high, and the ice shelf slope is steep. The BMB field agrees well with on-site measurements from phase-sensitive radar, although independent radar profiling indicates unresolved spatial variations in firn density. We show that an elliptical surface depression (10 m deep and with an extent of 0.7 km × 1.3 km) lowers by 0.5 to 1.4 m a−1, which we tentatively attribute to a transient adaptation to hydrostatic equilibrium. We find evidence for elevated melting beneath ice shelf channels (with melting being concentrated on the channel's flanks). However, farther downstream from the grounding line, the majority of ice shelf channels advect passively (i.e. no melting nor refreezing) toward the ice shelf front. Although the absolute, satellite-based BMB values remain uncertain, we have high confidence in the spatial variability on sub-kilometre scales. This study highlights expected challenges for a full coupling between ice and ocean models.


Introduction
Approximately 74 % of the Antarctic ice sheet is surrounded by floating ice shelves (Bindschadler et al., 2011a) providing the interface for interactions between ice and ocean. Marine ice sheets -characterized by a bed elevation below sea level and sloping down towards the interior -can be destabilized leading to a marine ice sheet instability (Mercer, 1978;Schoof, 2007;Tsai et al., 2015). However, ice shelves that are laterally constrained through embayments (or locally regrounded from below), mitigate the marine ice sheet instability , thus regulating the ice flux from the inland ice sheet through buttressing. Over the last decade, major advances in our understanding of the processes at this ice-ocean interface have emerged, both theoretically (e.g. Pattyn et al., 2013;Favier et al., 2012Favier et al., , 2014Ritz et al., 2015) and from observations (e.g. Rignot et al., 2014;Wouters et al., 2015). It is now established that ice shelf integrity plays an important part in explaining sea-level variations in the past (Golledge et al., 2014;DeConto and Pollard, 2016), enabling improved projections of future sea-level rise (Golledge et al., 2015;Ritz et al., 2015).

S. Berger et al.: Spatial variability of ice shelf basal mass balance
Ice shelf integrity can be compromised by atmosphericdriven surface melt ponding (Lenaerts et al., 2017) and hydrofracturing (Banwell et al., 2013;Scambos et al., 2004;Hulbe et al., 2004). From the ocean side, ice shelves may thin or thicken (Paolo et al., 2015) due to changes in basal mass balance (BMB), i.e. the difference between refreezing and melting. Point measurements with phase-sensitive radars (Marsh et al., 2016;Nicholls et al., 2015), global navigation satellite system (GNSS) receivers (Shean et al., 2017), observations from underwater vehicles (Dutrieux et al., 2014) and analysis from high-resolution satellites (Dutrieux et al., 2013;Wilson et al., 2017) have shown that BMB varies spatially on sub-kilometre scales. Ice shelf channels are one expression of localized basal melting (Stanton et al., 2013;Marsh et al., 2016) which, after hydrostatic adjustment, form curvilinear depressions visible at the ice shelf surface (Fig. 1). These surface depressions reflect basal incisions resulting in curvilinear tracts of thin ice. In some areas, ice shelf channels are twice as thin as their surroundings (Drews, 2015). However, the impact of ice shelf channels on ice shelf integrity is yet unclear because, on the one hand, excessive basal melting beneath ice shelf channels may prevent ice-shelf-wide thinning (Gladish et al., 2012;Millgate et al., 2013) but, on the other hand, increased crevassing due to channel carving may structurally weaken the ice shelf .
Here we attempt to derive the BMB of the Roi Baudouin Ice Shelf (RBIS), Dronning Maud Land, East Antarctica, at 10 m gridding, based on mass conservation in a Lagrangian framework. The RBIS (Fig. 1) is constrained by an ice promontory to the west and by Derwael Ice Rise in the east, blocking the tributary flow from West Ragnhild Glacier -one of the largest outlet glaciers in Dronning Maud Land . Analyses on Derwael Ice Rise Callens et al., 2016) and the larger catchment area  suggest that the RBIS is a relatively stable sheet-shelf system on millennial timescales. The RBIS contains a number of ice shelf channels (Drews, 2015, and arrows in Fig. 2e), many of which start at the grounding line and extend over 230 km to the ice shelf front.
We outline our approach of deriving the BMB, with focus on attaining high spatial resolution. Resolving BMB is challenging, because it is computed as the residual of several uncertain quantities and it relies on spatial derivatives, which amplify noise in the input data. The latter can be accounted for with spatial averaging (e.g. Neckel et al., 2012;Moholdt et al., 2014), which, however, may smear out the imprint of processes acting on sub-kilometre scales. Here, we use spatially well-resolved input data combined with total-variation regularization of the velocity gradients. This avoids spatial averaging, but still mitigates the noise in the input data. As a result, our BMB field shows high detail over different spatial scales that are validated with phase-sensitive radar, GNSS observations and ground-penetrating radar.  Fig. S2. Light blue and light red are the low-lying parts of the ice shelf, which are excluded from the GNSS-TanDEM-X comparison in Fig. 9. "Radar" denotes both ground-penetrating and phasesensitive radars. The background is from the Radarsat mosaic (Jezek and RAMP Product Team, 2002) and the black line delineates the grounding line (Depoorter et al., 2013b).
2 Data and methods

Basal mass balance from mass conservation
We derive the basal mass balance (Ṁ b ) from mass conservation, i.e.
whereṀ s is the surface mass balance (SMB, positive values for mass gain), H i is the ice thickness and u the columnaverage horizontal velocity of the ice.Ṁ b ,Ṁ s and H i are given in ice-equivalent units. ∂H i /∂t and DH i /Dt represent the Eulerian and Lagrangian thickness change, respectively, and ∇ · (H i u) denotes the flux divergence, which includes advection of thickness gradients (u · ∇H ) and ice flow divergence (H i ∇ · u). In principle, Eq. (1) does not depend on the reference frame and can be calculated in both a fixed coordinate system (i.e. Euler coordinates) or with a moving coordinate system that follows the ice flow (i.e. Lagrange coordinates). In practice, however, both approaches differ: Eulerian studies are often based on one thickness field and either assume a steady state (Rignot and Steffen, 2008;Neckel et al., 2012) or rely on an external dataset  (Jezek and RAMP Product Team, 2002) and the black line delineates the grounding line (Depoorter et al., 2013b). Key features regarding the input datasets are summarized in Table 1. (Depoorter et al., 2013a, Rignot et al., 2013 to account for the thickness changes ∂H i /∂t (e.g. Pritchard et al., 2012;Paolo et al., 2015). The Lagrangian approach, on the other hand, requires two thickness fields covering different time periods from which the Lagrangian thickness change is calculated implicitly (DH i /Dt). As shown below, the key difference between both approaches is how the advection of thickness gradients (u·∇H ) is accounted for. The Lagrangian  (Förste et al., 2014) (14.6;19.8) approach is best suited in areas with rough surface and significant advection (e.g. near ice shelf channels). We refer to previous publications (Dutrieux et al., 2013;Moholdt et al., 2014;Shean et al., 2017) that further explain differences between Eulerian and Lagrangian approaches.
In the following, we describe surface velocities in Sect. 2.2, surface mass balance in Sect. 2.3, the derivation of hydrostatic ice thickness in Sect. 2.4 and Lagrangian thickness change in Sect. 2.5. Key features of the input datasets are summarized in Table 1. As a novelty compared to previous studies, we base our hydrostatic thickness field on highresolution digital elevation models (DEMs) derived from TanDEM-X images from 2013 and 2014. Section 2.6 explains the implementation of spatial velocity gradients (∇ · u in Eq. 1), which is non-trivial when derivatives are taken over short distances with noisy input data. We compare the derived Lagrangian basal mass balance (LBMB) with field measurements of phase-sensitive radar and GNSS profiling (Sect. 2.7). Although this is not a direct validation, as the field data cover a different period, the comparison is insightful to understand the spatial variability in our LBMB field. The derived LBMB is only valid in freely floating areas, which excludes not only the grounding zone but also other small-scale features such as ice shelf channels where viscous inflow can occur (Humbert et al., 2015;Drews, 2015). (Examples where this may be the case are discussed in Sect. 5.)

Surface velocities from satellite radar remote sensing
Assuming that velocities do not vary with depth, we use surface velocities that were previously derived by combining interferometric synthetic aperture radar (InSAR) and speckle tracking . The velocities are mosaicked and gridded to a 125 m posting and are based on images from the European Remote Sensing satellites (ERS 1/2) from 1996 and the Advanced Land Observing System Phased Array type L-band Synthetic Aperture Radar (ALOS-PALSAR) from 2010. As shown in Berger et al. (2016), comparison with on-site measurements collected in 1965-1967 and 2012-2014 yields no evidence of prominent changes in the ice velocities over the last decades, which supports the combination of data from different dates. The velocity mosaic covers 75 % of our area of interest (dashed line in Fig. 2c). The remaining areas are filled with an Antarctic-wide flow field (Rignot et al., 2011b) gridded to 900 m postings (the 450 m gridded velocities are too noisy in our area of interest). We reduce seams -as high as 60 m a −1 in some places -using linear feathering (e.g. Joughin, 2002;Neckel et al., 2012) over 4.5 km.

Surface mass balance from atmospheric modelling
We use the surface mass balance from a high-resolution (5.5 km posting) simulation of the Regional Atmospheric Climate MOdel (RACMO) version 2.3, centred on Dronning Maud Land (25 and 45 • W) and averaged over the period 1979-2015 (Lenaerts et al., 2014(Lenaerts et al., , 2017. The SMB field correctly reproduces asymmetries across Derwael Ice Rise originating from orographic uplift and also simulates a corresponding shadowing effect on the Roi Baudouin Ice Shelf ( Fig. 2b and Lenaerts et al., 2014). Moreover, the simulation explains observed surface melting near the grounding zone due to a wind-albedo feedback caused by persistent katabatic winds in this area (Lenaerts et al., 2017).

Hydrostatic ice thickness
We calculate the ice thickness ( Fig. 2e) by imposing hydrostatic equilibrium on surface freeboard (Bindschadler et al., 2011b;Chuter and Bamber, 2015;Drews, 2015) derived from the TanDEM-X satellites. The details of hydrostatic inversion are presented in the following two sections.

Surface elevations
The digital elevation models are processed from 43 image pairs ( Fig. S1 in the Supplement) of the TanDEM-X mission (Krieger et al., 2007), in which the TerraSAR-X and TanDEM-X satellites image the surface simultaneously from different viewing angles. This allows us to infer topography interferometrically without the need to correct for ice flow. Images from the austral winters of 2013 and 2014 are processed to single-look complex scenes, using SARscape ® . After coregistration using the CryoSat-2 DEM (Helm et al., 2014), the pairs are differenced in phase. The resulting interferograms are then unwrapped and the phase difference is re-flattened before being geo-referenced in polar stereographic coordinates. The processing provides 43 single DEMs (32 from 2013 and 11 from 2014) gridded to 10 m. They cover time spans of June-October 2013 and June-July 2014 (Fig. S1), with time separations ranging from 231 to 379 days in areas where thickness rates where calculated between the 2 years. Digital elevation models from the same date and satellite path are concatenated together, with a linear taper on overlapping zones. Grounded areas are masked out using the composite grounding line from Depoorter et al. (2013a), based on differential InSAR with Radarsat and PALSAR (Rignot et al., 2011a) at RBIS. To correct for small elevation shifts between the different frames, which we assume to be uniform over the ice shelf, we tie the 2013 concatenated frames to each other and to the CryoSat-2 DEM (Helm et al., 2014), using constant offsets. We attribute these small shifts to tides, inverse barometric effects or different calibrations during the SAR processing.
All DEMs are smoothed with a Gaussian filter to remove small-scale surface roughness. The standard deviation of the filter is set to 7 pixels (or 70 m) in all directions. This means that points lying within that distance are weighted with 0.68. To determine the size of the Gaussian filter, we investigated standard deviations from 1 to 10 pixels and found that using 7 pixels minimizes the elevation discrepancy between 2012 GNSS and TanDEM-X surface elevation (Sect. 2.7). As shown in Fig. S2, the applied smoothing does not affect the shape of the surface depressions linked to ice shelf channels (with a typical width of 1-2 km).
The difference fields of the 2013-2014 overlapping DEMs exhibit a linear trend aligned with the satellite trajectory. We attribute this signal to the interferometric processing, which can leave a flawed elevation trend due to imprecise information about the satellite orbits or due to ill-constrained parameters during the SAR processing (Drews et al., 2009). To account for this effect, we subtract a plane from the 2014 DEM using the difference fields of 2013-2014 overlapping fields. The plane fit and the offset correction applied earlier mask absolute ∂H i /∂t changes, which we assume to be small in the following.
To assess the relative vertical accuracy of the final DEMs (Sect. 4.1) (i) we use the difference fields of overlapping, unconcatenated TanDEM-X frames from the same date and satellite path (Fig. S3), and (ii) we compare the DEMs to kinematic GNSS profiling. We estimate the relative vertical accuracy to be better than 1 m, although elevation differences in some areas are systematically higher (Sect. 2.7). The offset and plane fitting corrections are further discussed in Sect. 4.1, as they strongly impact the quality of our ice-thickness fields and the resulting LBMB rates.

Hydrostatic equilibrium
We invert hydrostatic thickness from freeboard heights (h asl ) with densities of ρ w = 1027 kg m −3 , ρ i = 910 kg m −3 and ρ a = 2 kg m −3 , for seawater, ice and firn air, respectively: The firn air content H a accounts for the lower firn and snow densities by subdividing the ice column in air-and iceequivalent layers. We use simulated values from the firndensification model "IMAU-FDM" (Fig. 2a and Ligtenberg et al., 2011;Lenaerts et al., 2017), which is forced by the SMB, exists on the same spatial grid (5.5 km, Sect. 2.3) and is averaged over the same time period . For converting ellipsoidal heights to freeboard elevations we employ the EIGEN-6C4 geoid (Förste et al., 2014) and the DTU12MDT mean dynamic topography model (Knudsen and Andersen, 2012). The hydrostatic ice thickness is most sensitive to the firn air content and the freeboard heights, resulting in an estimated uncertainty of at least ±25 m (Drews, 2015). However, as discussed in Sect. 4.1, uncertainties can be much higher in areas where firn density is ill-constrained.

Lagrangian thickness change
As the Lagrangian framework moves with the flow, computing the Lagrangian thickness change DH i /Dt requires us to shift one thickness field to match the geometry of the second one. Consequently, this approach implicitly accounts for advection of thickness gradients (u · ∇H i ). Here, the 2013 TanDEM-X frames are shifted forward with a normalized correlation coefficient matching algorithm from the computer vision library OpenCV (Bradski and Kaehler, 2008). Each 2013 concatenated frame is divided into 5 km × 5 km patches that are sampled every kilometre in both directions. Each 2013 patch is then compared with any possible 5 km × 5 km patch within a slightly bigger search region (6.6 km × 6.6 km) in the 2014 DEMs that overlap with the 2013 DEM. Comparison is based on normalized crosscorrelation coefficients technique, a more robust variant of 2-D normalized cross-correlation (Marengoni and Stringhini, 2011). The shift of the 2013 patches is found where the correlation coefficient is maximal. Mismatches are discarded when the correlation coefficient is smaller than 0.8, or when the detected offset is well beyond what would be expected from the available flow field. All the 2013 shifted patches are then mosaicked to construct a shifted 2013 frame that matches the geometry of its overlapping 2014 frame. The process is applied to each overlapping pair of 2013-2014 TandDEM-X frames before conversion to hydrostatic thickness. In Sect. 4.1, we investigate an alternative approach using observed surface velocities to shift the DEMs with a 10-day time step (as in Moholdt et al., 2014). We also apply this alternative approach to shift the 2016 GNSS profiles (Sects. 2.7, 3.2 and 5).

Spatial derivatives of noisy input data
Taking spatial gradients in Eq. (1) is not straightforward as naive discretization schemes (e.g. forward, backward or central differences) greatly amplify the signal-to-noise ratio if the input data are noisy. This issue can be accounted for by smoothing the input data (e.g. Moholdt et al., 2014) and/or by increasing the lateral distances over which the derivative is approximated (e.g. Neckel et al., 2012). However, smoothing prior to taking the derivative can lead to smearing out of the derivative in areas where the derived quantity changes abruptly (or discontinuously). We expect such abrupt changes in the surface velocities across ice shelf channels that experience strong basal melting (Drews, 2015). To circumnavigate this problem, we applied the total-variation regularization, a technique that suppresses noise from spatial derivatives while preserving abrupt changes (Chartrand, 2011). Noise is removed from the data by reducing the total variation in the signal to a certain degree controlled by a regularization parameter α (Chartrand, 2011). The α value we use (10 5 ) is given by the variance of the velocities, following the discrepancy principle (Chartrand, 2011). Figure 3 compares regularized derivatives with derivatives based on velocity fields that were smoothed to varying degrees prior to taking the derivatives using central differences. It should be noted that some ambiguity about the specific choice of α remains but this is inherent to regularization in general. We discuss the benefits and trade-offs of the different derivative schemes further in Sect. 4.2.

On-site geophysical measurements
Remote-sensing and modelling data are complemented by a series of geophysical measurements (ground-penetrating radar, GNSS profiling and phase-sensitive radar measurements) carried out in December 2012, December 2014 and January 2016 (Fig. 1).
The ground-penetrating radar profile shown in Fig. 8a (located in Fig. 1) was acquired in 2016 with a 20 MHz pulsed radar (Matsuoka et al., 2012). The data are geolocated with kinematic GNSS and migrated using Kirchoff depth migra- tion with a velocity-depth function that accounts for the low firn air content in this area. More details about acquisition and processing of the radar data are given in Drews et al. (2015). We use the radar ice thickness to validate the hydrostatic ice thickness (Sect. 4.1).
We use three sets of kinematic GNSS profiles that were recorded at 1 Hz intervals with geodetic, multi-channel receivers moving at a speed below 12 km h −1 . In December 2012, a 20 km × 25 km GNSS network was acquired at the front of the ice shelf (Drews, 2015). The profiles cross ice shelf channels multiple times. Two years later in December 2014, a 100 km long north-south GNSS transect was acquired (Lenaerts et al., 2017). The last GNSS dataset was acquired in January 2016, along and across an elliptical surface depression (Sect. 3.2). All GNSS elevations are de-tided using the circum-Antarctic tide model (CATS2008a_opt) from Padman et al. (2002Padman et al. ( , 2008. Datasets from 2012 and 2016 are processed differentially, relative to a non-moving base station , while data from 2014 are post-processed with Precise Point Positioning. Elevations from GNSS are used (i) to determine the size of the Gaussian filter applied to the TanDEM-X DEMs (2012 survey, Sect. 2.4), (ii) to assess the accuracy of the TanDEM-X DEMs (2012 and 2014 surveys, Sect. 4.1) and (iii) to extend the time period of surface elevation change detected by the TanDEM-X mission (2016 survey, Sects. 3.2 and 5; Figs. 7 and 8).
BMB was measured at point locations using a phasesensitive radar. Processing and acquisition schemes are as outlined previously (Nicholls et al., 2015;Marsh et al., 2016). The radar antennas were positioned at 22 sites. Each site was remeasured after 10 days at the same location at the surface (in a Lagrangian framework). This way, relative thickness changes due to strain thinning and basal melting can be detected within millimetres. Strain thinning is corrected using a linear approximation of the vertical strain rate with depth, based on tracking the relative displacement of internal reflectors. The strain correction of the BMB rates is small (6.6 × 10 −3 a −1 on average), because strain thinning is small.

Large-scale pattern of the basal mass balance
The LBMB rates range from −14.7 to 8.6 m a −1 (excluding outliers with 0.1 and 0.99 percentiles) and average −0.8 m a −1 (negative values signify melting, positive values refreezing). For the 9227 km 2 covered by the TanDEM-X DEMs, net mass loss at the ice shelf bottom is 6.7 Gt a −1 . Most melting occurs just seaward of the grounding zone where the western Ragnhild Glacier feeds into the Roi Baudouin Ice Shelf (Fig. 4, label A). This area corresponds to the thickest and fastest part of the grounding zone ( Fig. 2e and c). We also find elevated melting close to the western ice promontory (Fig. 4, label C) and on the southern side of Derwael Ice Rise (Fig. 4, label B).
The uncertainties of the absolute LBMB are typically higher than the LBMB itself, because errors unfavourably propagate in mass budgets (Moholdt et al., 2014). Here, we assess a lower bound of the LBMB errors by using the difference fields of the individual LBMB frames in overlapping areas. These show no systematic patterns and the standard deviation amounts to 2.3 m a −1 . Moreover, comparing the (yearly averaged) LBMB values with the 22 on-site phase-sensitive radar measurements, reveals differences of 1.1 ± 2.6 m a −1 in mean and standard deviation, respectively. We discuss this comparison in more detail in Sect. 3.2. Figure 2b, d and f illustrate the terms entering Eq. (1c), namely surface mass balance, ice flow divergence and Lagrangian thickness change, whereas Fig. 2a, c and e display the most critical input vari-ables needed to compute those different terms -i.e. firn air content, ice velocity and hydrostatic thickness. For the RBIS, the Lagrangian thickness change dominates the BMB (as in Shean et al., 2017), while ice flow divergence and SMB are both one order of magnitude lower. Qualitatively the largescale pattern agrees well with the results from Rignot et al. (2013) who also found the highest melt rates close to the grounding line, both for steady state or transient approximations.
To illustrate the advantages of the Lagrangian approach, Fig. 5 shows the Eulerian thickness change, flux divergence and Eulerian BMB. While the large-scale pattern of the Eulerian BMB agrees very well with that of the LBMB, the Eulerian approach fails in the vicinity of ice shelf channels (arrows in Fig. 5). Advecting topographic features imprint the Eulerian thickness changes (Fig. 5a), however, the Eulerian approach does not fully account for this advection of thickness gradients (u · ∇H i ) in the flux divergence (Fig. 5b). This results in spurious Eulerian BMB in the vicinity of ice shelf channels (Fig. 5c). These spurious signals in the Eulerian BMB become even stronger when thinning/thickening rates are taken from external datasets which are spatially less well resolved. Using ice-shelf-wide, average values (e.g. repeat satellite altimetry) does not account for the advection of ice shelf channels and other (transient) features in the ice shelf. It therefore introduces artifacts into the basal mass balance pattern.

Small-scale variability of the basal mass balance
The larger-scale LBMB pattern (> 10 km) is overlain by smaller-scale variability. Ice shelf channels appear most clearly in the DEMs and thus in the hydrostatic thickness fields (arrows in Fig. 2e). In some places, they also co-locate with areas of lateral inflow (i.e. negative flow divergence; arrows in Fig. 2d) and Lagrangian thinning (i.e. negative Lagrangian thickness change in Fig. 2f). In the LBMB, ice shelf channels appear partially as narrow bands of intense melting. Figure 6 shows one example where ice preferentially melts at the flanks of an ice shelf channel. LBMB rates drop to −5 m a −1 at both flanks, whereas outside the channel the LBMB is close to zero. The slight refreezing found at the channel's apex (1.5 m a −1 ) is very close to the detection limit and its magnitude is 3 times lower than what is observed at the flanks.
Another example of a small-scale feature is illustrated in Figs. 7 and 8. Here, we observe a 0.7 km × 1.3 km elliptical surface depression that is up to 10 m lower than its surroundings and located on the upstream end of an ice shelf channel. The surface topography also exhibits secondary elongated surface depressions that are shaped like fingers merging into the elliptical depression. We surveyed this area in 2016 with kinematic GNSS profiles, ground-penetrating radar and 22 point measurements of the BMB with phasesensitive radar (Sect. 2.7). Lenaerts et al. (2017)    this feature as one of the 55 features on the Roi Baudouin Ice Shelf, that can be linked to the formation of englacial lakes near the grounding line. They proposed that these features are initially formed as supra-glacial lakes in the grounding zone due to katabatic wind-albedo feedback. Freezing at the lake surface and subsequent burial by snowfalls form at first englacial lakes that are advected farther downstream. As a function of the advection time the liquid water then likely fully refreezes. For the elliptical surface depression considered here, the radar data show a bright reflector at approximately 30 m depth and no coherent signals appear at larger depths (Fig. 8a). We tentatively interpret the bright radar reflector as a refrozen surface of a former supra-glacial lake. The specularity of this interface hinders deeper penetration of the radar signal. However, a more detailed radar analysis is warranted to unambiguously clarify the origin and history of this feature. Here, we restrict ourselves to the elliptical surface depression where we observe significant surface lowering.
The elliptical depression appears prominently in our LBMB field with rates as low as −12 m a −1 (Figs. 7b and  8b). On the eastern side of the depression, the BMB from the phase-sensitive radar (Fig. 8b) agrees well with the LBMB estimate, both methods averaging about −0.5 m a −1 with little spatial variability. On the western side -which contains the finger-shaped surface features -larger differences and variability occur. The differences could reflect the more com-plex topography and/or temporal variations. The large negative LBMB rates in the elliptical depression reflect persistent surface lowering of 0.5 to 1.4 m a −1 . Ice flow divergence is negligible at that location. We extend the time series from the TanDEM-X DEMs to 2016 with the GNSS profiles (Fig. 8c) where we find the same localized lowering. This indicates that the high-resolution TanDEM-X DEMs reliably pick up surface elevation changes on sub-kilometre scales. Some of the finger-shaped surface depressions also show surface lowering, but less pronounced than what is seen in the elliptical depression itself. The flanks of the surface depression are significantly steeper on the eastern compared to the western side. Unlike the elliptical depression, the ice shelf channel located farther downstream does not actively experience melting or refreezing. Away from ice shelf channels or other surface depressions, our assumptions for the LBMB (such as hydrostatic equilibrium) likely hold explaining the comparatively good fit with the phase-sensitive radar measurements. Inside the elliptical depression, the observed surface lowering cannot unambiguously be attributed to basal melting. Regardless of the specific mechanisms causing the surface lowering, this example highlights that much of the small-scale variability seen in the resulting LBMB field can be used to investigate sub-kilometre-scale ice shelf processes that do not necessarily occur at the ice shelf base. a b c Figure 5. (a) Eulerian thickness change (∂H i /∂t), (b) Flux divergence (∇ · (H i u)) and (c) Eulerian basal mass balance (BMB). Arrows point to spurious signal due to advection of ice shelf channels. The background is from the Radarsat mosaic (Jezek and RAMP Product Team, 2002) and the black line delineates the grounding line (Depoorter et al., 2013b). 4 Error sources

Hydrostatic thickness and Lagrangian thickness change
The Lagrangian thickness change is the dominant error source of the LBMB for the Roi Baudouin Ice Shelf, since the magnitude of both ice flow divergence and SMB is 1 order of magnitude smaller (Fig. 2). The Lagrangian thickness change depends (i) on factors controlling the hydrostatic ice thickness, i.e. the surface elevation (above sea level), the seawater and ice densities, the depth of the firn pack and temporal variations thereof, and (ii) on the Lagrangian matching of the DEMs following the ice flow. It should also be clear that our approach is only able to detect basal changes reflected in the surface elevations, because ice thickness is derived from hydrostatic equilibrium.

Calibration and accuracy of TanDEM-X elevations
The interferometric DEMs provide excellent spatial resolution at the cost that they require calibration. It is straightforward to offset the DEMs to account for the relative phase unwrapping using Antarctic-wide DEMs based on altimetry. More challenging are residual phase trends that may originate from imprecise satellite orbits/SAR processing (Drews et al., 2009) or represent unaccounted tilting of the ice shelf surface due to tides. In our case, these trends are near-linear and become evident in the difference fields of overlapping DEMs from both different years and from the exact same date and satellite path. In the former, systematic biases extend in the azimuth direction with residual height differences typically ranging from −0.5 to +0.5 m. Such biases strongly imprint the corresponding LBMB fields resulting in a mosaic with linear trends typically ranging from −10 to +10 m a −1 in the azimuth direction and differences exceeding 13 m a −1 across seams. To account for this, we correct the 2014 DEMs with plane fitting (Sect. 2.4). We do not correct for systematic trends in individual TanDEM-X frames from the same dates, not only because the small overlapping areas would amplify plane-fitting errors dozens of kilometres away but also because the discrepancies are smaller (the standard deviation of the difference fields is 0.3 m, Fig. S3). An exception is the two northernmost difference fields, where a trend ranging from −0.8 to 0.8 m remains. In addition to residual phase trends, discrepancies of ∼ 0.5 m can occur in areas where surface slope is locally elevated (e.g ice shelf channels or surface ridges). Altogether, we therefore estimate the SAR processing uncertainties to be of the order of 0.5 m.
Next, we compare the 2013 and 2014 DEMs with kinematic GNSS profiles from 2012 and 2014, respectively. The time lag between the satellite data acquisition and the collection of ground-truth data is thus 8-10 months for the 2013 DEMs, and 5-6 months for the 2014 DEMs. For 2012-2013, differences are −0.44 ± 1.05 m, and for 2014 −0.04 ± 0.65 m. The largest discrepancies occur in both datasets near ice shelf channels were ice advection within the multiple months time lag is significant (Fig. 9). Removing those areas reduces the discrepancies to −0.37 ± 0.29 m in 2012-2013, and −0.07 ± 0.2 m in 2014. Ignoring the dynamic influence of ice shelf channels, the highest discrepan- cies are found in the most upstream part of the 2014 GNSS profile (Fig. 9d). There, TanDEM-X elevations are systematically overestimated by up to 2 m near the grounding line. We attribute this bias to decreasing penetration of the TanDEM-X signal, as the firn air content decreases towards the grounding zone (Fig. 2a). The X-band radar signal can penetrate up to 8-10 m in cold dry snow (Humbert and Steinhage, 2011), and the bulk part of such a signal penetration would be accounted for during our offset correction. However, errors due to spatial variations in signal penetration remain but affect both the 2013 and 2014 DEMs. To conclude, we estimate that in most areas the relative accuracy of the TanDEM-X DEMs is in the sub-metre range. Errors are slightly elevated in areas where the local surface slope is high, and surface elevation is systematically and significantly overestimated by up to 2 m in a narrow belt close to the grounding line.

Hydrostatic inversion
The main uncertainties for the hydrostatic inversion are referencing the surface elevation to height above sea level, and accounting for density variations. The former depends on the geoid, the mean dynamic topography, tides, atmospheric pressure variations and eustatic sea level. Drews (2015) estimates errors in the geoid and the dynamic topography for RBIS to be within ±1 m. We account for tides and atmospheric pressure variations implicitly by offsetting the TanDEM-X DEMs to the CryoSat-2 DEM, which contains these corrections. The smallest component in the error budget are changes in eustatic sea level rise, which we neglect.
Variations in firn air content are important because these propagate with a factor of 9 into the hydrostatically inverted ice thickness (Eq. 2). We illustrate this point along profile PP where the inferred thickness from radar profiling and from phase-sensitive radar agree closely, but the hydrostatic thickness is > 80 m thinner (Fig. 8a). Because surface elevation is well constrained by our kinematic GNSS profiles (Fig. 8c), we attribute this large, unphysical mismatch to an overestimation of the firn air content. The firn densification model predicts a value of 11 m at that location. However, in the field it became evident that this area is close to a spatially extensive blue-ice area where firn air content is negligible. Reducing the firn air content to 1 m reconciles the hydrostatic ice thickness with the observed radar ice thickness (Fig. 8a). Such a large deviation of the modelled firn air content may be site specific because it is located in the transition zone where turbulent mixing by the katabatic winds and a windalbedo feedback form a micro-climate that causes extensive surface melting with not yet fully understood effects on the firn densification (Lenaerts et al., 2017). The impact of the firn air content misestimation on the derivation of the hydrostatic ice thickness is further discussed in Lenaerts et al. (2017). Moreover, Drews et al. (2016) used wide-angle radar measurements in conjunction with ice coring and found that firn density varies spatially over tens of kilometres scales, in particular across ice shelf channels, where surface melt water collects in the corresponding surface depressions and locally refreezes. Therefore, we anticipate that at least some of the variability seen in the LBMB field is due to unresolved variations in firn density.
Because of unaccounted variations in firn density, and uncertainties in referencing the freeboard height, our ice thickness field has a lower bound error of at least ±25 m (Drews, 2015). In some areas the error can be considerably larger. However, the corresponding impact on the inferred LBMB rates is mitigated by the low ice flow divergence rendering the magnitude of ice thickness less important (Eq. 1c).

Lagrangian matching
Computing the Lagrangian thickness change, requires matching the DEMs to account for ice advection. We use a normalized cross-correlation to match 5 km × 5 km patches from 2013 to the 2014 geometry (Sect. 2.5). Alternatively, the matching can be based on the surface flow field (Moholdt et al., 2014). For the DEMs, this methods yields similar results in terms of the large-scale LBMB pattern, but introduces erroneous positive/negative patterns near ice shelf channels. This is because the flow velocities are not sufficiently constrained for the flow direction, and tilts by a few degrees cause a significant mismatch in areas where thickness gradients are larger. On the other hand, the 2016 GNSS has to be matched with the velocities, because 2-D crosscorrelation fails with profiles.

Ice flow divergence: the benefits of regularized derivatives
The high-resolution velocity field is too noisy in magnitude to approximate the derivatives in the flow divergence with finite differencing of neighbouring cells (gridded to 125 m posting). This can be accounted for by smoothing the velocity field prior to taking the derivative. However, this type of smoothing can blur abrupt changes in the flow velocities and corresponding strain rates. This is important, because we suspect that ice flow velocities change abruptly in ice shelf channels that experience strong basal melting (Drews, 2015). We therefore explore the use of total-variation regularization, which treats abrupt (and discontinuous) changes more accurately (Chartrand, 2011). Figure 3 illustrates a close-up of an ice shelf channel (inset " Fig. 3" in Fig. 4) where we compare the "normal" (unsmoothed) velocity divergence (b) with its regularized (c) and smoothed (c-e) versions. For the latter, we applied average filters of 375 m × 375 m, 1125 m × 1125 m and 1875 m × 1875 m (i.e. kernels of 3×3, 9 × 9 and 15 × 15 pixels, respectively) to the velocity field, before computing the gradients. The enhanced velocity divergence has a similar magnitude in the regularized and the smoothed version using a 375 m× 375 m window. However, the latter is noisier outside the ice shelf channel than the regularized version. In the regularized case, velocity divergence at the channel's apex is 8, 24 and 40 % lower than for the 375 m × 375 m, 1125 m × 1125 m and 1875 m × 1875 m kernels, respectively. However, the inferred LBMB rates are insensitive to the technical implementation of the derivatives, because the Lagrangian thickness change controls the signal at RBIS. Nevertheless, in order to study the dynamics of the smaller-scale ice shelf channels, efficiently denoising the derivatives becomes increasingly important, in particular for ice shelves where the dynamic thinning terms is more important.

Surface mass balance
Both the firn air content and the SMB are spatially less well resolved than our ice thickness and velocity fields. Consequently, we do not capture their spatial (and temporal) varia-tions on the length scales associated with ice shelf channels. Both Drews et al. (2016) and Langley et al. (2014) found evidence in the shallow radar stratigraphy that the SMB may be locally elevated in those areas, potentially reflecting the deposition of drifting snow at the bottom of surface slopes (Frezzotti et al., 2007). If this holds true, then the systematic underestimation of the SMB would result in a positive bias of the LBMB in those areas.

Discussion
The large-scale patches of enhanced basal melting (Sect. 3.1; labels A-C in Fig. 4) are sufficiently far away from the tidal bending zone so that we can safely assume hydrostatic equilibrium. These regions are also detected by Rignot et al. (2013), based on different input datasets (i.e. Eulerian thickness change based on ICESat-1). Patches A-C line up with deepest parts of the ice shelf base and the largest gradients in the hydrostatic ice thickness. A large ice draft fosters basal melting because the freezing point is lower with depth (e.g. Holland et al., 2008). The steep basal slopes facilitate entrainment of heat in the mixed layer beneath the ice shelf increasing basal melting (Jenkins and Doake, 1991;Little et al., 2009). The smaller-scale variations in LBMB are more difficult to interpret, because these are overlain by unaccounted variations in firn density, SMB, and ice that is not in hydrostatic equilibrium. Nevertheless, the comparison with the phasesensitive radar data and the kinematic GNSS profiling increases our confidence that much of the relative variability that we observe here is meaningful. The surface lowering of the elliptical surface depression is consistently observed over a 3-year time period marking this zone as dynamically active. Two other options are (i) a transient adjustment of the surface towards hydrostatic equilibrium (Humbert et al., 2015) as a response to some unknown event in the past which locally reduced the ice thickness, and (ii) the surface lowering may reflect vertical creeping of a liquid water body through the ice column. In any case, the surface lowering is restricted to a small area and the ice shelf channel farther downstream appears passive (i.e. does not show significant melting nor refreezing).
In most areas, ice shelf channels at RBIS seem to advect passively and basal melt rates there do not significantly stand out from those in the larger surrounding. Exceptions are the locally elevated basal melt rates in ice shelf channels in the interior of the RBIS (e.g. inset " Fig. 3" in Fig. 4) and close to the grounding zone ( Fig. 6 and its corresponding inset in Fig. 4). Almost all ice shelf channels at RBIS are connected to the grounding line and may arise from water-filled subglacial conduits injecting subglacial melt water into the ice shelf cavity, driving a spatially localized buoyant melt water plume (Jenkins, 2011;Le Brocq et al., 2013;Drews et al., 2017;Sergienko, 2013). Such localized melting near the grounding zone has been previously observed on Pine Island Ice Shelf using similar methods as done here (Dutrieux et al., 2013). However, on Pine Island Ice Shelf, background melt rates are an order of magnitude larger than what is observed here (Depoorter et al., 2013a;Rignot et al., 2013) and Dutrieux et al. (2013) analysed DEMs separated by 3 years (compared to the 1-year time period used here). This explains why locally elevated BMB values appear more clearly on other ice shelves. We find some evidence that basal melting is concentrated on the flanks, rather than on the apex (Fig. 6). This accords both with observations (Dutrieux et al., 2014) and modelling (Millgate et al., 2013). Dutrieux et al. (2014) suggest that the presence of a colder water blocks the heat flux from below near the apex of the channel. Alternatively, modelling suggests (Millgate et al., 2013) that a geostrophic current develops beneath the channels (if the channels are wide enough) which preferentially melts at the channel's flanks. This seems less likely here because ice shelf channels near the grounding line are narrow (i.e. a few hundred metres wide and high).
In summary, our observations suggest that the LBMB varies on multiple spatial scales which has several implications. First, point measurements with phase-sensitive radars are not necessarily representative for a larger area. Particularly in areas where thickness gradients are large, phasesensitive radar measurements are best understood in combination with satellite-based estimates covering larger spatial scales. On the other hand, on-site point measurements are crucial to estimate the quality of the satellite-based BMB estimates, which are uncertain in their magnitude. Second, this sub-kilometre variability in ice-ocean processes poses challenges for coupling ice flow with ocean models, because highly resolved ocean models and community efforts, such as the Marine Ice Sheet-Ocean Model Intercomparison Project (MISOMIP), are typically gridded with 1-2 km (Dinniman et al., 2016;Asay-Davis et al., 2016). This is too coarse to capture the spatial variability that we observe here.

Conclusions
We derived the Lagrangian basal mass balance (LBMB) of the Roi Baudouin Ice Shelf by combining TanDEM-X DEMs of 2013 and 2014 with high-resolution surface velocities and atmospheric modelling outputs. On a large scale, the LBMB shows the highest basal melt rates where the ice draft is deepest and steepest, i.e. close to the grounding line and near Derwael Ice Rise and the Western Ice Promontory. This pattern is overlain with significant sub-kilometre scale variability, as witnessed by localized surface lowering of an elliptical surface depression and large basal melting rates below some sections of ice shelf channels. For the latter, we find evidence that at least in some areas, basal melting is concentrated on the channel's flanks as opposed to its apex. Key advancements in our methodology to eluci-date this variability are (i) the calibration of the DEMs to account for residual trends from the interferometric processing, (ii) the quality of the matching procedure -using normalized cross-correlation coefficients -for calculating the Lagrangian thickness change, and (iii) the total-variation regularization of the spatial derivatives that preserves abrupt changes in flow velocities that are sometimes observed across ice shelf channels. New satellites (such as TanDEM-X or Sentinel 1) will continue to provide highly resolved datasets of surface elevation and ice velocity. In comparison, atmospheric modelling does not (yet) provide the required spatial resolution on firn density and SMB to solve the mass budget reliably on sub-kilometre scales. Although the uncertainty of the absolute LBMB values remains high, we find a good fit with on-site measurements from phase-sensitive radar, and we demonstrate that much of the spatial LBMB variability contains information about ice shelf processes occurring at sub-kilometre scales. This variability highlights the complexity of the ice-ocean and ice-atmosphere interactions on small spatial scales on ice shelves, which need to be accounted for by glaciologists, oceanographers and atmospheric scientists.
The Supplement related to this article is available online at https://doi.org/10.5194/tc-11-2675-2017-supplement. excellent logistic support by the Belgian Military, AntarctiQ and the International Polar Foundation during the field campaigns. Finally, we thank Jan Lenaerts and Stefan Ligtenberg for sharing results from atmospheric modelling, as well as David Shean and Geir Moholdt for their constructive comments on this paper.
Edited by: Andreas Vieli Reviewed by: Geir Moholdt and David Shean