2.1. Climatological observations

Meteorological information is available through the Greenland Climate Network (Reference Steffen, Box and AbdalatiSteffen and others, 1996) from an AWS at Crawford Point, in the vicinity of T05. The station has been in operation since April 1995, and the available record includes hourly mean values of surface air temperature, wind speed and direction as well as snow temperatures (Reference Steffen and BoxSteffen and Box, 2001). We extract mean surface air temperatures for the three warmest months: June, July and August (JJA). In addition, we calculate a modified positive degree-day factor for each year. This factor is the sum of daily mean temperatures above 0°C during a particular period, typically a year (Reference BraithwaiteBraithwaite, 1995). Daily mean temperatures generally do not exceed the melting point at this location, as warm air masses are effectively cooled over the ice surface. The factor used here is based on positive hourly mean temperatures and referred to as the positive degree-hour (PDH) factor to distinguish it from the classical positive degree-day factor. In order to yield comparable dimensions, hourly mean values are weighted by a factor of 1/24. Artefacts caused by insufficient ventilation of the temperature sensor are avoided by taking wind speeds above 1 m s⁻¹ as a constraint for the inclusion of temperatures. The PDH factor is only minimally altered by this constraint.

To determine the onset of melt we use the day of year (DOY) at which two threshold values of PDH factors are reached: DOY_0.1 (PDH factor = 0.1) indicates the DOY when the first short events with above-zero temperatures occur; DOY₁.₀ (PDH factor = 1.0) provides the DOY at which a stage of more persistent positive temperatures is present, likely causing some surface melt. For the time period between these days, the mean-snow temperature has been calculated using the record of a thermistor deployed at an initial depth of 1.1 m below the surface. This set of parameters describes the magnitude of meltwater generated at the surface as well as the cold content of the snow at the onset of melt. It gives a qualitative picture of the conditions to determine meltwater percolation and refreezing during a particular summer season. A detailed analysis of the AWS data is beyond the scope of the present study.

Data analysis indicates warm summers in 1995, 1999, 2002, 2003 and 2004 (Table 1). The years 1995, 1998, 2002 and 2003 were reported earlier as years of extreme melt throughout large parts of the Arctic, both in terms of the aerial extent and total amount of melt (ACIA, 2004; Reference Steffen, Nghiem, Huff and NeumannSteffen and others, 2004). Enhanced summer melting during these years covers the western slope of the Greenland ice sheet. The PDH factors from 1995 to 2004 range from 0.3 to 11.0. The highest values (>9.0) occurred in 1999, 2002, 2003 and 2004, a moderate value in 1995, and small values (<5.0) in the other years. The DOY_0.1 is in the range 142–188, and DOY₁.₀ varies from 152 to 224. The time period between these two threshold values ranges from 2 to 44 days. In 1996, the year with the coldest summer and lowest PDH factor, the second threshold value was not reached at all. The period between the threshold values is generally longer for cold summers than for warm summers. Snowpack temperatures at the onset of melt range from −15.9 to −19.0°C, indicating a large cold content of the snow for all years in the record.

2.2. Snow-pit and firn-core data

During spring (19 April–13 May) and late summer 2004 (28 August–21 September), stratigraphic records were retrieved from the survey area from nine snow pits and shallow firn cores. The sites were located at point T05 and at distances of 1, 10, 100 and 1000 m east-northeast and south-southeast from T05, parallel and perpendicular to the EGIG line. These points are referred to as T05-E1 to E4, and T05-S1 to S4, respectively. Observations and measurements comprise visual snow stratigraphy, a qualitative description of layers and snow density. In spring, the last-summer melt surface (LSS) was reached in snow pits at a mean depth of 1.43 m. The LSS was characterized by a flat and hard ice layer beneath the late-summer depth hoar (Reference ParryParry and others, 2007). The mean snowpack density ρ of all sites was 420 kg m⁻³. Density variations in the winter layer were small. From the bottom of each pit, a shallow core with an approximate length of 2 m was drilled. The cores penetrated into a heterogeneous zone of metamorphosed snow, containing ice lenses and discontinuous ice layers below the LSS. At point T05-E3, 100 m east of T05, a longer core down to 18.9 m depth below the surface was additionally retrieved.

