Spatial mapping of multi-year superimposed ice on the glacier Kongsvegen, Svalbard

Abstract

Ground-penetrating radar (GPR) and satellite ERS-2 synthetic aperture radar (SAR) are used to map the thickness and extent of the superimposed ice (SI) zone on the surge-type glacier Kongsvegen, Svalbard. GPR imagery shows sub-horizontal SI layers lying unconformably above a discrete boundary. Below this boundary, the ice has a GPR signature similar to that of ice further down-glacier in the ablation zone. This boundary is posited to represent the closing of crevasses that were created during the last surge of Kongsvegen in ∼1948. Open crevasses would have interrupted the formation of sheet layers of SI due to efficient vertical drainage of the snowpack. Aerial photographs suggest that the crevasses closed sometime in the period 1956–66. A classified SAR image from 2003 is used to delineate the extent of the SI zone. The SI extent in the SAR image agrees well with the SI zone mapped by GPR. Using the SI spatial depth distribution, we estimate the mean annual accumulation of superimposed ice to be 0.16 ± 0.06 m w.e. a⁻¹ (locally up to 0.43 m a⁻¹ w.e.). This corresponds to ∼15–33% of the local winter balance and ∼5–10% of the total winter balance measured since 1987.

Introduction

Retention of water in the accumulation zone of glaciers and ice sheets is of interest in the light of predicted climate change (Bøggild and others, 2005; Solomon and others, 2007). The water derives either from melting snow and ice or from rain. Refreezing occurs either when water reaches a cold underlying layer, or when atmospheric cooling is sufficiently strong to freeze the wet snow. The former mechanism generally occurs during spring and summer while the latter is associated with autumn (Wadham and Nuttall, 2002; Obleitner and Lehning, 2004). If the water refreezes above the previous summer surface, forming a continuous mass, it is denoted ‘superimposed ice’ (SI) (Paterson, 1994). The region where there is an annual increment of SI exposed at the surface is further defined as the ‘SI zone’ and will have an area smaller than the region where SI forms. SI formation is thus different from ‘internal accumulation’ which we define as the retention of water below the previous years’ summer surface.

The amount of SI formed depends on the amount of available water, the energy balance of the snowpack, the capability of water to percolate through the snowpack and the local topography (which dictates whether water channels or collects). SI formation is spatially and temporally variable, with the result that SI tends to form in events that create a signature of interleaving bubbly and clear ice layers (König and others, 2002; Hooke, 2005). Therefore, layers might not necessarily represent annual or seasonal signatures but rather a collection of short-term events.

Since the winter balance is generally small on arctic and low-precipitation glaciers, SI accumulation can be a relatively large percentage of the total accumulation. The formation of SI is therefore important when estimating the mass balance of these types of glaciers (Fujita and others, 1996; Wadham and Nuttall, 2002 and references therein) and has implications when estimating glacial run-off (Bøggild and others, 2005).

In this study, ground-penetrating radar (GPR) and satellite synthetic aperture radar (SAR) imagery is used to investigate the spatial depth distribution of accumulated SI and spatial extent of the SI zone on Kongsvegen, a 25 km long glacier situated in northwest Spitsbergen, Svalbard (Fig. 1). Kongsvegen is an arctic glacier, and the site of previous SI studies (König and others, 2002; Obleitner and Lehning, 2004; Langley and others, 2007). Obleitner and Lehning (2004) modeled SI formation and melting at the site of an unattended automatic weather station (AWS) at stake 6 on Kongsvegen (Fig. 1), using the one-dimensional, finite-element snowpack model ‘SNOWPACK’ (Bartelt and Lehning, 2002; Lehning and others, 2002a, b). The model successfully explained the temporal evolution of the snow-pack, matching the snowpack temperature observations reasonably well, but excess water needed to be added to the model to generate the amount of SI measured at the study site. Obleitner and Lehning concluded that the two-dimensional flow of water in the snowpack is important to explain the amounts of SI seen on Kongsvegen.

