Constraining the recent mass balance of Pine Island and Thwaites glaciers , West Antarctica , with airborne observations of snow accumulation

In Antarctica, uncertainties in mass input and output translate directly into uncertainty in glacier mass balance and thus in sea level impact. While remotely sensed observations of ice velocity and thickness over the major outlet glaciers have improved our understanding of ice loss to the ocean, snow accumulation over the vast Antarctic interior remains largely unmeasured. Here, we show that an airborne radar system, combined with ice-core glaciochemical analysis, provide the means necessary to measure the accumulation rate at the catchment-scale along the Amundsen Sea coast of West Antarctica. We used along-track radarderived accumulation to generate a 1985–2009 average accumulation grid that resolves moderateto large-scale features (> 25 km) over the Pine Island–Thwaites glacier drainage system. Comparisons with estimates from atmospheric models and gridded climatologies generally show our results as having less accumulation in the lower-elevation coastal zone but greater accumulation in the interior. Ice discharge, measured over discrete time intervals between 1994 and 2012, combined with our catchment-wide accumulation rates provide an 18-year mass balance history for the sector. While Thwaites Glacier lost the most ice in the mid-1990s, Pine Island Glacier’s losses increased substantially by 2006, overtaking Thwaites as the largest regional contributor to sealevel rise. The trend of increasing discharge for both glaciers, however, appears to have leveled off since 2008.


