Characterizing sub-glacial hydrology using radar simulations

The structure and distribution of sub-glacial water directly influences Antarctic ice mass loss by reducing or enhancing basal shear stress and accelerating grounding line retreat. A common technique for detecting sub-glacial water involves analyzing the spatial variation in reflectivity from an airborne radar echo sounding (RES) survey. Basic RES analysis exploits the high dielectric contrast between water and most other substrate materials, where a reflectivity increase ≥ 15 dB is frequently correlated with the presence of sub-glacial water. There are surprisingly few additional tools to further characterize the size, shape, or extent of hydrological systems beneath large ice masses.

We adapted an existing radar backscattering simulator to model RES reflections from sub-glacial water structures using the University of Texas Institute for Geophysics (UTIG) Multifrequency Airborne Radar Sounder with Full-phase Assessment (MARFA) instrument. Our series of hypothetical simulation cases modeled water structures from 5 to 50 m wide, surrounded by bed materials of varying roughness. We compared the relative reflectivity from rounded Röthlisberger channels and specular flat canals, showing both types of channels exhibit a positive correlation between size and reflectivity. Large (> 20 m), flat canals can increase reflectivity by more than 20 dB, while equivalent Röthlisberger channels show only modest reflectivity gains of 8–13 dB. Changes in substrate roughness may also alter observed reflectivity by 3–6 dB. All of these results indicate that a sophisticated approach to RES interpretation can be useful in constraining the size and shape of sub-glacial water features. However, a highly nuanced treatment of the geometric context is necessary.

Finally, we compared simulated outputs to actual reflectivity from a single RES flight line collected over Thwaites Glacier in 2022. The flight line crosses a previously proposed Röthlisberger channel route, with an obvious bright bed reflection in the radargram. Through multiple simulations comparing various water system geometries, such as canals and sub-glacial lakes, we demonstrated the important role that topography and water geometry can play in observed RES reflectivity. From the scenarios that we tested, we concluded the bright reflector from our RES flight line cannot be a Röthlisberger channel but could be consistent with a series of flat canals or a sub-glacial lake. However, we note our simulations were not exhaustive of all possible sub-glacial water configurations.

The approach outlined here has broad applicability for studying the basal environment of large glaciers. We expect to apply this technique when constraining the geometry and extent of many sub-glacial hydrologic structures in the future. Further research may also include comprehensive investigations of the impact of sub-glacial roughness, substrate heterogeneity, and computational efficiencies enabling more complex and complete simulations.

1 Introduction

The size, shape, and distribution of sub-glacial water is important to ice dynamics and remains a significant uncertainty in projecting sea level rise due to ice mass loss. Hydrological structures directly influence basal shear stress distribution, which defines the boundary condition for rheology at the ice–bed interface (Gilbert et al., 2022; Brinkerhoff et al., 2021; Hoffman et al., 2016). Widely distributed water has been shown to lubricate the base and weaken sediments (Dunse et al., 2015; Hoffman et al., 2016). Conversely, narrow water channels may have little impact on shear stress at the bed (Schroeder et al., 2013) but act as conduits for concentrating meltwater produced upglacier. The size and location of such channels can directly influence grounding line retreat by controlling the volume of water transported to the grounding line (Young et al., 2016; Wright et al., 2012; Schroeder et al., 2013).

Airborne radar echo sounding (RES) is an established technique for studying sub-glacial hydrology throughout Earth's cryosphere (Peters et al., 2005; Chu et al., 2016; Schroeder et al., 2013). Attenuation loss through ice at common radar wavelengths (∼ 2–5 m) is relatively low, enabling reliable imaging of bed surfaces beneath ice masses several thousand meters thick. The dielectric contrast at an ice–water interface is much higher than at an ice–rock interface, resulting in higher bed reflectivity when water is present. Additionally, some ice–water interfaces may be relatively smooth, depending on the geometry of the hydrological structure. These properties are commonly exploited in RES surveys to infer the location of sub-glacial water based on reflected power ≳ 15 dB higher than the surrounding area (Schroeder et al., 2015; Peters et al., 2005; Rutishauser et al., 2018; Schroeder et al., 2013; Young et al., 2016).

Deducing the presence of water beneath large ice sheets using RES is relatively common; however, methods for testing further hypotheses regarding size, geometry, or hydrological structure remain challenging. Direct observation of the glacier bed over any significant spatial extent is infeasible with current methods (e.g., drilling (Priscu et al., 2021)), limiting our ability to calibrate radar returns to observed hydrological features.

This problem has precedent in inter-planetary science, where radar experiments are designed to test hypotheses with limited in situ evidence about surface or sub-surface characteristics. Backscattering simulators have proven especially useful in modeling radar returns for celestial targets (Spagnuolo et al., 2011; Russo et al., 2008; Gerekos et al., 2018). Most simulators approximate a target surface as a series of flat facets acting as point backscatterers. The backscattered electric field strength at the radar antenna position is estimated with common mathematical approximations, such as the Stratton–Chu integral. In such a point-scatterer formulation, facets much smaller than the radar wavelength (typically < $\mathit{\lambda }/\mathrm{10}$) are required to approximate real-world instrument results (Gerekos et al., 2018). This constraint necessitates access to high-end computing resources, often making point-scattering radar simulators unrealistic for RES modeling.

However, as algorithms improve and computational costs decrease, it is increasingly attractive to attempt such a simulation method with ice-penetrating RES problems. Gerekos et al. (2018) described a simulation technique that is particularly intriguing for the study of sub-glacial hydrology. The methodology is unique in two distinct ways from other simulators, which are helpful in modeling ice penetrating RES. First, the simulator can estimate the strength and direction of signals transmitted through multiple layered material interfaces. This makes it conducive to targets such as the ice–bedrock system. Second, the algorithm allows phase to vary linearly across the facet (termed the linear phase approximation or LPA). This feature enables modeling with significantly larger facets (∼λ or larger), drastically reducing computational resources needed for accurate simulations.

This paper demonstrates the radar simulator's application to geometric scattering from common hydrological targets: flat canals and Röthlisberger channels. We discuss the relevant parameters to achieve accurate model results and illustrate its utility in interpreting radar signatures from hydrological features beneath the ice. Finally, as an example, we demonstrate simulations of a hydrological target beneath Thwaites Glacier. In the future, we anticipate the technique will have broad applicability to interpreting RES data in the context of the sub-glacial environment.

Figure 1 Two-dimensional conceptual model of a RES simulation, defined in the direction of flight. (a) The material and geometric model for a radar observation. (b) The radar's spherical wavefront is approximated as rays directed at facets within a footprint beneath the aircraft. The simulator calculates a reflected and transmitted field strength based on the dielectric constant at the surface (ϵ_ice). We do not show the surface reflection for simplicity. (c) The rays propagate through the ice, with a direction obeying Snell's Law. At the bed, reflected field strength and direction are controlled by dielectric contrast and incidence angle with the basal facet. (d) The backscattered electric field strength (E_b) is approximated with the Stratton–Chu integral. (e) The aircraft's position is incremented, and the process is repeated.

