Age–depth distribution in western Dronning Maud Land, East Antarctica, and Antarctic-wide comparisons of internal reflection horizons

Abstract

Radio-echo sounding provides the opportunity to study the internal architecture of ice sheets through imaging stratified englacial reflections, known as internal reflection horizons (IRHs). They represent consistent time horizons formed at the former ice-sheet surface and buried over time, thus reflecting the ice sheet's age–depth architecture. Their analysis allows crucial insights into past and present glaciological conditions, e.g. bed topography, surface and basal mass balance, and physical properties and ice dynamics. This study presents a comprehensive data set of IRHs and insight into the age–depth distribution in western Dronning Maud Land (DML), East Antarctica, spanning the Holocene to the Last Glacial Period (4.8–91.0 ka). Using data from various radar systems deployed by the Alfred Wegener Institute between 1996 and 2023, we traced and dated nine IRHs over an area of 450 000 km². A precise age could be assigned to the IRHs by two-way travel time to depth conversion and employing radar forward modelling based on conductivity peaks of the EPICA DML ice core. Six IRHs correlate with the timing of past volcanic eruptions, and our findings suggest that most IRHs correspond to IRHs of similar age in other regions of East and West Antarctica, thus likely originating from the same physical reflectors at depth, although some could not be physically connected. This work enhances understanding of the englacial architecture and relationships with snow accumulation and ice-dynamic processes of this sector of the Antarctic ice sheet and provides boundary conditions for numerical ice flow models and paleoclimatic studies.

1 Introduction

Studying the dynamics of the Antarctic ice sheet (AIS) through a combination of geophysical observations and ice-sheet modelling is crucial for understanding its response to climate change and better projecting future sea-level rise . Deducing its past evolution helps us to understand the processes driving ice flow, accumulation, and basal melting. As one unique observable property, the ice sheet's internal stratigraphy represents a valuable record of past kinematics and dynamics , offering insights into atmosphere–ice–ocean interactions and how the ice-sheet system responded to natural climatic variations over a broad range of timescales ranging from hundreds to hundreds of thousands of years e.g..

Figure 1

Radar profiles used for IRH tracing. (a) Radar profiles are sorted by system: EMR short pulse (white lines), UWB narrowband (orange lines), and UWB wideband (blue lines). (b) Radar data sorted by acquisition season. Although some profiles extend down to the ice shelf, the IRH tracing region is restricted upstream of the Maud Belt. The background map shows ice surface flow velocity superimposed on the hillshaded surface elevation . The IMBIE (Ice sheet Mass Balance Inter-comparison Exercise) drainage basins are marked by thin black lines, and deep ice core sites are marked with a white circle.

[Figure omitted. See PDF]

Radar data provide the opportunity to study the internal architecture of ice sheets (i.e. their stratigraphy) via continuous stratified englacial reflections . These internal reflection horizons (IRHs) occur when radar waves are reflected at boundary layers where the dielectric properties of the ice change . Dielectric contrasts in ice sheets are primarily caused by variations in density , electrical conductivity , and changes in the preferred orientation of ice crystals . IRHs that are located sufficiently deep below the ice surface, typically within a few hundred metres in Antarctica, are primarily linked to conductivity and density contrasts and can sometimes be detected thousands of kilometres across the ice sheet . Strong conductivity contrasts are often linked to deposits of acidic materials from major volcanic eruptions . Moreover, IRHs are often considered isochronous because they form simultaneously at the ice-sheet surface and are then buried under subsequent snow accumulation, carrying with them changes in dielectric properties which radar systems are sensitive to .

Information about the internal stratigraphy and age–depth architecture of the Antarctic ice sheet is crucial as it serves as a climatic and ice-dynamic record. For instance, it provides insights into past accumulation rates and basal melting, informing us about past surface and basal mass balance . Since many IRHs can be dated with ice cores, they can also serve as calibration points for ice-sheet models and to model basal conditions . Additionally, the geometry, depth, and continuity of IRHs contain information on the cumulative deformation due to ice flow, making them suitable passive markers for ice-dynamic activity and its changes over time, e.g. deciphering folding processes and ice stream activity , as well as ice-dynamic processes causing IRH discontinuity .

The generation and compilation of IRHs using radar data in Antarctica and Greenland remain largely a manual task . While automated methods can successfully trace IRHs to a certain extent (see e.g. review by ), these algorithms often fail to consistently follow the same reflection across multiple intersecting radar profiles and different radar systems, where reflections have a different radar signature due to their different wavelengths. Nevertheless, IRHs have been traced and dated across large parts of East Antarctica , West Antarctica , and Greenland . For Antarctica, in particular, these studies have been motivated by the AntArchitecture initiative, a Scientific Committee on Antarctic Research Action Group which aims to build a 3D age–depth model of Antarctica from IRHs .

Table 1

Alfred Wegener Institute (AWI) radar system specifications. The range resolution applies for a dielectric permittivity of ${\mathit{\epsilon }}_{\mathrm{r}}=\mathrm{3.15}$ in ice.

Radar system	Developer	Frequency	Transmit signal	Range resolution	Seasons
EMR	TU HHa	150 MHz	60 ns burst	$\sim \mathrm{5}$ m	1996/97, 1997/98, 1998/99
					2000/01, 2001/02, 2002/03
					2007/08, 2010/11, 2016/17
UWB	CReSISb	180–210 MHz	1, 3, and 10 $\mathrm{\mu }\mathrm{s}$ chirp	4.3 m	2018/19
		(narrowband)
UWB	CReSISb	150–520 MHz	1, 3, and 10 $\mathrm{\mu }\mathrm{s}$ chirp	0.35 m	2023/24
		(wideband)

^a Hamburg University of Technology. ^b Center for Remote Sensing and Integrated Systems, formerly known as Center for Remote Sensing of Ice Sheets.

In this study, we present a comprehensive and detailed insight into the dated radiostratigraphy of western Dronning Maud Land (DML; East Antarctica), spanning from the Holocene to the Last Glacial Period, using radar data collected over the past 3 decades by various recording systems of the Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research (AWI; Fig. ). Through radar forward modelling based on conductivity data from the EPICA DML (EDML) ice core, we can accurately date IRHs and, in most cases, link them to deposits from past volcanic eruptions. We discuss the age–depth distribution and links to glaciological parameters, such as ice flow velocity, accumulation rates, and bed topography. Moreover, we infer the extent to which our dated IRHs can be correlated with those from other Antarctic-wide studies and demonstrate the potential for extrapolating some of these IRHs to much larger areas where an AWI radar system with coarser vertical resolution was used.

2 Data and methods

2.1 Radar data

In this study, we utilize three types of radar products that provide sufficiently high range resolution (0.35–5 m) and penetration depth (Table ). These radar products were obtained during campaigns conducted by AWI between 1996 and 2023 with two different radar systems (Fig. ). The radar data were collected using AWI's polar aircraft Polar 2, Polar 5, and Polar 6 , operating between the stations Neumayer and Kohnen , Novolazarevskaya, Troll, and Dome Fuji.

2.1.1 AWI EMR short-pulse data

The largest pool of radar data is provided by AWI's EMR (Electromagnetic Reflection System), which has been operational in Antarctica and Greenland since 1994 . For IRH tracing, we utilize EMR radar data from nine Antarctic seasons spanning from 1996/97 to 2016/17 (Fig.  and Table ). The EMR system consists of two dipole antennas mounted underneath the wings of AWI's polar aircraft. The transmission signal is a burst with pulse lengths of 60 or 600 ns at a frequency of 150 MHz, with a pulse repetition frequency (PRF) of 20 kHz. For this study, we exclusively utilize EMR 60 ns pulse data (short-pulse data) due to their higher range resolution to detect IRHs compared to the 600 ns pulse data see. The short-pulse data have a range resolution of approximately 5 m and a vertical sample interval of 13.33 ns. EMR radar data processing comprises filtering and along-track stacking with a factor of 7 or 10, which results in a mean trace spacing of $\sim \mathrm{35}$ and $\sim \mathrm{50}$ m, respectively .