Introduction
Pine Island (PIG) and Thwaites (THW) glaciers are two of the largest Antarctic contributors to recent sea-level rise (SLR) (Rignot et al., 2008;Shepherd et al., 2012) and will likely continue contributing substantially over the next century (Joughin et al., 2010;Gladstone et al., 2012).Differences between snow accumulation and ice discharge (i.e., icebergs or ice shelf melting) to the ocean define the glacier mass balance.While measuring these processes at the catchment scale was once difficult, satellite observations in the B. Medley et al.: Constraining the recent mass balance of Pine Island and Thwaites glaciers vicinity of the grounding line have improved estimates of ice discharge.Remotely sensed measurements of ice-surface velocity over the past few decades revealed that the rate of ice discharge from Pine Island and Thwaites glaciers is increasing (Rignot, 2001(Rignot, , 2008;;Joughin et al., 2003), resulting in extensive thinning near their margins (Thomas et al., 2004;Pritchard et al., 2009).This rapid dynamical change is likely the consequence of warm ocean currents melting and thus thinning the buttressing ice shelves, an effect observed over much of West Antarctica (Shepherd et al., 2004;Joughin et al., 2012;Pritchard et al., 2012;Rignot et al., 2013;Depoorter et al., 2013).While our understanding of the dynamics of these glaciers has improved substantially over the past decade, snow accumulation over large areas of these glaciers has only been sparsely sampled (van de Berg et al., 2006) limiting a complete understanding of their overall mass change.
Determining catchment-wide snow accumulation using traditional methods is difficult because rates vary considerably in space and time, and field measurements of accumulation typically sample one dimension with exclusion of the other.For example, ice-core records of accumulation (e.g., Kaspari et al., 2004) capture the temporal signal but are sparsely distributed.Stake-farm accumulation measurements (e.g., Frezzotti et al., 2005;Kameda et al., 2008;Agosta et al., 2012) are collected over broader areas to capture the spatial variability yet typically span only a few years.In addition, these in situ measurements are inadequate for mass balance studies because recovery at the catchment scale is not possible.Accumulation measurements using ground-based radar systems overcome some of the disadvantages of the traditional in situ measurements: they capture the spatial variability in accumulation over discrete (i.e., annual to multidecadal) and consistent time horizons over hundreds of kilometers (Rotschky et al., 2004;Spikes et al., 2004).Groundbased systems, however, are insufficient for regional studies because of the scale issue discussed above.Recent work by Medley et al. (2013) indicates that accumulation rates derived from airborne radar provides spatial and temporal accumulation rates over large areas, highlighting their potential for more comprehensive and improved mass balance studies.
Here, we use data from two airborne radar systems to calculate the 1985-2009 average annual accumulation over the Pine Island-Thwaites drainage system along the Amundsen Sea coast of West Antarctica.The spatial coverage limitation that makes in situ accumulation measurements disadvantageous for regional mass balance studies is overcome by aerial surveys.Using two radar systems developed by the Center for Remote Sensing of Ice Sheets (CReSIS) (Rodriguez-Morales et al., 2013), we tracked a few near-surface horizons over hundreds of kilometers of the flight surveys.The radarderived accumulation survey was spatially extensive, which enabled us to generate a complete map of the recent accumulation rate.Combining these basin-wide accumulation measurements with flux-gate estimates of ice discharge, we de-termined the recent mass balance history of the Amundsen Sea coast glaciers and their contribution to SLR.

Study area
Located in West Antarctica along the Amundsen Sea coast, the Pine Island and Thwaites catchments cover areas of 167 × 10 3 and 176 × 10 3 km 2 , respectively.Their combined extent accounts for 3 % of the grounded ice-sheet area, but receives ∼ 7 % of the accumulation (Lenaerts et al., 2012).While Pine Island and Thwaites are the primary interest, smaller adjacent catchments are investigated as well (Fig. 1).Although the Crary Mountains in the Thwaites catchment reach over 3500 m above sea level (a.s.l.), the majority of both catchments lie below 2300 m a.s.l.(Fig. 2).
These glaciers receive large amounts of snowfall because their low-elevation coastal slopes allow moisture-rich cyclones to penetrate well into the interior (Nicolas and Bromwich, 2011).Until recently, few reliable measurements of snow accumulation existed from these glaciers (Favier et al., 2013).Kaspari et al. (2004) presented accumulation records from several ice cores collected during the International Trans-Antarctic Scientific Expedition (ITASE), but only four of these lie within the Pine Island-Thwaites drainage system.Based on three of these records (one is disregarded as it is just over 20 years in length), they found that recent accumulation (between 1970 and 2000) had increased relative to the 1922-1991 average.While the recent period is relatively high, radar-derived annual accumulation shows no significant trend over Thwaites Glacier between 1980 and 2009 (Medley et al., 2013).While recent work by Burgener et al. (2013) found a negative trend in several shallow cores straddling the divide between the Ross Sea and Amundsen Sea drainages, the three cores presented in this work indicate that accumulation rates within the Pine Island-Thwaites drainage system do not exhibit a recent trend, consistent with the snow radar record from Medley et al. (2013).

Data and methods
Ground-based radar imaging of both near-surface (Sinisalo et al., 2003;Rotschky et al., 2004;Spikes et al., 2004;Eisen et al., 2005;Anschutz et al., 2007Anschutz et al., , 2008;;Frezzotti et al., 2007;Urbini et al., 2008) and deep (Nereson et al., 2000;Siegert and Payne, 2004;Waddington et al., 2007;Huybrechts et al., 2009;MacGregor et al., 2009) internal horizons has provided the basis for calculating recent and historical spatiotemporal snow accumulation rates over Antarctica.Because radar-derived accumulation measurements capture the spatial variability better than widely spaced point measurements, they provide a more accurate representation of the spatial mean, and thus are more appropriate for mass balance studies (Richardson et al., 1997).While these ground-based studies capture the spatial variability over large areas, we With the exception of Byrd, density measurements from all the ice-core sites were used to create a regional density profile in Fig. 3.The inset map shows the location of our study area.
are unaware of any surveys that were designed to map accumulation rates over entire catchment areas for the purpose of determining mass balance, as we intend to do here.For this study, we recovered three intermediate-depth firn cores, which are connected by an airborne radar survey designed to capture regional variations in snow accumulation over the entire Pine Island-Thwaites drainage system (Fig. 1).

Accumulation radar
We use two CReSIS radars in this study.The first, referred to as the "accumulation radar," is an ultra-wideband steppedfrequency chirped pulse radar system that operated between 600 and 900 MHz and is designed to image horizons in the upper 300 m of the ice sheet (Lewis, 2010;Rodriguez-Morales et al., 2013).The near-surface horizons represent contrasts in dielectric permittivity, which are likely caused by seasonal variations in the physical and chemical properties of the firn (Arcone et al., 2005)  range resolution in ice is 50 and 62 cm in firn (for density of 550 kg m −3 ).The vertical resolution of this system is too coarse to image annual stratigraphy, and consequently, an independent ice-core depth-age scale is necessary to determine horizon ages.The radar data and full documentation are available at ftp://data.cresis.ku.edu/data/accum.We processed the echograms by (1) zeroing the ice-sheet surface under the assumption that it is represented by the strongest return from each trace, (2) stacking every 12 traces to reduce noise, (3) normalizing the range bin return (i.e., quantized radar-range value) relative to the average bin return from all traces in order to brighten deeper horizons, and (4) applying a horizontal Sobel edge-detection filter to enhance horizon contrast.The Twin Otter survey took place in December 2009 through January 2010 and was designed to maximize spatial coverage over the Pine Island-Thwaites drainage area (Fig. 1) with nearly 10 000 km of flight surveys covering an area of about 300 × 10 3 km 2 .The second radar system is discussed in Sect.3.5.We calculate depth as d = 0.5cτ ε −0.5 , where τ is measured two-way travel time and c is the wave speed in a vacuum (3 × 10 8 m s −1 ), and the dielectric permittivity ε is calculated using a mixture model (Looyenga, 1965) for ice and air and is dependent on density, which increases with depth.To estimate a density profile, we fit a steady-state density model (Herron and Langway, 1980) to the average of nine firn core density profiles from the region (Fig. 3a), which is integrated to obtain a cumulative mass profile (kg m −2 ; Fig. 3c).The density profiles should approximate the steadystate density profile from the late 2000s, even though the cores were collected between 2000 and 2010, because no recent trend in accumulation is found over the past ∼ 30 years.Any variations from steady state would be manifested in the near-surface firn, which is most susceptible to the interannual variability in accumulation.The d-τ profile is generated incrementally at 1 cm intervals throughout the firn column (Fig. 3b).Previous studies have interpolated the density profile between two ice cores (e.g., Spikes et al., 2004), thereby presuming that the dominant variation is a linear trend and that the cores represent the average conditions at either end (i.e., are not biased substantially in an area of higher/lower accumulation).Because accumulation is highly variable over short distances, both the former and latter inherent assumptions are challenged.Thus, we use the regional mean profile and attempt to capture the potential error in using a prescribed profile in our error estimates.Use of a regional density profile assumes that the d-τ and cumulative mass profiles are spatially invariant; this potential source of uncertainty is discussed further below.
In order to create a spatially complete and temporally consistent map of snow accumulation, we tracked a strong and continuous reflector (H1) over as much of the radar survey as possible (Fig. 4).The depth of H1 varied considerably, ranging from 4.3 to 36.9 m.Tracking such a shallow (and thus young) horizon is possible because any undulations in the stratigraphy due to accumulation rate variations have not yet been substantially steepened over the short time interval.Using a consistent horizon over multiple flight surveys is important to generate accumulation rates over the same temporal interval.In the few areas where H1 was not traceable with confidence, we mapped other brightly visible horizons H2 and H3, which respectively are deeper and shallower than H1 (Fig. 1).We next determined the ages of these mapped isochrones using depth-age scales derived from firn cores in order to estimate accumulation rate.All horizon tracking began at the PIG2010 site where the horizons were dated using the PIG2010 depth-age scale (see Sect. 3.2).

Firn cores
Data from the radar survey were used to select sites for three intermediate-depth firn cores that we drilled during the field season following the accumulation radar survey (December 2010 and January 2011).The cores were extracted in approximately 1 m sections (diameter: 81 mm) using a Badger-Eclipse drill provided by the US Ice Drilling Program.The Pine Island glacier (PIG2010) and Thwaites (THW2010) cores were ∼ 60 m long, while the core collected along the divide between the two catchments (DIV2010) was ∼ 110 m (Table 1).Specifically, the PIG2010 site was selected because of its location at the convergence of several flight surveys (see Fig. 1).The DIV2010 was selected based on its proximity to the coast and the presence of a sequence of relatively flat radar horizons.The THW2010 core was selected on the westernmost flight path to ensure we could date  the radar horizons in case we were unable to track horizons continuously from the PIG2010 site.Density profiles for DIV2010 and THW2010 were measured in the field, while the PIG2010 profile was measured at the US National Ice Core Laboratory in Denver, Colorado.Water isotope ratios and concentrations of more than 30 elements and chemical species were measured at high depth resolution (∼ 1 cm water equivalent) using a continuous ice-core melting system (McConnell et al., 2002(McConnell et al., , 2007;;Maselli et al., 2013).Most species exhibited pronounced annual cycles (e.g., Criscitiello et al., 2013), but here we used the summer maxima in hydrogen-peroxide concentration, water-isotope ratios, and non-sea-salt sulfur to sodium ratio to identify annual layers.Known volcanic horizons identified by marked increases in wintertime sulfur concentration provided verification of the annual layer counting, indicating a dating uncertainty of less than 1 year.
The PIG2010 core was selected to date horizons because of its optimal location at the intersection of several radar surveys (Fig. 1), and we evaluated the isochronal accuracy by dating H1 at the DIV2010 and THW2010 core sites.The depths of H1 at the PIG2010, DIV2010, and THW2010 cores were 19.75 ± 0.33, 18.95 ± 0.35, and 14.15 ± 0.34 m (error calculations described below), respectively, which correspond to ages of 25.4 ± 0.4, 25.7 ± 0.4 yr, and 25.2 ± 0.6 yr.Although the depth uncertainties at each core are comparable, the age uncertainty for THW2010 is larger because of the relatively low accumulation rate at this site.As a result, the THW2010 depth-age curve is shallower (Fig. 3d), which translates into a larger age uncertainty.Nonetheless, the H1 ages differ by 0.5 yr, which is remarkable given that the survey distance between cores is 750 km.This comparison confirms that radar-detected horizons are isochronous over large distances, consistent with others studies from this region (Arcone et al., 2004;Spikes et al., 2004).

Accumulation rate calculations
At shallow depths, the spatial variation in the depth to a given horizon is largely driven by variations in the accumulation rate and, to a lesser extent, the density profile and the firn compaction rate, which are also dependent on the accumulation rate.The depth variations are combined with firn density information to extract the accumulation rate along the radar survey.The long-term accumulation rate (between the horizon and surface) is determined by dividing the cumulative firn mass per unit area (Fig. 3c) above the horizon by the time since horizon burial (i.e., the horizon age in years).While the horizon age does not vary spatially, the horizon depth does vary, which results in variable firn mass above the horizon.Over the survey portions where we were unable to map H1, we measure accumulation using an alternate horizon (H2 or H3) and correct for the temporal bias.The bias corrections were based on accumulation measurements where both H1 and the alternate horizon were coincidently mapped.We completed a total of five robust regressions (Fig. 5), one for each flight survey segment where an alternate horizon was used to measure accumulation (see Fig. 1).The different relationships are the potential result of (1) different temporal biases from using two alternate horizons and (2) spatial variations in the bias.These corrected data make up only 8 % of the total accumulation measurements and are all located in the northern Pine Island Glacier catchment.

Radar-derived accumulation rate errors
Uncertainty in radar-derived accumulation rates arise from uncertainties in the regional density profile and horizon age.At any location, the deviation of the actual density profile from the regional mean translates into errors in the cumulative mass and d-τ profiles and ultimately the measured accumulation rate.To account for this error, we fit the aforementioned density model to the ±1σ deviation of the measured density profiles from the mean (Fig. 3a).We then calculate the error in the d-τ and cumulative mass profiles (Fig. 3b, where an alternate horizon was used to measure accumulation.Using a robust regression model, we correct accumulation measurements derived from the alternate horizon (either H2 or H3) to more appropriately represent its H1 measurement where both horizons were coincidently tracked.The relationship was then applied to the measurements derived from the alternate horizons where H1 was not tracked.Four of the corrections use the same alternate horizon H3, whereas one leg uses H2.Interestingly, H2 dates to 1939 and the resulting accumulation rates are found to be much lower than those from H1 (i.e., the regression model lies well above the 1 : 1 line).This result is consistent with ice-core observations that recent accumulation has increased relative to the long-term mean.Corrections for measurements derived from H3, dated to 2002, hover around the 1 : 1 line, especially in the areas of the majority of the measurements (0.3-0.5 m w.e.yr −1 ).
c) assuming that the density uncertainty could bias our results, which means that errors accumulate with depth.This assumption is conservative and reasonable based on evaluation of the individual core profiles relative to the regional mean.Finally, a digitization error of ±1 range bin is included in the depth error.Uncertainty in the age of the horizon also introduces an uncertainty into our accumulation measurements.Using the regional mean density profile and its uncertainty, we determined the H1 depth and error at each of the three core sites using their measured depth-age profiles shown in Fig. 3d (see above).While the measured age at each site falls within the error bounds of the other two sites, the range of values is large enough that we must consider its impact.We assign the error in the age of H1 at ±1 year, which is likely an overestimation based on the evaluation of the isochronal accuracy.
Finally, we must consider uncertainty in the bias correction for measurements based on the alternate horizons (H2-H3).These measurements are assigned errors equal to the root mean square error of the robust regression fits shown in Fig. 5, which vary from 0.018 to 0.069 m water equivalent (w.e.) yr −1 depending on survey leg.

Snow radar
The second CReSIS radar used in this study -referred to as the "snow radar" -is an ultra-wideband microwave radar that operated over the frequency range of 4-6 GHz in 2009 and 2-6.5 GHz in 2010 and 2011.This system is a frequencymodulated continuous-wave (FM-CW) radar that images the stratigraphy in the upper 20-30 m of the ice sheet with fine vertical resolution.The theoretical vertical range resolution for the 2009 (2010/11) survey is 8 cm (4 cm) in ice and 10 cm (5 cm) in firn (Panzer et al., 2013), which is much finer than the accumulation radar.The snow radar was flown as part of NASA's Operation IceBridge, which focused on areas of rapidly changing ice in and around the major outlet glaciers.While the survey was not designed to extend over the entire catchment, it provides additional measurements over ∼ 2000 km of the survey tracks primarily within the Thwaites catchment.To be temporally consistent with the accumulation radar measurements, here we use the 1985-2009 mean annual accumulation derived from the snow radar (Medley et al., 2013).

Interpolation
While the along-track resolution of the radar data is ∼ 10 s of meters, we extract horizon depths at a sampling interval of 500 m, which is small relative to the catchment size.Large data gaps still remain (up to 150 km between flight paths; Fig. 6a).The large gaps mean that the spatial resolution of an interpolated map will be substantially coarser than the alongtrack resolution: accumulation was not appropriately sampled to recover high-frequency (< 10 km) variations at the catchment scale.To minimize the high-frequency variability in the accumulation measurements, we applied a 25 km running average filter to the profile data (Fig. 6b).Approaching the ends of the surveys, the filtering length was tapered down to 5 km to maximize the spatial coverage for interpolation.Using the smoothed accumulation measurements, we generated a gridded map of accumulation using the geostatistical interpolation technique of kriging (Leuangthong et al., 2008).Prior to interpolation, we used an ordinary least squares (OLS) linear regression model with northing, easting, and elevation as explanatory variables to create an accumulation rate surface (Fig. 7a).The interpolation was then performed on the OLS model residuals in order to recover the moderate-scale features.The distribution of the residuals showed only a small degree of skewness (< 0.5), and, as a result, we did not perform a log transformation on the data as done by Arthern et al. (2006).We found that the best fit to the measured semivariogram was an isotropic spherical model with a range of 175 km.Sharp lines and edges in the radar survey result in unrealistic artifacts in the interpolated map.Therefore, we smoothed the 3 km grid with a 9 × 9 cell mean filter to minimize these high-frequency interpolation artifacts.Finally, the OLS surface is added to the kriged residuals to generate the final accumulation map shown in Fig. 7b.
We created an accumulation error grid that accounts for measurement and interpolation uncertainties (Fig. 7c).The kriging standard prediction error is based on the distances to the nearest measurements (i.e., cells farther from measurements have a greater error) and the spatial structure of the data as described by the semivariogram.We also investigate the impact of measurement error on the final accumulation map.Random error is added to each accumulation measurement and these perturbed values are then interpolated to a grid using the same parameters described above.The error added to each measurement point is randomly selected from a normal distribution with a mean of zero and standard deviation equal to the measurement error for that data point.This process was repeated 200 times, and the measurement error for each grid cell was taken as the standard deviation of these 200 realizations.The final accumulation error grid was generated by root-sum square (RSS) of the interpolation and measurement grids and was smoothed using a 9 × 9 cell mean filter.

Surface velocities and catchment discharge
We derived surface velocities from 1994 to 2012 using a combination of interferometric synthetic aperture radar (InSAR) and speckle-tracking techniques (Joughin, 2002).Velocities from 2000 and before were determined using data from the European Space Agency's ERS-1/2 (European Remote Sensing satellites) mission and later velocities were derived from a combination of data from the Japanese ALOS (Advanced Land Observing Satellite) and German TerraSAR-X (Terra synthetic aperture radar) mission (Joughin et al., 2003(Joughin et al., , 2010)).Our velocity measurements are made at the ice-sheet surface but we assume they are equivalent to column average velocity due to the high degree of sliding at near the grounding line.Based on analysis of estimated deformation velocities that are internal variables of a temperature model (Joughin et al., 2009), any biases introduced by this assumption are less than 1 % and are not included in our error analysis.System noise produces errors of ∼ 10 m yr −1 and there are additional velocity and slopedependent errors of ∼ 3 % (Joughin, 2002).
Catchment discharge was estimated using ice-surface velocities and a time-varying estimate of ice thickness along a transect that is roughly parallel to the grounding line.Ice-thickness estimates are sometimes confounded by the presence of basal crevasses downstream of the grounding line and, because grounding lines have retreated during the decade considered in this study (Joughin et al., 2010), the transect was displaced 5-10 km inland to ensure that highquality ice-thickness data were available.
We constructed a high-resolution model of the timevarying ice-surface height using ICESat (Ice, Cloud,and land Elevation Satellite) altimetry data (Zwally et al., 2012), and airborne scanning laser altimetry data supplied by NASA's IceBridge program (Blair and Hofton, 2012;Krabill, 2013).We combined our ice-surface height estimates with icethickness estimates derived from ice-penetrating radar (Holt et al., 2006;Vaughan et al., 2006;Allen, 2013) to estimate a set of bed elevations.We applied a minimum-curvature gridding technique that minimized, in a least-squares sense, the second spatial derivative of the bed elevation, while also minimizing the data misfit.The cost function on the derivatives was selected based on ice-velocity maps so that curvature in the along-ice-flow direction was penalized more heavily than the curvature in the across-flow direction, giving bed-elevation estimates that preserve channel structures while suppressing small-scale noise in the data.The rms (root mean square) misfits between data and the fit surface were better than 7 m over the smooth basal topography near the grounding line.This small misfit suggests that the fit surface adequately resolves the details of the bed topography; uncertainties in the data picking and in the location of radar footprints contribute substantially larger errors to the ice thickness, which we conservatively estimate at 50 m.
We derived ice-discharge estimates using surface-height, surface-velocity, and ice-thickness estimates assuming that ice flow is almost entirely due to sliding at the bed: (1) Here u•n is the component of the ice-surface velocity perpendicular to the transect, z s is the surface height at the time the velocity was measured, b is the bed elevation, and h air is the depth-integrated thickness of air contained in the firn, as estimated from van den Broeke (2008).We evaluate the flux integral on points spaced every 50 m along the transect.
When the velocity maps contain gaps, we interpolate spatially within the same map to close gaps smaller than 4 km, and interpolate in time between temporally adjacent maps to fill larger gaps.Velocity values so interpolated are assigned an additional error component of 100 m yr −1 in each direction.If temporally adjacent maps do not supply a valid velocity estimate, the velocity is estimated from the mean of all available velocities for that point, and the error estimate is set to 250 m yr −1 , which happens only for a few points at the north edge of PIG and a few points along the Wedge.Because no usable velocity data were available for the tributary of PIG (see Fig. 1) in 2006, we estimated the flux and its error for that part of the profile using a linear interpolation between the 2000 and 2009 values.Because the variation in flux for this part of the profile is on the order of 1 Gt between 2000 and 2009, the error incurred by this interpolation should be small.
We estimate errors in our flux estimate by propagating the measurement errors in Eq. ( 1), assuming components between 50 m grid points were independent.We then account for spatial correlation in the errors, which we conservatively

B. Medley et al.: Constraining the recent mass balance of Pine Island and Thwaites glaciers
assume to be spatially correlated on a 10 km scale, by multiplying these initial error estimates by (10 km/50 m) 1/2 , or about a factor of ∼ 14.Other reasonable choices of a correlation scale would reduce the error estimates roughly in half, but significantly larger errors are unlikely.

Radar-derived accumulation measurements
Mean annual radar-derived accumulation (including the bias correction) over the period 1985-2009 are shown in Fig. 6a.The ∼ 20 000 accumulation measurements span an order of magnitude ranging from 0.13 to 1.37 m w.e.yr −1 and have a mean (± standard deviation) of 0.41 ± 0.12 m w.e.yr −1 .The ∼ 6300 measurements within the Pine Island Glacier catchment vary from 0.18 to 1.37 m w.e.yr −1 with a mean of 0.43 ± 0.15 m w.e.yr −1 , and the 11 400 measurements within the Thwaites Glacier catchment area vary from 0.21 to 0.84 m w.e.yr −1 with a mean of 0.42 ± 0.08 m w.e.yr −1 .Rates exceeding 1.0 m w.e.yr −1 are found along coastal Pine Island glacier and in isolated surface depressions.Rates below 0.2 m w.e.yr −1 occur on bumps alongside those depressions as well as across the Thwaites southern divide toward WAIS (West Antarctic Ice Sheet) and Byrd camps.Measurements outside the Pine Island-Thwaites catchment area (n = 1568) are used in the OLS regression and interpolation.The average 1σ accumulation-measurement errors within the Pine Island and Thwaites catchments are 0.03 and 0.02 m w.e.yr −1 respectively.While the minimum error in each catchment is the same (0.01 m w.e.yr −1 ), the maximum is much greater within the Pine Island catchment (0.17 m w.e.yr −1 ) than over that of Thwaites (0.04 m w.e.yr −1 ), which is due to the larger accumulation rates and associated bias correction on Pine Island glacier.
Radar-derived accumulation rates (±1σ error) at the PIG2010, DIV2010, and THW2010 sites are 0.43 ± 0.02, 0.41 ± 0.02, and 0.29 ± 0.01 m w.e.yr −1 , respectively, which match those derived from the core records shown in Table 1 (0.42, 0.41, and 0.29 m w.e.yr −1 ).The nearly identical measurements at the PIG2010 site is not surprising because the core-derived depth-age scale was used to determine horizon ages and its density profile was one of nine used to determine a regional profile.The only information used from the DIV2010 and THW2010 cores was their density profiles; the radar-derived measurements at these cores are largely independent of the core-derived accumulation rates.At both sites, the core measurements fall within the radar-derived error interval.

Gridded accumulation rates
Not surprisingly, the OLS accumulation rate surface (Fig. 7a) shows relatively high accumulation at low elevations, which decreases progressively inland towards higher elevations.
While the general structure is correct, there are several moderate-scale (25-50 km) features that are not reproduced with the simple OLS model, which is clearly apparent when comparing the OLS model (Fig. 7a) with the smoothed measurements (Fig. 6b).The final gridded accumulation map (Fig. 7b) reproduces both the regional and moderate-scale features observed in the measurements and will provide the snow input values for our mass balance estimates.Unlike the OLS model, the final grid captures the precipitation shadow effect that is apparent, for example, along the northern slopes of the Pine Island catchment and is caused by the mountain ranges of Eights Coast.
The average (± standard deviation) gridded accumulation rates over Pine Island and Thwaites glaciers are 0.40 ± 0.13 and 0.43 ± 0.09 m w.e.yr −1 , values similar to those from the radar-only measurements (Table 2).The accumulation grid errors range from 2.6 to 32.7 % with an average of 8.6 %.Nearly 90 % of the cells have errors of less than 15 % (Fig. 7c).Maximum errors are found over the southern sector of the upper Pine Island catchment where the accumulation rates are very low and the radar coverage is sparse.The lowest errors are in central Thwaites where multiple overlapping flight paths provide better spatial coverage.
We compared several core-derived accumulation rate averages from within the Pine Island-Thwaites drainage system to their coincident gridded estimate.Specifically, we used the ITASE 01-1, 01-2, 01-3, and 01-6 average rates between 1985 and 2001 and the PIG2010, DIV2010, and THW2010 average rates between 1985 and 2009.At five of the seven sites, the average rate falls within the 1σ grid error.The other two sites (ITASE 01-3 and PIG2010), which interestingly are separated by only 20 km, fall within the 2σ error.Even though the grid has been smoothed removing small-scale accumulation features, it still matches isolated core measurements very well, which gives confidence that our gridded accumulation rates and errors are reasonable.

Accumulation distribution by elevation
We next investigate the elevation-dependent accumulation distribution for the Pine Island (Fig. 8) and Thwaites (Fig. 9) catchments by binning their accumulation grids over 100 m intervals.The average accumulation rates over both glaciers decrease consistently with increasing elevation, but the average rate over a given elevation bin is larger for Thwaites than for Pine Island (Figs. 8a,9a).Although accumulation decreases with elevation, the differences in catchment hypsometry indicate the elevations that contribute most to the total snow accumulation fall between 700 and 1500 m over the Pine Island catchment and between 1200 and 1900 m over that of Thwaites (Figs. 8b,9b).According to our results, totals of 67.3 ± 6.1 and 75.9 ± 5.2 Gt yr −1 accumulated on average between 1985 and 2009 over Pine Island and Thwaites, respectively, and 158.6 ± 12.5 Gt yr −1 accumulated over the entire region (Table 3).We assumed the gridded errors were not independent, and as a result, the errors were calculated by cell-by-cell summation (not RSS) of the grid errors.Therefore, the error bounds are likely a conservative estimate.

Comparison with climatologies and atmospheric models
The elevation-dependent accumulation distributions for various accumulation climatologies, reanalysis products, and a regional climate model output were also analyzed and compared.The model-and observation-based climatologies are available at the Antarctic Cryospheric Access Portal (A-CAP; Scambos et al. (2008), van de Berg et al. ( 2006), Arthern et al. (2006), and Monaghan et al. (2006)), which are henceforth referred to as VDB05, ART06, and MON06.While VDB05 and MON06 are model-derived, we separate them from the reanalysis and climate models described below because these products are provided by A-CAP as long-term averages and their temporal coverage is not consistent with our radar-derived measurements.The ART06 map was generated using field-based measurements of snow accumulation, which were gridded using remotely sensed microwave emission data to guide the interpolation.The MON06 mean annual (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)) simulated precipitationminus-sublimation (P-S) estimate is derived from the Polar MM5 atmospheric model forced by the European Centre for Medium-Range Weather Forecasts (ECMWF) 40-year Reanalysis (ERA-40).Finally, VDB05 is the 1958-2002 mean annual simulated surface mass balance (SMB) from the Regional Atmospheric Climate Model RACMO2 forced at its lateral boundaries by ERA-40.
The three global reanalysis P-S products include the ECMWF "Interim" Reanalysis (ERA-Interim) (Dee et al., 2011), the NASA Modern Era Retrospective Analysis for Research and Applications (MERRA) (Rienecker et al., 2011), and the National Centers for Environmental Prediction Climate Forecast System Reanalysis (CFSR) (Saha et al., 2010).Finally, we use SMB from a recent RACMO2 simulation that is forced at the lateral boundaries with the ERA-Interim (Fig. 7d) (Lenaerts et al., 2012).Even though these products do not all estimate precisely the same variable (i.e., accumulation, SMB, or P-S), they are all nearly equivalent to snow accumulation in this region (Medley et al., 2013).To ensure consistency, all grids were bilinearly resampled to the same 3 km grid and were binned as explained above (Figs.8, 9).Over both catchments, the high-elevation average accumulation rates for nearly all the products are lower than those from our grid.For Thwaites, all products have higher lowelevation accumulation rates than our grid, with the exception of ART06.Because our grid is not based on measurements below 800 m from Thwaites, we cannot confidently determine whether the products have higher accumulation at those elevations; however, they are higher than our gridded values between 800 and 1200 m.For Pine Island, a range of differences occurs at the low elevations: some are higher, some are lower, and others fall within our gridded values.In general, these accumulation products (except ART06) have steeper elevation-dependent accumulation gradients than our grid.
The largest spread in average accumulation between these products occurs at the lowest elevations (below ∼ 600 m), but these elevations occupy a relatively small area of the large basins and thus do not contribute substantially to the spread in their cumulative accumulation rates (Figs.8b, 9b).With the exception of MON06 and MERRA, the products fall within the error range of the our total cumulative accumulation rate for the Pine Island glacier catchment, albeit towards the low end (Fig. 8c).Only CFSR, RACMO2, and VDB05 fall within this range for Thwaites where the spread is much larger (Fig. 9c).Although several of these products generate values similar to our grid (Table 3), it is often the result of low-elevation regions of larger rates balancing high-elevation regions of lower rates, which is most apparent in Figs.8b and  9b.

Ice discharge and mass balance
The total flux of ice lost to the ocean from this region increased from 192.1 ± 6.0 Gt yr −1 in the mid-1990s to 257.4 ± 4.8 Gt yr −1 in 2010 (Table 4), which is a more than 30 % increase over ∼ 15 years and is consistent with earlier estimates (Rignot, 2008).The ice discharge of the Wedge from 2000 (9.7 ± 2.8 Gt yr −1 ) was used in the regional estimation of ice discharge for the mid-1990s.Over the www.the-cryosphere.net/8/1375/2014/The Cryosphere, 8, 1375-1392, 2014  same interval (mid-1990s-2010), discharge from Pine Island alone increased more than 50 % whereas Thwaites increased just under 20 %.Although between the mid-1990s and mid-2000s only a few data points exist, this is likely the period over which the discharge increased substantially and was followed by a period of relatively steady flow between 2008 and 2010 for Thwaites and between 2008 and 2012 for Pine Island (Table 4).From 2000 to 2010 the Wedge increased discharge by just over 20%, and from the mid-1990s to 2010 Haynes increased by more than 20 % and Thwaites East by 50 %.
The multiple ice discharge measurements presented show a strong trend towards more negative values, whereas no recent     2013), our accumulation estimates for a given discharge time interval could be biased by as much as 25 %.The integrated mass balance and sea-level measurements over the entire 1994-2010 interval should not be affected since no accumulation trend is observed.Therefore, the mass balance trends presented are entirely determined from the trends in discharge and we assume that, while year-to-year accumulation can vary, the interannual variability is not significant when considering the entire interval.
Based on those assumptions, the mass loss from the Pine Island-Thwaites drainage system nearly tripled between the mid-1990s and 2010, increasing from 33.5 ± 13.9 to 98.8 ± 13.4 Gt yr −1 , values which correspond to +0.09 ± 0.04 and +0.27 ± 0.04 mm SLR yr −1 (Table 5; Fig. 10).During the mid-1990s, the mass balance of the Pine Island glacier catchment was slightly negative (−6.0 ± 6.4 Gt yr −1 ) with errors large enough to suggest that it was at or near balance.Thwaites, Thwaites East, and Haynes glaciers all showed negative imbalances.Thwaites Glacier was farthest out of balance at −17.3 ± 7.1 Gt yr −1 followed by Haynes at −8.3 ± 1.3 Gt yr −1 .By 2010, Pine Island Glacier was significantly out of balance, losing 46.1 ± 7.1 Gt yr −1 and overtaking Thwaites (−35.7 ± 5.8 Gt yr −1 ) as the largest  3. The actual mass balance values are listed in Table 5.The points are offset slightly in time for clarity and do not represent actual differences in the data collection period.
contributor to SLR.The Thwaites East and Haynes mass balances decreased as well, but the Wedge remained essentially in balance.

