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, Reference Matsuoka2015). 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, Reference Favier, Gagliardini, Durand and Zwinger2012, Reference Favier2014; Favier and Pattyn, Reference Favier and Pattyn2015; Schannwell and others, Reference Schannwell2019; Henry and others, Reference Henry, Drews, Schannwell and Višnjević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, Reference Raymond1983). 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, Reference Davies, Hambrey, Smellie, Carrivick and Glasser2012), are unavailable for most of the Antarctic perimeter.
(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, Reference Howat2022). 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, Reference Howat2022), bed elevation and ice thickness (Morlighem, Reference Morlighem2022).
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, Reference Matsuoka2015). Clear signatures of transience are Raymond stacks that do not align with contemporary ice divides (Nereson and Waddington, Reference Nereson and Waddington2002), such as at Siple Dome (Nereson and others, Reference Nereson, Hindmarsh and Raymond1998). 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, Reference Drews2015). Cases between these two end members (Goel and others, Reference Goel2020) are more difficult to interpret and require advanced model-data comparison, including thermomechanically-coupled full Stokes models with anisotropic rheology (Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009a, Reference Martín, Hindmarsh and Navarro2009b; Martín and Gudmundsson, Reference Martín and Gudmundsson2012) and a dynamically evolving grounding line (Schannwell and others, Reference Schannwell2019, Reference Schannwell2020; Henry and others, Reference Henry, Drews, Schannwell and Višnjević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, Reference Martín and Gudmundsson2012; Drews and others, Reference Drews2015) 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, Reference Gillet-Chaulet, Hindmarsh, Corr, King and Jenkins2011; Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009a, Reference Martín, Hindmarsh and Navarro2009b; Drews and others, Reference Drews2015; Bons and others, Reference Bons2018). In contrast, the along-ridge flow component generally produces smaller arches (Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009a, Reference Martín, Hindmarsh and Navarro2009b). Similarly, variations in bed topography can lead to smaller arches when the bed is uneven (Kingslake and others, Reference Kingslake2014), while basal sliding also contributes to reduced arch sizes (Petit and others, Reference Petit2003). Additionally, localized factors such as surface mass balance and erosion at the crest can increase arch amplitudes (Conway and Wilbour, Reference Conway and Wilbour1999; Drews and others, Reference Drews2015). The historical thinning or thickening of the ice further translates to changes in arch sizes relative to their current geometry (Martín and others, Reference Martín, Hindmarsh and Navarro2006; Goel and others, Reference Goel, Martín and Matsuoka2018).
Ice-core analysis, in combination with shallow and deep radar, can constrain the three-dimensional ice geometry (Hindmarsh and others, Reference Hindmarsh2011) and the surface accumulation history (Philippe and others, Reference Philippe2016; Goel and others, Reference Goel, Brown and Matsuoka2017; Cavitte and others, Reference Cavitte2022). Strain measurements such as the coffee-can method (Hamilton and Whillans, Reference Hamilton and Whillans2000) and repeat surveys with phase-coherent radar can provide additional constraints on the vertical strain rates (Kingslake and others, Reference Kingslake2014). 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, Reference Drews2015). Consequently, so far, ice rises and their inferred dynamic history play a minor role in constraining larger-scale ice flow models (Bentley and others, Reference Bentley2014).
Phase coherent radar polarimetry using a ground-based phase-sensitive Radio Echo Sounder (pRES) (Brennan and others, Reference Brennan, Lok, Nicholls and Corr2014) 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, Reference Dall2010; Jordan and others, Reference Jordan, Schroeder, Castelletti, Li and Dall2019, Reference Jordan, Schroeder, Elsworth and Siegfried2020; Ershadi and others, Reference Ershadi2022; Rathmann and others, Reference Rathmann2022; Zeising and others, Reference Zeising2023). Anisotropic ice-flow models of steady-state ice rises, as detailed in studies by Martín and others (Reference Martín, Gudmundsson, Pritchard and Gagliardini2009a, Reference Martín, Hindmarsh and Navarro2009b); Martín and Gudmundsson (Reference Martín and Gudmundsson2012), 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, Reference Howat2022) 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−1 (Cavitte and others, Reference Cavitte2022) 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, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009a). 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, Reference Fromm, Oberdieck, Heitland and Köhler2019) 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, Reference Cavitte2022) (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.
A cross-sectional view along the extended pp’ profile (Fig. 1), illustrating (a) surface velocity (Shallow-ice approximation and Rignot and others (Reference Rignot, Mouginot and Scheuchl2017)), and (b) surface mass balance (Lenaerts and others, Reference Lenaerts2014; Cavitte and others, Reference Cavitte2022). 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 (Reference Koch2023a). 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, Reference Hargreaves1977, Reference Hargreaves1978; Fujita and others, Reference Fujita, Maeno and Matsuoka2006). 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 (λ 1, λ 2, λ 3 with λ 1 < λ 2 < λ 3 and λ 1 + λ 2 + λ 3 = 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, Reference Dall2010; Jordan and others, Reference Jordan, Schroeder, Castelletti, Li and Dall2019, Reference Jordan, Schroeder, Elsworth and Siegfried2020; Ershadi and others, Reference Ershadi2022; Zeising and others, Reference Zeising2023) 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, Reference Young2021; Ershadi and others, Reference Ershadi2022, Reference Ershadi2024). 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, Reference Young2021; Ershadi and others, Reference Ershadi2022, Reference Ershadi2024).
Here, the quad-polarimetric data (Pattyn and others, Reference Pattyn, Wauthy, Sun, Tison and Tsibulskaya2023) 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 = λ 2 − λ 1) and its georeferenced orientation as the direction of the strongest horizontal eigenvector ($\vec {v}_2$
) using a polarimetric forward model (Fujita and others, Reference Fujita, Maeno and Matsuoka2006) and an and inversion outlined in Ershadi and others (Reference Ershadi2022). 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 = λ 3 − λ 2) 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 (Reference Jordan, Schroeder, Castelletti, Li and Dall2019). 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 (Reference Woodcock1977). The key parameters in this scheme are K = ln(λ 3/λ 2)/ln(λ 2/λ 1) and C = ln(λ 3/λ 1), 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 (Reference Llorens2022) 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 (Reference Wauthy2024). 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 (Reference Langway1958). 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, Reference Wilson, Russell-Head and Sim2003). 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, Reference Kingslake2014). 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, Reference Nicholls2015; Sun and others, Reference Sun2019). 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, Reference Rodriguez-Morales2014; 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 (Reference Koch2023a) and are similar to those described by Franke and others (Reference Franke2021) and Franke and others (Reference Franke2022). 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, Reference Rodriguez-Morales2014; Hale and others, Reference Hale2016; Franke and others, Reference Franke2022). 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, Reference Hutter1983; Greve and Blatter, Reference Greve and Blatter2009) 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, Reference Morlighem2020). 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, Reference Fujita, Maeno and Matsuoka2006) and ice-fabric parameters resulting from a non-linear optimization method (Ershadi and others, Reference Ershadi2022).
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 λ 3 and a decrease of both λ 1 and λ 2 (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.
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 λ 1 and λ 2 (Fig. 4a) are about 0.07 and 0.03 larger than the measured values, respectively, and consequently, the estimated λ 3 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 λ 1 and λ 2 (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 λ 3 (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, Reference Koch2023a). Continuous tracking of the stratigraphy is also difficult in areas where internal layers are more inclined (i.e. near the divides) (Holschuh and others, Reference Holschuh, Christianson and Anandakrishnan2014), 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, Reference Drews, Martín, Steinhage and Eisen2013, Reference Drews2015; Goel and others, Reference Goel, Brown and Matsuoka2017, Reference Goel, Martín and Matsuoka2018; Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009a, Reference Martín, Hindmarsh and Navarro2009b; Martín and Gudmundsson, Reference Martín and Gudmundsson2012; Hindmarsh and others, Reference Hindmarsh2011). Two additional studies of a dome and ice rise, respectively, used the observed vertical strain rates (Gillet-Chaulet and others, Reference Gillet-Chaulet, Hindmarsh, Corr, King and Jenkins2011; Kingslake and others, Reference Kingslake2014). 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, Reference Ershadi2022), 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, Reference Rathmann2022), it leads to a more complicated forward model for which the inversion is not yet established. However, Rathmann and others, Reference Rathmann2022 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, Reference Budd and Jacka1989; Llorens and others, Reference Llorens2022). 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, Reference Ershadi2022). We investigated if using a scaling factor to the dielectric anisotropy for a single crystal (commonly assumed to be 0.034 Matsuoka and others, Reference Matsuoka, Fujita, Morishima and Mae1997) 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).
(a) Depth-averaged variation of Δλ H within a specific depth window. (b) Depth-averaged variation of λ 3 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 (Reference Alley1992) 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 (Reference Passchier1997) 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 λ 1 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, λ 1 and λ 2 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, Reference Drews2015). 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, Reference Hindmarsh2011), 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.
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, Reference Kingslake2014). 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−3 m a−1 are comparable to what has been observed at other triple junctions (i.e. Fletcher Promontory Kingslake and others, Reference Kingslake2014), 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, Reference Budd and Jacka1989). 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, Reference Durand2007). 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 λ 3 (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 (Reference Martín, Gudmundsson, Pritchard and Gagliardini2009a) 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 (Reference Martín, Gudmundsson, Pritchard and Gagliardini2009a) 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, Reference Ershadi2022). 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, Reference Favier and Pattyn2015). 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, Reference Drews2015; Favier and Pattyn, Reference Favier and Pattyn2015). 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.
Data
The source code used in this study for pRES fabric analysis, strain rate analysis, and SIA is available at https://github.com/RezaErshadi/HammarryggenIceRiseSourceCode_FabricInversion_Strainrates_SIA. The pRES and ice core data can be accessed at https://zenodo.org/record/8095508, and the UWB data is available in Franke and others (Reference Franke, Jansen, Drews and Eisen2020), and Koch and others (Reference Koch2023b).
Acknowledgements
We gratefully acknowledge the data acquisition of the UWB radar data by Daniela Jansen and Steven Franke and the contributions of Steven Franke to original data processing and the Ken Borek crew and PES staff for logistic support.
Author contributions
M. Reza Ershadi led the code development and writing of the manuscript. Frank Pattyn and Sainan Sun collected the pRES data at HIR. Veronica Tsibulskaya, Jean-Louis Tison, Sarah Wauthy led the ice core analysis. M. Reza Ershadi, Reinhard Drews and Carlos Martín designed the study outline. M. Reza Ershadi analyzed the pRES data. Reinhard Drews, Inka Koch and Olaf Eisen led the UWB data analysis and interpretation. A. Clara J. Henry designed the SIA model. Falk Marius Oraschewski inferred the vertical strain rate. Paul Dirk Bons and Jean-Louis Tison led the interpretation of the fabric type. All authors contributed to the writing of the final paper.
Financial support
We would like to mention:
– M. Reza Ershadi and Reinhard Drews were supported by Deutsche Forschungsgemeinschaft (DFG) Emmy Noether grant (grant no. DR 822/3-1).
– C. Henry was supported by DFG in the framework of the priority program 1158 ‘Antarctic Research with comparative investigations in Arctic ice areas’ by a grant SCHA 2139/1-1.
– F. Oraschewsik was supported by a scholarship from the Studienstiftung des Deutschen Volkes.
– Sarah Wauthy benefited from a Research Fellow grant of the F.R.S.-F.N.R.S.
– Veronica Tsibulskaya benefited from a Research Fellow grant of the F.N.R.S.-FRIA.
– The ice-core project was supported by The Belgian Research Action through Interdisciplinary Networks (BRAIN-be) from the Belgian Science Policy Office (BELSPO) in the framework of the project “East Antarctic surface mass balance in the Anthropocene: observations and multi-scale modeling (Mass2Ant) with contract no. BR/165/A2/Mass2Ant.
– J.-L. Tison, F. Pattyn and V. Tsibulskaya and ice core collection were supported by the Belgian Science Policy Office (BELSPO - Mass2Ant project – contract n° BR/165/A2/Mass2Ant)
Competing interest
Co-author Frank Pattyn is an associate chief editor in JOG.
Appendix A. Limitations in depth of investigation
In Fig. 7, we present the signal power (blue line) and coherence magnitude (red line) at the p0 radar site (located at the center of the profile). As detailed in Sec. 3.1, if the coherence magnitude falls below 0.4 (Fig. 7—red zone), the signal is considered unreliable for further phase analysis. As depicted in the figure, the coherence magnitude falls below 0.4 at approximately 400 m depth. As a result, for all pRES data analyzed in this study, only the top 400 m are used.
(a) Radar backscattered power (blue line) reveals the ice thickness. (b) The magnitude of complex polarimetric coherence between HH and VV signal (red line). The red zone is the area below 0.4 coherence magnitude.
Appendix B. 2D interpolated fabric spatial change
A 2D interpolated spatial distribution of fabric properties inferred from pRES data is provided in Fig. 8. The values depicted in Figs. 8a and 8b represent Δλ H and λ 3, respectively, directly estimated from the pRES data. On the other hand, Figs. 8c and 8d illustrate the deviation between the estimated ice fabric orientation $\vec {v}_2$
and the surface flow direction from SIA and between $\vec {v}_2$
and the maximum strain direction from SIA, respectively.
Showing the two dimensional interpolation of (a) horizontal ice fabric anisotropy. (b) Magnitude of the strongest eigenvalue (lambda 3). (c) deviation of $\vec {v}_2$
from surface flow direction. (d) deviation of $\vec {v}_2$
from maximum strain rate direction. Not that both X and Y axes are normalized by the mean ice thickness (H $\backsimeq 550$
m).
Appendix C. Woodcock plot (pRES and ice core)
Woodcock (Reference Woodcock1977) introduced the parameter K = ln(λ 3/λ 2)/ln(λ 2/λ 1) as a logarithmic ratio between the Eigenvalues, dividing the ice fabric type into the cluster zone (K > 1) and the girdle zone (K < 1). The extreme cases are the uniaxial girdle (K close to 0) and the uniaxial cluster (K close to infinity), with K = 1 representing the transition zone. Additionally, Woodcock introduced the parameter C = ln(λ 3/λ 1) as a measure of the preferred orientation strength. Higher C values indicate a greater concentration of the c-axis and a lower noise level. By using Woodcock's method, the ice fabric type obtained from estimated and measured Eigenvalues can be compared.
Here we regenerated the Fig. 1 from Woodcock (Reference Woodcock1977) and added some extra information to it. Hand-drawn Schmidt diagrams illustrate the shape of the ice fabric type in each zone, where the top left and bottom right show the uniaxial cluster and the uniaxial girdle, respectively. The isotropic ice fabric is situated at the origin of the figure. Not that the thin sections in Schmidt diagrams from the ice core analysis in Fig. 4c are vertical while the Schmidt diagrams shown in Fig. 9 are oblique. The estimated and measured ice fabric types are depicted as green squares and black circles, respectively, within the 50 to 260 m range. Both the estimated and measured ice fabric types suggest that the fabric is in the weak cluster zone, although the estimated fabric is slightly weaker compared to the measured fabric.
Regenerated Woodcock (Reference Woodcock1977), categorizing the ice fabric type according to Woodcock's parameters. The background color shows the change of Δλ H, green dashed contours show the Δλ V, blue dashed contours represent the K values, and red contours are the C values. The Schmidt diagrams are copied directly from Woodcock (Reference Woodcock1977). The green squares and black circles are estimated from radar data and measured from the ice core, respectively, between 50 to 260 m depth.
Appendix D. SIA results
The magnitude and orientation of surface velocity along with the magnitude and orientation of the maximum horizontal strain estimated from SIA as explained in section 3.5 are shown all over HIR in Fig. 10
Estimated from SIA, (a) magnitude and direction of surface velocity. (b) magnitude and direction of maximum horizontal strain rate.
Open access
