Chronostratigraphy of the Larsen blue-ice area in northern Victoria Land, East Antarctica, and its implications for paleoclimate

In blue-ice areas (BIAs), deep ice is directly exposed at the surface, allowing for the cost-effective collection of large-sized old-ice samples. However, chronostratigraphic studies on blue-ice areas are challenging owing to fold and fault structures. Here, we report on a surface transect of ice with an undisturbed horizontal stratigraphy from the Larsen BIA, northern Victoria Land, East Antarctica. Ice layers defined by dust bands and ground-penetrating radar (GPR) surveys indicate a monotonic increase in age along the ice flow direction on the downstream side, while the upstream ice exhibits a potential repetition of ages on scales of tens of meters, which result from a complicated fold structure. Stable water isotopes (δ¹⁸O_ice and δ²H_ice) and components of the occluded air (i.e., CO₂, N₂O, CH₄, δ¹⁵N–N₂, δ¹⁸O_atm (=δ¹⁸O-O₂), δ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, δ$\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$​​​​​​​, ⁸¹Kr, and ⁸⁵Kr) are analyzed for surface ice and shallow ice core samples. Correlating δ¹⁸O_ice, δ¹⁸O_atm, and CH₄ records from the Larsen BIA with ice from previously drilled ice cores indicates that the gas age at various shallow vertical coring sites ranges between 9.2–23.4 kyr BP, while the ice age sampled from the surface ranges from 5.6 to 24.7 kyr BP. Absolute radiometric ⁸¹Kr dating for the two vertical cores confirms ages within acceptable levels of analytical uncertainty. A tentative climate reconstruction suggests a large deglacial warming of 15 ± 5 ^∘C (1σ) and an increase in snow accumulation by a factor of 1.7–4.6 (from 24.3 to 10.6 kyr BP). Our study demonstrates that BIAs in northern Victoria Land may help to obtain high-quality records for paleoclimate and atmospheric greenhouse gas compositions through the last deglaciation, although in general climatic interpretation is complicated by the need for upstream flow corrections, evidence for strong surface sublimation during the last glacial period, and potential errors in the estimated gas age–ice age difference.

1 Introduction

Ice cores serve as very useful archives for paleoclimate and can provide information on, for instance, the relationship between greenhouse gas concentration and Earth's climate (Petit et al., 1999; EPICA Community Members, 2004; Siegenthaler et al., 2005; Lüthi et al., 2008; Bereiter et al., 2015), abrupt climate changes during the last glacial period (Dansgaard et al., 1989; Steffensen et al., 2008), and bipolar seesaw climate links on millennium timescales (Blunier and Brook, 2001; Landais et al., 2015). To date, continuous climate records from ice cores cover the last 800 kyr, while new deep-drilling projects in Antarctica now aim at recovering ice that is more than 1 Myr old (Fischer et al., 2013). However, the use of ice cores obtained from conventional deep-ice-core-drilling projects has some limitations, as many analyses require large amounts of ice, such as trace element isotope and trace gas analyses. In addition, deep-drilling projects incur high economic and labor costs. In contrast, coring in blue-ice areas (BIAs) has emerged as an alternative to obtain large ice samples in a cost-effective manner (Folco et al., 2006; Petrenko et al., 2006; Schaefer et al., 2006; Sinisalo et al., 2007; Korotkikh et al., 2011; Turney et al., 2013; Bauska et al., 2016; Aarons et al., 2017; Baggenstos et al., 2017; Zekollari et al., 2019; Yan et al., 2019; Fogwill et al., 2020).

Over the vast majority of the Antarctic ice sheet surface, ice forms through the compaction of snow. This ice then deforms and, through gravity, flows towards the margin of the ice sheet (Gardner et al., 2018). However, (basal) topographic obstacles in some cases redirect the ice flow towards the surface in some so-called blue-ice areas. In these BIAs, annually more surface snow is ablated than accumulated by katabatic winds and/or sublimation (Bintanja, 1999; Sinisalo and Moore, 2010). Because ice layers of the same age (isochrones) are extended on the surface, large amounts of old ice of a specific age can be obtained at the surface and/or relatively shallow depths, allowing for studying the paleoclimate without constraints related to limited sample sizes (Schaefer et al., 2006; Bauska et al., 2016; Fogwill et al., 2020).

In some BIAs, when ice flow and surface conditions combine favorably, large concentrations of meteorites can accumulate on the surface of the BIAs (Whillans and Cassidy, 1983; Cassidy et al., 1992; Harvey, 2003; Tollenaar et al., 2022). Recently, studies on BIAs have drawn particular attention to the possibility of identifying ice older than 800 kyr BP (to date, the longest continuous ice core record covers the last 800 kyr), which is of particular interest because glacial–interglacial cycles changed from 40 to 100 kyr around this time period, which is referred to as the Mid-Pleistocene Transition (MPT). It has been reported that ice at the Allan Hills BIA has ages of 90–250 kyr BP on the surface (Spaulding et al., 2013), reaching ∼ 2.7 Myr BP near the bedrock at depths of about 150–200 m (Yan et al., 2019). However, the use of ice samples taken from BIAs (hereafter referred to as blue ice) has several drawbacks. In most cases, the stratigraphy of the BIAs is complicated, as shown by the fold and fault structures on the surface (Folco et al., 2006; Petrenko et al., 2006; Curzio et al., 2008; Schaefer et al., 2009; Baggenstos et al., 2017), and the stratigraphy is discontinuous at deep depths near the bedrock where ice ages are similar to or older than the MPT (Higgins et al., 2015; Yan et al., 2019).

BIAs cover about 1.7 % of the Antarctic continent, of which a large fraction is concentrated in Victoria Land, the Transantarctic Mountains, Dronning Maud Land, and the Lambert Glacier basin (Hui et al., 2014) (Fig. 1). Several chronological studies of BIAs have been conducted in southern Victoria Land: on Taylor Glacier (Aciego et al., 2007; Buizert et al., 2014; Baggenstos et al., 2017, 2018; Menking et al., 2019), in the Allan Hills (Spaulding et al., 2013; Higgins et al., 2015; Yan et al., 2019, 2021), on Mullins Glacier (Yau et al., 2015), and in northern Victoria Land on Frontier Mountain blue-ice field (Folco et al., 2006; Curzio et al., 2008; Welten et al., 2008). Studies on Taylor Glacier and Allan Hills have constrained age along several transects and cores. However, chronologies of the other BIAs in Victoria Land remain largely unknown, prohibiting high-resolution paleoclimate studies. In addition, important paleoclimate proxies such as stable water isotopes, greenhouse gases, and isotopic ratios of oxygen gas have, to date, not been analyzed for the BIAs in northern Victoria Land.

Generally, because of the complicated stratigraphy, blue-ice samples are dated by combining multiple methods (Sinisalo and Moore, 2010). In previous studies, ice flow modeling, correlation analysis with stable isotopes of ice, and radiometric analysis of meteorites and tephra in ice have been used, but in most cases the age constraints have not been sufficiently precise to study the paleoclimate (Azuma et al., 1985; Nakawo et al., 1988; Reeh et al., 2002; Folco et al., 2006; Aciego et al., 2007; Curzio et al., 2008; Dunbar et al., 2008). One effective method for dating the gas age is to correlate globally well-mixed atmospheric gas records (e.g., CH₄, CO₂, and δ¹⁸O_atm) with existing well-dated ice core records (Spaulding et al., 2013; Baggenstos et al., 2017; Menking et al., 2019; Yan et al., 2021). Other methods include the use of stable Ar isotopes (Higgins et al., 2015; Yau et al., 2015; Yan et al., 2019) or radioactive ⁸¹Kr (Loosli and Oeschger 1969; Buizert et al., 2014; Tian et al., 2019; Crotti et al., 2021), both of which provide independent and absolute age constraints. However, the use of Ar and Kr isotopes has a certain limitation that allows for an accurate age constraint. The age uncertainty in Ar dating is ±180 kyr or 11 % of the age according to a regression line used by Bender et al. (2008). In addition, the ages can be corrupted by the injection of radiogenic ⁴⁰Ar from the continental crust (Bender et al., 2010). ⁸¹Kr is a cosmogenically produced radioactive isotope with a half-life of 229 ± 11 kyr, decaying to ⁸¹Br via electron capture. It is mixed in the atmosphere within 1–2 years and has no other significant sources or sinks, making ⁸¹Kr an ideal tracer (Oeschger, 1987; Zappala et al., 2020). Oeschger (1987) suggested the potential of ⁸¹Kr for radiometric ice core dating. However, at that time, 10⁵–10⁶ kg of ice was required. Owing to the development of Atom Trace Trap Analysis (ATTA), the required ice has continued to decrease (Lu et al., 2014; Tian et al., 2019; Jiang et al., 2020; Crotti et al., 2021). Buizert et al. (2014), for the first time, showed that ⁸¹Kr age dating is feasible for blue ice at Taylor Glacier, Antarctica. However, the age uncertainty in ⁸¹Kr dating ranges between 5 %–20 % of the age depending on the sample age and sample size (Jiang et al., 2020). It also has a systematic age uncertainty of ∼ 5 % due to the uncertainty in the ⁸¹Kr half-life. Given the uncertainties in using Ar and Kr isotopes for direct dating, here we opt to use correlation of atmospheric gas records with existing well-dated ice core records as the primary dating tool.

