Snow accumulation and compaction derived from GPR data near Ross Island, Antarctica

.


Introduction
In recent decades, satellite altimeters have been used to estimate and monitor the Antarctic mass balance by measuring surface height changes around Antarctica (Wingham et al., 1998;Davis and Ferguson, 2004;Nguyen and Herring, Correspondence to: N. C. Kruetzmann (nikolai.kruetzmann@pg.canterbury.ac.nz) 2005).The variability in the mass of polar ice sheets has important implications for sea-level rise and the global radiation balance (Davis et al., 2005).Key uncertainties for determining snow accumulation from changes in surface elevation are snow density and compaction (Arthern and Wingham, 1998) and the spatial variability thereof (Drinkwater et al., 2001).These uncertainties can only be quantified by means of ground truthing.
The amount of snow compaction near the surface is related to mechanical settling during and immediately after accumulation, the overburden pressure by additional snow deposition, and the complex mechanism of temperature metamorphosis ( Van den Broeke, 2008), while melt metamorphosis is mostly restricted to coastal areas.A change in surface height as measured by satellite altimeters, therefore does not necessarily reflect a mass imbalance but may instead be caused by meteorological conditions affecting snow compaction.Additionally, densification processes can change the snow morphology, which affects the height retrieval from reflected radar waveforms such as those from the CryoSat-2 radar altimeter (Wingham et al., 2006).In the present study we describe a method for using ground penetrating radar (GPR) measurements to estimate accumulation and compaction rates at two sites in, and close to, the dry-snow zone in Antarctica to reduce these uncertainties.
Studies which measure compaction of snow are rare.Recently, Arthern et al. (2010) presented an experimental setup for measuring snow compaction down to three distinct depth levels with very high temporal resolution.Zwally and Li (2002) developed a compaction model which shows good agreement with point measurements of compaction rates.Their results show that compaction of dry snow is a continuous process whose seasonal variability depends largely on temperature.In this study we investigate the feasibility of measuring compaction rates over larger areas using a ground based radar system.We believe that our GPR based methodology is complementary to point measurements like those Published by Copernicus Publications on behalf of the European Geosciences Union.
made by Arthern et al. (2010), in that it can measure compaction over larger regions with a higher vertical resolution, albeit on a longer time scale.
Radar has been utilised for glaciological analysis of ice thickness and the detection of internal layers for many years, e.g., Bentley et al. (1979), Bogorodsky et al. (1985), Arcone et al. (1995), Eisen et al. (2003), and Rotschky et al. (2006).Reflections seen in radar recordings can sometimes be associated with distinct accumulation or melt events, or depth hoar layers (Eisen et al., 2004;Arcone et al., 2005;Helm et al., 2007;Dunse et al., 2008).Analysis of internal layers in ice and snow with commercial GPR systems for estimating accumulation has been demonstrated by Arcone et al. (2004), Dunse et al. (2008), and Heilig et al. (2010), amongst others.Here we apply a more rigorous processing methodology based on deconvolution.Our study includes snow density and accumulation information derived from snow pits and firn cores.This is used to reference our processed GPR data and expand the point measurements to larger areas.Additionally, we show that it is possible to acquire estimates of the compaction rates of dry snow by tracking internal horizons in GPR data and comparing layer separations in different years.
The research was carried out at three land ice sites, the locations of which are described in Sect. 2. Section 3 details the deterministic Fourier deconvolution processing scheme used for enhancing the weak contrast of the GPR data.In Sect.4, the accumulation and compaction estimates are presented and results are discussed in Sect. 5. Section 6 summarises our findings.

