Is there 1.5-million-year-old ice near Dome C, Antarctica?

. Ice sheets provide exceptional archives of past changes in polar climate, regional environment and global atmospheric composition. The oldest dated deep ice core drilled in Antarctica has been retrieved at EPICA Dome C (EDC), reaching ∼ 800 000 years. Obtaining an older paleoclimatic record from Antarctica is one of the greatest challenges of the ice core community. Here, we use internal isochrones, identiﬁed from airborne radar coupled to ice-ﬂow modelling to estimate the age of basal ice along transects in the Dome C area. Three glaciological properties are inferred from isochrones: surface accumulation rate, geothermal ﬂux and the exponent of the Lliboutry velocity proﬁle. We ﬁnd that old ice ( > 1.5 Myr, 1.5 million years) likely exists in two regions: one ∼ 40 km south-west of Dome C along the ice divide to Vostok, close to a secondary dome that we name “Little Dome C” (LDC), and a second region named “North Patch” (NP) located 10–30 km north-east of Dome C, in a region where the geothermal ﬂux is apparently relatively low. Our work demonstrates the value of combining radar observations with ice ﬂow modelling to accurately represent the true nature of ice ﬂow, and understand the formation of ice-sheet architecture, in the centre of large ice sheets.


Introduction
Since around 800 000 years ago, glacial periods have been dominated by a ∼ 100 000-year cyclicity, as documented in multiple proxies from marine, terrestrial and ice core records (Elderfield et al., 2012;Jouzel et al., 2007;Lisiecki and Raymo, 2005;Loulergue et al., 2008;Lüthi et al., 2008;Wang et al., 2008;Wolff et al., 2006).These data have provided evidence of consistent changes in polar and tropical temperatures, continental aridity, aerosol deposition, atmospheric greenhouse gas concentrations and global mean sea level over numerous glacial cycles.Conceptual models (Imbrie et al., 2011) have been proposed to explain these asymmetric 100 000-year cycles in response to changes in the configuration of the Earth's orbit and obliquity (Laskar et al., 2004), and involve threshold behaviour between different climate states within the Earth system (Parrenin and Paillard, 2012).The asymmetry between glacial inceptions and terminations may, for example, be due to the slow build-up of ice sheets and their rapid collapse once fully developed due to glacial isostasy (Abe-Ouchi et al., 2013).Observed sequences of events and Earth system modelling studies (Fischer et al., 2010;Lüthi et al., 2008;Parrenin et al., 2013;Shakun et al., 2012) have shown that climate-carbon feed-F.Parrenin et al.: Is there 1.5-million-year-old ice near Dome C, Antarctica? backs also play a major role in the magnitude of glacialinterglacial transitions.
Critical to our understanding of these 100 000-year glacial cycles is the study of their onset, during the Mid-Pleistocene Transition (MPT; Jouzel and Masson-Delmotte, 2010), which occurred between 1250 and 700 kyr BP thousands of years before 1950; Clark et al., 2006), and most likely during Marine Isotope Stages (MIS) 22-24, around 900 kyr BP (Elderfield et al., 2012).Prior to the MPT, marine sediments (Lisiecki and Raymo, 2005) show glacial-interglacial cycles occurring at obliquity periodicities (40 kyr) and with a smaller amplitude.The exact cause for this MPT remains controversial and several mechanisms have been proposed, including the transition of the Antarctic ice sheet from a wholly terrestrial to a part-marine configuration (Raymo et al., 2006), a hypothesis which is, however, unsupported by long-term simulations (Pollard and DeConto, 2009); a nonlinear response to weak eccentricity changes (Imbrie et al., 2011); merging of North American ice sheets (Bintanja and Van de Wal, 2008); changes in sea ice extent (Tziperman and Gildor, 2003); a time varying insolation energy threshold (Tzedakis et al., 2017); a threshold effect related to the atmospheric dust load over the Southern Ocean (Martínez-Garcia et al., 2011); and a long-term decrease in atmospheric CO 2 concentrations (Berger et al., 1999), the latter hypothesis being challenged by indirect estimates of atmospheric CO 2 from marine sediments (Hönisch et al., 2009).
A continuous Antarctic ice core record extending back at least to 1.5 Myr BP would shed new light on the MPT reorganisation (Jouzel and Masson-Delmotte, 2010), by providing records of Antarctic temperature, atmospheric greenhouse gas concentrations and aerosol fluxes prior to and after the MPT.The opportunity to measure cosmogenic isotopes ( 10 Be) would also provide information on changes in the intensity of the Earth's magnetic field, especially during the Jaramillo transition (Singer and Brown, 2002).Retrieving Antarctica's "oldest ice" is therefore a major challenge of the ice core science community (Brook et al., 2006).A necessary first step towards this goal is to identify potential drilling sites based on available information on ice-sheet structure and accompanying age modelling (Fischer et al., 2013;Van Liefferinge and Pattyn, 2013).
The maximum age of a continuous ice core depends on several parameters (Fischer et al., 2013).Mathematically, the age χ of the ice at a level z above bedrock can be written as follows: where D(z) is the relative density of the material (< 1 for the firn and = 1 for the ice), a(z) is the accumulation rate (initial vertical thickness of a layer, in metres of ice yr −1 ), τ (z) is the vertical thinning function, i.e. the ratio of the vertical thickness of a layer in the ice core to its initial vertical thickness at the surface, and H is the total ice thickness.Increasing the maximum age χ max can be obtained by increasing H or by decreasing a or τ .At first glance, one might select a site where H is maximum and a is minimum, but this neglects the importance of τ , notably through basal melting.In general, τ decreases toward the bed and, in steady state, reaches the value µ = m/a, where m is the basal melting.m is therefore a crucial parameter of the problem, as it destroys the bottom of the ice record.As ice is a good insulator, H either increases the ice temperature towards melting for frozen basal ice conditions, or, when melting is present, m increases with H and with the geothermal flux underneath the ice sheet (Fischer et al., 2013).Consequently, "oldest-ice" sites have a better chance to exist where ice is not overly thick as to lead to basal melting (Seddik et al., 2011), yet thick enough to contain a continuous ancient accumulation.The distance of a site to the ice divide is also an important parameter.This distance influences the profile of τ , which is increasingly non-linear right at a dome.Therefore, χ max can be up to 10 times larger at a dome than a few kilometres downstream (Martín and Gudmundsson, 2012).Moreover, assuming a largely constant ice sheet configuration across glacial cycles, an ice record close to the divide has travelled a shorter horizontal distance and therefore has a better chance of being stratigraphically undisturbed (Fischer et al., 2013).
The depth-age profile in an ice sheet can be obtained using radar observations at VHF ranges to identify englacial reflections (e.g.Fujita et al., 1999) and trace them as isochrones across the ice sheet (Cavitte et al., 2016;Siegert et al., 1998).Until now, such analysis has been restricted to the top ∼ three-fourths of the ice thickness in East Antarctica.However, depth-age information from internal layers can be used in conjunction with ice flow models and age information from ice cores to extrapolate down to the bed.Radar observations allow estimates of poorly known ice-sheet parameters, such as the geothermal flux (Shapiro and Ritzwoller, 2004) and past changes in the position of ice domes and divides.
The Dome C sector is one of the target areas for the "oldest-ice" challenge and has a number of distinct benefits over other regions: it has already been heavily surveyed by geophysical techniques (Cavitte et al., 2016;Siegert et al., 1998;Tabacco et al., 1998), a reference age scale has been developed through the existing ice core work (Bazin et al., 2013;Veres et al., 2013) and it is logistically accessible from nearby Concordia Station.In this study, we concentrate on airborne radar transects (Fig. 1), which are all related to the EDC ice core.These data resolve the bed (Young et al., 2017) and internal isochrones (Cavitte et al., 2017) and are suitable for the oldest-ice search (Winter et al., 2017).The isochrones are dated up to about 366 kyr BP using the most recent AICC2012 chronology established for the EDC ice core (Bazin et al., 2013;Veres et al., 2013).We extrapolate the age of the isochrones toward the bed using an ice flow model in order to identify potential oldest-ice sites along these transects.We also build maps of surface accumulation rate, geothermal flux and of a linearity parameter of the ver-  (Fretwell et al., 2013) while the thin grey transparent lines represent the surface elevation (Fretwell et al., 2013).The red square in the inset show the location of the zoomed map around EDC.The red star is the location of the EDC drilling site.The orange squared areas are oldest-ice candidates from Van Liefferinge and Pattyn ( 2013).The red dotted line is the OIA/JKB2n/X45 radar line displayed in Fig. 3. tical velocity profile.The spatial and temporal variations of surface accumulation rates are discussed in detail in a companion paper (Cavitte et al., 2017).