Fig. 1. Map showing Svalbard and the location of Kongsvegen and midre Lovénbreen. Mass balance stakes are marked with dots and numbers.

Previous GPR studies on Kongsvegen focused on the thermal regime, firn–ice transition and firn make-up (Björnsson and others, 1996; Kohler, 2002; Kohler and others, 2003; Pälli and others, 2003). Recent work by Langley and others (2007) revealed differences in the GPR signature of the different glacier zones of Kongsvegen. They also showed ice with continuous sub-horizontal layers lying unconformably on ice with no apparent layering, as illustrated in Figure 2. The GPR signature of the latter ice is identical to that exposed at the surface down-glacier in which a distinctive repeated ‘herringbone’ pattern is apparent, with a spacing of 50–150 m. This discrete unconformable contact between the layered and non-layered ice is not seen in GPR imagery of SI on the nearby glacier, midre Lovénbreen (Wadham and others, 2006). The midre Lovénbreen behavior is what one would expect for a glacier in which SI is produced and preserved year after year. The appearance of a well-defined basal boundary of SI on Kongsvegen suggests that it has only accumulated for a limited period. This implies that on Kongsvegen, as expected (Hooke, 2005), the maximum multi-year SI depth occurs in close proximity to the firn–SI boundary zone, where maximum annual SI formation is expected, and that the depth subsequently decreases down-glacier towards the equilibrium line.

Fig. 2. (a) Hilbert-transformed 500 MHz GPR profile along center line (mass balance stakes 5–8 in Fig. 1) showing glacier ice, SI and firn. (b) Amplitude format close-up of SI zone. SI has a clear layered appearance, in contrast to the underlying glacier ice, which shows a distinctive repeated herringbone pattern with a spacing comparable to the typical crevasse spacing (50–150 m). A two-way travel time (TWT) of 100 ns corresponds to ∼9 m in the glacier ice and SI and ∼11 m in the firn.

Kongsvegen surged sometime just before 1948 (as seen in aerial photographs). This surge resulted in extensive surface crevassing, throughout today’s SI and ablation zone. Until the crevasses closed, they would have provided an effective vertical draining path for meltwater. This would ultimately limit the thickness of SI formation on the areas between crevasses and, more importantly, prevent the formation of SI layers as continuous sheets, at least until the crevasses had closed again.

The discrete boundary creates an opportunity to estimate the depth distribution of multi-year SI that has accumulated since the crevasses closed. Combining the GPR-derived spatial depth distribution with the SI zone extent determined from satellite SAR, an average annual SI zone accumulation rate is determined and compared to the annual mass balance measured since 1987.

Field Site

Kongsvegen covers an area of ∼100 km², ranging from 0 to 800 m a.s.l. The average equilibrium-line altitude (ELA) is ∼500 m a.s.l., with an average net mass balance of ∼0.11 m w.e. a⁻¹ (Lefauconnier and others, 1999; Hagen and others, 2003). The glacier is polythermal, with a cold surface layer partly reaching the bed in the ablation area (Hagen and Sætrang, 1991; Björnsson and others, 1996). Following the surge in ∼1948, Kongsvegen has been in a quiescent stage, with horizontal velocities less than 5 m a⁻¹ (Melvold and Hagen, 1998). Based on SAR imagery (König and others, 2002), Obleitner and Lehning (2004) estimated that 35% of the glacier area is covered by SI, comparable to the area covered by bare ice (36%).

Data And Methods

Satellite SAR

A European Remote-sensing Satellite-2 (ERS-2) SAR image from 23 March 2003 is used to determine the extent of the SI zone. Image acquisition details and processing methods are described by König and others (2004).

The SAR image shows three different backscatter zones: low, medium and high. The high-backscatter zone (> −2.5 dB) corresponds to the firn area, and the low-backscatter zone (<−8.5 dB) to the non-SI-covered ‘normal’ glacier ice area. The medium backscatter zone (>−8.5 and <−2.5 dB) has been interpreted as the surface area covered by SI (König and others, 2004).

