Skip to main content Accessibility help


  • Access
  • Cited by 8
  • Cited by
    This article has been cited by the following publications. This list is generated based on data provided by CrossRef.

    Eisen, Olaf Nixdorf, Uwe Wilhelms, Frank and Miller, Heinrich 2004. Age estimates of isochronous reflection horizons by combining ice core, survey, and synthetic radar data. Journal of Geophysical Research: Solid Earth, Vol. 109, Issue. B4,

    Brandt, Ola Björnsson, Helgi and Gjessing, Yngvar 2005. Mass-balance rates derived by mapping internal tephra layers in Mýrdalsjökull and Vatnajökull ice caps, Iceland. Annals of Glaciology, Vol. 42, Issue. , p. 284.

    Parrenin, F. Hindmarsh, R.C.A. and Rémy, F. 2006. Analytical solutions for the effect of topography, accumulation rate and lateral flow divergence on isochrone layer geometry. Journal of Glaciology, Vol. 52, Issue. 177, p. 191.

    Carter, Sasha P. Blankenship, Donald D. Peters, Matthew E. Young, Duncan A. Holt, John W. and Morse, David L. 2007. Radar-based subglacial lake classification in Antarctica. Geochemistry, Geophysics, Geosystems, Vol. 8, Issue. 3, p. n/a.

    Parrenin, Frédéric and Hindmarsh, Richard 2007. Influence of a non-uniform velocity field on isochrone geometry along a steady flowline of an ice sheet. Journal of Glaciology, Vol. 53, Issue. 183, p. 612.

    Rogozhina, I. Martinec, Z. Hagedoorn, J. M. Thomas, M. and Fleming, K. 2011. On the long-term memory of the Greenland Ice Sheet. Journal of Geophysical Research: Earth Surface, Vol. 116, Issue. F1, p. n/a.

    Panton, Christian 2014. Automated mapping of local layer slope and tracing of internal layers in radio echograms. Annals of Glaciology, Vol. 55, Issue. 67, p. 71.

    MacGregor, Joseph A. Fahnestock, Mark A. Catania, Ginny A. Paden, John D. Prasad Gogineni, S. Young, S. Keith Rybarski, Susan C. Mabrey, Alexandria N. Wagman, Benjamin M. and Morlighem, Mathieu 2015. Radiostratigraphy and age structure of the Greenland Ice Sheet. Journal of Geophysical Research: Earth Surface, Vol. 120, Issue. 2, p. 212.




      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Using internal layers from the Greenland ice sheet, identified from radio-echo sounding data, with numerical models
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Using internal layers from the Greenland ice sheet, identified from radio-echo sounding data, with numerical models
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Using internal layers from the Greenland ice sheet, identified from radio-echo sounding data, with numerical models
        Available formats
Export citation


Spatially extensive internal layers have been traced in airborne radio-echo sounding (RES) data collected over Greenland during the late 1990s. By linking internal layers within individual flight-lines at crossover points, it is possible to identify spatially continuous layers that are interpreted as isochronous surfaces. Several of the survey lines pass over the GRIP core site, and this allows us to use the published GRIP age–depth relationship to accurately date these surfaces. Two layers, with ages of 3891 and 6956 years BP, have been traced over a large part of North Greenland. Accurately dated and spatially continuous isochrones are valuable for both assimilation within, and verification of, numerical models. For example, comparison of isochronous surfaces from a numerical simulation with those layers observed in RES data can be used to inform the choice of parameters (e.g. rheology) and climate history used to force a numerical model. To demonstrate the potential of the RES data, two layers for North Greenland were used to determine palaeo-accumulation rates. The inversion from layer depth to accumulation rate requires a three-dimensional velocity field. This velocity field is constructed by combining a two-dimensional balance-velocity field with an assumed vertical structure for the horizontal velocity. The isochronous-layer derived accumulation rates were compared with the Bales and others (2001) rates. A larger east–west gradient was found across the central ice divide for the derived accumulation rate, suggesting a trend in the Holocene accumulation rates for this region. The layers were also compared with isochronous surfaces derived from simulations of a three-dimensional thermodynamic ice-sheet model. Using the isochronous-layer derived accumulation rates to force the model improved the match between modelled and observed layers.


Accumulation rate is an important forcing parameter for numerical models. Themost recent synthesis of accumulation data for Greenland is that of Bales and others (2001), subsequently referred to as the OSU (Ohio State University) dataset. They amassed numerous data from snow pits, automatic weather stations and ice cores to produce a dataset covering the whole of the Greenland ice sheet. Despite this substantial undertaking, there is still a large amount of interpolation required, and any technique that can provide accumulation rates over data-sparse regions is of value.

Airborne radio-echo sounding (RES) surveys over Greenland have been going on since the early 1990s as part of the Program for Arctic Regional Climate Assessment (Thomas and others, 2001). Numerous layers are visible within the echo returns. Assuming these layers are isochronous, it is possible to calculate accumulation-rate histories. Where flight-lines cross over each other, the layers can be extended, and, with enough flight-lines, relatively extensive coverage can be achieved. In this paper, we use a technique using the balance velocity to invert a distribution of these layers to obtain a spatial palaeo-accumulation pattern. The results demonstrate the use of RES data to obtain accumulation rates and histories, and the use these can be put to in constraining numerical ice-sheet models.