Discussion
The more than 20 000 radar-derived accumulation measurements reveal both regional and local features that have not been uncovered using the existing firn-core measurements alone.At the same time, the dated firn cores remain crucial for our analysis because they provide the depth-age scale necessary to date the radar horizons and the firn depthdensity profiles.Tracking the horizons over 100s of kilometers proved successful except over a few areas.Notably, horizons disappeared northward from the PIG2010 site into the area of enhanced ice flow and rougher surface undulations around the Pine Island trunk.Additionally, we were unable to differentiate with confidence between horizons moving westward from WAIS to Byrd because the accumulation rate is substantially lower at Byrd, resulting in the merging of horizons.The areas of more extreme surface undulations, as indicated by the tonal differences in the base map in Fig. 1, often coincide with data gaps where the horizons could not be tracked.Nonetheless, we were able to track horizons over the majority of the Pine Island-Thwaites drainage system and over a wide range of elevations and accumulation rates.Outside of regions with large accumulation gradients, the accumulation radar likely should image continuous and discretely trackable horizons that, when combined with ice cores and a well-defined survey, should provide catchment-wide accumulation measurements elsewhere in Antarctica.
Smoothing of the raw accumulation measurements and grid filtering indicates our map contains moderate-to large-scale accumulation features with scales on the order of 25 km or greater.Smaller-scale (< 25 km) features certainly exist as evidenced in the echograms and raw accumulation measurements, but a denser airborne survey would be required to capture these features over the large catchment areas.If we consider the small-scale features as high-frequency noise, the catchment-wide accumulation measurements are not negatively affected.Because of the moderate-scale data resolution (∼ 25 km), we do not expect ice-core accumulation measurements to match the coincident grid accumulation with high fidelity.For example, the PIG2010 core site is clearly located in a minor depression where accumulation rates are enhanced relative to the background (Fig. 4) and is one of the cores with the greatest mismatch from our grid.At the same time, accumulation rates from most cores are within the ±1σ grid error; thus, in the areas where the firn cores were retrieved, our final accumulation grid is a good representation of the moderate-scale accumulation rate with errors suitably large to encompass any smaller-scale features not resolved.
Beyond determining catchment-wide accumulation, our grid provides the means with which to test the abilities of various climatologies, reanalyses, and climate models to reproduce the accumulation rate over a large area.Rather than comparing the models to isolated point measurements of accumulation, which do not make for equal comparison with large grid cell values, our grid is of comparable data resolution (i.e., is more representative of the scales resolved by the models).A consistent feature amongst the reanalysis and climate models is the lower values of accumulation in the highelevation, low-accumulation interior.Whereas, our gridded product suggests lower values than those from the models in the low-elevation, high-accumulation coastal areas.Unlike in the interior, the discrepancy at low elevations does not lead to a substantial spread in the estimated basin-wide accumulation between the various products because the low elevations cover relatively small areas in our study area.We suggest the models are likely closer to the truth than our grid for Haynes and Thwaites East, as well as the lowest portions of Pine Island and Thwaites catchments.Not only do we lack measurements below 500 m a.s.l., the slopes are steepest at the low elevations, resulting in enhanced orographic lifting (Nicolas and Bromwich, 2011) and this effect is not incorporated in our interpolations scheme (i.e., slope was not used as a dependent variable in the OLS model).Additional accumulation measurements from lower elevations are necessary to properly assess model abilities in this high-accumulation area.Nevertheless, RACMO2 and CFSR are able to generate catchment-scale accumulation rates that fall within our grid error bounds for both Pine Island and Thwaites.
The ice discharge measurements presented here are in excellent agreement with those from Rignot (2008), showing a strong increase from Pine Island and a moderate increase from Thwaites between the mid-1990s and 2007.The additional values presented in this study from 2008 to 2012 indicate that the increasing trends have not continued for both glaciers and that ice discharge has, for the time being, leveled off.Several model simulations from this region (e.g., Joughin et al., 2010;Favier et al., 2014;Joughin et al., 2014) exhibit a similar pattern of accelerated loss followed by a period of sustained and stable high loss.While Thwaites had the largest mass imbalance between 1994 and 1996, the mass loss from Pine Island overtook that from Thwaites by 2006, making Pine Island the largest contributor to SLR from the region.Haynes, a relatively small glacier, shows large mass loss, which is likely overestimated because the interpolated accumulation rates from this sector are likely biased low.Nevertheless, Haynes Glacier would still show mass losses over the entire period even if we doubled our accumulation rates over the entire glacier.The same is true of Thwaites East but due to its small size, the imbalance is not as large.The Wedge, between Pine Island and Thwaites glacier basins, is relatively stable and according to our estimates in balance, receiving some of the highest accumulation rates in the region.While Thwaites East and Haynes glacier mass losses are dwarfed by those from Pine Island and Thwaites, their ice discharges have increased substantially since the mid-1990s (50 and 22%, respectively).If the increases continue into the future, these small glaciers are likely to have a substantial impact on SLR.In total, our gridded data show that the region receives on average 158.6 ± 12.5 Gt yr −1 of snow accumulation at an average accumulation rate of 0.43 m w.e.yr −1 for the 1985-2009interval. Between 1994and 1996, , total ice discharge was 192.1 ± 6.0 Gt yr −1 resulting in an annual mass change of −33.5 ± 13.9 Gt yr −1 .In 2010, ice discharge increased substantially to 257.4 ± 4.8 Gt yr −1 yielding a much larger mass change of −98.8 ± 13.4 Gt yr −1 .These mass imbalances translate into rates of SLR of 0.09 ± 0.04 (1994)(1995)(1996) and 0.27 ± 0.04 mm yr −1 (2010).To compare with the recent ice-sheet mass balance intercomparison exercise (IM-BIE) (Shepherd et al., 2012), which assessed the mass balance of Greenland and Antarctica from 1992 to 2011, we make the following assumptions about our discharge history: (1) the ice discharge from 1992 to 1998 is constant and equal to our 1994-1998 measurement; and (2) in subsequent years, the discharge steps up or down at the beginning of the next measurement interval.For example, the Thwaites Glacier discharge is set to 93.2 ± 4.8 Gt yr −1 from 1992 to 2005 and jumps to 104.3 ± 2.6 Gt yr −1 in 2006.These assumptions result in a conservatively low estimate of the 1992-2011 ice discharge because we assume all changes occur instantaneously.Subtracting the mean ice discharge over the entire record from our catchment-wide accumulation rates provides the 1992-2011 mean mass balances, which are −19.4 ± 6.1 and −22.0 ± 5.3 Gt yr −1 for Pine Island and Thwaites glaciers, respectively.These mass balance measurements match the IMBIE estimates well for Pine Island glacier, but our estimate for Thwaites Glacier is on the more negative end of the IMBIE range (i.e., shows greater mass loss).As a result, we have Thwaites Glacier losing more mass on average than Pine Island, which is reversed relative to the IMBIE measurements.This discrepancy could be the result of different catchment boundaries.Our results indicate the region as a whole has contributed ∼ 3 mm to SLR over the 1992-2011 period, which amounts to 27 % of the total contribution to SLR of 11.2 mm from both Greenland and Antarctica as determined by IMBIE.