Using these simple thresholds, the image is classified into three zones: ice, SI and firn (Fig. 3). Areas with steep slopes have been removed to avoid backscatter signals that are caused by differences in surface geometry rather than firn or ice properties.

Fig. 3. Classified SAR image showing the three classes: firn, SI and (glacier) ice. Steeper slopes are excluded and shown white. The SI depths along GPR profiles are indicated with colors. GPR profiles where no SI has been detected are marked with black dots.

GPR

A Geophysical Survey Systems Inc. (GSSI) SIR-2 control unit and a GSSI 3102A shielded 500 MHz antenna were towed behind a snowmobile in a regular grid on the glacier surface in April and May 2003. Traces were triggered either two or four times per meter using an odometer wheel. Survey lines were positioned using differential GPS (global positioning system) yielding positional accuracies of ∼0.1 m. The position of each trace is determined by distributing the traces equally along the survey line between calibration points (points where trace and GPS positions coincide). The accuracy of the trace position is therefore less than the GPS accuracy (due to variations in triggering distance) but should be in the range of 1 m. Nonetheless, this accuracy is nearly an order of magnitude more precise than the SAR pixel resolution of 30 m.

The GPR data have been post-processed by stacking, bandpass filtering and applying time-varying gain. The first reflection on the GPR images is from the bottom interface of the layer of winter snow upon which the snowmobile is driving (Fig. 2). Up-glacier from stake 7 (and below the winter snow), the GPR signature of the firn is characterized by high backscatter due to inhomogeneities in the firn (such as ice layers and lenses) as well as the presence of liquid water at depths greater than ∼15 m (Brandt and others, unpublished data). SI layers are located between stakes 6 and 7, which are characterized by a weaker backscatter than the firn, and appear as prograding sub-horizontal deposits that thicken up-glacier.

Down-glacier from stake 6 is an area composed of ‘glacier ice’, with no apparent continuous layering. However, this ice (when seen in the along-flow direction) has the distinctive herringbone structure mentioned above (Fig. 2). This is a wavy repeated pattern with spacing of ∼50–150 m, and wave crests aligned approximately cross-glacier. The spacing, appearance and orientation of the pattern suggests that it is the result of crevassing during the surge (leaving scars), but more detailed study is required to show the exact cause of the radar signature.

We calculate total SI zone accumulation along the GPR profiles by digitizing both the snow–SI and SI–glacier ice boundaries. The GPR-derived SI depths (excluding snow) are shown overlaying the classified SAR image in Figure 3.

The contact between snow and ice was digitized using an automatic layer-following routine (Heiberg, 2004). The two-way travel time (TWT) to the lower SI boundary is more irregular and thus was hand-digitized. To better visualize the SI layers, the data were Hilbert-transformed. The digitized TWT to the lower SI boundary was converted to depth based on GPR common-midpoint (CMP) profiles at stake 6.5 (midway between stakes 6 and 7 in the center of the SI zone) which revealed typical electromagnetic (EM) wave velocities in the SI ice from 1.7 to 1.9 × 10⁸ m s⁻¹. Conservatively, we assumed an electromagnetic (EM) wave velocity of 1.8 × 10⁸ m s⁻¹ ± 10% for SI ice and 2.24 × 10⁸ m s⁻¹ for snow, corresponding to a snow density of ∼400 kg m⁻³ (Eisen and others, 2002).

The contact between snow and ice has a well-defined reflection which is readily digitized without Hilbert-transforming the data. Snow depth is digitized to within ±1 ns. A bulk snow density variation of 300–450 kg m⁻³ gives an estimate of the uncertainty in the radar-wave velocity of roughly ±5%, leading to a snow depth uncertainty of ∼0.2 m. The total snow depth is therefore believed to be accurate to ±0.15 m.

