Investigating the dynamic history of a promontory ice rise using radar data

Abstract

Ice rises hold valuable records revealing the ice dynamics and climatic history of Antarctic coastal areas from the Last Glacial Maximum to today. This history is often reconstructed from isochrone radar stratigraphy and simulations focusing on Raymond arch evolution beneath the divides. However, this relies on complex ice-flow models where many parameters are unconstrained by observations. Our study explores quad-polarimetric, phase-coherent radar data to enhance understanding near ice divides and domes, using Hammarryggen Ice Rise (HIR) as a case study. Analysing a 5 km profile intersecting the dome, we derive vertical strain rates and ice-fabric properties. These align with ice core data near the summit, increasing confidence in tracing signatures from the dome to the flanks. The Raymond effect is evident, correlating with surface strain rates and radar stratigraphy. Stability is inferred over millennia for the saddle connecting HIR to the mainland, but dome ice-fabric appears relatively young compared to 2D model predictions. In a broader context, quad-polarimetric measurements provide valuable insights into ice-flow models, particularly for anisotropic rheology. Including quad-polarimetric data advances our ability to reconstruct past ice flow dynamics and climatic history in ice rises.

1. Introduction

Ice rises are grounded, locally elevated, ice features surrounded by ice streams or ice shelves. They form over regions with shallower bathymetry, enabling the accumulated ice to stay grounded in these areas. This then results in a locally different flow regime (Matsuoka and others, 2015). Promontory ice rises, such as Hammarryggen Ice Rise (HIR) (Fig. 1), are connected to the main ice sheet via a saddle in the surface topography. They may form triple junctions near their domes (Fig. 1—blue lines), from which three ridges extend into the ice-rise flanks. Ice rises have two main characteristics that make them of particular interest: Firstly, they decelerate ice flux from the main ice sheet toward the ocean and consequently delay grounding-line retreat (Favier and others, 2012, 2014; Favier and Pattyn, 2015; Schannwell and others, 2019; Henry and others, 2022). Secondly, they are an archive for the local atmospheric and ice-dynamic history. The latter is accessible through the englacial stratigraphy, which includes Raymond arches–anticlines in the ice stratigraphy that evolve once a local ice dome or ice divide has formed (Raymond, 1983). The presence or absence of Raymond arches provides insight into the ice-rise history, especially the temporal stability of the configuration, and can be used as a tie-point of the ice thickness to constrain continental ice-flow models. Such tie-points are important, as other constraints, such as exposure dating of rock outcrops (Davies and others, 2012), are unavailable for most of the Antarctic perimeter.

Figure 1. (a) The location of study area in Antarctica. (b) Hammarryggen ice rise, the white contour lines and satellite background represents the surface elevation derived from the REMA dataset (Howat and others, 2022). Two black dashed lines represent the UWB flight lines. The blue lines denote the approximate position of the ridges. The black dot represents the location of the ice core and the red line indicates the pRES profile. (c) The red shading corresponds to the location of the pRES profile. pRES measurement points depicted as red dots in the inset. (d) and (e) A cross-sectional view along the extended pp’ profile, illustrating surface elevation (Howat and others, 2022), bed elevation and ice thickness (Morlighem, 2022).

Much progress in previous studies has guided the interpretation of observed Raymond stacks (i.e. individual Raymond arches and their evolution with depth) in the context of the ice-dynamic history of a respective catchment (Matsuoka and others, 2015). Clear signatures of transience are Raymond stacks that do not align with contemporary ice divides (Nereson and Waddington, 2002), such as at Siple Dome (Nereson and others, 1998). Fully evolved Raymond stacks that align with the contemporary ice divide location are at the other end of the spectrum and indicate stability (e.g. Derwael Ice Rise; Drews and others, 2015). Cases between these two end members (Goel and others, 2020) are more difficult to interpret and require advanced model-data comparison, including thermomechanically-coupled full Stokes models with anisotropic rheology (Martín and others, 2009a, 2009b; Martín and Gudmundsson, 2012) and a dynamically evolving grounding line (Schannwell and others, 2019, 2020; Henry and others, 2022).

A drawback of the model-guided interpretation of observed Raymond stacks is that many unconstrained factors influence the arch amplitude. One of them being the ice anisotropy (Martín and Gudmundsson, 2012; Drews and others, 2015) for which so far virtually no observations away from ice cores were available. This is the main problem that we address in this paper using polarimetric radar as a main tool. Arch amplitude is influenced by multiple interrelated factors that affect the development of ice fabric. Firstly, the degree of non-linearity in Glen's flow law exponent significantly impacts arch size; a higher non-linearity typically results in larger arch amplitudes (Gillet-Chaulet and others, 2011; Martín and others, 2009a, 2009b; Drews and others, 2015; Bons and others, 2018). In contrast, the along-ridge flow component generally produces smaller arches (Martín and others, 2009a, 2009b). Similarly, variations in bed topography can lead to smaller arches when the bed is uneven (Kingslake and others, 2014), while basal sliding also contributes to reduced arch sizes (Petit and others, 2003). Additionally, localized factors such as surface mass balance and erosion at the crest can increase arch amplitudes (Conway and Wilbour, 1999; Drews and others, 2015). The historical thinning or thickening of the ice further translates to changes in arch sizes relative to their current geometry (Martín and others, 2006; Goel and others, 2018).

