1. Introduction

Radio-echo sounding (RES) is one of the most powerful tools for retrieving information on the subsurface properties of polar ice sheets. The technique has been used for several decades to determine the basal topography and ice thickness of Antarctica and Greenland (e.g. Reference DrewryDrewry, 1983; Reference Bamber, Layberry and GogineniBamber and others, 2001). However, the radar data also contain information about the englacial properties of the ice sheets, often manifested as internal layers (or internal radar reflectors) caused by changes in the dielectric properties of the ice. Studies have shown that the changes are induced by variations in density, impurity content, acidity or crystal fabric orientation. The latter often occurs in response to overburden ice pressure and therefore does not necessarily imply a continuous layer of the same age (Reference Eisen, Hamann, Kipfstuhl, Steinhage and WilhelmsEisen and others, 2007). Layers caused by variations in impurity content and acidity, however, are generally agreed to be isochrones (e.g. Reference SiegertSiegert, 1999; Reference Miners, Wolff, Moore, Jacobel and HempelMiners and others, 2002), and their stratigraphy provides information on the past mass-balance rate and ice flow dynamics (e.g. Reference Fahnestock, Abdalati, Joughin, Brozena and GogineniFahnestock and others, 2001a). The layers are thus a valuable tool for reconstructing the past of an ice sheet (e.g. Reference Hindmarsh, Leysinger Vieli and ParreninHindmarsh and others, 2009; Reference Leysinger Vieli, Hindmarsh, Siegert and BoLeysinger Vieli and others, 2011).

The Center for Remote Sensing of Ice Sheets (CReSIS) at the University of Kansas has conducted surveys of the Greenland ice sheet since 1993 (e.g. Reference GogineniGogineni and others, 2001), and this comprehensive RES dataset contains both ice thickness information and a wealth of other information, including internal layering. Comparisons with ice-core records obtained at the ice divide (e.g. the Greenland Ice Core Project (GRIP) and the North Greenland Ice Core Project (NorthGRIP) (GRIP members, 1993; NorthGRIP Members, 2004)) allow for dating of the internal layering. This paper presents one of the first attempts at mapping the internal stratigraphy of the Greenland ice sheet, specifically focusing on the climate transition at the end of the last glacial period 14.7 ka b2k (before year AD 2000, according to the GICC05 timescale with a dating error of ±0.3 ka (Reference RasmussenRasmussen and others, 2006)).

Sections 2 and 3 of this paper contain a description of the data and the method used for retrieving the depth of the transition. Section 4 presents the results, and Section 5 discusses the implications for reconstruction of past mass balance based on the transition depth.

2. Data

The CReSIS radars used for data collection operated at a centre frequency of 150 MHz and a bandwidth of 17 MHz with a peak transmit power of 200W from 1993 to 2008. In 2010 and 2011, they operated at a centre frequency of 195 MHz and a bandwidth of 30 MHz with a peak transmit power of ∼400W. Post-processing renders high-accuracy GPS positions, and the errors in the radar-derived ice thicknesses are typically within 10m in central Greenland. However, errors may be substantially larger close to the margin. More information on the radar system is provided by Reference Gogineni, Chuah, Allen, Jezek and MooreGogineni and others (1998, Reference Gogineni2001).

During extensive field campaigns in Greenland from 1993 to 2001, CReSIS acquired RES data along the ice divide intersecting several deep ice-core drilling sites (GRIP, NorthGRIP and NEEM (North Greenland Eemian Ice Drilling)). The internal layering in the RES data can thus be correlated with ice-core records and dated using ice-core stratigraphy. The reconstructed age-depth relationships show that a change in reflectivity often corresponds to a transition in climate (Fig. 1a); for example, the locations of several deep layers correlate with interstadials identified in the NorthGRIP core (NorthGRIP Members, 2004). One such very distinct change is the onset of the warm Bølling-Allerød interstadial at 14.7 ka b2k characterized by a warmer climate with higher accumulation than the glacial period preceding it (Reference Johnsen, Dahl-Jensen, Dansgaard and GundestrupJohnsen and others, 1995). At the NorthGRIP drill site this transition was found at 1655 m depth (Reference SteffensenSteffensen and others, 2008).