EMR survey lines cover the entire western DML (Fig. ). Surveys between 1997 and 1999 focused on the pre-site survey of the EDML ice core site , while surveys between 2001 and 2011 were used to investigate the ice internal structure as well as ice thickness for solid earth geophysical studies in DML . The survey in 2016/17 extensively surveyed ice thickness at the Dome Fuji ice core site .

2.1.2 AWI ultra-wideband data

In addition, we use two data sets acquired with AWI's ultra-wideband (AWI UWB) radar system. The system is an improved version of the Multichannel Coherent Radar Depth Sounder (MCoRDS, version 5), which was developed at the Center for Remote Sensing and Integrated Systems (CReSIS) at the University of Kansas and has been in operation since 2016 . The radar configuration for both campaigns in 2018/19 and 2023/24 consisted of an eight-element antenna array mounted under AWI's Polar 5 or Polar 6 Basler BT-67 aircraft's fuselage, which serves as a transmit and receive antenna array. The transmission signal is composed of several staged modulated chirp signals, which provide high resolution at different depths . Standard processing techniques were performed with the OPR Toolbox Open Polar Radar Toolbox, formerly termed CReSIS Toolbox;. The main steps include motion compensation, pulse compression, synthetic aperture radar (SAR) focusing in the fk (frequency–wavenumber) domain, and array processing. Aeroplane location and orientation are from Global Positioning System (GPS) precise point positioning (PPP) post-processed with a final estimated accuracy (commercial software package Waypoint 8.4) of better than 3 cm for latitude and longitude and better than 10 cm for altitude.

The two AWI UWB campaigns were flown in two different acquisition settings in the 2018/19 and 2023/24 seasons, respectively. In the 2018/19 season, data were acquired with a frequency range between 180–210 MHz (narrowband), which corresponds to a range resolution of 4.3 m. The acquisition settings were optimized to sound deep englacial reflections as well as the bed topography at the onset of Jutulstraumen Glacier . AWI UWB data in the 2023/24 season were acquired with a frequency range between 150–520 MHz (wideband), with a range resolution of $\sim$ 35 cm. The acquisition settings were chosen to map near-surface IRHs at high spatial resolution in the upstream part of the Jutulstraumen catchment as well as the central plateau south and west of the EDML ice core.

2.2 Internal reflection horizon tracing

For tracing IRHs, we utilize the commercial software Paradigm 19 from AspenTech Subsurface Science & Engineering. To enhance the visibility of isochrones at greater depths, we applied an automatic gain control (AGC) algorithm to the radargrams. IRHs were traced with the following rationale: (1) existing IRHs, such as the 38 and 74 ka isochrones previously traced by on segments of the same data we used here; (2) the shallowest and deepest clearly identifiable reflection in the EMR short-pulse data was traced; and (3) additional isochrones were traced at regular intervals between those defined in (1) and (2), which were clearly discernible in all radar products.

All isochrones were manually traced with the assistance of a semi-automatic picker, which searches for the local maximum in return power within a window width of 5 ns along two picks. Based on the intersection points of the radar profiles, IRHs were transferred between the different radar products. In cases where the assignment at the intersection point was not clear or the radar reflection was no longer traceable within a profile, no picks were created. We clarify that we use the existing IRHs from in DML and integrate these data for our representation, dating, and analyses.

2.3 Ice base determination

The determination of the ice base reflection for this study is composed of existing and newly traced data. We used existing ice base pick data from the AWI UWB 2018/19 campaign from . Ice base determination from the AWI UWB 2023/2024 season was generated using the OPR Toolbox and Paradigm. For the EMR data, we were able to partially rely on existing data , which were either directly traced in the 60 ns short-pulse data or, if flown in toggle mode, traced in the 600 ns long-pulse data and projected to the 60 ns data. For EMR short-pulse profiles used in this study that did not have existing ice base picks, the ice base reflection was traced at locations where IRHs were present and the base reflection was clearly visible.

Figure 2

Comparison between simulated radar data for 150 and 335 MHz (150 MHz for the EMR profile and 335 MHz for the AWI UWB wideband profile) based on measured conductivity at the EDML ice core and radar data collected near EDML. (a) Input parameters for dielectric conductivity and relative permittivity. (b) Comparison between the gprMax simulated radar trace at 150 MHz with overlaid conductivity peaks (yellow) and the EMR profile 20023154 (transmit frequency is 150 MHz). (c) Comparison between the gprMax simulated radar trace at 335 MHz with overlaid conductivity peaks and the UWB wideband profile 20231211_01_024 (centre frequency is 335 MHz). The correspondence between simulated reflections, conductivity peaks, and the reflections of our IRHs is marked with coloured arrows. Note that the simulated radargram is the 2-D representation of a single trace (without AGC), and the measured radargram has an along-track range of 20 km (with AGC) to enable a better comparison.

[Figure omitted. See PDF]

2.4 IRH depth and normalized depth

We derived the depth of all IRHs by calculating the two-way travel time (TWT), $t$, between the ice surface and the IRH. Then, we converted the TWT to depth e.g. using

1$$Z_{\mathrm{irh}}={\displaystyle \frac{t{c}_{\mathrm{air}}}{\mathrm{2}\sqrt{{\mathit{\epsilon }}_{\mathrm{r}}^{\prime }}}}+Z_{\mathrm{firn}},$$ where $c_{\mathrm{air}}$ is the electromagnetic wave speed in air ($c_{\mathrm{air}}=\mathrm{2.9971}\times {\mathrm{10}}^{\mathrm{8}}$ m s⁻¹). For the calculation of IRH depths, we use a constant dielectric permittivity of ${\mathit{\epsilon }}_{\mathrm{r}}^{\prime }=\mathrm{3.15}$. This value deviates slightly from the empirically determined value at EDML , but it is a realistic estimate for the spatial extrapolation within our study area and has been widely used in other studies. Additionally, we apply a firn correction of 13 m , which is the best-known value in this region. However, even though a 13 m firn correction was applied uniformly, snowfall and accumulation rates vary across western Dronning Maud Land e.g., potentially affecting depth calculations of IRHs slightly, i.e. on the order of a few metres.

Additionally, we calculated the normalized depth of the IRHs within the ice column, which provides a better comparison to determine whether IRHs lie deeper or shallower within the ice column in specific areas. For this, we divided the IRH depth below the ice surface by the respective ice thickness. Here, we only use ice thickness values at locations where ice base picks are available, as differences in ice thickness compared to gridded products like BedMachine Antarctica can result in significant discrepancies, which can strongly influence the relative position of IRHs within the ice column. These discrepancies are due to interpolation methods and the fact that not all ice thickness data used here are yet included in BedMachine Antarctica.

2.5 Internal reflection horizon dating

To determine the age of the IRHs, we used two different methods: (1) age dating by converting the IRH radar wave two-way travel time (TWT) from the ice surface to depth and correlating it with the age–depth scale of the EDML ice core e.g. and (2) matching reflections of our IRHs with simulated radar reflections based on conductivity peaks derived from measured DEP (dielectric profiling) data e.g.. Method (1) has the advantage of allowing the determination of an age for each IRH, albeit with the acceptance of a larger error in depth and age assignment. Method (2) offers greater accuracy and linkage to volcanic events because the conductivity peaks represent the physical origin of reflections ; however, this approach is not applicable for all IRHs due to some reflections not rising above background noise in the radar data.

Table 2

Overview of IRH key characteristics, such as the mean IRH TWT at EDML, mean IRH depth at EDML derived from TWT-to-depth conversion, IRH depths at EDML derived from matching conductivity peaks, IRH ages based on the TWT-to-depth converted depths, IRH ages based on the depth from matched conductivity peaks, cumulative length of IRHs (line kilometres), and link to possible volcanic eruption ages. The IRH name contains the region (DML) and age in ka (e.g. 4_8_ka for 4.8 ka). All TWT-to-depth converted depth values are relative to the ice surface reflection at first break. The conductivity peak depths are from the DEP data set and refer to the depth of the maximum value of the respective conductivity peak. The IRH ages highlighted in bold italics in the TWT-to-depth-based age and conductivity-peak-based age columns are the ones we have selected as the final IRH ages.