Method
We use a 1-D pseudo-steady (Parrenin et al., 2006) ice flow model, which assumes that the geometry, the shape of the vertical velocity profile, the ratio µ = m/a and the relative density profile are constant in time.Only a temporal factor R(t) is applied to both the accumulation rate a and basal melting m: where a(x) and m(x) are the temporally averaged accumulation and melting rates at a certain point x.Under the pseudosteady assumption, the vertical thinning function is given by where ω is the horizontal flux shape function (Parrenin et al., 2006).While there is no physical reason to assume co- .R(t) proportionality factor applied to accumulation and melting rates (see Eq. 2).The plot is cut at 1 Myr for better readability.R(t) is based on the accumulation record at EDC for the last 800 kyr (Bazin et al., 2013;Veres et al., 2013).
variance of basal melting and surface accumulation, comparison with a transient dating model (Parrenin et al., 2007) shows errors of only 6 % maximum in the evaluation of the thinning function.Moreover, the fact that there is an analytical expression for the thinning function allows one to drastically reduce the computation time, an important factor since the 1-D model needs to be applied on many locations and with many different sets of parameters.A steady age χ steady is first calculated assuming a steady accumulation a and a steady melting m.Then the real age χ is calculated using the following equation (Parrenin et al., 2006): R(t) (Fig. 2) is directly inferred from the accumulation record of the EDC ice core (Bazin et al., 2013;Veres et al., 2013).Beyond 800 kyr BP, it is assumed to be equal to 1; that is to say that the accumulation before 800 kyr BP is assumed equal to the average accumulation over the last 800 kyr.The horizontal flux shape function is determined using an analytical expression (Lliboutry, 1979;Parrenin et al., 2007): where ζ = z/H is the normalised vertical coordinate (0 at bedrock and 1 at surface) expressed in ice equivalent, and p a parameter modifying the non-linearity of ω (the smaller p, the more non-linear ω).This formulation assumes that there is a negligible basal sliding ratio, as is the case at EDC (Parrenin et al., 2007).This might not be the case elsewhere, but adding a basal sliding term has a similar effect as increasing the p parameter for the top ∼ three-fourths of the ice sheet.The age of the ice at any depth is deduced from Eq. (1) using the relative density profile at EDC (Bazin et al., 2013).
To compute the basal melting, we use a simple steady-state 1-D thermal model.We solve the heat equation (neglecting the heat production by deformation since there is minimal horizontal shear) as follows: where T is the temperature, u z is the vertical velocity, ρ i = 917 kg m −3 is the ice density (Cuffey and Paterson, 2010) and k T (W m −1 K −1 ) is the thermal conductivity (Cuffey and Paterson, 2010), is given by and c (J kg −1 K −1 ), the specific heat capacity (Cuffey and Paterson, 2010) is given by The boundary conditions are where T S = 212.74K is the average temperature at the surface assumed to be equal to the one at Dome C (Parrenin et al., 2013), G 0 is the geothermal flux and T f , the fusion temperature is given by Ritz (1992): where P = 10 6 Pa is the partial pressure of air and P , the pressure, is approximated by the hydrostatic pressure: where g = 9.81 m s −2 is the gravitational acceleration.We used this formula since it gives the best agreement to the measured temperature profile at EDC (Passalacqua et al., 2017).The basal melting is given by where G is the vertical heat flux at the base of the ice sheet and L f = 333.5 kJ kg −1 is the latent heat of fusion (Cuffey and Paterson, 2010).
To prevent p from being < −1 (Eq. 5 has a singularity for p = −1), we write The values of a, G 0 and p are reconstructed using a variational inverse method and using the radar isochrone constraints.The cost function to minimise is formulated using a least-squares expression: where N is the number of isochrones (3 ≤ N ≤ 18, see Table 1 and Fig. 3), d iso i and χ iso i are the depths and ages of the isochrones respectively, σ iso i is the confidence interval on their age and χ mod is the modelled age.p prior = ln(1 + 1.97) is the a priori estimates of p , inferred from the age-scale model of the EDC ice core (Parrenin et al., 2007) and σ p = 2 is its standard deviation, chosen to be sufficiently large to allow for a large range of p values.G 0,prior = 51 mW m −2 is the a priori estimate of the geothermal flux calculated from satellite magnetic data (Fox Maule et al., 2005;Purucker, 2013), and from analysis of the heat required to maintain melting above subglacial lakes (Siegert and Dowdeswell, 1996).σ G 0 = 25 mW m −2 is the uncertainty in the geothermal flux (Fox Maule et al., 2005;Purucker, 2013); the total uncertainty of the age of isochrones σ iso is composed of (1) the uncertainty in the depth of the traced isochrones (Cavitte et al., 2016), transferred in age, and (2) the uncertainty of the AICC2012 age of the isochrone at the EDC site.
To solve the least-squares problem formulated in Eq. ( 16), we used a standard Metropolis -Hastings algorithm (Hastings, 1970;Metropolis et al., 1953) with 1000 iterations.This allows one not only to obtain a most probable modelling scenario, but also to quantify the posterior probability distribu-  tion, in particular the confidence intervals or the modelled quantities.