CH₄ and O₂ are well mixed in the atmosphere (Blunier et al., 2007), and conducting correlations of both CH₄ concentration and δ¹⁸O of O₂ (δ¹⁸O_atm) is a well-known strategy for establishing chronologies of blue ice (Petrenko et al., 2006; Baggenstos et al., 2017; Yan et al., 2021). Variations in continental ice mass are known to be the main factor that controls δ¹⁸O_atm during glacial–interglacial cycles (Bender et al., 1985; Sowers et al., 1993). In addition, variations in the Intertropical Convergence Zone (ITCZ) are also known to influence δ¹⁸O_atm on millennial and orbital timescales (Severinghaus et al., 2009; Landais et al., 2010; Seltzer et al., 2017; Extier et al., 2018b). Because of the long lifetime of oxygen gas in the atmosphere (∼ 1 kyr), δ¹⁸O_atm varies more gradually compared to CH₄, limiting synchronization to millennial timescales. In contrast, atmospheric CH₄ concentration changes rapidly because of the short lifetime of ∼ 12 years, allowing for precise dating via stratigraphic matching.

For dating ice ages, glaciochemical records (e.g., nss-Ca²⁺ (nss signifies non-sea-salt), δ¹⁸O_ice, and δ²H_ice) can be used for correlation with existing well-dated ice core records (Baggenstos et al., 2018; Menking et al., 2019). Notably, at a given depth, the ice is older than the gas because the gas is only isolated and stops mixing with the atmospheric air when the firn completely transforms into ice (Schwander and Stauffer, 1984).

Figure 1 Distribution map of BIAs in Antarctica. BIAs are noted in dark blue (Hui et al., 2014). Orange dots represent the BIAs where the chronology has been studied. Deep-ice-core locations with brown dots are labeled as follows: TALos Dome Ice CorE (TALDICE), Taylor Dome (TD), EPICA (European Project for Ice Coring in Antarctica) Dome C (EDC), West Antarctic Ice Sheet (WAIS) Divide (WD), and Vostok. To keep a consistent orientation with the GPR profile in Fig. 3, we flipped the classical map of Antarctica (East Antarctica to the left-hand side). The Antarctic map was obtained from the QGIS (geographic information system) Quantarctica package.

Our study focuses on the chronostratigraphy of ice in the Larsen BIA, with the goal of facilitating future research in this region. We describe the local ice flow and structure of the ice layers using dust bands and ground-penetrating radar (GPR) surveys and assess the alterations of the measured stable water isotopes, greenhouse gases (CH₄ and CO₂), and gas isotopes (δ¹⁵N–N₂ and δ¹⁸O_atm). To constrain the unknown gas and ice ages, δ¹⁸O_atm, CH₄, and δ¹⁸O_ice are correlated with existing ice core records. We also independently confirm the ages using the radiometric ⁸¹Kr dating method. Finally, using δ¹⁵N–N₂ and the Δage (ice and gas age difference) results, we present a record of surface temperature and accumulation rate. For this, in contrast to previous studies, which used the Herron–Langway model (Herron and Langway, 1980; Baggenstos et al., 2018; Menking et al., 2019; Yan et al., 2021), we apply a recently developed analytical framework which does not require stable water isotope values to estimate past surface temperatures and accumulation rates (Buizert, 2021).

2 Study area and methods

2.1 Larsen blue-ice area (BIA)

