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

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 (IPCC, 2023). 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 (Siegert et al., 2004; Bons et al., 2016; Leysinger Vieli et al., 2018; Jansen et al., 2024), offering insights into atmosphere–ice–ocean interactions (Drews et al., 2020; Višnjević et al., 2022) 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. Leysinger Vieli et al., 2004; Sutter et al., 2021; Bodart et al., 2023).

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 (Mouginot et al., 2019a) superimposed on the hillshaded surface elevation (Howat et al., 2019). The IMBIE (Ice sheet Mass Balance Inter-comparison Exercise) drainage basins (Rignot et al., 2019) are marked by thin black lines, and deep ice core sites are marked with a white circle.

Radar data provide the opportunity to study the internal architecture of ice sheets (i.e. their stratigraphy) via continuous stratified englacial reflections (Bingham et al., 2024). These internal reflection horizons (IRHs) occur when radar waves are reflected at boundary layers where the dielectric properties of the ice change (Glen and Paren, 1975; Robin et al., 1977). Dielectric contrasts in ice sheets are primarily caused by variations in density (Robin et al., 1969), electrical conductivity (Paren and Robin, 1975), and changes in the preferred orientation of ice crystals (Harrison, 1973; Fujita et al., 1999). 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 (Millar, 1981; MacGregor et al., 2015; Winter et al., 2019). Strong conductivity contrasts are often linked to deposits of acidic materials from major volcanic eruptions (Hammer, 1980; Millar, 1981). 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 (Bingham et al., 2024).

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 (Leysinger Vieli et al., 2004; Eisen et al., 2005; Leysinger Vieli et al., 2011; Cavitte et al., 2018; Bodart et al., 2023; Koch et al., 2023). Since many IRHs can be dated with ice cores, they can also serve as calibration points for ice-sheet models (Sutter et al., 2021; Višnjević et al., 2022) and to model basal conditions (Leysinger Vieli et al., 2018; Fudge et al., 2023; Chung et al., 2023b; Wang et al., 2023a). 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 (Siegert et al., 2004; Bons et al., 2016; Franke et al., 2022a; Jansen et al., 2024), as well as ice-dynamic processes causing IRH discontinuity (Panton and Karlsson, 2015; Sanderson et al., 2023b).

The generation and compilation of IRHs using radar data in Antarctica and Greenland remain largely a manual task (Bingham et al., 2024). While automated methods can successfully trace IRHs to a certain extent (see e.g. review by Moqadam and Eisen, 2024), 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 (Leysinger Vieli et al., 2011; Winter et al., 2019; Cavitte et al., 2016, 2021; Beem et al., 2021b; Wang et al., 2023a; Sanderson et al., 2024), West Antarctica (Siegert et al., 2005; Jacobel and Welch, 2005; Ashmore et al., 2020a; Bodart et al., 2021b; Muldoon et al., 2018), and Greenland (Karlsson et al., 2013; MacGregor et al., 2015; Franke et al., 2023a; MacGregor et al., 2025). 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 (Bingham et al., 2024).

Table 1 Alfred Wegener Institute (AWI) radar system specifications. The range resolution applies for a dielectric permittivity of ε_r=3.15 in ice.

^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. 1). 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 1). These radar products were obtained during campaigns conducted by AWI between 1996 and 2023 with two different radar systems (Fig. 1). The radar data were collected using AWI's polar aircraft Polar 2, Polar 5, and Polar 6 (Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, 2016a), operating between the stations Neumayer and Kohnen (Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, 2016b), 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 (Nixdorf et al., 1999). For IRH tracing, we utilize EMR radar data from nine Antarctic seasons spanning from 1996/97 to 2016/17 (Fig. 1 and Table 1). 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 Wang et al., 2023a). 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 ∼35 and ∼50 m, respectively (Steinhage et al., 2001; Steinhage, 2001; Wang et al., 2023a).