Results and discussions
In our forward modelling, we used the 1-D pseudo-steady assumption.This assumption is very convenient numerically because in this case, there is an analytical expression for the thinning function (Eq.3).Therefore, there is no need to use a costly Lagrangian scheme, following the trajectories of ice particles.Of course, the reality is more complex than the pseudo-steady assumption because the temporal variations in melting and accumulation rates are not related and are not the same for each point in space.In Parrenin et al. (2007), we used a more complex age model with a ratio µ and with an ice thickness allowed to vary in time.The results are very similar with the pseudo-steady model.This is because melting is small compared to the accumulation, and the variations in ice thickness are small compared to the total ice thickness.Regarding the spatial pattern of accumulation, we assumed that it is stable in time, which is roughly confirmed by the inversion of internal layers (Cavitte et al., 2017).Moreover, the 1-D assumption dominates the uncertainty since we do not take into account horizontal advection and dome movement.Therefore, we suggest the pseudo-steady assumption is good enough for a 1-D model.
An example age profile along the OIA/JKB2n/X45 radar transect (see Fig. 1 for its position) is displayed in Fig. 3. From these profiles, maps of the modelled age at 60 m above the bed, minimum age at 60 m above the bed (at 85 % confidence level), the height above the bed of the 1.5 Ma isochrone and temporal resolution at 1.5 Myr are displayed in Fig. 4. We use 60 m above the bed as this is the height at EDC below which the ice becomes disturbed such that it cannot be interpreted stratigraphically (Tison et al., 2015).The modelled basal melting m and inferred steady accumulation rate a, geothermal flux G 0 and p parameter of the vertical velocity profile are displayed in Fig. 5.
The bottom age inferred at EDC at 3200 m is 785 kyr, which is remarkably close to the age of ∼ 820 kyr inferred from the analysis of the ice core (Bazin et al., 2013;Veres et al., 2013).This 35 kyr difference represent a depth mismatch of 24 m.This is a confirmation of the method used, despite its assumptions (i.e.1-D, pseudo-steady, Lliboutry velocity profile).
There are two main regions where the basal age is modelled to be older than 1.5 Myr.The first one is situated close to Little Dome C (LDC), ∼ 40 km south-west of EDC.In this region that we call LDC Patch (LDCP), the ice thickness is several hundred metres lower than at EDC, thus reducing the likelihood of basal melting.The second region is 10-30 km north-east of EDC in the direction of the coast, at a place where the ice thickness is comparable to the one at EDC but with a lower geothermal flux.We call this re-gion "North Patch" (NP).In those two oldest-ice spots, the height above the bed of the 1.5 Myr isochrone is modelled to be greater than 150 m.The temporal resolution at 1.5 Myr is ∼ 10 kyr m −1 , which is sufficient to resolve the main climatic periods (Fischer et al., 2013).
Our LDCP area is generally consistent with Candidate A of Van Liefferinge and Pattyn (Van Liefferinge and Pattyn, 2013) although our area is smaller and constrained to the subglacial highlands under LDC.Van Liefferinge and Pattyn (2013) did not find a candidate at NP.However, the geothermal heat flux maps they relied on have a lower spatial resolution than the details we examine here.Our model does not predict very old ages for Candidates B-E of Van Liefferinge and Pattyn (2013), although the 1-D assumption is problematic in those areas since ice particles experienced very different ice thickness conditions along their path.
One possible limitation of our simple ice sheet model is that it does not allow for a layer of accreted ice.We argue that there are no discernable accretion features in the UTIG radargrams, although it is possible that the accretion features do not show up in the basal layer which is difficult to interpret.
We now examine the other variables inferred from the inversion.Basal melting is of course negligible at these two oldest-ice spots.Melting is, however, significant around EDC (which is consistent with known basal melting at this place), on the other side of LDC and on the bed ridge adjacent to the Concordia Subglacial Trench (called here the Concordia Ridge), consistent with the observation of subglacial lakes (Wright and Siegert, 2012;Young et al., 2017).While it is surprising that basal melting is so large across the ridge of the bed, where the ice thickness is smaller, the 1-D assumption is probably invalid in this region, since the ice has been significantly advected horizontally over regions with very different basal conditions (i.e. over the wet-based Concordia Subglacial Trench and then over the adjacent Concordia Ridge which likely has a frozen base).The average surface accumulation rate shows a large-scale north-east-south-west gradient probably linked to a precipitation gradient.It also shows small-scale variations linked to surface features and probably due to snow redistribution by wind.The spatial and temporal variations of accumulation are the subject of a companion paper to this study (Cavitte et al., 2017).For the geothermal flux, it should be noted that its reconstruction is only relevant when there is some basal melting (i.e. a temperate base).When the base is cold, its evaluation mainly relies on the prior used for the least-squares cost function.Indeed, below the threshold of zero melting, further decreasing the geothermal flux has no effect on the basal melting, and therefore no effect on the modelled age.In the EDC region, the geothermal flux is estimated around 60 mW m −2 .A high geothermal flux of ∼ 80 mW m −2 is also estimated on the ridge adjacent to the Concordia Subglacial Trench.Again, these results should be taken with caution since they could be an artifact due to the 1-D assumption used.The p value inferred