Data acquisition
Ground penetrating radar data were acquired on the Mc-Murdo Ice Shelf and on Ross Island, Antarctica (Fig. 1a  and b) in November 2008 and 2009.A Sensors and Software pulseEKKO PRO GPR system emitting at a nominal frequency of 500 MHz was used to acquire reflection profiles of the subsurface.According to the manufacturer's specifications, the system has an effective isotropically radiated power (EIRP) just below 10 mW (Sensors and Software Inc., personal communication, 2010).With the settings used for this study the pulse repetition frequency (PRF) lies between 50 kHz and 60 kHz and the system has a duty cycle of 0.02 %.A frequency analysis shows that the system actually operates at an effective centre frequency of about 620 MHz (see Fig. 2b), which is equivalent to an approximate wavelength of 0.35 m in dry snow, assuming a density of 500 kg m −3 .Following Rial et al. (2009), the bandwidth of the system can be estimated from Fig. 2b to be f ≈330 MHz, giving a theoretical resolution v =c/(2• f • √ ε r )≈0.32 m, where ε r is the relative dielectric permittivity of snow (see Sect. 3) at a density of 500 kg m −3 .However, the relative resolution of the system was found to be considerably better.We tested the relative accuracy of the system by recording a GPR profile of a snow pit which had metal stakes inserted into one wall at 0.5 m intervals (not shown).From the apparent separations of the reflection hyperbolas we found that the average error of relative measurements within the snow is 11 %, as long as the distance between the reflectors is greater than the theoretical resolution.
The measurement sites were located within 30 km of New Zealand's Scott Base.Stake farms on land ice were established in 2008 at three locations in different climatic settings and named L1, L2, and L3 (Fig. 1b).The stake farms were set up to measure snow accumulation over a one year time period and to allow repeat GPR measurements along the same profiles.At L1 and L2, 81 stakes were installed on a regular 800 m × 800 m grid at 100 m intervals.The distances between stakes and the regularity of the grid were established with centimetre accuracy using a total station.The layout of a site is illustrated in Fig. 1c.As the study was conducted within a validation experiment for CryoSat-2, the sites were oriented along anticipated satellite ground tracks.The stake in the southwest corner of each farm is labelled A1.The stake in the northeast is I9, with numbers increasing from west to east.For topographic and safety reasons, site L3 was reduced to a 400 m × 400 m grid of 25 stakes.Only odd numbered labels were used in this case to maintain consistent nomenclature for the corners.In the following, The Cryosphere, 5, 391-404, 2011 www.the-cryosphere.net/5/391/2011/directions are relative to grid-north/east, unless stated otherwise.In addition to the stake farms, a profile of 20 stakes was established between L1 and L2, with a separation of approximately 1.5 km between the stakes.The transmit-receive system and the recording equipment were pulled along the grid lines on plastic sleds at a slow walking pace.The GPR data were acquired with a sampling interval of 0.1 ns.Recordings of radar traces (shots) were triggered at regular intervals with an odometer wheel.Using a 135 ns time window in 2008 allowed us to record one trace every 5 cm.In 2009 we attempted to image deeper reflections by using a time window of 205 ns.Due to system limitations, this required an increased horizontal step-size of 7 cm.An additional radar profile was recorded along the line between L1 and L2 using a skidoo, with a recording step-size of 0.4 m and a time window of 180 ns.
Based on previous studies, L1 is located in an area of low accumulation and frequent summer melting (Heine, 1967), on an almost stationary part of the ice shelf.L2, situated in Windless Bight, features considerably higher accumulation rates.Based on Scott Base temperature records we expected occasional summer melting at this site.However, no melt layers were observed in the upper 8 m of snow.GPS measurements corrected with base station data from Scott Base yielded an ice shelf movement of 58 m towards the southwest between November 2008 and October 2009.The third test site, L3, is located on the western slopes of Mt.Erebus (Ross Island), at an altitude of approximately 350 m a.s.l., in the dry snow zone and on undulating terrain.Here, GPS measurements show that this area had moved 3.4 m towards the Erebus Glacier Tongue.In both years, density profiles were taken in at least one snow pit at each site.Densities were measured by weighing known volumes of snow.Additionally, firn cores were drilled and logged in 2009 to obtain snow density profiles up to 8 m deep.L1 did not display coherent layers in radar or density profiles in either year.The irregular and low accumulation, high wind, and frequent summer melting at this site probably prevent the formation of clear stratification in snow pits and GPR images.Consequently, the data from this site are not discussed further.