The accuracy of the GPR-derived SI thickness depends on GPR resolution, the assumed EM wave velocity and the digitization accuracy. The assumed EM wave velocity should be accurate to ±10%. The minimum vertical resolution for an impulse GPR (i.e. the minimum distance required to differentiate two layers) is ∼1/4 of the wavelength. For our ∼500 MHz center-frequency system in ice and firn this is ∼0.10 m. The accuracy of the digitized Hilbert-transformed data is estimated to be ±5 ns, resulting in a depth accuracy of ±0.4 m. The EM wave velocity uncertainty (±10%) gives an additional error that increases with depth, varying from close to zero at the snow–ice interface to ±2 m at ∼20 m depth. The interpretation accuracy is difficult to quantify, but is probably similar in magnitude. Profiles acquired along the glacier flow direction are generally easier to interpret than the cross-glacier profiles, due to the orientation of the herringbone pattern. In two of the lower GPR cross-profiles, the contact between layered SI and underlying ice is therefore hard to differentiate. This affects the comparison of the lower SI zone boundary as inferred from GPR to that obtained with the SAR image. However, these profiles represent only a relatively small part (∼500m) of all the profiles, resulting in only minor uncertainty when calculating the annual average SI accumulation.

The total error for the SI thickness is therefore roughly estimated to be ±2 m for the deepest parts, decreasing to ±0.4 m close to the snow–ice interface.

DEM

The latest digital elevation model (DEM) available for Kongsvegen is based on aerial photographs taken on 18 August 1995, at a height of 6000 m a.s.l. (scale 1 :30 000), and has a pixel resolution of 10 m. Artifacts were found in the DEM, caused by both the stitching together of aerial photographs and problems of low-contrasting snow surfaces (which makes photogrammetry difficult) Therefore, a low-pass filter using linear interpolation was applied to the DEM to eliminate the discontinuities. The DEM was then re-sampled to contain a pixel resolution of 25 m.

Surface drainage

The GPR images show that SI thickness decreases close to surface meltwater channels. To determine if there is any large-scale connection between the drainage and the SI thickness, surface meltwater pathways are digitized over the SI area. The major melt channels were identified using aerial photographs, while the DEM was used to estimate the most likely drainage patterns.

Major melt channels were identified in a 1 m × 1 m resolution orthophotograph from 1995. Most of the digitized channels were wider than 3 m and thus covered several pixels. The digitized positioning of the channels is believed to be about the same as the pixel positioning accuracy, which is conservatively estimated to be 5 m.

Older aerial photographs from 1956, 1961, 1966 and 1990 were used to determine when the surge crevasses closed and the current surface drainage pattern started to evolve. Positioning of surface channels in these images is not as accurate as in the 1995 image since orthorectified images are not available.

While the 1995 orthophotograph shows the locations of the major streams, some channels were only identified for short distances due to image stitching, snowpatches on the surface and contrast problems. To complement these observations of the largest-scale drainage network, we use the DEM gradients, water mass conservation and a simple melt-rate approximation (based on the mass-balance measurements since 1987) to model the most likely water pathways and the relative differences in water flux through the different pathways. The model chosen comes from Lüthje and others (2006) and assumes an impermeable ice surface and that the water flow (vector) velocities, u, are governed by a simple gradient relation

(1)

where ᐁ ≡ (∂/∂x, ∂/∂y), x and y are the horizontal spatial coordinates, Ω is the hydraulic conductivity, assumed constant in space and time, and is the equipotential surface given by ψ = H + h, where H is the underlying ice topography and h the water level.

By discretizing the space and time domain using a finite-difference time-domain scheme, the water level in each space element is determined by the local melt rate and the local horizontal flux gradient (difference in flux going in and out of the box element) which is directly related to Equation (1) by multiplying u by the spatial discretization size and water level, h.

Water levels, h, are initially set to zero for the first time-step, and a daily melt-rate pattern as a function of elevation was estimated using the summer mass-balance measurements (1987 to present, as a function of elevation) and assuming an ablation season length of 60 days. Thus, given the elevation for each space element in the model, a melt rate was derived for that element by finding the corresponding summer balance for that elevation divided by 60 days. Further, the melt rate was assigned to zero in the firn area, implying that all melt in the firn percolates downward as internal accumulation. Meltwater that reaches the downslope boundary is removed, and closed boundaries are applied to the other three sides in the model. The extent of the firn area was determined as the area classified as firn in the SAR image. Melt rates were also assigned zero for surface slopes larger than 10%, which are mainly non-glaciated areas. Hydraulic conductivity, St, was set to 1 m s⁻¹ as no qualitative differences resulted from experimenting with other values.