We sampled ice from the Larsen BIA, an outlet glacier in northern Victoria Land, East Antarctica, during the austral summer of 2018/19. The Larsen BIA is located approximately 85 km southwest of the Korean Jang Bogo Station (Fig. 2a). The mean annual temperature is −24.4 ± 11.7 ^∘C (1σ) (for the year 2016; data and information were obtained from the Meteo-Climatological Observatory at MZS – Mario Zucchelli Station – and Victoria Land of PNRA – Programma Nazionale di Ricerche in Antartide; http://www.climantartide.it, last access: 1 February 2019), cold enough to prevent ice melting in the summer. Consequently, the ablation of ice at the surface occurs through sublimation. There were ∼ 20 cm wide dust bands with gentle folding structures in the mid- to downstream part, while we observed severely folded dust bands (e.g., S and Z folds) in the upstream part (Fig. 2b). Dust bands are frequently observed in BIAs in Victoria Land and can be used as isochrons (Sinisalo and Moore, 2010). To obtain ice samples with a simple stratigraphy, we avoided ice coring in areas with complex fold structures and sampled ice in the direction of ice flow identified by the Antarctic ice velocity map in the QGIS Quantarctica package (Rignot et al., 2011; Mouginot et al., 2012; Matsuoka et al., 2018, 2021). In contrast, the stratigraphy of the upstream ice, where we collected surface ice, may be inverted and repeated on a scale of tens of meters. Shallow ice cores were drilled along a 1 km long transect at intervals of 20–30 m (Fig. 2b). Most ice cores had lengths of approximately 2 m, but ice core no. 23 (74.9319^∘ S, 161.6018^∘ E) was ∼ 10.4 m (Fig. 2c). We also collected near-surface ice samples (∼ 500 g) along a 1.3 km transect at 20 m intervals to measure stable water isotopes of the ice (δ¹⁸O_ice and δ²H_ice) at depths of ∼ 5–10 cm (hereafter referred to as the surface ice sample). The average horizontal spacing of the surface ice samples and ice cores was approximately 10.3 m. The ice core TF (Tarn Flat) (74.93042^∘ S, 161.56975^∘ E) with a length of ∼ 12 m on the upstream side was sampled in 2016 for a preliminary study, as described by Jang et al. (2017). An imaginary line parallel to the ice flow direction was used to define the horizontal distance, while a perpendicular line from each sampling location to the line parallel to the ice flow direction was projected to identify the intersection point. Each intersection point was then used to measure the horizontal distance from the most upstream sampling site (Fig. S1).

Figure 2 Location of the Larsen BIA and sample collection. (a) Location of the Larsen BIA, EM-D (Elephant Moraine) core, and Jang Bogo Station. (b) Magnified map of the Larsen BIA, including sample locations. Orange lines represent dust bands observed in the blue-ice field. The line marks are derived from Google Earth imagery (https://earth.google.com/web/@-74.93,161.58,1080.8673325a,2487.72294087d,49.37163057y,347.89929343h,0t,0r, last access: 2 May 2022). Blue dots represent locations of surface ice samples. Red diamonds are locations of shallow ice cores. Six representative names are shown: red letters are for ∼ 10 m long cores (TF and no. 23), and black letters are for 2 m long cores (no. 306, no. 120, and no. 200) and surface ice (81w). Z and S folds are identified by the dust bands at upstream sites. The total length of the ice sample transect is approximately 1.4 km. (c) Schematized cross-section of the transect. Satellite photo of Antarctica is from the QGIS Quantarctica package (Bindschadler et al., 2008).

2.2 Ground-penetrating radar (GPR) survey

Approximately 17 km of GPR data was collected during 2 d in January 2019 (Fig. S1). A MALÅ ProEx (Professional Explorer) impulse radar system with a 50 MHz unshielded antenna for a larger penetration depth was used for data acquisition. Records were obtained at a sampling frequency of 559.5 MHz, with a time interval 0.1 s, and stacked four times. The survey was conducted at a speed of 2.5–3 km h⁻¹ to minimize noise caused by frictional vibrations between the antenna and the surface of the glacier. The position of traces was recorded using a single-frequency code-phase GPS with an accuracy of <3 m. Data processing was performed in Reflexw v9.5 in the following order: dewow, DC (direct current) filter, band-pass filter, time-zero drift, energy decay correction, background removal, static correction, migration, and stack. As the constant radar signal velocity v in the glacier was 0.17 m ns⁻¹ (Borgorodsky et al., 1985​​​​​​​; Reynolds, 1985) and the frequency f of the radar system was 50 MHz, the GPR wavelength λ was about 3.4 m ($\mathit{\lambda }=v{f}^{-\mathrm{1}}$).

2.3 Stable water isotope measurement

Stable water isotopes (δ¹⁸O_ice and δ²H_ice) were measured simultaneously at the Korea Polar Research Institute (KOPRI) using the cavity ring-down spectroscopy (CRDS) method using a Picarro L2130-i with a vaporizer (A0211). Surface ice samples (∼ 5–10 cm depth) and ice core samples (10–30 and 190–200 cm depths of ice cores) were used for the analyses. Ice from the ∼ 10 m long core no. 23 was measured at 20 cm vertical intervals. Ice samples were kept in Whirl-Pak bags and melted at room temperature. The melted sample was then injected into a 2 mL vial via a disposable syringe with a 0.45 µm filter. One batch of measurements consisted of 5 duplicate samples (10 samples in total) and a working standard. The first and second aliquots of the sample were measured 12 times each, while the working standard was measured 20 times. The last six measurements were used to remove memory effects from the previous samples. The precision was evaluated by measuring the working standard repeatedly (n=67); 1σ (standard deviation) was 0.07 ‰ for δ¹⁸O_ice and 0.90 ‰ for δ²H_ice. The working standards and samples were calibrated against the international standards for water isotopes, VSMOW2 (Vienna Standard Mean Ocean Water 2), SLAP2 (Standard Light Antarctic Precipitation 2), and GISP (Greenland Ice Sheet Precipitation).

2.4 Greenhouse gas (GHG) measurement

The CH₄ concentrations in the ice cores (no. 306, no. 23, no. 120, and no. 201) were analyzed at Seoul National University (SNU) using a melt–refreeze technique with a flame ionization detector gas chromatograph (FID-GC) at 10 cm vertical intervals. The CH₄ measurement process is described in detail by Yang (2019). Briefly, we trimmed the outermost side and the cracks of the ice samples by approximately 2 mm with a clean band saw to eliminate contamination by ambient air. The mass of the ice used for measurement was about 35–55 g. Then, we placed the ice inside a custom-made flask and evacuated the flask. We used a 740.6 ppb standard air sample from the National Oceanic and Atmospheric Administration (NOAA) Global Monitoring Division (GMD) on the WMO (World Meteorological Organization) X2004A scale to establish a daily calibration line (Dlugokencky et al., 2005). The NOAA standard air was injected into four flasks containing bubble-free ice samples, which served as the control group. The CH₄ mixing ratio of the control group was measured, and the daily average offset to the standard value of 740.6 ppb was considered a systematic error; thus, the daily average offset (average of the daily average offsets: 11.4 ± 3.2 ppb (1σ); average of the intra-day standard deviations of the control group: 3.7 ± 1.2 ppb (1σ)) was subtracted from each measurement result. The trapped gas in the ice was liberated by immersing the flask in a hot water and then refreezing the meltwater. In addition, as CH₄ is more soluble than the major air components (N₂, O₂, and Ar), the CH₄ concentration of the initially extracted air is lower than the original value. Therefore, a second gas extraction was conducted (refrozen meltwater was melted and refrozen once again) to correct the solubility effect.

The same ice cores (no. 306, no. 23, no. 120, and no. 201) were used for CO₂ measurements. Generally, CO₂ should be measured using a dry-extraction system rather than wet extraction because of its high solubility in water and the possibility of CO₂ production by carbonate–acid reactions in ice melt (Delmas et al., 1980). Therefore, for our measurement, air was collected using a needle-crusher dry-extraction instrument and detected by FID-GC at SNU. The ice was pre-treated in the same manner as that for CH₄, as described above, while the required amount of ice was 15–20 g. We placed the ice in the needle-crusher chamber, which was cooled to −35 ^∘C during ice preparation. The CO₂ results could be affected by adsorption and desorption by the instrument and the sample tube (Ahn et al., 2009). Therefore, before crushing the ice, NOAA standard air was released through the chamber and collected into the sample tube. The CO₂ mixing ratio was measured, and the average offset with respect to NOAA standard air (systematic error) was subtracted from each measurement results. The CO₂ mixing ratio of the standard air we used for making a daily calibration line was 285.66 ppm or 293.32 ppm (because gas in the tank had been exhausted), which is from NOAA GMD on the WMO X2019 calibration scale (Hall et al., 2021). Trapped air in the ice was collected in a sample tube cooled by a helium closed-cycle refrigerator (He-CCR) (∼ 12 K) after crossing the water trap (approximately −80 ^∘C). The entire process is described in detail by Shin (2014).

We also used ice samples from 35 cores at depths of 190–200 cm (hereafter referred to as horizontal measurement) and ice core no. 23 to measure CH₄, CO₂, and N₂O concentrations, along with other gas isotopes (see below), by a wet-extraction method at the National Institute of Polar Research (NIPR) in Japan. For CO₂ and CH₄, FID-GC was used, and for N₂O, an electron capture detector (ECD) GC was used. To determine the concentrations of CH₄ and CO₂, a calibration line was established using TU-2008 (Tohoku University) scale standard air; for N₂O, TU-2006 scale standard air was used. Details on the standard air used in this study (named STD 1, 3, and 5) are described by Oyabu et al. (2020). Three standard airs and a quadratic calibration line were used to determine the GHG concentration at NIPR, while a single standard air and a linear calibration line were used to determine the GHG concentration at SNU. Different calibration scales of standard air (TU and NOAA–WMO) and different calibration methods may contribute to an inter-laboratory offset of GHG concentrations. CO₂ was measured again at SNU using the dry-extraction method at 190–195 cm depth for the horizontal measurement (Table S7). The ice preparation process is briefly explained in Sect. 2.5.

Here, we also report CH₄ and CO₂ measurement results of the EM-D core for comparison (76.253^∘ S, 156.562^∘ E) (Fig. 2), which was sampled on 13 December 2016 at Elephant Moraine Texas Bowl for a preliminary study (Jang et al., 2017).

2.5 Analyses of δ¹⁵N–N₂, δ¹⁸O_atm, δ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, and δ$\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$

O₂, N₂, and Ar isotopes were measured at NIPR using a dual-inlet mass spectrometer (Thermo Fisher Scientific Delta V). The isotope results were evaluated with respect to the modern atmosphere. As mentioned previously, we measured 35 ice cores at depths of 190–200 cm to evaluate how ancient-air compositions change horizontally and also measured gas isotopes using ice core no. 23 at several depths with a range of 0–10 m to compare with the results from the neighboring cores at depths of 190–200 cm. Because gas loss fractionation by molecular diffusion during storage could have affected the isotope results (Ikeda-Fukazawa et al., 2005; Oyabu et al., 2021), we trimmed the surface and the cracks of the ice approximately 3–5 mm and then shaved away some blurry ice surface using a ceramic knife. The weight of the ice obtained was 55–80 g after trimming the surfaces. The ice inside the vessel was evacuated for ∼ 120 min. To extract the gas, each vessel was gradually immersed in a hot-water container. Simultaneously, the released air was collected into a sample tube, which was cooled by an He-CCR. After homogenizing the sample tube overnight, the gas was split into two aliquots, one for isotope analysis (1.5 mL) and the other for greenhouse gas measurement (5 mL). The measurement process is described in more detail by Oyabu et al. (2020).

Obtaining a true δ¹⁸O value in the past atmosphere from ice cores requires gravitational, thermal, and gas loss fractionation corrections. The gravitational factor is proportional to the difference in the mass number between isotopes (Craig et al., 1988; Severinghaus et al., 1998). Therefore, δ¹⁸O (${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$) was affected twice as much as δ¹⁵N (${}^{\mathrm{15}}\mathrm{N}{/}^{\mathrm{14}}\mathrm{N}$) by gravity. Hence, each gas isotope was gravity-corrected using Eq. (1):

$$\begin{array}{}\text{(1)}& {{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{atm}}={\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{2}\times {\mathit{\delta }}^{\mathrm{15}}\mathrm{N}.\end{array}$$

Following the same principle, δ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ is more affected than δ¹⁵N, by a factor of 4; thus it was corrected using Eq. (2) and used to assess the gas loss fractionation:

$$\begin{array}{}\text{(2)}& {\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/\mathrm{N}}_{\mathrm{2},\mathrm{gravcorr}}=\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}-\mathrm{4}\times {\mathit{\delta }}^{\mathrm{15}}\mathrm{N}.\end{array}$$

Along with gravitational correction, thermal fractionation should also be considered because temperature gradients in the firn column affect the distribution of the isotopes (Severinghaus et al., 1998; Goujon et al., 2003). However, thermal fractionation is typically small in Antarctic ice cores due to the relatively gradual nature of surface climate change (Goujon et al., 2003). Comparing ${\mathit{\delta }}^{\mathrm{40}}\mathrm{Ar}/\mathrm{4}$ (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$) and δ¹⁵N together allows for the discrimination of the contribution of thermal and gravity fractionation (Severinghaus et al., 1998). However, for our samples, thermal fractionation could not be considered because ³⁶Ar interfered with ³⁶O₂.

Following Capron et al. (2010), we did not correct for gas loss fractionation in our samples because the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/\mathrm{N}$_2,gravcorr values were significantly greater than −30 ‰, except for one measurement for sample no. 301 (Tables S4 and S5). $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/\mathrm{N}$_2,gravcorr of around −30 ‰ indicates that the ice is poorly preserved and has experienced considerable gas loss during either drilling or storage (Landais et al., 2003). We did not use the results of no. 301 when constraining the gas age using δ¹⁸O_atm. Gravity-corrected $\mathit{\delta }\mathrm{Ar}/\mathrm{N}$_2,gravcorr ($\mathit{\delta }\mathrm{Ar}/\mathrm{N}$${}_{\mathrm{2},\mathrm{gravcorr}}=\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}-\mathrm{12}\times {\mathit{\delta }}^{\mathrm{15}}\mathrm{N}$) is listed in Tables S4 and S5 but was not used in this study.

2.6 ⁸¹Kr dating

For the ⁸¹Kr measurement of ice core no. 23 and TF, 5.3 kg (depth: 711–1040 cm) and 5.4 kg (depth: 798–1192.5 cm) of ice were used, respectively. Air for the ⁸¹Kr measurement was extracted at SNU using the instrument provided by the University of Science and Technology of China (USTC). The ice was kept in a tank with an O-ring lid and evacuated using a dry scroll pump with a water trap, then melted by immersing the tank in hot water. The released gas was collected in containers and shipped to the USTC for Kr purification and ⁸¹Kr analysis using Atom Trap Trace Analysis (ATTA). The extraction procedure was described in detail by Tian et al. (2019), and the principle of ATTA is described by Jiang et al. (2012).

The anthropogenic ⁸⁵Kr was measured simultaneously with ⁸¹Kr to quantify any contamination with modern air. For both samples, the measured ⁸⁵Kr activity was below the detection limit; hence no correction for contamination with modern air was necessary. For the calculation of the ⁸¹Kr ages, the changes in the past atmospheric ⁸¹Kr abundance due to variations in the cosmic ray flux on Earth (Zappala et al., 2020) were considered. The age uncertainty calculation was based on the statistical error in atom counting.

2.7 Development of the WD2014 timescale for TALDICE

When constraining the ice age of the Larsen ice, the TALDICE record was selected to synchronize Larsen δ¹⁸O_ice because of its proximity to the Larsen BIA. For the TALDICE ice core, we used a timescale that was synchronized to the WAIS Divide (WD) WD2014 chronology (Buizert et al., 2015; Sigl et al., 2016) in the following way: the ice age scales were synchronized using volcanic deposits identified in sulfur/sulfate data from the WD and TALDICE cores (Severi et al., 2012; Buizert et al., 2018). Next, the TALDICE gas age–ice age difference (Δage) was established empirically by matching abrupt changes in atmospheric CH₄ to the WD core; at each match point, this provided one discrete Δage constraint. A dynamical firn densification model based on Herron–Langway firn physics (Herron and Langway, 1980) was then used to interpolate between these empirical Δage constraints in order to obtain a TALDICE gas age scale for all depths (Buizert et al., 2021).

3 Results

3.1 Ground-penetrating radar (GPR) survey

In the GPR survey, we identified ice layers (or isochrones) in the transect parallel to the ice flow direction (Fig. 3). The dips of the ice layers range from 1 to 6^∘ with a decreasing trend from the upstream to the downstream direction. The ice layers of the radargram were not clearly visible at a depth of <10 m because of the direct wave signal. We did not observe any stratigraphic folding structure in the ice layers (and associated age inversion) along the ice flow direction in the mid- to downstream areas (from ice cores no. 23 to no. 200). Therefore, we expect monotonic and continuous age changes along the ice flow direction. However, as shown in the dust bands with S and Z folds in the upstream area (Fig. 2b), the upstream stratigraphy might be repeated on a scale of tens of meters (from 81w to ice core no. 23). In addition, the subsurface ice layer in the upstream area (0–800 m from the most upstream side) was not well recognized from the GPR profile (Fig. 3c). It is possible that the noise caused by (subsurface) crevasses, cavities, or cracks obscured the signals. In addition, accurate data acquisition might have been hindered by antenna tremors or low battery power at severely cold temperatures. The basal topography is well defined from the GPR data; we observed an ice thickness variation of 200–380 m (Fig. 3a and b). The results of the bedrock elevation and ice thickness (Fig. 3a) were obtained using a kriging method. The ice thickness decreases as ice flows with increasing bedrock elevation, which is a favorable condition for old ice to be exposed at the surface.

Figure 3 Ground-penetrating radar (GPR) survey overview. (a) Bedrock elevation and ice thickness of the Larsen BIA. (b) GPR profile of the transect (A to A^′ in a) through the blue-ice field. Identifiable ice layers are indicated in orange lines. Ice layers are not well identified at a distance of <800 m. (c) Enlarged upper 80 m of the GPR profile. Ice layers at a distance of >2500 m have a relatively flat dip. Vertical axis is exaggerated, and the dip of the ice ranges between 1 and 6^∘. The average dips used for estimating the ice ages at 1.95 m depth are indicated by yellow boxes.

3.2 Stable water isotopes

Stable water isotopes in ice record surface temperatures of the past at snow deposition sites (Jouzel et al., 1997). The δ¹⁸O_ice and δ²H_ice records of the Larsen BIA are presented in Tables S8, S9, and S10. The horizontal δ²H_ice has a distinct local minimum around a horizontal distance of 600–800 m, indicating a transient cold event (Fig. 4b). The δ²H_ice of the transient cold event plunges by approximately 70 ‰, heading to the upstream ice. The most downstream sample, ice core no. 200, has a δ²H_ice of −369 ‰. Conversely, the most upstream sample, surface ice 81w, has a δ²H_ice of −245 ‰. Ice core no. 23, located in the middle of the transect, and the vertical δ²H_ice profile range between −353 ‰ and −291 ‰. The water isotope values from the Larsen BIA are scattered over a wider range than most other published ice core records (Petit et al., 1999; EPICA Community Members, 2004; Stenni et al., 2011). A highly variable δ¹⁸O_ice has also been reported at Taylor Glacier (Baggenstos et al., 2018; Menking et al., 2019). Severe scattering might indicate that the accumulation zone (original source of ice) of the Larsen BIA might have experienced more variability in temperature and/or the vapor source (variability in atmospheric conditions) than other sites. The wider range of water isotope values during the deglaciation as compared to other published ice core records might be due to a higher increase in temperature and/or because of the spatial difference between the original deposition site of upstream and downstream ices. Deuterium excess (d=δ²H${}_{\mathrm{ice}}-\mathrm{8}\times {\mathit{\delta }}^{\mathrm{18}}$O_ice) shows a wide range (8.44 to −8.69 ‰, Table S8) from the entirety of the near-surface ice samples. The negative d excess likely indicates that isotopic fractionation is attributed to the sublimation of ice in the accumulation zone (Hu et al., 2021). Negative d-excess values were also observed in the Allan Hills BIA (Hu et al., 2021). Meanwhile, sublimation can deplete ¹⁶O and ¹H in ice, thereby enriching isotopic ratios (δ¹⁸O_ice and δ²H_ice). Thus, a negative value of d-excess results indicates that stable water isotopes are not proper proxies for the changes in temperature and vapor sources.

To match the two δ²H_ice profiles (measurement of ice core no. 23 and horizontal measurement of near-surface ice), the depth of ice core no. 23 was converted to the horizontal distances by pinpointing the deepest result of no. 23 to the result of ice core no. 104. The horizontal distances for the rest of core no. 23 were assigned through linear interpolation. The two δ²H_ice profiles match well, as indicated by an r² value of 0.86 (p<0.001) (Fig. 4b). The similarity of δ²H_ice between the vertical and horizontal measurements demonstrates that the ice stratigraphy was not disturbed. The average difference in the values of δ²H_ice between vertical ice core no. 23 and the horizontal record is 9.8 ‰ and 0.6 ‰ for δ¹⁸O_ice. We assume that the average difference is the uncertainty (1σ) of stable water isotope values of the Larsen ice. The ∼ 10 m long no. 23 core covers ages that correspond to a surface distance of ∼ 117 m, from which we calculate an average dip of the ice layer of 4.96^∘, which is comparable with the average dip derived from the GPR profile (Fig. 3c). This depth–distance relationship is also supported by the comparison of gas isotope values (δ¹⁸O_atm and δ¹⁵N–N₂) from ice core no. 23 to the records from the horizontal measurement at a depth of 1.95 m (Fig. A2).

Figure 4 δ¹⁸O_atm and δ²H_ice records from the Larsen BIA. (a) δ¹⁸O_atm from each ice core at a 1.95 m depth (horizontal measurement). The δ¹⁸O_atm of the horizontal measurement was gravity-corrected using Eq. (1). (b) δ²H_ice from near-surface ice (∼ 5–10 cm depth surface ice samples and 10–30 cm depth ice core samples) and ice core no. 23. The surface ice record was matched with core no. 23 using two tie points (orange dots), with a corresponding r² value of 0.86 (p<0.001).

3.3 Analysis of gas entrapped in the ice

3.3.1 CH₄ and CO₂ mixing ratios

The measurement results of the CH₄ and CO₂ concentrations for the vertical cores are presented in Fig. A1. In the shallow ice cores (no. 306, no. 120, and no. 201) greenhouse gases are significantly altered for the top 2 m, showing out-of-range values of the range of natural greenhouse gas concentrations during the last 800 kyr: 340–800 ppb for CH₄ and 180–300 ppm for CO₂. The CH₄ records of no. 23 also fluctuate significantly at 0–4.6 m depth, which is reduced at >4.6 m. CO₂ records of no. 23, in contrast, gradually decrease and become steady at a depth of >4.6 m. A comparison of the results from NIPR and SNU for ice core no. 23 shows that the difference in the concentration decreases significantly at depths of >4 m. There are large differences (30–140 ppb) of CH₄ at 0–4 m depths, but these decreased considerably (5–10 ppb) at a depth of >4 m. The CO₂ difference was 10–20 ppm at <4 m but decreased to 2–10 ppm at a depth of >4 m (Tables S5 and S6). The CO₂ mixing ratio of the TF core, reported by Jang et al. (2017), showed that CO₂ was scattered even in the deeper part (>4 m) of the ice core. We speculate that this scattering is due to complicated ice stratigraphy, since we observed several folding structures identified by dust bands near the TF core site (Fig. 2b). Altered CH₄ in near-surface samples have also been identified at Pakitsoq (Greenland), Taylor Glacier, and Elephant Moraine (Petrenko et al., 2006; Baggenstos et al., 2017). By combining our new measurements with those from the literature, we find the depth at which greenhouse gas concentrations stabilized appears to be inversely proportional to mean annual temperature (Table 1).

Table 1 Depth of unaltered greenhouse gas compositions and mean annual temperature of BIAs.

^a Depth of unaltered greenhouse gas compositions from Elephant Moraine Texas Bowl remains uncertain due to the lack of data at depth of >10 m. Mean annual temperature of Elephant Moraine: provided by KOPRI's automatic weather station (AWS) record of the years 2020 and 2021. ^b Allan Hills: not provided by an AWS but by stable water isotopes; ^c deduced by vertical profile of δ¹⁵N–N₂ and δ¹⁸O_atm values in Allan Hills. ^d Taylor Glacier: assumed to be comparable with the mean annual temperature of nearby McMurdo Station (∼ 100 km away). ^e Pakitsoq (Greenland): assumed to be comparable with the mean annual temperature of the nearby town of Ilulissat (∼ 40 km away).

3.3.2 δ¹⁵N–N₂ and δ¹⁸O_atm

As δ¹⁵N–N₂ and δ¹⁸O_atm were measured at very shallow depths (∼ 1.95 m), we compared the horizontal results with the vertical distribution of ice core no. 23. The depths of ice core no. 23 were converted to horizontal distances at a depth of 1.95 m (Fig. A2) using Eq. (3):

$$\begin{array}{}\text{(3)}& {\mathrm{HD}}_{\mathrm{1.95}\mathrm{m}}=\left[\mathit{\right\{}\left(D_{\mathrm{no}.\mathrm{23}}\right)-\mathrm{1.95}\mathit{\}}/\tan (\mathrm{4.96}{}^{\circ }\left)\right]+\mathrm{663},\end{array}$$

where HD_1.95 m represents the horizontal distance in meters (at 1.95 m depth) and D_no. 23 represents the depth of ice core no. 23. The length of 663 m reflects the distance between ice core no. 23 and the most upstream sampling site. The average dip of the ice layer (4.96^∘) was calculated by matching the δ²H_ice of ice core no. 23 with the horizontal records from the neighboring cores (Fig. 4b). The δ¹⁵N–N₂ and δ¹⁸O_atm values obtained from the horizontal record at 1.95 m are comparable to those from ice core no. 23, except those from the very shallow depths of <0.5 m (Fig. A2), confirming that the gas isotope ratios are generally reliable at a depth of 1.95 m. We estimated the uncertainty in δ¹⁵N–N₂ and δ¹⁸O_atm values as ±0.05 ‰ (1σ) because the differences between no. 23 and the horizontal record were approximately 0.05 ‰. Petrenko et al. (2006) also reported for a site in western Greenland that as long as the ice is not affected by surface ice melting, nitrogen and oxygen gas isotopes are not altered even at a depth of 0.3 to 0.4 m. In contrast, δ¹⁵N–N₂ and δ¹⁸O_atm values in the Allan Hills are relatively stable below 7–10 m (Spaulding et al., 2013). No significant gaps, discontinuities, or anomalies were found within the results from horizontal measurements of δ¹⁸O_atm, indicating a relatively undisturbed stratigraphy (Fig. 4a).

4 Discussion

4.1 Glacial-termination identification in the Larsen ice

The increase in δ²H_ice and decrease in δ¹⁸O_atm from downstream (older ice) to upstream ice (younger ice) reflects warming atmospheric conditions over time (Fig. 4). More specifically, the δ¹⁸O_atm values of the Larsen ice (1.126 ‰ to −0.075 ‰) are typical of a glacial-termination period. Based on satellite imagery of the studied area, the Larsen BIA might have been located directly downstream of an accumulation zone, and the ice might have been originated solely from an accumulation zone because no other BIAs are identified upstream (refer to Fig. 1 in Zekollari et al., 2019). This expectation is supported by the GPR survey profile (Fig. 3b and c) and the δ¹⁸O_atm (Fig. 4a), which indicates no complex ice stratigraphy and no discontinuities in the downstream area of the Larsen BIA, respectively. In terms of this ice flow condition, we expect the age of the Larsen BIA corresponds to a relatively young one among the latest glacial terminations.

The maximum δ¹⁸O_atm value of >1 ‰ implies that the ice is younger than 800 kyr, as the glacial climate conditions became more extreme following the MPT (Lisiecki and Raymo, 2005; Elderfield et al., 2012; Chalk et al., 2017); values from the literature suggest that the maximum δ¹⁸O_atm values were <1.0 ‰ during pre-MPT glacial periods (Yan et al., 2019). The EPICA Dome C (EDC) record shows that Terminations I, II, IV, V, and VII are the only terminations that have both negative δ¹⁸O_atm values and values >1.0 ‰ over the last 800 kyr (Landais et al., 2013; Extier et al., 2018a) (Fig. 5). Hence, we exclude Terminations III and VI for the Larsen BIA ages. Among the candidates, Termination V should also be excluded because the maximum δ¹⁸O_atm value in the EDC record is ∼ 1.4 ‰, which is significantly higher than that of the Larsen ice. The δ²H_ice decrease in the middle of the glacial termination in the Larsen BIA (Fig. 4b) is similar to that during the Antarctic Cold Reversal (ACR, 12.7–14.6 kyr BP), which is a distinct feature only during Termination I among the candidate terminations (i.e., Terminations I, II, IV, and VII) as observed in the EDC δ²H_ice record (Fig. A4). To confirm the age of the Larsen ice, we compared the δ¹⁸O_atm–CH₄ relationship with the ice from core no. 23 at depths of 4.6–10 m with EDC and the WAIS Divide (Fig. 6) and found that the Larsen δ¹⁸O_atm–CH₄ distribution was well matched only at Termination I. The δ¹⁸O_atm of Larsen has a 0.05 ‰ offset with that of the WAIS Divide, which is within the same order of magnitude on the differences observed between core no. 23 and horizontal measurements at a depth of >1.95 m (Fig. A2). Due to the lack of a corresponding δ¹⁸O_atm record for the CH₄ of core no. 23, the difference may also come from the δ¹⁸O_atm records that were estimated by linear interpolation. The δ¹⁸O_atm–CO₂ and CO₂–CH₄ relationships do not clearly match with those of Termination I in the EDC and WAIS Divide records. The CO₂ concentration of the Larsen ice is typically about 10–20 ppm higher compared to EDC and the WAIS Divide (Fig. A3), warranting further investigations. Finally, we confirmed the ages of the Larsen ice with ⁸¹Kr dating, indicating 9–41 and 14–43 ka for ice from the TF and no. 23 cores, respectively (Table 2), and concluded that the Larsen ice covers the Last Glacial Termination (LGT; T1). The ice sample details and krypton dating results are presented in Table 2.

Table 2 Results of ⁸¹Kr and ⁸⁵Kr analysis. Errors in the fifth to seventh column have a 1σ confidence level, whereas upper limits in the fourth column have a 90 % confidence level. The ⁸¹Kr age in the sixth column is given with the statistic error due to atom counting. The systematic error in the ⁸¹Kr age in the seventh column is due to the uncertainty in the half-life of ⁸¹Kr (229 ± 11 kyr). pMKr: percent modern krypton; dpm cc⁻¹: decay per minute per cubic centimeter STP (standard temperature and pressure) of Kr. As ⁸⁵Kr in modern air in Antarctica is around 65 dpm cc⁻¹, the ⁸⁵Kr of samples indicates that alteration by modern air is negligible.

Figure 5 Comparison of Larsen δ¹⁸O_atm with the EDC record. δ¹⁸O_atm records are from Landais et al. (2013) for 0–40 ka and Extier et al. (2018a) for 100–800 ka, respectively. Green vertical bars represent the candidate age intervals for Larsen blue ice. T: termination. The range of Larsen δ¹⁸O_atm values is shown as a vertical grey bar. Horizontal grey bars with dashed lines represent the uncertainty in the measured δ¹⁸O_atm from the Larsen BIA (Fig. A2).

Figure 6 Comparison of the δ¹⁸O_atm–CH₄ relationship in Larsen ice core no. 23 with existing records over glacial terminations. δ¹⁸O_atm and CH₄ data of the WAIS Divide are from Severinghaus (2015) and Rhodes et al. (2017), respectively. δ¹⁸O_atm of EDC is from Landais et al. (2013) and Extier et al. (2018a). The CH₄ record of EDC is from Bazin et al. (2013d). The area within the black box is magnified to compare the Larsen no. 23 record with those from the WAIS Divide T1 and EDC T1.

4.2 Gas ages of the Larsen ice

The CH₄ was severely altered at a depth of 1.95 m, and the only available CH₄ record was from ice core no. 23 below a depth of 4.6 m. Therefore, we developed the tentative gas age by correlating the horizontal δ¹⁸O_atm values to those from the WAIS Divide record on the WD2014 timescale. Then, fine correlation was conducted with the CH₄ record from core no. 23, which corresponds to 13.3–14.1 kyr BP at depths of 4.6–10.4 m (Fig. 7d). We used the high-resolution records of CH₄ from both Larsen core no. 23 and WAIS Divide ice, which greatly improves the precision of age construction compared to the age construction only based on δ¹⁸O_atm.

Spline curves were created for both WAIS Divide and Larsen δ¹⁸O_atm records to reduce artifacts from insufficient sampling resolution and/or ice quality at a 1.95 m depth of the Larsen ice (Fig. 7). The spline curve was drawn after interpolating the original records at 5 m horizontal distance intervals. Three points were selected, which divide the upstream part of the Larsen ice into three equal parts to tie with those in the WAIS Divide record (pink dots in Fig. 7a and b). For the downstream side of the Larsen ice, we selected a local maximum and a minimum as the tie points (purple dots in Fig. 7a and b) because the slope of the δ¹⁸O_atm spline curve is small. When the horizontal distance was >1280 m, the age was constrained by extrapolation using the age–distance relationship at 18–22 ka. Based on the tentative horizontal gas ages and the depth–distance relationship of core no. 23 (obtained from Fig. 4b), we determined the gas age of core no. 23 (∼ 10 m vertical ice core). A small temporal difference (53.6 ± 38.5 years) existed between the CH₄ record of core no. 23 (>4.6 m) and the WAIS Divide record (difference between the red and dark blue lines in Fig. 7d). To eliminate this gap, four tie points (orange dots) were chosen and interpolated (light-blue line). The gas age, corresponding to 13.3–14.1 kyr BP, was obtained more accurately via CH₄ correlation. As a result, the gas age for the Larsen ice at 1.95 m depth was estimated to be 9.2–23.4 kyr BP in the WAIS Divide chronology in 2014 (WD2014). The ages on the AICC2012 (Antarctic Ice Core Chronology) scale were estimated to be 9.4–23.4 kyr BP using the depth and the WD2014 age of EDC (Bazin et al., 2013a; Buizert et al., 2021). The age–distance relationship was reasonable (Fig. A6), showing no abrupt change, as supported by the gradual change in the ice layer dip (Fig. 3c). Refer to Sect. S1​​​​​​​ in the Supplement for gas age uncertainty.

Figure 7 Correlation of Larsen δ¹⁸O_atm and CH₄ records with the WAIS Divide. (a) Interpolated δ¹⁸O_atm value at 5 m from the Larsen BIA with a spline curve. (b) δ¹⁸O_atm record from the WAIS Divide with a spline curve (Severinghaus, 2015​​​​​​​). Six tie points were used to correlate each other. (c) Comparison of synchronized Larsen δ¹⁸O_atm with the EDC (Landais et al., 2013) and WAIS Divide records. (d) Comparison of the CH₄ record from Larsen no. 23 (>4.6 m depth) with the WAIS Divide (Rhodes et al., 2017). Tentative gas age determined by δ¹⁸O_atm correlation with the WAIS Divide is tuned by correlating the CH₄ record using four tie points. The WD2014 timescale for the WAIS Divide and EDC is from Sigl et al. (2016) and Buizert et al. (2021), respectively.

4.3 Ice ages of the Larsen ice

The ice ages of the Larsen ice were determined based on the δ¹⁸O_ice correlation with the TALDICE record (Fig. 9). TALDICE was selected since the coring site is close to the Larsen BIA (∼ 240 km apart) and the direction of the TALDICE site from the Larsen BIA is similar to the ice flow direction of the Larsen ice (Rignot et al., 2011; Mouginot et al., 2012). It is likely that the trend of surface temperature changes in the snow accumulation zone of the Larsen BIA is comparable to that of the TALDICE site. Spline curves were constructed for the horizontal measurement of the surface ice. As described above, the horizontal resolution has been increased by interpolating the original δ¹⁸O_ice records at 5 m intervals. The correlation was established by visually selecting a similar inflection point. Then, to validate the ice age, we estimated the ice age at a 1.95 m depth (see below) and calculated the Δage; Δage for the Larsen ice is defined as (ice age) − (gas age) at a depth of 1.95 m.

To estimate the ice age at a depth of 1.95 m, the age–depth relationship of each ice core sample must be determined. The relationship was deduced by a simple calculation using Eqs. (4) and (5): α, β, and γ are the ice age at each location; D is the horizontal distance between two surface ice samples (at a depth of 20 cm); θ is the average dip of the ice layer; and H+0.2 m is the depth where the age is the same as β (Fig. 8). The average dip was estimated from the GPR profile (yellow box in Fig. 3c). The average dip of the ice layer located near core no. 23 was 4.96^∘, as inferred from Fig. 4b. Then, the ice age at a depth of 1.95 m was calculated using Eq. (6). Refer to Fig. 8 to understand the following equations:

$$\begin{array}{}\text{(4)}& {\displaystyle }H=D\times \tan \mathit{\theta },\text{(5)}& {\displaystyle }\text{Age/depth relation}=(\mathit{\beta }-\mathit{\alpha })/H,\text{(6)}& {\displaystyle }\mathit{\gamma }=\mathit{\left\{}\right(\mathit{\beta }-\mathit{\alpha })/H\mathit{\}}\times (\mathrm{1.95}-\mathrm{0.2})+\mathit{\alpha }.\end{array}$$

Figure 8 Schematic illustration for the mathematical relation between ice ages.

First, eight tie points were selected and linearly interpolated. In this case, Δage showed a minimum value at around 17.5 ka (Fig. A5), which is close to the period of the Last Glacial Maximum (LGM). Since Δage increases when the climate is colder (Schwander et al., 1997), we assumed that either the ice age is undervalued or the gas age is overvalued. Because the δ¹⁸O_atm value of the Larsen ice is comparable with the EDC and WAIS Divide records, the gas age seems to be reasonable and not likely to be overvalued (Fig. 7). Hence, we concluded that the ice age should be revised. In addition, the age–distance relationship of ice shows a rapid age increase near ∼ 1000 m (Fig. A6), which is not supported by the average dip of the ice layer deduced from the GPR profile (Fig. 3). As the dip decreases to the downstream ice, no abrupt increase in the slope of the age–distance relationship is expected. However, this is the case only when we assume no abrupt change in the snow accumulation rate.

After excluding one tie point (indicated by the red dotted line in Fig. 9), the Δage of the Larsen ice shows a maximum value at approximately 17.5 kyr BP, which appears more appropriate (Fig. A5). For records with horizontal distances less than 107 and greater than 1310 m, extrapolation was conducted to estimate the ice age using the age–distance relationship at 7–11 and 18–23.8 kyr BP, respectively. In conclusion, the ice age was estimated to be 5.6–24.7 kyr BP (WD2014) for the surface ice. However, as noted above, local age inversion may have occurred in the upstream ice areas. Hence, the age constraint for upstream ice should be considered cautiously. The ice age calculated for a 1.95 m depth also seems reliable because the water isotope record is comparable with the surface ice result when plotted with the ice age (Fig. 9c). Note that ice age synchronization via water isotopes is challenging in climatically stable periods such as the LGM; measurements of dust and ice chemistry could be applied in the future to refine the ice chronology presented here. Refer to Sect. S1 in the Supplement for ice age uncertainty and Δage uncertainty.

The established horizontal ice chronology shows that temporal resolutions of 10 and 17 yr m⁻¹ are available for the Holocene (5.6–12 kyr BP) and for the last deglaciation (12–24.6 kyr BP), respectively, at the surface ice from the Larsen BIA. These estimated resolutions are higher than those of the deep ice core TALDICE covering this time interval (∼ 18 yr m⁻¹) (Masson-Delmotte et al., 2011). In addition, the Larsen BIA allows for collecting large quantities of ice from the surface. Assuming that the average age–depth relationship of ice core no. 23 (∼ 189 yr m⁻¹) is maintained throughout the ice, the ice age may then reach ∼ 60 ka near the bedrock (∼ 240 m of ice), and thus, ice may be obtained for at least ∼ 60 ka via shallow coring at the Larsen BIA.

Figure 9 Correlation of Larsen δ¹⁸O_ice with the TALDICE record. (a) Interpolated δ¹⁸O_ice value at 5 m from the Larsen BIA with a spline curve. (b) δ¹⁸O_ice record from TALDICE with a spline curve (Stenni et al., 2011; Bazin et al., 2013c). The ice age of the surface ice from Larsen was estimated by correlating seven inflection points of the spline curve (indicated by black dotted lines) with TALDICE. (c) Comparison of near-surface δ¹⁸O_ice with those in 1.95 m depth at the Larsen BIA.

4.4 Estimation of past surface temperature and accumulation rate

As discussed in Sect. 3.2, δ¹⁸O_ice and δ²H_ice might have been enriched during the sublimation of ice. Therefore, it is inadequate to use stable water isotopes to estimate past surface temperatures. Instead, the analytical framework developed by Buizert (2021) allows us to estimate the past surface temperature using the δ¹⁵N–N₂ and Δage. The methods for estimating the surface temperature and accumulation rate are described in detail in the Appendix (Appendix D). The results are presented in Table S11 and Fig. 10. The accumulation rates for TALDICE, EDC, Taylor Dome, and Taylor Glacier are also shown for comparison (Veres et al., 2013; Baggenstos et al., 2018). The youngest three reconstructed values of the Larsen BIA were rejected (red line in Fig. 10) due to the large relative uncertainty in Δage (including the possibility of negative Δage) and because the implied 10.5 ^∘C increase is for just ∼ 1200 years. This might be due to a low bias in the reconstructed ice age, which leads to a lower Δage. Hence, the ice age in this part should be more strongly constrained by using dust concentrations in a future study.

During the last deglaciation, the surface temperature increased by approximately 15 ± 5^∘C (1σ) (from 24.3 to 10.6 kyr BP), which is greater than those for any other ice core sites including nearby Talos Dome, where a recent reconstruction suggests 7.1 ± 2 ^∘C of deglacial warming (Buizert et al., 2021). However, a cautious interpretation is recommended because of the large uncertainty. The large magnitude of the reconstructed temperature change is a consequence of Δage, which is around 300 % larger during the Last Glacial Maximum (LGM) than during the Holocene in our reconstruction. By comparison, most Antarctic sites show only a 60 % to 120 % Δage increase (Buizert et al., 2021). While the Larsen BIA is currently close to the open ocean, during the LGM, the Antarctic ice sheet reached the continental shelf break covering the entire Ross Sea Embayment (Conway et al., 1999). Climate models identified strong cooling over the Ross Sea sector, reflecting the increased elevation and increased albedo of the extended ice sheet (Buizert et al., 2021). We suggest that this enhanced cooling may have affected the original deposition site (upstream accumulation area) of the Larsen BIA. We note that the past ice sheet thickness at Larsen and the upstream distance traveled, since deposition is poorly constrained; both of these impact the magnitude of the reconstructed temperature change. The temporal isotope slope (α_T=0.58 ‰ K⁻¹) is lower than that for any other ice core site (Buizert et al., 2021) and smaller than the spatial regression slope of around 0.8 ‰ K⁻¹ (Masson-Delmotte et al., 2008​​​​​​​). However, the temporal isotope slope we reconstructed is likely to be a lower bound because the sublimation might have enriched the δ¹⁸O_ice of the glacial (downstream) part of the Larsen BIA. The reconstructed surface temperature change of 15 ± 5 ^∘C assumes that the selected tie points are correct. However, if the δ¹⁸O_ice features used for the matching are not climatic in origin but rather reflect local effects (such as sublimation intensity and accumulation controls by surface slope), then our tie points would be incorrect. Future measurements of ice chemistry and dust loading may improve our Δage estimates, which will allow for a refined estimate of past glacial cooling.

From 24.3 to 10.6 kyr BP, the accumulation rate increased by a factor of 1.7–4.6 (from 0.033 ± 0.007 to 0.103 ± 0.042 m ice per year). During the same time period, the accumulation rates at the Taylor Dome and TALDICE sites increased by a factor of ∼ 15.4 and ∼ 2.4, respectively (Veres et al., 2013; Baggenstos et al., 2018). Based on the three records (TALDICE, Larsen blue ice, and Taylor Dome), the accumulation rate increased more in southern Victoria Land than northern Victoria Land. The accumulation rate at TALDICE and EDC began to decrease transiently around 14.5 kyr BP following the Antarctic Cold Reversal (ACR), while the reconstructed accumulation rates at the deposition site of the Larsen BIA and Taylor Dome keep increasing across the ACR. We acknowledge that the accumulation rate of the deposition site of the Larsen ice younger than 14 kyr BP remains poorly constrained because of the large uncertainty but is well constrained for the older part (>14 kyr BP). The accumulation rate at the deposition site of the Larsen ice is lower than that of TALDICE during 14–21 kyr BP and exceeds the accumulation rate of TALDICE after 14 kyr BP. The Ross Ice Shelf (RIS) retreated during the last deglaciation (Ship et al., 1999; Yokoyama et al., 2016), and as the RIS retreats, the storm track migrates to southern Victoria Land from the northern part and increases precipitation in southern Victoria Land (Morse et al., 1998; Aarons et al., 2016; Yan et al., 2021). Therefore, we speculate that the storm track affects the original deposition site of the Larsen BIA more than the Talos Dome after 14 kyr BP. The accumulation rate at Taylor Dome (TD) during 14–21 kyr BP shows a lower value than that from TALDICE and the Larsen ice (Fig. 10). However, the accumulation rate increased by a factor of ∼ 5.2 at Taylor Dome, while those at Larsen and TALDICE increased by a factor of ∼ 2.1 and ∼ 1.8, respectively. Overall, it appears that as the RIS retreats, the accumulation rate at southern Victoria Land increases more than northern Victoria Land. This interpretation is in line with the previous studies (Morse et al., 1998; Aarons et al., 2016; Yan et al., 2021) and may help studies for reconstructing past atmospheric circulation associated with the retreat of the RIS. Likewise, a strong accumulation increase across the last deglaciation was seen at the coastal Law Dome site in the Indian Ocean sector (Van Ommen et al., 2004), also attributed to increases in storm-derived precipitation. However, a spatial difference in the original deposition site of the upstream and downstream Larsen ice is not constrained; this must be known to allow for a better interpretation of the accumulation rate.

Figure 10 Reconstructed surface temperature and accumulation rate over the Last Glacial Termination (Termination I). Accumulation rates of TALDICE and EDC are from Veres et al. (2013). Accumulation rates of Taylor Dome and Taylor Glacier are from Baggenstos et al. (2018). The three youngest values of Larsen blue ice are rejected (red line) due to the large uncertainty in the Δage that translates into the possibility of negative accumulation rates and an unexpected 10.5 ^∘C increase in ∼ 1200 years. ΔT is a value relative to 10.6 kyr BP; 10.6 kyr BP is the youngest time period available. The uncertainty for ΔT and the accumulation rate (1σ) was estimated through error propagation.

5 Conclusions

Based on the shapes of dust bands, GPR profiles, and analytical data (δ¹⁸O_atm, δ²H_ice, and CH₄) we conclude that a large ice section of the ice at the Larsen BIA (the downstream ice, from ice cores no. 23 to no. 200) monotonically increase along the ice flow direction. The vertical profile of δ²H_ice in core no. 23 correlated well with the record from the horizontal surface ice samples, showing that the stratigraphy of ice is not significantly disturbed and does not cause an age inversion at the site. The negative d-excess values indicate sublimation of the ice and the possibility of enrichment in the δ²H_ice and δ¹⁸O_ice values during sublimation. CH₄ is well preserved at depth in ice cores (>4.6 m for ice core no. 23), and gas isotopes (δ¹⁸O_atm and δ¹⁵N–N₂) are well preserved at 1.95 m depth.

Horizontal δ¹⁸O_atm values at a depth of 1.95 m in the Larsen BIA reveal a typical glacial termination (the Larsen δ¹⁸O_atm shows both negative and >1.0 values). Along with glaciological ice flow principles, stable water isotopes and the correlation of CH₄ concentration and δ¹⁸O_atm with existing ice core records, we suggest that the gas ages of the studied Larsen ice cover 9.2–23.4 kyr BP and that ice ages cover 5.6–24.7 kyr BP. This is additionally supported by ⁸¹Kr ages from local ice cores (no. 23 and TF). We provide high-precision ages for ice and fossil air trapped in blue ice at the Larsen BIA. Because well-constrained ages and high-quality proxies are essential for paleoclimate studies, our work may help future studies on BIAs across the Antarctic continent. For example, using the ice samples from the Larsen BIA, we suggest a tentative climate reconstruction of surface temperature (15 ± 5 ^∘C increase) and accumulation rate (increased by a factor of 1.7–4.6) at the deposition site of the Larsen ice during the last deglaciation.

Combining information derived from ice samples from Taylor Glacier, Allan Hills, and the Larsen BIA may also provide insights regarding the growth and/or retreat of the Ross Ice Shelf. In addition, various sites may be beneficial for studying the paleoclimate, since the ice quality for paleoclimate studies depends on local glaciological conditions. In particular, since a large amount of ice with the same age is exposed at the surface in the BIA, this study may facilitate further research on paleoclimate for which large quantities of ice are required (such as isotope analysis of trace gases or elements in the ice).

Appendix A: CH₄, CO₂, δ¹⁸O_atm, and δ¹⁵N–N₂ measurements of the Larsen ice samples

Figure A1 Vertical profiles of greenhouse gas concentrations at the Larsen BIA ice cores. (a) CH₄ records. (b) CO₂ records. CH₄ was measured using a wet-extraction method, while CO₂ was measured via both dry- and wet-extraction methods at SNU and NIPR, respectively. Results of the TF core are from Jang et al. (2017).

Figure A2 Nitrogen and oxygen gas isotope record of ice core no. 23 and horizontal measurement. Horizontal measurement is conducted using ice cores at 1.95 m depth. Depth of ice core no. 23 was converted to horizontal distance at 1.95 m depth by using Eq. (3). The uncertainty in a measured δ¹⁸O_atm at a depth of 1.95 m is assumed to be ± 0.05 ‰ (see Sect. 3.3.2).

Appendix B: Gas age estimation for the Larsen BIA

Figure B1 ​​​​​​​Comparison of the CO₂–CH₄ and δ¹⁸O_atm–CO₂ relationship in Larsen ice core no. 23 with existing records during glacial terminations. WAIS Divide: Marcott et al. (2014), Severinghaus (2015)​​​​​​​, and Rhodes et al. (2017). EDC: Bazin et al. (2013d), Monnin et al. (2001, 2004), Siegenthaler et al. (2005), Lourantou et al. (2010), Schmitt et al. (2012), Landais et al. (2013), Bereiter et al. (2015), and Extier et al. (2018a). Vostok: Petit et al. (1999) and Bender (2002).

Appendix C: Ice age estimation for the Larsen BIA

Figure C1 δ²H_ice records from EDC. δ²H_ice records for T1, T2, T4, and T7 on the AICC2012 scale are from Bazin et al. (2013b). The blue vertical bar represents the time interval of the ACR (12.7–14.6 kyr BP).

Figure C2 EDC, TALDICE, and Larsen BIA Δage. Δage for the Larsen BIA is defined as (ice age) − (gas age) at 1.95 m depth. The EDC WD2014 chronology is from Buizert et al. (2021).

Figure C3 Relation between horizontal distance and age.

Appendix D: Surface temperature and accumulation rate estimation

The surface temperature change can be estimated using the following equation (Buizert, 2021), where ${\stackrel{\mathrm{‾}}{T}}^{\mathrm{2}}=T_{\mathrm{1}}{T}_{\mathrm{2}}$, R is the ideal gas constant (8.314 J mol⁻¹ K⁻¹), L is the lock-in depth, γ=1, and E_T=23.5 kJ mol⁻¹.​​​​​​​

$$\begin{array}{}\text{(D1)}& T_{\mathrm{2}}-T_{\mathrm{1}}={\displaystyle \frac{R{\stackrel{\mathrm{‾}}{T}}^{\mathrm{2}}}{(\mathrm{1}+\mathit{\gamma })E_{\mathrm{T}}}}\ln \left[{\displaystyle \frac{\mathrm{\Delta }{\mathrm{age}}_{\mathrm{1}}}{\mathrm{\Delta }{\mathrm{age}}_{\mathrm{2}}}}\cdot {\left({\displaystyle \frac{L_{\mathrm{1}}}{L_{\mathrm{2}}}}\right)}^{\mathit{\gamma }}\right]\end{array}$$

To find the lock-in depth L, we first estimated the diffusive column height via the barometric equation (Eq. D2) (Sowers et al., 1992) then added the convective zone of the firn column; we assume the length of the convective zone to be 4 m. In Eq. (D2), Δm is the difference in the molecular weight between ²⁸N₂ and ²⁹N₂ (0.001 kg mol⁻¹), g is the gravitational acceleration constant (9.82 m s⁻²), H_diff is the diffusive column height, R is the ideal gas constant (8.314 J mol⁻¹ K⁻¹), and T is the temperature (K).

$$\begin{array}{}\text{(D2)}& {\mathit{\delta }}^{\mathrm{15}}\mathrm{N}=\left(e^{\frac{\mathrm{\Delta }m\times g\times H_{\mathrm{diff}}}{R\times T}}-\mathrm{1}\right)\times \mathrm{1000}\end{array}$$

We first applied the modern mean annual temperature of the Larsen BIA (−24.4 ^∘C) to calculate the lock-in depth at all times. Then we applied Eq. (D1) to calculate T₂ and then applied Eq. (D2) again to find the correct lock-in depth. Finally, we applied Eq. (D1) once again to estimate the correct surface temperature.

The ice-equivalent accumulation rate (A) can be estimated by using the following relationship between Δage and ice-equivalent lock-in depth (L_IE) (Buizert, 2021). The ratio of $L_{\mathrm{IE}}/L$ is known to be ∼ 0.70 (Parrenin et al., 2012) and was used to find L_IE.

$$\begin{array}{}\text{(D3)}& \mathrm{\Delta }\mathrm{age}={\displaystyle \frac{L_{\mathrm{IE}}}{A}}\end{array}$$

All of the uncertainties in the input parameters were assumed to be independent. Then, final uncertainty for the surface temperature and accumulation rate was estimated through error propagation.