Introduction
Surface mass balance is a key parameter for estimating the contribution of the Antarctic ice sheet to global sea-level changes and deriving palaeoclimatic conditions from ice cores (Reference Giovinetto and ZwallyGiovinetto and Zwally, 2000). A major problem for determining the net surface mass balance, hereafter referred to as accumulation rate, on continental scales is its irregular spatial characteristics. Regionally homogeneous precipitation is disturbed in the atmospheric boundary layer by small-scale variations in wind speed, leading to varying patterns of the surface mass balance on all scales.
Local changes in accumulation rate over time can be accurately reconstructed from ice cores retrieved at specific locations such as ice domes, provided their position was stationary over time. At other drilling locations, the surrounding setting of bedrock topography, ice thickness and surface slope determines whether observed local changes in annual layer thickness contain upstream effects (Reference HamiltonHamilton, 2004), ice dynamics (Reference Vaughan, Corr, Doake and WaddingtonVaughan and others, 1999) or correspond to real temporal changes in the past (Reference Spikes, Hamilton, Arcone, Kaspari and MayewskiSpikes and others, 2004; Reference Steinhage, Eisen and ClausenSteinhage and others, in press).
In order to estimate the spatial representativity of deep ice-core measurements, drilling activities are usually accompanied by repeated readings of surface stakes, snow-pit sampling and deployment of shallow firn cores. Although ice-dynamic numerical modelling improves the separation of small-scale temporal and spatial effects, the limited spatial resolution of input datasets as well as models makes reliable interpretation of results difficult (Reference KarlöfKarlöf, 2004).
*Present address: Versuchsanstalt für Wasserbau, Hydrologie und Glaziologie, Eidgenössische Technische Hochschule, ETH-Zentrum, CH-8092 Zürich, Switzerland.
Especially over the last decade, several methodological studies on the application of high-resolution ground-penetrating radar (GPR) have been carried out to improve the spatial resolution of accumulation maps (Reference Richardson-Näslund, Aarholt, Hamran, Holmlund and IsakssonRichardson-Näslund and others, 1997; Reference Nereson, Raymond, Jacobel and WaddingtonNereson and others, 2000; Reference Richardson-NäslundRichardson-Näslund, 2001; Reference Frezzotti, Gandolfi and UrbiniFrezzotti and others 2002; Reference PälliPälli and others, 2002; Reference KarlöfKarlöf, 2004; Rotschky and others, 2004). In most cases, internal horizons imaged with GPR represent surfaces of isochronal deposition (Reference GudmandsenGudmandsen, 1975; Reference Jacobel and HodgeJacobel and Hodge, 1995). With its ability to penetrate the upper hundreds of metres of the ice column, equivalent to a deposition history of several hundreds to thousands of years on the Antarctic inland plateau, GPR provides a means to accurately map accumulation rates and their spatial variations along continuous profiles.
In this study, we present the spatial and temporal variation of the accumulation rate in the vicinity of the EPICA (European Project for Ice Coring in Antarctica) deep-drilling site in Dronning Maud Land, Antarctica, (EDML) based on extensive GPR measurements (Fig. 1). The EDML drill site near Kohnen station (75.00˚ S, 0.07˚ E; 2893ma.s.l.) is located close to a topographic ice divide running approximately east-southeast–west-northwest. Repeated global positioning system (GPS) geodetic measurements during the period 2001–04 showed that surface velocities are ~0.7ma–1, with the main flow roughly parallel to the divide (personal communication from D. Schulte, 2004). The combination of continuous internal layers up to 110m depth, obtained during several EPICA pre-site surveys, with physical properties from the EDML ice core allowed us to calculate the spatial distribution of the accumulation rate for eight time intervals during the last 1000 years.
Fig. 1. Location of study area in Dronning Maud Land (top), and radar profiles surrounding the EPICA drill site (bottom). Contours (black) and greyscale refer to surface topography as derived from European Remote-sensing Satellite interferometry (see text for details). Superimposed as white lines are the 2000/01 grid- (thick) and star-like patterns (thin), with the 1998/99 traverse profiles as thick dashed white lines. The EDML drill site is in the centre of the star- and grid-profile patterns. Mean ice flow is from east-southeast to west-northwest at about 0.7ma–1 near the drill site; mean wind direction is from ~60˚ true north (Reference ReijmerReijmer, 2001). The black-and-white dashed line indicates the upwind profile 011203 discussed in the text and shown as a cross-section in Figure 4. The approximate position of the major topographic ice divide is shown as dotted-dashed line.
Our results show that the spatial variations in accumulation rate, related to undulations in surface topography, and temporal variations seen in the ice core over the period considered are similar in magnitude. Interpretation of the EDML ice core therefore requires an understanding of how the spatial and temporal variations interact to form the record preserved in the core.
Data and Methods
Ground-penetrating radar
We used a commercial RAMAC radar set (Malå Geosciences), operating with shielded 250MHz antennae, and simultaneous positioning measurements with a Trimble GPS receiver, to map the internal structure of the upper 120m of the ice sheet in 1m increments. The radar profiles recorded in 2000/01 formed two patterns: the first was a 10 × 10km2 grid with 2 and 3 km spacing roughly parallel to the ice divide and 1 km spacing across. A second, star-like pattern consisted of eight legs up to 25km long. Both patterns were centred on the EDML drilling location (Fig. 1). In addition, profiles recorded during the 1998/99 field season with the same system but unshielded 200MHz antennae were also used for our analysis.
A series of nine continuous internal reflection horizons (IRHs), which are isochrones (Reference Eisen, Nixdorf, Wilhelms and MillerEisen and others, 2004), were tracked in all the profiles throughout the region of interest. Details of data processing, the tracking procedure and examples of the data quality are available in Reference Eisen, Nixdorf, Wilhelms and MillerEisen and others (2004). After tracking, all reflectors were available in the two-way travel-time domain.
Ice-core properties
Dielectric profiling (DEP) along the EDML core formed the basis for converting the structure of the IRHs to accumulation rates. The upper 13m of the EDML core were missing because of the technical set-up of the deep drilling. They were replaced by DEP data from the nearby 150m ice core B32, located 1.5km to the west. The B32 and EDML records were matched by identifying unique volcanic peaks visible in both DEP conductivity records (personal communication from A. Lambrecht, 2003). This enabled us to transfer the age–depth profile from B32 to EDML, where the independent counting of annual layers is still underway. Dating uncertainty is ~5 years (Reference Traufetter, Oerter, Fischer, Weller and MillerTraufetter and others, 2004). Based on the age–depth model and density measurements along the EDML ice core, we determined the accumulation rates at the EDML drilling location used below.
Depth distribution of the electromagnetic wave speed and density was calculated from the DEP permittivity, allowing conversion of travel time to depth along the profiles. The integration of density yielded cumulative mass as a function of depth, thereby smoothing out small-scale variations in density. As we did not find significant spatial deviations of the mass–depth distribution from an analysis of eight ice and firn cores over >320km (Rotschky and others, 2004), we assumed that it was constant throughout the region. The age of each reflector was determined by merging reflector depth and the age–depth distribution at EDML (for details of the IRH dating procedure see Reference Eisen, Nixdorf, Wilhelms and MillerEisen and others, 2004).
Synthesizing accumulation rates
The above datasets provided cumulative mass and age for each point along each IRH. The accumulation rate along the profiles was derived from adjacent pairs of IRHs for different time intervals by calculating the ratio of difference in cumulative mass and difference in age. The age difference of neighbouring IRHs ranges from <100 to 200 years. Hence, the calculated accumulation rates represent a low-pass filtered mean of the annual accumulation rate. The cumulative effect of strain thinning is comparably small in the upper part of the ice column. Vertical strain rates were on the order of 1–5 ×10-5 a-1 (personal communication from H. Oerter and O. Rybak, 2004), two orders of magnitude smaller than the value found to be significant for internal deformation of an ice cap by Reference Vaughan, Corr, Doake and WaddingtonVaughan and others (1999).
The spatial distribution of accumulation rate was obtained by gridding the dataset with a spatial increment of 0.5km (Fig. 2). Prior to gridding with a tensioned minimum curvature technique, the data were spatially averaged within each gridcell to avoid aliasing of small wavelengths (Reference Smith and WesselSmith and Wessel, 1990). Within the 10km grid, this resulted in one interpolated data point between adjacent radar profiles. However, in the outermost parts of the star-like pattern, adjacent radar profiles were up to 20 km apart. The maximum distance of interpolated points to the nearest profile was therefore on the order of the typical correlation lengths of accumulation rates, which were ~10km in this area (Rotschky and others, 2004). Careful consideration of these values is therefore required.
Fig. 2. Regional distribution of accumulation rate (grey values and contours) averaged over the period AD 1848–2001, derived from the uppermost isochronous IRH along the radar profiles. Line types as in Figure 1.
The error in accumulation rates depends on the dating of the IRHs, conversion from travel time to depth, and integration of mass. From crossover point analyses within the grid we found a typical uncertainty of 5 ns in travel time, or 0.5m in depth. This equals a statistical error of ~1.5 kgm–2 a–1 for accumulation rates close to the radar profiles. Dating uncertainty results in a systematic error in accumulation rates, ~1kgm–2 a–1. For discussion of the present accumulation rate pattern below, we use the accumulation rates derived from the depth of the uppermost continuous horizon, dated to AD 1848. The resulting pattern of accumulation rates (Fig. 2) therefore presents the mean over the period AD 1848–2001.
Surface topography
Surface topography was available from kinematic GPS measurements accompanying each GPR profile. The static GPS base station at Kohnen, tied into the global GPS network, was used as reference for the GPS rover with a time interval of 1 s. The kinematic GPS data result in a relative accuracy of surface elevations of 0.1 m. Unfortunately, only within the 10km grid surrounding EDML were the GPS data spatially dense enough to interpolate a reasonably accurate surface topography. We therefore used satellite remote-sensing data to construct a digital elevation model (DEM) of the whole region, which provides an accurate datum for each gridpoint.
The topography was derived by European Remote-sensing Satellite (ERS-1 and -2) synthetic aperture radar (SAR) interferometry using two pairs of satellite data acquired in 1996 (Fig. 1). The ice motion and temporal changes in the investigation area were neglected. Precise orbits were refined using the GPS-based topography along our GPR profiles and ICESat’s Geoscience Laser Altimetry System (GLAS) data (http://nsidc.org/data/icesat/). A coherence adaptive filter was applied to reduce the noise of the interferogram. Noise in the DEM was further reduced by averaging the original 25m pixel spacing of the geocoded product with an estimated locational accuracy of 100m onto the same grid used for gridding of the GPR data with 0:5 × 0:5km2 resolution (Fig. 1). The rms difference of the DEM in respect to the reference topographies is 0.9 m.
Regional Distribution Of Accumulation Rate
Point sampling of accumulation rates by means of stake measurements, snow pits and ice cores during several field seasons in the 1990s yielded a generally decreasing trend of accumulation rate in DML away from the coast (Reference Isaksson, Karlén, Gundestrup, Mayewski, Whitlow and TwicklerIsaksson and others, 1996; Reference OerterOerter and others, 2000). The data synthesis of our ice-core and regional radar measurements confirms this trend, but with a much enhanced spatial resolution (Fig. 2).
General characteristics
The mean accumulation rate at EDML from AD 1848 to 2001 is ~64 kgm–2 a–1. The regional accumulation rate gradient is roughly perpendicular to the ice divide (Fig. 2). Based on the 153 year average, the accumulation rate decreases over a distance of 40km in the north-northeast–south-southwest direction from ~72 kgm–2 a–1 to 54 kgm–2 a–1, almost 30% of the accumulation rate at EDML. The average decreasing gradient perpendicular to the ice divide is ~0.5 kgm–2 a–1 km–1, whereas the mean spatial trend in accumulation rate is close to zero parallel to the ice divide.
The general pattern is overlain by small-scale variations in accumulation rate, on the order of 5–15% of the accumulation rate at EDML over a few kilometres, the origin of which is discussed further below. Within the 10km grid, where values reliably quantify the variations, maximum local accumulation rate gradients are up to 2.5 kgm–2 a–1km–1, about five times the regional gradient (Fig. 3). Even higher gradients occur on longer profiles outside the grid. The drill site is located on the eastern side of a southward-stretching tongue of higher accumulation rate, about 4 km wide (Fig. 3). Further upstream towards the east, accumulation rates decrease to <60 kgm–2 a–1 about 6 km away, before increasing again to 65 kgm–2 a–1 (Fig. 3).
Fig. 3. Close-up of Figure 2, focusing on the 10 km grid, with radar profiles as shown above. Note that the range of grey values has changed.
Relation to surface topography
Snow accumulation rate in this region is mainly determined by precipitation, wind erosion and deposition (Reference Van den BroekeVan den Broeke and others, 1999). Varying wind speeds and turbulences in the boundary layer lead to disturbances in the regionally homogeneous precipitation pattern, finally resulting in inhomogeneous spatial accumulation rate (Reference King, Anderson, Vaughan, Mann, Mobbs and VosperKing and others, 2004). To qualitatively identify the relation of surface topography and the observed accumulation rate pattern, we compare both values along profile 011203 (Fig. 4), running in the upwind direction of EDML, which is from ~60˚ true north (Reference ReijmerReijmer, 2001). The elevation along the profile rises gently east-northeastwards by an average of 10m over a distance of 15 km. The general increase is interrupted by two local elevation maxima at distances of 7.5 and 13.5 km, separated by a local minimum around 8.5 km. The first convex rise is 2 km wide and about 0.5m above the surrounding surface. The second rise is 4 km wide with a height of 2m (Fig. 4).
Fig. 4. Accumulation rate for the time interval AD 1848–2001, World Geodetic System 1984 (WGS84) surface elevation, and slope (dotted line) along upwind profile 011203 (Figs 1 and 2). Elevation is from kinematic GPS measurements; the magnitude of slope is derived from the DEM (Fig. 1) over segments of 1 km. Shaded areas mark the region of lower accumulation rates, as discussed in the text. The main wind direction is parallel to the profile from east-northeast to west-southwest (right to left).
The accumulation rate along this profile shows two pronounced minima along the same line (Fig. 4). For the first 3 km from EDML it is almost steady, then it drops by 3kgm–2 a–1 over 2 km and increases again to the same level as at EDML at 8.5 km. The second minimum has slightly different characteristics than the first one. From its lowest value at ~11km from EDML, the accumulation rate increases strongly over 5 km with a gradient of up to 2.5 kgm–2 a–1 km–1 and starts to level out around 16km (Fig. 4). A small local maximum occurs about 19km from EDML.
Comparison of the elevation and accumulation rates along this profile reveals two striking features: first, the local accumulation rate maximum at 8.5km coincides with the local elevation minimum; second, both accumulation rate minima are on the leeward side of two elevation maxima and occur at locations where the local slope has a maximum. Moreover, the drop in accumulation rate downwind along the profile at 11km is more than twice as deep (10 kgm–2 a–1) as the accumulation rate minimum at 6 km. The larger drop seems to be related to the height of the upwind surface undulation and local slope.
Similar relations between surface undulations and accumulation rate patterns have been observed in other parts of Antarctica (e.g. the East Antarctic megadune fields (Reference Frezzotti, Gandolfi and UrbiniFrezzotti and others, 2002) and Lyddan Ice Rise (Reference King, Anderson, Vaughan, Mann, Mobbs and VosperKing and others, 2004). Our interpretation of the observed features is in line with current understanding of the interaction of surface topography, wind field and accumulation rate. Undulations in the surface topography lead to variations in wind speed and turbulences, in turn resulting in varying depositional and erosional processes, as recently modelled by Reference King, Anderson, Vaughan, Mann, Mobbs and VosperKing and others (2004). Accumulation rate is lowest where wind speeds are likely to be higher than in the surrounding area, i.e. where surface gradients are larger. This leads to the ‘inphase’ relation of local slope maxima and minima of accumulation rate (Fig. 4). Concave depressions, on the other hand, are likely places of lower wind speeds or even turbulences, resulting in higher accumulation rate. This is the case at the accumulation rate maximum at 8.5 km. Likewise, the small maximum around 20km is also located in a small surface depression, which extends to about 23km from EDML, downwind of a slope with higher gradients.
The observed accumulation rate anomalies, unlike those of East Antarctic megadunes, seem to be stationary (Fig. 5). The internal structure is advected with the mean ice flow of <1ma–1. It does not provide a consistent picture of surface features moving in the main wind direction, such as the displacement of the apexes of internal layers over time found by Reference Frezzotti, Gandolfi and UrbiniFrezzotti and others (2002).
Fig. 5. Accumulation rates along the upwind profile 011203 (Fig. 2) for several time intervals (years AD). Based on the age–depth profile transferred from B32 to EDML and the density measurements at EDML, we calculated a time series of accumulation rate at EDML from AD 985 to 2001 (see text). The mean of accumulation rate over the complete period at EDML and the standard deviation of the 100 year filtered time series are indicated in error-bar type (grey) on the left.
Temporal Variations
Patterns of accumulation rate derived for other time intervals over the >1000 year period show qualitatively the same distribution as that discussed above for the most recent time interval (Fig. 2). We show an example of this for the upwind profile in Figure 5. Within each time interval, the position and value of spatial deviations from the mean accumulation rate are fairly constant. This means that the actual values of a distribution vary with the mean accumulation rate of the respective time interval. An exception is the region 19– 22km from EDML, where the four intervals before AD 1400 show lower accumulation.
The mean accumulation rate at EDML throughout the period amounts to 63 kgm–2 a–1, as calculated from matching with the age–depth model of ice core B32 described above. To compare the spatial variations derived from the radar data with the temporal variations of the ice core, both records need to be based on comparable timescales. We therefore filtered the ice-core time series of accumulation rate with a 100 year running mean. This equals the order of magnitude of the time difference of adjacent IRHs. The corresponding standard deviation that resulted was 2.1 kgm–2 a–1. As seen in Figure 5, the average accumulation rate within each time interval varies from the mean of the whole period by approximately the same amount. In contrast, the spatial standard deviation of accumulation rate for a single time interval along this profile is ~3.5 kgm–2 a–1, almost twice as large.
As the EDML ice core is not located at an ice dome, its signal is formed along the upstream particle trajectories. From comparison of the temporal changes in accumulation rate observed in the ice core and the regional variations indirectly derived from radar surveys, it is evident that temporal and spatial changes are of the same order of magnitude. If the accumulation rate pattern varies along these trajectories, the ice-core signal also contains a spatial component stemming from advection. This may be part of the reason for the slightly lower accumulation between 19 and 22km prior to AD 1400. An accurate interpretation of the ice-core record requires separation of temporal and spatial changes along the particle paths.
Implications and Conclusions
The accumulation rate pattern derived from joint GPR and ice-core data presented in this study reveals local variations in accumulation rate which have not previously been measured in a comparable spatial two-dimensional resolution in the vicinity of the EPICA deep-drilling site. Despite the low gradients in surface topography, tiny undulations alter accumulation rates by >10% over only a few kilometres.
Regarding climatic interpretations of ice-core records in terms of accumulation rates, these small-scale variations could lead to misinterpretations because of advection of the ice-core location through areas of different accumulation rates. Deconvolution of the ice-core record with respect to the spatial accumulation signal is therefore mandatory, and requires detailed knowledge of the upstream accumulation and velocity characteristics. Full understanding of the ice-dynamic mechanisms involved is also important when records from depth ranges subject to strain thinning are interpreted.
State-of-the-art satellite remote-sensing methods, such as passive microwave radiometers (Reference Winebrenner, Arthern and ShumanWinebrenner and others, 2001), cannot yet map the accumulation rate to the degree of resolution and accuracy that is possible with ground-based measurements, but progress is expected during the next few years, with more accumulation datasets for ground-truthing becoming available from the satellite-borne active radar sensor (Reference Drinkwater, Long and BinghamDrinkwater and others, 2001). The data presented here and in comparable studies (e.g. in this volume) will moreover make a valuable contribution to calibrating and validating satellite measurements, on board ICESat and CryoSat, aimed at determining the mass balance of the Antarctic continent. Once highly accurate topographies from those missions are available, their combination with numerical boundary-layer models and meteorological data will provide estimates of accumulation rates even in those areas of Antarctica where direct measurements are unavailable.
Acknowledgements
We acknowledge the contribution of the field parties during data acquisition, and especially G. Stoof for maintaining the system. Altimetry data were provided by the US National Snow and Ice Data Center, Boulder, CO (http://nsidc.org/data/icesat/). The ERS data were made available by the European Space Agency through the VECTRA project. Comments from C. Richardson-Näslund, D.G. Vaughan and an anonymous reviewer improved the manuscript. Preparation of this work was supported by the Deutsche Forschungsgemeinschaft grant Ni493/1 and two scholarships of the Studienstiftung des Deutschen Volkes. This work is a contribution to the ‘European Project for Ice Coring in Antarctica’ (EPICA), a joint European Science Foundation (ESF)/European Commission (EC) scientific programme, funded by the European Commission and by national contributions from Belgium, Denmark, France, Germany, Italy, the Netherlands, Norway, Sweden, Switzerland and the United Kingdom. This is EPICA publication No. 148.