GPR data processing methodology
In many cases the processing of GPR data has been adapted from the processing of seismic recordings.One frequently used procedure is to calculate the envelope of the received signal via the Hilbert Transform (e.g.Taner et al., 1979), thereby removing the phase information.The resultant trace gives a picture of the instantaneous amplitude of the received signal, but is still strongly influenced by the source signature.Deconvolution, the common remedy to this problem in seismics, has been shown to be more difficult for GPR data (e.g., Turner, 1994;Irving and Knight, 2003).Two key reasons for this difficulty are dispersion of the emitted waveform and its non-minimum-phase character (Belina et al., 2009).The former causes changes in the shape of the radar wave as it travels through the medium, which makes the task of removing one specific waveform inaccurate.The latter relates to the energy distribution of the waveform emitted by most commercial GPR systems, which has its maximum close to the centre of the time domain pulse rather than being frontloaded.This can lead to non-convergent deconvolution operators.Recently, Xia et al. (2004) and Belina et al. (2009) successfully tested deconvolution techniques for GPR recordings on low-dispersion soils.Additionally, Spikes et al. (2004) used a "spiking deconvolution in RADAN" (S.Arcone, personal communication) for simplifying firn radar profiles.In this study we use a similar method, the deterministic Fourier deconvolution, to analyse internal radar reflections of dry Antarctic snow, which is also a low-dispersion material.
A common assumption in analysing GPR data is that the distribution of dielectric contrasts in the ground is random.This assumption is also known as the whiteness hypothesis (Ulrych, 1999), because a (successfully) recovered reflectivity profile is expected to have a spectrum that is similar to that of white noise.While the whiteness hypothesis is widely accepted in a geological context, though sometimes modified to a "blueness hypothesis" (Walden and Hosken, 1985;Ulrych, 1999), it is not immediately evident that it should also apply to the reflectivity structure of stratified snow.Snow deposited in different weather conditions will have variable permittivity based on the thermodynamic properties of the environment at the time, and thus successive layers may be correlated.Nevertheless, the small-scale details of the contrast between these layers are still likely to be random in nature, even if there is some correlation.Hence, the assumption that the spectrum of the output of the deconvolution should be at least whiter than the recorded radargram is likely to be true.
The GPR data recorded for this study can therefore be assumed to measure a medium which consists of well-defined layers with variable dielectric properties.The conductivity of dry snow is very small and the imaginary part of the dielectric permittivity can be neglected (Kovacs et al., 1995).The real part of the relative dielectric permittivity, ε r , of dry snow can be related to its density, ρ (in kg m −3 ), using the empirical formula (Kovacs et al., 1995): The succession of snow layers with different ε r can be thought of as a reflectivity profile, r(t).The emitted signal is partially reflected at each interface between these layers, and the intensity of the reflection depends on the magnitude of the dielectric gradient.Therefore, the reflectivity profile can be directly related to snow density variations (Eisen et al., 2008).The received radar signal, s(t), can be described www.the-cryosphere.net/5/391/2011/The Cryosphere, 5, 391-404, 2011 as the convolution of the emitted waveform, e(t), with r(t) and a noise term, n(t): Fourier deconvolution provides a method to remove the effect of the high frequency carrier wave, e(t), and to recover the reflectivity profile from the radargram.According to the convolution theorem, the Fourier transform (FT) of the convolution of two or more functions is equal to the pointwise multiplication of the Fourier transforms of the individual functions (Bracewell, 2003): Thus, if the emitted waveform can be estimated, and assuming the noise term is negligible, dividing the FT of the recorded trace by the FT of the waveform results in the FT of the reflectivity profile.An inverse Fourier transform (IFT) can then be used to derive the reflectivity profile: If the source waveform is known and does not change with time (i.e.depth in the medium), this is referred to as deterministic deconvolution (Yilmaz, 1987), since it is a well defined mathematical problem that has a single solution.To avoid the problem of non-convergent filter functions mentioned above, we perform the deconvolution in the frequency-domain.
To successfully use Eq. ( 4) for recovering r(t), we measured the emitted waveform (see Fig. 2) by holding the transmitting and receiving antennae above the ground, directly facing each other.To avoid interference between the air-wave and the ground-reflection, the transmitter and receiver were placed 2.3 m above the surface, and 2 m apart.
Figure 2a shows the returned signal as a function of time for 100 superimposed shots.Two distinct returns are observed, the air-wave that travels directly from the transmitter to the receiver, and the ground-reflection.The small variation between shots displayed in Fig. 2a illustrates that the pulse emitted by the system is transmitted at a consistent phase and can be assumed to have the same form from one shot to another.We use an integrated version of the air-wave as our estimate of the emitted waveform for the deconvolution process.The frequency spectrum of the averaged trace, Fig. 2b, shows that the centre frequency of the pulse lies close to 620 MHz.Before performing the division in Eq. ( 4), the spectrum of the emitted wave is whitened by ten percent of the largest spectral peak in order to avoid over-amplification of frequencies of small amplitude (Yilmaz, 1987).
An overview of the processing steps and the order in which they are applied is given in the following list: 1. Dewow -a DC-correction to remove a localised offset caused by receiver saturation (Sensors and Software Inc., 2006).2. Linear gain -compensates for the loss of signal power due to the spherical spreading of the wave as it propagates through the medium.
3. Fourier deconvolution -removes the effects of the carrier wave.
4. Low-pass filter -applied in the frequency domain with a 1550 MHz cut-off based on the spectrum of the emitted wave.
5. Background subtraction -to remove interference due to ringing.This distorts the radargram in the first 6 to 8 ns which are therefore not analysed.
6. Envelope calculation using the Hilbert transform.
7. Integration (horizontal stacking) -used to improve the signal-to-noise ratio.
The Cryosphere, 5, 391-404, 2011 www.the-cryosphere.net/5/391/2011/Step ( 7) also allows us to match the horizontal sampling intervals of the data from the two seasons.After ten-fold stacking of the data from 2008 and seven-fold stacking of data from 2009, the horizontal sampling of both data sets reduces to approximately 0.5 m.Also, as the emitted waveform is not minimum-phase, the deconvolved data is offset by −1.8 ns (upwards).This corresponds to the time from the start of the waveform used for deconvolution to its peak amplitude (Yilmaz, 1987).Accordingly, the timezero for all data has to be adjusted by this amount after step (3).
Figure 3 illustrates the benefits of the Fourier deconvolution for a representative 200 m long radar profile from the L2 site.The simple application of dewow and gain in Fig. 3a shows that the carrier wave is still present in the data.After applying all processing steps except deconvolution, a reduced number of distinct reflections can be clearly identified in Fig. 3b.Including the deconvolution step (Fig. 3c) significantly sharpens these horizons.Comparison of the different panels in Fig. 3 shows that the deconvolution results in a more focussed radargram.
A more stringent test of the quality of the deconvolution processing methodology is to analyse both the autocorrelation and the frequency spectrum of traces before and after processing.Effective deconvolution should whiten the spectrum, reduce the characteristic correlation duration and flatten the tail of the autocorrelation function (ACF) of a radar trace (Ulrych, 1999).Figure 4a shows the ACF of a tracethere is a decrease in autocorrelation time after deconvolution (first minimum is shifted from 0.7 ns to 0.5 ns) and a flattened tail.As expected, the spectrum of the same trace (Fig. 4b) is also clearly whitened after decovolution (Fig. 4c).The ACFs and the spectra in Fig. 4 were computed from a truncated version of the trace that does not include the first 8 ns, for reasons mentioned above.
We use the processing detailed above to identify and track internal reflections in the radargrams.The tracking was performed using the KINGDOM Suite 8.2 software.From a starting point somewhere within the radargram (determined by the user), the program follows the reflection peak from trace to trace.The vertical guide window within which the algorithm searches for the amplitude maximum was set to 2 ns, corresponding to approximately 0.2 m in the vertical.If a reflection is weak or distorted in some part of the profile the selection is corrected manually.
The tracked internal reflection horizons are used for extrapolating accumulation measurements from snow pits and firn cores to larger areas and for deriving compaction rates as a function of depth, d, at the sites L2 and L3.In order to www.the-cryosphere.net/5/391/2011/The Cryosphere, 5, 391-404, 2011 analyse the GPR profiles it is necessary to establish a timedepth relationship.Using the density information from the snow pits and firn cores to estimate ε r with Eq. ( 1), we calculate the velocity, v, of the radar signal in snow: where c is the speed of light in vacuum.Using this analysis, the two-way-travel time (TWT) from a radargram can be converted to depth and vice versa, allowing accumulation estimates to be derived from the GPR data.