Repeated measurements in late summer reached the LSS at 1.51 m depth (5.3% increase compared to the spring depth) and a mean density of 530 kg m⁻³ (+26.2% compared to spring). Consequently, the water equivalent (w.e.) depth of the layer accumulated from the previous summer increased from 60.5 cm in spring to 79.6 cm in late summer, an increase of 31.6% (Reference ParryParry and others, 2007). The density profiles showed more complex variations in late summer than in spring due to the presence of solid ice structures. Ice layers showed little spatial continuity, even at length scales as short as 1 m. An exception was a spatially continuous, 10–15 cm thick ice layer located at 3–3.5 m depth, potentially related to extreme summer melting in 2002 (Reference Scott, Mair, Nienow, Parry and MorrisScott and others, 2006). Ice layers were present at all depths of the snow pits, more or less evenly distributed with no preferable depth of higher number density.

2.3. GPR data acquisition and processing

The GPR data were collected from 29 April to 11 May 2004 using a commercial radar system (RAMAC, Malå GeoScience) operating with shielded antennae at frequencies of 500 and 800 MHz. The GPR profiles form a grid of 100 m × 100 m with 10 m line spacing and were oriented parallel (west-southwest–east-northeast) and perpendicular (southsoutheast–north-northwest) to the EGIG line. Additional GPR measurements were performed along 1 km sections forming a square (Fig. 1). A global positioning system (GPS) receiver was operated together with the GPR for simultaneous kinematic positioning. Static GPS measurements were additionally performed at the beginning and end of each profile. Gaps in the GPS data, resulting from the kinematic operation, were filled by linear interpolation between nearest kinematic or static points. A GPS reference station was operated at T05 and used for post-processing differential GPS data.

GPR shots were triggered every 0.1 m by means of an odometer attached to the sledge, which was pulled by hand. However, some slip caused the actual trace interval to vary between 0.10 and 0.12 m, reducing the number of traces along the profiles by up to 15%. To increase the signal-to-noise ratio during data acquisition, eight traces were recorded at each shot point location at a repetition rate of 200 kHz (quasi-instantaneously) and stacked. The receiving time window was 130 ns, corresponding to a depth range of about 13.5 m with 1024 samples and a sample interval of ∼0.13 ns. The wavelength λ in dry firn (ρ = 500 kg m⁻³) is ∼0.42 m and ∼0.26 m for the 500 MHz and 800 MHz antennae, respectively. The GPR system is designed such that the centre frequency corresponds to the bandwidth. The theoretical resolution (λ/4) of the 500 MHz and 800 MHz antennae is therefore about 0.1 m and 0.06 m, respectively. This is limited by the length of the transmitted wavelets, however, comprising two to three cycles at 500 MHz and three to five cycles at 800 MHz, as interferences from partial reflections at inhomogeneities within the length of the wavelets occur.

Radar data were processed using the seismic software Focus/Disco (Paradigm Geophysical). Processing steps included static correction, bandpass filtering, complex trace analysis and gain adjustment (linear gain function and automatic gain control). For the static correction, individual traces were shifted such that the arrival times of the direct waves correspond to the theoretical arrival time in air, determined by the antennae spacing. Semi-automatic tracking routines in the Landmark software package SeisWorks2D (Halliburton) facilitated the tracking of continuous IRHs (Fig. 2a). To ensure that the same horizon was tracked throughout the record, intersections of previously tracked horizons were displayed as markers in the current profile. Cross-point errors in IRH depths from intersecting profiles are thus reduced.

Fig. 2. Sample 500 MHz GPR profile, aligned approximately east–west as indicated in the inset in (a). Panels on the left (a, b, c) cover 100 m; panels on the right (a′, b′, c′) provide a 10 m wide close-up. (a, a′) Statically corrected, filtered and gain-corrected radargram, displaying signal amplitude; (b, b′) as (a), displaying signal envelope; (c, c′) tracked bands of high reflectivity of IRHs (upper and lower edge as dashed lines, computed centre line as solid line).

Displaying the amplitude of the received signals, IRHs appear as bands of high reflectivity. Due to the above-mentioned length of the transmitter wavelet, the reflection horizons encompass multiple cycles of high amplitude (Fig. 2a). The uppermost 10–15 ns of the radargrams are masked by the arrival of the direct wave and cannot be used for analysis. Likewise, the lowermost 10 ns are blanked by the gain function.

