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.
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).
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.
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
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.
3.2. Brute force fitting
Based on the predicted transition depths, some parts of the radargram can now be disregarded and only the part where the transition is likely to be present is considered. The 2-D version of Eqn (1) is now fitted to the data using the following assumptions:
The value of the lower level of the ramp function (x 2 in Eqn (1)) is equal to the lowest mean value from the histogram constructed in the first processing step described above.
The value of the upper level of the ramp function (x 1 in Eqn (1)) is set equal to the average of t 1 and the 20 points above it.
The transition is sharp (t 1 and t 2 are closely spaced).
The best fit is found using a brute force search by constructing a ramp function with a horizontal width of 1000 m and moving it vertically up point-for-point while calculating the misfit between the ramp function and the RES data. Essentially, the search function is looking for the best fit of the shape shown in Figure 2c with the data shown in Figure 2d. Once the best fit is found, the fitting algorithm will move 250 m forward horizontally and start a new search using a ramp function of 250 m width. However, this new search is limited in that it can only move vertically five samples up or down from its previous position. This restriction is based on the observation that on small horizontal scales the layers are unlikely to move up or down substantially. Once the fitting algorithm reaches the end of the radargram it restarts at its original point and moves backwards to the start of the radargram. This process is described in more detail in the Appendix.
3.3. Selection of final estimate
A set of solutions will now exist (due to the different starting points). The final result is chosen to be the transition depth that is most frequently selected by the fitting routine. This final estimate is saved with the other estimates, and the next radar line is then processed.
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
In 5 1% of the data (35% of the full dataset), the transition was identified correctly as the final estimate.
In a further 29% (20% of the full dataset), the transition was identified but the final estimate was incorrect (see the description of the processing algorithm above). This was adjusted manually.
For the remaining 20% of the data (14% of the full dataset), the algorithm was unable to pick out the correct transition, mainly due to a very gradual increase in received relative power from interglacial to glacial ice. These datasets were discarded
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.
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.
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.
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.
N.B.K. is supported by European Research Council grant No. 246815 ‘Water Under the Ice’. The Centre for Ice and Climate is funded by the Danish National Research Foundation. We acknowledge the use of data and data products from CReSIS generated with support from US National Science Foundation grant ANT-0424589 and NASA grant NNX10AT68G, and the work of many people – too many to name individually – in collecting and processing the radar data. We thank J. MacGregor and an anonymous reviewer whose insightful comments helped improve the manuscript.
In the following, the method for locating the transition depth is described in more detail.
The radargram is smoothed, resampled if necessary, normalized and divided into subsamples as described in Section 3. The histogram of the relative received power for each radargram subsample is then calculated using, for example, seven bins as shown in Figure 2f. This returns seven mean bin values of relative received power, one from each bin. Each individual data point is then assigned the mean value that is closest to its original value. This leads to the construction of a radargram that only has seven unique relative received power values. This radargram is then simplified further, leading to the radargram shown in Figure 2g. This additional simplification is performed by assigning all points located at the same depth the value that most frequently occurs at that depth. From the simplified radargram shown in Figure 2g, the initial transition depth is estimated. This depth is assumed to be at the change from the lowest to the second-lowest value that occurs the highest up in the simplified radargram (black arrow in Fig. 2g).
The process described above is repeated for all sub-samples of a radargram, and the depth estimates are saved. The estimates provide the basis for discarding the parts of the radargram that are unlikely to contain the transition. In general, anything higher in the ice than the highest depth estimate and anything lower than the lowest estimate ±0.05 in normalized depth is discarded, although this range is adjusted slightly depending on the estimated depth range (e.g. if the depth estimate is so low that an additional 0.05 would set it below bedrock, the range is adjusted upwards).
Using the depth estimates from the previous processing step, parts of the radargram can now be disregarded and only the part where the transition is likely to be present is considered. The fitting algorithm is initiated at the start of every subsample (i.e. every few kilometres) using a ramp function with a horizontal width of 1000m and following the processing steps below:
1. Define t 2 as the deepest point in the (reduced) radargram.
2. Define the location of t 1 as five vertical sample points above t 2 .
3. The value of x 2 is set equal to the lowest mean bin value from the histogram constructed in step 1.
4. x 1 is set equal to the average of the relative received power in t 1 and the 20 points above in the vertical direction.
5. Calculate the misfit between the constructed ramp function and the RES data.
6. Move t2 one point up and repeat the process.
The misfits are saved and the values of t 2 and t 1 that return the lowest misfit are defined as possible transition depths. The fitting algorithm now moves 250 m forward horizontally and the fitting process restarts, only this time the search is limited vertically to only consider five samples up or down from its previous position and the width of the ramp function is 250 m. Once the fitting algorithm reaches the end of the entire radargram, it restarts at its original point and moves backwards to the start of the radargram.
Since the fitting processing step is applied at the start of every subsample, a set of solutions will now exist. The final transition depth estimate is chosen to be the transition depth that is most frequently selected by the fitting routine. This final estimate is then saved with the other estimates and the next radargram is processed.
To ensure that the transition depth returned by the algorithm is correct, every radargram is inspected visually. The final results presented in this paper have been generated following these guidelines:
In the cases where the transition is easily identifiable by a user and the algorithm predicts the correct depth, the transition data are saved (this returned >300 flight-lines).
In the cases where the transition is easily identifiable by a user but the algorithm chooses an incorrect depth, all the transition depth estimates are loaded. Frequently one of the estimates from the fitting algorithm will be the correct depth. If the radargram intersects another radargram with a correctly predicted transition depth, the intersection point is used as an additional checkpoint.
If none of the transition estimates is correct or if the transition is not easily identifiable by a user, the radargram is discarded.