Ice-core analysis, in combination with shallow and deep radar, can constrain the three-dimensional ice geometry (Hindmarsh and others, 2011) and the surface accumulation history (Philippe and others, 2016; Goel and others, 2017; Cavitte and others, 2022). Strain measurements such as the coffee-can method (Hamilton and Whillans, 2000) and repeat surveys with phase-coherent radar can provide additional constraints on the vertical strain rates (Kingslake and others, 2014). However, other factors, such as ice anisotropy, remain unconstrained, resulting in ambiguous matching of observed Raymond arch stacks with ice-flow models (Drews and others, 2015). Consequently, so far, ice rises and their inferred dynamic history play a minor role in constraining larger-scale ice flow models (Bentley and others, 2014).

Phase coherent radar polarimetry using a ground-based phase-sensitive Radio Echo Sounder (pRES) (Brennan and others, 2014) has seen much development in terms of inferring ice-fabric types for various flow regimes using the polarimetric coherence phase as a metric to extract information from the birefringent radar backscatter (Dall, 2010; Jordan and others, 2019, 2020; Ershadi and others, 2022; Rathmann and others, 2022; Zeising and others, 2023). Anisotropic ice-flow models of steady-state ice rises, as detailed in studies by Martín and others (2009a, 2009b); Martín and Gudmundsson (2012), predict significant gradients in ice-fabric types on either side of an ice divide. This prediction, highlights the impact of anisotropic rheology on ice dynamics. However, thus far, it has not been directly compared with observations, a gap that warrants attention in the field.

Here, we investigate to what extent ice-fabric properties can be derived from quad-polarimetric radar data near a triple junction of HIR in Dronning Maud Land, East Antarctica. We validate the inferred ice-fabric types with ice-core data near the summit and provide additional context in terms of variability in vertical strain rates and corresponding signatures in the radar stratigraphy.

2. Study area & data

HIR is a promontory ice rise located in eastern Dronning Maud Land (Fig. 1). It has a discernible dome at 367 meters above sea level (m a.s.l.) (Howat and others, 2022) that is co-located with a triple junction from which three ridge divides extend into the ice-rise flanks. The ice thickness at the dome is approximately 550 m (Fig. 7). The average accumulation rate and mean ice thickness within the 5 km pRES profile are reported as 0.4 m a⁻¹ (Cavitte and others, 2022) and 550 m, respectively. The ratio of both values (thickness/accumulation) provides a characteristic time scale (t _D), which is a reference of the time it takes for a change to advert through the system (Martín and others, 2009a). For HIR, t _D is approximately 1400 years. In this study we use three different dataset collected at HIR.

Phase coherent radar data: In 2019, 15 static, quad-polarimetric measurements were taken along a 5 km profile crossing the triple junction HIR in northwest-to-southeast direction (Fig. 1b—red line). At each site, we infer the magnitude and the orientation of ice fabric with depth (sec. 3.1). One static measurement (site name p0) was taken at a few tens of meters distance from the ice core, which validates our inference with values derived from ice-core data (sec. 3.2). In 2020, all static sites were revisited to determine the yearly-averaged vertical strain rates (sec. 3.3).

Airborne radar data: The airborne radar data were collected in December 2018 and January 2019 as part of CHIRP (Channel and Ice Rise Project; Jansen and others, in: Fromm and others, 2019) using the ultra-wideband radar system (UWB) of the Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung (2016) with a frequency range of 150–520 MHz. The system was deployed to survey the area providing ice thickness and internal ice stratigraphy data (sec. 3.4) at multiple cross sections roughly oriented along the East-West direction (Fig. 1b—black dashed lines).

Ice core fabric data: During the 2018–2019 austral summer field season, a 263 m long ice core was drilled at the summit of HIR (70.49960$^\circ$S, 21.88019$^\circ$E) (Fig. 1b—black dot). The ice core provided the age-depth relationship used to date near-surface radar stratigraphy imaged with a different ground-based radar in order to extrapolate the surface mass balance spatially (Cavitte and others, 2022) (Fig . 2b). The ice core was also analyzed to investigate ice crystal fabric. In this study we will only use the ice core fabric data to verify inferences drawn from 15 quad-polarimetric radar observations.

Figure 2. A cross-sectional view along the extended pp’ profile (Fig. 1), illustrating (a) surface velocity (Shallow-ice approximation and Rignot and others (2017)), and (b) surface mass balance (Lenaerts and others, 2014; Cavitte and others, 2022). The red shading corresponds to the location of the pRES profile.

This publication marks the first release of all the data presented here, with the exception of the AA’ UWB profile (illustrated by the black dashed lines in Fig. 1b), which was previously published by Koch and others (2023a). Additionally, an approximation of surface velocities and corresponding horizontal strain rates based on the shallow ice approximation (SIA) (sec. 3.5) is provided.

3. Methods

3.1. Ice-fabric derived from static, phase-coherent radar

