2.1 Technical background and functional principle of GPR

The principle of GPR consists of measuring the two-way travel times of parts of an electromagnetic pulse emitted toward the ground that are reflected on discontinuities within the observed medium. When pulses are emitted at points along a profile, lateral continuity of the physico-chemical transitions within the medium yields wiggles of similar travel time (hence depth) along the resulting radar trace. Once put side by side, these wiggles are expressed in the form of horizons, also called IRHs, visible on the recorded radargram. In this respect, GPR is analogous to seismic reflection but with an electromagnetic wave instead of an acoustic wave.

Figure 2Continuous post-processed radargram corresponding to the merged profiles as outlined in Fig. 1. The numerous observable IRHs result from as many reflections within the snowpack and are each characterized by a laterally varying depth and a specific age. The coloured (blue, green and pink) aligned dots emphasize the three most evident and continuous IRHs that are used in the present study. The radargram as depicted here (restricted to the first 35 m of firn) has undergone a time migration (see Sect. 2.3), meaning that the vertical axis represents true depths instead of two-way travel times. The apparently shifted vertical scale results from the time-zero correction. The horizontal axis is the cumulated distance along the profile as given by the GPS data. The white vertical stripe is due to a breakdown of the GPR leading to a gap of less than 2 km in the data. Also represented are the crossing points with the ice/firn cores DC, S1, S2 and S0 where the intersecting depths with the three selected IRHs have been inferred to check their isochronicity (see Sect. 4). The red and green upper banner refers to the position with respect to the topographic divide and correspond to the colour convention of the bottom right of Fig. 1 (red for “divide” and green for “off-divide” locations).

Download

Discontinuities responsible for these partial reflections result in changes in the complex dielectric constant ϵ^* of the medium. For the snowpack, essentially two main factors can be responsible for these changes as explained in Eisen et al. (2008) and references therein. The first concerns the real part, namely the permittivity, which varies with density, whereas the second acts on the imaginary part, which is sensitive to the conductivity of the medium. As a result, radar horizons within the first few hundred metres of firn and ice are essentially due to changes in density and acidity from volcanic deposits.

The huge advantage of the IRHs comes from their isochronous nature, although this has not been rigorously (physically) demonstrated, only heuristically by connecting ice core chemistry analysis of different drill sites (Eisen et al., 2004; Frezzotti et al., 2005) or by comparing accumulation rates obtained from both IRH depths and stake lines (Vaughan et al., 2004). This isochronicity is however understandable as it results from processes operating at the surface at the same time, such as acidic atmospheric deposits or density changes from the metamorphism of snow (at least the one initiated by surface interaction with the atmosphere like radiation crusts and wind redistribution). However, the physical processes leading to these IRHs still remain poorly known. For instance, as pointed out by Eisen et al. (2008), annual layers are still visible with radar wavelengths much larger than the annual thickness (a typical 100 MHz wave in the cold dry snow of the plateau has a wavelength of 2 m: much larger than the 10 cm annual snow thickness at Dome Concordia (DC) for instance), which contradicts the theory unless constructive interferences are considered as suggested by Palli et al. (2002) for instance.

Assuming isochronicity, the spatial variations of snow accumulation can be obtained from the varying depth of the IRHs. It thus implies transforming the two-way travel times of the reflected waves into actual depths. This is only possible if the wave velocity within the medium is known all the way from the surface down to the specific IRH and the time radargram is properly migrated (see Sect. 2.3). Inferring accumulation rates then requires the association of an age to the IRHs from independent means such as correct recognition of fallouts from identified volcanic and/or bomb testing events along firn/ice cores for instance (see Sect. 4).

2.2 GPR set-up

For the present project, GPR data were acquired with a MALÅ^® ProEx GPR equipment fitted with a 100 MHz “rough-terrain antenna” towed behind one of the tractors used during the traverse with a constant separation of 2.2 m between the emitter and receiver parts. The apparatus, operating frequency and set-up are very similar to those used on a similar traverse between base stations of Dumont d'Urville and DC during the 2008/2009 austral summer (Verfaillie et al., 2012). The triggering was set to 1 s, meaning a radar trace every 4 m or so given an average speed of 14 km h⁻¹ for the convoy. The time window was set to 1 µs (allowing for an investigation depth of some 100 m) and was sampled at 1.1 MHz, leading to roughly 1100 samples per trace. During acquisition a first improvement of the signal-to-noise ratio was obtained by an up to 64-fold stacking (which, given the time window of 1 µs, remained compatible with acquisitions made every second). Moreover, the system was connected to a GPS receiver mounted on the vehicle that recorded the geographic position of every single trace along the profiles.

2.3 Resulting GPR radargrams

Profiles outlined in Fig. 1 represent 5 days of measurements between 13 and 17 January 2012 and led to an almost continuous record of 630 km. The corresponding post-processed and concatenated radargram is shown in Fig. 2. The post-processing sequence consisted of time-zero corrections, a zero-phase low-cut filter (devow) to remove direct continuous currents and an “energy decay” gain to compensate for the volumetric spreading signal attenuation. Band pass filtering was reduced to a low-pass filter in the form of a spatial averaging over 50 traces which, given the huge number of traces (164 000 altogether), appeared more efficient than traditional finite impulse response (FIR) filters. Time-to-depth conversion was obtained by migrating the radargrams with the help of a vertical velocity profile for the radar wave. By investigating a similar region with similar firn properties, we followed the approach of Verfaillie et al. (2012), which consisted of only considering the effect of the firn density on the electromagnetic wave velocity as a result of dry and clean snow over the Antarctic Plateau. From the empirical relation of Kovacs et al. (1995) relating permittivity and density, the wave velocity c is obtained according to

