Brief communication: New radar constraints support presence of ice older than 1.5 Myr at Little Dome C

The area near Dome C, East Antarctica, is thought to be one of the most promising targets for recovering a continuous ice-core record spanning more than a million years. The European Beyond EPICA consortium has selected Little Dome C (LDC), an area ∼ 35 km southeast of Concordia Station, to attempt to recover such a record. Here, we present the results of the final ice-penetrating radar survey used to refine the exact drill site. These data were acquired during the 2019–2020 austral summer using a new, multi-channel high-resolution very high frequency (VHF) radar operating in the frequency range of 170–230 MHz. This new instrument is able to detect reflectors in the near-basal region, where previous surveys were largely unable to detect horizons. The radar stratigraphy is used to transfer the timescale of the EPICA Dome C ice core (EDC) to the area of Little Dome C, using radar isochrones dating back past 600 ka. We use these data to derive the expected depth–age relationship through the ice column at the now-chosen drill site, termed BELDC (Beyond EPICA LDC). These new data indicate that the ice at BELDC is considerably older than that at EDC at the same depth and that there is about 375 m of ice older than 600 kyr at BELDC. Stratigraphy is well preserved to 2565 m, ∼ 93 % of the ice thickness, below which there is a basal unit with unknown properties. An ice-flow model tuned to the isochrones suggests ages likely reach 1.5 Myr near 2500 m, ∼ 65 m above the basal unit and ∼ 265 m above the bed, with sufficient resolution (19 ± 2 kyr m⁻¹) to resolve 41 kyr glacial cycles.

1 Introduction

Ice cores provide one of the best records of paleoclimate on 100 kyr timescales, but to date no continuous ice core has been recovered that spans more than 800 kyr in stratigraphic order. There is great interest in extending ice-core records beyond the mid-Pleistocene transition (1.25 to 0.7 Ma), since this may provide unique insight in the mechanism which caused the switch between 41 and 100 kyr ice-age cycles. An ice core spanning the last ∼ 1.5 Myr would extend into the period characterized by regular 41 kyr cycles and would provide a more precise record of greenhouse gases through this transition than is currently available (Fischer et al., 2013). Several nations or consortia of nations are endeavoring to recover such cores in East Antarctica as part of the International Partnerships in Ice Core Sciences (IPICS; Beyond EPICA near Dome C, Australia near Dome C, China near Dome A, Japan near Dome F, Russia near Ridge B, and the USA in the Allan Hills and exploring other potential sites; see other articles in this special issue for details on these efforts).

The EPICA Dome C ice core (EDC; EPICA Community Members, 2004), drilled at the location now occupied by Concordia Station in East Antarctica, provides the oldest stratigraphic ice-core climate record recovered to date. The site's cold conditions, low accumulation, and thickness are conducive to preserving old ice. However, slight melting at the bed suggests that a nearby site with slightly thinner ice, and thus no basal melt, could preserve a longer record. Ideally, that site would have relatively smooth bed topography to prevent flow-induced disturbances. Modeling identified two candidate targets in the area (Parrenin et al., 2017), and subsequent work (Passalacqua et al., 2018; Young et al., 2017) narrowed the search to an area ∼ 35 km southwest termed Little Dome C (LDC, Fig. 1). To obtain the oldest ice at maximum resolution, the core would ideally be at a location where the ice is as thick as possible without allowing basal melting. While the minimum ice thickness to allow melting varies spatially with accumulation, ice flow, and geothermal heat flux, several constraints are available for the region. Analysis of airborne radar data identified a number of subglacial lakes, all of which lie beneath at least 2875 m of ice (Young et al., 2017). Though not framed specifically in terms of minimum ice thickness to cause melting, several thermal modeling studies suggest that parts of LDC with ice thickness around 2700 m are likely free of basal melt (Passalacqua et al., 2018; Parrenin et al., 2017; Van Liefferinge et al., 2018).