The bands of high reflectivity become more clearly visible and focused when instead of the signal amplitude the signal envelope is displayed (Fig. 2b). The envelope display emphasizes zones of high reflectivity (strong instantaneous signal magnitude) vs zones of low reflectivity as it omits the zero-amplitude transitions. As the recorded radar data represent the real part of a complex signal, the imaginary part can be reconstructed by the Hilbert transform. Complex trace analysis yields the envelope, also termed the instantaneous magnitude of the complex signal, while the phase information is lost (Reference Mari, Glangeaud and CoppensMari and others, 1999).

The upper and lower edges of each band were tracked by logging trace number and two-way travel time (TWT) of the sample exceeding a threshold magnitude. The mean TWT of both values determines the centre line of the IRH, which coincides approximately with the maximum signal magnitude (Fig. 2c). The difference in the TWT between the upper and lower edge is referred to as the IRH width.

For later analysis, we converted the recorded radar signals from the time domain to the depth domain. The wave speed v of the radar signal in snow is determined by the permittivity, , where c is the speed of light in vacuum. We used the empirical relation of Reference Kovacs, Gow and MoreyKovacs and others (1995), , to link density and permittivity ε′ where the snow density ρ is in kg m⁻³. Finally, the results for reflector depth are half the TWT multiplied by the average wave speed to that depth. A typically observed mean density of 500 kg m⁻³ corresponds to a wave speed of 2.1 × 10⁸ m s⁻¹. During data acquisition in spring, the snow was completely frozen (as evident from thermistor data at Crawford Point). Effects of liquid water on the propagation of the radar signal can therefore be safely neglected.

The error in the depth determination of a tracked reflection can be split into range-dependent and range-independent contributions. The former include inaccuracies in the correction of arrival times, assumed to be on average no larger than one sample (0.13 ns or ∼1.5 cm), and inaccuracies in IRH tracking. Tracking IRHs is a partly subjective process; quantifying the related error is therefore difficult. We estimate the related error to be five samples (0.6 ns or ∼6 cm).

Range-dependent errors are introduced in the conversion of TWT to depth, due to deviations from the snow density–depth model. Because of the depth-averaging of the velocity profile, the conversion from time to depth domain is less prone to errors from variations in density on short vertical length scales at a single location. Rather, it is a deviation of the mean or bulk density along the GPR profiles that introduces depth errors. We analyze density–depth profiles from seven snow pits and firn cores within 100 m of T05. The standard deviation of bulk densities from one pit to another rapidly decreases with depth and falls below 4% at 1 m depth. For simplicity, we assume a maximum error of 4% in the density–depth model, constant over the whole depth range of the GPR. The related maximum error in TWT–depth conversion increases from 0 cm at the surface to >15 cm at 10 m depth. The sum of the three error contributions adds up to a maximum error in IRH depth of ≤10 cm in the upper 2 m to almost 25 cm in 10 m depth.

3.1. Comparison of IRHs observed at 500 and 800 MHz

Altogether, six IRHs down to about 10 m depth could be clearly identified and tracked at both frequencies, within the grid as well as along the 1 km sections. These are referred to as IRH-15 to IRH-95, with numbers indicating their approximate depth in nanoseconds of TWT near T05 (Table 2). To validate the measurements, 500 and 800 MHz results were first treated independently and later combined. The data at both frequencies reveal the same number of clearly identifiable IRHs, located at approximately the same depths with respect to the IRH centre line. IRHs obtained at 800 MHz appear on average 4.8% deeper than in the 500 MHz record. The differences in IRH depths increase with depth and range from 0.00 m for IRH-15 up to 0.52 m for IRH-95.

To investigate this further, we quantify the spatial variation in depth of each IRH in terms of its standard deviation (depth-SD). This property will serve as a first-order approximation of the small-scale spatial variability in snow accumulation. We consider it to be a measure of the representativeness of layer depths derived from individual point measurements within a certain area. Applying it to the comparison of 500 and 800 MHz data, it turns out that the depth-SD is on average 7.0% smaller for the 800 MHz than for the 500 MHz measurements for corresponding IRHs. No trend with respect to depth is present for the ratios of the depth-SD at both frequencies. Regarding differences in IRH widths, the 800 MHz survey yields slightly larger values: +6.0% on average. The observed differences in both IRH depths and widths can be considered small. They are most likely caused by variations in constructive and destructive interference from inhomogeneities because of the different frequencies and wavelets. The small differences indicate that IRHs obtained at 500 and 800 MHz likely arise from the same cluster of physical subsurface features.