with c_v being the wave velocity in vacuum and ρ the firn density (relative to water density). The density profile comes from a recent core drilled at DC (Leduc-Leballeur et al., 2015) where high-resolution measurements allowed for a detailed profile (see Sect. 5.2). The true restitution of dipping reflectors theoretically requires a topographic migration, but the latter appears only necessary where the topographic variations are of the same order as the reflector depths. Given the extreme flatness of the investigated area of less than 2 m km⁻¹ (topography of the upper surface and hence that of the IRHs), no such correction has been applied here. Last, ice sheet dynamics potentially changing the geometry of the IRHs has to be accounted for, except in the upper part of the ice/firn column, especially over areas of slow motion such as over the Antarctic Plateau as in the present case (Eisen et al., 2008).

3.1 Identification of volcanic eruption signals and volcanic chronology

Although sulfate in Antarctic snow comes from sea-salt spray and to a lesser degree from crustal erosion (Cole-Dai et al., 2000; Maupetit and Delmas, 1992), and principally from atmospheric oxidation of biogenic dimethylsulfide (DMS) emitted by oceanic phytoplanktonic activity (Prospero et al., 1991; Saltzman, 1995), volcanic eruptions are also major sources of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ during active eruptions. To identify volcanic signals in ice cores, it is thus necessary to calculate the non-sea-salt sulfate, (nss${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) corresponding to total sulfate minus sea-salt sulfate, and set a threshold above which spikes can be attributed to volcanic deposition. Over the East Antarctic Plateau areas, part of the sulfate background that could be related to continuous emissions from non-explosive volcanic activity seems to be minor (Castellano et al., 2004; Cole-Dai et al., 2000; Patris et al., 2000). The sea-salt sulfate contribution to total sulfate budget, evaluated using Na⁺ as a specific marker (Palmer et al., 2001; Röthlisberger, 2002), is less than 10 %. The Holocene crustal contribution, calculated by non-sea-salt Ca²⁺ as a continental dust marker (Castellano et al., 2004; Plummer et al., 2012; Röthlisberger, 2002), is even lower (< 0.05 %). Since these contributions are of the same order of measurement reproducibility, we did not correct sulfate concentrations and will not distinguish between total sulfate and non-sea-salt sulfate in the following discussion. However, particular attention was paid to developing a reliable method for distinguishing volcanic signals from non-volcanic sulfate background, and an outline of the method used is described in Cole-Dai et al. (1997, 2000), Stenni (2002), Castellano et al. (2004) and Igarashi et al. (2011). The procedure employed to identify candidate signals of volcanic eruptions in the sulfate profiles of the different cores follows that of Gautier et al. (2016).

Between five and eight volcanic events were identified in each TASTE-IDEA firn core, and 10 events were detected in the VOLSOL-1 ice core from the surface to 20 and 30 m through the application of the above-described procedures (see summary in Table 1). In the considered depths for shallow TASTE-IDEA cores and the VOLSOL-1 ice core, the main volcanic events in the 1601–2012 CE period are Pinatubo (Philippines, 1992 CE), Agung (Indonesia, 1964 CE), Krakatau (Indonesia, 1885 CE), Coseguina (Nicaragua, 1836 CE), Tambora (Indonesia, 1816 CE), unknown (1809 CE), Jorullo-Taal (Philippines, 1758 CE), Serua (Indonesia, 1695 CE), Gomkonora (Indonesia, 1676 CE) and Huyanaputina (Peru, 1601 CE). The sulfate deposition of the Kuwae eruption (Vanuatu, 1457–1458 CE; Salzer and Hughes, 2007; Sigl et al., 2013, 2014, 2015) can now be dated as 1459 CE, i.e. 5 to 6 years later than previously assumed (Gao et al., 2006, 2008).

Table 1Depths of the dating events along the mentioned cores from known volcanoes as well as thermonuclear bomb testing.

^a The value in parentheses represents the depth shift to apply as the core was extracted a year earlier than radar measurements (see Sect. 4.2). ^b Time markers of the S0 and Explore cores have been merged to produce a hybrid depth relationship (see Appendix C). ^c A hyphen indicates an undetected volcano or an unmeasured nuclear fallout. ^d These two bomb test fallouts were actually measured on a nearby shallow core close to the VOLSOL-1 one in 2011/12 and served in the determination of its age–depth relationship. ^e Not accounted for in the resulting age curve since it was not clear to which volcano this sulfate peak should be attributed. ^f These two peaks, dated 1453 and 1460 CE, have recently been measured in the Explore core and are attributed to the Kuwae volcano (Robert Mulvaney, personal communication, 2017).

Download Print Version | Download XLSX

The temporal duration of volcanic signals in the ice has been evaluated by several studies (Castellano et al., 2004; Cole-Dai et al., 2000; Palmer et al., 2001), and general temporal durations range from 1 to 3 years. In this work, the temporal duration of volcanic signals lagged between 1.2 (3 consecutive samples) and 6 years (15 consecutive samples).

3.2 Artificial radionuclides deposition over Antarctica and associated chronology

Artificial radioisotopes resulting from atmospheric thermonuclear tests carried out between 1953 and 1980 were deposited in Antarctica after transport in the upper atmosphere and stratosphere, creating distinct radioactive reference levels in the snow. The dates of arrival and deposition in this polar region are well known and therefore provide a means to estimate Antarctic snow accumulation rates or describe air mass circulation patterns (Magand, 2009; Pourchet et al., 2003, and references therein). The 1955 and 1965 CE radioactivity peaks provide two very convenient horizons for dating snow and firn layers and thus measuring accumulation. Special techniques have been developed over the last 40 years to detect and measure artificial and natural radionuclides present in the ice sheets (Magand, 2009; Pourchet et al., 2003, and references therein). In Antarctica, ⁹⁰Sr, ²⁴¹Pu (deduced from ²⁴¹Am analysis) and ¹³⁷Cs radionuclides constitute the well-known debris products of atmospheric thermonuclear tests between 1953 and 1980 CE that we could still use to identify the two well-known reference layers in snow (1955 and 1965 CE peaks), corresponding to the arrival and deposition of artificial radionuclides in this region. Total beta counting and gamma spectrometry remain the most frequent radioactivity measurement device systems used to clearly detect artificial radionuclides and unambiguously determine the 1955 and 1965 CE peaks depths in the cores. The 1955±1 and 1965±1 CE peaks were both identified in each TASTE-IDEA firn core and in an extra-shallow core close to the VOLSOL-1 one (see Table 1).

4.1 Identification of major IRHs

Interpretation of radar data first consisted of identifying contrasting IRHs whose continuity could be tracked all along the entire merged profile. Three equally spaced IRHs were selected within the first 35 m of firn (no sufficiently clear reflectors could be used any deeper) and are emphasized with associated colours in Fig. 2. In order to properly capture the spatial variability, selection of these reflectors was done every kilometre or so, yielding each time a depth and the corresponding geographical position depicted as a coloured dot in the figure (the resulting number of dots leading to an almost continuous horizon). Of interest is the phasing of these three IRHs, which is characteristic of a stationary accumulation spatial pattern. It comes from the fact that local extrema in accumulation keep the same positions through time, allowing for cumulative effects with time which amplify IRH undulations with depth (Verfaillie et al., 2012).

4.2 Methodology for dating the IRHs

Dating is achieved by detecting the depths at which the reflectors intercept (or pass very close to) an ice/firn core where a depth–age relationship has been obtained (see Sect. 3). In the present case, as can be seen from Figs. 1 and 2, each of the three reflectors passes over several ice core locations and sometimes several times over the same one, thereby offering redundant dating possibilities. Due to its reliable dating curve including time–depth markers regularly spaced down to about 30 m, the DC-VOLSOL-1 core was selected as the reference and served for proposing an age ab initio for each of the three IRHs. The choice for the DC core is also justified by the numerous available ice core dating schemes for the site which gives a good level of confidence in the proposed age–depth relationships (Parrenin et al., 2007).

As for the other cores of the project, the dating is less reliable as the result of shallower depths, making the resulting fits questionable for the 2 deepest IRHs. A longer record has been artificially derived from the shallow S0 core supplemented by the nearby (70 m) 110 m deep Explore core from the companion Explore programme. Because of ambiguity with a sulfate peak at 10.48 m, only the very deepest points of the Explore core have been used to supplement the S0 age–depth relationship (see Fig. C1) used for estimating the age dispersion as described in Appendix C.

Having been extracted during the VOLSOL programme (2010/2011 field season), the DC-VOLSOL-1 core is a year older than the radar measurements. As a consequence, a depth correction has been applied to each IRH volcanic depth in the form of a downward shift of 8 cm (modern snow-equivalent accumulation for DC) weighted by the density ratio between the relevant depth and the surface (i.e. 8 cm for the surface, reducing to 4.53 cm for the Kuwae volcano at the depth of 31.20 m; see Table 1).

Figure 3Identified and dated deposits placed in an age–depth representation according to their depths within the DC-VOLSOL-1 ice core. Volcanic horizons were completed here by the integration of the two 1955 and 1965 thermonuclear horizons (depths of 4.20 and 4.60 m) from dedicated measurements along an extra core drilled at DC in 2012 during the scientific traverse. A degree-2 polynomial fit (black line, r²>0.99) was then used as a continuous dating curve permitting the dating of any intermediary depth.

Download

4.3 Age results and the isochronous character of the IRHs

As can be seen from Fig. 3, IRH depths of 15.08, 22.97 and 31.89 m at the VOLSOL-1 core location respectively give ages of 233.74, 390.54 and 591.71 years (later rounded to 234, 391 and 592 years).

The two main factors that can potentially alter the isochronous character of the selected IRHs are (i) errors in age–depth relationship at each coring site and (ii) errors in the depth estimates of the IRHs at these same sites deduced from the radargrams. As for the first, owing to the short-scale changes in snow deposition and erosion as expressed by the surface roughness, the archiving process can sometimes be significantly different for neighbouring cores, thereby raising the problem of the correct representativeness of a single core. For several volcanoes Gautier et al. (2016) found differences in the depths of the corresponding sulfate peaks by as much as ±20 cm from one core to the other among a set of five different cores all separated by 1 m. As a consequence, on top of errors in the correct depth measurements of the chemical peaks, surface processes lead to an uncertainty in the core-derived age–depth relationships. Figure B1 in Appendix B reports the time markers along the five cores from the study of Gautier et al. (2016). From quadratic fits along each core, five different ages are computed for our IRH depths of 15.08, 22.97 and 31.89 m, yielding respective root mean square (rms) deviations of 6, 7 and 7 years with respect to our reference ages of 234, 391 and 592 years. These deviations can be considered as good estimates of the uncertainty in the correct representativeness of the single core of our study. The residency time of volcanic aerosols in the stratosphere (from 2 to 4 years) is variable from one volcanic event to the other and therefore also contributes to some uncertainty in the depth–age relationship, which makes for an overall uncertainty of 8 years in the age for a given depth, a result very similar to that of Eisen et al. (2004).

As for the depth accuracy of the selected IRHs, it mainly results from the physics of radar and the chosen frequency for the electromagnetic wave. Eisen et al. (2004) gives a review of all possible resultant sources of errors, and their results are applicable to our case because of the use of the same commercial radar set from MALÅ^® Geoscience Sweden, with the only difference being that they use 200 and 250 MHz frequencies instead of our 100 MHz. One initial source of uncertainty lies in the thickness the 100 MHz wave can resolve. Although it theoretically amounts to λ∕4, one-fourth of the wavelength in the firn, λ∕2 is usually considered as more realistic, which in our present case gives 1 m. Additionally, the pulse width also matters (length of the energetic part of the source wavelet) and was determined by common midpoint analysis to be 12 ns with our system (Verfaillie et al., 2012), leading to a tracking accuracy of half this duration, which in the firn gives a shift of some 1.2 m. Since these two errors do not systematically add up, one generally considers the worse of the two. Finally, interpreting radargrams and picking reflectors are partly subjective and operator-dependent processes, which leads us to estimate the overall uncertainty of the IRH depth to be 1.5 m at the worst. Given a relatively constant slope of 20 yr m⁻¹ of the age–depth relationships (see Fig. 3), this depth uncertainty maps into an age uncertainty of about 30 years. This uncertainty in the depth positioning nevertheless has a strong systematic component (in other words the error in depth positioning remains mainly the same along a given IRH) and therefore only partly alters the isochronous character of the reflectors. The same can be said for depth errors resulting from a wrong assessment of the vertical velocity profile. Only lateral variations not properly accounted for would contribute, but it is well known that the plateau firn only undergoes minor lateral variations in its depth distribution of density (Muller et al., 2010) and hence velocity at the scale of a couple of hundred kilometres. As a result, one can reasonably argue that a true isochronous IRH can theoretically deviate by 10 to 15 years. Appendix C presents the different inferred ages when reporting our IRH depths along the other cores of the project (S1, S2 and S0–Explore). From the rms dispersions of respectively 3, 6 and 13 years for IRHs 1, 2 and 3 at the crossover points as revealed by Fig. C1 we come to the conclusion that our three IRHs can be considered as isochronous. We nevertheless have to keep in mind that the absolute error of the age (to be used in the budget error for mass balance computation; see Sect. 5.6) remains potentially higher than these 10 to 15 years by integrating the biases from systematic incorrect positioning and errors from time-to-depth conversion in the radargrams.

5.1 Accumulation computation

Assuming an IRH at a given depth H, the total mass M of the unit surface (1 m by 1 m) firn column above it can be expressed as

where M can be considered as a mass flux per square metre because of the implicit unit surface. It is therefore expressed in units of kilograms per square meter which, given the density of water, also represents millimetre water equivalent. The average accumulation rate over the period since the deposition of the IRH is finally obtained by dividing by the age of the IRH, leading to a final result in millimetre water equivalent per year.

5.2 Firn density profiles for accumulation rate computations

Similarly to the depth and age, proper estimation of the density profile down to a given IRH is crucial for a good assessment of the corresponding accumulation rate. If spatial variations of accumulation mainly result from differences in snow cumulative heights as revealed by undulating reflectors, one can reasonably wonder to what extent geographical changes in density may also have an impact. Unfortunately, over the Antarctic Plateau, density profiles are usually concentrated at limited sites with almost no reliable data in between. In the present case, density profiles from deep ice cores are only available at the two “extremities” of the radar line (DC and S0). Some extra-shallow cores were also drilled along the radar profile (S1, S2, S2B) but their maximum depths (about 20 m) are limiting at least for interpreting the two deepest IRHs. However, owing to the fact that over the plateau meteorological parameters controlling the densification process (e.g. air temperature, wind activity, solar radiation) exhibit small gradients at the regional scale, a relative uniformity in firn density profiles can be first assumed over distances of the order of 100 km as can be anticipated from the small divergence between the DC and S0 density profiles (see Fig. 4). The limited variability was confirmed by Fujita et al. (2011) from pits and surface snow measurements during the JASE traverse, and more specifically the resulting small impact on accumulation rates was demonstrated by Ruth et al. (2007) when comparing water equivalent depth profiles between Dome Fuji and EPICA (European Project for Ice Coring in Antarctica) DML (Dronning Maud Land). This led us in a first instance to consider only deep ice cores from the DC site along which precise and high-resolution density profiles have been performed (Leduc-Leballeur et al., 2015). A later sensitivity test is proposed (see Sect. 5.4) consisting of assessing changes in accumulation to expect from integrating the contribution of density profiles from different locations. These changes would then be compared to the uncertainties solely arising from the degree of representativeness (due to the small-scale variability) and measurements errors in density when using a single density depth profile for computing the accumulation rate all along the radar line (see Sect. 5.5).

Figure 4Density profiles (left scale) as mentioned in the text and detailed in Table 2. Also shown in the bottom part is the corresponding cumulative mass with depth (right scale). For the sake of clarity, only cores (and associated cumulative masses) exceeding 20 m are fully represented (DC5, DC3, Firetracc and ITASE-98 ice cores). Insets correspond to enlargements around the depths of IRHs 1, 2 and 3 where all mass curves are depicted and allow for an estimation of the dispersion in cumulative masses and in the resulting accumulation rates (to be used later in Sect. 5.5). The upper panel shows a comparison of the DC5 and S0 (Explore) density profiles with estimated measurement error bars that only served for the computation of the uncertainty in the wave velocity profile (see Table 3).

Download

Table 2Density cores at DC along with the characteristics of their quadratic fits. Cumulative masses (kg m⁻²) down to each of the three IRHs at respective depths of 15.08, 22.97 and 31.88 m are computed from the raw data as well as from the fitted curves.

^a In kg m⁻². ^b Leduc-Leballeur et al. (2015). ^c Depth not reached by the core.

Download Print Version | Download XLSX

5.3 Density profiles at DC

In the framework of the EPICA project (Augustin et al., 2004), along with the main core, several extra-shallow drilling projects have yielded as many density profiles over at least the first 50 m for the DC area (see Table 2 and Fig. 4). A first question arises as to which of these cores to use for computing accumulation rates; in other words, what can be the sensitivity of the results to the choice between different nearby cores at a given site like DC? The lower part of Fig. 4 shows the depth–density profiles for each of the cores listed in the table along with the corresponding cumulative mass with depth (right scale in the figure). The cumulative masses have been computed either from the raw data or from the quadratic fits (see corresponding results at the three IRH depths in the table). Despite a pronounced scatter in the individual data points around their respective fits (see Fig. 4 and associated root mean square error (RMSE) in the table), considering raw or fitted data does not significantly change the depth–cumulative mass with differences not exceeding 1.15 % at the most (ITASE-98 at the depth of IRH1). This illustrates a random scatter of data points above and below the fit, which leads to an almost exact compensation when summing up the cumulative mass. As a consequence, fitted density curves will hereafter be considered in the computations of accumulation rates. The second result comes from the limited dispersion in the cumulative mass from one core to the other as can be easily seen from the figure. The maximum discrepancy is to be found between the ITASE-98 and Firetracc ice cores with relative differences of 7, 5 and 4 % at the depths of respectively IRHs 1, 2 and 3. Of interest is the almost perfect overlap of the DC3 and DC5 curves, which come from the same project and have been measured according to the same protocol by the same operator. These results tend to show the following:

• a.

  Most of the difference between different cores results from measurement errors.

• b.

  The variability of the density expressed at the scale of the typical distance between these cores (from metre to kilometre scale) is of stochastic type (Libois et al., 2014) and is therefore rapidly cancelled out by the depth-averaging process of cumulative mass computation leading to accumulation rates. Indeed the DC3 and DC5 cores are separated by 1.5 km but still show remarkable consistency in their quadratic fits or in their associated cumulative mass distribution with depth.

• c.

  Provided that systematic error measurements can be minimized, a single density profile should be considered as representative of the local accumulation rate.

The very good similarity of the DC3 and DC5 cores attests to the quality of density measurements as confirmed by the novelty and strictness in the measurement protocols (Leduc-Leballeur et al., 2015). This led us to choose the DC5 as the reference for our accumulation rate computations.

5.4 Sensitivity of the accumulation rate to geographically different density profiles

Figure 5 shows the three accumulation rates derived along each of the three IRHs along the radar line from the density profile provided by the DC-5 reference core. Since the computation uses a single density profile, no lateral changes in density are accounted for, and the stationary accumulation pattern is therefore expressed as shifted IRHs of the same shape. Of importance is the pronounced gradient when following the ridge going from DC towards S2 (see map in Fig. 1) with a more than 20 % loss within 250 km.

In order to test the sensitivity to the use of density profiles originating from different distant places, accumulation rates are separately computed here from both the reference DC5 density profile and the Explore one (the latter being the only reliable and deep enough density profile at our disposal). The difference is expressed in the form of a uniform shift along a given IRH (as can be seen from the insets) but with a decreasing sensitivity with the depth of the IRH. The difference in accumulation rates from the use of the DC and S0 (Explore) densities amounts to 1 mm w.e. yr⁻¹ for the 234-year IRH1 (lower inset) and reduces to about half of it when considering the 592-year reflector (upper inset). The explanation is to be found in the rapidly converging fits of the DC and S0 density curves with depth (top of Fig. 4). Indeed, the increasing cumulative difference in mass is more than counterbalanced when dividing by the age of the IRH in the accumulation rate computation process. Using the less reliable S1, S2 and S2B (not represented) did not lead to a significantly higher sensitivity. As a consequence, a maximum deviation of the order of 1 mm w.e. yr⁻¹ is to be expected from accounting for the geographical changes in density along the profile.

Figure 5Computed accumulation rates since the deposition of the three IRHs according to the polynomial fit of the DC-5 density profile. The lower inset represents an enlargement over the dotted square and shows the sensitivity of the accumulation rate at the uppermost IRH to the use of the DC-5 and S0 (Explore) cores, whereas the upper one does the same but for the lowermost IRH. Error bars represent the estimated uncertainty resulting from the combined density error measurements and small-scale spatial variations. These error bars (see Sect. 5.5) are only depicted once since they are uniform all along the profile. Ranges for the accumulation rate at S2 inferred from the SR50 sensor (see Sect. 6.1) are represented by the black box along the S2 vertical line.

Download

Table 3Different source terms in the overall error budget for accumulation rate.

^a rms of $\mathrm{\Delta }{\dot{A}}_{\mathrm{1}}$ and $\mathrm{\Delta }{\dot{A}}_{\mathrm{2}}$. ^b rms of ΔT_c and ΔT_p. ^c Computed with averaged depths of respectively 13.87, 21.09 and 29.43 m and corresponding averaged densities of respectively 397, 498 and 553 kg m⁻³ for IRHs 1, 2 and 3. ^d rms of ΔH_r and ΔH_v.

Download Print Version | Download XLSX

5.5 Proper choice of the density core to use

Assessing the uncertainty in the accumulation rates solely arising from measurement errors in the density profile used is not straightforward since some of these errors are random and certainly compensate with depth, a least partially. Moreover, small-scale accumulation variability that influences the representativeness of a single core for a given location also potentially contributes to accumulation uncertainty, even if limited in amplitude as was shown in Sect. 5.3. Because they result from different projects, implying different measurements protocols, the numerous cores drilled in the DC area provide a means of estimating the order of magnitude of the combined effect of these two terms. In particular, the dispersion between all these cores should implicitly include most of the measurement errors and would only miss a systematic component (i.e. an identical error for all measurements), which in any case must remain small. Computation of the cumulative mass down to IRH1 for these eight cores listed in Table 2 gives an average of 6127 kg m⁻³ with a standard deviation of 141 kg m⁻³. When dividing by the age of IRH1 (234 years), one obtains 26.18±0.60 kg m⁻³ yr⁻¹ or 26.18±0.60 mm w.e. yr⁻¹. Similar computations down to the two deepest IRHs (but only concerning the first four cores) respectively give 25.75±0.54 and 25.49±0.40 mm w.e. yr⁻¹ for IRHs 2 and 3.

The resulting dispersion is then considered as representative of the total error due to the use of a local density curve and has been represented in the form of the error bars in the insets of Fig. 5. The same error bars have been applied to the accumulation profiles from the Explore core at S0, but they probably represent an underestimation due to a less stringent density protocol in the density measurement. From the systematic overlap and the fact that the Explore error bars represent lower bounds, we come to the conclusion that accounting for the available geographical spatial variability of density does not cause significant changes in the accumulation rates. In fact, properly accounting for the geographical distribution of density would normally require integration of all available density profiles evenly distributed along the radar line according to their respective distances (inverse distance weighting). However, the S1, S2 and S2B cores are of limited use because of a maximum depth of 20 m. Moreover, including the effects of density at S0 (potentially different from that at DC) according to an inverse distance weighting would not make much sense when considering the position of the S0 profile with respect to the entire radar profile (the DC–S0 gradient being perpendicular to the major axis of the radar line; see Fig. 1). As a consequence, our strategy thus consists in relying on a single but reliable density profile (DC5) rather than trying to integrate the non-significant geographical effects from limited or not-so-relevant extra data sets.

5.6 Overall error budget for surface accumulation rates

From the expression of snow accumulation rate $\dot{A}=\frac{\stackrel{\mathrm{‾}}{\mathit{\rho }}H}{T}$, the three error sources contributing to the overall uncertainty originate from (i) ΔH, the error in the depth of the IRH; (ii) ΔT, the error in its age; and finally (iii) $\mathrm{\Delta }\stackrel{\mathrm{‾}}{\mathit{\rho }}$, the density error resulting from the use of the selected density profile. All the resulting components are summarized in Table 3. The first error term ($\mathrm{\Delta }{\dot{A}}_{\mathrm{1}}$) due to density implicitly includes the uncertainty resulting from the degree of representativeness and the measurements errors of the selected core (DC5). It has already been assessed from the dispersion in accumulation rates of the numerous available density profiles at DC (see Sect. 5.5) and is directly expressed as an accumulation rate uncertainty (0.6, 0.54 and 0.4 mm w.e. yr⁻¹ for respectively IRHs 1, 2 and 3). As for the error due to the neglect of the geographical changes ($\mathrm{\Delta }{\dot{A}}_{\mathrm{2}}$), it can be estimated to be around the average along the radar line of the difference in accumulation rate when using the DC5 and the S0 core (Sect. 5.4), namely 0.98, 0.85 and 0.68 mm w.e. yr⁻¹ for respectively IRHs 1, 2 and 3. The overall error due to density is finally proposed as the rms of these two terms (see Table 3).

For the age, the two independent terms arising from (i) the relevance and associated uncertainty of the age–depth relationship of the used core ΔT_c (8 years) and from (ii) the positioning error due to radar resolution ΔT_p (30 years; see Sect. 4.3) give a quasi-systematic rms remaining around 30 years. This time error thus transforms into accumulation rate uncertainties of respectively 3.01, 2.06 and 1.39 mm w.e. yr⁻¹ for IRHs 1, 2 and 3.

Finally, errors in the depth (H) of the IRHs stem from the 1.5 m radar positioning accuracy as determined in Sect. 4.3 plus a systematic component due to the time–depth conversion of the two-way travel time. Proper derivation of Eq. (1), with an average error of 18 kg m⁻³ in density measurements (deduced from the estimated density measurement errors represented by the error bars of the DC5 density profile in Fig. 4), gives a relative error of less than 1.5 % for the average wave velocity. It leads to the resulting depth error ΔH_v as a function of the IRH depth as given by Table 3. The rms combination of these two terms finally gives an error component for the accumulation rate of respectively 2.54, 1.68 and 1.23 mm w.e. yr⁻¹ for IRHs 1, 2 and 3.

Following Muller et al. (2010), the cumulative uncertainty is finally proposed as the rms of these three main contributors ($\mathrm{\Delta }\stackrel{\mathrm{‾}}{\mathit{\rho }}$, ΔT and ΔH), which gives global maximum uncertainties of 4.12, 3.02 and 2.17 mm w.e. yr⁻¹ for IRHs 1, 2 and 3 respectively, later rounded to 4, 3 and 2 mm w.e. yr⁻¹. The noticeable feature is a strong decrease with IRH depth, in contrast to Muller et al. (2010), where an increasing error with depth seems to result from a rapid degradation of the signal-to-noise ratio due to the weak penetration depth and high sensitivity to water content of their high-frequency (2.3 GHz) system. It should also be noticed that the proposed uncertainties here pertain to the absolute errors in the accumulation rates and implicitly comprise a significant systematic part that does not come into play when interpreting spatial or time-dependent gradients as in the following section.

6.1 Spatial distribution

A geographical representation of averaged accumulation rates (and differences between the fields) over the three periods characteristic of the three selected IRHs is represented in Fig. 6, where significant trends can be observed. The relatively good phasing between the three IRHs already noticed on the radargram is observable from the difference in the represented fields (blue scale) not exceeding 1.2 mm w.e. yr⁻¹. It is also consistent with the computed SMB of Fig. 5 and confirms a quasi-stationary accumulation pattern over the past 600 years. In particular, an overall decrease in accumulation of some 20 % from DC to the south-west (towards S2) similarly appears for the three IRHs. More specifically, traversing towards the south-east leads to the minimum value between 19 and 20 mm w.e. yr⁻¹ (depending on the averaging period) at a place close to the S2B core site in the direction of a nearby megadune field (see Fig. 1).

Figure 6Geographical (polar stereographic, 71^∘ S) distribution of net accumulation rates according to the ages of the three IRHs. Panel (a) is the average over the last 234 years (1778–2012 CE), whereas panel (b) represents the latter field minus the 1621–1778 CE average. Last, panel (c) stands for the difference between the 1621–1778 and the 1420–1621 CE periods. Globally positive values in these differences confirm a steady overall increase through time as revealed by Fig. 9. Background is a RADARSAT image of the prospected area revealing some topographic features emphasized by the contour line of the surface digital elevation model (Bamber et al., 2009).

Download

Figure 7Snow surface level in millimetres above the reference level of 1 January 2012 obtained from acoustic distance measurements. Presented daily mean values were calculated from measurements recorded every 30 min and corrected for air temperature.

Download

It should be noted that the 2011/12 TASTE-IDEA traverse was the opportunity for the deployment of an ultrasonic range sensor (SR50, Campbell Scientific) allowing for a continuous monitoring of the surface height for more than 5 years at the S2 point. Corresponding surface heights since January 2012 are depicted in Fig. 7, where short-term variations revealing both precipitation and/or strong snow redistribution events are clearly visible. Despite these rapid changes, which eventually contribute to the metre-scale variability in surface accumulation (Libois et al., 2014), a significant trend emerges over these 5 years of measurements. When considering surface densities between 300 and 330 kg m⁻³ (Fig. 4), this trend maps into accumulation rates between 19.5 and 21.5 mm w.e. yr⁻¹, which appear to be fully consistent with our computed accumulation rate for the S2 point over the last 234 years as represented by the black box in Fig. 5.

Previous measurements of GPR-derived accumulation rates in the vicinity of Concordia Station were carried out in 2005 (Urbini et al., 2008). A radar profile through the station revealed an almost north-to-south gradient of $-\mathrm{0.02}\pm \mathrm{0.01}$ mm w.e. yr⁻¹ km⁻¹ in its southern part (see Fig. 6 in Urbini et al., 2008, and Fig. 8, where corresponding data have been summarized). Despite limited areas over which they overlap, and a slightly different averaging period (the last 266 years against our average of the last 234 years), these data nevertheless remain relatively consistent with ours as can be seen from the figure. Moreover, from our data, starting from around DC and considering the upper left blue points of the figure, a difference of about 1 mm w.e. yr⁻¹ can be derived over a distance of some 30 km, leading to a gradient of −0.033 mm w.e. yr⁻¹ km⁻¹: close to that of Urbini et al. (2008) along a similar direction. When assessing an uncertainty in our proposed gradient, only the non-systematic terms have to be considered, such as the effects of a laterally varying density, which in the present case should be insignificant because of the distance of only 30 km. Errors in positioning or age along IRH1 should also remain small. Indeed, as can be seen from the limited dispersion in both ages and depths at crossing points S1 and S0 (see Table 2), the spatial variations in positioning or in the age along a given IRH are minor (maximum differences of 18 cm and 3 years for IRH1, leading to a rms uncertainty of 0.43 mm w.e. yr⁻¹ according to the error budget of Table 3). It leads to an uncertainty of 0.014 mm w.e. yr⁻¹ km⁻¹ on the proposed gradient and makes the comparison still relevant.

Figure 8Comparison of our accumulation rates with those in Urbini et al. (2008) in the vicinity of Concordia Station over the area featured by the dashed inset of Fig. 6. Blue dots represent the rates as represented in Fig. 6 for the last 234 years, whereas the green ones are those from the radar transect of Urbini et al. (2008) after being equally sampled (5 km) and averaged over five measuring points. Also shown by red dots are accumulation rates inferred from an analysis of firn cores based on tritium ∕ β markers for the 1965–2000 period (Urbini et al., 2008) and for which an uncertainty of 10 % has been derived. Contours are altitude above sea level as given by Bamber et al. (2009). The Concordia Station is featured with a yellow star.

Download

It should be noted that this gradient has also been observed by Genthon et al. (2015) from snow accumulation at two stake networks 50 km apart along a north–south direction centred on the Concordia Station. Despite a high uncertainty for the uppermost snow density which limits the derivation of absolute values of accumulation rate, the use of a uniform snow density along the 50 km long section revealed a significant gradient similar to that of Urbini et al. (2008) and of the present study. Their study also shows that, despite an overall underestimation of snow accumulation, meteorological analyses from the European Centre for Medium-Range Weather Forecasts confirm this tendency (see the modelled precipitation field in Fig. 1). Moreover, from a careful observation of precipitation events and analysis of air mass back-trajectory pathways, Genthon et al. (2015) show that the core of snowfall results from relatively warm (and hence humid) air masses intrusions from the north. When passing over the Concordia dome on their way to the south, they undergo a significant temperature-driven depletion in moisture, which is fully consistent with the observed gradients.

Figure 9Time-dependent evolution of the accumulation rates spatially averaged over five sectors along the radar profile. The first sector corresponds to the DC vicinity (green squares on the map) within the firn core network as in Urbini et al. (2008) and featured by the red spots in Fig. 8. The second sector more or less follows the ridge from the end of the DC sector down to a topographic saddle (orange circles). The third (blue diamonds) continues along the ridge and climbs up to the S2 point. The fourth (purple stars) follows the perpendicular S2–S2B transect, whereas the fifth (brown triangles) comprises the vicinity of the S0 point. The time frames for the time averaging (a) correspond to the age difference from each of the three IRHs to the next above or to the surface, namely 1420–1621, 1621–1778 and 1778–2012 CE. Also presented with grey symbols are the accumulation rates as reported by Frezzotti et al. (2005) and Urbini et al. (2008) over the DC area (see text) and for which error bars of ±10 % occur (not shown for clarity). Surface topography with 2 m contours is that of Bamber et al. (2009).

Download

6.2 Time dependency of accumulation

Because of the limited number of IRHs being analysed in this study, only a coarse characterization of the time-dependent accumulation can be proposed as is the case in Urbini et al. (2008) for instance. This limitation comes from the stringent requirements for a reflector to be considered as properly isochronous (showing continuity over large distances) and assigned a relevant age (by passing through a coring site where an depth–age relationship can be independently derived). As for recent accumulation rates, it should be noticed that the multiple direct reflections in the air between the emitter and receiver parts of the antenna lead to a saturated signal screening potential sub-surface reflectors. The last decades of surface mass balance from GPR are therefore often missing, which prevents proper comparisons with “modern” accumulation values as derived from stake farms for example (only operational since the late 1990s for the DC area).

We therefore propose three periods over which mean accumulation rates can be derived (1420–1621, 1621–1778 and 1778–2012 CE) as can be seen from Fig. 9. Data were spatially averaged over five main sectors whose respective extensions are reported in the figure with appropriate symbols. As for the associated uncertainties (depicted in the form of the proposed error bars in the figure), amongst the various terms of Table 3, only the non-systematic components have been kept. For instance error due to density is now restricted to the $\mathrm{\Delta }{\dot{A}}_{\mathrm{1}}$ term since $\mathrm{\Delta }{\dot{A}}_{\mathrm{2}}$ results from potentially unaccounted-for geographical changes in density. It leads to respective accumulation rate uncertainties of 0.6, 0.54 and 0.4 mm w.e. yr⁻¹ for IRHs 1, 2 and 3 respectively. For the age uncertainty, we considered ΔT_c (error of 8 years resulting from the age–depth relationship) and the effects of a reduced positioning uncertainty ΔT_p of 0.5 m (since most of it is time-independent) leading to an age uncertainty of 10 years. Considering the average slope in the age–depth relationship, the global rms age error (13 years) maps into respective uncertainties of 1.3, 0.9 and 0.6 mm w.e. yr⁻¹. Finally, for the error solely due to positioning and its effect on the cumulative mass, we restricted to the above-mentioned 0.5 m and neglected errors stemming from the uncertainty in the velocity profile. It thus leads to accumulation uncertainties of respectively 0.85, 0.64 and 0.47 mm w.e. yr⁻¹ for IRHs 1, 2 and 3. The corresponding total errors finally result from the rms of these three terms yielding the ±1.66, ±1.23 and ±0.86 mm w.e. yr⁻¹ error bars reported in the figure.

The figure shows an overall accumulation increase which appears consistent from one sector to the other. Such a uniform increase is in line with the stationary character in the spatial distribution as already observed in Figs. 5 and 6 for example. Our results compare well with those mentioned in Urbini et al. (2008) and references therein over the DC area (depicted in grey in the figure). More specifically, from nss${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ volcanic spikes along the EPICA EDC96 core, an average accumulation rate of slightly less than 25 mm w.e. yr⁻¹ for the 1460 (Kuwae)–1816 (Tambora) period was proposed (Castellano et al., 2004). Then follows the 1816–1998 AD period with a 25.3 mm w.e. yr⁻¹ accumulation rate overlapped by a marked increase up to 28.3 mm w.e. yr⁻¹ between 1965 and 2000 deduced from nuclear test horizons along the firn cores represented by the red circles in Fig. 8 (Frezzotti et al., 2005). Unfortunately, these last results also suffer uncertainties (estimated to be 10 % of their value), which casts some doubt on the reliability of possible deduced time trends.

Accumulation rates from stake farms are also proposed by Frezzotti et al. (2005) with 32 mm w.e. yr⁻¹ for the 2004–2007 period and up to 39 mm w.e. yr⁻¹ between 1996 and 1999, but their relevance is also questionable because (i) the measuring period is extremely short with respect to a very high interannual variability; (ii) the derived accumulation rate requires proper knowledge of the sub-surface density, which is highly variable and difficult to accurately measure; and (iii) the intrinsically strong spatial variability in snow thickness at the scale of a stake farm also requires long integration periods to gain in significance. This is confirmed by the associated standard deviation of 14 mm w.e. yr⁻¹ for the 1996–1999 stake results. In any case, no comparison with radar data is possible because of corresponding reflectors too close to the upper surface to be detected. More generally, other studies along the East Antarctica Ice Divide (Dome C, Vostok and Dome Fuji) suggest similar increases of surface mass balance during the last century (Frezzotti et al., 2013; Fujita et al., 2011; Osipov et al., 2014; Urbini et al., 2008). If the results proposed here appear consistent with these studies, their associated uncertainties unfortunately remain too large (even after being diminished by the time-independent biases) to serve as incontrovertible evidence of a significant time increase in the accumulation rates.