In summary, this simple model does not attempt to model the flow in slush and snow in any detail. Instead, the most likely large-scale flow pattern is determined in order to estimate the water flux that can be expected in different areas.

Calculation of annual average SI accumulation

The spatial depth distribution of SI on Kongsvegen was estimated by combining the SI depths measured by GPR and the area classified as SI in the SAR image. Tributary glacier areas not covered by GPR profiles were excluded from this analysis. The GPR soundings were linearly interpolated to a regular grid. Positions outside the SI area (determined from the SAR imagery) were added (in the gridding) to ensure 0 m depth outside the SI area. The sensitivity of the gridding procedure was tested by altering the number of GPR depth soundings included in the interpolation and by altering the gridcell size. The cell size varied between 50 and 500 m, and the number of GPR soundings (included in the procedure) varied between 50 and 800. The resulting volume was insensitive to choice of cell size, and the average SI depth derived from the re-gridded data corresponded closely to the average depth in the digitized profiles. To estimate the mean annual SI accumulation during the period, the average SI depth is divided by the number of years since the SI started building up in continuous layers again, after the surge ended.

Results

GPR and SAR data

The firn–ice transition found using GPR coincides well with the transition from medium to high backscatter depicted in the SAR image (Fig. 3). Similarly, the area classified as SI in the SAR image corresponds well with the area of SI found from the GPR profiles. This suggests that SAR dB thresholds provide a crude but useful estimate of the SI zone.

Annual average SI accumulation

The 1956 aerial photographs show that some parts of the SI zone had already developed a surface drainage system reasonably similar to the one it has today. In 1956 and 1961, some crevasses were still open in the lower part of today’s SI zone, in which these areas were free of surface streams. Also visible in the image are old crevasse scars in other parts of the SI area. By 1966 a discernible surface drainage pattern had evolved over the entire SI zone. It was therefore concluded that the multi-year accumulation of SI started sometime between 1956 and 1966. The midpoint of these two years, 1961, is thus assumed as the starting date for closed crevasses and the initiation of sheet SI accumulation.

Dividing the average SI depth with the ∼42 year time frame (1961–2003), the average annual accumulation of SI is estimated to be 0.16 ± 0.06 m w.e. a⁻¹ (a SI density of 800 kg m⁻³ is assumed, taken from the GPR CMP-derived EM wave velocity of 1.8 × 10⁸ m s⁻¹ in the conversions discussed above and validated by Eisen and others (2002)). This corresponds on average to ∼15–33% of the specific winter balance and ∼5–10% of the total glacier-averaged winter balance which has been measured by the Norwegian Polar Institute since 1987. However, in the deepest multi-year SI, the SI accumulation rate must reach up to ∼0.43 mw.e. a⁻¹ (note the low horizontal ice velocity, less than 5 m a⁻¹,which affects these calculations).

Long-term SI depth distribution

The multi-year SI depth distribution accumulated over ∼42 years shows that the deepest SI is located close to the firn line and thinnest close to the ELA, forming a thinning ‘wedge’ shape down-glacier. This pattern, described by Paterson (1994) and Hooke (2005), is a consequence of the refrozen meltwater being less likely remelted at higher elevations and run off the glacier. Our measurements also show that the SI depth varies considerably within the SI zone (Figs 2 and 3). In addition, SI thickness decreases close to surface meltwater channels as seen by plotting the digitized surface channels and the modeled water drainage pattern together with the SI depth (Fig. 4). The SI depth is consistently shallower in the vicinity of model water pathways. These modeled pathways are mainly influenced by the topography, determined primarily by the DEM. However, the re-sampled DEM resolution is too coarse to resolve single meltwater channels, and, moreover, the simplicity of the model precludes obtaining a detailed drainage network. Model results are used to display likely positions for water pathways and where ponding can be expected due to low flow velocities.