Unless explicitly stated otherwise, we use mean values of both 500 and 800 MHz measurements for the discussion and interpretation in the rest of this paper. We will focus on the upper edges of IRHs, as they were the easiest to identify, usually occurring below a zone of low reflectivity and thus providing high tracking accuracy. The physical explanation for this observation is discussed in section 3.2.

3.2. Relation of GPR reflectivity to snow stratigraphy

Changes in density dominate the reflectivity of the firn pack, while changes in conductivity can be neglected. We now investigate which density features could cause the observed IRHs. The upper edge of IRH-15 is located at 1.48 ± 0.10 m. It corresponds well with the depth of the LSS (end of 2003 summer surface) of 1.43 ± 0.04 m (Reference ParryParry and others, 2007). We therefore consider it likely that IRH-15 arises from strong density contrast at the LSS formed during summer 2003. The IRH width and associated variability of IRH-15 is significantly smaller than for the deeper IRHs. We explain this with the time of data acquisition. As the data were acquired prior to melting conditions in summer, the overlying snowpack has not yet experienced metamorphosis through meltwater intrusion from above. The significantly larger reflector widths of deeper IRHs suggest that these may be influenced or even arise from ice-layer clusters formed during multiple summer seasons.

The depth range below IRH-15, down to approximately 4 m, is characterized by very complex and strong reflections. Two IRHs have been tracked in this depth range. The upper is IRH-25, at and below 2.29 m depth at T05. IRH-25 is laterally discontinuous and partly splits into two separate IRHs. IRH-35, at and below 3.45 m depth, is laterally well confined towards larger depths, where only weak signal returns are possible.

The AWS data at Crawford Point and the stratigraphic observations at TO5 in late summer 2004 provide an archetype for the development of the snowpack in warm summers, useful for interpreting deeper reflection patterns. The conditions during the melt period 2004 resemble those in the wet snow zone rather than the percolation zone. This is the year (2004) for which the AWS record yields the warmest summer months (June–August) with a mean temperature of −4.7°C, a high PDH factor and the highest initial snowpack temperature at the onset of melt (−15.9°C). Accordingly, after the melt season, discontinuous ice layers were found at all depths of the snow pits (Reference ParryParry and others, 2007).

Similar conditions can also be expected for the years 2003 and 2002, for which the AWS data show clearly above-average summer temperatures. An example is the 0.10–0.15 m thick ice layer at 3–3.5 m depth (Reference Scott, Mair, Nienow, Parry and MorrisScott and others, 2006; Reference ParryParry and others, 2007), which is probably related to the extreme summer melt in 2002. Strong percolation of surface-generated meltwater in 2002 eventually accreted as ice on top of the underlying summer surface. We therefore suggest that IRH-35 partly arises from ice which probably formed during summer 2002 on top of LSS 2001. The complex appearance of IRH-25 above IRH-35 may be a consequence of extensive downward meltwater flow and soaking of the previous accumulation layer during summers 2003 and 2002. IRH-25 could therefore represent both the 2002 summer surface and lower clusters of solid ice structures as well as the ice accreted on top during the following summer.

The 4–5.5 m depth range shows weak signal returns in the GPR data compared to the depth range above. The AWS data for the years 2000 and 2001 show significantly colder summer temperatures expressed in low PDH values as well as colder initial snow temperatures. The winter accumulation during this time period is expected to be less affected by meltwater infiltration and refreezing. To date, an independent age–depth model is not available to confirm these interpretations. However, we can utilize earlier estimates of accumulation to determine a tentative age–depth model.

3.3. Continuous IRHs

The detection of continuous IRHs is encouraging for deriving accumulation characteristics. Interestingly, other studies did not observe continuous IRHs. Reference Scott, Mair, Nienow, Parry and MorrisScott and others (2006) performed measurements at the same site in late summer 2004. Their data, collected along profiles of length 100 m, did not reveal continuous near-surface IRHs.