IRH name	Mean TWT	Mean depth	Depth of con-	TWT-to-depth-	Conductivity-	Volcanic	IRH
	at EDML	at EDML	ductivity peak	based age	peak-based age	eruption age	line km
DML_4_8_ka	3789 ns	347 m		4.8 $\mathbf{\pm }$ 0.62 ka			28 597 km
DML_7_2_ka	5503 ns	477 m	473 m	7.3 $\pm$ 0.76 ka	7.2 $\mathbf{\pm }$ 0.04 ka		29 872 km
DML_11_1_ka	7704 ns	662 m	650 m	11.2 $\pm$ 0.93 ka	11.1 $\mathbf{\pm }$ 0.03 ka	11.3 kaa	25 003 km
DML_15_4_ka	9527 ns	816 m		15.4 $\mathbf{\pm }$ 1.19 ka		15.6 kab	13 926 km
DML_25_6_ka	12 677 ns	1081 m	1070 m	25.9 $\pm$ 2.08 ka	25.6 $\mathbf{\pm }$ 0.03 ka	25.3 kac	24 212 km
DML_38_1_ka	15 661 ns	1332 m	1311 m	39.0 $\pm$ 2.69 ka	38.1 $\mathbf{\pm }$ 0.02 ka	38.1 kad	28 231 km
DML_52_2_ka	18 470 ns	1569 m	1552 m	53.5 $\pm$ 4.22 ka	52.2 $\mathbf{\pm }$ 0.07 ka	52.3 kae	20 169 km
DML_73_4_ka	22 217 ns	1884 m	1868 m	74.7 $\pm$ 5.34 ka	73.4 $\mathbf{\pm }$ 0.88 ka	73.9f/75.0 kag	25 361 km
DML_91_0_ka	24 723 ns	2095 m	2080 m	92.4 $\pm$ 7.52 ka	91.0 $\mathbf{\pm }$ 1.10 ka		8151 km

^a Large Northern Hemisphere higher-latitude eruption rank 18 in.^b Large lower-latitude or Southern Hemisphere eruption rank 32 in.^c Taupo Oruanui eruption .^d Large lower-latitude eruption rank 3 in.^e Large bipolar eruption rank 5 in.^f Toba eruption dated to 75.0 $\pm$ 0.9 ka see references in.^g Toba eruption dated to 73.88 $\pm$ 0.32 ka see references in.

2.5.1 IRH dating by TWT-to-depth conversion

For the IRH dating using the AICC2023 age–depth scale of the EDML ice core , we used three radar profiles from the seasons 1997/98, 2001/02, and 2023/24, which are located close to the drilling site. Based on the radar traces closest to the ice core for each radar profile (25 m for profile 19983101, 65 m for profile 20023154, and 280 m for profile 20231211_01_024), we determined the travel time between the ice surface reflection and the IRHs. For the depth conversion at EDML, we used the constant radar wave speed in ice suggested by of $\sim$ 1.69 $\times$ 10⁸ m s⁻¹, corresponding to a dielectric permittivity of ${\mathit{\epsilon }}_{\mathrm{r}}^{\prime }=\mathrm{3.145}$. Additionally, we added a value of 13 m to the IRH depths as a firn correction, calculated from measurements of the complex permittivity of shallow firn cores in the region . To account for the variation in acquisition years and radar products, we averaged the depth values for all three radar profiles. We disregarded the snow accumulation between the recording times of the radar profiles, as the average accumulation rate at EDML is approximately 5 cm ice equivalent per year, which amounts to about 1.25 m over 25 years and is significantly below the range resolution of 5 m of the EMR system.

For the depth determination error, we included the following parameters: (1) the standard deviation of the depth differences of the respective IRHs from the three radar profiles; (2) the smallest range resolution of 5 m from the EMR system; (3) the span of the tracing window of 5 ns, which corresponds to approximately 0.5 m; and (4) an error in the dielectric permittivity of 1 % . In addition to the age error resulting from the depth error range, the uncertainty of the age in the AICC2023 chronology is also included.

2.5.2 IRH dating by DEP-based radar forward modelling

For age assignment based on simulated radar reflections derived from the measured conductivity of the EDML ice core , we used the open-source software gprMax to simulate electromagnetic waves in three-dimensional space. This software solves Maxwell's equations using the finite-difference time-domain (FDTD) method and Yee cells , which has proven effective for 1-D and 2-D simulations in ice sheets and glaciers . gprMax uses the same algorithm as already successfully employed for ice-core-based forward modelling by but includes an efficient solution for model boundary conditions using perfectly matched layers PMLs; while still offering reasonable computation times for 2-D models. All simulations presented here were performed with gprMax version 3.1.5 (Big Smoke).

Figure 3

Representation of a subset of radargrams (AGC version) of the three different radar products used in this study and traced IRHs in the western DML. (a) Visualization of selected radar data in a 3D canvas with the bed topography and traced IRHs. Radar profiles with a black line at the top of the fence radargrams are EMR short-pulse profiles, those with a red line are narrow-bandwidth profiles, and those with a yellow line are wide-bandwidth profiles. (b) Representation of an EMR short-pulse radargram (profile ID 20023154), (c) a narrow-bandwidth radargram (profile ID 20181227_03_001), and (d) a wide-bandwidth radargram (profile ID 20231211_01_024) with IRHs. The radargrams in panels (b)–(d) are located close to the EDML ice core and are aligned with the surface reflection.

[Figure omitted. See PDF]

For our simulations, we defined a 2-D model domain extending 2.4 m in the $x$ direction (width) and 2774 m in the $z$ direction (depth) to represent the full depth of the ice core. To reduce the model size, the $y$-direction extension was limited to a single cell size, rendering the model in 2-D to effectively operate in the transverse-magnetic mode. The cell size in all dimensions is 0.02 m, which enables high resolution and numerical stability and avoids numerically induced dispersion by ensuring the cell size is at least 10 times smaller than the smallest wavelength in our model. Consequently, the time step $\mathrm{\Delta }t$ is 0.047 ns for our cell size. The transmit and receive antennas are placed 1 m below the upper limit of the model domain, in the centre of the $x$ domain (zero offset). The simulation window for all simulations is 30 $\mathrm{\mu }\mathrm{s}$, covering the entire ice thickness at EDML. The antenna represents a Hertzian dipole transmitting a Ricker wavelet polarized in the $y$ direction. For all simulations, we defined 30 PML cells at the outer boundary of the $x$ and $z$ directions.

We conducted simulations using two fixed transmission frequencies: 150 and 335 MHz. These frequencies represent the centre frequencies of the two radar profiles located close to the ice core (150 MHz for the EMR profile and 335 MHz for the UWB wideband profile). As input data, we used the conductivity measured by dielectric profiling DEP; along the EDML ice core. We used a constant dielectric permittivity value of ${\mathit{\epsilon }}_{\mathrm{r}}^{\prime }=\mathrm{3.145}$ for the depth range of solid ice. This corresponds to an EM wave propagation velocity of $\sim$ 1.69 $\times$ 10⁸ m s⁻¹, which is the average value determined by at EDML. For the upper $\sim$ 180 m, we fitted a function to the permittivity measured with the DEP until it reached a constant value. The relative magnetic permeability was set to a constant value of ${\mathit{\mu }}_{\mathrm{r}}=\mathrm{1}$, and the magnetic loss factor was set to ${\mathit{\sigma }}_{*}=\mathrm{0}$.

We applied the following signal-processing steps to the simulation output. We derived the reflected energy $P$ of the synthetic radar data following the methods described by and . This involved calculating the envelope of the dielectric field strength ($y$ component of the electric field, $E_y$) using the Hilbert magnitude transform. Additionally, the radar trace was smoothed using a 1-D Gaussian filter. This approach enhances the representation of the impulse response of a transmitted radar signal and reflections from a smooth reflector. The final signal strength was obtained by converting the electric field envelope ($E_y$) to decibel (dB) ($\mathrm{20}{\log }_{\mathrm{10}}\left(\right|E_y\left|\right)$). We did not apply corrections for geometric spreading or englacial attenuation.