Propagation of radio waves through ice is polarization dependent because ice is mechanically and dielectrically anisotropic (Hargreaves, 1977, 1978; Fujita and others, 2006). More specifically, radio wave speed depends on the orientation of the ice crystals relative to the radio-wave polarization which leads to variability in backscattered power (through birefringence and anisotropic reflections) as a function of antenna orientation at the surface. The degree and type of anisotropy in ice, in short the ice-fabric type, is often described using three eigenvectors ($\vec {v}_1$, $\vec {v}_2$, $\vec {v}_3$) and eigenvalues (λ ₁, λ ₂, λ ₃ with λ ₁ < λ ₂ < λ ₃ and λ ₁ + λ ₂ + λ ₃ = 1) which correspond to an ellipsoid best describing the bulk orientation of individual crystal c-axis at a given depth. The directions are locally defined, but can be georeferenced using the antenna orientation at the surface. Inferring the anisotropic ice properties from polarimetric radar data has been the subject of many previous studies (Dall, 2010; Jordan and others, 2019, 2020; Ershadi and others, 2022; Zeising and others, 2023) and some consensus has emerged that the polarization dependence can be fully captured using a quad-polarimetric setup in which four antennas are oriented perpendicuarily to each other (Young and others, 2021; Ershadi and others, 2022, 2024). Following the notation from the satellite remote sensing literature, the two orthogonal polarizations are referred to as horizontal (H) and vertical (V), although they are both situated in the horizontal plane. Each quad-polarimetric measurement consists of four individual measurements with co-polarized (HH, VV) and cross-polarized (HV, VH) orientations. The data can be synthesized to mimic a full azimuthal orientation of the antennas, and variations in backscatter power are displayed correspondingly (Young and others, 2021; Ershadi and others, 2022, 2024).

Here, the quad-polarimetric data (Pattyn and others, 2023) at each site were collected with a fixed antenna distance (5 m between centers), and the absolute, georeferenced orientation of the baseline connecting the two antennas is determined with a compass with approximately 15$^\circ$ uncertainty. We determine the horizontal ice fabric anisotropy (Δλ _H = λ ₂ − λ ₁) and its georeferenced orientation as the direction of the strongest horizontal eigenvector ($\vec {v}_2$) using a polarimetric forward model (Fujita and others, 2006) and an and inversion outlined in Ershadi and others (2022). This method employs HH and HV power anomaly data and the HHVV coherence phase, defined as the argument of the complex polarimetric coherence and its scaled phase derivative, which estimates the depth variability of Δλ _H and $\vec {v}_2$ assuming that one (in this case $\vec {v}_3$) of the eigenvectors is pointing vertically. Additionally, the method allows for the estimation of all three eigenvalues assuming that ice is isotropic at the surface. This enables the reconstruction of the vertical anisotropy (Δλ _V = λ ₃ − λ ₂) in a top-to-bottom approach. In this case a weak Δλ _H would be reflected in a smoothly varying coherence phase. A strong Δλ _H, on the other hand, would result in multiple nodes where the coherence phase is wrapped at the 2π boundaries. For HIR specifically, we limit our analysis to a magnitude of coherence of 0.4 following recommendations from Jordan and others (2019). This covers approximately the upper 400 m, corresponding to approximately 70% of the total ice thickness near the dome (Fig. 7).

To categorize the various observed ice fabric types and their development at different depths, we employ a classification method that uses the logarithmic ratios of the eigenvalues. This approach effectively distinguishes between cluster-type (point maximum) and girdle-type fabrics, as outlined by Woodcock (1977). The key parameters in this scheme are K = ln(λ ₃/λ ₂)/ln(λ ₂/λ ₁) and C = ln(λ ₃/λ ₁), where K serves to identify whether the fabric is a uniaxial girdle or cluster, and C measures the intensity of the identified ice fabric type. K and C are later referred to as ‘Woodcock parameters’. The evolution of fabric types in relation to flow regimes is well described by Llorens and others (2022) providing comprehensive models and visual representations that elucidate the relationship between ice deformation and the resultant fabric patterns.

3.2. Ice-fabric from ice-core data

The ice core was cut in 0.5 m sections on site, then packed, transported to and stored at the Laboratoire de Glaciologie (Université libre de Bruxelles (ULB), Belgium) respecting the cold chain (temperature below −25°C) at all times. Dating and interpretation of a series of environmental and climatic proxies for the upper 120 meters of the core are beyond the scope of this paper and are presented separately in Wauthy and others (2024). Here we will focus on the ice-fabric properties of the entire ice core, more specifically the eigenvalues of the eigenvectors, characterizing the ice-fabric anisotropy that we aim to reconstruct from the pRES measurements.

To determine the eigenvalues of the ice fabric from the ice core, 114 regularly spaced 8 cm high and 500 μm thick vertical thin sections of ice were produced following the standard procedure of Langway (1958). The thickness of the ice core sections, typically between 500 μm and 600 μm, ensures that there is no superposition of crystals, allowing for accurate 3D fabric analysis. The Automatic Fabric Analyzer effectively measures the orientation of individual pixels and uses image analysis to determine grain boundaries and calculate the mean orientation within each grain, providing robust data for deriving eigenvectors. Crystal (optic) c-axes orientations were measured using the G-50 Automated Ice Fabric Analyzer (Russell-Head Instruments, e.g. Wilson and others, 2003). Eigenvectors and eigenvalues were calculated using the FAME software (Hammes and Peternell, 2016). The same software was used to determine grain boundaries, to plot c-axis orientation density distributions in a lower hemisphere, equal-area or Schmidt diagram. Schmidt diagrams are a common representation in geology providing equi-areal 2D projections of the ice crystal's c-axes intersection with a lower hemisphere into the equatorial plane, chosen in the plane of the vertical thin sections in this study. Density diagrams are constructed by counting the number of c-axes falling in a reference counting circle displaced on a regular grid across the Schmidt diagram.