2 Simulation methodology

Figure 1 describes a generic conceptual model for a RES radar simulation. Parameters were chosen to emulate a typical helicopter-based airborne radar survey with the University of Texas Institute for Geophysics (UTIG) Multifrequency Airborne Radar Sounder with Full-phase Assessment (MARFA) instrument (Castelletti et al., 2017; Lindzey et al., 2020). An ice surface and a bed surface are defined in a three-dimensional Cartesian coordinate system, noted as S_ice and S_bed, respectively. Figure 1a illustrates a 2-D representation of these surfaces, divided into N facets using a Delaunay triangulation algorithm, with characteristic length l_f. The linear phase approximation (LPA) employed in Gerekos et al. (2018) allows the phase of the incident and reflected electric fields to vary linearly across each facet. LPA enables accurate simulations of coherent radar using relatively large facets (∼λ vs. $\mathit{\lambda }/\mathrm{10}$ with other methods). The facet size should be constrained according to Eq. (1), where h is the aircraft height and λ is the free-space radar wavelength (Gerekos et al., 2018). For simulations presented here, h=500 m and λ=5 m, consistent with a typical UTIG MARFA helicopter experiment. l_f=5 m was chosen, consistent with Eq. (1).

$$\begin{array}{}\text{(1)}& l_{\mathrm{f}}\le \mathrm{0.2}\sqrt{\mathit{\lambda }h/\mathrm{2}}\end{array}$$

The radar's spherical wavefront is simulated as a series of plane waves, with wavevectors ${\stackrel{\mathrm{‾}}{\mathit{k}}}_i$ and field strength vector ${\stackrel{\mathrm{‾}}{\mathit{E}}}_i$ directed at each ice surface facet within a radius R beneath the aircraft (Fig. 1b). A wider R will provide a more complete approximation of the radiated and returned electric field but comes at a computational cost proportional to R². Therefore, R is chosen to balance these competing priorities, as we will discuss in detail.

Each ray's path is traced through transmission at S_ice using Snell's Law, and then reflection is calculated from S_bed (Fig. 1c). In Fig. 1, the reflection from the ice surface is omitted for simplicity, since we are most concerned with reflections from the bed in this study. Reflected and transmitted field strengths (${\stackrel{\mathrm{‾}}{\mathit{E}}}_{\mathrm{t}}$ and ${\stackrel{\mathrm{‾}}{\mathit{E}}}_{\mathrm{r}}$) are calculated from the real component of the dielectric constant at S_ice and S_bed, as discussed below. Attenuation loss between the surfaces is dependent on the complex component of the dielectric constant for ice (Gerekos et al., 2020), which is assumed constant over the short simulation distances presented here.

The total field at the radar antenna is then approximated as a summation of the backscattered electric fields (${\stackrel{\mathrm{‾}}{\mathit{E}}}_{\mathrm{b}}$) from individual propagated waves using the Stratton–Chu integral (Fig. 1d). A detailed treatment of the simulation mathematics is presented in Gerekos et al. (2018). Because of higher ice–water dielectric contrast, rays reflecting from a facet within S_bed identified as water will have higher amplitude than facets comprised of rock.

The radar pulse repetition frequency (PRF) and aircraft velocity v determine the spatial resolution between radar observations in the azimuth direction, Δ_a, as shown in Fig. 1e. In the field, the MARFA instrument employs a native pulse repetition frequency (PRF) of 6400 Hz and peak power of 8 kW, then these observations are coherently stacked in post processing to achieve along-track resolution of 1 m before focusing (Peters et al., 2007). We chose to simulate the final (stacked) PRF and power instead of the native parameters for computational efficiency.

Table 1 Summary of RES simulation parameters. Parameters were chosen to replicate helicopter-based MARFA ice sounding experiments.

^∗ For fields marked with an asterisk, the native MARFA radar value differs from the parameter used in the simulation. Simulation parameters are in parentheses and represent values after coherent pre-summing of raw MARFA data.

2.1 Dielectric material model

It is instructive to consider energy returned to the radar as the superposition of dielectric and geometric effects. Dielectric effects result from material property changes at an interface, while geometric effects result from the orientation of the interface's topography and the radar antenna. We will begin by discussing the implications of dielectric parameters then move on to geometric considerations.