For the age assignment, we examined which of the simulated reflections originated from conductivity peaks and matched the IRHs . We then determined the depth of the corresponding conductivity peaks (if present and unambiguous to determine) and used this depth to assign the age via the AICC2023 age–depth chronology . This method is less dependent on the choice of radar wave propagation velocity in ice, firn correction factor, and ice surface determination, relying solely on the assignment of conductivity peaks, and provides a clearer age determination with a smaller associated error. Hence, for the error estimation, we only consider the dating error provided in the AICC2023 age–depth chronology. In addition, it can be used to better estimate the average relative permittivity of the ice at the ice core site.

Figure 4

Location, vertical overlap, and spatial distribution of IRHs. (a) Spatial distribution of IRHs from this study (dark blue) in DML and of IRHs from other studies in this region orange and purple, respectively. (b) Colour-coded overlap of IRHs. (c–e) Magnified views showing the spatial coverage of IRHs in key regions: (c) Jutulstraumen Glacier onset, (d) EDML ice core site, and (e) ice divide southwest of EDML (locations indicated by white boxes in panel b).

[Figure omitted. See PDF]

2.6 IRH depth and age uncertainty quantification

We follow and to quantify the depth and age uncertainty for all IRH depths and consider the following sources of depth uncertainties: (i) the range resolution of the radar system, (ii) the uncertainty associated with the value for the firn correction, and (iii) the uncertainty in the chosen dielectric permittivity. From the obtained total depth uncertainties we derive an age uncertainty by determining the lower and upper age bound of the respective IRH age at the EDML ice core chronology.

The range resolution of the EMR system in the short-pulse mode is approximately 5 m . For the UWB, the range resolution is determined by the bandwidth and the windowing factor in processing what results in a range resolution for the narrowband product of 4.3 m and for the broadband product of 0.35 m . For determining the uncertainty in the firn correction, we use the typically applied values in the range of 10–15 m. The values are uncertain due to the paucity of measurements and can vary spatially due to different rates of firn densifications over large regional scales. Previous studies used uncertainty values ranging from 0.5–2 m in East Antarctica and 3 m in West Antarctica . Therefore, we estimate an uncertainty of $\pm$ 2.5 m for our firn correction. For the uncertainty in the dielectric permittivity of dry ice, we apply an error of 1 % . The total IRH depth uncertainty is the root-sum-square error of the three uncertainties and is calculated for each IRH data point.

Figure 5

Colour-coded absolute depth representation of the nine IRHs from this study. The background map displays the ice surface elevation Reference Elevation Model of Antarctica (REMA) in hillshade. The fine black lines delineate the drainage basins according to .

[Figure omitted. See PDF]

2.7 Data validation

We validated our traced IRHs using a cross-point analysis for each IRH. To achieve this, we calculated crossover points along the traced horizons and determined the depth difference between two intersecting IRHs (see Appendix  for details). This ensured that the same IRHs were traced across intersecting transects across our study area. In total, we analysed more than 2000 crossover points, and the differences for all IRHs in the crossover analysis range between 0 and 74 m, with a mean standard deviation of 3.75 m. Since we traced the ice surface at the highest gradient of surface return power, while the IRHs were traced at their maximum return power, a certain fluctuation in depth at the cross-points is expected. The individual results of the cross-point analysis and details on the method are presented in Appendix .

3 Results

3.1 IRHs ages in EMR short-pulse and UWB data

Two different dating methods provided IRH ages ranging from 4.8 to 91.0 ka (Table ) covering the Holocene and Last Glacial Period. We use the age derived from the depth of conductivity peaks to date the IRHs due to their greater accuracy and only use the age from TWT–depth conversion in the radar profiles when no conductivity peak could be assigned. The radar reflections of seven of nine IRHs could be clearly associated with conductivity peaks in the EDML core using radar forward modelling (Fig. ).

Figure 6

Colour-coded normalized depth representation of the nine IRHs from this study. The background map displays the ice surface elevation REMA in hillshade. The fine black lines delineate the IMBIE drainage basins .

[Figure omitted. See PDF]

3.2 Spatial extent of IRHs in DML

The overall coverage of traceable IRHs in this study spans the entire western DML, covering an area of approximately 450 000 km². IRHs extend over 700 km south of the EDML ice core, $\sim$ 200 km north, and $\sim$ 400 km southwest and northeast. IRHs also extend along the divide between EDML and Dome Fuji ice core sites into a region where IRHs were traced extensively in AWI's EMR long-pulse data Fig. a;. The horizontal density of traced IRHs varies significantly (Fig. c–e) and is highest within a radius of about 30 km around the EDML ice core (Fig. d). Additionally, there is dense coverage in the Jutulstraumen drainage basin south of the Maud Belt, particularly at the onset of the Jutulstraumen ice stream (Fig. c). Furthermore, there is a high density of IRHs southwest of the EDML ice core at the divide between the Slessor, Bailey, and Recovery drainage basins and the drainage basin feeding the Riiser-Larsen and Brunt Ice Shelf (Fig. e). Generally, IRHs cover mostly slow-moving regions, in particular at the ice divides. Regions of higher flow velocities are covered at the onset of Slessor ice stream ($\sim$ 35 m yr⁻¹) and Jutulstraumen Glacier ($<\mathrm{100}$ m yr⁻¹).

IRHs could be consistently traced in most radar profiles of the different radar products (Fig. ). However, there are gaps, particularly in the area of the Jutulstraumen ice stream, which increase downstream as the internal layers become more folded due to ice flow convergence and the deep canyon system . Strongly dipping internal layers show a loss in return power depending on their dip angle at depth due to off-nadir ray path losses and destructive interference . Moreover, gaps in the 4.8 and sometimes in the 7.2 ka IRH occur because the surface multiple reflection overlays these IRHs, which makes continuous tracing difficult.

The IRH overlap is greatest in the vicinity of the EDML ice core, where all nine IRHs could be traced nearly continuously (Fig. b). Here, the layer continuity is least disturbed by ice flow, and both shallow and deep IRHs are well resolved in all radar systems. There is also high IRH overlap at the onset of the Jutulstraumen ice stream and along the ice divide towards Dome Fuji. The lowest overlap is generally found in the part of radar profiles closest to the Maud Belt, where deeper and older IRHs become absent.

In our region, IRHs with ages of 38.1 and 73.4 ka were previously traced and published by and . The most extensive of the two IRH archives, published by , is now significantly expanded by several hundreds of kilometres around EDML. The 38.1 and 73.4 ka IRHs are extended by $\sim$ 450 % and $\sim$ 480 %, respectively. Additionally, seven more IRHs were added with comparable coverage (Fig. b). Based on the existing data set of IRHs from , this represents a total increase in IRH data coverage by a factor of 20.

Figure 7

Depth distribution of the 38.1 ka IRH in the catchment area of the Jutulstraumen (a–c) and upstream of the Slessor ice stream (d–e). Background maps in panels (a) and (d) show ice flow velocity , panels (b) and (e) show bed elevation , and panels (c) and (f) show ice thickness with RACMO (Regional Atmospheric Climate Model) snow accumulation represented as contour lines .

[Figure omitted. See PDF]

3.3 Spatial variation in IRH depth distribution

