Annual Greenland accumulation rates ( 2009 – 2012 ) from airborne snow radar

Contemporary climate warming over the Arctic is accelerating mass loss from the Greenland Ice Sheet through increasing surface melt, emphasizing the need to closely monitor its surface mass balance in order to improve sea-level rise predictions. Snow accumulation is the largest component of the ice sheet’s surface mass balance, but in situ observations thereof are inherently sparse and models are difficult to evaluate at large scales. Here, we quantify recent Greenland accumulation rates using ultra-wideband (2–6.5 GHz) airborne snow radar data collected as part of NASA’s Operation IceBridge between 2009 and 2012. We use a semiautomated method to trace the observed radiostratigraphy and then derive annual net accumulation rates for 2009–2012. The uncertainty in these radar-derived accumulation rates is on average 14 %. A comparison of the radarderived accumulation rates and contemporaneous ice cores shows that snow radar captures both the annual and longterm mean accumulation rate accurately. A comparison with outputs from a regional climate model (MAR) shows that this model matches radar-derived accumulation rates in the ice sheet interior but produces higher values over southeastern Greenland. Our results demonstrate that snow radar can efficiently and accurately map patterns of snow accumulation across an ice sheet and that it is valuable for evaluating the accuracy of surface mass balance models.


Introduction
In the past 2 decades, climate warming over the Greenland Ice Sheet (GrIS) has accelerated its mass loss, nearly quadrupling from ∼ 55 Gt a −1 between 1993 and 1999 (Krabill et al., 2004) to ∼ 210 Gt a −1 , equivalent to ∼ 0.6 mm a −1 of sea-level rise, between 2003(Shepherd et al. 2012. As GrIS mass loss has accelerated, a fundamental change in the dominant mass loss process has occurred (e.g. Tedesco et al., 2015). It switched from ice dynamics to surface mass balance (SMB) processes, which include accumulation and runoff (van den Broeke, 2009;Enderlin et al., 2014). This recent shift emphasizes the need to monitor SMB which, over most of the GrIS, is dominated by net accumulation.
Here, we use the complete set of airborne snow radar data collected by NASA's Operation IceBridge (OIB) over the GrIS from 2009 to 2012 to produce net annual accumulation rates, hereafter called accumulation rates for simplicity, along those flight lines. The radar-derived accumulation rates are compared to both in situ data and model outputs from the Modèle Atmosphérique Régional (MAR).

Background
In situ accumulation-rate measurements are limited in number by the time and cost of acquiring ice cores, digging snow pits or monitoring stake measurements across large sectors of the ice sheet. Only two major accumulation-rate measurement campaigns have been undertaken across the GrIS: the first in the 1950s when the US Army collected pit data along long traverse routes (Benson, 1962) and the second in the 1990s when the Program on Arctic and Regional Climate Assessment (PARCA) collected an extensively distributed set of ice cores (e.g. Mosley-Thompson et al., 2001). A recent traverse and study by Hawley et al. (2014) reports a 10 % increase in accumulation rate since the 1950s and highlights the need to monitor how Greenland precipitation is evolving in the midst of ongoing climate change. Although many other accumulation-rate measurements exist, they are more limited in either space or time (e.g. Dibb and Fahnestock, 2004;Hawley et al., 2014).
To date there is no annually resolved satellite-retrieval algorithm for accumulation rate across ice sheets. Hence, the two primary methods used to generate large-scale (hundreds of kilometers) accumulation-rate patterns are model outputs and radar-derived accumulation rates . High-resolution, near-surface radar data have shown good fidelity at mapping spatial patterns of accumulation over ice sheets at decadal and annual resolutions from both airborne and ground-based radars (Kanagaratnam et al., 2001(Kanagaratnam et al., , 2004Spikes et al., 2004;Arcone et al., 2005;Anshütz et al., 2008;Müller et al., 2010;Medley et al., 2013;Hawley et al., 2006;2014;de la Peña et al., 2010;Miège et al., 2013). Radars detect the lateral persistence of isochronal layers within the firn. When these layers are either (1) dated in conjunction with ice cores or (2) annually resolved from the surface, they can be used to determine along-track accumulation rates.
Early studies by Spikes et al. (2004) in Antarctica and Kanagaratnam et al., (2001 and in Greenland used high/very-high-frequency (100 to 1000 MHz) ground-based and airborne radars, with vertical resolutions of ∼ 30 cm, to measure decadal-scale accumulation rates between dated ice cores. These high/very-high-frequency radars can penetrate hundreds of meters in the dry-snow zone and tens of meters in the ablation zone (Kanagaratnam et al., 2004). Subsequent studies utilized the larger bandwidths of ultra/super-highfrequency (2 to 20 GHz), frequency-modulated, continuouswave (FMCW) radars with centimeter-scale vertical resolutions capable of mapping annual layers within ice sheets (e.g. Legarsky 1999;Marshall and Koh, 2008;Medley et al., 2013). Ultra/super-high-frequency radars can penetrate tens of meters in the dry-snow zone and meters in the ablation zone. Legarsky (1999) was among the first to show that such radars could image annual layers, and Hawley et al. (2006) further demonstrated that a 13.2 GHz (Ku-band) airborne radar imaged annual layers in the dry-snow zone of the GrIS to depths of up to 12 m.
Most previous studies used radar data that overlapped spatially with ice cores or snow pits for both dating layers and density information. Medley et al. (2013) and Das et al. (2015) showed that accumulation rates could also be derived using density from a regional ice-core ensemble. Density end members are used to derive uncertainty limits, and the derived regional density profile is sufficient for radar studies of accumulation and SMB (Das et al, 2015). Additionally, Medley et al. (2013) showed that snow radar is capable of resolving annual layer in high accumulation regions where such layers were well preserved. Therefore, it was possible to date the layers by simply counting from the surface downwards.
Regional climate models (RCMs), general circulation models (GCMs) and reanalysis products provide the only spatially and temporally extensive estimates of accumulation-rate fields at ice sheet scales (e.g. Burgess et al., 2010;Hanna et al., 2011;Ettema et al., 2009;Fettweis, 2007;Cullather et al., 2014). In a comprehensive model intercomparison study, Vernon et al. (2013) found that modeled accumulation rates had the least spread across the RCMs considered but still had a ∼ 20 % variance. Chen et al. (2011) found the range in mean accumulation rate across the GrIS between five reanalysis models to be ∼ 15 to 30 cm a −1 , while Cullather and Bosilovich (2012) found the range in mean accumulation rate across the GrIS between reanalysis data and RCMs to be ∼ 34 to 42 cm a −1 . While these models continue to improve, there is clearly a continuing need for large-scale accumulation-rate measurements to evaluate their outputs.
3 Data, instruments and model description

Modeled accumulation rates and density
Accumulation rate and snow/firn density profiles were derived from the MAR RCM (v3.5.2; X. Fettweis, personal communication, 2015). MAR is a coupled surfaceatmosphere model that simulates fluxes of mass and energy in the atmosphere and between the atmosphere and the surface in three dimensions and is forced at the lateral boundaries with climate reanalysis outputs (Gallée, 1997;Gallée and Schayes, 1994;Lefebre et al., 2003). It incorporates the atmospheric model of Gallée and Schayes (1994) and the Soil Ice Snow Vegetation Atmosphere Transfer scheme (SISVAT) land surface model, which includes the multi-layer Crocus snow model of Brun et al. (1992). The MAR v3.5.2 simulation used here has a horizontal resolution of 25 km and utilizes outputs from the European Center for Medium Range Weather Forecasting (ECMWF) ERA-Interim global atmospheric reanalysis at the lateral boundaries (Dee et al., 2011). Additional details are described by Fettweis (2007), with updates described by Fettweis et al. (2011Fettweis et al. ( , 2013 and Alexander et al. (2014). MAR has been validated with in situ data and remote sensing data over the GrIS, including data from weather stations (e.g. Lefebre et al., 2003;Fettweis et al., 2011), in situ and remotely sensed albedo data (Alexander et al., 2014) and ice-core accumulation rates (Colgan et al., 2015), and it has been used to model both past and future SMB (Fettweis et al., 2005. We use accumulation rates and density profiles simulated by MAR for the period during which the radar data were collected (2009 to 2012). In MAR, the initial falling snow density (ρ s,0 ) is parameterized as a function of the temperature in the first model layer (T air ) in • C (at roughly 3 m above the surface) and 10 m wind speed (V ) in m s −1 . The parameterization differs depending on atmospheric temperature as follows.
If T air is greater than −5 • C, then If T air is less than −5 • C and V > 6 m s −1 , the parameterization of Kotlyakov (1961) is used: If V < 6 m s −1 the initial snow density is set to the fixed value of 200 kg m −3 .
After snow falls to the surface, snow compaction in MAR is described according to the scheme of Brun et al. (1989), where the compaction of a layer (dδz/dt) of thickness δz is given by where ρ is the dry-snow density (g cm −3 ), T s is the snow temperature (degrees Celsius) of the layer, σ is the vertical stress from the snow above (kg m −1 s −1 ) and C is a function of snow grain size and snowpack liquid water content.

In situ density and accumulation-rate data
The SUrface Mass balance and snow depth on sea ice working group (SUMup) dataset (July 2015 release) compiles publicly available accumulation-rate, snow depth and density measurements over both sea ice and ice sheets (Koenig et al., 2013). We use two subsets of these data. First, to characterize density across the GrIS, we extract the snow/firn density measurements ranging in depth from the snow surface to 15 m (the depth to which MAR predicts firn density), which contains over 1500 measurements from snow pits and cores at 62 sites. At each site, the number of measurements ranges in number between 8 and 170 and maximum depths range from 1 to 15 m. (Koenig et al., 2014;Koenig and SUMup, 2015;Miège et al., 2013;Mosley-Thompson et al., 2001;Hawley et al., 2014;Baker, 2012) (Fig. 1). Second, to compare radar-derived and in situ accumulation rates, we consider only accumulation-rate measurements within 5 km of OIB snow radar data, a criterion that includes 11 cores from the SUMup dataset (Mosley- Thompson et al., 2001). To expand this comparison, an additional dataset of 71 cores was included (J. McConnell, personal communication, 2015), providing 23 additional cores within 5 km of OIB snow radar data ( Fig. 1).

Determining the density profile and uncertainties
Because we seek to derive accumulation rates from nearsurface radars across large portions of the ice sheet, we require firn density profiles that cover the entire GrIS. Modeled snow/firn density profiles from MAR were investigated for use. However, a preliminary comparison of the SUMupmeasured density profiles to MAR-estimated density profiles showed that MAR simulated density values in the top 1 m of snow/firn were lower (0.280 ± 0.040 g cm −3 ) than observed (0.338 ± 0.039 g cm −3 ) (Fig. 2). The comparison of measured and modeled density was simultaneous in time. Specifically, the MAR density profile output on the day of the measurement was used in this comparison. We consider it beyond the scope of this study to investigate and explain why MAR underestimates near-surface density. Here we assume that the firn density in the top 1 m is 0.338 g cm −3 . Below 1 m, the model and observed densities are similar (4 % mean difference with the model generally overestimating measured density slightly), so the spatially varying modeled density profiles are used for 30 April of each year. Hence, a hybrid measured-modeled density profile is used to determine accumulation rates from the snow radar data (Fig. 2). Our assigned uncertainty in the top meter is the relative standard deviation in observed density (12 %), which we assume is due to the natural variability in surface density. This uncertainty is higher than the assumed mean measurement uncertainty of 2-5 % (Proksch et al., 2016) and smaller than the mean difference between the modeled and observed values within the top meter (16 %). No spatial bias is evident between the mean model used here in the top 1 m of snow/firn and the observed density.

Deriving accumulation rates from snow radar and uncertainties
The radar travel time is converted to depth (z) using the snow/firn density profile and the dielectric mixing model of Looyenga (1965). Errors in radar-derived depth come from two primary sources: (1) the dielectric mixing model chosen and (2) layer picking. The choice of the dielectric mixing model maximizes potential error at a density of ∼ 0.300 g cm −3 . The maximum possible difference in depth over 15 m is 3 % assuming a constant density of 0.320 g cm −3 and < 1 % assuming a constant density of 0.600 g cm −3 (Wiesmann and Matzler, 1999;Gubler and Hiller, 1984;Schneebeli et al., 1998;Looyenga, 1965;Tiuri et al., 1984). The second source of error occurs during manual adjustment of the picked layers (Sect. 4.3.4) and is estimated to be a maximum of ±3 range bins, or ∼ 8 cm. Given the relative standard deviation in accumulation rate, the range bin error contributes a mean uncertainty of 7 % with a range of 4 to 24 %. Lower accumulation rates have a higher relative error from layer picking. The water-equivalent accumulation rateḃ in m w.e. a −1 at along-track location where z is the depth of layer in m, and ρ is average snow/firn density to depth z in kg m −3 . Hence, the numerator is the mass in kg m −2 to depth z, a is age of the layer in years from the date of radar data collection and ρ w is the density of water in kg m −3 (e.g. Medley et al., 2013;Das et al., 2015). Depth z is calculated using the radar two-way travel time (TWT), the snow/firn density (ρ) and the Looyenga (1965) dielectric mixing relationship as follows: where TWT is the travel time to the dated layer in sec, c is the speed of light in m s −1 , ρ i is ice density in kg m −3 and ε i is the dielectric permittivity of pure ice. Combining these two equations giveṡ The cumulative mean snow/firn density (ρ) is determined by the density profile described in Sect. 4.1. The layers are picked in the radar data using a semiautomated approach described in Sect. 4.3. Layer ages are determined by assuming spatially continuous layers are annually resolved and dated accordingly from the year the radar data were collected. The radar data were collected during springtime (April-May) and the surface is assumed to be 30 April to align with the modeled accumulation rate, which was processed to monthly values. Subsurface picked layers are assumed to be 1 July ± 1 month, so the first layer represents 10 months and each subsequent layer is 12 months. Peaks in radar reflectivity are, assuming ice with no impurities, caused by the largest change in snow density. In the ablation and percolation zone, the peak in density difference occurs in the summer between the snow layer and ice or the snow/firn layer and the high-density melt/crust layer, respectively (e.g. Nghiem et al., 2005). In the dry-snow zone, the peak density contrast also occurs in the summer between the summer hoar layer and the denser snow/firn layer (e.g. Alley et al., 1990). While melt/crust and hoar layers can form at other times, it is assumed they will be smaller density contrasts and, therefore, cause a smaller radar reflection than the dominant layers which occur near 1 July. We assume picking errors are randomly distributed in space due to the heterogeneous nature of snow/firn. Density peaks have been shown to vary in depth along ice-core transects, likely due to smallscale microstructure differences (e.g. Machguth et al., 2016).
To calculate the total uncertainty on the radar-derived accumulation rate, the largest error is assumed for density (12 %) and age (10 %) and the error for the mean accumulation rate is assumed for layer picking (7 %). Equation (3) shows that the density profile is used for calculating both depth and water equivalent. The derivative of Eq. (3) determines the correlated error between depth and density. Assuming uncorrelated and normally distributed errors between density, age and layer picking the mean accumulation-rate uncertainty is 14 %, with a range of 13 % for the highest accumulation rates and 27 % for the lowest accumulation rates. This relative uncertainty is similar to previous studies by Medley et al. (2013) and Das et al. (2015) for radar-derived accumulation rates.

Semiautomated radar layer picker
A semiautomated layer detection algorithm is developed to process the large amounts of OIB snow radar data (> 10 4 km a −1 ), analogous to the challenges faced by Mac-Gregor et al. (2015) for analysis of very-high-frequency radar sounder data. While a fully automated method is ultimately desirable, we have found that it is necessary to manually check every automated pick, making adjustments as needed by an experienced analyst, to distinguish between spatially discontinuous radar reflections, caused by the natural heterogeneity of firn microstructure, and spatially consistent annual layers. Our algorithm processes the OIB snow radar data in four steps outlined below.

Surface alignment
The snow surface is detected by a threshold, set to 4 times the mean radar return from air, which is assumed to be the radar background noise level. A median filter is applied vertically to each radar trace to minimize noise. Any surface detection that is displaced by greater than 10 range bins (∼ 25 cm) from its adjacent traces is not used and that entire vertical trace is ignored in subsequent analysis. Data arrays are then aligned to the surface and truncated above and below the surface (200 and 800 range bins, respectively), equivalent to ∼ 25 m into the snow/firn, to reduce data volumes. Layer depths are measured relative to the snow surface. The radar data are then horizontally averaged (

Layer detection
The algorithm takes advantage of the difference between high-frequency and low-frequency spatial variability in the travel time/depth domain to identify peaks in returned power in the radar data. Such peaks are formed by the stratified accumulation layers of interest in this study, and they form over small spatial scales, equivalent to high frequency, in the travel time/depth domain. Our peak detection process is thus a type of high-pass filter, resulting in the set of disjointed points detected at radar reflection peaks in the time domain and in adjacent traces along the flight path. These points are connected into continuous layer segments using the halfmaximum width of each peak's waveform (Fig. 3, locations of radargrams shown in Fig. 1).

Layer indexing
Each along-track detected layer is indexed, with both a number and the corresponding year, counting down from the surface detection (Fig. 3). This indexing begins with the seg- Figure 3. Example snow radar radargrams from 2011 in the percolation zone (top), inland from Jakobshavn Isbrae, and dry-snow zone (bottom), near the ice divide ∼ 220 km south of Summit Station, showing automatically picked layers (black) resulting from the layer picking algorithm before any manual adjustments. Indexing by year is shown at the left end of each picked layer. Snow radar data frames represented are 20110422_01_218 to 20110422_01_244 (top) and 20110426_03_155 to 20110426_03_180 (bottom) (Leuschen, 2014). Locations of the radargrams are shown by the red lines in Fig. 1. mentation of the layers, so that each layer is uniquely identified with a layer number. The peak points within each segment are connected by spline fitting, resulting in a set of sharply defined along-track layers at different depths (Fig. 3). These layers represent 1 July in the appropriate year, counting from the surface and the year collected.

Manual adjustment with the layer editor
A graphical user interface (GUI) was developed to verify automated layer detections by displaying the snow radar radargram and the resulting automated-layer detections. An analyst uses the GUI to quickly visually compare the picked layers and the radargram. The GUI application allows for layer editing as needed including tools for layers or parts of layers to be added, deleted, gap-filled and re-indexed. The GUI accelerates layer picking by providing the ability to scroll through all the radargrams and picked layers, including the previous and subsequent along-track data, to detect errors. Scrolling allows for spatially continuous layers, which may not be datable at all locations, to be propagated and dated from a location where annually resolved layers are evident from the surface. Error statistics for the automatic algorithm were not kept but depend generally on the quality of the radar data, influenced by both radar and aircraft operations, and the regional characteristics of the snow/firn microstructure, which can either preserve or erode layering.

Radar-derived accumulation rates
Annual radar-derived accumulation rates and their uncertainties were calculated for all 2009-2012 OIB radar data that contained detected layers (Fig. 4). The increase in coverage from 2009 to 2012 is related to an increasing number of OIB flights over the GrIS and adjustments to the snow radar antenna and operations that improved overall data quality. These accumulation-rate patterns are consistent with observed and modeled large-scale spatial patterns for the GrIS: high accumulation rates in the southeast coastal sector and lower accumulation rates in the northeast (Fig. 5). Year-toyear variability in the accumulation rate is also evident, even at the ice sheet scale; e.g., in the southeast accumulation rates were lower in 2010 than in 2011.
The radar-derived accumulation rate in Fig. 4 represents only the first layer detected by the snow radar, or approximately the annual accumulation rate from the year prior to data collection. For simplicity, we refer to this quantity as the annual accumulation rate, but we caution that it does not strictly represent the calendar year. The values shown in Fig. 4 represent only 10 months of accumulation, based on our assumption that the radar layers date to 1 July (Sect. 4.2) and that the data collection date is 30 April for all OIB data, which may differ from the actual flight date by up to a month. When comparing the first layer of radar-derived accumulation to modeled estimates from MAR (Fig. 5) or other accumulation measurements, these timing differences must be considered. Although the first layer represents only a partial year, all deeper layers represent a full year, from 1 July to 30 June. We simultaneously compare the time represented by the layer to MAR estimates of accumulation. Figure 6 shows the number of detected layers, or previous years, discernable in the OIB radar data. For the majority of the GrIS, one to three annual layers are discernable. OIB flight lines are clustered in the ablation/percolation zones of the GrIS, where radar penetration depths are reduced by the increased density, englacial water and layering structure of the firn column (Fig. 3). In the GrIS interior, where dry-snow conditions allow deeper radar penetration, annual layering going back over 2 decades is detectable (Fig. 3).  Crossover points were assessed to determine the internal consistency of the radar-derived accumulation rates (Figs. 8 and 9). While no consistent spatial pattern is found in the crossover errors, the largest discrepancies were found in 2011 and 2012 in the northwest and southeast (Fig. 8). Other inconsistencies are likely due to snow deposition occurring between flights in the southeast and incorrectly picked layers that were either sub-or multiannual in the northwest. Figure 9 shows a scatter plot of crossover points. There are relatively few outliers, and those that are outlying are generally offset by a factor of 2, suggesting an error in layer detection/dating rather than a radar-system error. Crossover differences per year, including the mean, standard deviation and maximum, are given in Table 1. These differences are comparable (mean of 0.04 m w.e. a −1 or 4 range bins) to our inferred relative uncertainty of 14 %, emphasizing the overall validity of our methodology.

Comparison with modeled accumulation
The along-track radar-derived accumulation rates were gridded to the MAR grid for comparison. The mean-local, radarderived accumulation rate was used when gridding. Because OIB flight lines are not spatially homogenous, each MAR grid cell represents a different number of radar-derived values, so grid cells are not sampled equally. With this discrepancy noted, this gridding method is still the most straightforward approach for this comparison. Figure 10 shows the difference between the radar-derived and MAR accumulation rates. The mean difference for all years is low (0.02 m w.e. a −1 ). Table 1 shows the annual variability of the mean difference, which is low for every year except 2010, when large differences are seen over the southeastern coastal region of the GrIS (Fig. 10). Figure 10 shows that MAR generally reconstructs accumulation rates well in the GrIS interior (consistent with the comparison with ice-core estimates presented by Colgan et www.the-cryosphere.net/10/1739/2016/ The Cryosphere, 10, 1739-1752, 2016 Table 1. Radar-derived accumulation-rate crossover analysis. Columns include the year the radar data were collected and the number, the mean, the standard deviation and the maximum difference of radar-derived accumulation at crossover points. Minimum crossover values were 0 for all years. The final column shows the mean difference between the gridded radar-derived accumulation and the MAR estimates of accumulation from 1 July to 30 April. Year al., 2015) but has larger differences around the periphery, especially in the southeast and northwest in particular years. In the southeast, MAR generally has higher accumulation rates, except in 2011 when there is a mixed pattern of agreement and higher accumulation rates. Higher values in the southeast are not surprising and are likely due to the large changes in surface topography that are not resolved by the relatively coarse model grid (Burgess et al., 2010). In 2011, the northwestern coastal region of the GrIS was well sampled by OIB and MAR has lower accumulation rates there in contrast to 2010, when the area was sampled and had higher values. The origin of this anomaly in the northwest is less clear, but it may be related to forcing at the lateral boundaries of MAR that may not capture a relatively small storm track into this region. Figure 11 shows a scatter plot of the radar-derived and MAR-estimated accumulation rates. These values are not well correlated (Pearson correlation coefficient r 2 = 0.2) and have large root-mean-square errors (RMSEs) (0.24 m. w.e. a −1 ), emphasizing that further improvements in accumulation-rate modeling and measurements are needed, particularly over the southeastern and northwestern GrIS.

Comparison with annually resolved in situ data
Between 2009 and 2012, OIB flew within 5 km of 34 core locations but only two locations, NEEM and Camp Century (Fig. 1), were coincident in time with the layers we detected. Each of these locations has two cores, providing annual accumulation rates and a measure of spatial variabil- ity. Figure 12 compares the radar-derived to core measured accumulation rates. At NEEM, the two ice cores and radar data are nearly co-located, within 0.6 km of each other. The radar-derived accumulation rates are self-consistent between 2011 and 2012 and agree well with the ice cores (RMSE of 0.06 m w.e. a −1 ). For comparison, the two NEEM cores have an RMSE of 0.05 m w.e. a −1 for the period of overlap. A timing discrepancy arises with this comparison because the ice cores, with higher dating resolution from isotopic and chemical analysis, are dated and reported for calendar years, whereas the radar-derived accumulation is assumed 1 July-30 June (Sect. 4.2). This mismatch in the measurement is likely evident in Fig. 12 by the differences in the annual peaks between the cores and radar-derived accumulation having similar means yet differing magnitudes from year to year.
Near Camp Century, the cores and radar data are farther apart from each other. The radar data are located within 4.4 km of the Camp Century core and the GITS core is located ∼ 8.2 km from the Camp Century core. These separations are likely responsible for the poorer agreement at this  Figure 8 shows the spatial distribution of these crossover errors in (m w.e. a −1 ).
site of radar-derived accumulation rate to the Camp Century core (RMSE 0.10 m w.e. a −1 ) and the larger difference (RMSE 0.07 m w.e. a −1 ) in accumulation rate between the two cores for the period of overlap. At Camp Century, and throughout much of northern Greenland, two older, continuous layers were dated from the interior of the ice sheet and spatially traced. These layers, dated July 2000 and July 2001, could not be dated with the Camp Century data alone -hence the temporal gaps in annual accumulation rate at this location. While it is more difficult to analyze the results at Camp Century, with only three overlapping points and no continuous annual time series of radar-derived accumulation rates, our estimates are within the expected variability and capture the long-term mean value.

Discussion
This study is the first to derive annual accumulation rates from near-surface airborne radar data collected across large portions of the GrIS. The pattern of radar-derived accumulation rates compares well with known large-scale patterns and clearly shows that these accumulation-rate measurements have the potential for evaluating model estimates. At the two locations with contemporaneous cores, radar-derived rates agree well with the long-term mean. Additional cores, with direct overflights, are clearly needed to continue assessing the accuracy of the radar-derived accumulation rates.
The work shown here only incorporates layering detected in the radar data that is annual and continuously dated from the surface to depth at some location. We did not exhaustively trace all layering detected by the snow radar, i.e., there are still contiguous layers, not connected to a dated layer, in the dataset that were not utilized. For example, in the centralnorthern GrIS, there is a strongly reflecting layer varying be- tween 15 and 18 m that cannot be dated with the radar data alone. If ice cores were drilled to identify the age of this layer, techniques similar to those developed by MacGregor et al. (2015) or Das et al. (2015) could be used to determine multiannual accumulation rates in additional regions of the GrIS and extend the snow radar record. Further deconvolution processing of these data, currently ongoing, will likely also resolve additional deeper layers in the snow radar data.
Annual radar-derived accumulation rates are not extrapolated spatially here due to their sparseness relative to the scale of the entire ice sheet. Such extrapolation between flight lines, which vary from year to year, must be left for future work, as additional data are collected and fill in gaps.
The largest overall discrepancy between radar-derived and MAR estimates of accumulation is in 2010. In 2010 MAR has higher accumulation rates over most of the GrIS and particularly over the southeastern GrIS (Fig. 10). A previous study (Burgess et al., 2010) showed that modeling accumulation rate is difficult in this region. However, the discrepancy is also due, at least in part, to the fact that in 2010 there is a higher percentage of radar data collected over the lower Figure 11. Comparison between radar-derived and MAR-estimated accumulation rate (m w.e. a −1 ). Radar-derived accumulations (Fig. 4) were averaged within each MAR grid cell. Figure 10 shows the spatial distribution of the differences.
portions of the southeastern GrIS compared to other regions. This spatial sampling bias of OIB flight lines is amplifying the discrepancy in 2010. Because OIB data are not spatially consistent from year to year, caution must be used when extrapolating to ice sheet scales.
In 2011 MAR has lower accumulation rates over the northwestern GrIS in a region just to the south of Camp Century in contrast to higher values in 2010. This small region is known to receive more snowfall locally than the surrounding areas, because storms on Greenland's western coast are diverted as the land mass to the north protrudes farther west into Baffin Bay (K. Steffen, personal communication, 2015). MAR does show increased accumulation in this region (Fig. 5) but differs in magnitude from the radar-derived measurements in 2010 or 2011. It is possible that MAR is not reproducing this local pattern because it is close to MAR's lateral boundaries, where the coarser GCM may not adequately represent this phenomena. This discrepancy emphasizes the importance of understanding the possible effects of lateral forcing of RCMs on accumulation-rate fields and warrants further study.
Finally, the uncertainties in the radar-derived accumulation rate are approximately equally distributed between the layer picking, age and density. However, the layer picking is likely overestimated and in most cases likely much lower, leaving age and density uncertainties nearly equal (Medley et al., 2013). Age uncertainties could be better constrained with a better understanding of the timing of density peaks across the ice sheet. Our assumption that the surface date is 30 April could be adjusted to the flight date if the modeled accumulation rates were reprocessed to daily values.
With respect to density uncertainty, we assumed a constant and uniform density in the top meter of snow/firn as modeled outputs did not match measured values (Fig. 2). This assumption could lead to spatial bias in our analysis if regional density deviates significantly from the mean, though existing measurements do not show any clear evidence of such spatial bias. Spatially distributed density measurements and improved density models spanning the entire firn column are required to take full advantage of the layering detected by near-surface radars and to reduce errors in radarderived accumulation rates. The current sampling of in situ measurements has large spatial gaps over the southwestern, northern and northeastern GrIS and the majority of the measurements are located in the upper-percolation and dry-snow zones (Fig. 1). To further constrain and improve the density models required for radar-derived accumulation rates, these gaps must be filled. Additionally, the snow radar's signal penetration around the perimeter of the GrIS is relatively shal-low, resolving one to three annual layers only, with the majority of detected layers in the top meter of snow/firn (Figs. 6 and 7). Accumulation rates are calculated using measurement averages in this section of the snow/firn column, likely causing less error than the MAR-modeled density. Improvement to modeled near-surface density should be considered for improved snow radar analysis.

Conclusions
A semiautomated method was developed to process tens of thousands of kilometers of airborne snow radar data collected by OIB across the GrIS between 2009 and 2012. The resulting radar-derived accumulation-rate dataset represents the largest validation dataset for recent annual accu-mulation rates across the GrIS to date. This dataset captures the large-scale accumulation-rate patterns of the GrIS well. Over 2 decades of annual radiostratigraphy is observed in the dry-snow zone, near Summit Station, and 1 to 3 years are generally detectable in the ablation/percolation zones. Our estimated uncertainty in the radar-derived accumulation is 14 %. This study emphasizes the need for ice cores coincident in time with airborne overflights and, more importantly, for improved density profiles, particularly in the top 1 m of snow/firn. These radar-derived accumulation rates should be used to evaluate RCM/GCM and reanalysis products, as demonstrated here using the MAR model. MAR matches the radar-derived accumulation rates well for most of the interior of the GrIS but tends to have higher accumulation rates in the southeastern coastal region of the GrIS and, in at least 1 year, has lower accumulation rates in the northwestern coastal region of the GrIS. While determining the precise nature of these differences is left for future work, we have clearly demonstrated the usefulness of the ice-sheet-wide, radar-derived accumulation-rate datasets for improving SMB estimates. As the GrIS continues to lose mass through SMB processes, monitoring accumulation rates directly is vital.