Three factors are likely to cause this discrepancy. First, they used a ground-based step-frequency radar with a centre frequency of ∼13 GHz, one set-up with a bandwidth 1 GHz, and in a second operational mode with a bandwidth of 8 GHz. The higher resolution of their system, compared to that used in the present study, causes higher volume backscatter and a reduced penetration depth. Laterally coherent reflections from layers or similarly layered sequences are therefore less likely to be observed. Second, the trace interval of 0.5 m up to almost 1.5 m chosen by Reference Scott, Mair, Nienow, Parry and MorrisScott and others (2006) is very large with respect to the high antennae resolution. It is therefore possible that despite a small antenna footprint, laterally consecutive first Fresnel zones did not overlap, therefore disabling laterally coherent reflections. In comparison, we used a trace interval of 0.1 m combined with a significantly lower resolution. Future experiments in comparable areas should take this into account. Finally, the complexity of the radar returns is further enhanced by the change in physical snow properties in the transition from pre- to late-summer conditions.

3.4. Age–depth model

Cumulative mass within a depth section divided by the corresponding time of accumulation yields the accumulation rate. We can therefore estimate the age of the IRHs by dividing the ice-core-derived cumulative mass above a particular horizon, obtained by integrating the density down to this horizon, by an independently determined mean accumulation rate at the site (available from earlier studies). We focus on the upper edges of IRHs, interpreted to represent LSSs for which the overlying winter layer is relatively undisturbed by ice accreted on top of the LSS.

Reference Fischer, Wagenbach, Laternser and HaeberliFischer and others (1995) derived accumulation rates along the EGIG line based on isotopic and chemical analysis of shallow firn cores, drilled in 1992. They found the annual mean snow accumulation over an 8 year period at T05 to be 46.0 ± 10.2 cm w.e. We attribute the variation of 22% to interannual variability of accumulation. Guided by the hypothetical time-span of accumulation, we appoint age estimates to the upper edges of IRHs (Table 2). This yields mean accumulation rates of 37.6–63.6 cm w.e. for periods enclosed by adjacent horizons. The mean accumulation rate over the period 1994–2004 is found to be 49.6 cm w.e. We estimate a temporal variability of 9.6 cm w.e. by the standard deviation of the six values, corresponding to the six IRHs. This value has to be considered as a lower limit for the interannual variability, as three of the periods average over 2 and 3 years, respectively. The relative uncertainty in the time since accumulation, and consequently the age estimate, amounts to roughly 30%. The absolute uncertainties therefore increase significantly with depth.

3.5. Scenarios for IRH formation

The stratigraphy in the percolation zone is very complex and its interpretation is based more on similar layered sequences than on the identification of specific layers (Reference BensonBenson, 1962). The initial processes for formation of a lateral continuous IRH occur at the snow surface mainly during the summer season. For the subsurface processes, which finally lead to amplification of the initial characteristics to produce strong IRHs in the upper ∼10 m, we propose the following scenarios.

In general, the IRHs are based on the integral effect of solid ice clusters on GPR signals at 500 and 800 MHz, as individual lateral continuous ice layers may be absent. Occasional melt–refreeze cycles during summers with little to moderate surface melt form an ice crust at the surface. Some meltwater infiltrates the snowpack and refreezes in the upper fraction of the layer of the previous winter’s accumulation, immediately below the summer surface. Once this layer is submerged and overlain by new accumulation, these features probably appear in the GPR record as a continuous IRH with a defined upper boundary and a rather diffuse lower boundary consisting of isolated signals of high reflectivity. Extensive breakthrough of percolating meltwater during warm summers with severe surface melt produces isolated ice lenses throughout the snowpack, but percolating meltwater may accrete as thick ice on top of the LSS.

This was the case in late summer 2004 where accretion of ice was observed at the bottom of a snow pit at T04 on the EGIG line, at an elevation of 1860 m, just 80 m lower than T05 (Reference Scott, Mair, Nienow, Parry and MorrisScott and others, 2006). Because of the strong density contrast between the porous snow matrix and the solid ice layer, this constitutes an inhomogeneity of high reflectivity in the GPR record. However, due to its proximity to the underlying LSS, this feature is barely separable from the LSS below. In case of persistent warm temperatures, extreme surface melt and subsequent meltwater percolation could cause a wetting front to propagate from the surface downwards. This may lead temporarily to snowpack conditions similar to those in the wet-snow zone at lower elevations. It can then form a unit of uniform wetting in an upper fraction of the snowpack, containing a complex mixture of small but numerous ice structures surrounded by coarse-grained snow crystals. As the complete layer thickness down to the underlying LSS is affected, this unit likely appears in the GPR record as a band of high backscatter with no defined lower boundary. Therefore, the continuous reflection from the LSS may be lost as a result of high clutter.