A material model is applied to the simulation, where all facets on S_ice and S_bed are assigned a complex dielectric constant, ϵ. The refractive index η for each material is derived from ϵ (Eq. 2). For simplicity, we have chosen a three-material model consisting of ice, rock, and water, with dielectric constants as presented in Table 2. Facets across S_ice are assigned ϵ_ice. Facets on S_bed are assigned ϵ_rock unless they are part of a water structure, which are assigned ${\mathit{ϵ}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$.

Table 2 Complex dielectric constants for materials in radar simulations (Peters et al., 2005; Fujita et al., 2000; Midi et al., 2014; Glover, 2015). For each material, we chose a value within the range published in the literature. Further analysis of sensitivity to these material parameters is presented in the discussion.

For a facet at nadir, absolute reflectivity R_abs results from the contrast between the refractive indices of ice and the bed material (η_ice and η_rock, respectively). A horizontal facet assigned ϵ_rock will have a lower reflectivity than one assigned ${\mathit{ϵ}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ by about 15 dB (Eq. 2). Real RES instruments are rarely calibrated to measure absolute reflectivity, and thus changes in received power, corrected for ice attenuation loss and geometric scattering, are assumed proportional to relative changes in reflectivity at the bed (R_rel). This forms the basis for the widely accepted assumption that R_rel≥ 10–15 dB over surrounding reflections implies the presence of liquid water (e.g., Young et al., 2016; Peters et al., 2007; Rutishauser et al., 2022).

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& R_{\text{abs},m}={\left({\displaystyle \frac{{\mathit{\eta }}_{\text{ice}}-{\mathit{\eta }}_m}{{\mathit{\eta }}_{\text{ice}}+{\mathit{\eta }}_m}}\right)}^{\mathrm{2}},\\ & {\mathit{\eta }}_m=\sqrt{Re\left({\mathit{ϵ}}_m\right)},m=\mathrm{rock},{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\end{array}\end{array}$$

The simulations presented in this paper assume only the three materials described in Table 2, with well-defined boundaries between hydrological and bedrock features. In real-world sub-glacial environments, additional material heterogeneity such as clays or hydrated tills exist, in addition to ambiguity in hydrological boundaries. We address the implications of this relatively simplistic material model in Sect. 5.

Figure 2 Radar simulation geometry of the (a) flight path (red), aircraft or channel orientation, and simulation radius (R) across the model footprint. (b) A flat canal with width c_w. (c) Geometry of a Röthlisberger or round channel with width c_w.

2.2 Geometric model

We ultimately seek to develop the Stratton–Chu simulation method for constraining the extent and cross-sectional geometry of sub-glacial water features. Given this objective, it imperative to consider a menu of geometric constraints and how simulation parameters will emulate real-world RES returns from different targets. A basic simulation geometry for a hypothetical case is shown in Fig. 2a. These simulations consist of flat surfaces S_ice at elevation 0 m and S_bed at −d_ice, where d_ice is the nominal ice thickness in m. The flight path is defined along the y direction, with the radar's dipole antenna oriented along the x direction. On the bed surface, a straight channel of width c_w is oriented perpendicular to the flight path.

2.2.1 Surface roughness

Small-scale topography variation below the radar's detection limit, which we will refer to as roughness, can impact reflectivity by diffusely scattering incident radar energy. To account for this effect, random isotropic Gaussian variation is introduced to both S_ice and S_bed via Eq. (3). l_c is the correlation length in both the x and y directions. To capture scattering behavior due to topography changes at the radar wavelength scale, l_c should be at least a few times the facet length l_f, while remaining near λ. In our simulations, both l_f and λ are ∼ 5 m, leading us to select l_c = 15 m.

$$\begin{array}{}\text{(3)}& {\mathrm{cov}}_{x,y}={\mathit{\sigma }}_m{e}^{-(x^{\mathrm{2}}+y^{\mathrm{2}})/l_{\mathrm{c}}^{\mathrm{2}}},m=\mathrm{ice},\mathrm{bed}\end{array}$$

When we assume our surfaces are smooth relative to λ, we must add at least a negligible, non-zero roughness to avoid simulation artifacts, such as Bragg resonance (Gerekos et al., 2020). We always make these assumptions for S_ice; therefore, σ_ice = 0.2 m in all simulations presented here. We also use σ_bed = 0.2 m in “smooth” bed simulations.

Hubbard et al. (2000) measured topography of recently deglaciated bedrock at high resolution in the Swiss Alps. They showed variations near 1 m over horizontal distances of 15 m, which serves as our upper bound on σ_bed. We acknowledge this Alpine sub-glacial environment may differ considerably from common RES survey targets such as Thwaites Glacier in Antarctica. However, roughness studies from the Thwaites region integrate over long horizontal scales irrelevant to radar scattering (Bingham and Siegert, 2009; Hoffman et al., 2022). We proceed with the understanding that our range for σ_bed from 0.2 to 1 m may represent an imperfect but reasonable range of expected variation in small-scale topography.

2.2.2 Simulation radius

An appropriate choice of simulation radius R is vital to accurately simulate radar geometric scattering from the sub-glacial environment. R defines the radius of a vertical cylinder, bounding the simulation scope at each aircraft position increment. If R is sufficiently large to capture the entire antenna beam pattern cast on the bed surface, then the simulation will provide a complete representation of off-nadir clutter and target range migration. In many layered radar simulation experiments, R is chosen primarily to capture off-nadir clutter at the sub-surface target's apparent depth (Gerekos et al., 2018). Our experiments involve both a thick ice material layer and smooth S_ice relative to λ, limiting the impact of off-nadir clutter. Therefore, we base our choice for R on two alternative criteria.

• R must be greater than the pulse-limited radius R_pl at the glacier bed.

• R captures adequate range migration to facilitate along-track focusing.

The first criterion assumes the majority of returned energy from a nadir-directed radar will come from within the pulse-limited footprint (R_pl) beneath the aircraft. Equation (4) approximates an upper limit on R_pl for a radar instrument with bandwidth B_w, where c is the speed of light in a vacuum. Actual R_pl will always be smaller than Eq. (4) predicts, due to refraction at the air–ice interface. For our simulations, with d_ice∼1000 m, R_pl≈ 130 m (Eq. 4). We consider this the minimum acceptable simulation radius, although clearly changes in aircraft height or simulated ice thickness will alter this limit.

$$\begin{array}{}\text{(4)}& R>R_{\text{pl}}\approx \sqrt{{\displaystyle \frac{c(d_{\text{ice}}+h)}{B_{\mathrm{w}}{\mathit{\eta }}_{\text{ice}}}}}\end{array}$$

An appropriate choice for R must also consider the desired range cell migration (RCM) at the bed surface. RCM is proportional to the change in physical distance a signal travels through air (r_air) and ice (r_ice) to reach a target as the radar moves past (Eq. 5). In synthetic aperture radar (SAR) processing, an aperture length L_a is chosen with sufficient range cell migration to optimize along-track focusing. This process improves signal-to-noise ratio and along-track resolution (Cumming and Wong, 2005). In order to facilitate simulated data focusing, R must be greater than the aperture required for desired range migration. Selection of L_a and the focusing process is described at length in Sect. 2.3.

$$\begin{array}{}\text{(5)}& \mathrm{RCM}\left(y\right)={\displaystyle \frac{\mathrm{2}{f}_{\mathrm{s}}}{c}}\left[\right(r_{\text{air}}\left(y\right)-h)+{\mathit{\eta }}_{\text{ice}}(r_{\text{ice}}\left(y\right)-d_{\text{ice}}\left)\right]\end{array}$$

Our simulations target three cells of range migration at the bed surface (RCM(L_a)≥3). For d_ice = 1000 m and sampling frequency f_s = 50 MHz, this translates to L_a = 277 m. All simulations presented here use R=300 m, meeting both the pulse-limited and range migration criteria. Real airborne RES focusing typically includes more range cell migration (RCM(L_a)=5 fast-time samples for MARFA); however, such a simulation radius would be computationally unrealistic. Therefore, our choice of R represents a compromise, and the implications will be discussed in Sect. 5. Simulations in thicker ice may require significantly larger R or further compromises in range migration.

2.2.3 Channel geometry

We seek to distinguish between two-channel geometries common to sub-glacial hydrology. A canal-like structure (Fig. 2b) will have a flat cross-section and produce a specular reflection. This type of feature is common when the surrounding bed is comprised of sediment or other soft material and the water pressure is high. Conversely, in a Röthlisberger channel (Fig. 2c), sub-glacial water carves a path through the ice above the bed surface (Röthlisberger, 1972). This type of channel is likely to form in areas where the substrate is impermeable bedrock with low water pressure (Walder and Fowler, 1994; Schroeder et al., 2013). We have confined ourselves to Röthlisberger channels and flat canals in this study as examples of common large-scale hydrological structures. Other known hydrologic features, such as Nye channels, may appear radiometrically similar to small canals, depending on their size and the radar wavelength. As algorithms and computational power enable higher-resolution models, additional sub-glacial structures are a logical extension of this work.

A Röthlisberger channel has an elliptical cross-section, with height c_h in addition to channel width c_w. We can infer from basic geometry that it will reflect radar energy divergently (Fig. 2c); therefore, the actual radar signature of a sub-glacial channel will be the superposition of its geometry and dielectric contrast. When representing the Röthlisberger curvature with flat facets, we must approximate the divergent scattering behavior by capturing reflections from multiple facets in the channel cross-section at nadir. In our simulations, we seek to capture reflections from at least the upper six facets within the simulation footprint as a reasonable estimate of Röthlisberger channel scattering (Fig. 3a).

Figure 3 (a) Schematic of simulation geometry for an elliptical cross-section of the Röthlisberger channel at nadir. A second, smaller facet length scale l_f2 is introduced to facilitate accurate representation of curvature and divergent reflection character. (b) Relationship between R and l_f2 is defined for a range of c_w. We seek to limit R≤300 m for computational efficiency, imposing an upper limit on l_f2 when simulating Röthlisberger channels

Given our desired range of c_w and the constraint already imposed on R, it is clear that l_f=5 m is too large to capture the scattering effect of a Röthlisberger channel. For our smallest channels (c_w = 5 m), l_f≤ 0.25 m is required to capture the desired facet reflections (Fig. 3b). Setting such a high resolution over all of S_ice and S_bed is computationally unrealistic and negates many of the benefits of the Stratton–Chu simulation method (Gerekos et al., 2018). Therefore, we introduce a second facet length scale, l_f2, which defines the facet length on S_ice and S_bed only in a Röthlisberger channel location. l_f2 has a maximum of $\sim c_{\mathrm{w}}/\mathrm{20}$, for d_ice = 1000 m and $c_{\mathrm{h}}=c_{\mathrm{w}}/\mathrm{4}$ (Fig. 3b). When simulating flat canals, higher-resolution channel facets are unnecessary due to the specular nature of the reflection, and therefore l_f2=l_f.

2.3 Simulated data processing

The simulator outputs rangelines, which represent the electric field strength vs. fast-time τ returned from a single pulse, at azimuth time a. τ is discretized into range bins (index j) with increments of $\mathrm{1}/f_{\mathrm{s}}$, where f_s is the sampling frequency. The rangelines are compiled sequentially in azimuth, or slow time, into a 2-D raw radargram matrix (ξ_raw). Azimuth increments (index k) are equally spaced in time at $\mathrm{1}/\mathrm{PRF}$.

The rangelines are then focused using a version of the range Doppler algorithm (RDA) (Cumming and Wong, 2005; Hélière et al., 2007). We first perform range compression by convolving each raw rangeline with the complex conjugate of the radar chirp, g(τ) over the radar receiving window T_r. The operation is performed by multiplication in the fast-time frequency domain.

$$\begin{array}{}\text{(6a)}& {\displaystyle }& {\displaystyle }g\left(\mathit{\tau }\right)=e^{-i\mathit{\pi }{B}_{\mathrm{w}}{\mathit{\tau }}^{\mathrm{2}}/T_{\mathrm{s}}}\text{(6b)}& {\displaystyle }& {\displaystyle }{\mathit{\xi }}_{\text{RC}}(a,\mathit{\tau })=\underset{\mathrm{0}}{\overset{T_{\mathrm{r}}}{\int }}{\mathit{\xi }}_{\text{raw}}(a,\mathit{\tau }+{\mathit{\tau }}^{\prime })g(\mathit{\tau }{)}^{\ast }\mathrm{d}{\mathit{\tau }}^{\prime }\end{array}$$

We produce the focused radargram (ξ_f) by convolving the range-compressed radargram (ξ_RC) with a 1-D reference function (ϕ) in along-track blocks (Eq. 7). The block size, L_a, for each fast-time value of τ is chosen such that range migration equals three fast-time samples. Thus it is important to note that L_a increases with depth, and the simulation radius R must be greater than the maximum anticipated L_a.

$$\begin{array}{}\text{(7)}& {\mathit{\xi }}_{\mathrm{f}}(a,\mathit{\tau })=\underset{-L_{\mathrm{a}}/\mathrm{2}}{\overset{L_{\mathrm{a}}/\mathrm{2}}{\int }}{\mathit{\xi }}_{\text{RC}}(a+a^{\prime },\mathit{\tau })\mathit{\varphi }(a^{\prime },\mathit{\tau }{)}^{\ast }\mathrm{d}{a}^{\prime }\end{array}$$

A deeper discussion of mathematics and block processing required for along-track focusing is presented in Appendix A.

In all radargrams, we convert fast time to physical depth (d) and slow time to along-track distance (y) via Eqs. (8) and (9), where ${\mathit{\tau }}_{\mathrm{s}}=\mathrm{2}h/c$ is the fast-time value for the surface reflection and PRF is the radar pulse repetition frequency.

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{d}_j={\displaystyle \frac{c}{\mathrm{2}{\mathit{\eta }}_{\text{ice}}}}({\mathit{\tau }}_j-{\mathit{\tau }}_{\mathrm{s}}),{\mathit{\eta }}_{\text{ice}}=\sqrt{Re\left({\mathit{ϵ}}_{\text{ice}}\right)}\text{(9)}& {\displaystyle }& {\displaystyle }{y}_{\mathrm{k}}=k{\displaystyle \frac{v}{\mathrm{PRF}}}\end{array}$$

The focused field in ξ_f is converted to power in decibels. We define along-track bed reflectivity in absolute terms (R_abs) by taking the maximum reflected power beneath the ice surface (Eq. 10a). Relative reflectivity (R_rel) compares R_abs to the mean reflectivity for a simulation with the same σ_bed and σ_ice, but no channel is present (R_abs,0). R_rel therefore measures the relative reflectivity gain observed by the radar due to the channel's presence vs. surrounding frozen-bed material.

$$\begin{array}{}\text{(10a)}& {\displaystyle }& {\displaystyle }{R}_{\text{abs}}\left(y_{\mathrm{k}}\right)=max\mathit{\left\{}{\mathrm{RG}}_{f,\mathrm{dB}}\right(d,y_{\mathrm{k}}\left)\mathit{\right\}},d>\mathrm{0}\text{(10b)}& {\displaystyle }& {\displaystyle }{R}_{\text{rel}}\left(y_{\mathrm{k}}\right)=R_{\text{abs}}\left(y_{\mathrm{k}}\right)-\mathrm{mean}\mathit{\left\{}{R}_{\text{abs},\mathrm{0}}\mathit{\right\}}\end{array}$$

2.4 Hypothetical simulation cases

A basic simulation geometry is shown in Fig. 2. Our hypothetical simulation cases consisted of flat surfaces S_ice and S_bed with elevations of 0 m and −1000 m, respectively. Each surface had isotropic Gaussian roughness as defined above. On the bed surface a single channel of width c_w was oriented perpendicular to the flight path. Channel cross-sections were either flat canal-like structures (Fig. 2b) or round Röthlisberger channels with a channel height of $c_{\mathrm{w}}/\mathrm{4}$ (Fig. 2c).

A series of simulation experiments were run for both types of channels. Channel width and basal roughness were varied according to Table 3. Simulations involving flat canals used a single facet length of 5 m, while Röthlisberger channels necessitated a smaller l_f2 in the channel location to accurately represent geometric scattering. Rangelines from each simulation were processed as described above, and along-track bed reflectivity from the focused radargram R_rel was compared for various scenarios.

Table 3 Summary of RES simulation parameters for each channel geometry.

3 Results – hypothetical simulation cases

Basal roughness (σ_bed) impacted the absolute reflectivity of the solid bed material in our simulations (R_abs,0), which we use as a baseline for calculating the relative impact of channels on the radar echo (Eq. 10b). Simulations excluding liquid water, with σ_bed of 0.2 and 1 m, produced mean R_abs,0 equal to −107.9 ± 2.5 and −114.1 ± 2.4 dB, respectively.

Simulations of hypothetical Röthlisberger channels and flat canals with the same c_w demonstrate the distinctly different radar signatures expected for the two geometries. Figure 4 shows along-track relative reflectivity for the two channel types, with c_w = 20 m, compared to a frozen substrate with no liquid water. Flat canals have a distinct peak centered near the canal location, with maximum R_rel = 19.5 dB when σ_bed = 0.2 m and 22.9 dB when σ_bed = 1 m. Regardless of roughness, a 20 m flat canal influences the radar echo for a few hundred meters along track. The reduced gain of 3.4 dB for a canal surrounded by a smoother substrate is consistent with the higher absolute reflectivity of the surrounding surface as discussed above.

Figure 4 Simulated along-track R_rel for 20 m channels compared to a dry substrate with (a) σ_bed=0.2 m and (b) σ_bed=1 m. All simulations had the same Gaussian roughness correlation length of l_c=15 m.

When flowing through a rough bed (Fig. 4b), a 20 m simulated Röthlisberger channel increased along-track R_rel for more than 500 m. The maximum peak in R_rel for this channel is 10.3 dB, occurring when the channel is far from nadir. In our smoother bed simulations, the same 20 m Röthlisberger cross-section produced peak R_rel of only 7.6 dB (Fig. 4a), which could be nearly indistinguishable from fluctuations in reflectivity from the frozen substrate. Larger Röthlisberger channels produce only moderate gains in R_rel, with the largest (50 m) channels having peak R_rel = 13.2 dB (Fig. 5).

Figure 5 Positive correlation between c_w and R_rel for both Röthlisberger channels and flat canals. The increase in R_rel is much more pronounced for flat canals.

Figure 5 shows peak R_rel generally increases with width for both channels and canals, as an increasing proportion of the area within the radar's footprint contains liquid water. Flat canals exhibit a much stronger correlation between radar reflectivity response and c_w than Röthlisberger channels when c_w < 20 m. At c_w > 20 m, R_rel for flat canals increases more gradually as it appears to reach an asymptote of ∼ 24 dB.

4 Application to Thwaites Glacier

To demonstrate the simulator's utility in sub-glacial hypothesis testing, we compared simulated reflectivity to a single flight line collected with UTIG's 60 MHz MARFA instrument aboard an AS-350 B2 helicopter. The 16 km line (THW2/UBH0c/X243a) was part of an airborne radar survey of Thwaites Glacier conducted between 2020 and 2022. The helicopter was flown at a nominal height of 500 m above ground level with target velocity of 30 m s⁻¹. Precise aircraft positioning and orientation were recorded with an onboard Renishaw laser altimeter, Trimble Net-R9 dual-frequency GNSS, and Novatel SPAN IGM-1A inertial navigation, as described in Lindzey et al. (2020).

Figure 6 The grounding line region of Thwaites Glacier, with the region of channelized hydrology proposed by Schroeder et al. (2013) highlighted. UTIG MARFA radar survey lines collected between 2020 and 2022 are shown in grey. Likely high volume channel routes suggested by Hager et al. (2022) are shown in blue, and the location of THW2/UBH0c/X243a is highlighted for reference.

4.1 Study area

Schroeder et al. (2013) first postulated a hydrological system dominated by channelized drainage within ∼ 45–70 km of the Thwaites grounding line (Fig. 6) based on data interpretation from the extensive 2003–2004 Airborne Geophysics of the Amundsen Sea Embayment, Antarctica (AGASEA) radar survey. The hypothesis was largely derived from observations of high local bed reflectivity, combined with anisotropy in the radar specularity content. Specularity content measures the relative diffuse or specular nature of radar returned power. Schroeder et al. (2013) postulated that extensive Röthlisberger channels could produce such radiometric signatures due to their geometric orientation.

Hager et al. (2022) ran a suite of sub-glacial hydrology simulations to evaluate the probability of persistent channelization routes beneath Thwaites. Their analysis concluded that Thwaites' near terminus hydrology is most likely comprised of a few persistent, high volume channels flowing toward the central grounding zone. Two probable routes were proposed, also shown in Fig. 6.

Although evidence for sub-glacial hydrology dominated by channelization in this region of Thwaites Glacier is strong, structures such as canals or lakes are not unprecedented. For example, Smith et al. (2017) identified a series of four active sub-glacial lakes that drained between June 2013 and January 2014. As seen in Fig. 6, at least one of these lakes is within the channelized region as defined in Schroeder et al. (2013).

Our chosen flight line, THW2/UBH0c/X243a, transects one of the proposed Hager et al. (2022) channel routes, and the radargram shows an isolated bright reflection coincident with this location, as shown in Fig. 7. Data from THW2/UBH0c/X243a were range compressed, corrected for geometric spreading loss and aircraft position, and focused in azimuth as described in Peters et al. (2007). This azimuth focusing is analogous to the procedure described in Sect. 2.3, although a longer aperture sufficient for range migration of five cells is used.

Figure 7 Focused radargram for THW2/UBH0c/X243a. The bright bed reflection corresponding to the proposed channel location is boxed. Vertical lines represent the extent of the flight line reproduced in the simulations.

The along-track surface and bed profiles were picked from the radargram. One-way attenuation loss of 13.8 ± 1.4 dB km⁻¹ was estimated using a spatially constrained linear regression model as outlined in Schroeder et al. (2016).

4.2 THW2/UBH0c/X243a simulation methodology

We simulated a 4 km segment of THW2/UBH0c/X243a containing the proposed channel, as depicted in Fig. 7. Along-track ice and surface elevations at nadir were calculated from the THW2/UBH0c/X243a focused radargram. From these data, we built elevation matrices for both the ice and bed surfaces ($S_{\text{ice}}^{\text{rad}}$, $S_{\text{bed}}^{\text{rad}}$) which vary in y according to the respective radar profiles but have no x variation.

To build appropriate across-track (x) topography, separate 2-D topographic matrices were created from BedMachine V2 data ($S_{\text{ice}}^{\text{BM}}$, $S_{\text{bed}}^{\text{BM}}$) (Morlighem, 2020a, b). The radar- and BedMachine-derived topography were superimposed to create simulation surfaces S_ice and S_bed via Eq. (11), where w is a quadratic weighting function varying between 1 at the surface edges and 0 at nadir. S^rough is an isotropic Gaussian surface, adding random roughness with σ = 0.2 m over correlation length l_c = 15 m.

$$\begin{array}{}\text{(11)}& S=w{S}^{\text{BM}}+(\mathrm{1}-w)S^{\text{rad}}+S^{\text{rough}}-z_{\mathrm{0}}\left(y\right)\end{array}$$

In an ice-penetrating RES survey, the aircraft attempts to “drape” the ice surface by flying at a constant height above ground level (500 m for UTIG helicopter based surveys). We simulate this with a polynomial interpolation of the radar-derived ice elevation along track, z₀(y). For our simulations of THW2/UBH0c/X243a, z₀ is a seventh-order polynomial, although in practice the polynomial order is somewhat subjective. The fitted function should approximate major terrain features in $S_{\text{ice}}^{\text{rad}}$ with gentle elevation changes, consistent with actual aircraft operation.

z₀(y) is subtracted from S_ice and S_bed, setting the average ice surface elevation to ∼ 0 m (Eq. 11). The simulated aircraft elevation h is a constant 500 m. This approach preserves known ice geometry at nadir, and minor topographic features appear as variations in aircraft range to target.

Figure 8 (a) A 4 km region from actual THW2/UBH0c/X243a radargram containing the Hager et al. (2022) channel location. (b) Simulated THW2/UBH0c/X243a radargram with no water, σ_bed = 0.2 m. (c) Bed relative reflectivity from the actual THW2/UBH0c/X243a. The proposed channel location is highlighted in red as the area of interest. (d) Simulated THW2/UBH0c/X243a relative reflectivity with no water, σ_bed = 0.2 m. The along-track location corresponding to the area of interest is highlighted in red.

The dielectric material model was applied to S_ice and S_bed as described previously. We ran individual simulations using the same surfaces, varying the width and geometry of across-track-oriented channels in the location identified in Hager et al. (2022). Relative bed reflectivity from each simulation result was compared to the actual relative reflectivity from the focused THW2/UBH0c/X243a radargram (Fig. 8c). The comparison between simulated and real data provides constraints on the extent and geometry of any real hydrological features at this location.

4.3 THW2/UBH0c/X243a simulation results

Figure 8 compares a 4 km segment from THW2/UBH0c/X243a RES data with a simulation containing no water and σ_bed = 0.2 m. There are two bright reflections in this section of the radargram (Fig. 8a). The first is centered around 1100 m along track with a maximum R_rel = 11.6 dB. The second is a broad area of high reflectivity between 2300–2900 m, with R_rel peaks ranging from 13 to 15 dB (Fig. 8c). This reflector coincides with the location of persistent Röthlisberger channelization proposed in Hager et al. (2022) and is the primary area of interest along THW2/UBH0c/X243a for this study.

The radargram from our frozen-bed simulation (Fig. 8b) captures the basic bed topography well, but along-track reflectivity is not always aligned with the real THW2/UBH0c/X243a radargram. Simulated reflectivity near 1100 m is consistent with the real data, indicating that this reflectivity peak could be the result of geometric effects from topography, rather than dielectric contrast from a water feature. However, for the region between 2300–2900 m, mean simulated R_rel is 22.8 dB below the value observed in the real THW2/UBH0c/X243a data (Fig. 11). In our material model, a gain of this magnitude could only be consistent with a change in dielectric properties from rock to liquid water. It is also unlikely to be a Röthlisberger channel, since we have demonstrated such geometry is not conducive to increased reflectivity of more than 13.2 dB for very large channels. We therefore propose that this basal reflector at 2300–2900 m along-track could represent a specular hydrological structure, such as a flat canal, of unknown dimensions.

When a single flat canal with c_w = 20 m was added to the simulation, we observe a peak R_rel of 6.8 dB over a narrow 150 m range along track (Fig. 9). This reflectivity gain is consistent with our findings from hypothetical simulation cases described above. However, the gain in reflectivity and along-track extent are insufficient to match the reflectivity profile observed in the real THW2/UBH0c/X243a RES data (Figs. 9b and 8c). Therefore, we conclude that the reflector cannot be a narrow, isolated flat canal.

Figure 9 (a) Simulated radargram with a single 20 m canal at the bed. (b) Relative reflectivity for simulation with 20 m canal at the bed.

The above results imply the reflector between 2300–2900 m may be a broader hydrological feature, such as an area with multiple flat canals. This hypothesis is compatible with the topographic context, given that it is in a low-lying area with steep topography just down-glacier. This is an ideal location for till and liquid water to accumulate if the ice exists at its pressure melting point.

Figure 10 R_rel for distributed hydrology simulations at 2300–2900 m along track. (a) A 600 m wet bed surface (Weertman sliding), (b) 30 flat canals with c_w=10 m, (c) 8 flat canals with c_w=30 m, (d) 6 flat canals with c_w=50 m, and (e) the entire 600 m area covered with a flat, specular waterbody (e.g., sub-glacial lake).

Figure 11 compares the difference between THW2/UBH0c/X243a and simulated R_rel for five simulations of such hydrological features. We compare both the mean and peak difference in R_rel between 2300–2900 m and conclude that each more closely approximates THW2/UBH0c/X243a than our frozen bed or single 20 m channel simulation. First, we consider the glacier sliding condition originally proposed by Weertman (1964). In this scenario, a thin water film (perhaps a few centimeters thick) coats the basal interface. We model this by maintaining the basal topography with σ_bed = 0.2 m, but all facets from 2300 m $\le y\le$ 2900 m are assigned ${\mathit{ϵ}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$. This simulation has a broad increase in reflectivity in the area of interest (Fig. 10a); however, the mean R_rel between 2300–2900 m is 9.4 dB below the actual THW2/UBH0c/X243a value (Fig. 11).

Figure 11 Mean difference between actual and simulated R_rel over the area of interest at 2300–2900 m for all simulations of THW2/UBH0c/X243a. Positive error bars show the difference between mean and peak R_rel over 2300–2900 m.

Figure 10b shows simulated R_rel when 30 narrow, flat canals (c_w = 10 m) were placed over the same 600 m region along track. This creates a region where 50 % of the area is covered with 10 m flat canals between areas of bed material with smooth basal roughness (σ_bed = 0.2 m). Mean R_rel at 2300–2900 m for this simulation is slightly higher than the Weertman film scenario in Fig. 10a but still 8.2 dB below the target R_rel from THW2/UBH0c/X243a.

The third simulation in Fig. 10c has eight larger flat canals (c_w = 30 m) evenly spaced across the same 600 m region. This simulation has less water coverage than the 30 m × 10 m flat canal simulation (40 % vs. 50 %), yet the wider channels increased mean R_rel by 3.5 dB (Fig. 11). This result reinforces that the shape of hydrological feature, not just the extent, makes a significant difference to the resulting reflectivity profile. Mean R_rel for the 8 m × 30m simulation was 4.8 dB below actual values. When the geometry is changed to six canals of 50 m (Fig. 10d), mean R_rel improves to just 2.3 dB below THW2/UBH0c/X243a.

The final simulation includes a very broad area of specular (σ_bed = 0 m) water covering the bed between 2300–2900 m. Due to its size, we refer to this feature as the “600 m lake” in Figs. 10 and 11. Mean R_rel over 2300–2900 m for this simulation deviates by only 0.6 dB from the actual THW2/UBH0c/X243a data. However, this 600 m lake simulation includes several peaks as high as 18 dB, which is 3 dB higher than the maximum peaks observed in THW2/UBH0c/X243a. The maximum R_rel from both the 8 m × 30 m and 6 m × 50 m simulations more closely matched the peaks observed in THW2/UBH0c/X243a.

We infer our feature at 2300–2900 m could be a wide area of sub-glacial water. Of the scenarios we tested, simulated reflectivity was most consistent with canals averaging ∼ 30–50 m, with at least 50 % water coverage, or a sub-glacial lake. Although THW2/UBH0c/X243a and our simulations demonstrate similar reflectivity profiles, we have not explored a comprehensive range of bed materials or water geometries. In Sect. 5, we address these omissions and the resulting limitations.

5 Discussion

Our comparison of Röthlisberger vs. flat canal reflectivity implies that a real world Röthlisberger channel, even one of significant size, may not exhibit an obvious reflectivity increase. Röthlisberger channels will scatter energy divergently (Fig. 2c); therefore, we see a smaller magnitude impact to R_rel over a greater distance than an equivalent flat canal. The simulator also demonstrates the radar's sensitivity to relatively narrow, flat canals. A single canal ∼ 10 m can exhibit a reflectivity peak > 15 dB (Fig. 5, nFlat canal). Such a hydrological feature would cover < 5 % of the radar's pulse-limited footprint under 1000 m of ice.

In our simulations of hypothetical cases, increasing σ_bed from 0.2 to 1 m reduced radar echo power by 6.2 dB. These results are intuitive. Scattering loss from a rougher surface will be larger and more variable, resulting in lower absolute reflectivity. Increasing roughness also alters the geometric contrast between the hydrological feature and surrounding material. The effect is most pronounced for flat canals with c_w = 10–20 m. The impact of the bed roughness near flat canals diminishes as c_w exceeds 30 m, when reflected energy from the water feature dominates the radar return.

Our examination of roughness in this study was deliberately limited, as the relationship between roughness and scattering is complex and dependent on correlation length l_c. Bingham and Siegert (2009) and Peters et al. (2005) each calculated wide ranges in RES-derived roughness across Antarctica. Roughness in these studies was measured with l_c > 1 km. Our method for building simulated ice and bed surfaces directly incorporates RES observations of topography with resolution l_f = 5 m (see Sect. 4.2). Therefore, any relationship between scattering and large scale roughness, as defined in these previous studies, is captured explicitly in our simulations.

Smaller-scale sub-glacial roughness (l_c < a few wavelengths) is more challenging to measure directly; therefore, actual values are poorly constrained. Our hypothetical case simulations demonstrate bed roughness over l_c = 15 m may alter radar echoes without a change in dielectric properties. This is consistent with Peters et al. (2005), who calculated theoretical scattering loss due to roughness could exceed 20 dB when l_c≪ first Fresnel zone. A comprehensive examination of basal roughness, incorporating a broad range of both l_c and σ_bed, is possible with our simulation method. Such an investigation is an active area of interest for future work. Our limited current examination of roughness demonstrates that RES detection of sub-glacial water must be nuanced. Large changes in R_rel may constitute a liquid water signature but could also indicate spatial heterogeneity in roughness or substrate material.

Given the limited reflectivity response from Röthlisberger channels, more advanced analysis techniques may be required to detect and characterize them accurately. This could include examination of specularity content as described in Schroeder et al. (2015). Specularity content quantifies the diffuse or specular quality of radar reflections by comparing returned power over two different SAR apertures. Schroeder et al. (2013) argued that specularity content from Röthlisberger channels will be orientation dependent; high values are expected if the radar path is aligned along the channel, and lower values are expected if the path is transverse to it. By contrast, flat interfaces such as canals and lakes should exhibit high and isotropic specularity content. Although our simulation methodology produces a focused product, it requires a simulation radius R greater than the SAR aperture L_a. Current computational constraints limit our choice of R to well below the aperture length of 2 km recommended in Schroeder et al. (2015). Therefore, although simulating specularity content would be a valuable advance, we have left it to a future study.

Our THW2/UBH0c/X243a simulations from Thwaites Glacier demonstrate the reflectivity magnitude may be consistent with a wide, flat water feature surrounded by bedrock. Although we did not simulate specularity content, we observe a prominent peak in actual specularity content from THW2/UBH0c/X243a (Fig. 12). This observation could support the hypothesis that the reflector at 2300–2900 m is a series of flat canals or a lake, since THW2/UBH0c/X243a is transverse to the channelized flow hypothesized in Hager et al. (2022) (Fig. 6). An orthogonally oriented radar transect could enable a more definitive diagnosis of the hydrological feature.

Figure 12 Actual specularity content from THW2/UBH0c/X243a.

Further, while our flat canal or lake configurations appear likely, we cannot preclude all competing hypotheses without exploring many additional scenarios beyond the scope of this work. Other water geometries, such as Nye channels or water sheets (Creyts and Schoof, 2009), combined with a broader range of roughness length scales, could fit the original THW2/UBH0c/X243a reflectivity profile. Future work should identify computational or process efficiencies to allow comprehensive coverage of all possible simulation variables.

The material model in this study was also deliberately simple. Many possible bed materials exist, ranging from frozen bedrock to clays with real dielectric constants between 2.7–36 (Fig. 13) (Christianson et al., 2016; Tulaczyk and Foley, 2020; Peters et al., 2005). As we demonstrated in the frozen-bed simulation, the feature at 2300–2900 m requires an apparent reflectivity gain of 22.8 dB over the surrounding substrate. Roughness may contribute some of this difference; however, the observation is most consistent with a significant material transition. Figure 13 shows that the dielectric contrast between sub-glacial water and most clays or hydrated tills is less than 10 dB, making them unlikely bed materials. Conversely, substrates with the lowest dielectric constants are generally frozen, and they are therefore inconsistent with the presence of liquid water. These constraints limit the likely substrates to a narrow range of materials with ϵ∼ 5, including saturated bedrock or hydrated tills with low permittivity.

Figure 13 Dielectric contrast between liquid water (ϵ=78) and substrates with permittivity from ϵ = 2.7–36. Example bed material permittivity ranges from Tulaczyk and Foley (2020) are provided as reference.

We also assumed constant radar attenuation in the ice. This assumption is reasonable given the short physical distances of our radar simulations. Simulations over greater distances may require introducing variation in the imaginary component of ϵ_ice to account for heterogeneous attenuation loss.

Several notable inconsistencies between our best simulations and the original THW2/UBH0c/X243a data remain. First, R_rel over the beginning and final 500 m of all the simulations are about ∼ 5–8 dB higher than the THW2/UBH0c/X243a RES data. This may indicate a change in σ_bed near the edges of the simulated region. The low-lying area in the region is a perfect topographical feature to accumulate sediments, which likely have low roughness. The steeper and elevated topography near the edges of the simulated region may be exposed bedrock, with greater roughness, which we have shown in our hypothetical case simulations can reduce R_rel. In future iterations of the simulator, we would enable heterogeneity in σ_bed in order to capture this type of variation explicitly.

R_rel near 2000 m along track in all simulations was ∼ 10 dB lower than actual reflectivity from THW2/UBH0c/X243a (Fig. 10). This is coincident with a steep topographical feature in the THW2/UBH0c/X243a radargram. Steep slopes such as this likely represent a limitation of our simulation approach. The surface representation using 5 m flat facets will inherently direct more reflected energy away from the antenna position than a real surface. Therefore, caution is imperative when interpreting results near steep topography.

There are several additional limitations of our simulation method for testing sub-glacial hypotheses. The simulation radius of R = 300 m explicitly limits the impact of range migration and clutter. A more complete simulation incorporating a larger R would be beneficial to simulations in thicker ice or focusing over longer apertures. Future work should combine additional computing power with simulated facet roughness (Gerekos et al., 2023) for additional realism and efficiency.

Our three-dimensional model for topography was derived from a single along-track RES flight line, which enables high confidence and sufficient resolution for ice geometry at nadir. However, this has the obvious limitation of requiring a previous flight line to build our simulation. We lack similar observation density across-track, leaving low-resolution (500 m) open-source DEMs (such as BedMachine V2 (Morlighem, 2020a, b)) as our best option for approximating off-nadir features. This asymmetry is not problematic for replicating basic topography at nadir, but the lack of realistic off-nadir features reduces the sharpness of the focusing algorithm. We can see this effect by comparing the image quality in Fig. 8a to simulated results in Fig. 8b or Fig. 9a. This limitation also clearly reduces our ability to assess hypotheses involving any off-nadir targets.

It is also important to note that the simulation in Fig. 8d exhibits along-track R_rel variation > 30 dB without any change in dielectric properties at the bed. This demonstrates that significant changes in bed echo strength are possible due to topography alone. When inferring the presence of sub-glacial hydrological features from RES data, care must be taken to consider reflectivity within the context of bed topography, hydraulic potential, attenuation, and other factors that may influence the radar echo strength. This observation also demonstrates the value of our simulation methodology for confirming the presence and extent of sub-glacial water.

6 Conclusions

In the exercise presented here, we optimized a radar simulation technique developed by Gerekos et al. (2018) to study the theoretical RES response from sub-glacial systems. The simulator incorporates the Stratton–Chu integral and linear phase approximation to efficiently estimate backscattered radar signal from simulated targets. Through a series of hypothetical simulation cases, we demonstrated the impact on relative reflectivity from rounded Röthlisberger channels or specular flat canals surrounded by bed materials of varying roughness. These simulations confirmed that reflectivity is highly dependent on both the size and cross-sectional shape of the sub-glacial water structure. Our results can be applied broadly to infer the presence, size, and structure of sub-glacial waterbodies from RES surveys in a more robust and sophisticated way than previous methods.

In our example flight line THW2/UBH0c/X243a from Thwaites Glacier, we demonstrated the simulator's utility in testing relevant hypotheses in sub-glacial hydrology. A large water structure could produce the elevated reflectivity beneath Thwaites Glacier, in a region coinciding with a channel route proposed by Hager et al. (2022). Of the scenarios we tested, the radar signature was most consistent with canals averaging > 30 m and covering at least half the area for 600 m along track. Additional simulations including a broader range of hydrological structures and bed materials may be useful to conclusively characterize the feature. Further, although we do not see a Röthlisberger channel at this precise location, our findings do not preclude Röthlisberger channelization further upglacier. Simulations of new and existing RES survey data across the Thwaites catchment could characterize the extent of upglacier channelization.

The method we outline has broad applicability for studying the basal environment of large glaciers. As we have shown, the simulation methodology can offer useful constraints when testing sub-glacial hypotheses. Scientific intuition, additional data inputs, and more computational power will improve the promise of this technique. Future work may include additional computational efficiency to enable more extensive simulations, as well as comprehensive investigations on the impact of substrate roughness, material properties, and additional water geometries.

RES data collection is logistically challenging and expensive. A forward model capable of predicting optimal locations for sub-glacial survey targets, instead of modeling existing flight lines, is an area of interest for future work. In such a forward model, computing resources must be optimized by strategically constraining parameter sets, and performing sensitivity tests for many of the variables considered here.

Appendix A: Along-track focusing

We produce the focused radargram (ξ_f) by processing in along-track blocks. For a given fast-time range bin τ_j, the block size, L_a, is chosen such that range migration for a target equals three fast-time samples, as depicted in Fig. A1. Thus, it is important to note that the block size increases with depth, and the simulation radius R must be greater than the maximum anticipated L_a.

To process a block centered at slow-time a₀, with depth τ_j, we begin with a block of length 2L_a from the range-compressed radargram (ξ_RC) as shown in Fig. A1. We calculated a 1-D reference function (ϕ), representing the Doppler phase modulation as the antenna travels across the aperture in slow-time a (Eq. A1) Peters et al. (2007); Legarsky et al. (2001); Hélière et al. (2007). The amplitude term, b, in Eq. (A1) is used in real-world RES processing to account for along-track variations in instrument gain and aircraft motion and to attenuate high Doppler frequency contributions at long apertures Legarsky et al. (2001); Peters et al. (2007). Our simulations do not contend with non-ideal flight or instrumentation variables, and we therefore use a simple Hamming window of width L_a for suppression of higher-frequency sidelobes.

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_j\left(a\right)=b\left(a\right)e^{-i\mathrm{2}{f}_{\mathrm{c}}\mathit{\tau }\left(a\right)}{|}_{a=a_{\mathrm{0}}-L_{\mathrm{a}}/\mathrm{2}}^{a=a_{\mathrm{0}}+L_{\mathrm{a}}/\mathrm{2}}\text{(A2)}& {\displaystyle }& {\displaystyle }{\mathit{\xi }}_{\mathrm{f}}(a,\mathit{\tau })=\underset{-L_{\mathrm{a}}/\mathrm{2}}{\overset{L_{\mathrm{a}}/\mathrm{2}}{\int }}{\mathit{\xi }}_{\text{RC}}(a+a^{\prime },\mathit{\tau })\mathrm{\Phi }(a^{\prime },\mathit{\tau }{)}^{\ast }\mathrm{d}{a}^{\prime }\end{array}$$

The data block and reference function are Fourier transformed and convolved in the frequency domain (Fig. A2). The result is transformed back to the slow-time domain via inverse Fourier transform. Because the reference function ϕ is tuned at the center of the original block, the middle of the final block is better focused than the edges Hélière et al. (2007). For this reason, only half of the final block (length L_a) is written to the focused radargram ξ_f (Fig. A1). This process is repeated for each along-track block at all fast-time range bins until a complete focused radargram, ξ_f, is formed.

Figure A1 Schematic representation of along-track (azimuth) focusing for a single block at discrete fast-time increment τ_j and azimuth block k in a simulated radargram. The orange hyperbola superimposed on the range-compressed radargram illustrates the theoretical range migration of a target with shortest fast-time range τ_j. We select an aperture length L_a such that range migration spans three sample cells (τ_j, τ_j+1, and τ_j+2).

Appendix B: List of variables used in this paper

Table B1 List of variables used in this paper.