3.3. Vertical strain rate

The sites used for the polarimetric surveys (sect. 3.1) were marked with bamboo stakes and revisited one year later. The phase-coherent repeat measurements enable tracking of the submergence of internal reflections relative to the bed (Kingslake and others, 2014). This allows us to infer yearly averaged vertical strain rates, a method which is commonly applied to ice shelves in order to isolate the basal melt rate signal from observed thickness change (e.g. Nicholls and others, 2015; Sun and others, 2019). For HIR specifically, we calculated depth-averaged values of vertical strain rate for ice thickness intervals over tens of meters in order to highlight signatures of the Raymond effect.

3.4. Airborne radar data

The UWB radar is an improved version of the Multichannel Coherent Radar Depth Sounder (MCoRDS 5) developed at the University of Kansas, Center for Remote Sensing and Integrated Systems (Rodriguez-Morales and others, 2014; CReSIS, 2021), operated on AWI's Polar6 BT-67 aircraft (Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, 2016). The radar system consists of an eight-element antenna array polarized in HH, which serves as a transmitter and receiver unit for radar signals. Data acquisition and processing methods are detailed in Koch and others (2023a) and are similar to those described by Franke and others (2021) and Franke and others (2022). During CHIRP, the radar transmitted three-stage linear modulated chirp signals (1 μs low-gain, 1 μs high-gain and 3 μs high-gain to sound the upper, middle and deeper part of the ice column in high quality) in a frequency range of 150–520 MHz and at an acquisition height of ~ 360 m above the ice surface. Radar data processing was conducted with the CReSIS Toolbox (CReSIS, 2021) and comprises pulse compression, synthetic aperture radar (SAR) processing with a wide angular range, and array processing (Rodriguez-Morales and others, 2014; Hale and others, 2016; Franke and others, 2022). The processed radar data have a range resolution of ~ $0.35$ m and an along-track trace spacing of approximately 6 m. Here, we use selected sections of the airborne radar data to analyze signatures of the Raymond arches beneath the dome and the landward-oriented ice divide (Fig. 1).

3.5. Shallow ice approximation: surface velocities and strain rates

Surface velocities at HIR are too low to be reliably measured up by remote sensing data (Fig. 2a). Therefore, we use the shallow-ice approximation (SIA; Hutter, 1983; Greve and Blatter, 2009) as a rough estimate of the surface velocity and maximum horizontal strain rate (${\dot {\bf \varepsilon }}_{max}$), whilst being aware that a higher-order ice flow model would be more accurate in the region. We use the calculated surface flow direction and the maximum strain rate direction, ${\dot {\bf \varepsilon }}_{max}$, to compare with the estimated strongest horizontal anisotropy eigenvector, $\vec {v}_2$. The map of HIR with the estimated magnitude and orientation of the surface velocity and maximum horizontal strain rate is shown in Appendix D.

Our calculation of velocities using SIA is not without uncertainty. Although bed elevation errors in BedMachine data are relatively low at Hammarryggen Ice Rise, there are some error estimates of up to 100 m on the southern side of the ice rise away from the radar profile (Morlighem and others, 2020). Furthermore, we have made the assumption that ice is isothermal, but given that we are most interested in comparing strain rate directions with the observational anisotropy data rather than strain rate magnitudes, errors due to this assumption are likely to be small.

4. Results

4.1. Inference of ice-fabric parameters from pRES measurements

We use the pRES measurement site closest to the ice-core site (marked p0 in Fig. 1) to illustrate results from the quad-polarimetric analysis. The observations from the quad-polarimetric measurements are displayed using multiple metrics. The HH power anomaly (Fig. 3a) represents the backscatter dependence as a function of antenna orientation and is indicative of anisotropic reflections, e.g. due to vertical variability in ice-fabric strength. The HHVV coherence phase (Fig. 3b) shows the phase correlation between the HH and VV directions. Stronger vertical gradients correspond to a stronger Δλ _H. The HV power anomaly (Fig. 3c) is an analog to the HH power anomaly but for the depolarization component and is a proxy for the ice-fabric orientation (marked with green dots). The scaled phase derivative (Fig. 3d) of the ice-fabric orientation for a given depth interval (marked with green dots) is defined as Δλ _H. Figure 3e–f show the same metrics based on a radio-wave propagation model (Fujita and others, 2006) and ice-fabric parameters resulting from a non-linear optimization method (Ershadi and others, 2022).

Figure 3. Results for the p0 radar site: (a) to (d) pRES observations, with green dots in (c) and (d) marking the minima in PHV. (e) to (h) Optimized model output capturing the principal patterns of the observations.

