The regional scale surface mass balance of Pine Island Glacier, West Antarctica over the period 2005–2014, derived from airborne radar soundings and neutron probe measurements

. We derive recent surface mass balance (SMB) estimates from airborne radar observations along the iSTAR traverse (2013,2014) at Pine Island Glacier (PIG), West Antarctica. Ground based neutron probe measurements of snow density at 22 locations allow us to derive SMB from the annual internal radar reﬂection layers. The 2005 layer was traced for a total distance of 2367 km allowing us to determine annual mean SMB for the period 2005–2014. Using complementary SMB estimates from the RACMO2.3p2 regional climate model and a geostatistical kriging scheme we determine a regional scale 5 SMB distribution with the same main characteristics as that determined for the period 1985–2009 in previous studies. Local departures exist for the northern PIG slopes, where the orographic precipitation shadow effect appears to be more pronounced in our observations, and the southward interior, where the SMB gradient is more pronounced in previous studies. For the PIG basin, we derive a total mass input of 79 . 9 ± 19 . 2 Gt yr − 1 . This is not signiﬁcantly different to the value of 78 . 3 ± 6 . 8 Gt yr − 1 for the period 1985–2009. Thus there is no evidence of a secular trend in mass input to the PIG basin. We note, however, 10 that our estimated uncertainty is more than twice the uncertainty for the 1985–2009 estimates. Our error analysis indicates that uncertainty estimates on total mass input are highly sensitive to the selected krige methodology and assumptions made on the interpolation error, which we identify as the main cause for the increased uncertainty range compared to the 1985–2009 estimates. we lack measurements from the southward interior. Earlier observations from ME14 suggest


15
The stability of the West Antarctic Ice Sheet (WAIS) is a major concern for scientists seeking to predict global sea level rise.
Transport of heat from the upwelling circumpolar deep water has proved to be a critical driver of Antarctic ice shelf thinning and grounding line retreat, thus stimulating the acceleration of marine-terminating outlet glaciers (e.g. Hillenbrand et al., 2017).
In particular the Amundsen Sea sector has experienced an unprecedented acceleration in ice discharge since the beginning of satellite based ice flow observations in the 1970s. Three quarters of this ice discharge stem from the Thwaites and Pine Island 20 glaciers with both showing evidence of rapid acceleration since the 1970s (Mouginot et al., 2014) and spreading of surface lowering along their tributaries over the past two decades (Konrad et al., 2017). While spaceborne observations indicate that this acceleration has levelled off recently (Rignot et al., 2019), they also support model projections suggesting modest changes  Table 1. Approximate sample bin resolution and maximum depth of: SAR level_ 1b processed ASIRAS data with indicated vertical range bin standard deviation based on our two-way-travel time (TWT) to depth conversion, CreSIS Accumulation radar according to Medley et al. (2014), and pulseEKKO PRO GPR discussed in Konrad et al. (2019). For the GPR system we estimate the maximum sampling depths based on shared radargrams, which resolve the internal stratigraphy at PIG for TWTs up to (1000-1200) ns.
that the SMB decreases towards the interior so the contribution from this area to the total mass input will be less than that from 60 the rest of the basin.
Additional SMB measurements were made with a ground penetrating radar during traverse T1 and published in Konrad et al. (2019). They selected the ∼ 1986 reflection layer, which approximately coincides with the observed main reflector by ME14, and traced the layer along sections of the 900 km traverse, amounting to a total of a 613 km distance covered by GPR observations. Due to the limited maximum sampling depth of the ASIRAS and NP measurements, the 1985/86 reflection layer used by ME14 and Konrad et al. (2019) is not contained in most of our data. To benefit from the ASIRAS coverage while simultaneously accounting for its limited depth range, we manually traced the continuous 2005 reflection layer over a distance of 2367 km to derive mean annual SMB estimates for the 2005-2014 period. Due to the reported consistency between the GPR and airborne SMB measurements in Konrad et al. (2019), we limit the comparison of our results to the basin wide estimates by ME14. In addition, we assume that the effect of strain history, which could affect our SMB estimates at the fast flowing sections 70 of PIG, is negligible. Konrad et al. (2019) conclude that the total effect over the whole catchment is small, even though it can have a very significant effect at some sites. However, this effect is expected to be further reduced for the shallower reflection layer depths from the ASIRAS measurements.