Fig. 1. (a) Example of a radargram (Z-scope) from central Greenland. Notice the upper part has an overall higher background reflectivity than the deeper part. The Bølling–Allerød transition that corresponds to a change in the radargram from higher to lower background reflectivity is indicated with a red line. The age–depth scale is based on the NorthGRIP ice core (NorthGRIP Members, 2004). (b) Radar acquisition lines over Greenland as black lines overlain on surface elevation (Reference Bamber, Layberry and GogineniBamber and others, 2001). The location of the radargram in (a) is marked with a thick, white line, and the NorthGRIP drill site with a blue dot.

3. Method

The method in this study aims at automating the processing and tracing of internal layers as much as possible. The method is based on the observation that the majority of the CReSIS RES data from central Greenland display a characteristic appearance: in their upper half, the data contain numerous internal layers and have a higher background reflectivity than the lower, older part, where only a few visible layers can be seen. This is clear both from the radargrams (Fig. 1a) and in the individual RES records (A-scopes; Fig. 2b). Compared to the depth–age relationship from the NorthGRIP ice core, this change in the RES data coincides with the transition from the glacial period to Bølling–Allerød. In other words, the upper part is mainly ice from the Younger Dryas and the Holocene while the deeper ice with lower background reflectivity is from the last glacial period.

Fig. 2. (a) Illustration of a ramp function (Eqn (1)). (b) A RES data record (A-scope) from Figure 1 with the transition between interglacial and glacial ice marked with a dashed, black box. (c) The 2-D version of Eqn (1). (d) Zoom-in on the area in the dashed, black box from (b) and all the RES records 250 m horizontally on each side of it. (e) Excerpt of radargram from Figure 1. (f) Histogram of the relative received power values from the radargram in (e). (g) Reconstructed radargram from (f) with a black arrow indicating the first estimate of the transition depth.

Previous studies have attempted to locate transitions in time series automatically. One approach is described by Reference MudelseeMudelsee (2000) using a ramp function (Fig. 2a) defined

(1)

as and is based on the assumption that the standard deviation of the time-series data also takes the shape of a ramp function. The ramp function is then fitted to the climate data using a search function to minimize the misfit. The RES data take a very similar shape when considering the individual RES records (Fig. 2b). However, the method proposed by Reference MudelseeMudelsee (2000) presents some problems for the RES data at hand, primarily for two reasons: (1) It requires a manual selection of the size of the smoothing window, which is not feasible for the RES data due to the vast data volume; and (2) by considering only the individual RES records, information about the horizontal coherence between each RES record is easily lost. Thus, while the algorithm proposed by Reference MudelseeMudelsee (2000) is unsuitable in the context of RES data, the idea of fitting a ramp function to the RES records is sound. The method presented here uses a two-dimensional (2-D) version of Eqn (1) as shown in Figure 2c (compare the similarity in shape with the RES data shown in Fig. 2d). An overview of the processing algorithm is presented below, and a more detailed description is provided in the Appendix. The results will be made available on the Centre for Ice and Climate’s website (www.iceandclimate.nbi.ku.dk).

3.1. Initial identification of the transition

The radargram is first horizontally smoothed (incoherently averaged) over a window size of 20 records. For the datasets acquired after 2002, samples are spaced 3 m apart; these datasets are resampled to a horizontal resolution of ∼50m. In older datasets, sample spacing is 130-150 m, so no resampling is performed. The radargram is then normalized with respect to ice thickness using

(2)

where is the normalized depth, z is the depth, and s(x), b(x) are the elevations of the surface and bed, respectively. Equation (2) is equivalent to setting the surface equal to 0 and the bed equal to 1. The normalization serves to negate any large-scale fluctuations of the layer stratigraphy caused by an undulating bedrock. The upper 500 m of the radargram is discarded to avoid additional noise from the firn layer, since this layer is often not resolved very well in the depth-sounding radar. This is especially true for the older RES data, because these data were collected with an analogue sensitivity time control that was set to prevent receiver saturation by strong surface echoes.

The radargram is then subsampled into data slices of a few kilometres width (Fig. 2e), and the histogram of the relative received power for each subsample is calculated (Fig. 2f). The purpose of the histogram is to sort the data points into bins, in this case seven bins, and the mean value of each bin is then used to construct a new representation of the radargram using only those seven mean values (Fig. 2g). In most cases, the interglacial/glacial transition will be located at the change between the lowest and second-lowest mean data values (black arrow in Fig. 2g) because of the decreasing relative received power with depth and the overall low data values in the glacial data. The number of bins will partly determine the location of the estimate, so the method is repeated with three, five and ten bins. Tests showed that using more bins leads to overestimation of the transition depth, while using fewer than three bins does not capture the change in relative received power at the transition. This process is described in more detail in the Appendix.