The characteristic signatures (e.g. nodes, location of maxima, etc.) in the observations (Figs. 3a–d) are well reproduced by the optimized forward model output (Figs. 3e–h) demonstrating that the inferred ice-fabric eigenvalues and their changes with depth are adequately captured by the inversion. The gradient in the polarimetric phase coherence indicates a gradual strengthening of Δλ _H with depth (Fig. 3d), and the minima in the HV power anomaly suggest that the ice-fabric orientation changes are small with depth (Fig. 3b). An exception occurs in the depth interval between 150 and 200 m, where a cross-polarization extinction node suggests a rotation of the $\vec {v}_2$ eigenvector of several degrees (Figs. 3c,g). We first substantiate the inferred ice-fabric parameters from the radar polarimetry by comparing them to ice-core measurements in the following section, and then continue by tracing the ice-fabric parameters away from the ice core into the ice-rise flanks.

4.2. Ice core validation

The fabric data measured from ice core samples show an increase with depth of λ ₃ and a decrease of both λ ₁ and λ ₂ (Fig. 4). The measured Δλ _H indicates a weak horizontal anisotropy within the ice column and remains almost constant with depth. In contrast the measured Δλ _V increases with depth. This behavior of eigenvalues results in Woodcock parameters K > 1 and C < 2, categorizing the fabric type into a weak uniaxial cluster. This pattern (increasing areal concentration of crystals’ c-axes from white to red) is evident in the density Schmidt diagrams (Fig. 4c), which are directly measured from the ice core. Additionally, Fig. 9 allows for a better comparison between the observed and estimated fabric types, demonstrating that the fabric is nearly isotropic and evolves toward a weak uniaxial cluster.

Figure 4. Comparison between estimated and measured (a) eigenvalues, (b) horizontal and vertical ice fabric anisotropy as Δλ _H and Δλ _V, respectively and (c) Woodcock values K and C with density Schmidt diagrams measured from the ice core. Note that the estimated values are the results from the inverted radar data, and the measured values are from the ice-core laboratory analysis.

The estimated eigenvalues from the quad-polarimetric radar measurement at site p0 are compared with the measured ice-core eigenvalues (Fig. 4). The estimated eigenvalues and anisotropy in both the horizontal and vertical directions exhibit the same behavior as the measured ones. However the estimated λ ₁ and λ ₂ (Fig. 4a) are about 0.07 and 0.03 larger than the measured values, respectively, and consequently, the estimated λ ₃ is systematically smaller than the measured value. Both estimated and measured Δλ _H are weak (approximately 5 % of the maximum possible horizontal anisotropy Δλ _H = 1 (Fig. 4b), with the estimated one being slightly weaker than the measured one). In contrast, both the estimated and measured vertical ice fabric anisotropy Δλ _V increase with depth (Fig. 4c). Similar to Δλ _H, the estimated Δλ _V is also weaker than the measured Δλ _V.

Similarly to the ice core data, the radar-derived fabric shows a tendency to form clusters which increase in strength with increasing depth (Fig. 4c). The differences seen in the eigenvalue magnitudes correspondingly translate into the K and C classification: The estimated C values (color of marks in Fig. 4c) are weaker than the measured ones, particularly on the shallower part of the ice column. The minimum C value estimated from radar at site p0 is 0.19, and the maximum is 1.81. In contrast, the ice core values are 0.36 and 2.35, respectively. The estimated Δλ _H between 350 and 380 increases to 0.12 (Fig. 4b), resulting from the corresponding change in λ ₁ and λ ₂ (Fig. 4a). This behavior does not affect Δλ _V (Fig. 4c), but it does affect the K value (Fig. 4c) which is close to unity. However, no ice-core data are available at that depth to validate this behavior. It is important to note that Fig. 4c shows a limited range of fabric types, while Fig. 9 in the appendix provides a fuller context for better comparison between measured and estimated fabric. Although the fabric type is broadly captured, its depth variability is overestimated by the pRES data. This overestimation stems from systematically low horizontal anisotropy values (Fig. 4b—blue line), which are disproportionally amplified because the low horizontal anisotropy values appear in the denominator of K.

4.3. Spatial changes in ice-fabric and vertical strain rates along the 5 km transect

After comparing the consistency between the estimated eigenvalues derived from polarimetric radar data at the p0 site and the measured ice core eigenvalues, we reconstruct ice-fabrics for all sites p1 to p14 along the 5 km long transect. To interpret our results, we normalize distances and elevation with the ice thickness at the dome (H ≃ 550 m). The distance of the pRES points from the dome denoted as X is normalized as x = X/H. Additionally, elevation is expressed as the normalized ice height above the bed, denoted as z = (H − Z)/H, where Z represents the depth. In this context, the mean bed elevation and mean surface elevation along the pRES profile are designated as z = 0 and z = 1, respectively. Subsequently, we employ linear interpolation to obtain the spatial variation of the fabric parameters along the 2D transect (Fig. 8).