Accumulation estimates from stake farms, snow pits, and firn cores
Stake farm and firn core measurements provided independent sources of accumulation data (Table 1).The stake farms at L2 and L3 were installed in November 2008 and revisited in November 2009.Recordings of the stake heights above the snow in both years allow the measurement of one year of accumulation in these areas using the conversion detailed in Takahashi and Kameda (2007).Taking the mean and standard deviation of the snow depths at the 81 stakes from L2 gives an average reduction in stake height of 60.2 ± 8.1 cm or, using the snow densities recorded in a snow pit, 224 ± 21 kg m −2 a −1 .At L3, only 25 stakes were installed.The measurement from one stake was omitted because of a likely error in the height recording.The average decrease in stake height measured at the remaining 24 stakes was 70.1 ± 11.1 cm, equivalent to 304 ± 83 kg m −2 a −1 of accumulation.
A previously identified dust layer (Dunbar et al., 2009) originating from a severe storm (with maximum southerly wind speeds exceeding 55 m s −1 ) that occurred on 16 May 2004 (Xiao et al., 2008) served as a reference point for dating.In 2008, we found the dust layer at a depth of 2.93 m in a snow pit near stake C3 (not shown) at L2.Using the dust layer for dating, we calculate an average annual accumulation of 251 kg m −2 a −1 .In 2009, we were unable to clearly identify the dust layer in a snow pit.However, a firn core was drilled one metre east of the stake G7, to a depth of 8.46 m.Approximately at the depth at which we expected to find the dust layer, the core contained an unusually coarse grained low-density layer (starting at 3.42 m).As we found a similar low-density layer right below the dust layer in the previous year, we believe this could be a depth hoar layer that formed either as surface hoar prior to the storm or underneath the wind crust afterwards.Considering this as a marker for the storm in May 2004 gives an accumulation rate of 245 kg m −2 a −1 .Figure 5a shows the density profile of the firn core and some of the snow pit data from 2009.
At L3, no particularly distinct features could be found in a two metre snow pit in either year (not shown).Slight changes in snow grain size and hardness at about one metre depth are indicative of the previous year's summer layer, but it was not possible to determine this with any certainty.Generally, the snow at L3 was found to be very homogeneous, with a slightly higher average density than at L2.In 2009, we drilled a firn core down to about 7.5 m near stake C3 at L3.The resulting density profile is shown in Fig. 5b.At 5.3 m depth we were able to identify a dust layer similar to that found at L2 in the previous year.Assuming that this is also associated with the May 2004 storm yields an average accumulation of 437 kg m −2 a −1 .This value is considerably higher than the one measured by the stake farm and could either indicate an error in the dating of the dust layer, or a high inter-annual variability in snow accumulation in this area.
The firn core data can be used to determine a depth-density relationship for the two sites, as a reference for converting the vertical scale of the radargrams to depth.Following Alley et al. (1982), we determine an empirical depth(d)-density(ρ) relationship for both sites by fitting an exponential model of the form ρ(d) = ρ i − a • e −c•d to the firn core data, where ρ i = 917 kg m −3 is the density of ice, and a and c are the constants to be fitted.For the upper 0.6 m at L2, where the snow was loose and density difficult to measure by coring, snow pit data were used in order to obtain a fitted depth-density curve.The core at L3 was drilled two weeks after the snow pit and GPR data had been recorded, as weather conditions prevented earlier attempts to return to the site.The core was started at the previous surface level by removing the fresh snow.However, the additional overburden will have caused The Cryosphere, 5, 391-404, 2011 www.the-cryosphere.net/5/391/2011/GPR/dust layer some densification.Therefore, the top 1.9 m of the firn core were replaced by snow pit data when determining the depthdensity relationship.The resultant equations for L2 and L3 are shown in Fig. 5a and b, respectively.As the radar profiles extend below the maximum depth of the firn cores, we also use these equations to extrapolate the density to greater depths in the following analysis.Integrating these empirical relationships allows the estimation of the total snow mass to a certain depth, which can then be used to determine the mean column density, and therefore the TWT, to this depth.