EMR survey lines cover the entire western DML (Fig. 1). Surveys between 1997 and 1999 focused on the pre-site survey of the EDML ice core site (Steinhage et al., 2001; Steinhage, 2001), while surveys between 2001 and 2011 were used to investigate the ice internal structure (Eisen et al., 2006, 2007; Drews et al., 2009; Steinhage et al., 2013) as well as ice thickness for solid earth geophysical studies in DML (Riedel et al., 2012; Eisermann et al., 2020). The survey in 2016/17 extensively surveyed ice thickness at the Dome Fuji ice core site (Karlsson et al., 2018; Wang et al., 2023a).

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 (Rodriguez-Morales et al., 2013; Hale et al., 2016) and has been in operation since 2016 (Kjær et al., 2018). 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 (Franke et al., 2022b). Standard processing techniques were performed with the OPR Toolbox (Open Polar Radar Toolbox, formerly termed CReSIS Toolbox; Open Polar Radar, 2023). 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 (Franke et al., 2021, 2024b). AWI UWB data in the 2023/24 season were acquired with a frequency range between 150–520 MHz (wideband), with a range resolution of ∼ 35 cm. The acquisition settings were chosen to map near-surface IRHs at high spatial resolution (Koch et al., 2023) 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 Winter et al. (2019) 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 Winter et al. (2019) 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 Franke et al. (2021). Ice base determination from the AWI UWB 2023/2024 season was generated using the OPR Toolbox (Open Polar Radar, 2023) and Paradigm. For the EMR data, we were able to partially rely on existing data (Steinhage et al., 2001; Riedel et al., 2012), 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.

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. Steinhage et al., 2001; MacGregor et al., 2015; Cavitte et al., 2016; Winter et al., 2017; Ashmore et al., 2020a; Bodart et al., 2021b) using

$$\begin{array}{}\text{(1)}& Z_{\mathrm{irh}}={\displaystyle \frac{t{c}_{\mathrm{air}}}{\mathrm{2}\sqrt{{\mathit{\epsilon }}_{\mathrm{r}}^{\prime }}}}+Z_{\mathrm{firn}},\end{array}$$

where c_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 (Eisen et al., 2006), 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 (Steinhage, 2001), 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. Rotschky et al., 2007), 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 (Morlighem et al., 2020) 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. Cavitte et al., 2016; Bodart et al., 2021b) and (2) matching reflections of our IRHs with simulated radar reflections based on conductivity peaks derived from measured DEP (dielectric profiling) data (e.g. Eisen et al., 2003, 2006; Winter et al., 2017). 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 (Eisen et al., 2003); 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.

^a Large Northern Hemisphere higher-latitude eruption (rank 18 in Lin et al., 2022).
^b Large lower-latitude or Southern Hemisphere eruption (rank 32 in Lin et al., 2022).
^c Taupo Oruanui eruption (Dunbar et al., 2017).
^d Large lower-latitude eruption (rank 3 in Lin et al., 2022).
^e Large bipolar eruption (rank 5 in Lin et al., 2022).
^f Toba eruption dated to 75.0 ± 0.9 ka (see references in Svensson et al., 2013).
^g Toba eruption dated to 73.88 ± 0.32 ka (see references in Svensson et al., 2013).

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 (Bouchet et al., 2023), 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 Eisen et al. (2006) of ∼ 1.69 × 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 (Steinhage, 2001; Steinhage et al., 2001). 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 % (Bohleber et al., 2012). 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 (Eisen et al., 2006; Winter et al., 2017; Mojtabavi et al., 2022), we used the open-source software gprMax (Giannopoulos, 2005; Warren et al., 2016) to simulate electromagnetic waves in three-dimensional space. This software solves Maxwell's equations using the finite-difference time-domain (FDTD) method (Taflove and Hagness, 2005) and Yee cells (Yee, 1966), which has proven effective for 1-D and 2-D simulations in ice sheets and glaciers (Franke et al., 2023b; Santin et al., 2023). gprMax uses the same algorithm as already successfully employed for ice-core-based forward modelling by Eisen et al. (2006) but includes an efficient solution for model boundary conditions using perfectly matched layers (PMLs; Giannopoulos, 2012) 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 (Morlighem et al., 2020) 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.

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 Δ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 µ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; Moore, 1993; Wilhelms et al., 1998) 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 ∼ 1.69 × 10⁸ m s⁻¹, which is the average value determined by Eisen et al. (2006) at EDML. For the upper ∼ 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 μ_r=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 Winter et al. (2017) and Franke et al. (2023b). 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 (Eisen et al., 2006). 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 (Bouchet et al., 2023). 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 (Winter et al., 2019; Wang et al., 2023a, 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).