Depth-averaged values of the horizontal anisotropy Δλ _H show differences on both sides of the divide (Fig. 8a). On the south-eastern side, where ice is thicker, values of Δλ _H are in general larger and more variable than on the north-western side. In the 30–35% depth-interval, the averaged Δλ _H exhibits a local maximum beneath the summit that is approximately one ice thickness wide and is asymmetrical. The north-western side also exhibits slightly smaller maxima beneath the ice-rise flanks. The spatial distribution of the magnitude of the strongest estimated eigenvector λ ₃ (Fig. 8b) exhibits a similar pattern in terms of a local maximum beneath the divide and has generally larger values on the north-western side. The depth-average orientation of $\vec {v}_2$, aligns within 10° with the North-South direction (Fig. 8c). This direction is ~40° offset to the mean flow direction in the ice-rise flanks and ~81° offset to the direction of maximum horizontal strain inferred from the SIA-based velocity field. The magnitude of the depth-averaged vertical strain rates (Fig. 8d) is highest in the top 20° of the ice thickness (80 to 100° depth interval), where the densification of firn is strongest. Vertical strain rates are also overall smaller in absolute value in the thinner north-western flank than the thicker south-eastern flank. At approximately 50° of the ice thickness, the vertical strain rates exhibit a pronounced (weakly double-peaked) minima beneath the divide which extends laterally for 1–2 ice thickness into the ice-rise flanks.

4.4. Internal stratigraphy

The airborne UWB radar profiles (Fig. 6) image ice thickness and internal radar stratigraphy in profiles located nearly perpendicular to the local ice divides (Fig. 1). The average ice thickness is between 500 and 600 m beneath the divides. The bed increases in elevation toward the west and deepens from the triple junction into the landward direction. The bed beneath the saddle (Profile B-B’) appears distinctly rougher than beneath the dome area (Profile AA’). The internal radar stratigraphy is clearly visible in both profiles but cannot be identified unambiguously at depths deeper than the surface multiple (Koch and others, 2023a). Continuous tracking of the stratigraphy is also difficult in areas where internal layers are more inclined (i.e. near the divides) (Holschuh and others, 2014), and in areas where the flight track is curved (Fig. 1). Nevertheless, internal radar stratigraphy close to the surface appear deeper in the south-eastern flanks compared to the north-western flanks, and their syncline arching beneath the divide is clearly visible in B-B’ (i.e. beneath the saddle) and to a lesser extent also along AA’ (just north-west of the dome). The arches increase in amplitude with increasing depth and are vertically aligned with today's divide position (Fig. 1).

5. Discussion

Previous studies have investigated ice-rise evolution using flow-line modeling in combination with the internal isochronal radar stratigraphy as principal observations (Drews and others, 2013, 2015; Goel and others, 2017, 2018; Martín and others, 2009a, 2009b; Martín and Gudmundsson, 2012; Hindmarsh and others, 2011). Two additional studies of a dome and ice rise, respectively, used the observed vertical strain rates (Gillet-Chaulet and others, 2011; Kingslake and others, 2014). Here we use all of the previous observations and add quad-polarimetric radar measurement as another possible observational constraint. We now investigate whether those observations capture signatures of the Raymond effect and, if so, how these can be contextualized with other geophysical observations of the contemporary flow regime. This may guide the application of a future 3D model (incl. thermo-mechanical coupling and anisotropic rheology) which is capable of simulating the complex dynamics occurring at triple junction ice rises. Given that the extraction of ice-fabric parameters from quad-polarimetric data using non-linear inversion has so far only once been compared with direct ice-core measurements (Ershadi and others, 2022), we first discuss the benefits and limitations of this method in general before moving on to investigate the flow history of HIR.

5.1. Applicability of the inferred the ice-fabric eigenvalues

The quad-polarimetric analysis has a limitation in that it assumes one of the principal eigenvectors points upwards. Although this assumption can be relaxed (Rathmann and others, 2022), it leads to a more complicated forward model for which the inversion is not yet established. However, Rathmann and others, 2022 show that the polarimetric radar response of nadir-looking radars is comparatively insensitive to ice fabrics that are vertically tilted. However, beneath ice domes vertical compression is assumed to dominate, which is expected to lead to a vertical point maximum in the c-axes distribution (Budd and Jacka, 1989; Llorens and others, 2022). We, therefore, consider the assumption of horizontal and vertical eigenvectors to be justified, and not likely a cause for the systematic mismatch in magnitude that we observe between the eigenvalues from the quad-polarimetric method and the ice-core-based values (Fig. 4a).

The systematic underestimation of Δλ _H and Δλ _V compared to ice-core values has to a lesser extent also been observed at the ice-core site of the European Project for Ice Coring in Antarctica (EPICA) at Concordia Dome C (EDC), however, it does not occur at the EPICA site in Dronning Maud Land (EDML) (Ershadi and others, 2022). We investigated if using a scaling factor to the dielectric anisotropy for a single crystal (commonly assumed to be 0.034 Matsuoka and others, 1997) can explain the underestimation. However, the mismatch did not significantly improve when changing dielectric anisotropy within the reported uncertainties. The inversion is also sensitive to the fabric orientation and backscatter ratio. The latter in turn varies according to the ice-core data on shorter spatial scales than what the inversion can currently resolve, particularly because it involves vertical averaging to smooth the phase gradient. The reason for the underestimation of Δλ _H and Δλ _V therefore requires further investigation, but given that the gradients are well reproduced, this does not hinder the interpretation of lateral ice-fabric variability.

5.2. PRES detects geo-referenced fabric orientation