Conclusion
We find that a well-designed accumulation radar survey combined with glaciochemical analysis of one or more well-sited firn cores is sufficient to generate a catchment-wide accumulation map that resolves moderate-to large-scale features.We found that various climatologies and reanalysis and climate models have lower accumulation rates than our gridded values in the high-elevation interior and potentially higher rates in the low-elevation coastal areas, consistent with our prior finding (Medley et al., 2013).These discrepancies often cancel out each other, resulting in reasonable estimate of catchment-wide accumulation.Between the mid-1990s and 2010, the mass balance of the region decreased from −33.5 ± 13.9 to −98.8 ± 13.4 Gt yr −1 , a near tripling of its imbalance and associated contribution to SLR.Although contribution to SLR from Pine Island Glacier exceeded that from Thwaites Glacier in 2006, Thwaites showed greater mass loss on average between 1992 and 2011.Although both glaciers experienced a substantial increase in ice discharge between the mid-1990s and 2008, the trend of increasing ice discharge ended, at least for now, around 2008.

Figure 1 .
Figure 1.The Amundsen Sea coast glaciers and locations of the radar flight surveys.Here, the MODIS mosaic is overlaid transparently by measured ice velocities from Joughin et al. (2010) and elevation contours (200 m intervals).The catchments are outlined and labeled.The complete accumulation radar survey consists of the white and grey lines.A dashed white line indicates no horizon was mapped, a solid white line indicates H1 was mapped, and light (dark) grey indicates that H2 (H3) was mapped.The solid black lines show where snow radar accumulation measurements were taken from Medley et al. (2013).The three 2010 ice cores are labeled and indicated by a yellow circles.The ITASE cores are displayed as blue triangles, and the WAIS divide and Byrd camps are labeled and indicated by light purple squares.With the exception of Byrd, density measurements from all the ice-core sites were used to create a regional density profile in Fig.3.The inset map shows the location of our study area.

