1. Introduction
In Antarctica, the accumulation rate is the primary, and in East Antarctica the only, component in the calculation of the surface mass balance. There is some ablation in particular areas of the East Antarctic ice sheet, mainly due to sublimation from blue-ice fields close to outcropping rocks, but this plays only a minor role in the overall surface mass balance of the ice sheet. Post-depositional redistribution of snow (snowdrift), on the other hand, is an integral component of the net accumulation rate and cannot be distinguished from the local atmospheric snow input from glaciological field data alone. In general, accumulation patterns are not smooth, either in time or in space (Reference EisenEisen and others, 2008).
Direct measurements of precipitation in Antarctica are carried out at a number of stations dispersed over the continent, including automatic weather stations. Annual accumulation rates can be obtained in situ from stake readings, snow pits and shallow cores (e.g. Reference Parrenin, Rémy, Ritz, Siegert and JouzelOerter and others, 1999, Reference Reeh, Oeschger, Langway and Jr2000; Frezzotti and others, 2005), or can be derived from ice-penetrating radar profiles (e.g. Reference Siegert and FujitaRotschky and others, 2004; Reference Waddington, Neumann, Koutnik, Marshall and MorseSpikes and others, 2004). More recently, new methods based on satellite observations have been introduced (e.g. Reference Arthern, Winebrenner and VaughanArthern and others, 2006). Compilations of available accumulation estimates combining different methods have resulted in gridded datasets of the present-day accumulation rate over large parts of Antarctica (Reference Wesche, Eisen, Oerter, Schulte and SteinhageVaughan and others, 1999; Reference Hindmarsh, G.J.M., Raymond and GudmundssonHuybrechts and others, 2000; Reference Arthern, Winebrenner and VaughanArthern and others, 2006).
A more challenging task is to estimate the palaeo-accumulation rate. At ice-core drilling sites, past accumulation rates can be evaluated from established chronologies provided the thinning function (ratio of a layer to its initial thickness) is known. For a spatial coverage, ground- or ice-penetrating radar provides continuous measurements of internal reflections in the near-surface firn layer where ice ages may be up to ~1000 years. Deeper layers can be identified by airborne radio-echo sounding (RES). Internal layers in ice are caused by ice density variations in the upper part of the ice sheet (e.g. Reference Steinhage, Nixdorf, Meyer and MillerSiegert, 2003) and by the enclosure of volcanic ashes and acids (Fujita and others, 1999; Reference Vaughan, Bamber, M.B., Russell and CooperSiegert and Fujita, 2001). They can be assumed to be isochronous surfaces. Theoretical aspects of inferring accumulation-rate patterns from such deep layers in ice sheets were discussed by Reference WilhelmsWaddington and others (2007) based on geophysical inverse theory.
The depth of an internal layer contains information about two variables intimately linked to each other: the past accumulation rate and the age of the ice. Neither of these variables can be estimated or calculated independently directly from field data. To separate the age of the ice and the accumulation rate, an additional relation is needed that specifies how the thickness of the original annual layer has changed by both firn densification and strain thinning. If independent dating of a specific layer is carried out (e.g. at a place of intersection with an ice core with established ice chronology), the past accumulation rate can be inferred from the layer-thinning function. This is most often done locally using an approach originally suggested by Reference Dansgaard and JohnsenDansgaard and Johnsen (1969) (Reference Steinhage, Nixdorf, Meyer and MillerSiegert, 2003; Reference Oerter, Graf, Wilhelms, Minikin and MillerLeysinger Vieli and others, 2004). However, such an approach cannot deal with the effects of complex flow regimes conditioned by bedrock topography, spatial and temporal variations of flow velocity, and spatial variations of the accumulation rate. In addition, in order to link palaeo-accumulation rates to the present-day accumulation regime, it is often necessary to take into account the effect of horizontal advection and thus to trace the layer particles back to their time and place of deposition.
Our approach to reconstructing the past and present accumulation rate from internal reflection horizons is to derive the strain-thinning function from a flow history calculation with a comprehensive nested high-resolution higher-order ice-flow model applied to eastern Dronning Maud Land, East Antarctica (Reference HuybrechtsHuybrechts and others, 2007). The main assumption to make our method work is that the thinning function is determined mainly by ice-dynamic features (bedrock topography, spatio-temporal structure of the ice-flow field) and not by the specific details of the age– depth distribution or temporal variations of the accumulation rate as produced by the ice-flow and mass-balance models themselves. The method of tuning the temperature–accumulation-rate relation further assumes that the present-day spatial pattern of accumulation rate is preserved in time. The dating of the internal ice layers, on the other hand, is taken directly from the Dome Fuji and eastern Dronning Maud Land ice cores bordering the RES profile under consideration at either end.
2. The Res Profile between Kohnen and Dome Fuji
2.1. Field data
Since the mid-1990s the Alfred-Wegener-Institut (AWI) has flown numerous airborne RES profiles in eastern Dronning Maud Land, over an area covering approximately 106km2 (Reference Watanabe, Jouzel, Johnsen, Parrenin, Shoji and YoshidaSteinhage and others, 2001). The survey produced detailed maps of the subglacial topography at spatial resolutions down to a few km. Of interest in this study is the continuous RES profile of ~1300km length connecting the ice-core drilling sites of Dome Fuji and Kohnen (Fig. 1). The Kohnen ice core was retrieved under the European Project for Ice Coring in Antarctica (EPICA), and is also known as ‘EDML’ (EPICA Dronning Maud Land). The profile between the two ice-core sites approximately follows the crest-like ice divide connecting a point located nearly 25km west of Kohnen station with Dome Fuji station which is located several km from the apparent dome. Along this profile, internal layers can be traced with confidence to about half the ice thickness, as can the bed topography with only minor interruptions. Moreover, these internal layers can be independently dated from the ice-core chronologies at either end of the profile.
Since the RES profile follows the ice flank, horizontal flow velocities are significant and moreover exhibit a complex pattern dictated by the bed topography and the ice-sheet geometry in this part of Dronning Maud Land. The magnitude of the surface velocity as predicted by the ice-dynamic model in the nested domain (see further below) is shown in Figure 2. Along the profile, these velocities vary from <0.1ma–1 close to Dome Fuji to >2.5ma–1 a few hundred km upstream from Kohnen. At Kohnen itself, the model predicts a surface velocity of 0.73 ma–1, close to the observed surface velocity of 0.76 ma–1 (Reference Wesche, Eisen, Oerter, Schulte and SteinhageWesche and others, 2007). Likewise, horizontal advection of deeper ice layers in at least the first 400 km of the RES profile upstream from Kohnen should be taken into account to be able to relate the information from these deeper layers to present-day surface characteristics.
2.2. Data preparation
Preparation of the observational data is needed to minimize the influence of unwanted layer features. Original RES measurements contain features that are not relevant to past accumulation rates. In the first place, vertical positions of the layers are contaminated with noise caused by various instrumental errors of measurement. Noisy data make it difficult to trace a particular internal layer (Reference Oerter, Graf, Wilhelms, Minikin and MillerLeysinger Vieli and others, 2004). Secondly, the layering has a wave-like structure due to the ice flow over peaks and depressions of the bedrock. Reference GudmundssonHindmarsh and others (2006) showed that the isochrones follow bed topography (drape) at long wavelengths of bed undulations and tend to override at short wavelength (less than or equal to one ice thickness), while single mountains deform layers whatever their spatial scale. They also showed that overriding can also occur over bedrock undulations exceeding one ice thickness when the undulation of the topographic feature orthogonal to the ice flow is smaller than or equal to one ice thickness. Besides, the translation of the basal signal to the surface is wavelength-dependent (Reference FujitaGudmundsson, 2003). Clearly, a rather complex topography along our RES profile must have caused many features in the ice-layer stratigraphy that are likely to obscure the use of the layers for our purpose of deriving accumulation rates.
In our approach we scaled the absolute depths of the layers to the ice thickness to obtain relative depths varying from 1 at the bottom of the ice sheet to 0 at the surface. In this way the direct effect of bedrock undulations on the layer geometry is largely, though not completely, removed from the layers. From theoretical considerations, it is known that the extremes of the waves of the internal ice layers are phase-shifted with respect to extremes of the bedrock relief. Such a phase shift depends on the the layer depth and the undulation wavelength (Reference FujitaGudmundsson, 2003), so a complex pattern of phase shifts is expected. For instance, bedrock undulations with a wavelength more than ~2.5 km will influence both the vertical and horizontal position of the layer. The minimal horizontal shift will therefore be of ~600m (0.25 of the minimal wavelength in topography overridden by the ice flow).
Secondly, we had to fill the gaps in the bedrock profile at sections along the RES profile where internal layers were visible but a clear bedrock echo was missing. These missing sections were filled by projecting the deepest internal layer according to its relative depth in an iterative way. While prone to a degree of circular reasoning, this approach is certainly superior to any other linear or non-linear interpolation that could have been attempted.
Finally the bedrock profile was smoothed with a 20 km running average. This distance appeared to give the best results and is moreover in accord with the 20 km resolution of the surface topography used to calculate surface slopes in the ice-flow calculations. In this way we smoothed out all basal signals less than 20km wavelength. Light smoothing (as a running average over a distance of ~1 km, i.e. smaller than the Fine Scale Model resolution) was also performed for the individual layers, mainly to remove the observational noise. The final result after normalizing for ice thickness is displayed in Figure 3b. It shows a layer architecture smoothly varying with respect to relative depth. The remaining features are now believed to represent mainly surface accumulation patterns rather than bedrock topography patterns (although there is some relation between them, as surface depressions and ridges that influence accumulation patterns may be caused by bedrock topography).
2.3. Dating the internal layers
To date the internal layers, we simply took the average of the independent datings from the ice cores drilled at each end of our profile. We used the chronology from EPICA Community Members (2006) for Kohnen, and from Reference WilhelmsWatanabe and others (2003) for Dome Fuji. In general, the dating of the same layers did not differ by more than a few hundred years at either end, except for layer 6 (Table 1). We therefore suspect this layer to have discontinuities, especially since it is located in the brittle ice zone where interpretations may be unreliable. To exclude the unclear segment between x = 560 and x = 800 km along the RES profile, layer 6 was split into two separate layers and dated at one side only. Layer 4 also has numerous interruptions for x > 560 km and therefore is interpreted only for the segment between 0 and 560 km.
3. Calculation of the Accumulation Rates
The core of our approach consists in the calculation of the thinning function along the RES profile with an ice-dynamic model. Once this function has been established, the original average accumulation rate of a specific layer of certain thickness can be found by integration provided the chronology at the layer interfaces is known. The following relations hold between annual layer thickness λ(z), ice age A(z), accumulation rate P A(z) and thinning function R(z) (e.g. Reference SiegertReeh, 1989; Reference RotschkyParrenin and others, 2004): where λ 0 = P A(t) is the initial annual layer thickness at the surface (or accumulation rate) we want to reconstruct for a certain time interval. Here, λ 0 and P A(t) are expressed in ice-equivalent units, whereas λ(z) is expressed in actual thickness of the annual layer. Hence, R(z) represents both firn compaction and dynamic strain associated with ice flow and therefore will be larger than unity in the firn compaction layer. This unconventional representation of the thinning function avoids the need to artificially reduce the ice thickness and enables the layer depths to be taken directly from the RES records.
To reconstruct the thinning function, we use the relations given above, but the ice age and time-dependent accumulation history are now provided from an internally consistent calculation within the ice-dynamic model itself. The background of this approach is that the thinning function is determined mainly by ice dynamics and the geometry of the ice flow and is much less sensitive to the details of the temporal variations of the accumulation rate and the exact modelled chronology linked therewith, as was, for example, demonstrated by Reference RotschkyParrenin and others (2004). This enables us to use the modelled thinning function to infer the measured distance between two isochrones into an average accumulation rate for any period in the past. As it turns out, the chronology of the internal layers as simulated by the ice-dynamic model is also very close (to within several thousand years or better) to the inferred dating based on both ice cores (Reference HuybrechtsHuybrechts and others, 2007). We consider this to be a good validation of the modelled ice-flow history in the region of interest because the effect of horizontal and vertical advection of individual layer particles has been captured well. We use this feature to relate accumulation rates inferred for different time periods to the same location on the surface of the ice sheet, thereby correcting for the effect of horizontal advection of the deeper layers with respect to shallower layers.
The nested ice-dynamic model consists of two components which respectively provide the large-scale evolution of the whole Antarctic ice sheet and the details of the local flow history within the nested domain (Reference HuybrechtsHuybrechts and others, 2007). The Large Scale Model (LSM) is a comprehensive three-dimensional (3-D) thermomechanical model run on a 20 km grid with 30 layers in the vertical (Reference HamiltonHuybrechts, 2002). It provides boundary conditions for the Fine Scale Model (FSM), which considers higher-order stress gradients in the force balance and is implemented over þan area of 960 ±640 km in eastern Dronning Maud Land at a horizontal resolution of 2.5 km and 101 equally spaced layers in the vertical. All prognostic calculations take place in the LSM. The main output of the FSM is the diagnostic 3-D velocity field in the nested domain as a function of time. The initial field of accumulation rate used in the nested model is based on the Reference FrezzottiGiovinetto and Zwally (2000) compilation used in Reference Hindmarsh, G.J.M., Raymond and GudmundssonHuybrechts and others (2000), and updated with data obtained from shallow ice cores during the EPICA presite surveys (Reference Parrenin, Rémy, Ritz, Siegert and JouzelOerter and others, 1999, Reference Reeh, Oeschger, Langway and Jr2000; Reference Spikes, Hamilton, Arcone, Kaspari and MayewskiRotschky and others, 2007).
The procedure to derive the thinning function R(z) along the RES profile is identical to that described in Reference HuybrechtsHuybrechts and others (2007), to which the interested reader is referred for more details. In short, the regional ice-flow history is obtained from a forward experiment with the nested model over several glacial cycles forced by climatic data from ice-core and oceanic sediment records. Subsequently the 3-D velocity field obtained by the FSM is used in a Lagrangian back-tracing routine to reconstruct particle trajectories along
the RES profile to their place and time of deposition. Initially, ice particles were placed in the vertical at every 0.1% of ice-equivalent thickness in locations spaced 25 km apart along the RES profile in the region with relatively high horizontal flow velocities (0–350 km) and for every 50 km further along the RES profile. The procedure directly yields the depth–age distribution, surface conditions at particle origin, and a suite of relevant parameters such as initial layer thickness. The strain-thinning function is straightforwardly derived from the latter parameter. Finally, a correction is made for firn compaction in the upper 180m by fitting a relation for firn density ±(z) to the observed density profile of the B32 core (H. Oerter and F. Wilhelms, http://doi.pangaea.de/10.1594/PANGAEA.58815) obtained by dielectric profiling (Wil-helms, 2005): where ± 0 = 910 kgm–3 is the ice density, ± sur = 350 kgm–3 is the snow density at the surface, a = –0.0212m–1 and ρ c = 12.329 kgm–3 is a fitting constant to obtain the constant value for ice density at 180m real depth. The depth z is taken positive downwards.
Examples of thinning functions obtained in this way for locations spaced 100 km apart along the Kohnen–Dome Fuji RES profile are shown in Figure 4.
4. Inferred Accumulation Rates
The average palaeo-accumulation rate P A ð z 1, z 2 Þin each individual layer along our RES profile is found from numerical integration of the following expression:
which can directly be deduced from Equation (1) assuming a linear age profile within a layer. The thinning function R(z ') between the layer boundaries z 1 and z 2 comes from the ice-dynamic model as described above. A 1 and A 2 are the ages at the respective layer boundaries taken from the established ice-core chronologies.
The inferred average accumulation rates for the different internal layers are shown in Figure 5. The most conspicuous feature is the apparent difference between ‘Holocene’ (0–13.1 ka BP) and ‘glacial’ (14.7/15.8–72.7 ka BP) accumulation rates, the latter being a factor 1.5–2 times lower. Also the spatial variability of the interglacial accumulation rates is higher than that of the glacial accumulation rates by about the same factor, i.e. the relative spatial variability for these two groups is similar. In between these two groups of curves is an intermediate accumulation-rate reconstruction characterized by the highest spatial variability. This corresponds to the time period 13.1–14.7 ka BP (blue curve) encompassing the Antarctic Cold Reversal. Particularly remarkable are the inferred accumulation rates close to Kohnen station, which are lower than the interglacial accumulation rates but with much higher spatial variability, whereas further east along the RES profile the inferred accumulation rates take on characteristics more similar to the interglacial accumulation rates.
The inferred accumulation rates obtained here contain spatial variability features on a scale of 20–40km that appear, moreover, quite constant in time. Such a pattern of spatial variability is confirmed by many studies (e.g. Reference Giovinetto and ZwallyHamilton, 2004; Reference Eisen, Rack, Nixdorf and WilhelmsEisen and others, 2005; Reference LoriusKarlöf and others, 2005; Reference Spikes, Hamilton, Arcone, Kaspari and MayewskiRotschky and others, 2007) and is often attributed to systematic snow redistribution in accordance with surface micro-relief variability caused in turn by the transfer of bedrock undulations to the surface of the ice sheet.
In Figure 6 we compare the average accumulation rate (in ice equivalent) reconstructed for the period 0–5 ka BP with published large-scale accumulation maps of Antarctica (Reference Wesche, Eisen, Oerter, Schulte and SteinhageVaughan and others, 1999; Reference Hindmarsh, G.J.M., Raymond and GudmundssonHuybrechts and others, 2000; Reference Arthern, Winebrenner and VaughanArthern and others, 2006). Assuming that the late-Holocene accumulation rate is close to the present-day distribution, we find a good fit with the accumulation values established for both the Dome Fuji and Kohnen ice-core locations (inferred values of respectively 2.8 and 6.8 cma–1 compared to the observed values of 2.7 cma–1 (Reference WilhelmsWatanabe and others, 2003) and 7.0 cm a–1 (Reference Reeh, Oeschger, Langway and JrOerter and others, 2000)), but substantially lower accumulation rates between x = 600 and x = 1200 km in an area where direct measurements are lacking. In particular, the Reference Wesche, Eisen, Oerter, Schulte and SteinhageVaughan and others (1999) distribution systematically overestimates our reconstruction by 30–70% all along the RES profile. Also the inferred accumulation rates obtained in this study exhibit a much larger spatial variability that is absent in the large-scale datasets because of the smoothing and interpolation procedures employed in those datasets.
5. Improving the Temperature–Accumulation-Rate Relation
The relation between changes in accumulation rate P A and changes in mean annual air temperature above the surface inversion layer T I (in K) is widely used in ice-sheet modelling to reconstruct palaeo-accumulation rates. It owes its success to the fact that changes in the accumulation rate can be estimated directly from temperature variations obtained from ice cores, without the need to invoke more complicated physical processes. The relation is usually expressed as a modification of the equation based on the original approach suggested by Reference OerterLorius and others (1985): where P A[T IðtÞα is the local accumulation rate at time t, P A[TIð0Þα is the present-day accumulation rate at the same place, T 0 = 273.15 K, T I(t) = 0:67T S(t) þ 88:9 (in K) is the temperature at the top of the inversion layer at time t, T I(0) is the same inversion temperature at the present day, and T S(t) is the surface temperature at time t. The factor β[T I(t) ρ TI(0)] takes into account glacial–interglacial changes of the accumulation rate that are not explained by the original formula based on the Clausius–Clapeyron relation . In the surroundings of Dome C ρ = 0.035 ±0.012 (EPICA Community Members, 2004). To fine-tune this relation for eastern Dronning Maud Land and establish a value for ±, we assume that P A[TI(0)] can be equated with the average late-Holocene accumulation rate inferred from the upper layer of the ice sheet (0–5 ka BP).
To estimate β in Equation (5), we need to evaluate the ratio between the average accumulation rate during a particular time period and the average present-day accumulation rate at the position where the snow was originally deposited. An approximate way to account for the effect of horizontal advection of the ice, i.e. for the effect that accumulation patterns in deeper layers have been transported downstream with respect to their time of deposition, is to introduce a scaling factor C for each internal layer point in the nested model. This is defined as the present-day ratio between surface accumulation rates at the current horizontal location xi of the layer point and the upstream location xj where the layer point was originally deposited: We thenmultiply P A[TI (t)] xi by the factor C to be able to relate it directly to P A[T I(0)] xi at the same location along the RES profile. This assumes that the spatial pattern of accumulation is stable in time at distances comparable to the maximum advection distance (~100 km; see Figs 2 and 5), which seems reasonable. Note, however, that the back-tracing algorithm directly provides the horizontal distance over which deeper layers have been advected, without the need to introduce the factor C. However, since this operation is only performed at discrete locations spaced 25–50km apart along the radar profile, the use of an interpolated scaling factor C enables us to account for horizontal advection for every internal layer point of the profile in a continuous way.
Figure 7 shows the ratios between the derived palaeo-accumulation rates and the average accumulation rate for the period 0–5 ka BP including the correction for the effect of horizontal advection. To estimate ±, we use these ratios and the ratios of present-day and palaeo-temperatures at location xi (Fig. 7b) reconstructed from using the d18O record obtained from the Kohnen ice core (EPICA Community Members, 2006) corrected for surface elevation changes simulated by the flow model.
Since the multiplier β is introduced in Equation (5) to scale the glacial–interglacial magnitude of temperature, it is reasonable to carry out its estimations only for the glacial (14.9/15.8–73.3 ka BP) period using the period 0–5 ka BP as the interglacial reference. Because of the discontinuity in the original layer 6 (see Table 1), we estimated the average β separately for two segments of the RES profile, namely between 10 and 560 km and between 800 and 1280 km. For every point along the RES profile within these two spatial intervals, β is calculated as a weighted average (i.e. proportional to the duration of the time period between two layer boundaries). A moving average was then applied to β over an interval of 20 km and its standard deviation was calculated.
The result is shown in Figure 8. For the segment between 10 and 560 km, closest to Kohnen station, we find an average β = 0.046l0.024, and for the segment between 800 and 1280 km, closest to Dome Fuji, β = 0.052l0.010. The large variation in the value of β is due to the fact that the ratios of palaeo-accumulation rates include a component from spatial variability that obviously cannot be explained by a corresponding spatial variability of temperature. It may also indicate that introducing a linear β-correction, while being the simplest thing to do, may not capture well all of the processes not described by the dependence of precipitation on the saturation vapour pressure at the top of the surface inversion layer. It therefore makes sense to consider the value of β over a much larger area than the scale of its spatial variability. Clearly, the average values of β around Kohnen are higher than the estimate of β = 0.035l0.012 derived for Dome C (EPICA Community Members, 2004), so the glacial–interglacial contrast is larger than that derived from the temperature-dependent change of the saturation vapour pressure alone. The reason for this difference is not immediately clear. One possible explanation could be a relatively higher influence of synoptical events in the present-day accumulation rate regime in Dronning Maud Land as compared to the surroundings of Dome C, and perhaps a more abrupt reduction of this influence during the last glacial period, but this remains very speculative.
6. Concluding Remarks
The model experiments discussed in this paper showed how internal reflection horizons can be used to infer current and past accumulation rates over a large distance and how detailed flow modelling can help to estimate the thinning functions required by our approach. The field data are often noisy and contaminated with information that is not relevant to extraction of accumulation rates. Smoothing and averaging procedures are required to prepare datasets and make the results amenable for interpretation. Nevertheless, we have shown that a consistent picture can be distilled. Our main finding is that accumulation rates in central Dronning Maud Land, where surface control points are lacking, are substantially lower than shown in available large-scale gridded datasets. Moreover, the spatial variability of glacial accumulation rates appears much lower than that inferred for the late Holocene. Using inferred palaeo-accumulation rates along our RES profile, we were able to calculate the enhancement factor for the glacial–interglacial contrast in accumulation rates. In the area surrounding Kohnen station, glacial accumulation rates were found to be 25–35% lower than predictions based on the gradient of the saturation vapour pressure above the surface inversion layer.
Acknowledgements
We greatly appreciate suggestions and detailed comments by E. Waddington, G. Leysinger Vieli, F. Parrenin and an anonymous reviewer, which led us to significantly improve the manuscript. Part of the work was supported by the Belgian Research Programmes on the Antarctic and on Science for a Sustainable Development (Belgian Federal Science Policy Office) under contracts EV/11/08A (AMICS) and SD/CA/02A (ASPI). This work is a contribution to the European Project for Ice Coring in Antarctica (EPICA), a joint European Science Foundation/European Commission scientific programme, funded by the European Union (EPICA-MIS) and by national contributions from Belgium, Denmark, France, Germany, Italy, the Netherlands, Norway, Sweden, Switzerland and the United Kingdom. The main logistic support was provided by Institut Polaire Franc¸ais– Emile Victor (IPEV) and Programma Nazionale di Ricerche in Antartide (PNRA) (at Dome C) and the AWI (at Dronning Maud Land). This is EPICA publication No. 215.