The estimated $\vec {v}_2$, as depicted in Fig. 5c, is derived solely from pRES data, without validation from field datasets. To overcome this limitation, we used surface flow direction data obtained from SIA modeling to compute the eigenvectors and eigenvalues of the strain rate tensor ${\dot {\bf \varepsilon }}$. When comparing $\vec {v}_2$ to the surface flow direction, a deviation of ~40° is observed (Fig. 5c blue vs. dashed red). In contrast, when compared to the direction of maximum horizontal strain rate, $\vec {v}_2$ shows a deviation of ~81° (Fig. 5c blue vs. solid red).

Figure 5. (a) Depth-averaged variation of Δλ _H within a specific depth window. (b) Depth-averaged variation of λ ₃ within a specific depth window. (c) Depth-averaged horizontal ice fabric orientation (blue line), surface flow direction derived from SIA (dashed red), and maximum strain direction derived from SIA (red line). (d) Vertical strain rates measured at each pRES site averaged over different depth intervals. Note that more negative strain rates indicate stronger deformation. The x-axis is the distance from the dome normalized by H.

It was established by Alley (1992) that during ice deformation, c-axes consistently rotate toward compressional axes and away from tensional axes. Also, the principles of fabric orientation under vertical shortening is discussed by Passchier (1997) where the theory explains that basal planes rotate toward the horizontal plane, which serves as the fabric attractor. Consequently, the perpendicular c-axes rotate toward the vertical direction. The rotation is most rapid in the plane containing the direction of maximum shortening (vertical) and maximum stretching. As a result, the variation in the horizontal c-axes, described by λ ₁ in the direction $\vec {v}_1$, is narrowest in this plane. $\vec {v}_2$ is perpendicular to this direction in the horizontal plane, hence it is expected to be oriented at 90° to the direction of maximum stretching, which indeed corresponds to our observations in Fig. 5c (~81°). Also, as suggested by the pRES measurements, λ ₁ and λ ₂ exhibit similar intensities (weak Δλ _H), it follows that the same might hold true for $\dot {\varepsilon }_1$ and $\dot {\varepsilon }_2$. Their combination would then yield maximum horizontal strain at approximately 45° from $\vec {v}_1$ and $\vec {v}_2$. This explains why $\vec {v}_2$ appears at approximately 45° from the flow direction in Fig. 5c (~40°).

5.3. Synthesis of radar observations within the ice-dynamic setting of HIR

The radar stratigraphy, the strain rates and the ice-fabric properties are all jointly influenced by the ice-dynamic evolution of HIR and encode parts of its history, even though it is not yet clear how rapidly ice-fabric parameters change with the ice-dynamic flow regime. Here we synthesize the different datasets with a particular focus on the Raymond effect and contextualize our findings with available modeling and observational studies of other ice rises.

The upward arches observed beneath the saddle (Fig. 6, BB’ flight line) are typical of ice rises in the sense that they are located beneath today's divides and that they are asymmetrical in shape. For example, a syncline as on the western side has also been observed at Derwael Ice Rise and explained with persistent accumulation patterns including erosion of snow at the crest and re-deposition in the flanks (Drews and others, 2015). Erosion of snow at the crest increases the amplitudes of the upward arches at larger depths which are, however, primarily formed by the Raymond effect. Both mechanisms require a stable ridge divide position and therefore testify that the saddle connecting HIR with the main ice sheet was stationary, probably for several t _D, i.e. several thousands of years. Upward arching also occurs beneath the dome (Fig. 6, at kilometer 10 in profile AA’), but the amplitudes are smaller compared to the saddle. The eastern side of AA’ is near-parallel to the eastern arm of the triple junction and hence strong upward arching is not expected in the stratigraphy here. It is therefore unclear if the triple junction of HIR exhibits a Raymond cupola as modeling would suggest (Hindmarsh and others, 2011), but if it does, the lower arch amplitudes could suggest that the dome position is younger than the saddle, although three-dimensional effects may be responsible.

Figure 6. Airborne UWB radargrams crossing two ridges of the triple junction dome (AA’) and the saddle ridge (BB’). Red curves highlight laterally coherent internal reflection horizons, and red dashed lines contain in parts data gaps, particularly in areas where the layers are more inclined.

For a two-dimensional, plane strain flow regime, it is well understood that lateral differences in vertical velocities that accompany the formation of Raymond arches, are expressed by corresponding patterns in the vertical strain rates. More specifically, the vertical strain rates are expected to be smaller in magnitude for 100 to 300 m depths (z from ~0.8 to 0.5) beneath the divide compared to the flanks (Kingslake and others, 2014). Our observations (Fig. 5d) comply with these predictions, particularly for 330 to 380 m depths (z from ~0.4 to 0.3). At shallow depths (top 100 m) the vertical strain rates are dominated by firn compaction, and deeper depth intervals could not be resolved. Observed magnitudes of approximately 1.0 × 10⁻³ m a⁻¹ are comparable to what has been observed at other triple junctions (i.e. Fletcher Promontory Kingslake and others, 2014), although the amplitude of the vertical strain rate anomaly across the dome is smaller. Once a local divide or dome has formed, the effect on the velocity field is instantaneous, and hence the vertical strain rates do not contain information about the ice-rise history per se. However, ice-fabric types are strain-induced and develop over time (Budd and Jacka, 1989). Consequently, if the dome position was temporally stable, corresponding signatures should appear in the derived ice-fabric types from the quad-polarimetric analysis, indicating a temporarily stable dynamic regime.