Conclusions
We developed a simple 1-D thermo-mechanical model constrained by radar observations to infer the age in an ice sheet.We identified two areas where the age of basal ice should exceed 1.5 Myr.They are located only a few tens of kilometres away from the French-Italian Concordia station, which could provide excellent logistical support for deep drilling.
The first area, LDCP, is close to a secondary dome and on a bedrock massif where ice thickness is only ∼ 2700 m.It is located only ∼ 40 km away from the Concordia station in south-westerly direction.The second area, NP, is 10-30 km north-east of Concordia in the direction of the coast.These "oldest-ice" candidates will be subject to further field studies to verify their suitability.A 3-D model approach would be necessary to study the effect of horizontal advection.Using the shape of the isochrones, which is better constrained than their absolute age, would shed more light on this problem.The possibility of a layer of stagnant ice should also be investigated.Ultimately, in situ study of the age of the bottommost ice at these sites will soon be feasible at minimal operational costs using new rapid access drilling technologies (Chappellaz et al., 2012;Schwander et al., 2014), which will provide in situ measurements to further assess the age of the basal ice and the integrity of the ice core stratigraphy.If successful, this next step will open an exciting opportunity for expanding the EDC records ∼ 700 kyr further back in time, which could help reveal the mechanisms controlling the last major climate reorganisation across the MPT.