Neutron Probe measurements
NP measurements of snow density were performed at all stations during both traverses as described in Morris et al. (2017). To 75 evaluate the effect of densification, they repeated the density profiling in the same boreholes during traverse T2. Because the most recent accumulation is missing in these profiles, they drilled an additional borehole of less than 6 m depth and a nearby distance of about 1 m to capture it during traverse T2. The only exception is site 2, where the ground team decided to auger a completely new 14 m borehole for the density profiling due to poor data from the T1 hole.
The deep firn cores (∼ 50 m) denoted in Fig. 1   modulating processes of density and H 2 O 2 allow for an independent determination of annual snow accumulation at the 10 deep core sites. Morris et al. (2017) applied an automatic annual layer identification routine to their snow density profiles and used 85 the annual H 2 O 2 peak depths as an additional guidance for the annual layer dating. Thus, the depth-age scales from both annual markers are consistent.
We use a single regional density-depth profile for the two-way-travel time (TWT) to depth conversion of the ASIRAS soundings, which we derive from the 43 NP profiles of traverse T2. First, we merge the ∼ 13 m and nearby ∼ 6 m density-depth profiles from the NP measurements during traverse T2 by linearly relaxing their overlapping segments. To reduce the effect 90 of lateral noise convolution, we limit the relaxation length to the overlapping segments that correlate well with each other.
Then we align the intercepting depth-age scales to create a consistent depth-age scale for each compiled profile. The resulting merged profiles are shown by the grey lines in Fig. 2. From these profiles, we then determine a smoothed regional mean profile, which is denoted by the black line. Finally we fit an exponential function to the regional mean profile (red dashed line), which we apply to the TWT to depth conversion. The blue dashed lines show the fitted standard deviation of the density as  a function of depth. Following Medley et al. (2013), we consider the fitted standard deviation to be representative for the spatial uncertainty of the regional scale density-depth profile.

ASIRAS soundings
ASIRAS is a Ku-band radar altimeter which operates at a carrier frequency of 13.5 Ghz and a bandwidth of 1 GHz (Mavrocordatos et al., 2004). As in a previous campaign in Greenland (Overly et al., 2016), it was operated in Low Altitude Mode at 100 PIG (i.e. less than 1500 m above ground). A Synthetic Aperture Radar (SAR)-processing of the collected data was performed, which yields the spatial resolution of the SAR level_1b data shown in Tab. 1. The associated cross track footprint is ∼ 15 m.
We use the electromagnetic wave speed v = c/ √ to convert TWT to depth, where c is the vacuum speed of light and is the real part of the dielectric permittivity of the firn column. For the latter, we apply the commonly used empirical relation by Kovacs et al. (1995): where ρ s = ρ/ρ w is the specific gravity of snow at current depth with regard to the water density ρ w = 1000 kg m −3 . An alternative model by Looyenga (1965) is with ice = 3.17 (Evans, 1965) and ρ ice = 917 kg m −3 . Sinisalo et al. (2013), who consider a similar depth range to this study, conclude that the difference between wave speeds based on Eq.
(1) and (2) has a negligible impact on their SMB estimates.
This is also the case for our estimates (see Section 3). The maximum depth of the radargrams is ∼ 30 m based on the TWT to depth conversion from the fitted regional mean profile of density with depth and substituted Kovacs relation. The depth range of resolved internal stratigraphy varies along the flight track, but the layering remains visible for most of the upper 13 m depth covered by the NP measurements. Using the mean depth density profile, we determine the w.e. depth value for each 115 waveform bin and calculate the mass per unit area between the selected reflection layer (magenta line in Fig. 3) and surface.
We assume that the annual peaks in density lead to reflection layers because of the associated dielectric constant (Eisen et al., 2003). However, we also find that the ASIRAS soundings appear to resolve sub-annual layers at locations where the snow accumulation is high, e.g. around iSTAR site 21 as shown in Fig. 3. These layers may be created by intra-annual events, e.g. hoar events that form thin ice layers (Arcone et al., 2004), which generate a noticable dielectric contrast with regard to 120 the ASIRAS soundings. Hence, we rely on the independent measurements at the iSTAR sites to distinguish intra-annual from annual layers. Before the layer tracing, we apply an automatic-gain-control filter to all waveforms and limit their dynamic range to twice the standard deviation centred around the mean amplitude of each waveform. This improved the signal contrast of the radargram. Initially we tested a phase following algorithm of the Paradigm EPOS geophysical processing software to trace semi-automatically the selected reflection layer. However, this method became unstable for lower contrast and cases with close 125 layer spacing. Furthermore, remaining SAR-processing artefacts were interfering with the phase following algorithm. Because of the complex nature of the observed stratigraphy, as has been also reported by Konrad et al. (2019), a manual layer tracing was used. Following Richardson et al. (1997), we attempted to bridge distorted or merged layer segments whenever distinct characteristics of a vertical layer sequence could be identified with confidence before and after the bridging. Different processes can lead to distortion of the reflection layers, e.g. processes changing deposition of the annual snow layers or excessive rolling 130 angles of the airplane, while merging layers can result from low snow precipitation and ablative processes (e.g. wind-scouring) or a combination of both. We compared the traced layer at 34 intersections and 8 nearby flight track segments. In case of layer mismatches, the layer tracing was iteratively adjusted to find a consistent layer for the entire flight track. In this sense, the layer tracing is performed independently from the annual layer dating at each iSTAR site.

