Edinburgh Explorer Spatially extensive estimates of annual accumulation in the dry snow zone of the Greenland Ice Sheet determined from radar altimetry

. We present estimates of accumulation rate along a 200 km transect ranging in elevation from 2750 to 3150 m in the dry snow zone on the western slope of the Greenland Ice Sheet. An airborne radar altimeter is used to estimate the thickness of annual internal layers and, in conjunction with ground based snow/ﬁrn density proﬁles, annual accumulation rates between 1998 and 2003 are derived. A clear gradient in the thickness of each layer observed by the radar altimeter and in the associated estimates of annual accumulation is seen along the transect, with a 33.6% ± 16% mean decrease in accumulation from west to east. The observed inter-annual variability is high, with the annual mean accumulation rate estimated at 0.359 m.w.e. yr − 1 (s.d. ± 0.049 m.w.e. yr − 1 ) . Mean accumulation rates modelled using meteorological models overestimate our results by 16% on average, but by 32% and 42% in the years 2001 and 2002. The methodology presented here demonstrates the potential to obtain accurate and spatially extensive accumulation rates from radar altimeters in regions of ice sheets where ﬁeld observations are sparse, and accumulation rates greater than several tens of cm.


Introduction
It is known that the arctic region is warming and that the margins of the Greenland Ice Sheet are experiencing substantial thinning (Hassol, 2004;Steffen et al., 2004;Zwally et al., 2005;Pritchard et al., 2009). Increasing fresh water input from the Greenland Ice Sheet (GrIS) to the oceans has created a growing concern about the effects it could have on climate and sea level rise. Although changes in the volume of an ice sheet may be estimated by observing variations in its surface elevation (e.g. Wingham et al., 1998;Shepherd et al., 2001 andZwally et al., 2005) and from gravity measurements (e.g. Velicogna, 2009), these methods do not distinguish ice dynamics and accumulation fluctuations which may introduce errors in mass balance estimates.
Observations show that surface melt extent and duration have increased dramatically since the late 1970s (Abdalati and Steffen, 1997;Steffen et al., 2004;Mote, 2007;Hanna et al., 2008). Meanwhile, there are reports of an increase in the surface elevation at higher altitudes since 1990 (Krabill et al., 2000;Nghiem et al., 2005;Hanna et al., 2006), supporting calibrated climate models which show that some regions of the ice sheet have experienced higher than average accumulation rates, especially in the south (Burgess et al., 2010). A recent modelling study found larger than average accumulation rates across the GrIS between 1958 and2007, and that the surface mass balance over the same period is between 32-63% higher than previous estimates (Ettema et al, 2009). However, Greenland's weather and accumulation patterns have high spatial and temporal variability (Mosley-Thompson et al., 2001;Steffen et al., 2004;Wang et al., 2007), and elevation changes may be caused by factors other than accumulation and melt such as snow compaction and densification processes and by the redistribution of drifting snow . For this reason, accurate estimates of mass balance made from elevation change observations require continuous measurements of accumulation and compaction rates across different regions of an ice sheet. Previous estimates of annual accumulation rate have been made using data collected by scatterometers (Drinkwater et al., 2001), and from Rayleigh-scattering microwave model inversion techniques applied to radar and passive microwave sensors data (Munk et al., 2003;Flach et al., 2005). However, direct observations of snow accumulation in Greenland are limited, and measurements are affected by wind patterns and for the most part made at stations located in coastal regions (Dethloff et al., 2002), while accumulation estimates derived from ice cores and snow density measurements are also limited spatially.
Here, we use airborne radar altimetry data and snow density measurements derived from snow/firn cores and neutronprobe observations obtained over the dry snow zone of the GrIS to estimate annual accumulation rates for the years 1998-2003 along a 200 km transect ranging in elevation from 2750 m to 3150 m on the western slope of the GrIS (Fig. 1). The estimates obtained are compared with snow accumulation rates derived from a weather prediction numerical model run retrospectively using data from the European Centre for Medium-Range Weather Forecast (ECMWF) 45-year (1957-2002) reanalysis ERA-40 and ECMWF operational analysis from 2002 (Hanna et al., 2005(Hanna et al., , 2006(Hanna et al., and 2008. This model has been recently corrected and calibrated for Greenland with recent field measurements . Technical advances in radar altimetry have made the European Space Agency's Airborne SAR/Interferometric Radar Altimeter System (ASIRAS) capable of resolving internal layers in the snowpack and firn (Hawley et al., 2006;Helm et al., 2007). We extend the previous ASIRAS analyses in order to resolve spatial and temporal trends in the estimated accumulation rates presented. ASIRAS is designed to emulate the operation of the SAR/Interferometric Radar Altimeter (SIRAL) onboard the recently launched CryoSat-2 , thus one of the objectives of this paper is to demonstrate the viability of estimating accumulation rates from CryoSat-2 or other space borne altimeters with similar operating characteristics.

Data and Methodology
The data reported here were derived from observations made by ASIRAS during the 2004 CryoSat Validation and Calibration field experiment. On 17 September 2004, the ASIRAS radar altimeter was flown along the Expedition Glaciologique au Groenland (EGIG) line (Fischer et al., 1995) over the GrIS operating in High Altitude Mode (HAM). We present a 200 km transect, along the EGIG line, on the western slope of the dry snow region of the ice sheet ranging in elevation from 2750 m to 3150 m. The system is a pulse-width limited, coherent, and phase sensitive radar altimeter that operates in the Ku band, at a central frequency of 13.5 GHz and with a bandwidth of 1 GHz. The system has a footprint of 4.5 m along track (after SAR processing) and the large bandwidth provides high vertical resolution, so it is possible to analyze in detail the structure of the observed snowpack.
Radar signals can penetrate snow, and are partially backscattered by changes encountered in both snow density and crystal structure. Thus stratigraphically distinct layers in the snowpack, as are commonly formed early in the fall in the dry snow region of the GrIS (Autumn hoar), make it possible to identify widespread internal layers reflecting annual accumulation (Hawley et al., 2006). An advantage in using the Autumn hoar to identify annual layers is that it coincides with the commonly used concept of a mass balance year for an ice sheet, lasting typically from the onset of accumulation (in Autumn) in one year to the cessation of ablation Table 1. Snow and firn density values (g cm −3 ) for each annual accumulation period and for each altitude where accumulation rate is estimated. Observed annually averaged densities are presented from field data collected at two sites located above (Summit) and below (T21) our transect.
Year/Height (m) 2650 T21 2750  2810  2870  2930  2990  3050  3090  3110  3150  (end-of-summer) the next year (Anklin and Stauffer, 1994). Thus, when referring to the observed accumulation rate in 1999, we are referring to the period from the end-of-summer of 1998 to the end-of-summer of 1999.
To estimate the distance travelled, we calculate the speed of light in snow based on the refractive index of snow. By neglecting any water content (dry snow), we assume we are working on a lossless media so we can obtain the refractive index of snow n e directly from the snow density, for which we use mean values for each of the layers detected based on snow density profiles obtained from snow/firn cores and neutron-probe observations. The refractive index for dry snow with a density under 500 kg m −3 is given by (Yankielun et al., 2004): where ρ s is the density of snow in kg m −3 . Using equation 1, we can estimate the speed of light v c through the snowpack from its linear relationship with the refractive index of snow, v c = c/n e , where c is the speed of light in vacuum. Two snow density profiles obtained from shallow firn cores and neutron probe measurements at the lowermost and the highest points of the survey are used. Our lower density profile from T21 (2650 m) was obtained from neutron probe measurements made in May 2004 (Hawley et al, 2006). Our upper density profile, from Summit (3200 m) was obtained by averaging the density estimates derived from both a shallow firn core and from neutron probe observations made in June 2004 (Hawley et al., 2008). Summit is located 150 km north from the surveyed path, but along the ice divide and only 50 m higher than our easternmost ASIRAS observation. For each layer (year), at both the upper and lower sites, a mean annual snow density was estimated. The annual layer density values at each of our nine sites were then determined by linear interpolation between the upper and lower observed values according to elevation, and the estimated densities at each point along the transect were used to derive refractive index and accumulation (Table 1). Although this is a simplistic approximation, climatic conditions in the dry snow regions of the GrIS above 2500 m are similar and there is no reason to suspect that changes in density in the area of study are not directly dependant on elevation.
We estimated the depth of each observed layer at 25 km intervals along the 200 km transect. At each point, 100 individual waveforms are averaged to minimize small scale variations in layer thickness. There are clear differences in the degree of variability in waveform structure with elevation. Between 2800 m and 2950 m elevation (50 to 110 km along the transect), the depth of the observed layers vary as much as ± 9 cm between the different waveforms, while at higher elevations no differences larger than the system range accuracy, which is estimated to fall within 3 cm (Hawley et al., 2006), are observed. In a previous study, winter accumulation rate estimated from ASIRAS at one site over the percolation zone was estimated at 0.63 ± 0.06 m.w.e., which compared favourably with the 0.61 m.w.e. measured in the field (Helm et al., 2007). Uncertainty is introduced from the assumption that each layer delimits exactly one annual accumulation period. In reality, the snow density changes that cause radar reflections are not formed at exactly the same time every year. For example, it is known that the 2002 summer melt season was unusually long with areas experiencing 2 weeks more melt than on average (Steffen et al., 2004), so if temperatures remained unusually high across the whole ice sheet, the formation of the autumn hoar may have been delayed, slightly modifying the length of the period for which the accumulation rate estimates are made. However, this is unlikely due to the predominantly low temperatures in the dry snow zone of the GrIS and the "annual" layering will still reflect an accumulation "year".

Results and Discussion
The radar altimeter observations along the 200 km transect are shown in Fig. 2. Individual waveforms are aligned with respect to the signal in each waveform estimated to be the surface, which reveals horizontal continuity in the firn structure with clear layers observable to a depth of about 10 m (travel time = 80 ns). A strong return from the surface is evident, with several weaker reflections from the volume of  the snowpack along the whole transect. These weaker reflections occur at interfaces within the snowpack where sharp changes in the firn density/structure occur, as described earlier, and represent annual accumulation periods. The distance between the annual layers decreases with increasing elevation; this is expected, since the mean snow accumulation rate increases in a westerly direction away from the ice divide along the EGIG line (Thomas, 2004). Seven distinct layers are clearly visible along the whole transect, with several more appearing at depth at higher elevations in the east where snow densities are lower. It has been observed elsewhere that radar signals penetrate deeper into the snowpack as snow density decreases, thereby revealing additional layers with increasing elevation and decreasing snow density (Scott et al., 2006 Figure 3a shows the estimated thickness of the annual layers for 1998-2003. The decrease in layer thickness with increasing elevation is evident in all years, and there is substantial inter-annual variability observed in the rate at which thickness decreases with elevation. A decrease in thickness of 18% and 14% is observed in 2000 and 2003; 27% and 23% in 1998 and 2002; and 45% and 50% in 1999 and 2001 respectively. Figure 3b combines the estimated annual layer thickness (Fig. 3a) with derived snow density to estimate the total accumulation for each annual layer down to 1998. There are large inter-annual variations in accumulation ranging from 0.267 ± 0.054 m.w.e. yr −1 in 1999 to 0.447 ± 0.037 m.w.e. yr −1 in 2000 (Table 2). Accumulation rate decreases from west to east along the 200 km transect Table 3. Mean accumulation rates (in m.w.e.) with elevation derived from ASIRAS every 25 km along the transect in conjunction with field data (Fischer et al., 1995;Anklin and Stauffer, 1994) and model results. The observed % change in accumulation are estimated from a comparison between the ASIRAS data and ground observations made by Anklin and Stauffer. (1994).
Elevation 1998-2003Fischer et al. Anklin and Stauffer 1998-2003 ASIRAS derived (1984-1989) (1977-1989 Figure 3a show the thickness for each layer. Figure 3b show the annual snow accumulation estimates obtained from the different layers identified (m.w.e.). The estimates were made every 25 km. reducing by 14% and 18% for 2000 and 2003; 36% and 30% for 1998 and 2002; and 50% and 55% for 1999 and 2001 (percentages calculated from the trend shown by a linear regression estimated for each year, with an average correlation coefficient r 2 of 0.853). These gradients in accumulation rate confirm that the dominant trajectory of moist air masses follows an easterly path along the western slope of the GrIS in the vicinity of the EGIG line. Furthermore, strong katabatic winds may further contribute to the transfer of snow from higher to lower elevations away from the ice divide. It is noteworthy that while the thickness in the annual layers between the lowest and highest elevations along the transect decreased on average by 25%, the net drop in mean accumulation is estimated to be around 33%, clearly indicating that any loss (or gain) in the ice sheet mass may not be completely reflected by proportional changes in surface elevation .
The calculated mean accumulation rates at 25 km intervals along the transect are presented (Table 3) in conjunction with both the mean annual estimates derived from the ERA-40 data-based model output and estimates made from two previous field based studies (Fischer et al., 1995;Anklin and Stauffer, 1994). Our mean accumulation estimates are higher than between 1977-89 (Anklin and Stauffer, 1994) at sites above 3100 m, with increases in accumulation of 15% and 22%, but with no significant change at sites below 3000 m ( Table 3). The ERA-40 model over-predicts the ASIRAS derived accumulation rates by an average of 16.2%. As stated earlier, there is a clear spatial gradient in accumulation rate, with an annual mean of 0.435 m.w.e. ± 0.065 at 2750 m elevation decreasing to 0.291 m.w.e. ± 0.095 at 3150 m. To emphasize that the assumption of a linear decreasing density from 2650 m to Summit does not introduce large errors, we made a second estimate of the accumulation rates using a constant density for each annual layer derived from the observations at T21 (Table 1) be greatest, the estimate for mean annual accumulation is 0.0298 m.w.e. ± 0.014 higher than observed, equivalent to just a ∼10% increase.
To further illustrate the importance of monitoring short scale spatial and temporal variability, we compare our results with accumulation rates derived from an ERA-40 data-based meteorological model (Hanna et al., 2005). The model runs with meteorological data and thus does not depend on density measurements. The model was downscaled to 5×5 km grid and calibrated against kriged ice-core data mapped to the same grid . Figure 4 illustrates the spatial and the temporal variability of the estimates derived both from the ASIRAS observations (Fig. 4a), and the model run retrospectively (Fig. 4b). The different approaches reveal very similar patterns with decreasing snow accumulation from west to east and similar inter-annual variability, with the exception of the year 2002, in which the model predictions overestimate our results by 42%. Modelled accumulation for a given year may be subject to unusual meteorological conditions that affect the model output. It is known that the GrIS experienced unusually high temperatures in 2002 (Steffen et al., 2004;Wang et al., 2007), which increased precipitation in the coastal regions but may not have significantly enhanced accumulation rates in the higher regions of the ice sheet.

Summary and conclusions
We present estimates of annual snow accumulation rates derived from observations made by the ASIRAS airborne radar altimeter along a 200 km transect of the dry snow region on the western slope of the Greenland Ice Sheet. Different layers are identified within the snowpack, each associated with the annual snow accumulation over six mass balance years for the period 1998-2003. The average annual accumulation is estimated at 0.359 m.w.e. ± 0.049, but the rates exhibit both high spatial and temporal variability, with the mean accumulation rate estimated at 0.267 m.w.e. in 1999 and at 0.447 m.w.e. in 2000. A clear gradient in the thickness of each layer and in the associated estimates of annual accumulation is seen along the transect, with a 33.6% ± 16% mean decrease in accumulation from west to east each year. This trend has been observed in previous studies (Fischer et al., 1995;Anklin and Stauffer, 1994;Burgess et al., 2010) and reflects the dominant transport of water vapour in the region from west to east. Inter-annual variability in accumulation has also been observed previously (Fischer et al., 1995). Comparison with historical records (Anklin and Stauffer, 1994) suggests that the accumulation patterns have increased by 15-20% above 3000 m over the last 20-25 years.
A new generation of radar altimeters have made it possible to obtain accurate estimates of snow accumulation from altimetry by combining the radar observations with field measurements of snow/firn density. Comparison with ac- cumulation rates derived from a meteorological model are favourable although our estimates reveal small scale variability that cannot be predicted by the model. The differences seen between the model and our estimates provide insight on uncertainties in surface mass balance estimates. Assuming similar accumulation patterns, the difference would mean an uncertainty in the mean annual mass gain due to accumulation of approximately 44.5 km 3 (275 km 3 for ASIRAS and 320 km 3 for the model) if applied to a dry snow area of 7.7×10 5 km 2 (approximately 55% of the GrIS), significant if compared to a 1958-2007 estimate of total surface mass balance for the whole GrIS of of 422 ± 37 km 3 yr −1 (Ettema et al., 2009). In the context of the ESA CryoSat-2 mission, designed to obtain very accurate measurements of ice sheet elevation change, the methodology we present should prove particularly useful for improving the accuracy of mass balance estimates of the GrIS. Preliminary radar returns from the mission indicate that CryoSat-2 will have the potential to identify subsurface accumulation layers in the dry snow zone. Due to CryoSat-2's larger vertical resolution compared with ASIRAS (47 cm in free space ); approximately 32 cm in snow with a density of ∼0.4 g cm 3 (Scott et al, 2006)), it may be difficult to identify annual accumulation layers in areas with very low accumulation rates, e.g. in certain regions of Antarctica. However, in dry snow zone areas with accumulation greater than several tens of cm such as in the GrIS, CryoSat-2 should be able to resolve these annual layers. The methodology presented here, using a radar altimeter with similar characteristics to CryoSat-2, therefore shows the potential that these satellites have to yield reliable estimates of snow accumulation over large regions of the ice sheets which are otherwise difficult to monitor.