2.6 IRH depth and age uncertainty quantification

We follow Cavitte et al. (2016, 2021) and Winter et al. (2019) 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 (Nixdorf et al., 1999). 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 (Franke et al., 2022b) and for the broadband product of 0.35 m (Koch et al., 2023). 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 (Cavitte et al., 2016; Winter et al., 2019; Wang et al., 2023a) and 3 m in West Antarctica (Ashmore et al., 2020a; Bodart et al., 2021b). Therefore, we estimate an uncertainty of ± 2.5 m for our firn correction. For the uncertainty in the dielectric permittivity of dry ice, we apply an error of 1 % (Bohleber et al., 2012). 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) (Howat et al., 2019) in hillshade. The fine black lines delineate the drainage basins according to Rignot et al. (2019).

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 A 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 A.

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 2) 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. 2).

Figure 6 Colour-coded normalized depth representation of the nine IRHs from this study. The background map displays the ice surface elevation REMA (Howat et al., 2019) in hillshade. The fine black lines delineate the IMBIE drainage basins (Rignot et al., 2019).

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, ∼ 200 km north, and ∼ 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. 4a; Wang et al., 2023a). The horizontal density of traced IRHs varies significantly (Fig. 4c–e) and is highest within a radius of about 30 km around the EDML ice core (Fig. 4d). 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. 4c). 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. 4e). 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 (∼ 35 m yr⁻¹) and Jutulstraumen Glacier (<100 m yr⁻¹).

IRHs could be consistently traced in most radar profiles of the different radar products (Fig. 3). 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 (Franke et al., 2021). 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 (Holschuh et al., 2014). 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. 4b). 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 Winter et al. (2019) and Steinhage et al. (2013). The most extensive of the two IRH archives, published by Winter et al. (2019), is now significantly expanded by several hundreds of kilometres around EDML. The 38.1 and 73.4 ka IRHs are extended by ∼ 450 % and ∼ 480 %, respectively. Additionally, seven more IRHs were added with comparable coverage (Fig. 4b). Based on the existing data set of IRHs from Winter et al. (2019), 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 (Mouginot et al., 2019a), panels (b) and (e) show bed elevation (Morlighem et al., 2020), and panels (c) and (f) show ice thickness (Morlighem et al., 2020) with RACMO (Regional Atmospheric Climate Model) snow accumulation represented as contour lines (Lenaerts et al., 2013).

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. 5). 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. 5 and 6). 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. 5). 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. 6).

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. 5 and 6) 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. 5). 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. 5). 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 (Lenaerts et al., 2013), velocity field (Mouginot et al., 2019a), and ice thickness distribution (Morlighem et al., 2020) 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 (Sanderson et al., 2024) or near Jutulstraumen Glacier (Franke et al., 2021), IRHs in DML likely show variability in depth where bedrock is uneven, particularly near subglacial ridges and troughs.

Feature A (Figs. 5, 6, and 7) 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; Lenaerts et al., 2013) snow accumulation data (Fig. 7c) and by firn core analysis (Rotschky et al., 2007). 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. 5 and 6), which could be linked to the ice stream activity and bed topography (Fig. 7a 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 (Franke et al., 2021, 2024b).

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 (Rotschky et al., 2007; Oyabu et al., 2023). 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 (Fujita et al., 2011; Van Liefferinge et al., 2021; Cavitte et al., 2023). 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 (Morlighem et al., 2020) 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. 7d–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. 7d), as previously suggested by Rippin et al. (2006). 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 ∼ 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 (Millar, 1981). 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 (Lin et al., 2022; Svensson et al., 2013). 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 (Lin et al., 2022). 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; Lin et al., 2022).

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 (Lin et al., 2022) and confirmed by tephra deposits in Antarctica (Dunbar et al., 2017). 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 (Lin et al., 2022). 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 Lin et al. (2022). The deposits from this eruption can be traced in three Greenland and three Antarctic ice cores (Lin et al., 2022), 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 Winter et al. (2019), 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 ± 0.32 ka and 75.0 ± 0.9 ka, depending on the location of the deposits (see references in Svensson et al., 2013). The Toba eruption is considered the largest known supereruption of the last 2.5 million years (Chesner, 2012) and likely occurred during the cooling transition from Greenland Interstadial 20 to Greenland Stadial 20 (Lin et al., 2023, 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; Svensson et al., 2013) 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 (Steinhage et al., 2013; Winter et al., 2019).

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 3.

