INTRODUCTION
Extracting the record of past climate, simulating the evolution of ice sheets and ice-sheet mass balance all require data on accumulation rates (Petit and others, Reference Petit1999; Lenaerts and others, Reference Lenaerts, van den Broeke, van de Berg, van Meijgaard and Munneke2012). Recent accumulation rates can be obtained in situ from stake readings (Ding and others, Reference Ding2011), shallow cores (Oerter and others, Reference Oerter2000), or can be derived from satellite observations (Arthern and others, Reference Arthern, Winebrenner and Vaughan2006). A more challenging task is estimating the palaeo-accumulation rate. Past accumulation rates can be evaluated from established chronologies at ice-core drilling sites. Drilling a deep ice core, however, is expensive and labor intensive, and data are spatially restricted. Radio-echo sounding (RES) surveys provide a good method to acquire spatially extensive accumulation rate data where no ice cores exist by continuous measurements of internal reflections within the ice sheet (Paren and others, Reference Paren and Robin1975). Internal layers in ice are caused by (1) firn density variations in the upper part of the ice sheet, (2) crystal-orientation fabrics and (3) the presence of impurities such as volcanic ash, acids and sea salt (Corr and others, Reference Corr, Moore and Nicholls1993; Fujita and others, Reference Fujita1999). These layers are believed to be isochronous surfaces (Millar and others, Reference Millar1981). The depth of an isochronous layer contains information about two variables intimately linked to each other: the past accumulation rate and the age of the ice. If independent dating of a specific layer is carried out, the past accumulation rate can be inferred from a layer-thinning function.
One of the earliest models that describe the thinning of the annual layers at the summit of an ice sheet is the Dansgaard-Johnsen model (D-J model; Dansgaard and Johnsen, Reference Dansgaard and Johnsen1969). This model was successfully used by Dahl-Jensen and others (Reference Dahl-Jensen, Johnsen, Hammer, Clausen and Jouzel1993) and Fahnestock and others (Reference Fahnestock, Abdalati, Luo and Gogineni2001) to reconstruct past accumulation rates from the sequences of dated layers for the Greenland ice sheet. In Antarctica, Siegert (Reference Siegert2003) used the D-J model to determine the spatial and temporal variation in East Antarctic ice accumulation over the Last Glacial cycle. Similarly, Leysinger-Vieli and others (Reference Leysinger-Vieli, Siegert and Payne2004) reconstructed the ice accumulation history for the region from Ridge B to Vostok station. Huybrechts and others (Reference Huybrechts, Rybak, Steinhage and Pattyn2009) derived new strain-thinning functions from established chronologies and used internal ice layers to infer the spatio-temporal pattern of accumulation rate in Dronning Maud Land. Furthermore, these palaeo-accumulation rates were used to fine-tune the relation between accumulation rate and temperature. In addition to the above studies, the spatial variation of accumulation rate along a flowband can be also recovered from analysis of the shapes of internal layers by simulating a range of accumulation patterns with an ice flow model to calculate sets of internal layer patterns. Inverse techniques can then find which accumulation pattern reproduces the observed shape of internal layers from RES measurements (Nereson and others, Reference Nereson, Raymond, Jacobel and Waddington2000; Waddington and others, Reference Waddington, Neumann, Koutnik, Marshall and Morse2007; Hindmarsh and others, Reference Hindmarsh, Leysinger-Vieli and Parrenin2009; Macgregor and others, Reference Macgregor2009; Leysinger-Vieli and others, Reference Leysinger-Vieli, Hindmarsh, Siegert and Sun2011).
Dome A (80°22′63″S, 77°22′90″E) is the highest plateau (4092.46 m a.s.l.) of the Antarctic ice sheet, and it has been identified as an important location to drill a deep ice core (Xiao and others, Reference Xiao2008; Sun and others, Reference Sun2009, Reference Sun2014; Zhang and others, Reference Zhang2014). There are several publications about recent accumulation rates at Dome A (Hou and others, Reference Hou, Li, Xiao and Ren2007; Ding and others, Reference Ding2011; Wang and others, Reference Wang2013), but little is known of accumulation rates deeper in the past. Sun and others (Reference Sun2014) modeled the age of the ice beneath Dome A, and we use these age-depth scales as part of the basis for the dating here of six layers (spanning 30–161 ka BP) that are linked to the Vostok ice core along the continuous RES profile between Dome A and Vostok. We then infer the spatial and temporal pattern of past accumulation rates (ice equivalent accumulation) at Dome A using the depth-age measurements of the six isochrones.
FIELD DATA
The radar data profiles used in our study are parts of two surveys: one is part of the data collected using a dual-frequency radar system over the Dome A region by the Polar Research Institute of China (PRIC) during the 21st Chinese National Antarctic Research Expedition (CHINARE) (the black line in Fig. 1) in 2004/05; the other is the continuous RES profile of ~1300 km length connecting the ice-core drilling sites of Vostok and Dome A, which was flown by the Alfred-Wegener-Institut (AWI) in 2007/08 as part of the Dome Connection East Antarctica (DoCo) project (the red line in Fig. 1). The PRIC survey intersects the AWI survey at two points: S (80°21′27.41″S, 77°22′46.54″E) and N (80°20′31.30″S, 77°25′56.60″E), as marked in Figure 1.
MODELS
The D-J model was used to calculate ice accumulation rates from dated radar layers in the Dome A region over the last ~161 ka. The D-J model assumes a constant strain rate from the surface down to some height z = h above the bed, and from there decreases linearly to zero at the base. The ice thickness (H) is assumed fixed over time. The modeled age t of the ice layer at depth d = H − z is given by Eqn (1):
where z and H are the present and initial distance (equivalent to the ice thickness) from the bed (z is positive upwards from the bed). $\bar b$ is the average accumulation rate we are interested in.
The accumulation for each time slice is given by:
where $\bar b_1 $ and $\overline {b_2} $ are the average accumulations obtained by the D-J model for layers to the depth of the respective ages t 1 and t 2.
The D-J model requires a value for h, which here is unknown. Different values of h have been used in different applications of the D-J model. Dahl-Jensen and others (Reference Dahl-Jensen, Johnsen, Hammer, Clausen and Jouzel1993) determined h to 1255 m, Fahnestock and others (Reference Fahnestock, Abdalati, Luo and Gogineni2001) estimated h at 400 m, Leysinger-Vieli and others (Reference Leysinger-Vieli, Siegert and Payne2004) took a mean value of h = 484 m. Here, we choose the parameter h = 800 m. We defer the details of this to the methods section.
METHODS
Internal radio-echo layering can be identified and traced across the ice sheet from analysis of radargrams. The depth of a given internal layer can be calculated from the time delay and the velocity of radio-waves. The wave velocity varies according to the state of the ice (density, temperature, fabric, etc.) at any location. Kovacs and others (Reference Kovacs, Gow and Morey1995) provide an empirical formula for radar velocity in firn of various densities taking permittivity variations of ice into consideration. An estimate of the vertical density distribution appropriate for a given site enables a ‘firn correction’ Z f, to be added to a preliminary calculation of the depth based on the assumption of a constant velocity appropriate for solid ice (V i = 168.9 m μs−1 in our study), see Appendix for details. The firn correction for PRIC data is 15 m, which is generally applicable over the central Antarctica (Dowdeswell and Evans, Reference Dowdeswell and Evans2004). The firn correction for AWI data, however, is 16.22 m derived from the density profiles measured at the Vostok drill site (D. Steinhage, 2015, personal communication). Hence, there is a small difference (1.2 m) in firn corrections for the two datasets.
Six layers can be continuously traced from Vostok along the AWI radar profile to the PRIC intersection points (S and N). Thus the depth-age scale of each layer at the location of S and N can be calculated (Fig. 2a). We traced all the clearly identifiable internal layers from PRIC radargrams, to ensure that crossover errors between these two radar lines were minimized. We calculated depth of each layer in the PRIC survey and selected the layers based on their proximity to AWI layers at the locations of S and N. The different depths of layering found in the PRIC and AWI radar surveys may be explained by considering the complex nature of the radar echo, which is a superposition of echoes from various physical impedance contrasts in the ice as well as from various points within the fresnel zone of the radar beam (Moore, Reference Moore1988). This means differences in radar centre frequency, bandwidth and pulselength as well the radar hardware such as beam width, antenna focusing and receiver hardware processing of the received radar signal can all affect the appearance of layering in the final radargram. The six layers are numbered according to their age, and the depth differences of each layer between PRIC and AWI survey are shown in Figure 2b. Moreover, we give the range of ages for each layer with the depth error. The depth-age reference data is Vostok Ice Core Data for 420 000 a (Petit and others, Reference Petit1999).
Previous research shows that the appropriate value of h lies between ~0.2H and 0.5H (Cuffey and Paterson, Reference Cuffey and Paterson2011). Choosing the summit point of Dome A, H = 2200 m as an example, and limiting 0.2H < h < 0.5H, we show the effects on reconstructed accumulation rates in Figure 3.
Figure 3 shows that the reconstructed accumulation rates increase gradually as h is increase from 400 to 1200 m. As h increases by 100 m, the calculated accumulation rate for six layers will increase on average by 0.19, 0.28, 0.35, 0.45, 0.86 and 0.99 mm a−1, respectively (shown with different coloured lines in Fig. 3). These variations are relatively small compared with the values of inferred accumulation rates (0.90%, 1.20%, 1.45%, 2.50%, 4.00%, 4.80% percentages of the inferred accumulation rates, respectively). Thus over a reasonable range of h there is minimal impact on reconstructed accumulation rates. We also find that as h increases from 800 to 1200 m, the rate of accumulation change increases, especially for the two oldest layers. Therefore a choosing h = 800 m is reasonable in our study given the steady-state assumptions in the D-J model.
RESULTS AND DISUSSION
Spatial and temporal changes in ice accumulation
Six layers are traced from each radargram in the Dome A region and dated by linking them to the Vostok ice core site. The ages of these isochrones from shallow to deep are 34.3 ± 1.3, 39.6 ± 0.1, 47.5 ± 1.7, 93.3 ± 0.4, 123.5 ± 1.5 and 161.4 ± 1.0 ka, respectively. Using the depth-age scale data of dated layers as input with a parameter of h = 800 m, we reconstructed the pattern of past accumulation rates at Dome A using Eqn (1). Here, we use the reconstructed results along the radar profile with a length of 216 km at Dome A to analyze the spatial distribution of accumulation for each time period (Fig. 4).
The mean accumulation rates along the transect are 2.21, 2.44, 2.62, 1.95, 2.46, 2.07 cm a−1 for TP1, TP2, TP3, TP4, TP5 and TP6, respectively. The inferred accumulation rates vary spatially over 40–50 km, so we smooth past accumulation rates with a 40 km length filter from south to north along the transect. There is an average increase of 3.6 mm a−1 between the first parts (0–40 km) and second parts (40–80 km) of the profile, and an average decrease of 2.2 mm a−1 from the second part to the third part (80–120 km). Further north (130–180 km), past accumulation rates increase slightly again. However, spatial variability across the whole transect is relatively consistent over time as seen by the similarity in variability of all the layers. The relatively consistent spatial patterns at Dome A suggest that the geometry of the ice-sheet surface has not changed in gross form over the last ~161 ka (Siegert and Payne, Reference Siegert and Payne2004).
To investigate if the spatial pattern in ice accumulation has changed with time, we calculated the averaged accumulation rates for five different time slices along the RES profile (Fig. 5). The accumulation for each time slice is determined using Eqn (2). The mean accumulation rates are 4.20, 3.41, 1.29, 3.74, 1.10 cm a−1 for TP12, TP23, TP34, TP45 and TP56, respectively. The conspicuous feature shown in Figure 5 is the different accumulation rates between TP12, TP23, TP45 and TP34, TP56, in which the latter two are lower than former three. Calculated results for TP12, TP23 and TP34 show how the accumulation changed through time since the last glacial period (13–116 ka ago). From the comparison of mean accumulation rates between TP2, TP3 and TP12, TP23, it is clear that the latter two are larger than the former, indicating that there were relatively high accumulation rates at Dome A for the time period of ~34–47 ka ago. We compare our results to Lorius and others (Reference Lorius1985), a study that derived a 150 ka climatic record from Vostok ice core. The successive warm and cold stages are designated by letters A-H consecutively from the top of the record downwards in Figure 6. As can be seen in the Figure 6, stage C (30–58 ka ago), which is roughly equivalent to Marine Isotope Stage 3 (MIS3: 29–59 ka ago), was a relatively warm interstadial, an intermediate stage between full glacial and interglacial conditions, which may explain why Dome A experienced relatively high accumulation for the time period ~34–47 ka. While temperature might play a role in accumulation rates, it might not be the only story since the most recent warm interval does not appear to have high accumulation rates.
Comparison with previous studies of Dome A
There are several studies on present-day accumulation rates at Dome A. Hou and others (Reference Hou, Li, Xiao and Ren2007) calculated recent accumulation rates at Dome A as 2.3 cm w.e. a−1 based on the β radioactivity horizon and the density profile. Ding and others (Reference Ding2011) studied the overall spatial characteristics of the surface mass balance (SMB) using stake measurements along a traverse route that starts at Zhongshan and ends at the Dome A summit. They find that the Dome A zone has the lowest snow accumulation rate (3.5 cm a−1, 2005–08). Wang and others (Reference Wang2013) inferred the spatial and temporal variability of snow accumulation near Dome A from snow pit and stake measurement data. The results suggest that the snow accumulation rate shows large interannual variations, but stable multi-decadal averages over the last 7 centuries. The calculated past accumulation rates at Dome A in our study are within a factor of 2 of these recent accumulation rates, although the lowest past accumulation rate for layer 4 was below 2 cm a−1. Accumulation rate variability between glacial and interglacial periods varied by a factor of 2 in the Vostok ice core (Lorius and others, Reference Lorius1985).
Studies about past accumulation rate at Dome A are nearly absent. Siegert (Reference Siegert2003) calculated the mean accumulation rates (1.6 cm a−1, Fig. 6 red dashed line) for the past 94 ka along a radar line ~100 km from the summit of Dome A, which suggests an accumulation rate ~0.27 cm a−1 lower than the result in our study (1.87 cm a−1 for TP4. Figure 6 red solid line). They also inferred accumulation rates using radar layer depths along a radar line of length 200 km that crossed over Ridge B and Vostok station (the corresponding results are indicated in Fig. 6 by the blue solid line). The study used a layer of a similar age (45.9 ka) and time period (45.9–83.5 ka) to layer 3 and TP34 in our study, respectively. These results suggest that the mean accumulation rate along the transect in our study is 0.38 cm a−1 more than the average rate over the transect between Ridge B and Vostok for the last ~46 ka. The rate of ice accumulation in our study (1.28 cm a−1 at Dome A zone), however, is lower than that in Siegert's study for the time period ~46–90 ka (1.50 cm a−1). Both Siegert's and our study show a similar temporal pattern that the accumulation rates for the past ~60 ka are higher than that for ~46–90 ka ago.
Discussion
The D-J model is a simple one-dimensional flow model and by using it to calculate the accumulation rate, we minimize the effect of the horizontal flow. The use of D-J model is most appropriate near an ice divide, where slow or negligible horizontal flow means that all the ice at depth originates from near the same location on the surface. The horizontal surface velocity values are close to zero at the summit of Dome A (11.1 ± 2.4 cm a−1; Yang and others, Reference Yang2014) and increases with distance from the summit. The horizontal velocity increases from 1.3 m a−1 at site of 150 km distance from the summit of Dome A to 3.0 m a−1 at site of 230 km distance from the summit of Dome A (Zhang and others, Reference Zhang2008). Macgregor and others (Reference Macgregor2009) infer accumulation rates from three radar layers (26, 35 and 41 ka old) in the Vostok Subglacial Lake region using two methods: one is the same as the method used in our study and the other is a combination of steady-state flowband modeling (model inputs include ice-surface elevations, ice velocities and temperatures) and formal inverse methods. By comparing the results of two methods, for all three layers along one flowband, which crosses the Vostok station (surface-velocity is ~2 m a−1), the mean difference between two results is 4%, which is a relatively small considering other uncertainties in the method. Sun and others (Reference Sun2014) calculate the basal temperatures around Dome A; for a geothermal heat flux of 50 mWm −2, only a small fraction of the basal ice in the 30 × 30 km2 domain centered around Kunlun station is at the pressure melting point. Consequently, we do not expect the absence of horizontal shear strain and basal melt in the D-J model to affect the results of spatial and temporal pattern of accumulation rate at Dome A in our study. However heat flux is largely unknown in the region, and larger areas of the bed are at pressure-melting point for higher heat fluxes.
Bell and others (Reference Bell2011) showed that around Dome A, 24% of the area contains frozen-on ice features, which drive substantial mass redistribution from localized basal melt and subsequent re-freezing at the bottom of the ice sheet. These features clearly show that regions of the basal ice must be at pressure-melting point. A freeze-on ice structure was also found in the PRIC RES profile along the traverse from Zhongshan station to Kunlun station (Tang and others, Reference Tang2015) located 1044–1056 km from the coast. No freeze-on ice structures were however found within the transect in our study. Sun and others (Reference Sun2014) also showed the importance of the ice fabric in determining the ice vertical velocity, and hence the age of layers with depth. This effect is not included in the D-J model, however, the effects of fabric are typically significant only much closer to the base of the ice than in the relatively shallow ice containing the 6 RES isochrones analyzed here.
Waddington and others (Reference Waddington, Neumann, Koutnik, Marshall and Morse2007) defined a non-dimensional depth number D, to express the relationship between the horizontal distance that a particle has traveled from the surface to any layer, and the characteristic length scales of accumulation-rate and ice-thickness variability. The equation for calculating D is defined as follows:
where $\bar U$ is the temporally averaged horizontal velocity of a particle as it following a particle path, ${\rm {\cal A}}$ is the age of the layer. $1{\rm /}L_{\dot b} = \left\vert {(1{\rm /}\dot b)({\rm d}\dot b{\rm /d}x)} \right\vert$ and $L_{\dot b} $ is the length over which $\dot b$ changes by an appreciable fraction; Similarly, $1{\rm /}L_H = \left\vert {(1{\rm /}H)({\rm d}H{\rm /d}x)} \right\vert$ and L H is the length over which ice thickness changes by an appreciable fraction.
If D ≪ 1 the spatial gradients in accumulation rates and ice thickness do not significantly affect the layer depths, and hence the D-J model is sufficiently accurate to infer accumulation rates. The values of D for the upper four layers in our study are both <0.09, and the D values of the other two layers are ~0.1. This suggests that even the deepest and oldest layers studied here are not significantly affected by spatial variability in accumulation and that the D-J model is acceptable to infer these accumulation rates at Dome A region.
CONCLUSIONS
Depth-age relationships are extended from Vostok ice core to Dome A. Six dated radar isochrones layers are traced from a radar profile collected in the Dome A region. The D-J model, using the depth-age measurements as input, is then used to infer the spatial and temporal variations in Dome A accumulation rates. The main result of this investigation is a generalized view of the spatial and temporal variation in ice accumulation over the last ~161 ka.
The lowest average ice accumulation rates for each time window were calculated at Dome A (smoothed over 40 km), which are 2.09, 2.30, 2.41, 1.79, 2.21 and 1.88 cm a−1 for TP1, TP2, TP3, TP4, TP5 and TP6, respectively. The spatial variability shows an increase or decrease of 0.7–4.0 mm a−1 on a scale of 40 km (Fig. 4). As a whole, ice-accumulation rates increase slightly from south to north along the radar profile over a distance of 216 km, starting from the summit of Dome A.
At Dome A, the mean accumulation for TP2 and TP3 is higher than that for TP1 and TP4, which suggests that Dome A experienced relatively high accumulation 34–47 ka ago, roughly corresponding to the warm interstadial Marine Isotope Stage 3. The variability of accumulation rates is almost a factor of 4 on scales of tens of ka (see Fig. 5, TP12 and TP56). This is probably more than would be expected from climatically driven changes in accumulation rates. This suggests that other factors such as the steady-state assumptions in the modeling may be questioned. However, despite complicating factors known to exist in ice around Dome A, such as basal accumulation of ice and varying ice fabric with depth, tests on the suitability of the simple 1 dimensional flow modeling method for the relatively shallow radar isochrones suggest that it can provide reasonable estimates of accumulation rate variability over the past 160 ka.
This paper is a first attempt to calculate past accumulation rates for Dome A, East Antarctica. A sophisticated 3-D flow model for the Dome A region has already been developed (Sun and others, Reference Sun2014), and with advances in basal hydrology incorporated within it, a more accurate and longer time series of past accumulation rates at Dome A will be derived. This work can also be independently verified with data from the Dome A deep ice core in the near future.
ACKNOWLEDGEMENTS
Work was supported by China's National Key Science Program (2012CB957702, 2013CBA01804), the National High-tech R&D Program of China (2011AA040202) and the Chinese Polar Environmental Comprehensive Investigation and Assessment Programs (CHINARE-02-02), Foreign Cooperation Support of Chinese Arctic and Antarctic Administration, SOA, China (IC201214), National Science Foundation of Shanghai, China (13ZR1445300), National Natural Science Foundation of China (41376192).
Appendix
FIRN CORRECTION
Since the conductivity of dry polar firn and ice is extremely low, the effective phase velocity of high-frequency electromagnetic waves can be determined from
where c is the velocity of an electromagnetic wave in a vacuum, 0.3 m ns−1. ε′e is the relative effective dielectric constant.
Kovacs and others (Reference Kovacs, Gow and Morey1995) combined published field data for electromagnetic wave velocities in polar firn. They found that an empirical linear relationship between the refractive index n and density ρ, that had been developed much earlier (Robin and others, Reference Robin, Evans and Bailey1969), produced the best statistical fit:
since n 2 = ε′e in the radio-radar frequency range of interest, then
Robin and others (Reference Robin, Evans and Bailey1969) adopted the coefficient K = 0.851 × 10−3 m3 kg−1, and Kovacs and others (Reference Kovacs, Gow and Morey1995) made a small adjustment to the coefficient K, K = 0.845 × 10−3 m3 kg−1.
An estimate of the vertical density distribution appropriate for a given site enables a correction (firn correction z f) that will be added to a preliminary calculation of the depth based on the assumption of a constant velocity appropriate for solid ice (Dowdeswell and others, Reference Dowdeswell and Evans2004).
where t is the two-way vertical travel time, i – ice, f – firn. Here the density of solid ice is ρ i = 917 kg m−3. Kovacs and others adopted a value of ${\epsilon} _{\rm i} {\rm =} n_{\rm i}^2 {\rm =} 3.15$ for solid ice.