Measurement error estimation 135
We attempt to trace the 2005 reflection layer, which is covered by all NP depth-age scales. Following Morris et al. (2017), we define mass balance years between the density peaks in the NP profiles (nominally 1st of July). For instance, the mass balance year 2013 begins at the second annual density peak below the surface (nominally 1st July 2013 and ends at the first peak (1st July 2014). At N points along the closest approach of the ASIRAS track to each iSTAR site we estimated the depth of the reflection layer from its TWT, using the fitted density profile from that site. We also tested the use of the measured density 140 profile for the TWT conversion instead, but we find that the impact on our results is negligible similar to Morris et al. (2017).
The estimated depths do not necessarily coincide with the depth of a peak in the density profile i.e with the start of a mass balance year. We assign an estimated date to the mean value of N points by interpolation using the depth-age scale from the  Table 2. Dating (Year) with associated uncertainty of closest ASIRAS reflection layer to i th iSTAR site. "Track number" refers to the ASIRAS flight track naming convention and "∆D" is the closest distance between the ASIRAS track and the iSTAR site. Years in brackets are discarded from the regional layer age estimation for the reasons given in the Comment column. uncertainty as the standard deviation σ x in the N measurements. In this sense, our error estimate is more conservative than the standard error of the mean. Furthermore, we assume that the uncertainty due to local variation in the stratigraphy is isotropic, which does not generally need to be true. However, according to Tab. 2 in most cases the overall impact of this effect is one 150 order of magnitude smaller than the variability of dated years among all iSTAR sites. As indicated in Tab. 2, we excluded dating estimates around iSTAR sites 2 and 19. In both cases, our layer tracing revealed a large offset contrary to the neighbouring iSTAR sites. Possible reasons for these offsets could be systematic errors in the layer dating from the NP profiles, the variability of internal stratigraphy between the ASIRAS measurements and their closest approach to both iSTAR sites or systematic errors in the manual reflection layer tracing. The remaining exclusion of dating estimates at iSTAR sites 7, 12, and 18 is either due to 155 high noise levels of the radargram or reflection layer depths exceeding the NP depth-age scales. Following Konrad et al. (2019) we estimate the reflection layer year by the mean of dating values at each site with an uncertainty of ∆t = δt 2 + δt 2 + t 2 x , with the annual layer tracing uncertainty δt = ±1 years, the standard deviation of dating estimates δt, and the propagated error t x = 1/n n i σ x (i) 2 from the n lateral error estimates around each iSTAR site (i = site number), which we introduced in addition. The resulting reflection layer dating estimate is T = 2004.8±1.4, which corresponds to a layer age of a = 10.1±1.4.

160
The associated average surface accumulation rateḃ in terms of w.e. depth per year iṡ where δz i is the i th depth increment of the radar waveform and ρ i is the associated density. Substitution of the wave propagation

165
where t s = 0.37 ns is the ASIRAS vertical bin sampling time (i.e. 0.5 × TWT per bin), and i refers to the permittivity value at the i th bin. To avoid any confusion with previous summations, the final index m refers to the traced waveform bin at the reflection layer depth. It is evident from Eq. (4) that the spatial uncertainty of the density profile affects both the integration depth and incremental mass. Medley et al. (2013) and ME14 estimated the spatial uncertainty from the resulting SMB change by directly applying the standard deviation fits of their regional density profiles to the TWT to SMB conversion. Instead, we 170 may propagate the error in Eq. (4), assuming that errors are uncorrelated and normally distributed. Based on the Kovacs relation according to Eq. (1) we account for the temporal, spatial, and digitization: where ∆a = ±1.4 years is the temporal uncertainty and ∆ρ i are the standard deviation intervals according to Fig. 2. Due to the small incremental density change of < 0.7% along the entire profile, we approximate the digitization error by the mean SMB 175 value of three consecutive bins centred at the final profile bin of the current integration depth. Figure 4 (a) displays the integrated individual measurement error components as well as the combined error according to Eq. (5) as a function of geometric depth.
In addition, we include the error partitioning in (b). The grey background shades highlight the distribution of layer depths to visualise the relevant error range of our SMB estimates, which peaks around (5, 8, and 10) m. By comparison with Medley et al. (2013) and ME14, we find that our spatial error estimate based on Eq. (5) is reduced by about one order of magnitude 180 while the standard deviation fits of their regional density profiles cover a similar range compared to ours. We may ignore the spatial error compensation in Eq. (5) by replacing the root-sum-of-squares (RSS) with absolute values: Hence, to comply with the studies above, we consider the more conservative spatial error propagation based on the sum of absolute values, but we keep the RSS of individual error components for the combined measurement error estimate as shown in Fig. 4 (c-d). Following these assumptions, we find that our measurement error estimate is still dominated 185 by the temporal layer dating uncertainty for most of the traced layer depths, but the spatial error reaches a similar range as reported in Medley et al. (2013).

Kriging scheme
We focus on the regional scale variability of the SMB distribution on PIG.    window length. We initially tested the same interpolation scheme as described in ME14 to estimate a regional scale SMB field for the PIG basin from our smoothed SMB points. This scheme is based on the ordinary kriging (OK) algorithm, a widely used geostatistical interpolation technique (e.g. Isaaks and Srivastava, 1991). Instead of a direct OK interpolation of smoothed SMB observations, ME14 consider the residual SMB values with regard to an ordinary least squares linear regression model for the 195 Thwaites-PIG basin area with northing, easting, and elevation as explanatory variables. This, in turn, yields a small degree of skewness < 0.5 of their new SMB distribution. However, we failed to reduce the skewness of residual SMB values effectively by the same method, which may by due to the different aerial coverage considered in our regression model. Examination of the DEM contour lines in Fig. 5 reveals that a simple relation between surface elevation and SMB is not evident, which may hint that the prevailing synoptic scale weather conditions at the Amundsen and Bellingshausen Sea sectors in combination with the 200 precipitation shadowing effect of the Eights Coast mountain range (Fig. 1) require a more sophisticated model to capture the SMB at the PIG basin scale. An alternative approach, which is also mentioned in ME14, is a logarithmic transformation of the SMB observations prior to the OK interpolation: where C is an arbitrary constant and x 0 represents the current interpolation location. After the OK interpolation of transformed 205 SMB observations, the estimates must be transformed back into the original measurement scale. This backtransformation requires the addition of a nonbias term for each OK estimate to ensure that the expected value is equal to the sample mean and that the smoothing effect is adequately compensated (i.e. resulting estimates reproduce the sample histogram and sample mean [Yamamoto, 2007]). We implemented such ordinary logarithmic kriging (OLK) method in our analysis by adopting the 4-step post-processing algorithm proposed by Yamamoto (2007) for the estimation of nonbias terms. According to Yamamoto (2008),

210
OLK does not necessarily require a log-normal sample distribution to produce improved estimates in terms of local accuracy.
Furthermore, Yamamoto (2007) tested the impact of constant C according to Eq. (6) and found that a data translation towards higher values yields an approximation from OLK to OK estimates, thus, eliminating the advantage of improved sample mean reproduction and local accuracy of OLK estimates. Indeed, we find that adding a negative constant C to all SMB values, such that the lowest SMB value reaches 0.1 kg m −2 yr −1 , yields an improved reproduction of the observation data characteristics.
215 Figure 6 shows the experimental isotropic semivariogram of our log-transformed SMB observations from Fig. 5 (b) together with a Gaussian model fit with a practical range of ∼ 190 km, i.e. the range at which the spatial autocorrelation of sample points is vanishing. Following Yamamoto (2005Yamamoto ( , 2007, we investigate the reproduction of observational data characteristics by means of PP-plots (i.e. percentiles of cumulative distributions of observations and estimates against each other). Figure 7 shows the PP-plots for our OLK and OK interpolation constrained to a maximum estimation range with regard to the closest ASIRAS However, after extending the estimation range, it is evident from Fig. 7 that the best match exists between the observation and OLK estimation values. Hence, we limit our analysis to these values in the following.  (1996) to correct for negative kriging weights (Yamamoto, 2000) and constrain all processing steps of the OLK estimation to the 16 nearest neighbours for each estimate according to the quadrant criterion. Depending on the considered neighbourhood 230 the effect of smoothing as well as local stationarity of observation data is affected. As a guidance for our final setting, we aimed at generating an optimal PP-relation according to Fig. 7, but also considered potential artefacts, which may arise from the OLK procedure.
In addition to each OLK estimate, we calculate the associated interpolation error. While ME14 choose the kriging standard deviation as a measure of interpolation error, our error estimation is based on the introduced interpolation standard deviation 235 S 0 by Yamamoto (2000) for two reasons. Firstly, as shown by the author, S 0 represents a more complete measure of local accuracy and has, therefore, been implemented in the post-processing algorithm in Yamamoto (2007). Secondly, for the OLK method we need a corresponding backtransformation of the interpolation error from the logarithmic to the measurement scale, which has already been investigated for S 0 in Yamamoto (2008). Thus, we adopted the proposed backtransformation of S 0 in this study.

240
Following ME14, we estimate the total error of each SMB estimate by the RSS of the measurement error and transform S 0 back.
The measurement error is estimated by generating 500 realisations of OLK SMB estimates with add noise to the smoothed SMB observations, which follows a normal distribution with a mean of zero and standard deviation equal to the measurement error of the SMB observation at x 0 .
For the basin wide SMB estimation we have to keep in mind that our OLK estimation is limited in terms of the practical 245 range according to Fig. 6. By comparison with the flight track shown in Fig. 1, even when considering the practical range as a maximum threshold for the spatial SMB estimation, we do not cover the entire PIG basin. Hence, for the calculation of total mass input for the PIG basin, we replace our SMB estimates with modelled SMB from the RACMO2.3p2 (van Wessem et al., 2018) regional climate model(in the following abbreviated as RACMO) at distances where the spatial autocorrelation of measurements is low.

Regional scale SMB distribution
Based on the adopted OLK interpolation scheme, we produced the mean annual SMB map for the PIG basin from the ASIRAS observations in Fig. 8 (a). SMB observations and estimates are colour coded with the same scale. Each estimate covers a pixel By comparison with ME14 and the extracted RACMO estimates in Fig. 8 (  the inland, and a region of low SMB in response to the shadowing effect from the Eights Coast mountain range. The RACMO estimates indicate a sharp boundary of this effect at the ice divide at the northern tip of PIG. This boundary is also captured by the ASIRAS observations at its northernmost flight track near iSTAR site 10, which is best seen in the high resolution observations according to Fig. 5 (a). the RACMO estimates. Despite the missing ASIRAS observation towards the southern interior, the latter is already apparent inside the 100 km maximum distance margin (inner dashed line). By comparison with the SMB distribution according to ME14, we find similar departures to our results with respect to the northern slopes of PIG and its southward interior.
The ME14 results indicate an elevation dependent drift between observed and simulated SMB values, which we also find 270 between our results and RACMO as shown in Fig. 9. Here, we first calculated the average SMB for each RACMO pixel from collocated ASIRAS estimates. Then we subtracted the RACMO SMB estimates from the ASIRAS averages and plot the resulting difference against their associated elevation, which we determined from the DEM model in Fig. 5. Due to the limited practical range of the ASIRAS observations, we only consider ASIRAS-RACMO estimates within the 100 km maximum range limit to the observations inside the PIG outlines. Like ME14, we find for this version of RACMO that simulated SMB estimates 275 appear to be lower than the ASIRAS estimates at higher elevation levels and vice versa. Despite the scatter between elevation and SMB difference, a Kendall's rank correlation test with α = 0.05 significance level suggests that both are significantly correlated with each other.

Total mass input
For the calculation of spatially integrated SMB, denoted by Σ + in the following, we produced the hybrid SMB map in Fig. 8 (d), 280 where the ASIRAS SMB estimates are linearly relaxed into the RACMO SMB estimates between the 100 km and 190 km range limits to account for the spatial autocorrelation in terms of Fig.6. Table 3 summarizes Σ + and further statistical SMB characteristics for different data sets and basin definitions according to Fig. 10 (d). Here, we replaced the interpolation artefact highlighted in Fig. 8 (a) with averaged values from neighbouring pixels. Σ + uncertainty estimates refer to the RSS of the  (408) 151 (157) 147 (144) 935 (989) PIG (Mouginot et al., 2017) 70.6 (72.1) 400 (409) 150 (157) 132 (131) 935 (989) PIG (Zwally et al., 2012) 81.6 (83.6) 392 (401) 151 (158)  interpolation and measurement error-grids ( Fig. 10 d) in accordance with ME14. Because RACMO is missing an error-grid, 285 we consider the total combined error for the entire PIG basin as a conservative error estimate. In this sense, we are augmenting the missing model error estimation. To quantify the relative contribution of ASIRAS to the hybrid SMB estimates, OLK area and OLK Σ + denote the relative contribution in terms of covered land area and integrated SMB respectively. For comparison with our hybrid based estimates, we include results from RACMO and ME14, which we converted from w.e. depth to SI units.
Because of the different averaging periods between this study and ME14, we add RACMO estimates in brackets, which we 290 extracted based on the same averaging period as the ME14 results.

Pine Island and Wedge Zone
The Pine Island Σ + values are in agreement between all data sets within the estimated error margins. This is different for the definitions. The impact on the mean annual SMB is limited to 10 kg m −2 yr −1 (∼ 2%) between all catchment definitions.

Discussion
We discuss first the pronounced differences between annual layer dating from ASIRAS reflection and neutron probe density profiles at some sites and then secondly the systematic differences in SMB distribution between our results and those of ME14 and RACMO.

SMB departures
Key to the evaluation of our selected internal reflection layer is its isochronic nature, which we assume based on matched depth-age relations from the iSTAR ground truthing measurements. One may argue that these measurements can be subject to local noise in the density profile, which would challenge any comparison with nearby radar observations. For instance, Laepple et al. (2016) observed dominating stratigraphic noise at single pit density profiles near Kohnen station (East Antarctic plateau,

315
Dronning Maud Land). Stacking of multiple profiles is one possibility to filter out noise. While this is not possible for the single iSTAR sites, the estimated dating uncertainty of ±1.4 years according to this study suggests that iSTAR ground truthing measurements at PIG are less prone to stratigraphic noise, which is most likely to be related to the higher SMB compared to ∼ 70 kg m −2 yr −1 nearby Kohnen station (Laepple et al., 2016). However, on a few occasions we identified larger departures in the annual layer dating, as it is the case for iSTAR site 2 (Tab. 2). While the layer tracing appears to be in agreement  Site 19 is directly located at the centre of a pronounced accumulation trough of ∼ 2.5 km width, which adds to the uncertainty in the layer matching because of the spatial displacement between the iSTAR site and point of closest approach. Because the traced reflection layer exceeds the depth range of the NP density profiles at site 18 and 19, we discarded both sites for the layer dating.

335
Additional local departures between our results and those from ME14 were identified for the northern slopes and southward interior of PIG. Because of difficulties in the layer tracing at the northern slopes, they had to augment their SMB estimates with results from a diferent layer, which they dated back to 2002 and corrected for a temporal bias to the 1985 layer based on overlapping segments. Thus, one possible explanation for the observed differences is that the true local temporal bias correction may be different from the regional scale bias correction, which they estimate from regression models. Other possible 340 explanations are differences in the observational coverage and local accuracy from the different interpolation methods. With regard to the southward interior, the spatial coverage is superior for the ME14 results. Despite our limit of a maximum range between 100 km and 190 km for the ASIRAS based estimates, the missing observational constraints towards the interior may still yield an underestimation of the southward SMB gradient. However, we also cannot rule out that the smaller gradient in our observations is due to a local increase in SMB between the different observational periods among both studies. Additional 345 observational constraints of the selected reflection layer may resolve the cause on the observed difference.

Impact on recent mass balance estimates
Despite the local differences in the SMB distribution, the difference between the Σ + estimates for the PIG catchment (including Wedge) between this study and ME14 is small, i.e. our Σ + is greater by 1.7 Gt yr −1 , which corresponds to 2% of the ME14 value. This indicates that the local differences in the SMB estimates between both studies cancel out. If we take 350 into account that the temporal averaging time used by ME14 is about a factor of 2.7 larger than that used in this study, this indicates that the observed total mass input for the PIG basin is not subject to any secular trend throughout the observational period, a finding which is further corroborated by the iSTAR measurements and RACMO simulations. This provides additional evidence to Medley et al. (2013) that the recent temporal evolution of the PIG mass balance is primarily driven by dynamic ice loss into the Amundsen Sea.

355
With regard to existing mass balance estimates for PIG, we have to take into account that the definition of basin can differ significantly as illustrated in Fig. 10 (d). To evaluate the impact of our hybrid SMB estimates on recent mass balance inventories, we extracted results from the literature in Tab. 4 and added updated mass balance estimates Σ − + by replacing the Σ + estimates from the literature with the Σ + estimates of this study. For the periods shown, we assume that the SMB remains stationary for the mass balance calculation. In addition, we linearly interpolated the estimated ice discharge measurements in ME14 for the 360 missing periods before 2007. Furthermore, we assume that the unspecified basin definitions in ME14 are in close agreement with the basin definitions based on Fretwell et al. (2013).
The small difference between the Σ + estimates of this study and ME14 directly translates into the Σ − + mass balance estimates. The largest impact of our results is on the Σ − + estimate by Gardner et al. (2018). After replacing their Σ + estimate from RACMO2.3 simulations with our hybrid Σ + estimate, the mass balance increases by ∼ 9 Gt yr −1 .

SMB uncertainty
While the agreement in Σ + estimates between this study and ME14 supports the hypothesis that the regional SMB of PIG is stationary at decadal scales, our uncertainty estimates are much larger. The temporal error according to Fig. 4, which is ∼ 5% larger than Medley et al. (2013) and ME14, cannot fully explain the difference between both uncertainty estimates. We also do not expect any major differences with regard to the spatial uncertainty of the density profiles. According to the error-grid 370 statistics in Tab. 5, we identify the backtransformed interpolation standard deviation S 0 from the OLK scheme as the dominating error source of our results, while the combined error in ME14 is slightly above our measurement error. The dominating S 0 uncertainty is also evident in Fig. 10 (a,b,c), where the spatial features of the combined error-grid are predominately determined by the S 0 grid. We find that the low accumulation zone at the northern slopes of PIG, which is next to the main trunk between iSTAR site 1 to 6, shows combined S 0 patches that considerably exceed 100%. In contrast, combined error estimates in ME14 375 do not exceed 20% at the same location.
Initial tests on our OLK setting revealed that the choice of the negative kriging weights correction method has a noticeable impact on the uncertainty estimates, a finding, which according to our knowledge, has not been reported before. However, our applied method by Deutsch (1996) already yields the minimum uncertainty estimates for our results, whereas the additional methods cited in Yamamoto (2000) yields an additional uncertainty increase between 20% (Froidevaux, 1993) and 50% (Jour-380 nel and Rao, 1996).
Additional tests, where we used the kriging standard deviation based on non-transformed OK estimates, did not improve our interpolation uncertainty. Therefore the different choice of the interpolation uncertainty measure is not the source of the larger uncertainty range of this study. We hypothesize that despite the homoscedastic (i.e. data value independent) nature of the krige standard deviation, the reduction of data variance after subtracting the regression surface according to ME14 is most likely the cause of their significantly lower uncertainty estimates.
In addition to the larger uncertainty range of this study, we note that the choice between cell-by-cell summation and RSS of grid errors has a quite substantial impact on the Σ + uncertainty estimates. If we make the optimistic assumption that gridded errors are independent and choose the calculation of RSS instead, Σ + uncertainty estimates would reduce to ±0.5 Gt yr −1 (i.e. ∼ 97% less) for the combined Pine Island and Wedge basin.

Systematic retrieval impacts
In addition to the uncertainty assessment in section 4.3, we evaluated the impact of artificial cluster removal, the choice of the permittivity model, and the non-transformed OK scheme.

Artifical cluster removal
Inspection of the artificial cluster highlighted in Fig. 5 revealed that it is centred around the location with the lowest observed 395 SMB and is essentially generated by the local nonbias terms of the OLK procedure. Owing to its steep contrast with the surroundings, it appears to be plausible to replace this cluster by averaged values of its nearest neighbours. However, due to the limited extend of this cluster, its additional contribution to the Σ + estimates would be less than 0.8%. Similarly, the impact on the PP-plot is negligible. Increasing the translational constant C helps removing this cluster, but at the cost of statistical agreement between observations and estimates. 400

Looyenga based results
Defining by Eq. (2) instead of Eq. (1) yields a minor reduction of Σ + for the PIG catchment of 0.6%, which we expect from Sinisalo et al. (2013). However, despite the minor impact of the alternative definition for , we noticed an additional small impact on the layer dating, which shifted our estimated layer formation from November to September 2004. Thus, we had to adjust the time range in the RACMO SMB extraction for the calculation of hybrid SMB estimates. While the choice of 405 the model only has a minor impact on our total mass input estimates, it is worth noting that the effect on our annual layer dating is detectable.

Non-transformed kriging results
If we choose the OK procedure instead, Σ + increases by 4% for the Pine Island and 12% for the Wedge zone outlines, which would further increase the offset between this study and ME14. However, inspection of the SMB distribution (not shown) 410 indicates that estimates tend to overshoot near the coastline of the Amundsen Sea, which becomes particularly evident for the Wedge outlines. Hence, the OK procedure appears to be more sensitive to the limited observational constraints near the Wedge outlines. In addition, S 0 based uncertainty estimates increase by 27% and 88% for the Pine Island and Wedge outlines, which highlights the improved performance of the OLK procedure.
Our analysis provides updated mean annual SMB estimates for the PIG basin and 2005-2014 averaging period based on a comprehensive airborne radar and ground truthing survey and complementary model simulations. Based on these estimates, we calculated a total mass input of 79.9 ± 19.2 Gt yr −1 for the PIG basin area. In comparison with earlier estimates from airborne radar observations, which consider the 1985-2009 averaging period, our result shows 2% greater total mass input.
Given the uncertainty in both values, there is no significant difference between them. Hence, no distinct secular trend is visible 420 between both averaging periods. We conclude that our results provide further evidence that the recent total mass input can be considered stationary at decadal scales. This implies that the increased dynamic ice loss over past decades remains the driving force in the recent mass balance evolution of PIG. However, departures between both observations at the northern slopes and southward interior of PIG may indicate temporal changes in the local SMB distribution, which cancel out for the estimates on total mass input. Furthermore, our radar based observations can resolve a discrepancy between strainrate and SMB 425 measurement at iSTAR site 2, which highlights the benefit of such complementary SMB measurements for future missions.
Despite the minor changes in total mass input between both studies, the more than twofold uncertainty range of our results remains striking. Neither can the applied model for the wave propagation speed of radar soundings nor the uncertainty related to the regional density profile explain the larger uncertainty of this study. The same also applies for the reduced temporal averaging time. A comprehensive evaluation of our uncertainty estimation revealed that assumptions on the geostatistical interpolation 430 error as well as grid-error dependences can have a substantial impact on the uncertainty estimation. In terms of the error partitioning, our interpolation error is the dominating source of combined grid-errors. Moreover, varying basin definitions have an impact on our total mass input estimate by up to 19%. This highlights the importance of a thorough documentation of uncertainty estimates and basin definitions to improve future intercomparisons between different SMB and mass balance inventories.
calibration, wrote the manuscript with input from all authors, V. Helm performed the SAR level_1b ASIRAS data processing, provided the digital elevation model, and established access to the RACMO2.3p2 data, E. Morris delivered the Neutron Probe density profiles, and O.
Eisen contributed to layer analysis and interpretation. All authors discussed the results and contributed to revising the manuscript. and discussions in an early stage. The authors gratefully acknowledge the provision of RACMO2.3p2 model output by S. Ligtenberg and M.
van Wessem (IMAU, Utrecht University, NL). The authors would like to thank Emerson E&P Software, Emerson Automation Solutions, for providing licenses in the scope of the Emerson Academic Program.