Beem et al. (2021b)Beem et al. (2021a)Steinhage et al. (2013)Ashmore et al. (2020a)Ashmore et al. (2020b)Bodart et al. (2021b)Bodart et al. (2021a)Muldoon et al. (2018)Muldoon et al. (2023)Steinhage et al. (2013)Ashmore et al. (2020a)Ashmore et al. (2020b)Bodart et al. (2021b)Bodart et al. (2021a)Siegert and Payne (2004)Beem et al. (2021b)Beem et al. (2021a)Chung et al. (2023b)Mulvaney et al. (2023)Muldoon et al. (2018)Steinhage et al. (2013)Beem et al. (2021b)Beem et al. (2021a)Steinhage et al. (2013)Cavitte et al. (2021)Cavitte et al. (2020)Chung et al. (2023b)Mulvaney et al. (2023)Leysinger Vieli et al. (2011)Sanderson et al. (2024)Sanderson et al. (2023a)Wang et al. (2023a)Wang et al. (2023b)Winter et al. (2019)Winter et al. (2018)Beem et al. (2021b)Beem et al. (2021a)Beem et al. (2021b)Beem et al. (2021a)Steinhage et al. (2013)Cavitte et al. (2021)Cavitte et al. (2020)Chung et al. (2023b)Chung et al. (2023a)Chung et al. (2023b)Mulvaney et al. (2023)Wang et al. (2023a)Wang et al. (2023b)Leysinger Vieli et al. (2011)Beem et al. (2021b)Beem et al. (2021a)Chung et al. (2023b)Chung et al. (2023a)Sanderson et al. (2024)Sanderson et al. (2023a)Winter et al. (2019)Winter et al. (2018)

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).

^a The age and error for the respective IRHs from Ashmore et al. (2020a) are based on the dating at the WAIS Divide ice core by Bodart et al. (2021b), with vertical conformity between the two data sets shown in Bodart et al. (2021b). ^b The age and error for the respective IRH from Winter et al. (2019) 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 Wang et al. (2023a) 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. 10) where IRHs from this study and IRHs from Wang et al. (2023a) 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 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 Wang et al. (2023a) IRHs intersect this study and Winter et al.'s 2019 IRHs at the Dome Fuji ice core site. (a) The radargram on the left is an EMR long-pulse profile with IRHs from Wang et al. (2023a), while the radargram on the right is an EMR short-pulse profile with IRHs from this study together with the IRHs from Winter et al. (2019). (b) The location at the Dome Fuji core, where the 36.3 and 75.3 ka IRH from Wang et al. (2023a) intersects with the 38.1 and 73.4 ka IRH from this study and Winter et al. (2019) are marked with arrows. The colours of the IRHs in the left radargram correspond to those in Fig. 2 in Wang et al. (2023a). The location of the two profiles is shown in Fig. 9a.

4.3 Comparison with IRHs from other studies

Steinhage et al. (2013) 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 (Bazin et al., 2013). 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 by Steinhage et al., 2013, shown 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 3). 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 Steinhage et al. (2013) were also dated using the Dome Fuji age–depth chronology, resulting in different ages. In addition to the study from Steinhage et al. (2013), the 38.1 and 73.4 ka IRHs had also already been traced extensively in our study area by Winter et al. (2019) (Fig. 8 and Table 3), 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 (Leysinger Vieli et al., 2011; Cavitte et al., 2016, 2021; Winter et al., 2019; Chung et al., 2023b), suggests that some of these IRHs likely correspond to those identified in our study (Table 3). 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 (Cavitte et al., 2021). The 38.1 and 73.4 ka IRHs (dated to the same ages by Cavitte et al., 2021) 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 (Lin et al., 2022; Svensson et al., 2013). 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 Chung et al. (2023b) 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 (Lilien et al., 2021) at Little Dome C (Chung et al., 2023b, Table 3).