Figure 2 .
Figure 2. The hypsometric distributions of the Pine Island and Thwaites glacier catchments.The sub-basins included in PIG are Pine Island glacier and the Wedge and in THW the Thwaites, Thwaites East, and Haynes glaciers.The median elevations within the PIG and THW are 1210 and 1540 m a.s.l., respectively.

Figure 3 .
Figure 3. Profiles of (a) density, (b) two-way travel time, (c) cumulative mass, and (d) age with depth.(a) The density profiles from nine firn cores from the region are plotted in light grey and the model fit to their mean is shown as a solid black line.The dashed black lines show ±1 standard deviation (σ ) from the mean.(b) The solid line was produced using the formula, d = 0.5cτ ε −0.5 , for conversion between two-way travel time and depth using the density model in (a), and the dashed lines were generated from the ±1σ deviations in (a).(c) The cumulative mass profiles were created by integrating the mean (solid) density profile with depth as well as the ±1σ deviations (dashed) from (a).(d) The depth-age profiles for the three 2010 cores determined from glaciochemical analysis.

Figure 4 .
Figure 4. Echograms at each 2010 core site and the tracked H1 horizon.From left to right, we display echograms from the PIG2010, DIV2010, and THW2010 sites along with H1 mapped in dashed red and H2 and H3 in dashed green.The vertical white line shows the location closest to each core, which extends from the surface to the actual recovery depth.