Fig. 4. The log₁₀ (to enhance small-scale features) of the modeled water depth (grayscale), the SI depth along the GPR profiles (color) and digitized surface streams (×). Circles indicate places where the SI depth is distinctly shallower in connection with meltwater channels.

Local accumulation of SI is, to a large extent, indirectly dictated by the topography, which governs the water pathways. It is also apparent, though, that the long-term SI depth increases with elevation (Figs 3 and 4) and therefore must also be related to snow depth and distance to the firn line. This implies a strong relation between SI depth and the mass balance. Finally, SI depth also shows a marked cross-glacier gradient, with more SI on the southern side of the glacier. Despite the trends in SI depth, clear mismatches and very abrupt changes are apparent in SI depth (∼5 m) over short distances (∼20 m) within the GPR data (see, e.g., the northeast end of the line passing through stake 7 in Fig. 4). The cause of these very abrupt changes is still not well understood and can, presumably, vary from case to case.

Discussion

Implications for mass balance

The uncertainty in the calculated average annual SI accumulation is mainly dependent on the uncertainty of the start time of the build-up and the estimated SI average depth. A 10 year uncertainty (start of accumulation in 1956 or 1966) gives the same magnitude of uncertainty (of average annual accumulation) as that in the calculation of the average depth.

However, the average annual SI accumulation of 0.16 ± 0.06 m w.e. a⁻¹ is comparable to the SI contribution to the mass balance reported by Wadham and Nuttall (2002) for the nearby glacier midre Lovénbreen (Fig. 1). They estimated that SI at an individual stake could account for 20% of the specific winter balance or 16–25% of the annual accumulation. Our estimate (15–33%) is in the same range for specific balance, but our estimate of SI contribution to the total winter balance (5–10%) is lower. Obleitner and Lehning (2004) reported up to 60 cm of SI accumulation (locally) during the summer of 2000, but most of it melted the following summer.

SI area mapped by GPR and SAR

The area covered by SI was more extensive in the classified SAR image than in the results from the GPR. This mismatch can be explained by choice of the SAR backscatter dB threshold value, by the threshold method in general or by GPR interpretation difficulties. In the lowest profiles, the SI extent interpreted from the GPR profiles could potentially have been confused with ice comprising a herringbone pattern (Fig. 2b). This would be most likely to have happened in the cross-glacier profiles, where it was difficult to resolve the discrete boundary between SI layers and the underlying glacier ice. More profiles in the ice-flow direction could resolve this uncertainty and more clearly determine the boundary.

The discrete unconformable contact between layered and non-layered ice suggests that GPR measurements can ideally reveal information of previous surges and be used to classify whether a glacier is of surge type or not. Whether or not a glacier that surges results in such a discrete boundary between SI and ‘glacier ice’ depends on the degree of surface crevassing in the SI zone, and how fast the glacier develops a reasonably impermeable surface on which SI can start to build up. A discrete boundary of the type identified on Kongsvegen may not develop on all surge-type glaciers, but, if one is observed, the glacier can most probably be classified as surge type. The frequency of this phenomenon still requires clarification.

Drainage model compared to observed surface streams

The modeled drainage pattern reached a steady state after ∼25 days. As the goal was to gain insight into the most probable water pathways at the end of the melt season, which is predicted by Obleitner and Lehning (2004) to be an important period for formation of SI, the temporal evolution of the drainage pattern was not of importance.

The three most important factors that determine the final model pattern are the water inputs used to drive the model, the assumed value of the snow hydraulic conductivity and the ice topography. A simple parameterization of constant daily melt rates was used, but more realistic daily-varying input signals should not significantly alter the steady-state drainage pattern. Testing the model with other values of Ω resulted in little influence on the water-pathway pattern. It would have been interesting to test different ice topographies to investigate the model’s sensitivity to the fact that the topography may evolve. However, only one DEM was available and thus no sensitivity test was performed.