Our IRHs, which range between 4.8 and 91.0 ka in age, are found at minimum and maximum depths of 200 and 2200 m below the ice surface, respectively (Fig. ). Both the absolute depth below the ice surface and normalized depth distribution of the IRHs show significant variation. We observe a general pattern where IRHs around the ice divide near EDML and to the west are deeper compared to IRHs further south and towards Dome Fuji (Figs.  and ). This pattern is even more pronounced in the normalized IRH depths along the ice divide between EDML and Dome Fuji. Additionally, we note a contrasting trend when comparing absolute and relative IRH depths at the edges of the Maud Belt and on the plateau. For the 73.4 ka IRH, for example, it is evident that northeast of EDML, this IRH lies at depths of approximately 1000 to 1200 m, similar to the southeastern part near Dome Fuji (comparison of Feature A and Feature B in Fig. ). However, the relative depth reveals that the 73.4 ka IRH is about 65 % below the surface northeast of EDML but only about 45 % near Dome Fuji. This pattern is also evident in the 11.1, 25.6, and 38.1 ka IRHs but with smaller differences as the IRHs become younger (comparison of Feature A and Feature B in Fig. ).

The IRH depths within the Jutulstraumen drainage basin are particularly variable and show the following pattern in their relative depths. In the southwestern part of the Jutulstraumen onset region, the IRHs are significantly deeper (absolute depth) than in the northeastern area. The IRHs in this region (denoted as Feature A in Figs.  and ) are significantly shallower compared to the southwest and east but show similar relative depths to EDML. This general pattern is evident in all IRHs that cover this region. A different pattern emerges when considering absolute depths. We find a greater variation with very shallow IRHs in the northeast and very deep IRHs in the Jutulstraumen Trough area. In the 38.1 ka IRH, these depth differences amount to up to 1500 m.

Additionally, the variation in depth of the IRHs is striking in two profiles extending from EDML towards the south-southeast. Around 79.9° S, there is a significant increase in absolute depth followed by a decrease to more shallow depths, observable in all IRHs from 7.2 ka onward (Feature C in Fig. ). This pattern is again more pronounced in older IRHs compared to younger IRHs. However, there is no information on the relative IRH depth at this location since the ice base reflection is not visible in this section of the radar profile. Moreover, we observe an increase in absolute IRH depth at the southwestern edge of our IRH data, which is particularly pronounced in the 25.6 and 38.1 ka IRHs (Feature D in Fig. ). Notably, older IRHs of this particular region could not be traced in the radar data.

4 Discussion

4.1 Interpretation of IRH depth distribution variations

Bed topography undulations, the overall accumulation patterns and complex ice flow dynamics are among the most likely reasons why IRH depths in DML vary widely across different subregions. The overall snow accumulation pattern in East Antarctica , velocity field , and ice thickness distribution agree on the broader scale with IRH depths, e.g. deep IRHs on the plateau and shallower IRHs towards the coast. Similar to the behaviour seen, for example, over the Gamburtsev Mountains or near Jutulstraumen Glacier , IRHs in DML likely show variability in depth where bedrock is uneven, particularly near subglacial ridges and troughs.

Feature A (Figs. , , and ) represents a notable anomaly in a region characterized by particularly shallow IRHs in both absolute and relative depths. This region deviates from the general trend of having a gradient from shallower IRHs on the plateau to deeper IRHs (especially in relative depth) near the Maud Belt. This general trend is likely linked to the overall snow accumulation pattern, with a higher accumulation rate towards the Antarctic coast and a lower accumulation rate on the plateau, as evidenced by the RACMO2 Regional Atmospheric Climate Model; snow accumulation data (Fig. c) and by firn core analysis . The fact that this pattern is also evident in older and thus deeper IRHs suggests that this low-accumulation anomaly may have persisted in this region not only during the Holocene but also during the Last Glacial Period, given that the ice in this region shows very low flow velocities and relatively undisturbed layering throughout the ice column.

Additionally, IRHs are particularly deep (both absolute and relative) in the southwestern area of the Jutulstraumen onset (Figs.  and ), which could be linked to the ice stream activity and bed topography (Fig. a and b). Besides dynamic thinning due to acceleration of ice flow, increased friction at the ice base could lead to basal melting, causing the IRHs to subside. However, our work also shows that ice older than 73.4 ka is present in the main trunk of Jutulstraumen, at least at its onset. Further northeast, the shallow depths of the IRHs might be explained not only by the low-accumulation anomaly but also by reduced basal sliding, less basal melting, and even a positive basal mass balance due to basal freeze-on of subglacial water .

A comparison of IRH depths between Feature A (south of the Maud Belt) and Feature B (near Dome Fuji) highlights the importance of considering the relative depth of IRHs in relation to ice thickness. The fact that all IRHs traced near Dome Fuji are significantly shallower relative to ice thickness than those further northwest can be explained by the ice sheet's geometry and surface mass balance. Ice and firn cores indicate that both regions are characterized by low accumulation rates . This low-accumulation pattern is confirmed in the wider Dome Fuji area via shallow-sounding radar measurements in combination with firn cores and snow pits . However, the ice at Dome Fuji is much thicker and significantly older at greater depths compared to further north. This results in lower relative depths for the IRHs. The low-accumulation pattern at Dome Fuji likely persisted over time (hundreds of thousands of years), impacting the long-term stratigraphy.

The sudden IRH depth increase of all ages along the southward-trending radar profiles (Feature C) and some IRHs in the southwest (Feature D) could be explained by changes in bed elevation and ice dynamics, respectively. Although there are no visible bed reflections in the radar profiles where the IRHs were traced, ice thickness data from nearby regions suggest a connection with the underlying topography. The bed topography data from BedMachine indicate that the IRH depth increase in Feature C is mostly associated with a subglacial valley with a width between 15–25 km and not with snow accumulation or flow velocity anomalies (Fig. d–f). In contrast, the IRHs at Feature D are located at the onset of one of the tributaries of the Slessor ice stream, where faster flow dominates. Higher flow velocity at this point would dynamically thin the ice, which must be compensated by lateral mass influx or increased snow accumulation, which would increase IRH depths. Hence, the relative drop in elevations of the IRHs in this area could indicate the long-term activity of this tributary (Fig. d), as previously suggested by . The fact that older IRHs, such as the 52.2 and 73.4 ka, could not be traced in this area further supports the hypothesis of a significant dynamic component causing IRH subsidence or layer disruption, which could suggest that this ice stream system has been stable at least for the past $\sim$ 70 ka.

4.2 Link between IRHs and volcanic eruptions

The dominant origin of IRHs was attributed to changes in conductivity due to variations in acidity originating from deposits of large volcanic events . Six of the IRHs that we dated using conductivity peaks and one IRH dated using the EDML ice core stratigraphy can be associated with major bipolar volcanic eruptions during the Holocene and Last Glacial Period . Our 11.1 ka IRH, associated with a moderate conductivity peak, can be correlated with the deposits of a bipolar eruption event dated to 11.3 ka. These deposits have been identified in both Antarctic and Greenland ice cores . Additionally, we consider the possibility that our 15.4 ka IRH could be linked to an event dated to 15.6 ka. Although there is no clear and strong conductivity peak in the EDML ice core for this event, one explanation could be that this event is classified as a relatively weaker bipolar eruption compared to others rank 32;.

The radar reflection of the 25.6 ka IRH can be linked to sulfate deposits from the New Zealand Taupo Oruanui eruption, 25.32 ka, which were detected in two Greenland (although undetected in the GISP2 ice core) and three Antarctic ice cores and confirmed by tephra deposits in Antarctica . The IRH we dated to 38.1 ka is very likely associated with the third strongest bipolar eruption in the Last Glacial Period known to date, which occurred at 38.13 ka. The exact eruption location is unknown, but it is likely a low-latitude or Southern Hemisphere eruption with an estimated average climate forcing twice as strong as the Tambora eruption in 1815 CE . Our 52.2 ka IRH could correspond to an eruption that occurred at 52.3 ka and is listed as the fifth strongest eruption in . The deposits from this eruption can be traced in three Greenland and three Antarctic ice cores , and the possible eruption site is estimated to be in the high-latitude Northern Hemisphere (above 40° N).

From tracing IRHs in the radar data in this region, we find that the reflection of the 73.4 ka IRH is the most discernible and extensive in the deeper third of the ice sheet. As per , we also speculate that this reflection most likely corresponds to the deposits of the Toba eruption in northwest Indonesia, dated using tephra and Ar–Ar dating to 73.88 $\pm$ 0.32 ka and 75.0 $\pm$ 0.9 ka, depending on the location of the deposits see references in. The Toba eruption is considered the largest known supereruption of the last 2.5 million years and likely occurred during the cooling transition from Greenland Interstadial 20 to Greenland Stadial 20 and references therein.