Figure 5 .
Figure 5. Bias correction regression models for five survey legs where an alternate horizon was used to measure accumulation.Using a robust regression model, we correct accumulation measurements derived from the alternate horizon (either H2 or H3) to more appropriately represent its H1 measurement where both horizons were coincidently tracked.The relationship was then applied to the measurements derived from the alternate horizons where H1 was not tracked.Four of the corrections use the same alternate horizon H3, whereas one leg uses H2.Interestingly, H2 dates to 1939 and the resulting accumulation rates are found to be much lower than those from H1 (i.e., the regression model lies well above the 1 : 1 line).This result is consistent with ice-core observations that recent accumulation has increased relative to the long-term mean.Corrections for measurements derived from H3, dated to 2002, hover around the 1 : 1 line, especially in the areas of the majority of the measurements (0.3-0.5 m w.e.yr −1 ).

Figure 6 .
Figure 6.The (a) raw and (b) smoothed radar-derived accumulation rates.(a) The raw accumulation rates derived from the accumulation radar, including the bias-corrected rates as well as the 1985-2009 mean annual accumulation from the snow radar.(b) Same as from (a) except a 25 km (tapered to 5 km approaching the ends) running average has been applied.
These give irregular spatial and temporal coverage between the late fall of 2002 and 2012, which we integrated into an estimate of surface elevation and elevation change by fitting a DEM (digital elevation model) surface for 2010 and a series of correction surfaces giving height differences between 2002 and 2010.The solution was selected to minimize, in a least-squares sense, the difference between the model surface and the observed surface heights, while also minimizing the second derivative of the DEM surface and the second derivative of the ice-surface change rate between any pair of years.This technique provides an estimate of surface heights at any time and at any position on the grounded ice.The accuracy of any estimate depends on the temporal and spatial sampling of the input data; typically, surface-elevation error estimates for points within 2-4 km of a flight line are less than 10 m.Ice-surface elevations for flux estimates before 2003 are calculated by linear extrapolation of the 2003-2007 elevation rate of change.This extrapolated elevation difference is assigned an error of 100 %.