Extensive radar work has been conducted in the area of LDC in the framework of Beyond EPICA, which greatly narrowed the area this present work examined. The initial aerial survey (Young et al., 2017) mapped the bedrock extensively, greatly improving the knowledge of the bed compared to the single Operation IceBridge flight line in the area. These results also allowed further inference of basal conditions (Passalacqua et al., 2017; Van Liefferinge et al., 2018) as well as the accumulation rate in the area over the last 73 kyr (Cavitte et al., 2018). Those results led to a targeted ground-based survey using an impulse radar operating in the 1–5 MHz range (Cavitte et al., 2020), which narrowed the core location to an ∼ 8 km² area. While modeling shows that LDC is likely to have old ice (e.g., Van Liefferinge and Pattyn, 2013), those previous radar surveys were not able to connect any isochrones older than ∼ 400 kyr throughout the area to the EDC ice core, leaving large uncertainties in the ages nearer the bed. This inability to trace older horizons stems from a near-basal region that is common in Antarctica, previously described as the “echo-free zone”. The cause of the echo-free zone is unclear and has been variously attributed to a sharp thermal transition, folding or buckling, and re-circulation and re-crystallization (see Drews et al., 2009, for a detailed discussion). Indeed, the existence of an echo-free zone is disputed, as it may simply be an artifact of radar-system detection limits. Regardless of the cause, the lack of reflections near the bed in prior surveys of LDC limited constraints on very old ice in the region. In addition to the echo-free zone, some radargrams showed a diffuse horizon near the bed, hinting that basal ice in the area may have different physical properties than the overlying ice (Cavitte, 2017), though the implications for the depth relationship were unclear.

Here, we present the results of an additional radar survey designed to connect the stratigraphy of that site with the EDC core and to identify the area at LDC with the highest potential for old ice. This survey utilized a new, highly sensitive radar, which allowed for detection of older horizons nearer the bed. The exact location for the ice core, Beyond EPICA LDC (BELDC; 75.29917^∘ S, 122.4451^∘ E; 3230 m above the WGS84 ellipsoid as of 2020), has deep, flat, visible stratigraphy and lies within the region of LDC identified by previous studies to be free of melt and likely to contain old ice. Here, we present the age constraints provided by this new radar survey for that chosen site by comparing the stratigraphy to EDC.

Figure 1 Map of study area, with inset showing location in Antarctica. Black lines show radar profiles acquired in 2019, with bold line showing the profile evaluated here. EDC core coincides with Concordia. Contours show surface elevation from Helm et al. (2014). Background colors are bed elevation relative to sea level from BedMachine Antarctica v2 (Morlighem et al., 2020).

2 Methods

2.1 Data collection and processing

Data were collected using a new very high frequency (VHF) radar, built by the Remote Sensing Center at the University of Alabama (Yan et al., 2020). The system was configured to transmit 8 µs chirps, with 200 MHz center frequency and 60 MHz bandwidth. Peak transmit power was varied from 125–250 W through the campaign to maximize the signal-to-noise ratio while limiting problems with radio-frequency interference. The system has eight transmit and receive channels, paired with eight monostatic antennas. Due to the logistical challenges of the operating environment, the number of channels in use varied from five to eight. The system was pulled behind a tracked vehicle, with controlling electronics in the rear passenger compartment and antennas approximately 12 m behind. The antennas were set up such that the electric field polarization was oriented across track, above a single sheet of plywood for stiffness and thin PVC mat for slipperiness. Data were collected at travel speeds of 2–3.5 m s⁻¹ over the course of a week in November and December 2019.

Data processing consisted of coherent integration (i.e., unfocused SAR), pulse compression, motion compensation (by tracking internal horizons), coherent channel combination, and de-speckling using a median filter. Two-way travel time was converted to depth assuming a correction of 10 m of firn air and a constant radar wave speed in ice of 168.5 m µs⁻¹ (e.g., Winter et al., 2017). After other processing was complete, different radargrams were spliced together to create a continuous profile extending from EDC to BELDC, and then the data were interpolated to have constant 10 m horizontal spacing. The re-interpolated data were used for horizon tracing, which was done semi-automatically to follow amplitude peaks between user-defined clicks. For the bed reflection, there were often weak, diffuse events shallower than a clear return. We always picked the first notable return in the region of the bed, so ice-thickness estimates are likely biased shallow; the number of hyperbolic and diffuse events, likely originating from roughness of the ice–bed interface or physical changes in the lower parts of the ice sheet, would cause a high risk of misinterpretation with other approaches.

2.2 Horizon dating and depth–age reconstruction