Regarding the ice fabric, below 150 m where the distribution of the orientation of the fabric slowly (low values on the density scale) evolves from a more random distribution in the top of the ice core toward a single maximum closely centered on the vertical as expected from dominant uniaxial compression at ice domes (e.g. Durand and others, 2007). The gradual strengthening of the fabric anisotropy is clearly seen in the evolution of the measured eigenvalues (Fig. 4a). The small but increasing horizontal anisotropy (Fig. 4b), indicating that the strain is not purely uniaxial flattening (compaction) but includes differential deformation (such as lateral extension) in the horizontal plane which is coherent with the complexity of the geomorphological setting (triple junction). The ice-fabric reconstruction from the quad-polarimetric data shows that minima in the vertical strain rates (Fig. 5d) are accompanied by corresponding maxima in Δλ _H and λ ₃ (Figs. 5a,b) in 330 to 380 m depth interval (z from ~0.4 to 0.3). This is in line with measured ice fabric and two-dimensional model predictions of Martín and others (2009a) which predicts a single maximum fabric which is stronger beneath the divide compared to the flanks.

A quantitative comparison in terms of timing between our observations and the model predictions of Martín and others (2009a) is hampered in several ways: first, the assumed two-dimensional geometry does not include the triple junction geometry of HIR, and second, the model predictions assume an evolution from fully isotropic to fully anisotropic ice. The latter is unlikely to be the case for HIR as demonstrated by the measured ice fabric data. Notwithstanding, in steady-state (i.e. at approximately 10 times t _D) the predicted degree in ice-fabric anisotropy is larger than what is reconstructed from quad-polarimetric data here. The reconstructed Δλ _H consistently remains below 0.1 which is comparable to other domes such as Dome C, but is much weaker than what has been observed in flank flow regimes such as the transient divide at the EDML drill site (Δλ _H > 0.3, Ershadi and others, 2022). Based on these comparisons, it appears that HIR in terms of its ice-fabric characteristic is not older than 4 times t _D (i.e. not older than approximately 5600 years). However, given the discrepancies between the model assumptions and observations, this time interval is not well constrained.

Taken together, the UWB radar profile across the saddle suggests a temporally stable divide position. The data at the dome are less conclusive in that sense, because arch amplitudes are smaller and because the ice fabric is only weakly developed. One plausible scenario uniting this would be that HIR undergoes a transition from a promontory toward an isle-type ice rise, which is a feature of deglaciation scenarios in this particular region (Favier and Pattyn, 2015). Thinning in the saddle area would then result in comparatively large arches relative to today's ice thickness in this area. The good match to the ice-core data reinforces that quad-polarimetric surveys can be a reliable tool to further constrain ice-rise evolution, in particular the influence of ice-anisotropy on Raymond arch evolution. For HIR, the comparatively weak ice -fabric suggests a comparatively young dome. However, a single two-dimensional profile heavily simplifies the dynamic complexity and modeling should account for these three-dimensional effects in the future.

6. Conclusion

We have investigated radar-derived properties of Hammarryggen Ice Rise (HIR): radar stratigraphy, strain rates, and ice-fabrics. HIR is a representative triple junction promontory ice-rise, making it an excellent laboratory to study ice dynamic processes, where we additionally, had access to both the ice core for c-axes measurements and the corresponding radar data.

Upward arching in the stratigraphy indicates a stable ice divide in the saddle area over, at least, several thousands of years. Upward arching beneath the dome is also observed but is less clear. Vertical strain rates are dominated by firn compaction near the surface, and exhibit a minimum closer to the bed indicative for the Raymond effect. The derived ice-fabric properties from quad-polarimetric radar fit ice-core-based values. The horizontal anisotropy is weak and thus young compared to steady-state, ice-dynamically evolved ice-fabric types predicted from two-dimensional models in comparable settings. This is perhaps indicative of thinning of the saddle connecting the dome to the mainland. There are also signatures of the Raymond effect in the ice-fabric. However, it is unclear how the triple junction geometry of Hammarryggen Ice Rise impacts both the vertical strain rates and the ice-fabric development. Previous studies have indicated that the region is ice-dynamically stable and comparatively resilient to sea-level changes (Drews and others, 2015; Favier and Pattyn, 2015). Our study on Hammarryggen Ice Rise provides further evidence for this stability, although it is the first instance where we suspect the dome position may have a younger history compared to the connected saddle. This could be an important consideration when using ice rises as proxies for ice-dynamic changes in their respective catchments.

Overall, the synthesis of the different radar observations has the potential to constrain unknown parameters like the ice fabric in future ice-flow modeling, particularly if measurements cover larger areas. We suggest that these additional geophysical constraints provide another step forward toward a quantitative interpretation of Raymond arch amplitudes using observationally constrained, anisotropic, three-dimensional ice-flow models of triple junctions, flow regimes common to many ice rises around Antarctica. To better understand the ice fabric and the dynamics of a triple junction ice rise, it is advised that future pRES measurement campaigns have profiles perpendicular to each ridge.