Introduction
Arctic glacier mass balances contain a record of current climatic change and, by extrapolation, should yield a warning of future changes (Reference Jania and HagenJania and Hagen, 1996). For the overall balance of the arctic ice masses, accurate information on snow accumulation is required. The traditional methods used, such as repeated stake measurements, snow probing, or stratigraphic studies in snow pits and firn-core studies, provide insufficient information about the spatial variability of snow accumulation.
The aim of this study is to resolve the spatial variability of snow accumulation on the upper parts of Nordenskjöldbreen on Lomonosovfonna, Svalbard (Fig. 1), using ground-penetrating radar (GPR) data and data from three ice cores along the radar profile. Ice cores provide excellent archives of past climate and environmental change, but data are limited to the drill site. Radar is routinely used to map spatial distribution of glacier properties for which there is an associated change in dielectric impedance, such as the interface between glacier ice and bed (ice thickness), the interface between cold and warm ice in polythermal glaciers, or internal layering in snow or ice. Layering in radar images of glacier snow or firn can be attributed to differences in ice density (Reference Paren and RobinParen and Robin, 1975), liquid-water content (Reference Macheret, Moskalevsky and VasilenkoMacheret and others, 1993), chemical composition, microparticle concentration (Reference DanielsDaniels, 1996) and crystal fabric (Reference Glen and ParenGlen and Paren, 1975; Reference Fujita and MaeFujita and Mae, 1994).
(a) Map of Svalbard and the locations of Lomonosovfonna, Longyearbyen, Ny-Ålesund and Hornsund. (b) Lomonosovfonna and Nordenskjöldbreen with the position of the deep core (10) and shallow cores 1–9 and Uranusfjellet. The radar profile from core 7 to core 10 is indicated by the black solid line.
In addition to being able to quickly image long profiles, radar, with its relatively wide beamwidth and large spatial footprint, images ice properties averaged over greater horizontal areas than do the equivalent data obtained from ice cores. Thus, a marker horizon resulting from a dielectric contrast, such as a bubbly ice layer or a sulfate-rich layer from a volcanic event, which may be absent in an individual core due to local accumulation variations, can nevertheless be detected in a radar image taken over the core site.
Using radar, data from cores separated from each other by distances of tens of kilometers can be correlated and the spatial behavior of snow layers, and thereby accumulation between layers, mapped between the cores. For example, Reference Siegert and HodgkinsSiegert and Hodgkins (2000) used 60 MHz airborne radar data from Antarctica to trace isochronous internal ice-sheet layering between Lake Vostok and Titan Ice Dome. Other examples can be found in Reference Forster, Davis, Rand and MooreForster and others (1991), Reference WeertmanWeertman (1993) and Reference Richardson, Aarholt, Hamran, Holmlund and IsakssonRichardson and others (1997), and in Reference Kohler, Moore, Kennett, Engeset and ElvehøyKohler and others (1997) who use GPR to image previous years’ summer surfaces for comparison with mass-balance measurements on a temperate Norwegian glacier.
We have used a 50 MHz antenna frequency in this study. This antenna was chosen because we wanted to map the bedrock and the internal structure of Nordenskjöldbreen at the same time. We did obtain some profiles with 800 MHz antennas, but the penetration was only 10 m. 50 MHz is a relatively low frequency for mapping layers in snow and firn because the pulse waveform (about 3 m in firn) gives a vertical resolution of 1.5 m, which is too long to resolve individual layers. What appears to be a discrete layer is more likely in fact to be a phase front resulting from interference of two or more reflectors. However, the lower frequency permits the deeper penetration required to reach the desired depths at our site.
To map accumulation layers outward from ice-core sites using GPR, the following steps are taken. First, reflecting horizons are dated using ice cores. Depth vs age relations are established in ice cores by such methods as annual-layer counting, radioactive layers or volcanic reference horizons. Depths of known ages can then be converted to radar times. In dry firn this is straightforward and is accomplished via a two-component mixing relation (e.g. Reference RobinRobin, 1975). In temperate or polythermal glaciers, the relation between the electromagnetic wave speed and the medium’s depth is more problematic. The large difference in dielectric constant is ε between ice (ε = 3.15–3.17) and water (ε = 88 at 0°C) means that the radar properties of ice containing water can be especially variable and sensitive to changes in water content of a few per cent. Water pockets in the firn have a strong influence on the absorption and scattering of radar waves, depending on frequency and bandwidth (Reference Smith and EvansSmith and Evans, 1972; Reference DanielsDaniels, 1996).
Next, the dated radar times are associated with the appropriate reflecting horizon, as closely as possible. The reflecting horizons must then be tracked through the GPR profiles. Ideally, horizons are continuous in the profile, but they can become difficult to track in low-accumulation areas where wind scour disrupts or removes layering, in areas with radar noise from crevasses or water and ice inclusions, or in areas in which there has been alteration of the dielectric properties through melting.
The final step is to convert the two-way radar travel times back into depth units, so as to calculate accumulation rates. In this step, the layer depth relative to the total ice thickness must be considered, as layer thinning does not permit the straight division of depth by age to arrive at accumulation.
An additional challenge for investigation is that our GPR data were obtained 2 years after the ice cores were taken. This means that the depths of the reference horizons in the ice cores had to be adjusted to compensate for the additional accumulation.
Equipment and Methods
GPR
In May 1999 a Ramac GPR (Malå Geoscience) with 50 MHz antennas was used to obtain a profile along Nordenskjöldbreen (Fig. 1). The 11.4 km profile follows the central flowline, and passes within 100 m of four ice-core drilling sites, 4, 7, 8 and 10 (Fig. 1), all of which were drilled in May 1997 (Reference IsakssonIsaksson and others, 2001). The equilibrium-line altitude was estimated from core 4 and from GPR stratigraphy to lie at about 660 m a.s.l. (Reference IsakssonIsaksson and others, 2001), at about the location of ice core 4 (Fig. 1).
Stakes set out in 1997 to mark the core sites were not found in 1999, so the GPR profile does not go directly over the cores. The route was determined in the field with handheld global positioning system (GPS), but precise positions along the profile were obtained using differential GPS, with a rover station attached to the snowmobile towing the GPR, and with the base station at Uranusfjellet (Fig. 1). Thus, GPS baselines were never more than 15 km, and the post-processed positions are good to about 1 m. The radar transmitter and receiver antennas were mounted on a non-metallic sledge 2 m apart from each other, which was pulled 7 m behind the snowmobile. The radar control unit and computer were mounted on the snowmobile together with the rover GPS station. Data were collected and stored on a laptop computer. Scans comprise 2048 samples in a time window of 5.11 μs, and were collected at 0.2 m intervals with triggering from a sledge wheel mounted on the snowmobile. Scans were not stacked during data collection. Post-processing of the radar data was done in the program Haescan (Road-scanners Oy). The amplitude zero-level correction was applied, background noise was removed and vertical highpass and low-pass filtering was done.
GPR data were of a good quality, although hard snow surfaces and sastrugi led to some noise in the data. The phase fronts appear fairly constant throughout the profile. Stacking 2–4 pulses helped tracking through areas of noise. Strong reflections were seen down to 25 m at core 7, down to 30 m at core 8 and down to 27 m at core 10. The disappearance of strong layering at these depths is probably due to the firn–ice transition, below which there are significantly weaker density contrasts to produce radar echoes. The firn–ice transition was estimated to occur at 29 m depth at core 10 (Reference IsakssonIsaksson and others, 2001).
Layering was disturbed in several places. These disturbances can be caused by wind scouring, crevasses, ice lenses, refrozen vertical water channels or large ice glands, all of which exist in the percolation zone (Reference PatersonPaterson, 1994), through which most of our profile went.
Ice cores
Core 10 was drilled to 122 m at the summit (1250 m a.s.l.) of Lomonosovfonna in May 1997 (Reference IsakssonIsaksson and others, 2001). The core was obtained using a 105 mm diameter electromechanical drill owned by the Institute of Marine and Atmospheric Research in Utrecht,The Netherlands. The core was drilled almost to bedrock, which is estimated from the radar data to lie at 127 m depth. Temperature measurements from the borehole suggested that the ice was below melting point. The stratigraphy did not seem to have been severely altered by meltwater percolation (Reference IsakssonIsaksson and others, 2001). Core density was measured at the core site. Subsampling work was performed in a cold laboratory at Stockholm University, with further analysis of the core done at the Laboratoire de Glaciologie et Géophysique de l’Environnement, France, the Finnish Forest Research Institute, Rovaniemi, and Tallinn Technical University, Estonia (Reference IsakssonIsaksson and others, 2001).
The uppermost 36 m of core 10 was dated using three different methods: annual-layer counting (Reference IsakssonIsaksson and others, 2001), identification of the 1963 bomb horizon from beta activity (Reference PinglotPinglot and others, 1999), and a one-dimensional ice-flow relation that accounts for layer thinning due to horizontal shear (Reference NyeNye, 1963).
Layer counting using peaks of Na and SO4 2− puts the 1986 layer at 9.21 m core depth (Reference IsakssonIsaksson and others, 2001). Reference PohjolaPohjola and others (2002) use pseudo-annual layers in δ 18O, Cl, NO3 and NH4 to put the 1986 layer at 7.78 m core depth. A third layer-counting method using principal component analysis (J. C. Kohler, unpublished information) yields depths which vary depending on the chemistry components used, but which are intermediate between 7.78 and 9.21 m core depth. We use these two depths as bracketing values for the 1986 layer.
The layer-thinning model predicts the age T at a particular core depth D we (in m w.e.) via
where a c is the mean accumulation rate and H is the ice thickness in m w.e.We use a c = 0.38 m w.e. a−1, a value which correctly puts the depth of the 1963 bomb horizon at the observed core interval 18.5–18.95 m (Reference PinglotPinglot and others, 1999). The 1986 horizon is then predicted to lie at 7.82 m core depth, at the low end of the bracketing values obtained from layer counting.
Shallow cores 7 and 8 were situated at 1044 and 1173 m a.s.l., respectively (Fig. 1), and were drilled using a Polar Ice Coring Office (PICO) drill with a diameter of 75 mm (Fig. 1). Maximum depths were 12.0 and 14.9 m, respectively. Ice thickness at these core sites is unknown, but deeper than the maximum radar penetration depth of 240 m. The 1986 Chernobyl event was seen as a clear 137Cs peak at cores 7 and 8 (Table 1) (Reference PinglotPinglot and others, 1999).
Neither core 7 nor 8 was deep enough to detect the 1963 bomb horizon. At core 10, despite the fact that it penetrated deep enough, there is no evidence of Chernobyl in the 137Cs record. This is attributed to wind scour since the accumulation rate at the summit of Lomonosovfonna is relatively low, and the time-span during which snow contaminated by the Chernobyl event was deposited was relatively short (Reference PinglotPinglot and others, 1999).
We calculate a composite density profile from selected cores rather than assuming variations of the density profile along the GPR profile. This composite density profile is based on cores 7, 8 and 10, and the 1999 core closest to the summit core 10S (Fig. 2).
Density records for 1997 cores 8 and 10 and 1999 core 10S, and the calculated average density.
Core 10 is the best dated of the cores and provides the starting point for layer tracking, with 1963 and 1986 as the target dates, and the results from cores 7 and 8 are used to corroborate the results at the appropriate location along the profile. Because the Chernobyl layer is absent at core 10, the depth of the 1986 layer can only be determined from layer counting. If this dating and the GPR layer tracking are done properly then we should arrive at the appropriate depths at cores 7 and 8, which are observed to contain the Chernobyl peaks. The depth of the 1963 layer at core 10 is better defined than the 1986 layer, but cores 7 and 8 can only be used to ascertain that the 1963 layers are deeper than the lowermost point drilled at each core.
Adjusted dates
Because the core data are from 1997 and the radar data are from 1999, the accumulation for the additional 2 years has to be considered. In May 1999, three new shallow ice cores 10N, 10C and 10S were drilled to depths of 15, 10 and 8 m, respectively, all within 100 m of the core 10 drill site. The density profile and a visual stratigraphy were determined for each core. These data show an average additional accumulation of 2.5 m snow and firn (1.0 m w.e.) for the period 1997–99 (Reference SamuelssonSamuelsson, 2001).
At core 10, the 1986 layer is found from the layer counting to lie at 7.78–9.21 m core depth or 4.4–4.99 m w.e. Considering densification alone, the additional accumulation of 1.0 m w.e. for 1997–99 might simply be added to push the 1986 layer to 5.4–5.99 m w.e. depth or, applying the composite density profile (Fig. 2), 8.9–10.5 m core depth. Since we expect some layer thinning, we use Equation (1) and the long-term average value for ac to determine a resultant change in depth for the 1986 layer of 0.74 m w.e. This results in the 1986 layer lying at 5.14–5.73 m w.e. or 8.6–10.05 m core depth (Table 1). The change in depth over 2 years for the 1963 layer is estimated from Equation (1) to be 0.68 m w.e. The additional accumulation results in a depth of 12.8–13.2 m w.e. for the 1963 layer in 1999, or 19.3–19.7 m core depth (Table 1).
Accumulation for the period 1997–99 was not measured at cores 7 and 8, so we must rely on the core 10 amount adjusted by the 1986–97 ratio of accumulations between the core 7 and 8 sites and core 10. Accordingly, the adjustment is 1.4 m w.e. at core 7 and 2.2 m w.e. at core 8 (Table 1).
Conversion to two-way travel time and accumulation rates
Reference RobinRobin’s (1975) expression for converting core depth D c to two-way travel time t t of a radar wave is given by:
where c = 3 × 108 m s−1 is the speed of light in vacuum,
is the mean core density at the core depth D c, and l a = 2 m is the antenna separation.
We use the composite density profile (Fig. 2) to give D we as a function of core depth D c, and interpolate using a lookup table to calculate two-way travel times for any value of D c (Fig. 3).
Travel time vs depth in m.
After yielding a record of profile distance against t t we then interpolate with a look-up table to determine the water equivalent depth as a function of two-way travel time. Finally, we calculate accumulation ak by simply dividing the calculated target depths Dk by the number of years Tk :
for the periods k = 1963–86, 1963–99 and 1986–99.
Ideally, we would account for layer thinning in calculating accumulation rates, as we do for the dating of the cores. Only in the vicinity of core 10 do we register a bed echo, however, and ice depth is >240 m over most of the GPR profile. Otherwise, there are no bedrock topography maps available, and ice thickness is unknown over most of the GPR profile. In calculating accumulation, therefore, we neglect layer thinning.
Results
Two reflecting layers that gave strong, consistent echoes, and that were within 20 ns (one wavelength) of the appropriate times for the 1986 and 1963 layers, were picked at the core 10 end of the GPR profile (Fig. 4). The two reflecting layers were then tracked and digitized manually down-glacier in the radar profile, passing nearby core sites 8 and 7, to yield a record of profile distance against t t. The GPR-determined accumulation rates 1986–99, 1963–99 and 1963–86 along the profile were calculated (Equation (4)) and are presented in Figure 5.
GPR data from core 10. The 1986 layers (two layers within 20 ns) and the 1963 layer are indicated by white cut lines, and the firn–ice transition depth from GPR data by the black arrow.
(a) Annual accumulation rates (calculated from accumulation in 1986–99 and 1963–99) vs distance. (b) Surface topography along the track indicated in (a), measured using differential GPS.
The distance from core sites 7 and 8 to the GPR profile track is roughly 100 m, so the accumulation from these cores does not permit a direct comparison. However, the spatial variation over 100 m distances is relatively small over much of the profile, and in particular near these core sites, so such a comparison is justifiable. At cores 7 (Fig. 6) and 8, the GPR-derived accumulation rates for the interval 1986–99 (Table 1) compare favorably to those values obtained from the depth to the radioactive layers. The accumulation rates for the interval 1963–99 are not verifiable at cores 7 and 8, but the GPR-tracked layers are deeper than the cores, offering some measure of the consistency (Table 1).
GPR data from core 7. The 1986 and 1963 layers are indicated by white cut lines, and hyperbolas marking water inclusions in the firn are indicated by black arrows.
Errors
Errors in the accumulation rates determined from Equation (4) are given by:
and can be attributed to at least four sources:
-
I error in assigning a date to a GPR target layer (δTk );
-
II error in finding and following a GPR target layer (δTk ), and converting from two-way travel time to depth (δDk = f (δTk ));
-
III error in neglecting layer thinning in the accumulation calculation; and
-
IV error in assuming dry firn, implicit in the simple two-phase mixing concept embodied by Equation (2).
Error I varies according to the dating method used. For the 1986 layer at cores 7 and 8, and the 1963 layer at core 10, error I is taken as half of the length of the core piece containing the Chernobyl peak, 0.13, 0.22 and 0.19 m w.e., respectively, as well as the error for the 2 years additional accumulation, which we assume is 20%, or 0.20 m w.e. Using the Nye time-scale to convert these depth errors to dating errors results in δT = 1.3, 1.6 and 1.6 years or δT/|T| = 12%, 15% and 5% for the 1986 layer at cores 7 and 8, and the 1963 layer at core 10, respectively.
The 1986 layer in core 10 is identified by annual-layer counting using chemistry peaks, and a reasonable error estimate is more elusive. We estimate the error for an individual layer δDi to be half of the average annual-layer thicknesses determined from peak counting in the 1986–97 interval, in this case 0.20 m w.e. The error for a particular layer is then 
, where n is the number of years. For the interval 1986–99, the error would be 0.68 m w.e. As above, we apply the Nye time-scale to obtain the dating error δT = 3.9 years or δT/|T| = 35% for the 1986 layer at core 10. A dating error δT = 3.9 years is fairly conservative since it is greater than the range in values obtained for the mean 1986 depth using different means of layer counting.
Error II for going from core data to flat-lying horizons and for tracking through speckle might be estimated as half of the radar wavelength, or 7.5 ns. However, the error in exact positioning of the cores relative to the radar profile must also be included, since no stakes marking the 1997 core locations remained in 1999 and the placement of the core data within the GPS profile could be in error.
The radar profile was made using geodetic GPS (Javad Positioning System, JPS), which is accurate to about 1 m. The location of core 10 was measured with a hand-held single-phase Garmin GPS receiver in 1997, to an accuracy of about 30 m. The coordinates measured for core 10 in 1999 disagree with the coordinates measured in 1997 after the core drilling by about 20 m, which is within the error limits of the single-phase GPS. The layers at core 10 dip with a slope of about 20 ns over 30 m, so we estimate the positioning error to be about 20 ns. Combined with the half-wavelength error, this results in an error of 21 ns.
Error II in terms of depth can be roughly estimated by combining Equations (2) and (3), neglecting antenna separation (l a = 0 m):
Expanding in terms of errors and assuming that the error in D c is similar to that in D we yields:
Thus, an error in assigning and tracking a layer of 21 ns corresponds to an error in thickness of 1.7 m w.e., or δD we/|D we| = 24%, 16% and 30% for the 1986 layer at cores 7, 8 and 10, respectively, and 13% for the 1963 layer at core 10.
The errors accumulated from Equation (5) are therefore 27%, 22% and 46% for the 1986 layer at cores 7, 8 and 10, respectively, and 14% for the 1963 layer at core 10.
Error III can be estimated by combining expressions for accumulation estimated by Equations (1) and (4), and making an approximation for small values of D/H to yield
Radar measurements show a minimum thickness of 240 m throughout the profile. Given that the measured layers are less than 12 m for the 1986 layer and 24 m for the 1963 layer, neglecting layer thinning in the accumulation calculation leads to a maximum error of 3% for the 1986 layer and 5% for the 1963 layer.
Error IV may produce a bias because water inclusions reduce radar wave velocities significantly (Reference MooreMoore and others, 1999). At lower altitudes, the firn layers are potentially wetter, and there are indeed a few places in the lowest third of the GPR profile in which there are several strong hyperbolic reflectors. The computed depths may therefore be in error at these specific spots. However, a spatially consistent error from wet firn is unlikely from both the evidence of the cores, which did not contain any wet snow or firn, and the GPR profile, which is mostly free of hyperbolas.
Discussion and Conclusion
Temporal variability of accumulation
The annual accumulation-rate values calculated from the GPR data for the time periods 1963–86, 1963–99 and 1986–99 averaged over the entire radar profile were 0.61, 0.58 and 0.67 m w.e., respectively. The 1986–99 accumulation rate A 86–99 shows an increase of roughly 12% compared to A 63–86
However, the error for GPR tracking of the 1986 layer is around 30–35% throughout the profile, while the error for the 1963 layer is roughly 15%. So while our results suggest an increase in accumulation along the GPR profile during the late 1980s and 1990s, the difference between the calculated accumulation rates is less than the errors.
Nevertheless, a precipitation increase in this period has been reported elsewhere.There are instrumental records covering various time periods from Ny-Ålesund, Longyearbyen and Hornsund (Fig. 1) The precipitation record from Ny-Ålesund shows an increase of 14% for the period 1986–98 and a further 6% increase for 1990–98 (Norwegian Meteorological Institute, unpublished information). The precipitation record from Hornsund shows that the precipitation in the 1990s was 50% higher than in the 1980s (Reference MietusMietus, 1999). Since the beginning of the record in 1912 the precipitation in Long-yearbyen has increased, with higher precipitation from the early 1980s to 1996 (Reference Hanssen-Bauer and FørlandHanssen-Bauer and Førland, 1998).
Part, but not all, of these positive precipitation trends is apparently caused by reduced undercatch in the precipitation gauges due to a concurrent warming and thus an increasing amount of precipitation falling as rain (Reference Førland and Hanssen-BauerFørland and Hanssen-Bauer, 2000). Nevertheless, applying a temperature-dependent correction factor, the increase in the annual precipitation at Longyearbyen Airport during the period 1964–97 is still 1.7% per decade (Reference Førland and Hanssen-BauerFørland and Hanssen-Bauer, 2000).
Ny-Ålesund, Longyearbyen and Hornsund are all situated on the west coast (Fig. 1). Lomonosovfonna is located in the central part of northern Svalbard, and has in general a more continental precipitation pattern, with 45–60% less accumulation than the western part (Reference Winther, Bruland, Sand, Killingtveit and MarechalWinther and others, 1998). But while accumulation is less at Lomonosovfonna, it is not unlikely that the temporal trends in precipitation records from the west coast would apply to Lomonosovfonna as well. Our results would suggest that this is the case, that there has been an increase in precipitation over the past decade, although the increase is within the conservatively estimated error limits.
Spatial variability of accumulation
A 86–99 shows large spatial variations over short distances along the measured radar profile (Fig. 5). The first peak in A 86–99 occurs at 3750 m along the profile, with an accumulation rate of 0.90 m w.e. a−1, and is situated in a wet firn area. There are hyperbolic reflections seen in the radar image at this location, indicating the presence of point reflectors in the ice. There is no apparent topographical change, a shallower surface gradient associated with a bedrock overdeepening for example, to explain the sudden increase in accumulation, so it is possible that the inferred depths are too large in this area because of water lenses in the firn.
However, the three other significant peaks in A 86–99 (at 7300, 9200 and 10 050 m along the profile), with values 0.82, 0.91 and 1.0 m w.e. a−1, all occur at higher altitudes, where water in the firn is less likely, and where there are no hyperbolas in the radar record. These accumulation peaks are all located in places where the topography is relatively flat, indicating increased ice thickness within a bedrock depression, and thereby increased accumulation.
Finally, we can compare the trend in increase of precipitation as a function of elevation with values reported elsewhere. The snow water equivalent gradient as a function of elevation is calculated to be 80–150 mm (100 m)−1 at altitudes of 1000–1200 m a.s.l. This is consistent with other observations (Tveit and Killingtveit, 1996; Reference Winther, Bruland, Sand, Killingtveit and MarechalWinther and others, 1998).
In conclusion, we successfully used GPR and ice cores in combination to estimate spatial distribution of snow accumulation between cores, and to estimate temporal changes in accumulation more accurately than at the core site, whose values may be dependent on specific local conditions.
Comparing the annual average accumulation rates for 1986–99 and 1963–99 indicates an increased accumulation during the past decade over the entire length of the profile. This illustrates the importance of combining GPR and core data, as there is essentially little change in accumulation rates for the two periods at core 10, at the summit of Lomonosovfonna. While the increase is within the conservatively estimated error limits, it is consistent with recent changes in precipitation in Svalbard detected in meteorological records at Hornsund, Longyearbyen and Ny-Ålesund.
Acknowledgements
We appreciate the helpful comments of the Scientific Editor, J. B. Johnson, and the two anonymous referees. We are grateful to A. Sinisalo from the Arctic Centre for assistance in the field, and A. Igesund from the Norwegian Polar Institute for construction of the map in Figure 1b. Wihuri Physics Laboratory, University of Turku, Finland, kindly allowed the first author to do her work at their office. Logistical support came from the Norwegian Polar Institute in Longyearbyen. Financial support came from The Finnish Academy, The Finnish Graduate School of Snow and Ice, and The Nordic Council of Ministers.