If no melting takes place during the summer, the conditions in the uppermost snow layer resemble those of the dry snow zone. The density, and thus the dielectric contrast, of layers originating from different accumulation events is small compared to normal percolation conditions. Only the formation of distinct layers of hoar or wind crust could provide a large enough density contrast to produce continuous IRHs. However, if this is not the case, then no continuous IRH might be observed with GPR, making it difficult to separate accumulation from subsequent years.

3.6. Spatio-temporal variation of accumulation

Here we discuss the spatial distribution of IRHs and the underlying accumulation pattern on different spatial scales. The results are first discussed independently of an age–depth scale in terms of depth and cumulative mass. We then apply the age–depth model established above to discuss the variation of accumulation in the larger regional context.

Over the area of the 100 m × 100 m grid, IRHs appear as interfaces aligned parallel to the reference surface (Fig. 2). The depth distribution of each IRH is basically of Gaussian shape. The depth-SD of the IRHs varies from 0.10 to 0.21 m, with a mean of 0.17 m, with the smallest value for the uppermost IRH-15, but shows no further systematic trend with depth (Table 2). As all depth-SD values are of the same order of magnitude for all horizons, we assume that a consistent spatial gradient in snow accumulation cannot be detected at scales of 100 m at this location. For this length scale we can therefore estimate the accumulation variation to be 8.5 cm w.e., derived by multiplying the mean depth-SD with a typical density of 500 kg m⁻³.

To identify a possible spatial gradient in snow accumulation over the area enclosed by the 1 km × 1 km square, we investigate the distribution of IRH-50. This IRH produces strong and relatively undisturbed return signals, providing reliable results. We compute inverse-distance weighted mean values of cumulative mass for measurements within a search radius of 100 m from mid- and end-points of the 1 km sections. Deviations are expressed as percentages of the average mass along all 1 km sections, indicated as bars in Figure 3. The analysis reveals increasing accumulation towards the northeast, a trend not recognized over the area of the 100 m × 100 m grid. Values range from −6.9% at the mid-point of the line connecting point T05 and point T05-S4 in the southwest to +6.1% at point T05-E4 in the northeast. IRHs along the 1 km sections show an increasing depth-SD with depth, varying between 0.11 m for IRH-15 and 0.34 m for IRH-95, with an overall mean of 0.22 m.

Fig. 3. Digital elevation model with 5 m contour lines, based on laser scanning during the ASIRAS campaign (Reference HelmHelm and others, 2007). Also shown is the location of the 1 km × 1 km square and the variation of cumulative mass above IRH-50 along its profiles. Bars visualize percentage deviations of inverse-distance weighted cumulative mass within a search radius of 100 m with respect to overall mean.

Spatial gradients in accumulation of such magnitude have been observed in many areas on ice sheets. They usually represent local effects caused by topographic undulations on small scales (few kilometres) that affect wind-driven redistribution of snow (Reference King, Anderson, Vaughan, Mann, Mobbs and VosperKing and others, 2004). To investigate if this process is also of relevance here, we use the new digital elevation model of the survey area, acquired during the ASIRAS campaign in spring 2004 by laser scanning (Reference HelmHelm and others, 2007).

The topography data indicate that our survey area lies on a plateau, with a nearby slope break to the southwest (Fig. 3). ASIRAS detected a slope-dependent, bimodal distribution of the uppermost-layer snow thickness along a 2.7 km crossover flight of T05 (Reference HelmHelm and others, 2007). A layer thickness of 1.30 ± 0.05 m was found on the slope (eastwards-rising elevation of 1920–1940 m), and 1.50 ± 0.13 m on the plateau (elevations of 1940–1945 m). Our 1 km GPR profiles are situated on the plateau, close to the edge of the area covered by ASIRAS, and obviously fall in a transition zone between the two accumulation modes. The ASIRAS-based layer thickness for the plateau corresponds very well with the depth of IRH-15 and the depth of the LSS from snow-pit analysis. While the strong response of the LSS causes complete power loss of the ASIRAS signal for the depth range below, we showed above that the GPR system is capable of detecting continuous IRHs down to 10 m depth. This enables us to investigate the spatial variability for earlier time periods.

