Layer tracing and dating

The methodology for tracing layers is similar to work previously presented by Reference Fahnestock, Abdalati, Luo and GogineniFahnestock and others (2001b). A point of the layer to be traced is selected from an echogram (Fig. 2). A weighted window centred on this point is compared with a windowed portion of the neighbouring vertical profile. By identifying the maximum cross-correlation between the two windowed profiles, the layer being traced is advanced. Figure 2 shows a sample echogram from a portion of a flight-line, highlighted in Figure1, with a traced layer overlain in black. Also shown (in white) is the return from the bedrock. The layer in one echogram is extended to other echograms by matching layers at flight-line crossover points.

Fig. 2. Example echogram through the Greenland ice sheet. This ∼161km west-to-east section is part of the flight-line shown in bold in Figure 1. The distances are from the start of the western end of that flight-line. The bedrock return is highlighted by the white line, and the 6956 year layer is shown by the black line. This isochrone is located (in this section) at around 750 m below the ice surface. The depths below the ice surface were calculated assuming a constant refractive index of n = 1.78 for ice (see text for more details).

Figure 1 shows all the flight-lines flown during 1997, 1998 and1999. The flight-lines used in this study are all from surveys flown in May 1999. The location of the traced layers is shown in black. The spatial extent of these layers is less than that of the flight-lines themselves. This is due to gaps in the echograms and to regions of faded or noisy IRHs which make the tracing of particular layers impossible. In some cases, it is possible to intervene manually in order to jump over such regions and to restart a layer, but this is subjective and care needs to be exercised. Several jumps have been made to extend the dataset.

The layer tracing was done using the binned two-way travel-time data. To convert this travel time to a depth, a constant velocity of 1.68610₈ ms–¹ (Gogineni and others, 2001) was used for the radar velocity through ice. This results in each bin corresponding to 4.494m. The two layers used were selected for their spatial coverage rather than their relation to any specific time horizon. The depths of these two layers, at the Greenland Icecore Project (GRIP) core site, were found to be 783.5±4.5 and 1203.6±4.5m. The error estimate assumes that errors result from the layer-tracing technique only, and that 4.494m is the depth resolution of the radar, as flight-lines were never more than one pixel out at crossover points. In fact, the rms depth difference at crossover points was 0.16 m. We use the larger error estimate, however, as area-sonable estimate of other sources of error such as the vertical variation of ice density which is discussed later. The depths are converted to ages by linear interpolation of the published GRIP age/depth data (NSIDC, 1997). The shallower layer was dated to 3891±31BP and the deeper layer to 6956±42BP.

Palaeo-accumulation rates

A common approach to calculating accumulation rates from IRHs is to use so-called “sandwich” models such as those of Reference NyeNye (1957) and Reference Dansgaard and JohnsenDansgaard and Johnsen (1969). These models assume that the layer depth at a given location has resulted from a constant ice thickness and accumulation rate. In other words, the ice thickness and accumulation rate along a flow-line have remained unchanged.

If a velocity field is available, it is then possible to move beyond these basic “sandwich”models. By explicitly following a flowline, the ice thickness and accumulation rates can vary. In order to calculate flowlines, a spatially extensive velocity field is required. Balance velocities are the only practical means of providing such a synoptic dataset for the whole ice sheet.

Balance velocities are an estimate of the depth-averaged velocity, assuming that the ice sheet is in a state of balance. Comparison between this steady-state “balance” field and observed or modelled velocities allows insight into the dynamics of an ice sheet (e.g. Reference Budd and WarnerBudd and Warner, 1996; Reference Joughin, Fahnestock, Ekholm and KwokJoughin and others, 1997; Reference Bamber, Hardy, Huybrechts and JoughinBamber and others, 2000; Reference Huybrechts, Steinhage, Wilhelms and BamberHuybrechts and others, 2000).

Three datasets are required to calculate balance velocity (see, e.g., Reference Budd and WarnerBudd and Warner, 1996): ice surface slope, ice thickness and ice accumulation rate. Ice surface slope is calculated from a digital elevation model (DEM) of the Greenland ice sheet (Reference Bamber, Ekholm and KrabillBamber and others, 2001a), and the ice thickness used is that described by Reference Bamber, Layberry and GogineniBamber and others (2001b). Surface slope and ice thickness are smoothed on a scale of 20 times the local ice thickness.