Although our IRH ages align closely with major volcanic events, minor age discrepancies are inherent to the dating methods for the IRHs and volcanic events, respectively. Inconsistencies can arise depending on the specific dating method for volcanic events e.g. the different ages for the Toba eruption; as well as due to atmospheric transport and deposition processes of the volcanic ash and sulfate deposits. Additionally, ages assigned to IRHs are core-specific and vary between ice cores .

Figure 8

IRHs from published data sets whose ages match with the (a) DML 4.8 ka horizon, (b) DML 38.1 ka horizon, (c) DML 73.4 ka horizon, and (d) DML 91.0 ka horizon. References to the scientific studies and data publications are provided in Table .

[Figure omitted. See PDF]

Table 3

IRHs from published data sets whose ages match with the IRHs of this study in western DML and likely represent the same physical reflector. The following abbreviations apply to the “Region” column: East Antarctica (EA), Central East Antarctica (CEA), Central West Antarctica (CWA), Dronning Maud Land (DML), Institute and Möller (I&M), and South Pole (SP).

DML IRH age	IRH age	IRH name	Region	Scientific reference	Reference data
4.8 ka	4.7 ka	SPC_IRH0	Titan Dome
	4.7 ka	I1	DML
	4.7 $\pm$ 0.3 kaa	H2	I&M catchment
	4.7 $\pm$ 0.3 ka	R2	Pine Island Glacier
	4.7 $\pm$ 0.2 ka	LM9	Marie Byrd Land
7.2 ka	7.4 ka	I2	DML
	6.9 $\pm$ 0.3 kaa	H3	I&M catchment
	6.9 $\pm$ 0.3 ka	R3	Pine Island Glacier
	6.4 $\pm$ 0.2 ka	Layer no. 4	CWA
11.1 ka	10.7 ka	SPC_IRH1	Titan Dome
	10.5 $\pm$ 0.5 ka	DELORES_IRH_1	Dome C
25.6 ka	25.0 $\pm$ 0.3 ka		Marie Byrd Land
	25.1 ka	I5	DML
38.1 ka	37.8 ka	SPC_IRH4	Titan Dome
	38.1 ka	I6	DML
	38.1 $\pm$ 0.7 ka	OIA_EDC_IRH3	CEA
	39.2 $\pm$ 1.1 ka	DELORES_IRH_3	Dome C
	37.1 ka		CEA
	38.2 $\pm$ 2.2 ka	H1	SP & Dome A
	36.3 $\pm$ 3.6 ka	H2	Dome F
	38.2 $\pm$ 0.6 kab	H1	CEA
52.2 ka	51.4 ka	SPC_IRH5	Titan Dome
73.4 ka	72.5 ka	SPC_IRH6	Titan Dome
	72.4 ka	I8	DML
	73.4 $\pm$ 2.1 ka	OIA_EDC_IRH5	CEA
	73.3 $\pm$ 2.4 ka	LDC-VHF_IRH_1	Little Dome C
	75.1 $\pm$ 2.2 ka	DELORES_IRH_6	Dome C
	75.3 $\pm$ 7.0 ka	H4	Dome F
	73.7 ka		EA
91.0 ka	93.9 ka	SPC_IRH7	Titan Dome
	90.0 $\pm$ 2.1 ka	LDC-VHF_IRH_3	Little Dome C
	90.4 $\pm$ 3.6 ka	H2	SP & Dome A
	90.2 $\pm$ 1.6 kab	H5	CEA

^a The age and error for the respective IRHs from are based on the dating at the WAIS Divide ice core by , with vertical conformity between the two data sets shown in . ^b The age and error for the respective IRH from are based on the dating at the Dome C ice core.

Figure 9

Potential to extrapolate IRHs to AWI EMR long-pulse data. (a) Overview of EMR short-pulse data (red) superimposed with the IRHs traced in this study (black) and in the background showing the EMR long-pulse data (dark blue). The IRHs from are depicted in purple. The lines A–A^′ (black) and B–B^′ (red) in panel (a) indicate a long-pulse and short-pulse radargram at Dome Fuji (Fig. ) where IRHs from this study and IRHs from intersect. (b) Composite radargram (automatic gain control; AGC version) from EMR short-pulse and long-pulse data (profile ID 19992116, long-pulse, right; 19993116, short-pulse, left). At the intersection, transitions where IRHs from this study are also represented in the long-pulse data are marked with arrows.

[Figure omitted. See PDF]

Figure 10

Illustration of two radargrams (automatic gain control; AGC version) and the IRHs intersecting at the Dome Fuji ice core. Panel (a) shows the radargrams with IRHs, and panel (b) shows the radargrams without IRHs but where the IRHs intersect this study and 's IRHs at the Dome Fuji ice core site. (a) The radargram on the left is an EMR long-pulse profile with IRHs from , while the radargram on the right is an EMR short-pulse profile with IRHs from this study together with the IRHs from . (b) The location at the Dome Fuji core, where the 36.3 and 75.3 ka IRH from intersects with the 38.1 and 73.4 ka IRH from this study and are marked with arrows. The colours of the IRHs in the left radargram correspond to those in Fig. 2 in . The location of the two profiles is shown in Fig. a.

[Figure omitted. See PDF]

4.3 Comparison with IRHs from other studies

traced eight IRHs along the ice divide connecting EMR short-pulse profiles between EDML and Dome Fuji and dated them via TWT-to-depth conversion using the EDML AICC2012 age–depth chronology . Several of these isochrones correspond with those we have traced over a much larger spatial extent in this study corresponding IRH ages dated at the EDML core byshown in parentheses: 4.8 ka (4.7 ka), 7.2 ka (7.4 ka), 25.6 ka (25.1 ka), 38.1 ka (38.1 ka), and 73.4 ka (72.4 ka) (Table ). The discrepancies arise due to a combination of the choice of radar profiles, different dating methods, the type of surface pick, and the selected propagation speed for radar waves in ice, all of which are within the stated error margins. Moreover, the IRHs in were also dated using the Dome Fuji age–depth chronology, resulting in different ages. In addition to the study from , the 38.1 and 73.4 ka IRHs had also already been traced extensively in our study area by (Fig.  and Table ), and their paper also discusses briefly the spatially limited presence of IRHs with ages of 4.8, 7.6, 15.4, and 25.0 ka. Together, these data provide an essential foundation for this study and enable a direct and independent comparison of both tracing and dating methods.

A comparison with traced IRHs in East Antarctica, in the region of the Vostok ice core and EPICA Dome C ice core , suggests that some of these IRHs likely correspond to those identified in our study (Table ). Most of these isochrones extend far beyond an age of 100 ka, reaching up to approximately 700 ka, with the dating of the oldest IRHs being particularly challenging due to the inherent uncertainties associated with the 1-D age–depth model used to date . The 38.1 and 73.4 ka IRHs dated to the same ages by are almost certainly the same reflections in this region, as they appear as very strong reflections in our radar data and correlate with significant past volcanic eruptions in ice cores . Furthermore our 11.1, 38.1 and 74.3 ka IRHs also likely correspond to the 10.5, 39.2, and 75.2 ka IRHs from at Dome C traced in ground-based DELORES (Deep Looking Radio Echo Sounder) data acquired in the 2016/17 and 2017/18 seasons as well as two IRHs (dated to 73.7 and 90.3 ka) traced in radar data acquired with the UA LDC-VHF system at Little Dome C Table .