Accumulation estimates from GPR measurements
The snow pit logged at L2 in 2008 is approximately located at the centre of Fig. 3c.The yellow arrow corresponds to the dust layer depth and coincides with a horizon which is more undulating than other horizons in its vicinity.The particularly strong roughness might be related to buried sastrugis caused by the storm event in May 2004 (Steinhoff et al., 2008;Dunbar et al., 2009).Figure 6a and b show processed radargrams from line E4 to E6 at L2, recorded in 2008 and 2009, respectively.The more undulating horizon, as well as several other reflections, are observed throughout the entire survey grid in both years (not shown).The result of tracking the dust layer horizon and eight other distinct reflections in the radar lines in Fig. 6a and b, is shown in Fig. 6c and d, respectively.The tracked reflection horizons are numbered with roman numerals from top to bottom for easier referencing.The vertical scale in Fig. 6d is shifted by 4.7 ns to align the first horizon (I) in both years to facilitate comparison.The yellow line (II) in Fig. 6c is the reflection we associate with the dust horizon.Comparing its path with the neighbouring horizons further illustrates that it is unusually variable; its standard deviation from the linear trend is 0.65 ns, as opposed to 0.30 ns, 0.44 ns, and 0.34 ns for the red (I), purple (III), and black (IV) horizons, respectively.Assuming that the undulating horizon is related to the storm-event, we can calculate the accumulation over the whole grid since May 2004.Tracking this reflection along all 18 grid lines gives an average TWT of 27.5 ± 0.8 ns, which is equivalent to an average accumulation of approximately 269 ± 9 kg m −2 a −1 .The error term is the standard error of the measured depths, reflecting the geophysical variability of the accumulation over the whole site rather than www.the-cryosphere.net/5/391/2011/The Cryosphere, 5, 391-404, 2011 measurement error.The latter is estimated by assuming that the actual depth of the tracked layer is 16 cm (half of the theoretical resolution) above or below the measured value, giving an error of ± 14 kg m −2 a −1 .Figure 7a illustrates the variability in this reflection's depth over the L2 area.While there is no particularly distinct pattern, it appears that accumulation is slightly higher in the north-east part of the grid.Figure 8 displays a radargram of the C-line at L3 from 2009.The location of the firn core is indicated by the red box.Using the dust layer depth (5.3 m) and the densities from the firn core at stake C3, the expected TWT to the dust layer is calculated to be 49 ns (yellow arrow in Fig. 8).Comparing this with the radargram shows that there is a clear reflection horizon at approximately this TWT (tracked in yellow in Fig. 8).Analogous to the approach at L2, we consider this as a marker for the May 2004 dust storm.Accordingly, the average accumulation over this site is calculated to be 404 ± 22 kg m −2 a −1 .In this case the measurement error due to the resolution of the system is approximately The internal horizons in Fig. 8 are noticeably deeper in the middle of the profile, indicating more accumulation at the centre of the grid than at the edges.This inhomogeneity is probably related to local topography since L3 is situated on sloping terrain.The elevation difference between the lowest (I1) and the highest point (A9) is almost 20 m, with A9 located on a local crest and the terrain sloping down towards the north-west.The dip in the observed reflection horizons is on the leeward side of this crest where more drift snow accumulates.Hence, the accumulation rate calculated above needs to be considered as a large scale average for the whole area, rather than an accurate estimate at any particular location.The interpolated accumulation grid for L3 (Fig. 7b), calculated from the tracked dust layer reflection, illustrates the overall pattern.The central dip in the terrain clearly captures more snow than the surrounding areas.
The decreasing trend in the depth of the dust horizon from north to south observed at L2 (Fig. 7a) can be followed south along a GPR transect (Fig. 9) that was recorded as an extension of the line going from I1 to A1.Along this profile, the internal horizons gradually migrate upwards (not shown).After about 14 km, the horizon we associate with the dust layer becomes too indistinct to be reliably identified.At this point, the dust layer reflection is found at a depth of about 1.9 m, which is equivalent to a reduced average accumulation of approximately 165 kg m −2 a −1 (see Fig. 9).This has to be considered a crude estimate, since it assumes the same average snow density between the surface and the dust layer as at L2.A qualitatively similar trend was previously reported by both Heine (1967) and the McMurdo Ice Shelf Project (McCrae, 1984).
Unfortunately, we were unable to reliably date other layers within the snow pit or firn core profiles and associate them with discrete reflections in the radargram.However, due to the high precision of the system, the vertical separation of other apparent horizons can be used to estimate compaction, which is the focus of the next section.