4. Results

In total, 939 radar flight-lines acquired between 1993 and 2011 over central and northern Greenland were processed. Of these, the transition was not visible or its presence was inconclusive in 301 cases (approximately one-third), and these datasets were discarded. While we cannot completely exclude that poor radar sensitivity may have caused the absence of clear layering, the majority of the discarded flight-lines were from areas of high surface velocity and/or shallow ice thickness. It therefore seems likely that the lack of a visible transition can mainly be attributed to either disruption of the layering by ice flow, as has been observed elsewhere (e.g. Reference Siegert, Payne and JoughinSiegert and others, 2003), or simply to the fact that layering cannot be discerned in thin ice (Reference Fahnestock, Abdalati, Joughin, Brozena and GogineniFahnestock and others, 2001a). From the remaining 638 lines, the success of the algorithm was assessed by visual inspection. It was found that

The obtained transition depths from >400 flight-lines are linearly interpolated onto a 5 km grid (Fig. 3a) that also imposes a degree of smoothness on the resulting map. The normalized depths are shown with a value of 0 for the ice surface. Blue indicates that the transition is situated in the upper part of the ice, and red indicates that it is in the deeper part of the ice.

Fig. 3. A gridded map of the normalized depth of the transition in colours, and the coast of Greenland outlined as a thick grey line. White contours outline the transition depth from 0.4 to 0.85 with intervals of 0.05. All data points >100km from a flight-line have been masked out. (a) The 50ma–¹ surface velocity contour from Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others (2010) is shown as a thick black line, and radar flight-lines as thin grey lines. (b) Black contours show the surface mass balance (Reference EttemaEttema and others, 2009). The red and blue arrows indicate areas where the general trend of increasing transition depth as the ice thickness decreases is less clear.

The most distinct feature in Figure 3a is seen in the central part of North Greenland, where the depth of the transition is at a minimum. At this minimum, the transition is only slightly more than one-third of the way down into the ice column, and its location does not coincide with the location of the present ice divide which is at least 100 km farther east. West of the depth minimum there is a distinct gradient in the transition depth where the contours follow the ice divide. Closer to the margin, the transition is very deep down in the ice column: over two-thirds of the way down in places. While the general trend is one of uniformly increasing transition depth as ice thickness decreases, there are areas where this relationship is less clear, for example the area east of the GRIP drill site and on the northeast coast close to 80° N (red and blue arrows in Fig. 3b). Here the change in transition depth appears not to follow the decrease in ice thickness as closely.

The normalized depth of the transition based on the smoothed, gridded data is 0.566 for GRIP and 0.527 at NorthGRIP. High-resolution layer counting at the drill sites has determined that the onset of the Bølling period is located at 0.558 for GRIP (Reference Seierstad, Johnsen, Vinther and OlsenSeierstad and others, 2005) and at 0.52 for NorthGRIP (Reference RasmussenRasmussen and others, 2006). Considering that the results presented here show the transition depth on a low-resolution, smoothed grid, we find that the agreement is reasonable.

Discussion

The depth of an isochrone is influenced by the surface accumulation, basal melting, spatial gradients in ice flow velocity, and ice-sheet evolution. These processes work on different spatial and temporal scales; as such, it cannot be expected that all the processes that have influenced the isochrone stratigraphy are clearly reflected in the results presented here.

To the first order, the surface mass balance controls the vertical velocity of ice particles. Thus, the depth of the transition is expected to be correlated with the current surface mass balance, provided the surface mass balance has not undergone large changes since the deposition of the transition. The correlation is expected to be positive since higher surface accumulation leads to faster burial rates and thereby greater transition depths. For the data at hand, the correlation coefficient between the transition depth and the surface mass balance is 0.8, and Figure 3b makes it clear that there is an overall agreement between the two. The increase in transition depth along the western coast is well correlated with areas of high accumulation, while the low accumulation values on the northwest coast fit with comparatively lower transition depths. The spatial gradient in accumulation west of the ice divide mirrors the changes in transition depth very well. In other words, these findings indicate that the large-scale accumulation pattern has not undergone substantial changes since the onset of the Holocene. The part of North Greenland that is situated east of the ice divide appears to have been, on average, a low-accumulation area throughout the Holocene, while the west coast appears to have been, on average, a high-accumulation area throughout the Holocene. Our findings agree with the results of Reference Fahnestock, Abdalati, Joughin, Brozena and GogineniFahnestock and others (2001a) who focused on a single flight-line from 1999 but traced several internal layers in that flight-line. Reference Fahnestock, Abdalati, Joughin, Brozena and GogineniFahnestock and others (2001a) found that the accumulation reconstructed from the internal layers is consistent with present-day values.