Our isochrones also extend into the area around Dome Fuji, for which there are dated IRHs from (Fig. ). However, these were traced using EMR long-pulse data, which have a vertical resolution of $\sim$ 50 m, making the direct comparison more difficult. Nonetheless, they also identified IRHs dated to 36.3 $\pm$ 3.6 and 75.3 $\pm$ 7.0 ka, which, when examined at the intersection points with our EMR short-pulse radar data, most likely represent the same IRHs (Figs.  and ). A large amount of our IRH ages correlate with IRHs from at Titan Dome, particularly their 4.7, 10.7, 37.8, 51.4, 72.5, and 93.9 ka IRHs. Moreover, we see a connection of our 38.1 and 91.0 ka IRHs to the 38.2 $\pm$ 2.0 and 90.4 $\pm$ 3.6 ka IRHs observed in radar data between South Pole and Dome A acquired with the British Antarctic Survey's Polarimetric radar Airborne Science Instrument (PASIN and PASIN 2) during the AGAP and PolarGAP surveys (Fig.  and Table ).

Dated IRHs with similar ages to ours are not confined to East Antarctica. Several studies have traced, dated, and published IRHs in West Antarctica . In particular, all of these studies identified a significantly widespread IRH dated at $\sim$4.7 ka that is present across most of West Antarctica e.g. and which matches our shallowest IRH over DML (Fig.  and Table ). For instance, identified a 4.7 ka IRH across West Antarctica's ice divide and throughout the Thwaites Glacier catchment, and identified a 4.72 ka IRH across the divide and throughout the Pine Island Glacier catchment, with direct links showing vertical conformity with another IRH widely traced across the Institute and Möller ice stream catchment by . Additionally, there may be further links between our 7.2 ka IRH and the 6.94 ka IRH from and , respectively, and our 25.6 ka IRH and the 24.9 ka IRH from (Fig.  and Table ). However, to accurately determine if these dated isochrones indeed originate from the same reflector Antarctic-wide, it would be necessary to link not only the IRH ages (which can vary depending on the ice core, method, and determined TWT between the ice surface and the IRH) but also the backscatter patterns in the radar data. Here, radar lines connecting data sets between East and West Antarctica would be beneficial (Fig. ). However, this presents challenges, such as long flights over the Antarctic Plateau, and geographical difficulties, such as the loss of layer continuity across the Transantarctic Mountains. For a definitive connection of IRHs that relate to the same reflectors and, therefore, the same deposits, geochemical analyses of ice cores could provide additional constraints.

The AntArchitecture initiative has been aiming to compile all such IRH data sets into a continent-wide radiostratigraphic database, which will be used to create a three-dimensional age–depth model of the ice sheet by interpolating between known IRH along-track data sets similar to already existing gridded topographic models such as Bedmap and BedMachine Antarctica . For such a task to be successful in the future, additional IRH tracing will need to take place in poorly surveyed areas of the ice sheet or away from ice divides and ice core locations using existing RES coverage to avoid over-fitting the interpolation algorithms. By utilizing 3 decades of AWI RES data from multiple radar systems in the area, this study connects a significant gap in the East Antarctic radiostratigraphic record for IRHs older than 7.2 ka and provides additional evidence of the presence of West Antarctica's 4.7–4.8 ka IRH deep in the East Antarctic ice sheet. The success of identifying similar IRHs across large swaths of both East and West Antarctica as demonstrated here (Fig. ) invites further opportunities to make use of the available historical record of yet unused radar data across Antarctica to accelerate the retrieval of a continent-wide age–depth model of the ice sheet, as envisaged by AntArchitecture .

4.4 Potential to expand IRH coverage in and beyond DML

In our study, we have not fully utilized the complete archive of AWI radar data. We focused on EMR short-pulse data, UWB narrow data, and UWB wideband data because these have similar vertical resolutions (5 m or better), enabling the capture of a higher number of IRHs as well as shallow IRHs that are comparably represented across these systems. However, there are additional areas in DML covered exclusively by EMR long-pulse data (approximately 50 m vertical resolution but deep sounding). At intersections between EMR short-pulse and long-pulse data where we were able to trace IRHs, it becomes apparent that the 25.6, 38.1, 52.2, and 73.4 ka IRHs are also represented in the long-pulse data (Fig. ), which aligns with the IRHs of similar age in (Fig.  and Table ). However, the range resolution of the long-pulse system used in is significantly lower, resulting in a larger uncertainty in age dating compared to our IRHs. A direct comparison of two radargrams intersecting at Dome Fuji confirms that particularly the 36.3 and 75.3 ka IRHs from correspond to the same reflections identified in our study, which we date as 38.1 and 73.4 ka, respectively (Fig. ). Additionally, it is likely that the reflection of the 95.5 ka IRH from corresponds to our 91.0 ka IRH (Fig. ). However, the offset in IRH depth and age is relatively large.

The representation of the same reflectors in EMR long-pulse and short-pulse data does not only account for reflections caused by changes in electrical conductivity but also for reflections caused by changes in ice crystal lattice orientation . Considering the AWI EMR radar line coverage in Fig. a, there is additional potential, particularly for regions upstream of Princess Ragnhild and Princess Astrid coasts, as well as the onsets of the Slessor and Recovery ice streams. Additionally, the spatial coverage of the East Antarctic Plateau between EDML and Dome Fuji could be significantly improved. It may even be possible to establish connections between our IRHs and the western Coats Land at the ice divide to the Bailey and Slessor basins using AWI's and the British Antarctic Survey's radar data , which are of similar range resolution.

4.5 Significance of dated IRHs for reconstructing the ice-sheet history in western DML

Reconstructing past ice-sheet configurations is crucial for understanding ice-sheet processes and their impacts on the Earth system, as well as providing essential constraints for ice-sheet models. The spatial depth distribution of our nine IRHs in western DML can significantly contribute to understanding ice-sheet evolution. By comparing IRH depths and spatial anomalies, we can test hypotheses about underlying mechanisms, such as surface or basal mass balance and ice-dynamic processes, as well as how they may have changed over time. Our dated IRHs also provide an opportunity to validate or expand existing findings and hypotheses on ice surface changes based on cosmogenic nuclides about past ice dynamics in DML . In addition, the absolute and relative depth of our IRHs offer the opportunity to complement studies on, for example, preserved paleo-geomorphological structures e.g. or regions of basal freeze-on of subglacial water e.g. providing further insights into past flow behaviour and ice-sheet configurations.

Furthermore, we see significant potential for improving regional and Antarctic-wide ice-sheet models to provide better projections for future sea-level rise . Indeed, IRHs have previously been used as a tool for calibrating results in model simulations and would be well suited to testing hypotheses of the stability of the dynamic Slessor and Jutulstraumen ice stream onset regions over the last interglacial period . Conversely, IRH data from less dynamically active regions, such as the ice divide southwest from EDML, along the divide between EDML and Dome Fuji, the region south of EDML, and the shallow-IRH anomaly region (Feature A), are well suited for inferring the spatial distribution of past accumulation rates .

An additional approach to decipher the past ice dynamics of a region using IRHs lies in the structural analysis of three-dimensional englacial structures, particularly the geometries of folds and fold axial planes, as well as the development of fold amplitudes with depth . Especially for Greenland, it has been shown in various dynamic ice regions that a three-dimensional representation aiming to depict deformation structures like folds accurately can contribute to understanding ice mechanical properties and past flow patterns . Notably, the profiles in the upstream region of the Jutulstraumen ice stream, as well as the region around EDML, where the high density of IRHs allows resolving small-scale englacial structures, and the profiles at the ice divides have the potential to provide insights into past ice dynamics, such as ice flow direction changes or shifts in the location of ice divides.

5 Conclusions

We have mapped nine internal reflection horizons (IRHs) using radar data from various systems of the Alfred Wegener Institute, collected over the last 3 decades. Our IRHs cover the western DML, south of the Maud Belt, and range from 4.8 to 91.0 ka, covering the Holocene and Last Glacial Period. Accurate dating of the IRHs was achieved by combining the age–depth chronology of the EPICA DML ice core and DEP-based radar forward modelling. Additionally, six of the new IRHs could be linked to deposits from significant historic bipolar volcanic eruptions, facilitating synchronization of these reflectors across Antarctica and potentially Greenland. Many of the IRHs mapped in this study likely correspond to the same englacial reflectors found in extensive regions over East and West Antarctica. A comparison with AWI EMR long-pulse data suggests that some of the IRHs identified here could be extended to even larger sectors of East Antarctica.