Balance velocity is a two-dimensional depth-averaged quantity, and for the flowline calculation a three-dimensional velocity field is required. The vertical velocity is obtained by making some assumptions about the vertical distribution of the horizontal velocity and using continuity. A shape function is used to define the vertical profile of the velocity. The shape function and flowline model are described in the Appendix.

This framework of a velocity field and a flowline model allows an IRH surface (an isochronous layer) to be inverted and a mean accumulation rate over the age ofthe layerto be obtained. An initial balance velocity is calculated using the ice-sheet geometry datasets, already mentioned, together with the latest synthesis of accumulation data for Greenland (the OSU dataset). Using a flowline model and a prescribed shape function, particle paths are traced backwards from starting points defined by the isochronous layer towards the surface. The particle being followed starts off with an age equal to the age of the layer. Its trajectory is followed until the particle has an age of zero. If the accumulation history along the path is correct, the particle will be at the surface of the ice, as defined by the DEM. After particles are traced for the whole isochronous layer, a modelled ice surface is obtained. Any mismatch between the modelled ice surface and that of the DEM is assumed to be due to changes in accumulation rates. We use this mismatch to calculate an updated accumulation field. However, because the balance velocity depends on the accumulation rate, the velocity used to update this rate is no longer the velocity that would have resulted given the updated accumulation field.

To overcome this circularity, the updated accumulation rate is used as input to the balance-velocity calculation to obtain a new balance velocity. This velocity can then be used to recalculate the accumulation rates. This continues until the accumulation and velocity converge. Typically five to six iterations are required to obtain differences of <1% for the whole region. A 1% error corresponds to a maximum absolute accumulation-rate difference of ∼0.0025m a–1 for the areas of higher accumulation rate in the region. Across most of the region, the accumulation rate and balance-velocity fields converge after two iterations. The last three to four iterations reduce the errors towards the edges of the region associated, primarily, with the paucity of data points available for interpolation (Fig.3) in those areas.

Fig. 3. Derived accumulation anomaly (in m a_–1) compared to the OSU dataset. The dot-dashed contour line represents no anomaly. The dashed line marks the –0.02 m a_–1 contour representing -15% difference for the region of the contour line. The solid contour lines indicate where OSU accumulation is lower than the derived rate. The 0.02 ma_{– 1}represents a difference of -15% in the northwest and -10% in the south of the area. The grey points mark the data points used for the interpolation. There is a reduction in accumulation rate from west to east between 76° and 78° N.

The balance-velocity code requires accumulation rates at each point of a two-dimensional grid. (The grid resolution for all the balance-velocity datasets is 2.5 km 62.5 km.) The accumulation rates calculated from the particle paths, however, result in a non-uniform distribution mirroring the flight-line distribution (Fig. 3). Since accumulation rate is a relatively smooth field, the non-uniform data are gridded to the balance-velocity grid using bilinear interpolation and a Gaussian smoothing filter. The magnitude of the resulting accumulation-rate anomaly calculated for the deeper layer is shown in Figure 3. This shows the difference between the updated accumulation rate and the reference (OSU) accumulation rate.

Figure 3 shows that around the GRIP site there is reasonably good agreement between the OSU dataset and the mean accumulation rate for the last 6956 years. The 0.02 ma–1 contour represents ∼10% difference in the southern part of the region where accumulation rates are higher. The GRIP age/depth relationship was used to date the layers, but during the iterative calculation of accumulation rates the layers are no longer tied to the GRIP data. In theory, the calculated accumulation rate at the GRIP site could differ considerably from the ice-core derived rate. The agreement in this region therefore gives us confidence in our methodology. To the north, there is a distinct east–west gradient, with a 0.02 ma–1 lower palaeo-accumulation rate to the east of the ice divide and a region of 0.04 ma–1 higher rates to the west. These patches coincide with regions that have little data coverage in the OSU dataset. The difference between the accumulation-rate fields could therefore be a consequence of data interpolation. If these differences are purely due to climate changes, then the results suggest that for North Greenland the east–west gradient in accumulation rates has been reducing during the Holocene.