Figure 7 .
Figure 7. Accumulation and error grids including (a) the accumulation surface derived from an OLS linear regression model with northing, easting, and elevation as dependent variables; (b) our final accumulation surface derived by adding the kriged OLS model residuals to the OLS regression surface from (a), (c) the combined measurement and interpolation errors displayed as a percentage of our gridded accumulation map shown in (b), and (d) the 1985-2009 average surface mass balance from RACMO2, widely used in mass balance studies, is shown for comparison to our gridded map in (b).

Figure 8 .
Figure 8.The elevation-dependent accumulation distribution for Pine Island glacier (including the Wedge) and comparison with climatologies and reanalysis and climate models.For each part, the grey shaded area shows the quantity of interest from our final accumulation grid including its ±1σ deviation.The climatologies and reanalysis and climate models do not have errors because the products do not provide error grids.The three part figure shows (a) the average accumulation rates over 100 m elevation bins, (b) the bin-summed accumulation rates scaled by the cell size of 9 km 2 with bars representing the bin size, and (c) the cumulative bin-summed accumulation rates from (b).The solid vertical black lines on each plot display the elevation limits of our radar-derived accumulation measurements, and thus bound the area with high data integrity.

Figure 9 .
Figure 9.The elevation-dependent accumulation distribution for Thwaites glacier catchment (including the Thwaites East and Haynes) and comparison with climatologies and reanalysis and climate models.For each part, the grey shaded area shows the quantity of interest from our final accumulation grid including its ±1σ deviation.The climatologies and reanalysis and climate models do not have errors because the products do not provide error grids.The three part figure shows (a) the average accumulation rates over 100 m elevation bins, (b) the bin-summed accumulation rates scaled by the cell size of 9 km 2 with bars representing the bin size, and (c) the cumulative bin-summed accumulation rates from (b).The solid vertical black lines on each plot display the elevation limits of our radar-derived accumulation measurements, and thus bound the area with high data integrity.

Figure 10 .
Figure 10.Mass balance history from 1994 to 2012 for each catchment.The mass balance was measured by subtracting the ice discharge measurements from various time intervals from Table 4 from the 1985-2009 mean catchment-wide accumulation from Table3.The actual mass balance values are listed in Table5.The points are offset slightly in time for clarity and do not represent actual differences in the data collection period.

Table 1 .
Summary of ice-core accumulation records.

Table 3 .
Catchment-wide annual accumulation rate from climatologies, reanalysis products, and a climate model compared to this study.

Table 4 .
Flux-gate discharge measurements and errors from 1994 to 2012.