Radar reflectors were dated by interpolating from the AICC2012 timescale (Bazin et al., 2013) at the point of closest approach to the EDC drill site. The radar line ended approximately 100 m horizontally from the EDC borehole. The depth of the bed reflection there is 3238 m (38.32 µs), within the depths at closest approach to EDC found by other radar systems (approximately 3220–3286 m in Winter et al., 2017) but shallower than the minimum actual thickness implied by the 3260 m EDC borehole (Parrenin et al., 2007). This offset is likely due to some combination of off-nadir reflection, debris in the ice, small differences in topography over the 100 m offset, and uncertainty in firn-air content and wave speed. Regardless of the cause, we must either re-scale the thickness to match EDC or leave it as measured. We choose the latter since any re-scaling would be highly uncertain.

The dating uncertainty has two primary components: uncertainty in the ice-core timescale and uncertainty in the radar-horizon depth. The horizon-depth uncertainty can be further subdivided into the component caused by the radargram not extending exactly to the EDC core site and the component caused by the firn correction and dielectric constant (see Winter et al., 2017, and references therein for a detailed discussion of the components of the error). For the ice-core uncertainty, we use the previously published estimates from the chronology (Bazin et al., 2013). We estimated slope-induced uncertainty from the ∼ 100 m offset of the radargram from the core using each horizon's average slope; slopes ranged from 10 to 60 m km⁻¹, resulting in depth uncertainty of 1 to 6 m, increasing with depth. The depth uncertainty introduced by anisotropy and temperature affecting the dielectric constant is taken to be 1 %, and we assume an additional 3 m uncertainty in the firn-air correction. The formal quarter-wavelength uncertainty of the horizon position is small (0.2 m) compared to other terms. Thus, total depth uncertainties range from 11 m for the upper horizons to 31 m for the lower horizons, introducing age uncertainties of 1 to 33 kyr (found using the depth gradient of age following Winter et al., 2017). Combining with uncertainties in the timescale itself, total age uncertainties increase from 2 kyr for shallow horizons to 34 kyr for the deepest ones. The uncertainties of the horizons' ages are correlated with each other, since an incorrect firn-air correction or dielectric constant, or an incorrect age scale at EDC, affects the inferred age of all these horizons, though we are unable to quantify the extent of this correlation.

2.3 Modeling the depth–age scale

While previous work has used sophisticated models to make estimates of the depth–age scale at LDC (Parrenin et al., 2017), here we seek a more simple constraint relying on the tighter age bounds of the new radar data. We fit a modified Lliboutry model (Lliboutry, 1979) to the horizons at the chosen core site. The model provides an analytic solution for the vertical thinning function, λ(d), where d is depth, the inverse of which can be integrated to find the steady-state age, $\overline{t}$, at a given depth. That is,

$$\begin{array}{}\text{(1)}& \overline{t}\left(d\right)=\underset{\mathrm{0}}{\overset{d}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }\left(d^{\prime }\right)}}\mathrm{d}{d}^{\prime }.\end{array}$$

Assuming no basal melt or accretion, and ignoring firn, the Lliboutry model approximates layer thickness as

$$\begin{array}{}\text{(2)}& \mathit{\lambda }\left(d\right)=a\left(\mathrm{1}\mathrm{-}{\displaystyle \frac{p+\mathrm{2}}{p+\mathrm{1}}}\left({\displaystyle \frac{d}{H}}\right)+{\displaystyle \frac{\mathrm{1}}{p+\mathrm{1}}}{\left({\displaystyle \frac{d}{H}}\right)}^{p+\mathrm{2}}\right),\end{array}$$

where a is ice-equivalent accumulation, H ice thickness, and p a shape factor controlling the vertical profile of deformation.

As in Parrenin et al. (2017), we used a temporally variable accumulation rate and solved for the depth–age relationship using a pseudo-steady method which permits analytical solutions for λ even with temporally variable accumulation (Parrenin et al., 2006). This involves a simple change of variable between time, t, and steady time, $\overline{t}$, of the form

$$\begin{array}{}\text{(3)}& \overline{t}=\underset{\mathrm{0}}{\overset{t}{\int }}R\left(t^{\prime }\right)\mathrm{d}{t}^{\prime },\end{array}$$