There are some regions where our calculated accumulation rates are the result of the interpolation. To overcome this, further flight-lines need to be processed to fill in the gaps. Despite this, the accumulation pattern to the north, associated with the denser flight-line coverage, is robust. The same spatial distribution of accumulation is obtained from both layers and for a range of shape-function parameters. For example, the magnitude of the accumulation rate is reduced by 7% when the shape-function parameter increases, from 10 to 20, as a result of changing velocity magnitudes in the constant-strain region (see Appendix), but the spatial pattern remains the same.

The sensitivity of the results to the age estimate was also investigated by calculating the accumulation rate assuming the layer to be 50 years older. The accumulation rate should be fractionally smaller if the layer is older, and over the flowlines themselves (grey lines in Fig. 3) this is the case, with the rate being ∼0.002 ma–1 lower. However, the rates differ by a larger amount (up to ±0.005m a–1, ∼5%) over the two regions where there are no flight-lines. It is unclear how much of this larger difference is due to the differences in layer age and how much is due to the interpolation scheme.

Modelled isochronous surfaces

Ideally, the observed IRHs wouldbe compared with isochronous surfaces calculated within an evolving three-dimensional ice-sheet model. To date, efforts to incorporate such a tracer transport scheme have been largely unsuccessful. A notable exception is the recent paper by Reference Clarke and MarshallClarke and Marshall (2002). Their novel approach and use of a Lagrangian transport scheme appears to be successful. Developing such a tracer transport model is, however, beyond this study.

Instead the velocity fields from two numerical simulations of the Greenland ice sheet, run to steady state, are used to calculate flowlines and, hence, isochronous layers. This mirrors the use of balance velocities earlier. The numerical model used is a newly developed three-dimensional thermomechanical ice-sheet model based on the model of Reference PaynePayne (1999). The climate forcing employed by the model uses a constant accumulation-rate field together with a lapse rate and a positive degree-day model to obtain modelled accumulation/ablation rates. Thus, in terms of model forcing, it is the spatial pattern of accumulation rate, rather than its magnitude, that is fundamental. Based on ideas presented by Reference Ritz, Fabre and LetréguillyRitz and others (1997), the degree-day model was tuned so that the modelled steady-state ice-sheet extent, simulated using the OSU accumulation rate, was in reasonable agreement with the current Greenland ice sheet.

Using the modelled steady-state velocity field, an isochronous surface with an age of 6956 years is calculated and compared with the observations. The isochronous surface is simply found by following particle paths for 6956 years. The particle paths are selected so that they terminate under the flight-lines used to create the observed layers. This allows a direct comparison to be made.

Figure 4 presents the observed 6956year IRH together with two modelled 6956 year layers. The observed layer (OL; solid line) is for the flight-line section highlighted in Figure 1, a section of which is shown in Figure 2. This section passes through the two regions of largest fractional accumulation-rate differences found between the OSU dataset and the IRH-derived accumulation rate (Fig. 3). The first model layer (ML1; dotted line) results from a model run using the OSUdataset. Theother layer (ML2; dashed line) is a result of using the IRH-derived palaeo-accumulation for the numerical simulation.

Fig. 4. Modelled isochronous-layer depths compared with observed IRH. The solid line (OL) is the result of tracing a layer through the data collected along the flight-line highlighted in Figure 1. By linking with flight-lines passing over the GRIP ice-core site, this layer is dated to 6956 BP. The dotted line (ML1) is the 6956 isochrone simulated from an ice-sheet model run using the OSU accumulation rate. The dashed line (ML2) demonstrates the marked improvement in the modelled isochrone when the updated accumulation rate is used to drive the model.

Using the updated accumulation rate results in a significant improvement in the fit between modelled and observed data. Comparing ML1 with ML2 clearly shows the east–west nature of the anomalous accumulation corrected for. The oscillations in both modelled profiles at the eastern end of the section are associated with the start on the northeast Greenland ice stream within the model. This improved fit is not limited to the one flight-line segment presented in Figure 4, as demonstrated by the following statistics. The mean difference between the observed layer and the modelled layer using the OSU accumulation rate is –70.9m with a standard deviation of 91.5 m. For the modelled layer using the updated accumulation, the mean difference with the observed layer drops to –33.1m and the standard deviation is 63.8 m.