To this end we use selected 800 MHz sections oriented parallel and perpendicular to the EGIG line. IRHs appear to be parallel to the surface along the section perpendicular to the EGIG line. The section following the EGIG line from T05 in a direction east-northeast reveals that IRHs are clearly dipping in that direction (Fig. 4a). To analyze the temporal variability of the spatial gradient in snow accumulation, we determine the thicknesses of units enclosed by adjacent horizons (see unit numbering in Fig. 4). The units were normalized by their corresponding mean thickness. We then performed a linear regression for all units. The fitted lines pass through the centroid (x mean, y mean; Fig. 4b). It is evident that the trends along this GPR profile are well represented by linear ap-proximations, indicating that the underlying gradient in snow accumulation is also linear. In general, the gradient appears less strong for deeper and older units than for younger units. The positive gradient in accumulation towards the east-northeast persists for most of the accumulation periods covered by our data.

Fig. 4. (a) IRHs along the 1 km profile east of T05 (see inset in (b)). The positions of IRH centre lines and corresponding linear regression lines are plotted. (b) Normalized thickness for layers enclosed by adjacent horizons, with profile centres crossing the mean (normalized unit thickness ≡ 1).

Layer thicknesses increase by about 10% from southwest to northeast for units 1, 2, 3 and 5. The trend is reversed for the deepest unit 6, for which normalized thickness is about 3% greater at the northeast than at the southwest end. The thickness of unit 4 is almost constant. The observed weakening of the gradient with depth could be explained by the fact that we neglect ice advection. The surface velocity of the ice flow at the survey site is about 100 m a⁻¹. Ten-year-old firn originates about 1 km upstream from the present location, further up on the plateau, where ASIRAS detected little variability in snow accumulation (Reference HelmHelm and others, 2007).

The age estimate of the six IRHs enables us to convert the observed IRH depths and cumulative mass discussed above to accumulation values. The previous year’s accumulation (summer 2003–spring 2004) yields 63.6 ± 4.3 cm w.e. This value is about 40% larger than previously derived accumulation rates (Reference Fischer, Wagenbach, Laternser and HaeberliFischer and others, 1995), but is in agreement with the snow-pit derived accumulation of 60.5 ± 3.4 cm w.e. (Reference ParryParry and others, 2007). It should be noted that this is a preliminary lower-limit result, since the mass-balance year was not fully completed by the time of data acquisition. The observed northeast gradient in snow accumulation is of the order 5.0 cm w.e. km⁻¹, roughly 10% of the mean per km. This value is one order of magnitude larger than the regional trend on a 100 km scale on existing maps of accumulation rates over Greenland (Reference Bales, McConnell, Mosley-Thompson and LamoreyBales and others, 2001). We explain this discrepancy as the result of the difference between local and regional elevation trends influencing the katabatic wind speed and thus post-depositional redistribution of snow. Whereas the regional slope inland of the study site is of the order 1–5 × 10⁻³ (i.e. 100–500 m over 100 km; Fig. 1), local slopes in the plateau-like area of the 1 km square amount to ∼10⁻² (Fig. 3). The ASIRAS observation of further decreasing accumulation on the larger slope in the west confirms this relation.

The redistribution of mass by percolation in years of strong melt could cause a higher uncertainty for the estimated accumulation values. Percolating meltwater during summer 2003 could have caused ice accretion on top of the LSS of 2002 (IRH-25). This mass would therefore contribute to the value computed for the previous period 2001/02, but would be lacking in the period 2002/03. Thus, the minimum accumulation rate of 37.6 ±12.5 cm w.e. for the period 2002/03 (Table 2) is likely an underestimate of the full balance year. Mass loss by percolation may also apply for the the period 2001/02. If the amount of loss by percolation and gain by internal accumulation were comparable in both periods, the downward loss of mass and the gain from above would be compensated.