Nonetheless, despite the drainage model’s simplicity and relatively crude parameterization, the steady-state drainage pattern corresponds well with the surface meltwater channels identified in the orthophotographs (Fig. 4). In order to emphasize the flow pattern, the steady-state water flow is plotted using a logarithmic scale. This is done to minimize deep water ponds deriving from inconsistencies in the DEM and enhance the small-scale differences in flow pattern. The hydrological conductivity is, in reality, not a single value for all types of flows as water will flow at different speeds, depending on the channel size and surrounding media. This causes non-perfect matching between the observed and modeled water pathways. Nonetheless, the pattern corresponds well with the drainage systems identified in the photographs from both 1990 and 1966. This also suggests that the positions of the water pathways have been fairly constant during most of the period investigated.

Factors governing the long-term distribution of SI

The long-term SI spatial depth distribution varies considerably within the SI zone (Figs 2 and 3), which leads to the question: what parameters best explain and govern the spatial distribution of multi-year SI depth? The main difficulty in answering this question is that many of the parameters suggested in other studies (i.e. Wadham and Nuttall, 2002; Obleitner and Lehning, 2004) correlate with each other, thus hindering a quantitative analysis using regression of observed and derived parameters. For this reason we restrict ourselves to stating that the long-term SI depth distribution shows a strong connection to: (1) the water pathways and thus topography; (2) elevation, snow depth and distance to the firn, all of which are highly correlated; and (3) cross-glacier gradient. This approach does not explain the physical reason for the differences in depth but instead confirms what Obleitner and Lehning (2004) and Wadham and Nuttall (2002), for example, have postulated.

Our approach of simply mapping the depth distribution differs from the studies of Wadham and Nuttall (2002) and Obleitner and Lehning (2004). These studies cover shorter time periods, are performed at smaller spatial scales and use physically based models forced by meteorological data to explain observed SI on a sub-annual timescale in order to determine the most prominent parameters. These studies concluded that SI formation in the short term is highly variable, both spatially and temporally.

Interestingly, while Obleitner and Lehning (2004) assert that SI accumulation is highly variable in space, the layers observed in our GPR profiles are trackable over several hundreds of meters, indicating that major events build up the long-term accumulation that occurs over larger areas.

It is shown that over long periods the water pathways, and implicitly the topography, explain some of the smaller spatial-scale variability, while elevation, snow depth and distance to firn, which all correlate to the mass balance, explain the larger-scale variability. Thus, essentially, the mass-balance gradient determines the overall wedge shape of the SI area, increasing in thickness in the up-glacier direction, while the water flow influences the local SI accumulation pattern.

Conclusions

The surge of Kongsvegen in ∼1948 makes it possible to map the spatial depth distribution of multi-year SI built up during a ∼42 year period using GPR. The multi-year SI mapped by GPR in 2003 agrees well with a zone classified as SI in a SAR image from spring 2003.

Over the ∼42 year period, the mean annual SI accumulation is 0.16 ± 0.06 m w.e.a⁻¹, but varies locally up to 0.43 m `w.e.a⁻¹. The mean annual estimate corresponds to ∼15–33% of the local winter balance and ∼5–10% of the total winter balance.

The long-term SI depth distribution shows a strong relation to (1) the water pathways and thus topography, (2) elevation, snow depth and distance to the firn, all of which are highly correlated with the mass balance, and (3) cross-glacier gradient. The first two confirm results by Obleitner and Lehning (2004) and Wadham and Nuttall (2002). Furthermore, the GPR layers in the SI zone are often trackable over several hundreds of meters, indicating that at least the major SI-forming events occur over large areas. Therefore, the high spatial variation found by Obleitner and Lehning (2004) might not be of importance for long-term accumulation.

Finally we suggest that the discrete unconformable contact between layered and non-layered ice ultimately could be used to reveal whether a glacier has surged and whether it is of surge type or not.