Our results show that the influence of bed topography, accumulation rates, and ice dynamics is engraved in the age–depth stratigraphy of the ice in the western DML. Thus, our study significantly contributes to a broader understanding of the englacial age architecture of the DML region in East Antarctica and highlights the potential for linking individual IRHs to other regions in Antarctica. Furthermore, the findings presented here are fundamental to enhancing our comprehension of past ice-sheet processes and are crucial for numerical ice flow models aimed at improving our understanding of the paleoclimatic history of Antarctica.

Appendix A Cross-point analysis

We performed an analysis of IRH depth differences at crossover points of intersecting radar lines to validate that we assigned the same reflections for the respective IRHs. For the analysis of crossover differences in the respective IRHs, we calculate the exact intersection points between the geometric lines formed by connecting the individual IRH picks. We create line subsets of one radar profile where IRH gaps exist to avoid creating lines, and therefore potential fake intersection points, where no IRHs are traced within a radar profile. Using these lines, we determine the location of the resulting intersection points. We then create a circular buffer around these intersection points and capture all IRH picks within this buffer. We select a buffer radius of 50 m and average the IRH depths of all picks within a profile to minimize small-scale variations introduced by the semi-automatic tracker. The difference is derived from the mean depths of picks from both intersecting lines within the buffer.

Examining the histograms in Fig. , we observe that most intersection point differences are below 10 m depth difference, and the depth difference generally decreases for higher values. For the 4.8 and 7.2 ka IRHs, we see a small local maximum around a 5 m intersection point difference. Moreover, the differences become larger for older IRHs, leading to a higher number of larger cross-point errors compared to younger IRHs.

Table A1

IRH data set description for column-separated text files.

Column name	Unit	Description
Season		The year in the season column represents the second year of the respective season (e.g. antr2001 represents the Antarctic season 2000/01)
Profile ID		Radar profile ID
Paradigm ID		Radar profile ID in the Paradigm system (AWI picking software)
Device		EMR short-pulse, UWB narrowband, or UWB wideband
Trace		Trace number
Longitude	decimal degree	EPSG:4326
Latitude	decimal degree	EPSG:4326
IRH age [ka]	thousand years
IRH age uncertainty [ka]	thousand years	Estimated IRH age uncertainty at the respective IRH depth
TWT IRH [s]	seconds	Absolute TWT in radargram
TWT surface [s]	seconds	Absolute TWT in radargram
TWT ice base [s]∗	seconds	Absolute TWT in radargram
Depth IRH [m]	metres	Depth below ice surface
Depth IRH uncertainty [m]	metres
Surface elevation [m]	metres	Ice surface elevation based on REMA;
Ice base depth [m]∗	metres	Depth below ice surface

^∗ The ice base reflection is not present in all radar profiles.

Figure A1

Differences of IRH depths up to 50 m at crossing points for the respective IRHs. The histograms in the lower-left corner of each panel show the total count of differences of IRH depths at crossing points up to 50 m ($x$ axis). The $y$ axis shows the number of crossing points, and the bin width is 1 m.

[Figure omitted. See PDF]

Appendix B IRH data set description

We publish our dated IRHs on the PANGAEA data publisher in a column-separated text file format, and the key information is summarized in Table . It is worth noting that the vertical position of the IRHs is provided in different formats: (1) as absolute TWT in the radargram, representing the direct reference to the radargram, and (2) as TWT (in seconds) and depth (in metres) beneath our ice surface pick. To ensure this reference is clear for future users, we also provide the absolute TWT of the ice surface pick, the elevation of the ice surface based on the REMA surface DEM version 2; at each location, and the depth of the ice base if the bed reflection is visible in the radargram. For conversion between TWT to depth, we used a dielectric permittivity of ${\mathit{\epsilon }}_{\mathrm{r}}^{\prime }=\mathrm{3.15}$ and a firn correction term of $Z_{\mathrm{firn}}=\mathrm{13}$ m (see Eq. ). We also provide the depth uncertainty of each IRH depth as well as an estimate for the resulting IRH age error (see “Data and methods” section). The “profile_id” is the ID of the radargram, which also appears in the filename. The “paradigm_id” is an internal AWI ID for referencing the AWI radar database and radar processing software where the picks are created and archived.

Data availability

Our IRHs are available on PANGAEA : 10.1594/PANGAEA.973266. The AICC2023 EDML ice core age–depth chronology is available on PANGAEA: 10.1594/PANGAEA.961019. Ice surface velocities from are available from the National Snow and Ice Data Center (NSIDC), 10.7280/D10D4Z . The IMBIE drainage system boundaries can be obtained from http://imbie.org/imbie-3/drainage-basins/ (last access: 17 January 2025). The BedMachine Antarctica v03 bed topography data from are available at 10.5067/FPSU0V1MWUB6 . The Reference Elevation Model of Antarctica (REMA) is available from the Polar Geospatial Center at https://www.pgc.umn.edu/data/rema/ (last access: 17 January 2025) and 10.7910/DVN/EBW8UC . RACMO2.3 snow accumulation data can be obtained at the Institute for Marine and Atmospheric research Utrecht at Utrecht University: https://www.projects.science.uu.nl/iceclimate/models/racmo-data.php (last access: 17 January 2025) and 10.5281/zenodo.10854319 . For the data availability of IRHs used for comparison, we refer to Table .

Author contributions

SF conceptualized the study, wrote the manuscript, traced and dated all IRHs, performed all data analysis, created all figures, and validated the data. AMZ, DS, and SF acquired radar data from the 2023/24 AWI UWB campaign, while AMZ processed the data from this campaign with contributions of SF and VH. AMZ, VH, and SF contributed to tracing the ice–bed interface. JAB provided context of comparable studies and for the discussion about connecting IRHs in this study over Antarctica. VH developed an algorithm for standard ice surface tracking and processed GPS data. OE, DS, SF, and VH implemented access to and usability of AWI's radar data archive. All authors jointly revised and edited the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We thank the AWI polar aircraft technicians Eduard Gebhard and Christoph Petersen for their support in the field during the 2023/24 radar campaign, the Kenn Borek crew of AWI's polar research aircraft, and all supporters of previous radar campaigns mentioned in acknowledgements of related publications. We acknowledge support via airborne radar campaign funding grants AWI_PA_02004, AWI_PA_02020, AWI_PA_02084, AWI_PA_02105, and AWI_PA_02146. Logistical support in the field over the past 3 decades has been provided by Neumayer III station (Germany), Troll station (Norway), Kohnen station (Germany), Princess Elisabeth station (Belgium), and Novolazarevskaya airbase (Russia). We acknowledge the use of software from Open Polar Radar generated with support from the University of Kansas; NASA grants 80NSSC20K1242 and 80NSSC21K0753; and NSF grants OPP-2027615, OPP-2019719, OPP-1739003, IIS-1838230, RISE-2126503, RISE-2127606, and RISE-2126468. The authors would like to thank Aspen Technology, Inc. for providing software licenses and support. This study was motivated by the AntArchitecture Scientific Committee on Antarctic Research Action Group (https://scar.org/science/cross/antarchitecture, last access: 17 January 2025).

Steven Franke was funded by the Walter Benjamin Programme of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation; project no. 506043073). Alexandra M. Zuhr was funded by the DFG in the framework of the priority programme SPP 1158 “Antarctic Research with Comparative Investigations in Arctic Ice Areas” (grant no. 522419679). Julien A. Bodart acknowledges funding from the Swiss National Science Foundation (grant no. 211542).

Financial support

This research has been supported by the Deutsche Forschungsgemeinschaft (grant nos. 506043073 and 522419679) and the Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (grant no. 211542).This open-access publication was funded by the Open Access Publication Fund of the University of Tübingen.

Review statement

This paper was edited by Joseph MacGregor and reviewed by Rebecca Sanderson and Marie G. P. Cavitte.