1. Introduction
Over the past few decades much effort has been put into recovering possible changes in the mass budget and dynamics of polar ice sheets in response to climate change. Mass-balance estimates, which still suffer from large uncertainties, require knowledge of the amount and spatial distribution of snow accumulation. Traditionally, shallow firn and ice cores, supplemented by snow pits, are used to derive the accumulation and climate history throughout the past few centuries. Resulting from variability in precipitation and wind-driven redistribution of snow, surface accumulation rate is highly variable spatially and temporally. Local values depend on the overall climate, regional weather conditions and terrain characteristics at macro- to microscale. For example, in accumulation averages for 4–7 years, Reference Reijmer and van den BroekeReijmer and Van den Broeke (2003) found standard deviations to range from 20% to 70% of the average annual accumulation, based on data from several automatic weather stations in Dronning Maud Land (DML), Antarctica. On an annual scale, large variability in accumulation has been observed from stake measurements by Reference Melvold, Hagen, Pinglot and GundestrupMelvold and others (1998), reaching up to 500 kg m–2 a–1 over a distance of <3 km on Jutulstraumen, DML. In this respect ice cores serve as a useful tool because they provide long-term averaged accumulation rates which more likely represent a typical value for their surrounding area. However, due to the time-consuming and often costly procedure of firn-core drilling, their numbers remain limited.
Sampling density of accumulation records varies strongly over the Antarctic ice sheet (Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others, 1999). Drilling sites are usually concentrated along traverse routes and easily accessible coastal regions. The sparse and uneven data distribution sets high demands on interpolation techniques used to map the broad-scale accumulation distribution over wider regions. Furthermore it is questionable whether firn-core records covering different accumulation periods can be combined in a single map if trends in accumulation history cannot be precluded. Nevertheless, for most regions, firn cores still represent the best in situ source of information that can be used for accumulation compilations.
DML has been under investigation by several countries over the past 50 years, resulting in a large body of available data. This region is particularly interesting because of its sensitivity to conditions over the Atlantic sector of the Southern Ocean, possibly connecting it with the climate record of the Greenland ice sheet. Field activities increased in the 1990s in two regards: (i) within the framework of the International Trans-Antarctic Scientific Expedition (ITASE), established in order to gain a greater understanding of Southern Hemisphere climate within the past 200 years;and (ii) as a result of pre-site surveys for the European Project for Ice Coring in Antarctica (EPICA), aimed at locating a suitable deep ice-core drilling site in the Atlantic sector of East Antarctica.
As more data on accumulation rates have been compiled successively each season, it becomes necessary to summarize and combine the information for a region-wide picture and estimation of snow accumulation in DML. To date, there is no accumulation map specifically for western DML that is based entirely on in situ data. In this paper, we present a complete set of reliable accumulation records available at the time of this study for the western DML region. For spatial interpolation we apply algorithms based on geostatistical methods and exploratory data analysis. In doing so, we aim to recover the large-scale accumulation pattern resulting from overall climate, weather and topographic conditions. Clearly, small-scale variations in accumulation cannot be resolved using point data. We discuss that dataset in terms of its spatial and temporal validity, justifying the combination of unevenly distributed field data representing a variety of accumulation periods at each location.
2. Geographical Description of Western Dml
2.1. Morphology
The area of investigation stretches inland from the Antarctic coast at 70° S to about 80° S, bounded by the 20° W and 15° E meridians (Fig. 1). Based on its topography, western DML can be subdivided into four regions, three of which are: (1) the major ice shelves Riiser-Larsenisen, Ekströmisen and Fimbulisen; (2) the lower inland ice regions Ritscherflya and Hellehallet;and (3) the high-altitude plateau area with an elevation >2500 m a.s.l., i.e. the Wegener Inlandeis made up by Amundsenisen and Wegenerisen. A clear borderline between coastal/lower inland ice regions and the plateau consists of (4) the nunatak areas of the Heimefrontfjella, Kirwanveggen and Mühlig-Hofmann-Gebirge mountain ranges, which penetrate the ice sheet, reaching heights of about 2700m a.s.l. Acting as a barrier to air masses approaching from the Weddell Sea and the Atlantic Ocean to the north, this chain of mountains has a large impact on the accumulation distribution, separating the area into two distinct accumulation regimes.
Fig. 1. Area of investigation and sampling network. Varying symbols for point records were chosen in order to distinguish between different institutions involved in data acquisition. Elevation contours in 200 m intervals are shown as solid black lines.
The surface topography of the plateau is smooth, with wind-shaped features, i.e. small-size sastrugi (10–15 cm high) (Reference StenbergStenberg and others, 1998), indicating that the mean wind velocities in this part of Antarctica are moderate. Undulations with amplitudes on the metre scale and wavelengths on the kilometre scale have been linked to bedrock morphology (e.g. Reference Budd and CarterBudd and Carter, 1971; Reference Rotschky, Eisen, Wilhelms, Nixdorf and OerterRotschky and others, 2004). Large areas of Amundsenisen have been delineated as ice drainage basins feeding the Filchner Ice Shelf with the ice divide as far north as 75° S (Reference Giovinetto and BentleyGiovinetto and Bentley, 1985). Major outlets from Amundsenisen are Stancomb-Wills Ice Stream, draining into the Brunt Ice Shelf, and Jutulstraumen, feeding Fimbulisen. With an area of approximately 124 000 km2, Jutulstraumen is the largest ice stream within our area of investigation. Topography varies to a greater degree coastward of the nunataks, where drainage basins alternate with mountain ridges. Varying slope gradients and aspects as well as a variable surface-wind field introduce complex patterns in the accumulation distribution (Reference Van den BroekeVan den Broeke and others, 1999).
2.2. General characteristics of snow accumulation
DML is situated in the Atlantic sector of Antarctica, and in the sphere of influence of cyclonic systems that move along the coast;the strongest winds blow from the east to northeast (Reference LundeLunde, 1961). Field observations reveal a typical continental precipitation distribution with a general decrease in accumulation with increasing elevation and distance from the open ocean, and with decreasing mean annual air temperature (Reference Giovinetto, Waters and BentleyGiovinetto and others, 1990; Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others, 1999). In DML, the nunataks interrupt this steady modification by forcing the humid air masses to rise, resulting in accumulation extremes of up to 780 kgm–2 a–1 on Jutulstraumen, as reported by Reference Melvold, Hagen, Pinglot and GundestrupMelvold and others (1998) based on stake measurements over a 1 year period. Large regional variations in accumulation occur depending on predominant cyclone pathways, topographic disturbances and katabatic wind activity (Reference Richardson and HolmlundRichardson and Holmlund, 1999; Reference Van den BroekeVan den Broeke and others, 1999; Reference King, Anderson, Vaughan, Mann, Mobbs and VosperKing and others, 2004). A detailed discussion of precipitation and related snow accumulation characteristics in DML was presented by Reference Noone, Turner and MulvaneyNoone and others (1999).
3. Sampling Framework
Our dataset includes accumulation data from snow-pit studies as well as from firn- and ice-core drillings carried out by six nations during 20 different Antarctic field campaigns, held over a period of about four decades. Compiled records and their sources are listed in the Appendix, tabulated by site name, campaign name and year of observation. Because of the susceptibility to errors we did not consider records based on snow-pit stratigraphy alone. Results from this method have been criticized as being generally unreliable due to subjective interpretation. Total obliteration of the annual layering, resulting from strong metamorphism or wind erosion, is considered responsible for the usual overestimations of accumulation rates derived from firn stratigraphy, particularly where accumulation is low (e.g. Reference Picciotto, Crozaz, De Breuck and CraryPicciotto and others, 1971; Reference Giovinetto and ZwallyGiovinetto and Zwally, 2000; Reference FrezzottiFrezzotti and others, 2004).
Sporadic visits to DML began with the Norwegian– British–Swedish Antarctic Expedition in the years 1949–52; however, systematic accumulation data acquisition did not start until the late 1960s. The earliest reliable dataset originates from the third leg of the US South Pole–Queen Maud Land traverse (SPQMLT 1964–68) covering the southernmost part of the area of interest between 75° S and 80° S and between 10° W and 40° E on Amundsenisen (Reference Picciotto, Crozaz, De Breuck and CraryPicciotto and others, 1971). Here the dating of 17 snow pits was based on total β activity deposited as a result of nuclear weapon tests. Data acquisition continued with German activities resulting in nine firn cores drilled within two seasons (1981/82 and 1986/87) at locations near the Georg von Neumayer station on Ekströmisen (Reference Schlosser, Oerter and GrafSchlosser and others, 1999) and on Ritscherflya (Reference Oerter, Graf, Wilhelms, Minikin and MillerOerter and others, 1999). Glaciological work conducted in the framework of the Swedish Antarctic Research Program (SWEDARP) in the seasons 1988/89 (Reference Isaksson and KarlénIsaksson and Karlén, 1994) and 1991/92 (Reference Isaksson, Karlén, Gundestrup, Mayewski, Whitlow and TwicklerIsaksson and others, 1996; this paper) included 11 shallow firn cores along a traverse route starting at Riiser-Larsenisen across Ritscherflya towards Amundsenisen east of the Heimefrontfjella mountains. The latter area was visited again during the Swedish ITASE campaigns in 1993/94, providing accumulation data based on snow pits dated by multi-parameter chemical records (Reference StenbergStenberg and others, 1998). Collaborative Swedish–Norwegian–Dutch activities continued in 1997/98 within EPICA, including snow-pit studies (Reference Stenberg, Hansson and HolmlundStenberg and others, 1999) as well as medium-depth ice-core drillings at Camp Victoria (Reference KarlöfKarlöf and others, 2000) and a depot of the British Antarctic Survey (BAS) (Reference HofstedeHofstede and others, 2004). Within five consecutive seasons from 1995/96 to 1999/2000, EPICA traverse work under the leadership of the German Alfred Wegener Institute (AWI) retrieved more than 30 shallow firn cores and three medium-depth ice cores, with the majority taken on Amundsenisen, but also close to Neumayer station, in the Heimefrontfjella, and on Ritscherflya (Reference Oerter, Graf, Wilhelms, Minikin and MillerOerter and others, 1999, Reference Oerter2000; this paper).
In an extensive mass-balance and glacier research project, the Fimbulisen and Jutulstraumen region was traversed by several Norwegian Antarctic Research Expeditions (NARE) in the seasons 1992/93 (Reference Melvold, Hagen, Pinglot and GundestrupMelvold and others, 1998; Reference MelvoldMelvold, 1999) and 1993/94. They recovered 12 shallow firn cores. In 1996/97 a Norwegian–Swedish–Dutch ground traverse as part of EPICA drilled 15 firn cores along a traverse from Fimbulisen passing Hellehallet and Mühlig-Hofmann- Gebirge and on the plateau along the ice divide leading to the Japanese deep-drilling site at Valkyrjedomen (Dome F). Those records define the eastern boundary of our investigation area (Reference Isaksson, van den Broeke, Winther, Karl¨f, Pinglot and GundestrupIsaksson and others, 1999; Reference Van den BroekeVan den Broeke and others, 1999). NARE fieldworkers in collaboration with Dutch scientists also recovered two medium-depth ice cores on Amundsenisen in 2000/01 (Reference HofstedeHofstede and others, 2004) and on Fimbulisen (Reference KaczmarskaKaczmarska and others, 2004). BAS expeditions retrieved two shallow firn cores in 1986/87 (Reference Wolff and SuttieWolff and Suttie, 1994) and 1987/88 (Reference Mulvaney and WolffMulvaney and Wolff, 1993), the only records available in Coats Land, the southwestern corner of the area of investigation. Ten years later, in 1997/98, BAS in a joint venture with NARE participated in EPICA studies recovering three shallow firn cores (this paper) and one medium-depth ice core on Amundsenisen (Reference HofstedeHofstede and others, 2004). The latest available record, at the time of writing, results from AWI firn-core drilling near Heimefrontfjella in 2003/04 (this paper).
4. Spatial and Temporal Validity of Data
4.1. Annual local noise and small-scale spatial variability
Spatial and temporal variations in accumulation are attributed to the combined effects of: (1) air-mass transport pathways; (2) interaction of airflow with topography;and (3) local effects, i.e. sastrugi and redistribution by snowdrift. Post-depositional processes driven by wind reshape the primary precipitation pattern leading to preferential accumulation of material in topographic depressions, in contradiction to the lowered accumulation rates found on crests (e.g. Reference SommerSommer and others, 2000; Reference Arnold and ReesArnold and Rees, 2003; Reference FrezzottiFrezzotti and others, 2005). Complex local wind erosion and redeposition processes occur due to undulations in surface topography at the kilometre scale and more specifically on variations of the downslope surface gradient (e.g. Reference Melvold, Hagen, Pinglot and GundestrupMelvold and others, 1998). Reference King, Anderson, Vaughan, Mann, Mobbs and VosperKing and others (2004) observed large spatial variability in snow accumulation (20–30% on a scale of about 1 km) even over gentle topography, which was attributed to a highly non-linear relationship between wind speed and snow transport. Stake-line readings and ice-penetrating radar (IPR) studies provide continuous accumulation information at high spatial resolution (Reference Richardson-Näslund, Aarholt, Hamran, Holmlund and IsakssonRichardson-Näslund and others, 1997; Reference Richardson-NäslundRichardson-Näslund 2004; Reference Rotschky, Eisen, Wilhelms, Nixdorf and OerterRotschky and others, 2004; Reference Spikes, Hamilton, Arcone, Kaspari and MayewskiSpikes and others, 2004). For instance, Reference Isaksson and KarlénIsaksson and Karlén (1994) found a standard deviation of 70% of the average annual accumulation using stake measurements at Ritscherflya. Such results raised questions on how representative single-point records are at a regional scale. Their spatial scale of significance has been widely discussed and is a critical point when interpolating over large distances.
Snow deposition patterns that are highly variable locally are preserved in the form of considerable interannual variability generally detected in firn-core stratigraphy. This phenomenon has been discussed in a large body of literature (e.g. Reference GoodwinGoodwin, 1991; Reference FrezzottiFrezzotti and others 2004, Reference Frezzotti2005). Specifications for year-to-year variations range from 30% on Amundsenisen (Reference SommerSommer and others, 2000) to as high as 250% for firn cores drilled on Jutulstraumen (Reference Melvold, Hagen, Pinglot and GundestrupMelvold and others, 1998). High-frequency variability in annual-layer thickness within ice cores is associated in the first place with wind redistribution controlled by micro-relief, rather than with climatically forced precipitation changes (Reference Reijmer and van den BroekeReijmer and Van den Broeke, 2003; Reference FrezzottiFrezzotti and others, 2005). As a consequence, accumulation time series generally miss a clear correspondence even for nearby firn-core sites (Reference SommerSommer and others, 2000). Such findings indicate that only multi-year accumulation averages should be used for regional studies. Five to ten years is generally considered sufficient for providing reliable estimates of local snow accumulation (Reference Petit, Jouzel, Pourchet and MerlivatPetit and others, 1982; Reference MorganMorgan, 1985; Reference Mosley-Thompson, Kruss, Thompson, Pourchet and GrootesMosley-Thompson and others, l985). We account for this by excluding all records that cover <5 years of accumulation, in order to avoid biases in our interpolation map resulting from annual local noise and small-scale spatial variability. When a choice was available between several measurements taken at the same location, we decided on the one with the longer average time interval.
4.2. Multi-decadal accumulation trends
Both positive and negative long-term trends in snow accumulation have been reported over the Antarctic ice sheet (e.g. Reference Petit, Jouzel, Pourchet and MerlivatPetit and others, 1982; Reference Morgan, Goodwin, Etheridge and WookeyMorgan and others, 1991; Reference Mosley-Thompson, Bradley and JonesMosley-Thompson, 1992). Our dataset consists of records that are sampled irregularly in space and time, covering a wide variety of accumulation time periods. Average time intervals range from several years to several hundreds of years, with the majority (78 out of 111) <50 years. Sampling was carried out over a time period of roughly 40 years (see Appendix). Therefore, differences in accumulation values might partly originate from a changing accumulation rate over time. This could induce misinterpretations when modelling spatial variations statistically in order to define optimal interpolation criteria for applying geostatistical surface prediction methods, i.e. kriging (see section 6). Whether this dataset is suitable for compiling an accumulation map must therefore be investigated.
All investigators involved in data sampling agree that time series of accumulation rates and stable isotopes do not indicate a drastic change in atmospheric moisture transport in DML over the covered time periods of their studies (e.g. Reference Isaksson, Karlén, Gundestrup, Mayewski, Whitlow and TwicklerIsaksson and others, 1996; Reference MelvoldMelvold, 1999; Reference Oerter, Graf, Wilhelms, Minikin and MillerOerter and others, 1999; Reference KarlöfKarlöf and others, 2000). At decadal timescales, the records show no clear trend or different magnitude and direction in changes between the single sites. Reference Oerter, Graf, Wilhelms, Minikin and MillerOerter and others (1999) and Reference GrafGraf and others (2002) observed temporal variations over the last two centuries in central DML that were linked to temperature variations derived from records of stable isotopes. However, for the same region, constant accumulation rates characterized the last 2000 years in three intermediate-depth ice cores (Reference SommerSommer and others, 2000). Slight increases or decreases in accumulation rates over several decades were found to be small compared with interannual variability. This is in accordance with the findings of Reference Isaksson, Karlén, Gundestrup, Mayewski, Whitlow and TwicklerIsaksson and others (1996) and Reference HofstedeHofstede and others (2004) who combined information from different sites on Amundsenisen. We found a stable plateau climate by comparing cores M and M150 (see Appendix). Both were drilled at 75°00’S,15°00’E and they cover a time period of 31 and 1033 years, respectively. Despite the different average time intervals, similar accumulation values of 43 and 45kg m–2a–1, respectively, indicate stable accumulation rates for the longest time interval included in this study.
Coastal areas experience stronger climatic fluctuations since they are exposed to the combined effects of changing sea-ice extent and cyclonic activity (Reference Isaksson and KarlénIsaksson and Karle´n, 1994; Reference KaczmarskaKaczmarska and others, 2004). However, since observed changes (e.g. Reference Isaksson, Karlén, Gundestrup, Mayewski, Whitlow and TwicklerIsaksson and others, 1996; Reference MelvoldMelvold, 1999) are generally restricted in time and show alternating signs and no inter-site correlation, a likely explanation is the movement of a borehole site through a complex undulation pattern related to bedrock topography (Reference Melvold, Hagen, Pinglot and GundestrupMelvold and others, 1998; Reference SommerSommer and others, 2000). Depending on ice-flow velocity and wavelength of surface undulations (in the order of some kilometres) (Reference Richardson-Näslund, Aarholt, Hamran, Holmlund and IsakssonRichardson-Näslund and others, 1997; Reference Rotschky, Eisen, Wilhelms, Nixdorf and OerterRotschky and others, 2004; Reference Eisen, Rack, Nixdorf and WilhelmsEisen and others, 2005), a coring site requires several decades to travel from a low-accumulation surface crest to a high-accumulation surface hollow. For this reason an apparent trend may be detected from changes in annual-layer thickness of firn cores, even in stable climatic periods. However, the average accumulation derived from the affected firn core would still be representative for a larger surrounding area. Based on these data, we assume that no major change in accumulation occurred in western DML during the past 1000 years that would significantly corrupt the quality of a spatially interpolated picture of accumulation rates over our area of interest.
5. Exploratory Spatial Data Analysis
Nearly all 111 selected accumulation records result from traverse work conducted by different nations as described above. Therefore, data collection did not follow an optimal sampling procedure. The resulting sampling density varies and remains relatively small considering the well-known high spatial variability of snow accumulation, particularly in mountainous areas. In order to define suitable parameters for a spatial interpolation we first examined the statistical properties of our dataset, i.e. value distribution, spatial variability and global spatial trends.
5.1. Value distribution
Accumulation records range from as low as 19kg m–2a–1 at the south eastern corner of our study area on Amundenisen to a maximum of 491kgm–2a–1 on Jutulstraumen. By looking at the number of records measured within defined accumulation intervals in steps of 20 kg m–2 a–1, we find our dataset is not normally distributed (Fig. 2). A bell-shaped histogram would be required for an optimal interpolation result. Instead, the majority of records are lower than 200 kgm–2 a–1, causing the histogram bars to be clearly shifted to the left. The highest number of records (36) fall into the accumulation interval between 45 and 65 kgm–2a–1, while medium and higher values are underrepresented in our dataset. This skewness is also expressed by a large difference between the median (63.5 kg m–2 a–1) and the mean value (129.49 kg m–2 a–1), which in turn is of the same order as the standard deviation of 124.67 kg m–2a–1. The considerably higher proportion of low data values within the dataset is clearly due to the preferential sampling on the DML plateau area (75 out of 111 records). Data points with values higher than 200 kg m–2a–1 are, without exception, located near the coast and the adjacent lower inland regions north of the nunataks, while nearly all records below this threshold belong to the inland plateau. Already from this frequency distribution two separate accumulation regimes can be distinguished.
Fig. 2. Frequency distribution of accumulation rates. Grey: records sampled on the plateau (generally below 200 kgm–2 a–1); white: coastal records; black: total number.
By looking at the accumulation values separately for coastal and plateau regions, we find a distribution close to normal for the coast, expressed by similar mean and median values (Table 1). The coastal records account for the high standard deviation in our dataset, with a value of 128.2 compared with only 30.0 kg m–2 a–1 for the plateau. The latter is in agreement with the findings of Reference Oerter, Graf, Wilhelms, Minikin and MillerOerter and others (1999). The higher spread in accumulation values for the coastal data is not surprising, considering the complex geographical nature described earlier. However, the relative standard deviation is similar for both regions. For the subset of records located on the plateau a distribution close to normal can be found by performing a logarithmic transformation of our data values. By this we can adjust the coefficient of skewness, a measure of the symmetry of a distribution, which should be near zero for optimal interpolation results.
Table 1. Statistical distribution of accumulation value in kgm–2 a–1
For individual regions within the area of investigation we find distinct differences in accumulation mean and spread (Fig. 3). We are aware that the number of records is a critical factor for such comparisons, particularly for the nunatak areas. Amundsenisen and Wegenerisen exhibit smooth topography and remote location with low spread and values generally lower than 120kg m–2a–1. In contrast, the data are highly variable for coastal regions. The most striking region is Jutulstraumen, which covers the entire range of accumulation values within our dataset. This large spatial variation in accumulation has been attributed to strong erosion and redeposition processes on the ice stream controlled by surface undulations and a complex near-surface wind field (Reference Holmlund and NäslundHolmlund and Na¨slund, 1994). Accumulation extremes also result from increased precipitation between 900 and 1200m elevation due to orographic uplift and adiabatic cooling of humid air masses (Reference Melvold, Hagen, Pinglot and GundestrupMelvold and others, 1998).
Fig. 3. Minima, mean and maxima of accumulation rates for individual regions. The number of respective records is given in brackets.
5.2. Outliers
Outliers within a dataset can affect the performance of a surface prediction dramatically by distorting the variogram modelling associated with the kriging procedure later used for the interpolation. They should therefore be identified prior to the interpolation (Reference BurroughBurrough, 1986). There are no global outliers within our dataset, i.e. no records are very high or very low relative to all of the values. However, at four locations we find sample points within the normal range of the entire dataset, but exhibiting unusually low values relative to their neighbours. Major accumulation deficits are found for core sites Jut_E (72°13’ S, 00°43’ W) and Jut_I (71°31’ S, 01 ° 110 W), both located on Jutulstraumen, which have been linked to strong katabatic wind activity as indicated by wind crusts, sastrugi and nearby blue-ice fields (Reference Melvold, Hagen, Pinglot and GundestrupMelvold and others, 1998). Anomalous low values also characterize core sites B10_11 (72°300 S, 09°06’ W) and ‘km 270’ (72°560 S, 09°410 W), situated on Ritscherflya at about 1000 m a.s.l. From field observations and isotope analysis it was found that these low values are also due to erosion by wind activity and were not caused by lower precipitation rate (Reference Oerter, Graf, Wilhelms, Minikin and MillerOerter and others, 1999). From stake-line readings (Reference Rotschky, Rack, Dierking and OerterRotschky and others, 2006) and IPR measurements (Reference Richardson-Näslund, Aarholt, Hamran, Holmlund and IsakssonRichardson-Näslund and others, 1997) we know that these core sites are situated within a wider zone of very low accumulation rates, which is likely a consequence of strong katabatic winds originating in the Heimefrontfjella range. We can therefore conclude that all the local outliers have been measured correctly and represent real abnormalities, reflecting the high spatial accumulation variability in complex terrain. We expect that measuring errors are not a major problem for interpolation of the dataset.
5.3. Global spatial trends
To identify global spatial trends we examined possible relationships between snow accumulation and distance to the nearest coast, elevation as well as slope gradient and orientation. These variables are often used to explain the general distribution of snow accumulation rates (e.g. Reference Giovinetto, Waters and BentleyGiovinetto and others, 1990; Reference Bernhard, Weibel, Dikau and SaurerBernhard and Weibel, 1999). Terrain information was extracted from the RADARSAT-1 Antarctic Mapping Project (RAMP) digital elevation model (H. Liu and others, http://nsidc.org/data/nsdic-0082.html). The ice front was taken from the Antarctic Digital Database (ADD) version 4, provided by the Scientific Committee on Antarctic Research (SCAR). Using the entire dataset, we did not find a correlation between accumulation rates and the chosen parameters. For instance, both the lowest (46 kg m–2 a–1) and highest (491 kgm–2a–1) coastal values are found within the elevation interval 900–1100 m a.s.l. at a distance of 200–250 km from the coast. Even on the flat ice shelves, accumulation ranges from 228 to 451 kgm–2a–1. Clear cross-correlations between slope gradient and accumulation reported by Reference Melvold, Hagen, Pinglot and GundestrupMelvold and others (1998) or Reference King, Anderson, Vaughan, Mann, Mobbs and VosperKing and others (2004) apply only for limited areas.
Excluding all coastal records from the analysis yields an expected negative cross-correlation between accumulation and elevation as well as distance to the nearest coast each with R = –0.6. As with Reference Giovinetto, Waters and BentleyGiovinetto and others (1990), we can describe those relations using first- and second-order regression functions (Fig. 4). Concerning continentality, the correlation improves when moving gradually away from the coast (Table 2). The best correlation of R = –0.8 was found for the distance interval 700–1000km from the coast, indicating that other reasons for spatial accumulation variations (e.g. surface undulations) become less important inland.
Fig. 4. Global trends. Correlations between accumulation rates and (a) elevation and (b) distance to the coast. Coastal records are shown as black dots; records sampled on plateau regions are shown as triangles.
Table 2. Global trends: correlation coefficients (R) between accumulation rates and elevation as well as distance to the coast
6. Surface Prediction
6.1. Semivariogram modelling
The process of variogram modelling, i.e. variography, is an instrument for characterizing the spatial variability of a dataset. It gives useful information for selecting optimal interpolation parameters (e.g. the mathematical model to be used for computing efficient interpolation weights, the size and shape of the search window, as well as the number of neighbouring measurement sites to be included in the prediction of unsampled sites). In order to estimate accumulation rates for the western DML region, we used kriging as the geostatistical interpolation method. Fundamental to kriging is the assumption that a relationship exists between the location of sites and the measured values, which can be modelled statistically. In theory the difference in values should be the same between any two points that are the same distance and direction apart. Naturally, dissimilarities between records increase with distance since nearby samples tend to be more alike than those farther apart. The degree of spatial autocorrelation is estimated from the data based on a scatter plot of half of the squared difference between pairs of sample values against their spatial separation, h, called the semivariogram cloud. An average of dissimilarities, y*(h), is calculated for all sample pairs falling into a certain distance bin of defined size, i.e. the lag size. The theory concerning kriging is discussed in more detail by Reference BurroughBurrough (1986) and Reference WackernagelWackernagel (1998).
6.2. Interpolation settings
Following the exploratory spatial data analysis, the dataset was separated into coastal and inland sites in order to define a set of domain-specific interpolation criteria. Different statistical properties suggest using two local models, rather than a single global model, for the interpolation, even though for the coastal domain the number of records becomes a critical factor for achieving a stable semi-variogram. The Heimefrontfjella, Kirwanveggen and Mühlig-Hofmann- Gebirge mark the borderline between the data subsets. To improve the western boundary conditions, where data are missing over large distances, we used data from Coats Land for both domains.
The selection of a proper set of interpolation parameters is highly subjective. However, we found that the resulting patterns remained quite stable for reasonable values of tested interpolation parameters. Even the specific variogram model type, used for computing interpolation weights, seems to be of minor influence. Over the smooth plateau region, where accumulation variability is relatively small, we selected a spherical model type which shows a progressive decrease of spatial autocorrelation until some distance, beyond which the autocorrelation becomes zero. For the coastal domain, a circular model type with a steeper curve near the origin, giving more weight to the closest neighbours for each prediction, is more appropriate. Variations in lag size and number did not significantly affect the output surface. A lag size of 20 km was chosen for both domains in order to get representative averages for semivariogram bins.
Concerning the number of measured values to be used in a prediction, stable results are computed including at least two to five neighbours in all four directions. Any further increase in the number of neighbours has no impact on the resulting pattern because of the decreasing influence of remote points. The specific shape of the neighbourhood determines where to look for data points. A clear indication for selecting optimal search radii is given by the range of the semivariogram cloud defining the distance to which data are statistically correlated (Reference BurroughBurrough, 1986). Higher autocorrelation in one direction than in another might result from global trends and/or anisotropy. Our earlier spatial data analysis showed a first-order trend on the plateau running from northwest to southeast towards increasing elevation and distance to coast, which was removed before variogram modelling and added back before predictions were made. From variography we still find anisotropic autocorrelation for this region, i.e. values in northeast–southwest direction change more slowly than perpendicular to this. Such a phenomenon can be explained from topography and resultant predominant wind direction, as modelled for example by Reference Van Lipzig, Turner, Colwell and van den BroekeVan Lipzig and others (2004). Therefore, we assume that these directional differences are not a measurement artefact resulting from the preferential alignment of records along traverse lines which could lead to insufficient sampling density for robust estimates in all directions. Consequently, we accounted for this anisotropic spatial autocorrelation by applying varying search radii as described below.
In general, we found a single plateau record to be representative for a larger area than coastal records. This is in agreement with other findings (Reference Stenberg, Hansson and HolmlundStenberg and others, 1999). The spatial autocorrelations suggest that the search radii should be 200–250 km from northeast to southwest and 100–150km perpendicular to this direction. For the coastal domain, where major variation occurs over shorter distances, the search radius was kept as small as possible (80–125 km in any direction). Little variations in search radii led to changes in the predicted surface, particularly in the coastal area. Since the most efficient search radius remains a subject of speculation, a good solution is to produce an average composite of ten realizations, each time adjusting the search window by slightly enlarging its two axes (i.e. the search radii) in 5 km steps. The differences between these predictions are used as a measure of uncertainty in the predicted surface.
Having computed the experimental variograms and fitted model parameters, our accumulation records were interpolated using ordinary kriging for the coast and universal kriging for the plateau domain. The latter should be applied only if an apparent spatial trend in the data can be explained from a priori knowledge (Reference BurroughBurrough, 1986). Using a logarithmic transformation makes the data more normally distributed in case of the plateau domain. Here, the interpolation was carried out using the logarithmically transformed data which were then transformed back. Our interpolation settings are summarized in Table 3.
Table 3. Interpolation settings
7. Results and Discussion
7.1. Surface prediction
The resulting accumulation distribution is mapped in Figure 5. Clearly, predicting a surface from sparse point data gives only a smoothed picture of the underlying reality. Details of spatial accumulation variability cannot be reproduced and the resulting pattern is highly dependent on the random distribution of sample points. For the coastal area, high accumulation alternates with accumulation lows reflecting a spatial variability in precipitation and snow transport that is topographically forced. Isolated circular patterns around single records result from values that deviate significantly from their surrounding neighbours. On the plateau, a clear trend of decreasing accumulation with increasing distance to coast and elevation from northwest to southeast (as described above) is reproduced. Values range from 150 kg m–2 a–1 just south of the mountain ridges down to as little as 25kg m–2a–1 near the southeastern corner of our study area. Tongue-shaped patterns of higher accumulation advance from coastal regions further inland, where nunatak barriers are absent over larger distances (e.g. south of Jutulstraumen and Coats Land). Some problems are visible near Mühlig-Hofmann-Gebirge where predicted values at the northern margin of the plateau domain increase rapidly over short distances. In this area, computed accumulation rates are unreliable. Further interpolation difficulties occur where in situ data are missing over large distances between traverse lines running inland from the coast. The resulting striped patterns are clearly an artefact of the interpolator and not representative of the true accumulation.
Fig. 5. Contoured accumulation pattern for western DML (data are available from the PANGAEA® website http://doi.pangaea.de/10.1594/PANGAEA.472297). A greyscale as well as contours at 50 (coastal domain) and 25 kgm–2 a–1 (plateau domain) are used to indicate the broad spatial variability in snow accumulation across the predicted area. Elevation contours in 500m intervals are shown as dashed white curves. The borderline between coastal and plateau domain follows the mountain chains and is shown as solid grey line.
Due to interpolation uncertainties as a result of the sparse and uneven data-sampling scheme, estimating the total surface mass balance in western DML is difficult. In spite of the larger area, the total accumulation on the plateau (58 × 104 km2) sums up to only 38 Gta–1 vs 62 Gta–1 for the coastal region with a spatial extent of 24 × 104km2. Discrepancies in total accumulation between our ten slightly varied realizations are marginal (generally <1%). In Table 4 we compare our results with previous accumulation compilations by Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others (1999) and Arthern and others (2006), who both guided the interpolation using background fields of satellite microwave observations, as well as by Reference Huybrechts, Steinhage, Wilhelms and BamberHuybrechts and others (2000), who applied a spline interpolation to nearly the same input data as in the present study. While the results of Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others (1999) and Reference Huybrechts, Steinhage, Wilhelms and BamberHuybrechts and others (2000) are 18% and 16% larger, respectively, the recently published findings by Arthern and others (2006) yield about 13% less accumulation for the western DML area than this study. These differences result from significantly higher estimates for the plateau in the case of Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others (1999) and for the coastal area in the case of Reference Huybrechts, Steinhage, Wilhelms and BamberHuybrechts and others (2000). Conversely, Arthern and others (2006) yield similar estimates for the plateau;however, their coastal estimates are significantly lower than ours. Generally the same input data have been used in all of the studies, which suggests that differences occur due to approaches to interpolation. In general, spline interpolation may produce local artefacts of excessively high or low values if input values change dramatically over short distances. Therefore, the results of Reference Huybrechts, Steinhage, Wilhelms and BamberHuybrechts and others (2000) are unreliable for the DML coastal area. The main reason for the disagreements is the absence of a barrier between the coastal and plateau accumulation regimes in the earlier studies, which results in an overestimation of the area of validity of high-accumulation coastal records in the case of Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others (1999) and of low-accumulation plateau records in the case of Arthern and others (2006). Like those of Arthern and others (2006), our results indicate that snow accumulation in western DML is probably smaller than previously assumed.
Table 4. Total accumulation in western DML: comparison between this study and earlier published results
7.2. Interpolation accuracy
In general, accuracies vary across a predicted area with data-sampling density. To evaluate the uncertainty of the interpolated surface, an error estimate is required. In the absence of an extra set of independent observations randomly distributed across the predicted area, we cannot quantify the accuracy directly. From the relatively small number of available data, a subset for validation purposes could not be created before interpolation. Each point in the dataset is essential for the variogram modelling and original prediction.
Some measure of certainty is provided as a by-product of the kriging procedure, known as the kriging variance. In the literature, this quantity has been criticized as being ineffective and a poor substitute for a true error (e.g. Reference JournelJournel, 1986; Reference Chainey, Stuart and CarverChainey and Stuart, 1998). The kriging variance, depending only on the geometrical arrangement of the sample data points, simply states that accuracy decreases with growing distance from input data.
As a measure of confidence we instead used the differences between the ten interpolation realizations, each with slight adjustments in the applied search window, as described above. The lower the standard deviation of the different predictions, the more accurate the interpolation is likely to be in terms of self-consistency of the model applied to the data. The standard deviations across the predicted area are mapped in Figure 6, expressing a clear difference between coastal and plateau domains. On the plateau, data density seems sufficient for achieving stable interpolation results, with the exception of the Jutulstraumen inlet area. In contrast, for coastal regions, exhibiting a much larger spatial accumulation variation, major interpolation uncertainties occur for unsampled locations between measurement sites. Clearly, small adjustments in the shape of the search window impact mostly where single values deviate greatly from their surrounding neighbours (local singularities, as described above). This result leads back to the question of how representative single-point records are for their wider surrounding. In order to catch the full range of accumulation variability occurring due to the complex coastal terrain and wind-field conditions, further data acquisition should be carried out in directions other than from north to south. Regional climate modelling may also help in this respect (Reference Van den BroekeVan den Broeke and others, 2006).
Fig. 6. Map of standard deviation between ten predictions with stepwise increase of search radii (data available from the PANGEA website http://doi.pangaea.de/10.1594/PANGAEA.472298). A greyscale and contours at 10 kgm–2 a–1 (black lines) indicate the interpolation inaccuracy. Elevation contours in 500m intervals are shown as dashed white curves. The borderline between coastal and plateau domain follows the mountain chains and is shown as solid grey line.
8. Summary and Conclusions
A set of 111 accumulation measurements from western DML is compiled and presented. Most data result from fieldwork carried out over the last 15 years in the framework of EPICA and ITASE. Older records from the literature are used to supplement the data. In spite of a large variety of covered accumulation periods, ranging between 5 and 1000 years, we find the datasets temporally representative, as no significant trend in accumulation history over the last few centuries has been reported from ice-core records of this region. Uneven data distribution resulting from preferential alignment of measurement sites along traverse lines, which lead mostly inland from the coast, is a more critical factor for judging the spatial accumulation variability. Nevertheless, an interpolated map of mean annual snow-accumulation rate could be derived that is in accordance with presumptions made from known overall spatial trends and topographically forced changes in accumulation. The interpolation was carried out using the geostatistical interpolation method of kriging, based on two local models and domain-specific interpolation settings for coastal and inland sites. Our results suggest that the total accumulation rate for the predicted area is about 16–18% smaller than published in earlier studies. For a quantitative evaluation of the uncertainty in predicted accumulation rates, we calculated the standard deviation of ten interpolation estimations, each with a slight stepwise increase in search radii. A map of the standard deviation thus obtained could help in locating future sampling sites to further minimize interpolation uncertainties. The existing sampling network should be improved with respect to surface topography and associated surface wind-field patterns which in combination have the potential to change accumulation conditions dramatically over short distances, especially in the complex terrain of the coastal regions. These results are of interest for the interpretation of climate records from the study area, as well as for ongoing surface mass-balance modelling aimed at resolving the fine-scale accumulation pattern across western DML.
Acknowledgements
This work is a contribution to the European Project for Ice Coring in Antarctica (EPICA), a joint European Science Foundation (ESF))/European Commission (EC) scientific programme, funded by the European Union and by national contributions from Belgium, Denmark, France, Germany, Italy, The Netherlands, Norway, Sweden, Switzerland and the United Kingdom. The main logistical suppoert was provided by Institut Polaire Français–Emile Victor (IPEV) and Programma Nazionale di Richerche in Antartide (PNRA) (at Dome C) and AWI (at DML). This is EPICA publication No. 175. We thank the scientific editor R. Bindschadler and two anonymous reviewers for valuable comments on the manuscript.
Appendix
Table 5. List of firn- and ice-core (C) and snow-pit (P) measurements sampled over the area of interest (also available from the PANGAEA® website http://doi.pangaea.de/10.1594/pangaea.407654). Records included (i) in the interpolation are marked with a ‘y’ in the first row
Table 5. (continued)
Table 5. (continued)
