Introduction
The growth and decay of ice sheets is driven by a balance between accumulation of snow on the surface, primarily in the high-elevation interiors, and the melting, runoff, evaporation, sublimation and iceberg calving that takes place primarily along the lower-elevation margins. The mass balance of the Greenland ice sheet, in particular, is of increasing importance to scientists and policymakers, as rising air and ocean temperatures have increased the rate of melting and the velocity of calving outlet glaciers, contributing to a rising sea level (Reference Lemke and SolomonLemke and others, 2007). The net mass balance of the Greenland ice sheet has been the subject of much recent work (e.g. Reference Alley, Spencer and AnandakrishnanAlley and others, 2007; Reference Rignot, Box, Burgess and HannaRignot and others, 2008; Reference Wouters, Chambers and SchramaWouters and others, 2008; Reference ZwallyZwally and others, 2011). Additionally, investigators have studied the surface mass balance, primarily accumulation, using point measurements (Reference BensonBenson, 1962; Reference Ohmura and ReehOhmura and Reeh, 1991; Reference Dibb and FahnestockDibb and Fahnestock, 2004; Reference BalesBales and others, 2009), automated weather stations (e.g. GC-Net: http://cires.colorado.edu/science/groups/steffen/gcnet/), remote sensing (Reference Drinkwater, Long and BinghamDrinkwater and others, 2001; Reference Kanagaratnam, Gogineni, Gundestrup and LarsenKanagaratnam and others, 2001; Reference Hawley, Morris, Cullen, Nixdorf, Shepherd and WinghamHawley and others, 2006), and more recently, using climate models constrained by observations (Reference Box, Bromwich and BaiBox and others, 2004, Reference Box2013; Reference BalesBales and others, 2009; Reference BurgessBurgess and others, 2010; Reference HannaHanna and others, 2011). Spatially extensive ground-based profiles of accumulation rates, however, are rare (Reference BensonBenson, 1962). The spatial variability of accumulation on varying length scales is difficult to determine with point measurements or climate models, especially when local topography is considered.
Here we derive spatially extensive measurements of accumulation along a 1009 km traverse from ~1700 m a.s.l. in northwest Greenland to the ~3200 m a.s.l. center of the ice sheet. We trace internal reflecting horizons (IRHs) in 400 MHz ground-penetrating radar (GPR) data. We use density profiles collected along the traverse to calculate electromagnetic wave speed through the snow and firn in order to determine the depth of the IRHs, and date the IRHs using snow/firn chemistry data from cores collected at each end of the traverse.
The Greenland Inland Traverse
During the spring of 2011, a science field team from Dartmouth College joined the Greenland Inland Traverse (GrIT), an effort funded by the US National Science Foundation to take heavy and bulky resupply items to Summit station, at the center of the Greenland ice sheet. The GrIT 2011 convoy consisted of three Case tractors, each hauling cargo, and a single Pisten Bully tractor hauling the science team and their equipment. Along the traverse route (Fig. 1), the science team collected continuous geodetic-quality GPS profiles and GPR data. At selected points along the route, the team made point measurements: snow pits (five pits, <2 m depth) for physical and chemical strati-graphic analyses, and shallow cores (six cores, <10 m depth) for density profiling. The locations for these stops were a compromise between high sample frequency (every ~40 km) along the route and the expeditious transport of cargo to Summit. Thus the science team sampled as frequently as feasible without impeding progress. The traverse team spent 24 days traveling from Thule to Summit, hauling cargo and making science stops, and 11 days traversing back to Thule. While GPS and GPR profiles were collected in both directions, here we report only the profiles collected from Thule to Summit. We use point measurements of density from both legs.
Methods
Positioning using GPS
We collected differential GPS data using a Trimble NetR8 reference receiver, with a Zephyr Geodetic antenna mounted to the roof of the Pisten Bully. In addition to our geodetic-grade GPS receiver, we also operated a navigation-grade Garmin GPS receiver, recording positions along the traverse route as a regular part of navigation.
We processed the GPS data using TRACK, the kinematic processing module of the GAMIT/GLOBK suite (Reference King and BockKing and Bock, 2009). Estimated root-mean-square horizontal values were generally between 13 and 18 cm, largely due to long baselines. Occasionally there were gaps in the TRACK output, due to extensive phase ambiguities, multipath signals or unknown causes. If these gaps spanned <5 min, we inserted navigation-grade GPS data recorded from the Garmin GPS. If these gaps spanned >5 min, we discarded the corresponding GPR data.
Ground-penetrating radar
We collected GPR data using a Geophysical Survey Systems Inc. (GSSI) SIR-3000 radar unit with a 400 MHz antenna. The choice of frequency was a trade-off between range resolution in discriminating IRHs and total range; lower frequencies can image deeper IRHs at the expense of lower depth resolution. The antenna was towed on the snow surface in a small plastic sled, affixed to a wooden ‘outrigger’ that held the sled out of the impressions left by the tracks and kept the antenna ~3 m away from the metal of the vehicle. We recorded 2048 samples per trace over a range window of 600 ns. At a relative permittivity of ~2.4, typical of dry-firn conditions, the range was ~58 m. The 400 MHz short-pulse radar we used has a range resolution (ability to resolve distinct features) of ~0.35 m in firn (Reference Spikes, Hamilton, Arcone, Kaspari and MayewskiSpikes and others, 2004). We recorded 16 traces per second, which at the Pisten Bully’s average travel speed of ~3.2 m s-1 results in 5 traces recorded per meter. Note that this spacing between traces varies with the vehicle speed.
We computed the position of each trace using our differential GPS measurements. To co-register GPR and GPS data, we used time stamps embedded in the two data streams. To reduce data volume and increase the signal-to-noise ratio, we ‘stacked’ the radar data by averaging adjacent traces together. Rather than choosing a number of traces to stack, we stacked using a fixed distance, based on our GPS measurements. We stacked each 20 m stretch of traces together, generally resulting in a stack of between 70 and 120 traces. The final radargram of the entire GrIT route is shown in Figure 2. We did not apply any time-variable gain or filtering.
Snow/firn density profiles
We measured density profiles along the traverse using snow-pit and core samples. In each snow pit, we used a cutter of a known volume to excise a snow sample, and weighed the sample. Using the cutter we measured density at 0.1 m resolution throughout the depth of each snow pit (~1 m). At many locations we were able to extend the density profile by drilling shallow cores for density only. We determined densities from these cores by measuring the diameter, length and mass of 10–20 cm core sections. Finally, we measured densities from cores extracted on the traverse for chemical analysis, as well as a published density profile from Summit (Reference Hawley, Morris and McConnellHawley and others, 2008).
Travel-time to depth conversion
In order to estimate a density profile for each radar trace to determine electromagnetic wave velocity, and thus depth, we interpolate linearly along constant depths between density profiles. At any given depth, we interpolate between all of the profiles that exist at that depth. Thus our interpolation becomes spatially coarser as fewer points reached the greatest depth. Our scheme is limited by the depth of the shallower (‘2Barrel’) of the two end profiles, and beyond this depth we use the density from the deeper Summit core. The final density array we use is shown in Figure 3. We use these densities to convert from two-way travel time to depth.
The electromagnetic wave speed, ν (m s-1), is related to the dielectric permittivity, εr (dimensionless), and the speed of light in a vacuum, c (~3 × 108 m s–1), by
The relative permittivity of snow and firn is related to density, ρ (g cm-3), following Reference Kovacs, Gow and MoreyKovacs and others, (1995):
We identified the direct-coupling wave in the radar profiles, and set the depth to zero at that point. We then calculated the depth for each subsequent GPR sample, consecutively, using the velocity profile of the firn above that sample. In this way we avoided using an average velocity for depth determination.
Core depth/age scales
We collected three cores for glaciochemical analysis: one close to each end of the traverse route, and one along the route (Fig. 1). The ‘Owen’ core was drilled to 100 m depth at Summit in summer 2010. The ‘2Barrel’ core was drilled to 20 m depth near Camp Century (76.935358 N, 63.153868 W) in summer 2011. The ‘Galen’ core was drilled to 10 m depth at 74.422338 N, 39.294338 W at ~800 km along the traverse route, also in 2011. To account for the year between cores, we adjusted the Owen depth using a published mean accumulation rate from Summit (Reference Dibb and FahnestockDibb and Fahnestock, 2004). We processed and sampled all three cores in the Dartmouth ice-core chemistry laboratory, using a continuous melter system with discrete sampling (cf. Reference Osterberg, Handley, Sneed, Mayewski and KreutzOsterberg and others, 2006). The melter system collected samples continuously with a sampling interval of 4–10 cm, producing 6–20 samples each year for all three records. We measured stable water isotope ratios (δD and δ18O) using a Picarro L2310-i isotope analyzer with a typical analytical uncertainty of ±0.4% for δD and ±0.2% for δ18O. We measured major and trace element sample concentrations (Na, Mg, Ca, Al, Fe, S, K, Ti, Sr, Cd, Cs, Ba, Pb, Bi, U, As, V, Cr, Mn, Co, Cu, Zn, La, Ce, Pr) by inductively coupled plasma mass spectrometry (ICP-MS) using a Thermo Element2 at the Climate Change Institute, University of Maine, Orono, USA (Reference Osterberg, Handley, Sneed, Mayewski and KreutzOsterberg and others, 2006). The seasonal variations in stable water isotopes, sea-salt Na and dust (represented by many trace elements, including Al, Fe, Ti, U and rare earth elements) are well defined and can be used to annually date the cores (Reference McConnellMcConnell and others, 2000; Reference Mosley-ThompsonMosley-Thompson and others, 2001). We used multiple chemical time series to determine annual demarcations in each core record, and different analytes were favored for timescale development at each site, depending on local site characteristics. For example, stable water isotopes were strongly diffused at the low-accumulation Galen core site (0.12 m w.e. a-1), so spring dust peaks and winter sea-salt Na peaks were preferentially used to develop the timescale. In contrast, stable water isotopes and dust were preferred for the 2Barrel depth/age scale, because meltwater percolation at the site more strongly eluted salts than dust particles (Reference Wong, Hawley, Lutz and OsterbergWong and others, 2013). The uncertainty of the timescales is estimated to be ±0.5 years, based on repeated independent annual-layer counts.
Identifying internal reflecting horizons
In the radargram of Figure 2, many identifiable IRHs are visible. We manually selected clear, strong IRHs that could be consistently traced for the full length of the traverse route. For the purposes of our manual picking, we chose to follow local maxima in the radar return. Where a layer appeared to bifurcate (occasionally this happens when accumulation increases and existing layers become discernible within the ability of our radar), we continued the layers based on the dip and amplitude of surrounding layers (cf. Reference Arcone, Spikes and HamiltonArcone and others, 2005). If this was not possible, or if the layer was discontinuous for more than several hundred meters along-track, we discarded the layer.
Estimating uncertainty
Uncertainty in our accumulation rates can arise from errors in correctly locating layers, errors in dating the layers and/or errors in the density used for converting to water equivalent. It is important to note that these error sources are independent of one another.
To assess error in layer picking, we note that our GPR settings result in a ~0.03 m spacing between radar samples. Close inspection of the layers reveals that peaks defining IRHs are clear within ±2 samples, i.e. within ±0.06 m. With our average epoch length of 2.4 years, at a density of 0.5 g cm-3, this will introduce an error in accumulation of 0.0125 m a-1.
The estimated uncertainty in dating is 0.5 year. At the lowest-accumulation core (Galen) the smallest distance between layers we used was 0.26 m w.e. over an epoch of 2.2 years. This gives an uncertainty in accumulation due to dating of ~±0.02 m a-1.
Reference KarlöfKarlöf and others, (2005) estimated the error in measuring density using similar techniques as 1.4%. If this is an optimistic estimate and our measurements resulted in twice as much error, the resulting 2.8%, even at higher-accumulation sites, results in errors in accumulation of ~±0.014 m a-1. As this scales with accumulation, error in accumulation measurement due to uncertainty in density is much lower over most of the traverse.
These errors are all random, not systematic, and thus are not additive. It is then most reasonable to average the errors together; thus we estimate an uncertainty in our accumulation rates due to all causes at ~0.015 m a–1. Because these errors are random and not systematic, they are unlikely to contribute to a bias.
Results
The depth/age scales from the cores allow us to demonstrate the isochronous nature of the IRHs. The depth and age of the IRHs are summarized in Table 1; the IRHs we have delineated are clearly isochrones, as ages of the IRHs are nearly identical at all core sites. We take the date of any given traced IRH as the date of the Owen core from the depth at which the layer intersects the core. The final dated layers with depth and age can be seen in Figure 4.
To calculate an accumulation rate between two dated IRHs at each 20 m stacked GPR trace, we subtract the depth of one layer from the other to determine a firn thickness. We then convert the firn thickness to water equivalent by multiplying the thickness by the average density of the firn between the two IRHs. We then divide by the difference in the ages of the two IRHs to determine an average accumulation rate for the period between the two IRHs. Our accumulation rate estimates are presented in Figure 5.
Several spatial features are evident in these accumulation profiles. The broad trend in accumulation is higher values near the coast giving way to lower accumulation inland as the traverse route gains altitude; this trend is consistent with orographic lifting causing the most precipitation closer to the coast. Accumulation rates then rise again just before approaching Summit, as the traverse curves back south.
On a smaller scale, the effect of small-scale (5–10 km wavelength) local topography is evident in the measured accumulation rates; as noted by Reference Black and BuddBlack and Budd, (1964) in the Antarctic, local surface curvature can drive short-scale accumulation patterns, in which concave regions tend to receive greater accumulation and convex areas receive lower accumulation. This effect can be seen at ~820 km in both Figures 4 and 6, with the variability increasing with depth in Figure 4, as expected. Reference MiègeMiège and others, (2013) noted similar accumulation features in radar measurements collected in southeast Greenland.
Discussion
Comparison with modeled accumulation rates
Reference BurgessBurgess and others, (2010) used firn cores and coastal precipitation records to calibrate the precipitation output of the Polar MM5 numerical regional climate model, producing annual modeled accumulation grids for Greenland for the years 1958–2008, accessible for download at http://bprc.osu.edu/Greenlandaccumulation/. We use nearest-neighbor interpolation along the traverse to determine a MM5 modeled accumulation rate for each 20 m GPR trace for which we have measured accumulation, and summarize the comparison in Figure 6.
On average, differences are small (median –0: 012 m; mean 0.001 m). However, biases are highly correlated into two regions. In the lower-accumulation interior of the ice sheet (~250–1000 km), at higher elevations, MM5 typically overestimates accumulation rates slightly (by 0.014 m), while in the higher-accumulation, lower-elevation marginal area near the coast, MM5 typically underestimates average accumulation by 0.0386–0.0466 m w.e. (median and mean, respectively). As is evident in the measured accumulation rates, this is a region of high spatial variability. The spatial accumulation maps of Reference BurgessBurgess and others, (2010) show a high-accumulation region in this area, probably related to a combination of factors: southerly flow from Baffin Bay during storms combined with orographic lifting combines with the close proximity to the North Water Polynya (the largest polynya in the world, 8: 5 × 104 km2), increasing moisture content of these storms in winter. As this northern end of the traverse route runs almost parallel to the accumulation contours, subtle variations in the route cause large changes in measured accumulation.
It is worth noting that the region that MM5 underestimates by up to 25%, within 200 km of the coast, is the region that shows the most change since 1955 (Fig. 9, below), and also the area with the highest accumulation. Thus, accumulation estimates from MM5 do not capture the increase we see since the 1950s. It is possible that this is a resolution issue in a region of high accumulation gradient, or it may be that MM5 does not capture the other possible climate changes contributing to the increase. The Wilcoxon signed-rank test finds these differences to be significant (p = 2: 2 × 10–16). This test was chosen over a t –test because the Shapiro–Wilk normality test showed both datasets to be non-normally distributed (p –valueMM5 = 2: 2 × 10-16; p –valueGPR = 2: 2 × 10-16), as are their differences (p –valueGPR MM5 = 2: 2 × 10-16).
At a shorter length scale (~10 km), the measured accumulation rates show subtle spatial changes related to surface topography, as noted above. Figure 7 illustrates the relationship between accumulation and small-scale topography along a 30 km segment of the traverse close to Summit. The highs in accumulation (e.g. 845 km in Fig. 7) correspond to local dips in topography, creating concave areas that accumulate more snow. The modeled accumulation rate does not show these small-scale effects, as its 24 km gridscale is independent of small-scale surface topography.
To illustrate temporal changes in accumulation over our measurement period, we plot time series of accumulation near each end of the traverse in Figure 8. At the northern, higher-accumulation, lower-elevation end of the traverse, though there is an offset in average accumulation rate (as seen in Fig. 6), the changes in accumulation rate correlate closely (r 2 = 0: 5158, p = 0: 0128), indicating that the calibrated MM5 model captures the interannual trends in accumulation well. At the southern, lower-accumulation, higher-elevation end of the traverse, the interannual changes are more subtle, and the correlation between measured and modeled series is lower and not statistically significant (r 2 = 0: 2610, p = 0: 1083).
Comparison with 1950s in situ measurements
Reference BensonBenson, (1962) measured accumulation rates in a series of snow pits along a similar route to ours in 1952–55 (Fig. 1). Our modern measurements thus allow us to investigate changes in accumulation rates over the 52 years between traverses. The different routes taken by the two traverses (Fig. 1) require that we disentangle spatial changes in accumulation from temporal changes. In order to compare the two datasets, we corrected the accumulation rates of Reference BensonBenson, (1962) using the spatial pattern of accumulation for 1958 (Fig. 1), the earliest year of MM5 output presented by Reference BurgessBurgess and others, (2010). For each accumulation location presented by Reference BensonBenson, (1962), we calculated the closest point along the GrIT. We then interpolated the 1958 Polar MM5 accumulation rates between the two points (Benson’s and GrIT’s) to find a spatial correction value to be added to the Reference BensonBenson, (1962) data for comparison. It is possible that one or more of Benson’s snow pits were located in areas of topographical roughness, which can affect accumulation rates (Fig. 7), and it is important to note that we cannot account for this effect.
The adjusted accumulations of Reference BensonBenson, (1962) are shown in Figure 9, together with our 10 year average accumulation rate and 1σ range from 1997–2007. A subtle but persistent positive difference is evident along the traverse (Fig. 9). On average the GPR-derived modern rate was 0.022–0.031 m w.e. a-1 greater than the Benson values, based on the mean and median differences, respectively. The one-sided non-parametric Wilcoxon signed-rank test shows that the positive difference between paired observations along the traverse is statistically significant (p = 0. 0000076). This test was chosen over a t –test because both datasets were non-normally distributed, based on the Shapiro–Wilk normality test (p –valueBenson = 0: 001455; p –valueGPR = 0: 000223), as were their differences (p –valueGPR-Benson = 0. 0179). It is important to consider the possibility of differing measurement techniques contributing to this result. Reference BensonBenson, (1962) used 500 cm3 SIPRE tubes to measure density and reports his overall accuracy to be ±0.005 g cm3 or ~1.5% at the densities he encountered. Thus it is unlikely that technique contributes to the difference in measured accumulation rates.
At the lower-elevation, coastal end of the traverse, the difference between Reference BensonBenson, (1962) and GrIT accumulation values is more pronounced. This could be the result of decreasing sea-ice extent in Baffin Bay (Reference RaynerRayner and others, 2003), which would increase the moisture available to storms in this region, and potentially increase the strength of these storms.
The 0.022 m w.e. a-1 increase we find implies a ~10% change in average accumulation in the dry-snow zone of the ice sheet in the past 52 years. The significance of this result is that increased accumulation, as seen on the East Antarctic ice sheet by Reference Davis, Li, McConnell, Frey and HannaDavis and others, (2005), driven by the increased ability of warmer air to hold moisture, is also evident inland on the Greenland ice sheet. Reference BoxBox and others, (2013) found a ~1.2% decade-1 increase over the 20th century, and our result is consistent with this, if slightly higher. The significance of this for the mass balance of the Greenland ice sheet, is that at least some of the increased mass loss from melting at the lower-elevation margins of the ice sheet is balanced by the small increases in mass gain from increased accumulation in the higher-elevation interior.
Conclusions
We have traced internal reflecting horizons in 400 MHz radar data along ~1000 km of the GrIT. Using densities measured along the traverse, we determined the depth of these reflecting horizons. We have shown these horizons to be isochrones using depth/age scales at each end of the traverse, with the ages of the layers at each end agreeing within an average of 2.5 months.
Using the isochrones and density measurements, we determined accumulation rates along the traverse for the period 1980–2007. These accumulation rates differ slightly from those modeled by Reference BurgessBurgess and others, (2010), both spatially and temporally over the coincident period. Though the modeled accumulations accurately reflect the large-scale patterns in accumulation, they understandably fail to capture small-spatial-scale (tens of kilometers) accumulation anomalies (~10%) driven by surface topography.
Comparison of our measured accumulation rates with those measured in the 1950s by Reference BensonBenson, (1962) indicates a ~2% decade-1 increase in accumulation between the periods 1945–55 and 1997–2007. Thus, due to a warmer atmosphere driving an increased capacity for moisture, and in common with the findings of Reference Davis, Li, McConnell, Frey and HannaDavis and others, (2005) in East Antarctica, accumulation in the interior of the Greenland ice sheet has increased slightly in the currently warming climate.
Acknowledgements
This project was supported by the US National Science Foundation under grant NSF-OPP 0909265. The Owen core was drilled by Terry Gacke of the Ice Drilling and Design Office (IDDO), University of Wisconsin–Madison, USA. Logistical support was provided by Ch2M HILL Polar Services. Mike Handley was instrumental in processing our samples at the University of Maine. This work would not have been possible without the support of the traverse ‘swing’ crew, for which we are grateful. The authors thank and commend Reference BurgessBurgess and others, (2010) for making their results freely available. We also thank Clément Miège and an anonymous reviewer whose comments improved the manuscript. Finally, we thank Scientific Editor John Woodward and Chief Editor Jo Jacka for their handling of the manuscript.