To test the sensitivity of the results to the age estimate of the layer, the IRH depths were also calculated assuming ages of 6896 and 7016 BP. At this depth a ±60 year error in the age estimate results in a difference of ±6±1m to the modelled IRH depth.

Discussion and Conclusion

This paper has demonstrated the potential of RES data to generate palaeo-accumulation maps. The methodology proposed uses balance velocities together with a shape function enabling a flowline model to be used. This has advantages over the usual“sandwich”-type models. Not only are accumulation rate and geometry allowed to vary along the flowline for the age calculation at a point, but with balance velocities flow-direction information is also available. Reference Fahnestock, Abdalati, Joughin, Brozena and GogineniFahnestock and others (2001a) used an approach based on a combination of Reference NyeNye (1957) and Reference Dansgaard and JohnsenDansgaard and Johnsen (1969) models to calculate accumulation and basal melt rates required to match the observed IRHs. They reach some startling conclusions with regard to the basal melt rates required.

To interpret the changing separation of IRHs purely based on their one-dimensional model, requires that the radar transects used follow flowlines. For transects that are not parallel to flowlines (as is the case for some of the transects that they present) the upstream flow history, which is likely to be important, is unknown. Neglecting the three-dimensional nature of the problem, especially with regard to flow in the vicinity of large basal topography, may affect the interpretation of their results. Though more computationally intensive, our methodology goes some way to overcoming this restriction.

Further flight-line processing is obviously required to improve the current spatial coverage and limit the influence of the interpolation scheme. This should now be possible for most of the Greenland ice sheet above the ∼2000m ice-thickness contour. Closer to the margin, IRHs become harder to follow. IRH-derived accumulation rates are a useful additional source of data for the compilation of accumulation-rate datasets, as they offer large spatial coverage over inaccessible regions to complement more marginal and coastal datasets.

No correction was applied to our depth estimates to account for the variation of density of ice that occurs for the top 100m or so. One common technique to compensate for this is the addition of a few O(10) m to the calculated depth. Reference Fahnestock, Abdalati, Luo and GogineniFahnestock and others (2001b), for example, used a value of 8 m. Without supporting data on the vertical variation of ice density and how it varies across the region, the best that can be done is to use such a constant value for the correction. This is straightforward to apply. Applying any correction globally will change the absolute accumulation rate, but the spatial pattern, accumulation-rate gradients and their interpretation remain unchanged. An 8m difference in depth results in a difference in dating of 68 years at 1200 m. Making a crude estimate for the error in accumulation rate that represents, we obtain 0.0011 ma–1 (8m/(6956+68)). For numerical modelling studies, it is the prescription of the accumulation-rate distribution that is fundamental, and such small differences in value can currently be considered negligible.

For this study a spatially constant shape function was used. The IRHs used lay no deeper than one-third of the ice thickness below the surface, and are thus within the constant-strain region of the shape-function profile. In addition, the region of study is the central part of North Greenland, and only small variations of the vertical velocity profile are expected along a particle path. It is therefore reasonable to use a constant shape function. However, allowing the vertical profile for horizontal velocity to change along the flow-line could offer potential improvements to the model used to invert for accumulation rate, particularly when studying older layers. Following Reference Nereson and RaymondNereson and Raymond (2001), two profiles could be used: a vertical profile for flow proximal to the ice margin, and a distal profile. The shape function at a point along a flowline is then a weighted combination of the two. Reference Nereson and RaymondNereson and Raymond (2001) used model-derived profiles along individual survey lines. Whether it is possible to develop two generic end-member profiles for the whole Greenland ice sheet is questionable, due to the spatial dependence of many parameters that affect the profile (e.g. the latitudinal dependence of ablation rates). It may be possible, however, at the scale of ice-divide separated basins.

The numerical-model results for IRH depth show the improvements that are possible with improved accumulation-rate estimates and indicate how IRHs could provide an additional way to verify numerical-model results in future. Numerical-model runs also provide the information required for selecting and defining shape functions. As previously mentioned, using regionally based shape functions could improve the results. It should be noted that these are still relatively crude comparisons requiring the model to be run to steady state. A better approach would be to include a tracer transport model (Reference Clarke and MarshallClarke and Marshall, 2002) and this is currently under development.