A second factor that controls isochrone depth is basal melting which is often observed as a dragging down of internal layers (e.g. Reference Fahnestock, Abdalati, Luo and GogineniFahnestock and others, 2001b; Reference Buchardt and Dahl-JensenBuchardt and Dahl-Jensen, 2007). Ice-core drilling (e.g. NorthGRIP Members, 2004), as well as properties of the bedrock from RES reflections (Reference Oswald and GogineniOswald and Gogineni, 2012), have shown that a significant part of the base of the ice sheet in North Greenland is at the pressure-melting point. In the results presented here, the influence of basal melting is obscured by the smoothing of the gridded data. Localized areas of increased basal melt were often observed in the radargrams, but the basal melting usually only affected a comparatively small area downstream of the basal melt spot. Modelling studies have shown that depending on the basal melt rate, surface accumulation and ice velocity, the isochrones recover their elevation as they move away from the basal melt spot (e.g. Reference Leysinger Vieli, Hindmarsh and SiegertLeysinger Vieli and others, 2007); this is also what we observe. We therefore suggest that future studies use only the transition depths from the original transects rather than the gridded versions shown in Figure 3.

While the surface mass balance explains the large-scale pattern of the transition depth, there are several features that might be influenced by other processes. In some areas, drawdown of the transition apparently occurs in areas of fast ice flow (Fig. 3a), for example in the area east of the GRIP drill site (red arrow in Fig. 3b). It also seems likely that the presence of the large ‘Northeast Greenland Ice Stream’ would influence the depth of the transition and potentially impact the extent of the observed transition high. Thus, the fact that the contours of the transition depth do not precisely match those of the low-accumulation area might be attributed to its proximity to the Northeast Greenland Ice Stream. However, the transition depth might also have been influenced by the changing surface elevation of the Greenland ice sheet since the end of the last glacial period (Reference VintherVinther and others, 2009). It is unclear how the change in surface elevation would affect the depth of the transition since the change might have been accompanied by changes in accumulation, thinning rates and ice flow dynamics.

Figure 3 shows that less than half of the ice column in the central part of North Greenland is Holocene ice; therefore, a substantial amount of older, glacial ice must be present in this region. This contradicts results of some modelling studies that suggest that this part of Greenland might have been ice-free during the Eemian (e.g. Reference HuybrechtsHuybrechts, 2002; Reference Born and NisanciogluBorn and Nisancioglu, 2012). If North Greenland was ice-free during the Eemian, then either a substantial build-up of ice has occurred between the Eemian and the end of the last glacial period, or there is a significant volume of accreted ice pushing the layers upwards. This phenomenon has been observed in Antarctica (Reference BellBell and others, 2011), but the timescales under which the accretion operates are uncertain at best.

Conclusion

The transition between ice from the Bølling–Allerød period and ice from the glacial period is easily identifiable in a large portion of the RES data collected over northern and central Greenland. This study presents an automated method for retrieving the depth of the transition. While the method successfully identifies the transition in a large portion of the data, it requires a final visual assessment to ensure the accuracy of the data presented here.

Results reveal that the transition is located relatively high up in the ice in the central part of North Greenland, with increasing depths towards the margins. The transition depth mirrors the present-day surface accumulation pattern with some exceptions, which are thought to be caused by high-velocity areas dragging the layers down. Furthermore, the results indicate that either a substantial amount of ice from the last glacial period is present in the central part of North Greenland or that large rates of basal ice accretion are taking place over a substantial area of North Greenland.

The depth of the transition between glacial ice and Bølling–Allerød will assist future modelling efforts in understanding the evolution of the Greenland ice sheet. The data presented here will be made available to the wider glaciological community on the website of the Centre for Ice and Climate (www.iceandclimate.nbi.ku.dk) or upon request to the lead author.