Snow compaction
Most of the apparent internal horizons found in the 2008 GPR data can also be identified in the following year's records.Figure 6c shows nine reflections tracked between stakes E4 and E6 at L2 in 2008, and Fig. 6d shows the same horizons tracked in the 2009 data.One might expect that compaction, caused by temperature metamorphosis and the additional overburden on the surface, will reduce the separation between horizon pairs.This is confirmed by the poor alignment of the bottom horizon (IX) in Fig. 6d compared to Fig. 6c.Clearly, the total separation between the top (I) and the bottom (IX) reflection has been reduced.Assuming that the same horizons have been identified, calculating the average separation between two successive horizons in 2008 and comparing it with that of the same horizon pair in the 2009 The Cryosphere, 5, 391-404, 2011 www.the-cryosphere.net/5/391/2011/data, allows us to estimate the compaction of the snow in the intervening time period.Not all horizons in Fig. 6c and d are suitable for compaction calculations because the high variability of the undulating horizon, for example, does not allow it to be tracked reliably enough to be confident that the same horizon was selected in both years.Similarly, the relatively strong doublehorizon visible between 41 ns and 46 ns in 2008 (Fig. 6a) and between 45 ns and 49 ns in 2009 (Fig. 6b) may be easy to recognise visually, but the noise between both horizons makes them unreliable for automated tracking.In order to establish which of the reflections are likely to be reliably tracked in both years, we calculate the correlation between the TWT profile of each horizon in 2008 and its counterpart in 2009.Additionally, we calculate the distance between pairs of horizons in terms of TWT for each year and the correlation of these relative TWT profiles.Only those horizons that show a correlation greater than 0.5 in all cases are used for compaction calculations.To ensure accuracy, we also require a minimum average TWT difference between horizons of 10 ns, corresponding to approximately 1 m vertical separation.
Fig. 10.Snow compaction along the line from stake E1 to E9 at L2, calculated from the changes in the average separations between some of the horizons shown in Fig. 6c and d.Horizontal error bars are a combination of the relative accuracy of the radar system and spatial variability of the horizons along the line.
For example, the average separation of horizon III and horizon IV in 2008 is 18.7 ns (Fig. 6c), but only 18.4 ns in 2009 (Fig. 6d).This means that the snow between these two horizons was compressed by about 1.6 % during the intervening time period, assuming a constant wave speed.Using the density information from the firn core, this translates into a compaction of about 1.6 cm m −1 .The same calculations for five other horizon pairs result in the compaction diagram shown in Fig. 10.Clearly, the compaction decreases with depth.The horizontal error bars are a combination of the standard errors of the respective horizons' depths along the whole line (spatial variability) and the relative accuracy of the system.The vertical scale in Fig. 10 relates to the depth of the horizons in 2008, i.e. before the compaction has occurred.Since the time between the acquisition of the GPR datasets was nearly one year (355 days), the values in Fig. 10 can be considered estimates of annual compaction rates.
The TWT difference between two horizons is converted to a physical separation by calculating the mean depth between the two horizons for the whole profile in 2009, determining the density at this depth with the formula shown in Fig. 5a and then using this density to estimate the radar velocity between the two horizons via Eqs.( 1) and ( 5).The density of the snow between two horizons is determined from the TWT in 2009 (post-compaction) and assumed to be the same in 2008 (pre-compaction), as the empirical depth-density relationship is purely based on data from 2009.Therefore, we likely underestimate the radar velocity, and hence the distance, between the horizons in 2008, making our compaction rates conservative estimates.For example, assuming a constant density profile with time (Sorge's law; Bader, 1954), the average snow density in a layer initially between 5 to 10 m depth at L2 would increase by 2.5 % between 2008 and 2009.This translates to an underestimation of the wave velocity www.the-cryosphere.net/5/391/2011/The Cryosphere, 5, 391-404, 2011 in 2008 by approximately 1 % and an increase in the compaction rates of about 0.5 cm m −1 .The apparent expansion observed below 10 m in Fig. 10 probably indicates the error in our measurements, and the actual amount of compaction over a one year time period at this depth is too small to be measured with our system.Furthermore, below 8.4 m the conversion from TWT to depth is based solely on an extrapolation of the exponential fit to the density data from the firn core, causing additional uncertainties in the calculations.
The same calculations were performed for five more lines at L2 (A2-I2, A5-I5, A7-I7, B1-B9, and H1-H9) and are summarised in Fig. 11a.The vertical error bars indicate the average separation of the two horizons in 2008.Figure 11a shows that the variability in layer depths and compaction between these six lines is small and thus tracking was not repeated for the remaining lines.The thicker, black data points present the average compaction between each of the horizon pairs and the connecting line illustrates the decreasing trend.
Figure 11b was calculated using the same methodology on seven internal horizons at L3.At this site, the variability of the compaction rates and the average depths is significantly higher.Therefore, the calculations were performed for all grid lines, except for A9-I9 due to corrupt data.The larger spread is in accordance with the observed spatial variability in the accumulation pattern and the resulting dipping of the internal horizons (see Figs. 7b and 8).Nevertheless, the mean compaction rate (black line) shows a clear trend of reduced compaction with depth, similar to that observed at L2.
Assuming constant density profiles with time (Sorge's law) and using the measured average accumulation from the stake farms to determine the initial offsets, we can calculate expected compaction curves.The results are the dashed blue lines in Fig. 11a and b.In both cases the general trend matches that of the compaction measurements, but above about 4 m measured compaction rates are higher than those predicted by the model, and lower at greater depths.

Discussion
Table 1 summarises our results relating to accumulation.At both L2 and L3, the measurements derived from the stake farms showed significantly less accumulation than the combined firn core and GPR measurements.This discrepancy is probably largely due to temporal variability.If the readout of the stake depths had occurred one week later, the measurements at both sites would have been noticeably higher due to the occurrence of a high precipitation event from the 13 to 15 November 2009.The accumulation rates derived from the snow pit and firn core observations should therefore be considered more reliable, since they cover a longer time period.
The average accumulation at L2 -a site that is relatively typical for coastal areas -was found to be 269 ± 9 kg m −2 a −1 , which is much lower than the value reported by Heine (1967)   at station 200, about 6 km north east of our site).There is a strong accumulation gradient in this area and L2 is close to the 320 kg m −2 a −1 contour suggested by McCrae (1984), a value which still lies above our estimate.
Generally, the conversion of TWT to depth based on measured densities is critical for our method and a potential source of error.A bias in the density data would lead to over-or under-estimation of the wave velocity in the snow and therefore the accumulation estimates, but good agreement between core-and snow pit densities from the two seasons at both sites (not shown) allows us to be confident in the density measurements.As the dust layer is a reliable reference point for dating, the difference between our results and previous studies could indicate an overall reduction in annual accumulation or high natural variability in this area.
The Cryosphere, 5, 391-404, 2011 www.the-cryosphere.net/5/391/2011/Moreover, historic accumulation rates derived from stake measurements are very likely to be underestimates, as the consideration of snow compaction -as suggested by Takahashi and Kameda (2007) -is not reported.The gradient in the accumulation map for L2 (Fig. 7a) is found to continue along the transect toward L1 (Fig. 9).The southernmost point up to which we were able to track the dust horizon lies relatively close to the location of a 20 m firn core analysed in Dunbar et al. (2009).They estimate an accumulation rate of 53 ± 20 cm of snow per year.If we assume an average density of 600 kg m −3 for the whole length of their core, this corresponds to 318 ± 120 kg m −2 a −1 .Again, our estimate from the tracked reflection lies significantly lower at 165 kg m −2 a −1 .However, the two points are approximately 3 km apart and the conversion from TWT to depth uses the average column density down to the dust layer determined at L2, which is probably an underestimate for this location.
Only one reflection horizon could be associated with a layer in the snow stratigraphy acquired via simple glaciological tools and visual observation.Such difficulties in linking snow pit and radar observations are quite common (Harper and Bradford, 2003) and are due to the limited resolution of our density profiles (Eisen et al., 2003).Furthermore, the study sites investigated here were located in areas of low dielectric variability.The snow was dry and homogeneous and therefore contained few major dielectric contrasts, such as might be caused by occasional melt layers (e.g., Dunse et al. 2008).Alley (1988) and Arcone et al. (2004) suggest that a combination of thin layers and depth hoar are likely sources for GPR reflections in dry snow.While we did not observe many distinct hoar layers in the snow pit and core data, some may have been overlooked due to the coarseness of the recorded density profile.Accurate identification of the origin of the radar reflections would require high-resolution data from e.g.dielectric profiling, as suggested by Eisen et al. (2004) and Hawley et al. (2008).
The locations of the various horizons at the crossover points of the grid are consistent between perpendicular profiles.In most cases, the same reflection tracked along different profiles can be found within two time samples (2 × 0.1 ns which corresponds to 4 cm) at the point of intersection.The concurrence of the horizon depths is equally high for all reflections at both L2 and L3, showing that the precision of the processed radar data is higher than the theoretical resolution might suggest.This high precision allowed us to estimate snow compaction with depth from changes in the separations of internal horizons from one year to the next.
At both L2 and L3, the average compaction measured between 2 and 13 m depth over a one year time period ranges from 0 cm m −1 to 7 cm m −1 .The sudden drop in compaction observed below about 4 m at both sites (see Fig. 11) could be related to a change in the compaction mechanism.Fresh snow largely compacts via settling, but above a density of 550 kg m −3 sintering is usually considered to become dom-inant (Maeno, 1982;Van den Broeke, 2008).According to the density profiles in Fig. 5, the 550 kg m −3 level lies around 6 m depth at both sites, which compares well with the model estimate by Van den Broeke (2008), who gives a range of 5 to 8 m depth for the 550 kg m −3 level in this area.However, the depth at which we observe a change in compaction rate (4 m) is more shallow than this theoretical threshold density.A recent study by Hörhold et al. (2011) suggests that the "classic" picture of snow compaction is too simple, but better resolved density measurements would be required to explain the origin of the observed step in compaction rate.Similarly, the origin of the dip in compaction rate between 7.5 m and 8.5 m at L2 (Fig. 11a) could be a result of above average snow density at this depth, possibly related to one or two particularly warm summers, but we do not have sufficient data to test whether this is a measurement error or a subsurface feature.
High resolution measurements of densification rates to directly compare with our results are sparse.As we cannot identify the 2008 surface in the radargrams from 2009, it is not possible to calculate the total compaction of the whole snow column.However, total compaction between 5 m and 10 m depth can be estimated from Fig. 11 and compared to the measurements from Arthern et al. (2010).Summing up the average compaction (black lines in Fig. 11a and b) between 5 m and 10 m, gives a total compaction of approximately 4.3 cm at L2 (assuming zero compaction for the deepest 30 cm, where the results are negative) and 5.5 cm at L3, over a time period of 355 days.Using the daily rates given in Table 3 of Arthern et al. (2010) to calculate the total compaction over the same time span and depth range gives 5.7 cm for their "Berkner Island" site, the only site at which their strainmeters worked throughout the whole trial time.While our results qualitatively agree with the "Berkner Island" data, the other sites detailed in Arthern et al. (2010) show considerably higher compaction rates.As mentioned above, one reason for this could be that we do not take into account the change in snow density due to compaction and its effect on the velocity of the radar signal, since this would require additional density data for 2008.Ultimately, the discrepancies between Arthern et al. (2010) and our results could also be due to a difference in climatic conditions, since most of their sites are located in regions with lower mean annual temperature, lower latitude, higher elevation and higher annual accumulation.Therefore, considerable differences in the compaction behaviour of the snow could be expected.

Summary and conclusion
The deterministic Fourier deconvolution scheme suggested here can be used to remove some of the effects of the carrier signal from GPR data of dry firn, provided that it was recorded with a system that has a very stable output.The result is a more focussed radargram with improved contrast.This type of processing improves the identification and precise tracking of weak internal reflection horizons www.the-cryosphere.net/5/391/2011/The Cryosphere, 5, 391-404, 2011 associated with density variations.The processing also facilitates recognition of the same horizons in follow-up surveys.Our methodology also has the potential to work successfully for recordings in areas that are subject to sporadic melt events, although one should keep in mind that an essential assumption of the deconvolution is that there is no -or only very little -frequency dependent absorption and dispersion in the medium.
Using the thusly processed radargrams we extrapolated point measurements of average accumulation from a snow pit at site L2 (251 kg m −2 a −1 ) and a firn core at site L3 (437 kg m −2 a −1 ), to a larger area by identifying a dateable dust layer horizon in the radargram.From the GPR data the extrapolated average accumulation over the 800 m × 800 m site on the ice shelf in Windless Bight (L2) was found to be 269 ± 9 kg m −2 a −1 .The 400 m × 400 m grid on Ross Island (L3) showed higher variability in the internal horizons with an overall average accumulation of 404 ±22 kg m −2 a −1 .Stake farm readings at both sites, maintained over approximately a one year time period, measured an accumulation of 224 ±21 kg m −2 a −1 at L2 and 304 ± 83 kg m −2 a −1 at L3.The discrepancy between these values and the combined firn core and GPR measurements was probably caused by the short time period spanned by the stake observations and the high temporal variability of precipitation events.Additionally, we measured a decreasing accumulation trend along a 14 km long GPR transect heading south from L2.At the southernmost point, the accumulation was about 165 kg m −2 a −1 .
By comparing vertical separations of internal reflection horizons from one year to the next, we were able to estimate compaction rates from GPR measurements down to 13 m depth.This technique might have implications for the validation of the CryoSat-2 satellite altimeter which measures the surface height of ice sheets and shelves, in order to monitor the polar mass balance.
Our results show that internal reflectors found in GPR data, combined with density information, can be used for estimating compaction rates of dry snow.However, estimating densification in percolation areas is probably much more difficult due to the more complex snow morphology (Parry et al., 2007).Frequent repetition of GPR measurements over a longer time period, combined with high resolution dielectric profiling of firn cores, could be used to establish a more detailed representation of time-dependent firn densification for the validation of current firn densification models.The suggested method is applicable over large areas in an efficient and non-invasive manner and is complementary to point measurements of snow compaction at a higher temporal resolution, such as those performed by Arthern et al. (2010) and Heilig et al. (2010).Using a higher frequency system, it might also be possible to improve the vertical resolution of snow compaction data from GPR measurements.

Fig. 1 .
Fig. 1.(a) The Ross Sea region.(b) Measurement area in the western Ross Sea region and corner points of stake farms on the Mc-Murdo Ice Shelf (L1 and L2) and on Ross Island (L3).The circles between L1 and L2 indicate the locations of additional stakes installed for accumulation and radar measurements.(c) Outline of the measurement grid and numbering of the stake farms.(Envisat ASAR image courtesy ESA)

Fig. 2 .
Fig. 2. (a) One hundred superimposed shots recorded with facing antennae.The air-wave and the ground reflection can be clearly distinguished.(b) The amplitude spectrum of the averaged air-wave shows the centre frequency to be around 620 MHz.

Fig. 3 .
Fig. 3. Radar profile at L2 between the stakes C2 and C4 (2008) after applying (a) dewow and gain filters, (b) the full processing excludingand (c) including deconvolution.The yellow arrow in (c) indicates the depth and location of the dust layer identified in the snow pit.The vertical scale of the image is exaggerated, representing 200 m in the horizontal direction and about 13 m vertically.

Fig. 4 .
Fig. 4. (a) ACF of a radar trace at L2 before (dashed line) and after deconvolution (solid line).(b) Amplitude spectrum of the same trace before deconvolution and (c) afterwards.The spectra are normalised to a maximum value of one.

Fig. 5 .
Fig. 5. Snow density profiles in 2009 from snow pits (grey shaded area) and firn cores at sites (a) L2 and (b) L3.In both cases, the snow pit was logged according to visually identified stratigraphy, while the core was largely logged in 10 cm intervals.

Fig. 6 .
Fig. 6.Processed radargrams from stake E4 to E6 at site L2, recorded in (a) 2008 and (b) 2009.Nine tracked reflection horizons (I-IX) are shown for (c) 2008 and (d) 2009.The vertical scale in (d) is shifted up by 4.7 ns, aligning horizon I.

Fig. 7 .
Fig. 7. Accumulation maps for (a) L2 and (b) L3, based on the reflection horizon associated with the dust layer tracked along each line.A four metre interpolation grid resolution is used in both cases.

Fig. 8 .
Fig. 8. Processed radargram from C1 to C9 at site L3 (2009).The red box marks the approximate location of the firn core; the yellow arrow indicates the expected TWT to the dust layer.

Fig. 9 .
Fig. 9. Accumulation along a transect from L2 to L1, estimated by tracking the dust layer reflection for 14 km.
at the closest station (510 kg m −2 a −1

Fig. 11 .
Fig. 11.Snow compaction vs. depth for (a) L2 and (b) L3.(a) shows the compaction calculated from six lines at L2, while (b) includes all lines from L3 except one.The vertical axis represents depth before compaction has occurred.The horizontal error bars illustrate only the spatial variability along each of the lines, not measurement error.The vertical error bars show the thickness of the layers before compaction.The thick black line represents the average compaction for each horizon pair and the dashed blue line represents the expected compaction when assuming a constant density profile with time (Sorge's law).