Figure 1 .
Figure1.Radar transects used in this study (dotted blue and red lines).The light colour scale represents the bedrock elevation(Fretwell et al., 2013) while the thin grey transparent lines represent the surface elevation(Fretwell et al., 2013).The red square in the inset show the location of the zoomed map around EDC.The red star is the location of the EDC drilling site.The orange squared areas are oldest-ice candidates from Van Liefferinge and Pattyn (2013).The red dotted line is the OIA/JKB2n/X45 radar line displayed in Fig.3.

Figure 2
Figure 2. R(t) proportionality factor applied to accumulation and melting rates (see Eq. 2).The plot is cut at 1 Myr for better readability.R(t) is based on the accumulation record at EDC for the last 800 kyr(Bazin et al., 2013;Veres et al., 2013).

Figure 3 .
Figure 3. One-dimensional ice flow simulation along the OIA/JKB2n/X45 radar transect (see red dotted line in Fig. 1 for location).(a) Various inferred parameters (plain lines) as well as their 15th and 85th percentiles (dashed lines).From top to bottom panels: average surface accumulation rate, geothermal heat flux, p + 1 parameter of the velocity profile, average basal melting, bottom age 60 m above bedrock, height above bed of the 1.5 Myr isochrone and resolution of the 1.5 Myr isochrone.(b) Modelled age (in colour scale; white is for ages older than 1.5 Myr), together with observed isochrones (in white) and bed (in thick black).Note the two main oldest-ice candidates at distance 25 km (North Patch, NP) and at distance 75 km (Little Dome C Patch, LDCP).

F
.Parrenin et al.:  Is there 1.5-million-year-old ice near Dome C, Antarctica?

Figure 4 .
Figure 4. Various bottom-age-related variables along the radar transects, in vivid colours.The bedrock and surface elevations (greyscale and isolines respectively) are shown as in Fig. 1.LDCP and NP are the two old ice patches that we discuss in this study.(a) Modelled bottom age at 60 m above bedrock.(b) Minimum bottom age at 60 m above bedrock with 85 % confidence.(c) Height above bed of the 1.5 Myr isochrone.(d) Temporal resolution for the 1.5 Myr modelled isochrone.

Figure 5 .
Figure 5. Various variables reconstructed by the inverse method along the radar transects, in vivid colour scale.The bedrock and surface elevations (greyscale and isolines respectively) are shown as in Fig. 1.(a) Modelled temporally averaged basal melting.(b) Inferred temporally averaged surface accumulation rate.(c) Inferred geothermal flux.(d) Inferred p vertical velocity parameter.

Table 1 .
Age and total age uncertainty of the 18 isochrones used in this study.