where $R\left(t\right)=a_{\mathrm{E}}\left(t\right)/{\overline{a}}_{\mathrm{E}}$ is the normalized accumulation at a given time of the EDC record (Bazin et al., 2013), and we assume R=1 for ages older than the extent of the EDC record. Equation (3) defines a bijection between t and $\overline{t}$, so we can first find the steady-state age profile using Eq. (1) and then convert to the equivalent profile incorporating the EDC accumulation variations using Eq. (3). In this formulation, the temporally variable accumulation enters only as the non-dimensional scaling, R(t), while a in Eqs. (1)–(2) is treated as a constant.

We used a Markov chain Monte Carlo method, implemented with PyMC3 (Salvatier et al., 2016), to find the probability distribution of the resulting depth–age scale by varying a, H, and p. We allowed H to vary to account for the possibility of stagnant ice which does not affect the deformation of the overlying ice column. Aside from the minor difference of using a newer Monte Carlo sampler and fitting for the effective thickness, H, the essential difference between this model and that used by Parrenin et al. (2017) is that we have excluded the possibility of basal melt and thus all thermal modeling; this is justified by the previous thermal modeling in the area (Passalacqua et al., 2018; Parrenin et al., 2017; Van Liefferinge et al., 2018) as well as radar evidence from other surveys (Young et al., 2017). This model strikes a balance between numerical requirements and a realistic, singularity-free profile of strain at different depths (Parrenin et al., 2006).

3 Results

The processed radargram shows a clear bed reflection and a number of horizons that can be continuously traced from EDC to LDC (Fig. 2). The data indicate that the thickness at the chosen site at LDC is 2764 ± 20 m, with the top 2565 ± 20 m showing continuous stratigraphy. Around 2565 m there is a change in the amplitude of the radar returns that suggests ice below this depth has different properties; this feature is discussed further in Sect. 4.1. The ice thickness is >100 m less than the minimum thickness over any subglacial lake observed in the area (Young et al., 2017), and there is no non-bed-parallel down-warping of englacial horizons, both indicating that the site is free of basal melt. To the east of LDC (near kilometer 14 and 27 in Fig. 2) there is some down-warping of englacial horizons, but the ice in that area is thicker than at BELDC as the bed deepens in subglacial valleys.

We traced a subset of the visible horizons, selected to span all depths with a concentration in the deepest areas. In addition to the continuous horizons, there were three deep horizons that could be identified near both ends of the radargram, where horizon slopes are relatively flat, but not in the middle of the radargram, where horizon slopes were steeper. We also traced these three partial horizons where possible, to attempt a more complete connection of the core's timescale to the BELDC site. Despite the gaps in the middle of the profile, the relatively distinctive pattern of these horizons gives us confidence that they are the same isochrone at each end of the profile (see Fig. 2b–e for zoom-ins on these horizons). However, since it is possible that we have misidentified these horizons, the depth–age analysis in Sect. 4.2 is repeated with the three deepest horizons excluded.

At the intersection with EDC, horizons were dated from 71 to 565 ka (Fig. 2). Of these, four isochrones were older than 399 kyr, the age of the oldest previously isochrone previously dated and traced from EDC (Parrenin et al., 2017). The oldest isochrone that could be continuously traced to BELDC was dated to 465 ka. This increase in the age of dated isochrones introduces a significantly tighter constraint upon the age and age resolution at the LDC site than was previously available. These data indicate that ice reaches 600 kyr at a relatively shallow 2373 ± 20 m at the chosen BELDC site. The older radar ages lead us to infer ice older than 1.5 Myr, with ∼ 19 kyr m⁻¹ resolution as detailed in Sect. 4.2.

Figure 2 (a) Radargram extending from EDC to the chosen ice-core site, BELDC. Blue and orange lines show horizons traced by Winter et al. (2017) and newly traced here, respectively, with ages on the right. Marks at the top show turns in the profile, as noted on Fig. 1. (b–e) Zoom-ins of discontinuous horizons at locations in panel (a).

4 Discussion

We first discuss a thick basal unit at LDC that may affect the depth–age relationship, before focusing on the depth–age distribution at the BELDC site.

4.1 Basal unit