Our isochrones also extend into the area around Dome Fuji, for which there are dated IRHs from Wang et al. (2023a) (Fig. 8). However, these were traced using EMR long-pulse data, which have a vertical resolution of ∼ 50 m, making the direct comparison more difficult. Nonetheless, they also identified IRHs dated to 36.3 ± 3.6 and 75.3 ± 7.0 ka, which, when examined at the intersection points with our EMR short-pulse radar data, most likely represent the same IRHs (Figs. 9 and 10). A large amount of our IRH ages correlate with IRHs from Beem et al. (2021b) 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 ± 2.0 and 90.4 ± 3.6 ka IRHs observed in radar data between South Pole and Dome A (Sanderson et al., 2024) acquired with the British Antarctic Survey's Polarimetric radar Airborne Science Instrument (PASIN and PASIN 2) during the AGAP and PolarGAP surveys (Fig. 8 and Table 3).

Dated IRHs with similar ages to ours are not confined to East Antarctica. Several studies have traced, dated, and published IRHs in West Antarctica (Muldoon et al., 2018; Ashmore et al., 2020a; Bodart et al., 2021b, 2023). In particular, all of these studies identified a significantly widespread IRH dated at ∼4.7 ka that is present across most of West Antarctica (e.g. Bodart et al., 2023) and which matches our shallowest IRH over DML (Fig. 8 and Table 3). For instance, Muldoon et al. (2018) identified a 4.7 ka IRH across West Antarctica's ice divide and throughout the Thwaites Glacier catchment, and Bodart et al. (2021b, 2023) 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 Ashmore et al. (2020a). Additionally, there may be further links between our 7.2 ka IRH and the 6.94 ka IRH from Ashmore et al. (2020a) and Bodart et al. (2021b), respectively, and our 25.6 ka IRH and the 24.9 ka IRH from Muldoon et al. (2018) (Fig. 8 and Table 3). 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. 8). 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 (Bingham et al., 2024) 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 (Lythe and Vaughan, 2001; Fretwell et al., 2013; Frémand et al., 2022; Pritchard et al., 2025) and BedMachine Antarctica (Morlighem et al., 2020). 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. 8) 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 (Bingham et al., 2024).

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. 9), which aligns with the IRHs of similar age in Wang et al. (2023a) (Fig. 8 and Table 3). However, the range resolution of the long-pulse system used in Wang et al. (2023a) 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 Wang et al. (2023a) correspond to the same reflections identified in our study, which we date as 38.1 and 73.4 ka, respectively (Fig. 10). Additionally, it is likely that the reflection of the 95.5 ka IRH from Wang et al. (2023a) corresponds to our 91.0 ka IRH (Fig. 10). 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 (Eisen et al., 2007). Considering the AWI EMR radar line coverage in Fig. 9a, 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 (Frémand et al., 2022), 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 (Andersen et al., 2020) about past ice dynamics in DML (Andersen et al., 2023; Braga et al., 2023). 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. Näslund, 1997; Rose et al., 2015; Franke et al., 2021; Carter et al., 2024) or regions of basal freeze-on of subglacial water (e.g. Bell et al., 2011; Leysinger Vieli et al., 2018; Franke et al., 2024b) 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 (Sutter et al., 2021). Indeed, IRHs have previously been used as a tool for calibrating results in model simulations (Sutter et al., 2021; Višnjević et al., 2022) 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 (Rippin et al., 2003, 2006). 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 (Eisen et al., 2005; Huybrechts et al., 2009; Leysinger Vieli et al., 2011; Bodart et al., 2023).

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 (Bons et al., 2016; Franke et al., 2023a; Jansen et al., 2024). 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 (Bons et al., 2016; Zhang et al., 2024) and past flow patterns (Franke et al., 2022a; Jansen et al., 2024). 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. A1, 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.

(based on REMA; Howat et al., 2019)

Table A1 IRH data set description for column-separated text files.

^∗ 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.

Appendix B: IRH data set description

We publish our dated IRHs on the PANGAEA data publisher (Felden et al., 2023) in a column-separated text file format, and the key information is summarized in Table A1. 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; Howat et al., 2019) 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_firn=13 m (see Eq. 1). 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.