RES Data

The ongoing RES surveys by the University of Kansas have now covered much of the Greenland ice sheet (Fig. 1). The radar system used to acquire these data has been described by Gogineni and others (1998,2001). In addition to the ice surface and bedrock returns, numerous internal reflecting horizons (IRHs) are observable as a result of englacial dielectric inhomogeneities. The nature of these changes in dielectric properties has been studied by several authors (e.g. Robin and others, 1969; Millar, 1981; Bogorodsky and others, 1985; Fujita and others, 1999; Hempel and others, 2000) and is still under discussion. This study, however, does not concern itself with the exact nature of the IRHs, merely noting that, regardless of the reflection mechanism, they can be identified as isochronous layers.

Fig. 1. Flight-line coverage of the Greenland ice sheet. The 1997, 1998 and 1999 flight-lines are shown in grey. The flight-line sections in black are those used to track the deeper internal reflecting horizon (IRH). Layer tracing starts from the GRIP core site, marked by the grey circle, and extends along flight-lines using crossover points to continue the layer. The flight-line picked out in bold is used for a comparison of modelled and observed IRHs (Fig. 4).

By following a particular IRH through the RES data and matching it at flight-line crossovers, it is possible to produce isochronous surfaces (Fahnestock and others, 2001b). By selecting flight-lines that pass over well-dated ice-core sites, the age/depth relationship determined for the core site can be extended for as far as it is possible to follow an I RH. It is then possible to invert the dated isochronous layer to obtain mean palaeo-accumulation rates.

Layer tracing and dating

The methodology for tracing layers is similar to work previously presented by Fahnestock 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.686108 ms–1 (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.


Two uses associated with improving numerical ice-sheet models are presented here. The first is the calculation of palaeo-accumulation rate, an important forcing function for numerical ice-sheet simulations. The second use is in assessing numerical ice-sheet model performance by directly comparing measured IRH depth with isochronous surfaces generated by a numerical model.

Palaeo-accumulation rates

A common approach to calculating accumulation rates from IRHs is to use so-called “sandwich” models such as those of Nye (1957) and Dansgaard 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. Budd and Warner, 1996; Joughin and others, 1997; Bamber and others, 2000; Huybrechts and others, 2000).