At LDC, there is a notable event (i.e., change in the return power) in the radargram around 2565 m depth (pink line in Fig. 2). While there are some distinct peaks in return power that arrive later than this event, and thus most likely originate deeper, there are no continuous or coherent reflecting horizons below this depth. However, this event is at a depth where the radar is capable of imaging continuous reflections throughout most of the radar profiles, suggesting a change in ice properties rather than a system detection limit. Other radars noted a diffuse event at this depth (Cavitte, 2017), but the additional sensitivity of this new survey helps establish that the origin of the event is in fact a physical property of the ice and not an artifact. It may be caused by a sharp transition in crystal fabric, heterogeneous small-scale roughness, stagnant ice, tightly spaced or disrupted isochrones, or perhaps some other relic feature; regardless of its origin, the lack of stratigraphy suggests that recovery of a climate record below this depth may be difficult, potentially impossible. However, despite the lack of continuous radar isochrones, studies at other ice-core sites (e.g., EDC and EDML) showed that climate information might still be retrievable at least in the top part of basal units with similar characteristics (e.g., Tison et al., 2015).

4.2 Depth–age scale at the BELDC site

The depth–age scale from a random sample of parameters drawn from the posterior distributions of a, H, and p is shown by the gray bands in Fig. 3b. The green bars show the results with only continuous horizons used to constrain the model. The 95 % confidence interval with the discontinuous horizons included is a subset of the 95 % confidence interval using continuous horizons only; including discontinuous horizons narrows the distribution to be slightly younger and to have slightly higher depth–age resolution. For the rest of this section, stated values include the discontinuous horizons as a constraint, while parentheticals use continuous horizons only. Uncertainties in this section are the standard deviation of the distribution of modeled depth–age profiles.

The best-fit ice thickness is 2579 ± 22 m (2549 ± 32 m), within uncertainty of the top of the basal unit. The mean accumulation is 14.5 ± 0.2 kg m⁻² a⁻¹ (14.4 ± 0.2 kg m⁻² a⁻¹), close to the 13.9 kg m⁻² a⁻¹ average found in the EDC core (Bazin et al., 2013). The shape factor is 5.5 ± 1 (6.5 ± 1.5). The relatively small effective ice thickness suggests that the basal unit is partially stagnant or flows much more slowly than the overlying ice, such that the deformation of the overlying ice column is unaffected by this deeper ice and ages reach an asymptote near the top the basal unit. The corresponding depth–age scales reach 1.5 Myr at 2498 ± 14 m (2476 ± 25 m), with average age resolution of 19 ± 2 kyr m⁻¹ (21 ± 3 kyr m⁻¹). Including the discontinuous horizons in the age modeling thus suggests that resolution is sufficient for measuring 41 kyr glacial cycles at this age, while using continuous horizons only suggests that the resolution may be marginal at 1.5 Myr. The original Beyond EPICA resolution target, 20 kyr m⁻¹ (Fischer et al., 2013), is passed at age 1552 ± 83 kyr (1468 ± 101 kyr). Regardless of whether discontinuous horizons are used, our modeling suggests that conditions are close to the desired depth–age properties, though the deepest isochrones add confidence that the site fully meets the age and resolution targets. Further refinement of the depth–age scale with more sophisticated models and elaborate assumptions could refine the estimated age and age resolution but will likely result in broadly similar values.

Figure 3 (a) Radargram near BELDC (as in Fig. 2). (b) Depth–age scale at BELDC and EDC. Dots with error bars show traced horizons with age uncertainties; blue line shows AICC2012 chronology (Bazin et al., 2013). The gray shaded regions show model results utilizing constraints from all horizons at 95 % confidence (light) and 1σ (dark). Green shaded region shows the model results with the discontinuous horizons excluded from the analysis.

5 Conclusions

Newly collected radar data provide a tighter constraint on the depth–age scale at LDC. These data reveal traceable stratigraphy in ice >500 kyr old in the region, with continuous horizons to 465 kyr. Near LDC, they also indicate a unit of basal ice in which few events are visible; the origin of this basal unit requires further investigations as its flow properties and composition are unknown. The stratigraphy indicates that old ice lies much shallower at BELDC than at EDC. A Lliboutry model of the depth–age scale, fitted to the isochrone data, indicates that 1.5 Myr old ice lies at ∼ 2500 m depth, where stratigraphy is still intact and preserved with ∼ 19 kyr m⁻¹ resolution. Very old (>1.5 Myr) ice could exist atop the basal unit, which appears partially stagnant and presumably also contains >1.5 Myr old ice, though the lack of stratigraphy does not allow firm conclusions regarding the extent to which useful climatic information may or may not be preserved below 2565 m depth.