Three datasets are required to calculate balance velocity (see, e.g., Budd 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 (Bamber and others, 2001a), and the ice thickness used is that described by Bamber 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– 1represents 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 Clarke 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 Payne (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 Ritz 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. Fahnestock and others (2001a) used an approach based on a combination of Nye (1957) and Dansgaard 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. Fahnestock 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 Nereson 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. Nereson 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 (Clarke and Marshall, 2002) and this is currently under development.


This work was carried out with funding from the Leverhulme Trust (grant F/182/BJ) and the U.K. Natural Environment Research Council (grant GR3/EOCE/4). We gratefully acknowledge S. P. Gogineni for supplying the RES data required to undertake this study.


Bales, R. C., Mcconnell, J. R., Mosley-Thompson, E. and Csatho, B. 2001. Accumulation over the Greenland ice sheet from historical and recent records. J. Geophys. Res., 106(D4), 33, 813–33, 825.
Bamber, J. L., Hardy, R.J., Huybrechts, P. and Joughin, I. 2000. A comparison of balance velocities, measured velocities and thermomechanically modelledvelocities for the Greenland ice sheet. Ann. Glaciol., 30, 211–216.
Bamber, J. L., Ekholm, S. and Krabill, W. B. 2001a. A new, high-resolution digitalelevation model of Greenland fully validated with airborne laser altimeter data. J. Geophys. Res., 106(B4), 6733–6746.
Bamber, J.L., Layberry, R.L. and Gogineni, S. P. 2001b. A new ice thickness and bed data set for the Greenland ice sheet. 1. Measurement, datareduction, and errors. J. Geophys. Res., 106(D24), 33, 773–33, 780.
Bolzan, J. F., Waddington, E. D., Alley, R. B. and Meese, D. A. 1995. Constraints on Holocene ice-thickness changes in central Greenland from the GISP2 ice-core data. Ann. Glaciol., 21, 33–39.
Bogorodsky, V. V., Bentley, C.R. and Gudmandsen, P. E. 1985. Radioglaciolog y. Dordrecht, etc., D. Reidel Publishing Co.
Budd, W.F. and Warner, R. C. 1996. Acomputerscheme for rapidcalculations of balance-flux distributions. Ann. Glaciol., 23, 21–27.
Clarke, G. K. C. and Marshall, S. J. 2002. Isotopicbalance of the Greenland ice sheet: modelled concentrations of water isotopes from 30, 000 BP to present. Quat. Sci. Rev., 21(1–3), 419–430.
Dansgaard, W. and Johnsen, S. J. 1969. A flow model and a time scale for the ice core from Camp Century, Greenland. J. Glaciol., 8(53), 215–223.
Fahnestock, M., Abdalati, W., Joughin, I., Brozena, J. and Gogineni, P. 2001a. High geothermal heat flow, basal melt, and the origin of rapid ice flow in central Greenland. Science, 294(5550), 2338–2342.
Fahnestock, M. A., Abdalati, W., Luo, S. and Gogineni, S. 2001b. Internal layer tracing and age–depth–accumulation relation ships for the northern Greenland ice sheet. J. Geophys. Res., 106(D24), 33, 789–33, 797.
Fujita, S. and 6 others. 1999. Nature of radio-echo layering in the Antarctic ice sheet detected by a two-frequency experiment. J. Geophys. Res., 104(B6), 13, 013–13, 024.
Gogineni, S., Chuah, T., Allen, C., Jezek, K. and Moore, R. K. 1998. An improved coherent radar depth sounder. J. Glaciol., 44(148), 659–669.
Gogineni, S. and 9 others. 2001. Coherent radar ice thickness measurements over the Greenland ice sheet. J. Geophys. Res., 106(D24), 33, 761–33, 772.
Hempel, L., Thyssen, F., Gundestrup, N., Clausen, H. B. and Miller, H. 2000. A comparison of radio-echo sounding data and electrical conductivity of the GRIP ice core. J. Glaciol., 46(154), 369–374.
Huybrechts, P., Steinhage, D., Wilhelms, F. and Bamber, J. 2000. Balance velocities and measured properties of the Antarctic ice sheet from a new compilation of gridded data for modelling. Ann. Glaciol., 30, 52–60.
Joughin, I., Fahnestock, M., Ekholm, S. and Kwok, R. 1997. Balancevelocities of the Greenland ice sheet. Geophys. Res. Lett., 24(23), 3045–3048.
Millar, D. H. M. 1981. Radio-echo layering in polar ice sheets and past volcanic activity. Nature, 292(5822), 441–443.
National Snow and Ice Data Center (NSIDC). 1997. The Greenland Summit ice cores. Boulder, CO, University of Colorado. Cooperative Institute for Research in Environmental Sciences. National Snow and Ice Data Center, Distributed Active Archive Center.
Nereson, N. A. and Raymond, C. F. 2001. The elevation history of ice streams and the spatial accumulation pattern along the Siple Coast of West Antarctica inferred from ground-based radar data from three inter-ice-stream ridges. J. Glaciol., 47(157), 303–313.
Nye, J. F. 1957. The distribution of stress and velocity in glaciers and ice-sheets. Proc. R. Soc. London, Ser. A, 239(1216), 113–133.
Payne, A. J. 1999. Athermomechanical model of ice flow in West Antarctica. Climate Dyn., 15(2), 115–125.
Reeh, N. 1988. A flow-line model for calculating the surface profile and the velocity, strain-rate, and stress fields in an ice sheet. J. Glaciol., 34(116), 46–54.
Ritz, C., Fabre, A. and Letréguilly, A. 1997. Sensitivity of a Greenland ice sheet model to ice flow and ablation parameters: consequences for the evolutionthrough the last glacialcycle. Climate Dyn., 13(1), 11–24.
Robin, G. de Q., Evans, S. and Bailey, J.T. 1969. Interpretation of radio echo sounding in polar ice sheets. Philos. Trans. R. Soc. London, Ser. A, 265(1166), 437–505.
Thomas, R. H. and PARCA Investigators. 2001. Program for Arctic Regional Climate Assessment (PARCA): goals, key findings, and future directions. J. Geophys. Res., 106(D24), 33, 691–33, 705.


It is convenient to work with a non-dimensionalized vertical coordinate σ which is defined as


where s is the surface elevation and h the ice thickness. Thus at the surface σ = 0 and at the bed σ = 1. A“shape function”, ϕ, is defined as the ratio of velocity, u, to the depth-averaged velocity,


Following Reehs (1988) methodology it can be shown that the vertical strain is given by


where w is the vertical velocity and a is the net accumulation. Note this equation is slightly different to that derived by Reeh (1988) due to the use of a different non-dimensionalization of the vertical velocity (Equation (1)). Vertically integrating Equation (3), assuming no basal motion, the vertical velocity is found to be


Thus, once the shape function has been specified, all three components of velocity are known (using balance velocity for u), allowing particle paths to be tracked. The shape function used in this study has the general form


It can be shown for an isothermal ice sheet that p = 4, but for much of Greenland a more realistic value for p is in the range 10–15 (Bolzan and others, 1995). Values around p = 15 are used for this work. The value of p affects the depth of the basal shear layer. Lower values of p have deeper shear layers. At the surface (σ = 0), ϕ = ( p+1)/p so for smaller values of p the surface velocity is higher. A value of p = 10 gives the surface velocity as 1.1 u, whilst p = 20 gives it as 1.05ū. For the values of p used, the constant-strain region of the profile, for which the horizontal component of velocity is still given by surface velocity, extends down at least as far as the deepest (6956 years) layer investigated. This explains the qualitatively similar results across such a large range of p values.