Introduction
Thwaites Glacier, the largest glacier of the Amundsen Sea Embayment in West Antarctica, contains enough ice to raise global sea levels by approximately 60 cm (Holt and others, Reference Holt2006; Scambos and others, Reference Scambos2017; Yu and others, Reference Yu, Rignot, Seroussi and Morlighem2018; Miles and others, Reference Miles, Stokes, Jenkins, Jordan, Jamieson and Gudmundsson2020). Thwaites has exhibited rapidly accelerating ice loss over the past several decades, losing 634 Gt of ice and contributing 1.8 mm of sea level rise between 1979 and 2017 (Mouginot and others, Reference Mouginot, Rignot and Scheuchl2014; Rignot and others, Reference Rignot, Mouginot, Scheuchl, Van Den Broeke, Van Wessem and Morlighem2019). Its retrograde bed geometry is recognized as inherently unstable, and recent evidence suggests irreversible collapse could be underway (Joughin and others, Reference Joughin, Smith and Medley2014; Rignot and others, Reference Rignot, Mouginot, Morlighem, Seroussi and Scheuchl2014; Seroussi and others, Reference Seroussi2017). Significant retreat of the Thwaites grounding line could initiate a widespread collapse of the West Antarctic Ice Sheet, triggering rapid sea level rise on the order of several meters (Scambos and others, Reference Scambos2017).
The subglacial environmental characteristics near the current grounding line will exert direct control on Thwaites’ future. Substrate materials such as sediments, clays, till and bedrock all impact basal friction (Rippin and others, Reference Rippin, Vaughan and Corr2011; Kazmierczak and others, Reference Kazmierczak, Gregov, Coulon and Pattyn2024). At steady-state, deformable bed materials (e.g. hydrated tills, clays and sediments) and inefficient subglacial hydrological structures (e.g. lakes, films and cavities) reduce effective bed pressure. This will reduce basal shear stress and enable ice acceleration (Creyts and Schoof, Reference Creyts and Schoof2009; Hoffman and others, Reference Hoffman2016; Brinkerhoff and others, Reference Brinkerhoff, Aschwanden and Fahnestock2021; Gilbert and others, Reference Gilbert, Gimbert, Thøgersen, Schuler and Kääb2022; Kazmierczak and others, Reference Kazmierczak, Gregov, Coulon and Pattyn2024). Conversely, efficient structures such as Röthlisberger or Nye channels are associated with rigid bedrock. These structures tend to allow for high effective bed pressures and have little impact on basal friction (Röthlisberger, Reference Röthlisberger1972; Hewitt, Reference Hewitt2011; Schroeder and others, Reference Schroeder, Blankenship and Young2013). As the name suggests, efficient basal water structures can transport large volumes of water toward the grounding line, where subglacial discharge enhances mixing and accelerates basal melting (Walder and Fowler, Reference Walder and Fowler1994; Schroeder and others, Reference Schroeder, Blankenship and Young2013; Young and others, Reference Young, Schroeder, Blankenship, Kempf and Quartini2016; Gourmelen and others, Reference Gourmelen2025).
Kazmierczak and others (Reference Kazmierczak, Gregov, Coulon and Pattyn2024) recently projected the impact of these structures’ distribution on Thwaites’ retreat. They demonstrated that specific spatial variations in bed materials (hard vs soft bed rheology) and hydrological structures (channels vs canals) could have a dramatic effect on the extent of retreat in the next century. Projected ice loss varied widely depending on the particular patterns of bed rheology and hydrological structures. These findings are consistent with other research, which routinely cites ambiguity in bed conditions as a primary uncertainty for predicting Thwaites’ future (Scambos and others, Reference Scambos2017).
Multiple studies have attempted to map these spatial variations of subglacial features across Thwaites Glacier using a variety of methods. Jordan and others (Reference Jordan, Thompson, Kulessa and Ferraccioli2023); Borthwick and others (Reference Borthwick2025) inferred several soft sedimentary basins upglacier, with predominantly igneous bedrock near the grounding line, using gravity, magnetic and seismic observations. Others have used high-resolution swath radar and active seismic techniques to distinguish between heterogeneities in bedrock formations and sediment deposition (Muto and others, Reference Muto2019b; Clyne and others, Reference Clyne, Anandakrishnan, Muto, Alley and Voigt2020; Holschuh and others, Reference Holschuh, Christianson, Paden, Alley and Anandakrishnan2020; Alley and others, Reference Alley2021; Borthwick and others, Reference Borthwick2025). Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017) identified a series of connected active subglacial lakes approximately 70–170 km from the grounding line, which were interpreted from satellite-based radar altimetry data. A second drainage cycle from these active lakes in 2017 was linked to temporary ice acceleration (Hoffman and others, Reference Hoffman, Christianson, Shapero, Smith and Joughin2020; Malczyk and others, Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020). Recently, Chartrand and others (Reference Chartrand, Howat, Joughin and Smith2024) also documented enhanced retreat in locations where ice shelf channels intersect the current grounding zone, also utilizing satellite altimetry data.
Airborne radar-echo sounding (RES) is another common technique that has been used extensively to characterize the Thwaites subglacial environment. In hydrologic investigations using RES, the high dielectric contrast between glacier ice, liquid water and other substrate materials (e.g. bedrock, clay and frozen tills) is exploited. Observations of elevated RES reflectivity are commonly interpreted as information about the subglacial hydrological network (Carter and others, Reference Carter, Blankenship, Peters, Young, Holt and Morse2007; Chu and others, Reference Chu, Schroeder, Seroussi, Creyts, Palmer and Bell2016; Rutishauser and others, Reference Rutishauser2018). Similarly, high specularity content—a measure of the mirror-like character of radar reflections—provides an additional metric for inferring the distribution of basal water (Schroeder and others, Reference Schroeder, Blankenship, Raney and Grima2015). However, these interpretations are subject to inherent ambiguities, which are discussed later in this paper.
In 2004–05, the University of Texas Institute for Geophysics (UTIG) conducted the Airborne Geophysical Survey of the Amundsen Embayment (AGASEA), providing RES observations throughout the entire Thwaites catchment. From these data, Schroeder and others (Reference Schroeder, Blankenship and Young2013) identified a pattern of elevated reflectivity and specularity content upstream from a bedrock ridge feature 45–70 km from the grounding line. Downstream of this ridge, radar returns exhibited more asymmetric specularity content. Schroeder and others (Reference Schroeder, Blankenship and Young2013) hypothesized that this radiometric change indicated a hydrological transition beneath Thwaites Glacier, with distributed water structures dominating the upstream hydrology, and largely channelized features below the bedrock ridge. Subsequent studies correlated Schroeder and others (Reference Schroeder, Blankenship and Young2013)’s reflectivity and specularity observations with basal shear stress variations (Haris and others, Reference Haris, Chu and Robel2024) and modeled likely persistent channelized drainage routes through the lower Thwaites region (Fig. 1) (Hager and others, Reference Hager, Hoffman, Price and Schroeder2022).
The downstream region of Thwaites Glacier, with contours of hydraulic potential mapped in brown. Contours intervals are 50 m
$_{H_2O}$, with index lines every 200 m
$_{H_2O}$. Locations of all UTIG survey lines from 2020 and 2022 are shown, with a colormap indicating specularity content. The 30 simulated flight segments in this study are outlined in grey. Segments A–A
$^{\prime}$, B–B
$^{\prime}$, C–C
$^{\prime}$ and D–D
$^{\prime}$ are described in greater detail in the results section. Channel routes from Hager and others (Reference Hager, Hoffman, Price and Schroeder2022); Chartrand and others (Reference Chartrand, Howat, Joughin and Smith2024), subglacial lake locations from Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017) and hydrological transition from Schroeder and others (Reference Schroeder, Blankenship and Young2013) are also shown. The estimated grounding line location from Bindschadler and others (Reference Bindschadler2011) is shown for reference, and light blue shading indicates the Thwaites ice shelf. The East/West dividing line at 
$106.5^{\circ}$ longitude is discussed in the results as a boundary between eastern and western Thwaites.
While RES has provided valuable information about Thwaites’ subglacial environment, it is challenging to verify its conclusions. The most reliable method remains direct borehole access, but complex and expensive logistics mean relatively few subglacial water bodies are verified this way (e.g. subglacial lakes Vostok, Whillans and Mercer) (Engelhardt and Kamb, Reference Engelhardt and Kamb1997; Tulaczyk and others, Reference Tulaczyk, Kamb, Scherer and Engelhardt1998; Jouzel and others, Reference Jouzel1999; Christner and others, Reference Christner2014; Priscu and others, Reference Priscu2021). Multiple ambiguities can obfuscate the basal reflectivity analysis of RES data. Surface roughness can increase scattering losses (Grima and others, Reference Grima, Schroeder, Blankenship and Young2014b; Schroeder and others, Reference Schroeder, Grima and Blankenship2016a), although this surface effect is often minor compared to uncertainty in englacial attenuation (Schroeder and others, Reference Schroeder, Grima and Blankenship2016a; Hills and others, Reference Hills, Siegfried and Schroeder2024; Killingbeck and others, Reference Killingbeck2024), and therefore commonly ignored. Additional complications at the basal interface include the dielectric properties of the substrate material (Tulaczyk and Foley, Reference Tulaczyk and Foley2020) and volumetric scattering associated with water or sediment entrained in basal ice (Hills and others, Reference Hills, Siegfried and Schroeder2024). Pierce and others (Reference Pierce, Skidmore, Beem, Blankenship, Adams and Gerekos2024b) also demonstrated that basal topography can introduce large variations in reflected power. As such, RES-based validation of basal conditions will benefit most from isolating heterogeneities in the basal environment vs englacial or topographic effects.
The simulation method first demonstrated in Pierce and others (Reference Pierce2024a) may prove useful in addressing some concerns around RES validation. Radar returns from a three-dimensional ice/substrate domain are modeled, and simulated results are compared to actual RES reflectivity data. The simulator independently accounts for attenuation, subglacial topography and material properties, making it possible to isolate geometric, material and attenuation impacts on the radar signal.
In this study, we enhance our understanding of Thwaites’ bed heterogeneities in the context of the hypothesized Schroeder and others (Reference Schroeder, Blankenship and Young2013) drainage network. New RES data, collected in 2020 and 2022 with a helicopter-based platform, provided dense coverage of the downstream Thwaites region, including locations coincident with the proposed subglacial water structures in Schroeder and others (Reference Schroeder, Blankenship and Young2013); Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017) and Hager and others (Reference Hager, Hoffman, Price and Schroeder2022) (Fig. 1). We expand the use of the Pierce and others (Reference Pierce2024a) modeling approach to simulate more than 400 km of survey transects. Those simulated outputs are compared to actual reflectivity from the 2020 and 2022 UTIG surveys. Using this tool, we can assess the relative homogeneity of the glacier bed and identify possible locations for widespread subglacial water within the Thwaites downstream environment. We analyze our findings geospatially within the context of previous Thwaites hydrological studies and additional evidence in the 2020 and 2022 RES data.
Methodology
RES data collection and processing
The 2020 and 2022 Thwaites RES surveys were conducted using UTIG’s 60 MHz Multifrequency Airborne Radar Sounder with Full-phase Assessment (MARFA) instrument on an AS-350 B2 helicopter platform (Lindzey and others, Reference Lindzey2020). The airframe and instrumentation were transported to the Thwaites region on the South Korean icebreaking research vessel, RV Araon, where survey flights were launched from the ship’s deck or a field camp on the Thwaites ice shelf. The helicopter was flown at a nominal height of 500 m above the ice surface with a target velocity of 
$30\,\text{m s}^{-1}$. A Renishaw laser altimeter tracked aircraft height above ground level, while a Trimble Net-R9 dual frequency Global Navigation Satellite System and Novatel SPAN IGM-1A recorded precise aircraft positioning and orientation, as described in Lindzey and others (Reference Lindzey2020).
Data were range compressed, corrected for geometric spreading loss and aircraft position, and focused in azimuth (Peters and others, Reference Peters, Blankenship, Carter, Kempf, Young and Holt2007). We estimated englacial attenuation separately for the 2020 and 2022 survey data, using a linear regression model similar to the methods described in Jacobel and others (Reference Jacobel, Welch, Osterhouse, Pettersson and Macgregor2009); Schroeder and others (Reference Schroeder, Blankenship and Young2013; Reference Schroeder, Seroussi, Chu and Young2016b). In our regression model, a linear fit between geometrically corrected bed returned power 
$\left[ P_r^G \right]_{dB}$ and ice thickness 
$d_{ice}$ was identified independently for each field season. The resulting slope is equal to the one-way attenuation rate 
$\left \lt N_{ice}\right \gt $ in 
$\mathrm{dB km}^{-1}$ (Eqn (1)). Bed reflectivity 
$\left[ R_{bed} \right]_{dB}$ and system gain 
$\left[ S \right]_{dB}$ are assumed homogeneous in this method. We ignore other possible sources of power loss (
$\left[L \right]_{dB}$) such as birefringence and volumetric scattering in our model, although we will revisit this assumption in the Discussion. All quantities denoted by 
$\left[ \right]_{dB}$ are expressed on the decibel scale, and 
$\left \lt \right \gt $ indicates a depth-averaged quantity.
\begin{equation}
\left[P_r^G \right]_{dB} = \left[R_{bed} \right]_{dB} - 2\,d_{ice}\,\left \lt N_{ice} \right \gt + \left[S \right]_{dB} + \left[L \right]_{dB}
\end{equation} The resulting one-way attenuation rates were 14.5 and 
$14.0\,\text{dB km}^{-1}$ for the 2020 and 2022 survey data, respectively. The 2020 root-mean-squared error (RMSE) between corrected bed power 
$\left[ P_r^G \right]_{dB}$ and the regression fit is 9.7 dB, while for the 2022 survey, the RMSE was 6.4 dB. These values are in line with Schroeder and others (Reference Schroeder, Seroussi, Chu and Young2016b), which identified localized attenuation rates between 10 and 16
$\,\text{dB km}^{-1}$ within 100 km of the grounding line.
Basal returns are identified in the radar data through manually determined limits (one above and one below), encapsulating the identified basal reflector. The magnitude of the basal return is taken as the sample with the maximum returned power between those bounds, at each vertical trace along-track. This value is modified as described above, through geometric and attenuation corrections, to determine the corrected bed reflectivity 
$R_{bed}$ (Blankenship and others, Reference Blankenship2001). We note that this value does not include any correction for volumetric scattering of basal ice, which can impact apparent reflectivity (Hills and others, Reference Hills, Siegfried and Schroeder2024).
We also consider the specularity content of the basal reflection, as described in Schroeder and others (Reference Schroeder, Blankenship, Raney and Grima2015). Specularity content measures the diffuse or specular nature of radar reflections by comparing the ratio of returned bed power over two different focusing apertures (nominally 700 m and 2 km, Schroeder and others (Reference Schroeder, Blankenship, Raney and Grima2015)). Specularity content close to 0 implies that reflected energy is scattered diffusely, as may be expected with a rough bed surface. By contrast, specularity content closer to 1 implies a smooth surface relative to the radar wavelength 
$\lambda$, consistent with a range of distributed subglacial water structures.
RES simulations
The radar simulation method described in Pierce and others (Reference Pierce2024a) estimates variations in returned power from the subglacial interface by tracing the path and field strength of many individual rays transmitted through the ice surface then reflected from the bed. The simulated ice geometry mimics observations from the RES flight path. Bed roughness (
$\sigma_{bed}$) and dielectric permittivity for ice (
$\epsilon_{ice}$), liquid water (
$\epsilon_{H_2O}$) and primary substrate material (
$\epsilon_{bed}$) are defined by the user. A conceptual model for the process, along with simulation parameters used in this study, is described in the Supplemental Materials and in Pierce and others (Reference Pierce2024a).
We modeled 30 segments of RES transects from the 2020 and 2022 Thwaites surveys (Fig. 1). Each simulated segment was between 10 and 15 km, totaling 408 km. Many segments were chosen for their proximity to hypothesized hydrological features (Smith and others, Reference Smith, Gourmelen, Huth and Joughin2017; Hager and others, Reference Hager, Hoffman, Price and Schroeder2022), while providing sampling across the entire survey footprint. This effort encompasses over 
$20\%$ of the total RES data collected over grounded ice, and is a 100-fold increase in the geographic extent simulated in Pierce and others (Reference Pierce2024a).
For each flight segment, we first compared the actual variation in along-track reflectivity to a simulation with bed roughness 
$\sigma_{bed} = 0.2$ m, substrate permittivity 
$\epsilon_{sub} = 5$, and no liquid water. This scenario, henceforth referred to as the ‘homogeneous’ case, serves as the experimental control. This null hypothesis helps to diagnose whether variations in bed returned power could result from bed topography alone, without other material or geometric variations in the subglacial environment.
We can deviate from the homogeneous case by manipulating two primary independent variables over part of the segment. Introducing spatial heterogeneity in substrate permittivity and/or small-scale roughness will influence the resulting along-track simulated reflectivity:
• Substrate permittivity—changing material properties along the flight line will invoke a response in reflected power, because reflected power depends on the dielectric contrast between ice (
$\epsilon_{ice}=3.18$) and the substrate (Peters and others, Reference Peters, Blankenship and Morse2005; Tulaczyk and Foley, Reference Tulaczyk and Foley2020). In this study, we often model transitions between bedrock (
$\epsilon_{sub}=5$) and liquid water (
$\epsilon_{H_2O} = 78$), in locations where the observed RES reflectivity increases by 10–15 dB over a broad area (
$ \gt 500$ m along-track), but the homogeneous simulation does not exhibit the same reflectivity response. Additional transitions (such as from bedrock to clay or till with
$\epsilon_{sub}=10$) are applied in this study as appropriate to fit the RES observations.• Small-scale roughness—we impose a random Gaussian roughness
$\sigma_{bed}$ over correlation length
$l_c=15$ m. Diffuse radar scattering due to roughness is highest over correlation lengths near the radar wavelength (
$\lambda =5$ m) (Haynes and others, Reference Haynes, Chapin and Schroeder2018; Gerekos and others, Reference Gerekos, Haynes, Schroeder and Blankenship2023). Therefore,
$l_c=15$ m was chosen as the minimum distance over which we could introduce roughness for a DEM with 5 m resolution. Realistic estimates for
$\sigma_{bed}$ over this length scale are somewhat speculative. Many studies have attempted to capture basal roughness magnitude beneath Antarctic ice sheets, but current methods limit analysis to km-scale correlation lengths (Taylor and others, Reference Taylor, Siegert, Payne and Hubbard2004; Bingham and Siegert, Reference Bingham and Siegert2009; Rippin and others, Reference Rippin, Vaughan and Corr2011) or evaluate small-scale roughness over a very limited geographic extent (Hoffman and others, Reference Hoffman, Christianson, Holschuh, Case, Kingslake and Arthern2022). Without strong observational data, we chose roughness magnitude
$\sigma_{bed}=0.2$ m for most subglacial DEMs in this study, as in Pierce and others (Reference Pierce2024a). For locations where we seek to model a rougher bed,
$\sigma_{bed}$ was increased to 1 m. Further analysis of these constraints will be offered in the Discussion section.
It is important to note that our segment-based approach introduces an important constraint when evaluating fit to observed RES reflectivity. A short segment relative to the along-track resolution of the radar data will always appear homogeneous, whereas a very long segment will generally exhibit multiple heterogeneous features. We assume that our simulation approach can successfully evaluate basal features more than an order of magnitude greater than the along-track resolution of the synthetic aperture (SAR) focused RES data (10–15 m) and significantly less than the overall segment length (10–15 km). Therefore, simulated substrate permittivity or small-scale roughness is manipulated only over distances between 300 and 5000 m. This aligns with Muto and others (Reference Muto, Alley, Parizek and Anandakrishnan2019a; Reference Muto2019b), which identified variations in Thwaites’ substrate material and basal shear stress over distances of hundreds of meters to a few km.
Evaluating simulation fit quality
Simulations were performed in 10–15 km segments along the flight paths depicted in Fig. 1. This simulation length was chosen largely on the basis of computational processing constraints. We compared the simulated relative variation in bed reflectivity to actual data from the 2020 and 2022 Thwaites RES surveys. We assess the fit between the real and simulated datasets using three different metrics:
• Pearson correlation coefficient (
$\Phi_{RS}$) is defined in Eqn (2), where
$\sigma_R$ and
$\sigma_S$ are the standard deviations for the RES and simulated reflectivity data (in dB). The covariance between the real and simulated datasets (
$cov_{R,S}$) is calculated after smoothing each with a Savitzky–Golay filter (720 m window) to reduce high-frequency noise.
$\Phi_{RS}$ ranges between
$-$1 and 1. Positive values indicate positive correlation, with higher magnitudes implying the correlation strength.
(2)\begin{equation}
\Phi_{RS} = \frac{cov_{R,S}}{\sigma_R \: \sigma_S}
\end{equation}• Standard deviation ratio (
$\sigma_{RS}$) is the ratio in standard deviation of along-track reflectivity between the real and simulated data. Values near one indicate a strong match between the magnitude of fluctuations in along-track reflectivity.
(3)\begin{equation}
\sigma_{RS} = \sigma_S / \sigma_R
\end{equation}• Range ratio (
$\chi_{RS}$) is the ratio of the total range in reflectivity (in dB) between the real and simulated data. A value near one indicates that both datasets have similar magnitudes.
(4)\begin{equation}
\chi_{RS} = \frac{\max{\left[R_{bed} \right]_{dB}} - \min{\left[R_{bed} \right]_{dB}}}{\max{\left[R_{bed, S} \right]_{dB}} - \min{\left[R_{bed,S} \right]_{dB}}}
\end{equation}
It is important to note that these metrics differ in their significance to fit quality assessment. Only 
$\Phi_{RS}$ accounts for spatial patterns in the radar return, whereas the other two metrics indicate whether fluctuations have consistent magnitudes. Therefore, 
$\sigma_{RS}$ and 
$\chi_{RS}$ are only meaningful in simulations where a minimum 
$\Phi_{RS}$ threshold is met.
Results
RES specularity content
The 2020 and 2022 RES surveys included nearly 1800 km of radar transects over Thwaites’ grounded ice. Nearly all of these data were located downglacier from the bedrock ridge identified in Schroeder and others (Reference Schroeder, Blankenship and Young2013). The data exhibit some patterns in basal reflectivity and specularity content similar to the 2004–05 UTIG RES survey. We are primarily interested in the subglacial hydrological system, and therefore RES observations over the Thwaites ice shelf are excluded from our analysis (Fig. 1). Further, we divide the data along 
$106.5^{\circ}$ W longitude when comparing observations between eastern and western Thwaites.
Observations of asymmetry in the specularity content are noteworthy for this downstream region of Thwaites Glacier, as this formed the basis for the channelized drainage hypothesis outlined in Schroeder and others (Reference Schroeder, Blankenship and Young2013). Table 1 summarizes key findings in the specularity content data from the 2020/2022 Thwaites survey. Mean specularity content across all included data from these surveys was 
$0.19\pm0.18$. East of the dividing line, mean specularity content was more than 
$50\%$ higher than the value for western Thwaites (
$0.23$ vs 
$0.15$, Table 1). Also similar to observations from Schroeder and others (Reference Schroeder, Blankenship and Young2013) in the 2004–05 AGASEA survey, radar transects oriented north/south have a higher mean specularity than transects oriented east/west. However, this orientation dependence is primarily observed east of the dividing line, where mean specularity content for north/south transects is 
$0.24\pm0.20$ and east/west transects is 
$0.16$. West of the dividing line, specularity content of north/south and east/west oriented flight paths is statistically equivalent. Despite relatively high standard deviations, two-tailed 
$p$-values for each of these observations indicate greater than 
$99\%$ confidence that specularity content patterns are correlated across the east/west divide (Box and others, Reference Box, Hunter and Hunter1978).
Specularity content for various subsets of the 2020/2022 Thwaites survey data.
SC = specularity content.
Comparing RES and simulated reflectivity
The simulation method assumes constant englacial attenuation and power losses at the ice surface are governed by the DEM, which mimics the real ice surface. Eliminating these common uncertainties ensures that variations in simulated bed reflectivity are the result of basal topography and substrate materials. Comparisons between simulations and actual RES reflectivity are, therefore, helpful in diagnosing subglacial environmental conditions.
In the field, calibrating the returned power of the ice-penetrating RES to absolute reflectivity is often infeasible. Therefore, we consider relative reflectivity along each simulated segment by centering the datasets about the mean value. This approach enables us to compare the magnitude of reflection variations along the segment without requiring calibration.
We introduce a simple framework for characterizing the fit quality of our simulations to the actual RES along-track reflectivity:
(1) Homogeneous fit: 12 out of our 30 simulated flight segments (40%) exhibited similarity to the real RES data across all three fit quality metrics, with positive correlation
$\Phi_{RS} \gt 0.2$, standard deviation ratio
$\sigma_{RS} \gt 0.65$ and range ratio
$\chi_{RS} \gt 0.65$. We interpret these locations as having relatively homogeneous bed material and subglacial roughness along-track.(2) Hydrology fit: for nine of the remaining segments (30%), the fit could be significantly improved by introducing a water feature extending across the simulated bed DEM. In these scenarios, parameters for the remainder of the bed DEM remain unchanged from the homogeneous simulation.
(3) Substrate transition: for three segments (10%), fit quality improvement was achieved by introducing a coupled change in both
$\sigma_{bed}$ and
$\epsilon_{sub}$ over an area extending across the DEM. Such a coupled change could represent a material transition from a deformable substrate such as clay or till to bedrock.(4) Heterogeneous: for the remaining six segments (20%), the homogeneous simulation demonstrated a poor fit to the RES data. Improving the fit quality would require complex manipulations of the bed DEM variables in multiple locations along the segment.
Figure 2 displays the location of all 30 simulated segments, with color coding indicating the bed diagnostic taxonomy described above. We will evaluate representative case studies from each category to build intuition for the above framework. This should establish an understanding of the simulation method’s potential and limitations.
Fit quality for all simulated flight segments. Hydraulic potential contours shown in brown with increments of 50 m
$_{H_2O}$ (calculated using Goff bed map (Goff and others, Reference Goff, Powell, Young and Blankenship2014) and REMA surface data (Howat and others, Reference Howat, Porter, Smith, Noh and Morin2019)). Grayscale indicates MEaSURES ice speed, with dark gray 
$\sim 0$ m yr
$^{-1}$ and white 
$ \gt 1\,\mathrm{km yr}^{-1}$ (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017; Matsuoka and others, Reference Matsuoka, Skoglund and Roth2018). Locations where water or material transition was modeled to generate a fit are shown with darker shading for hydrology and substrate transition fits. Insets 1 and 2 are provided for clarity.
Homogeneous fits
Figure 3a shows the focused radargram for a 14 km segment between A and A 
$^{\prime}$ in Figs. 1 and 2, with the along-track distance measured from the center point of the simulated segment. Between 
$-$650 and 1250 m along-track, we observe a 15 dB reflectivity peak (Fig. 3c). This bright reflection is coincident with the location of a hypothesized active subglacial lake from Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017). From the RES surface elevation (
$z_s$) and ice thickness (
$d_{ice}$) measurements, we can estimate the hydraulic potential along-track (Eqn (5)), assuming zero effective pressure at the bed (Brocq and others, Reference Brocq, Payne, Siegert and Alley2009). In Eqn (5), 
$\rho_{ice}$ is ice density, and 
$\rho_{H_2O}$ is the density of water. Water will flow along hydraulic potential gradients and pool in locations that are hydraulically flat. Hydraulic potential for the same flight segment is plotted in Fig. 3b, with color corresponding to variation in 
$R_{bed}$. This shows the bright reflector is also hydraulically flat, and coincides with a distinct increase in specularity content (Fig. 3b)
\begin{equation}
\psi = d_{ice} \left(\frac{\rho_{ice}}{\rho_{H_2O}} - 1 \right) + z_s
\end{equation}(a) Focused radargram for the segment between A and A 
$^{\prime}$. The bright reflector is boxed. (b) Graph of hydraulic potential for the same simulated segment, with color indicating 
$R_{bed}$ variation. Specularity content is shown in gray. (c) Along-track reflectivity variation 
$R_{bed}$ and simulated homogeneous 
$R_{bed,S}$. Locations A and A 
$^{\prime}$ correspond to the locations in Figs. 1 and 2. The hypothesized lake boundaries from Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017) intersect this segment from 0 to 7000 m along-track.
Figure 3c compares the along-track variation in 
$R_{bed}$ vs 
$R_{bed,S}$ for the homogeneous simulation of the same 
$14km$ segment. We can visually confirm the correlation of peaks and valleys between the two datasets, resulting in 
$\Phi_{RS} = 0.78$. Both datasets exhibit similar standard deviations and peak heights, as indicated by 
$\sigma_{RS}=1.18$ and 
$\chi_{RS}=1.17$, respectively. Although the bed is hydraulically flat and radar-bright, the strong fit with the homogeneous simulation indicates the 15 dB increase in along-track 
$R_{bed}$ may result entirely from topography. A material change with a large dielectric contrast (such as a transition from rock to water) is not required to induce the reflectivity increase between 
$-$650 and 1250 m.
Table 2 summarizes the fit quality metrics for all simulated segments identified as homogeneous.
Hydrologic fits
Eighteen of our 30 simulated flights do not exhibit a strong homogeneous fit as described above. Of these simulated segments, the fit quality for nine can be improved by changing the permittivity over part of the bed DEM from the default value 
$\epsilon_{sub}=5$ to 
$\epsilon_{H2O}=78$. We impose this change over an area at least 300 m along-track and extending across the simulated bed DEM, which could represent an extensive area of distributed water (e.g. a subglacial lake or seawater intrusion). Our modeling choice is justified by the geographic location or alternative evidence in the radar data, such as increased specularity content.
Figure 4a shows the radargram for an 11 km segment located from B to B 
$^{\prime}$ in Figs. 1 and 2. This segment intersects the grounding line near several of the channel locations referenced in Hager and others (Reference Hager, Hoffman, Price and Schroeder2022); Chartrand and others (Reference Chartrand, Howat, Joughin and Smith2024). An abrupt reflectivity transition of nearly 20 dB occurs around 
$-$1800 m along-track. Figure 4b compares the along-track variation in 
$R_{bed}$ to the homogeneous simulation case. Fit quality metrics 
$\Phi_{RS} = 0.26$, 
$\sigma_{RS} = 0.54$ and 
$\chi_{RS} = 0.50$ imply the bed has lower homogeneity than the previous example from A to A 
$^{\prime}$.
(a) 1-D focused radargram for the 11 km segment from B to B 
$^{\prime}$ (see Figs. 1 and 2). The bright basal reflector is highlighted. (b) Comparison of 
$R_{bed}$ and 
$R_{bed,S}$ for the homogeneous simulation. (c) Hydrologic simulation comparison of 
$R_{bed}$ vs 
$R_{bed,S}$, where the basal material is changed to liquid water (
$\epsilon=78$) from 
$-$4100 to 
$-$1800 m along-track.
The 20 dB change in reflectivity is approximately coincident with the grounding line, and therefore, we consider this change in reflectivity could mark the transition between grounded and floating ice. In Fig. 4c, we simulate the same segment with a substrate permittivity change from 
$\epsilon_{rock}=5$ to 
$\epsilon_{H_2O}=78$ at 
$-$1800 m along-track. The material transition induces a strong reflectivity response, and fit quality improves to 
$\Phi_{RS} = 0.91$, 
$\sigma_{RS} = 1.11$ and 
$\chi_{RS} = 0.90$. We therefore consider this bed material model as a plausible interpretation of the RES data.
This pattern is repeatable for four additional flight segments near the Thwaites grounding zone (locations 15, 16, 17 and 19 in Table 3 and Fig. 2). Four others occur further inland (locations 13, 14, 20 and 21, Table 3 and Fig. 2), where the modeled water feature is consistent with a peak in specularity content. For each of these simulated segments, fit improvements are outlined in Table 3. Along-track comparisons of 
$R_{bed}$ vs 
$R_{bed,S}$ can be found in the Supplemental Materials.
Substrate transitions
Three additional simulated segment locations from the 2020 and 2022 RES survey exhibit 
$R_{bed}$ patterns consistent with a simple substrate transition over a single area. A 14 km segment from C to C 
$^{\prime}$ (Figs. 1 and 2), on the western edge of the lower Thwaites catchment, demonstrates this type of heterogeneity. There are no obvious bright reflectors in the radargram (Fig. 5a), though 
$R_{bed}$ spans a range of 25 dB over the segment. The homogeneous simulation (Fig. 5c) has fit metrics 
$\Phi_{RS} = 0.15$, 
$\sigma_{RS} = 0.57$ and 
$\chi_{RS} = 0.67$, indicating the homogeneous model does not fit the RES data well. However, the fit is binary. The section from 
$-$4000 to 
$-$1500 m demonstrates poor coherence with the actual variation in 
$R_{bed}$, while the rest of the homogeneous simulation appears to correlate.
(a) 1-D focused radargram for the 14 km segment from C to C 
$^{\prime}$ (Figs. 1 and 2). (b) Specularity content from the RES data for the simulated segment. (c) Comparison of along-track 
$R_{bed}$ and 
$R_{bed,S}$ for the homogeneous scenario. (d) 
$R_{bed}$ vs 
$R_{bed,S}$ comparison for a roughness only transition (constant 
$\epsilon_{sub}$). (e) 
$R_{bed}$ vs 
$R_{bed,S}$ comparison for a substrate material transition (coupled 
$\epsilon_{sub}$ and roughness change).
This shift in fit quality is paired with an apparent reduction in specularity content to near zero (Fig. 5b). Low specularity indicates that incident radar energy is scattered more diffusely, and therefore may imply a change in basal roughness (Schroeder and others, Reference Schroeder, Blankenship, Raney and Grima2015). When the basal roughness 
$\sigma_{bed}$ is increased to 1 m over the segment from 
$-$4000 to 
$-$1500 m along-track (Fig. 5d), we see improved correlation (
$\Phi_{RS} = 0.68$). Standard deviation ratio and range ratio agreement also increase, although 
$\sigma_{RS} = 0.61$ remains below our fit quality threshold.
The results from Fig. 5d indicate that a change in roughness alone may not explain the along-track variation in 
$R_{bed}$ from C to C 
$^{\prime}$. We could reasonably speculate that the change in roughness indicates a material transition, with a corresponding change in the bed permittivity. Figure 5e shows results from such a coupled scenario, where the lower roughness value 
$\sigma_{bed}=0.2$ m is paired with 
$\epsilon_{bed}=10$ (e.g. clay bearing or sandy till), and higher roughness 
$\sigma_{bed}=1.0$ m is associated with 
$\epsilon_{bed}=5$ (e.g. rough bedrock) (Tulaczyk and Foley, Reference Tulaczyk and Foley2020). Fit quality metrics 
$\Phi_{RS} = 0.65$, 
$\sigma_{RS} = 0.91$ and 
$\chi_{RS} = 0.98$ all imply this is a possible interpretation of the subglacial environment consistent with the RES data at this location.
Two additional simulations exhibit improvement in fit quality if we assume similar coupling of bed material and roughness (segments 21, 22 and 23 in Fig. 2). Fit improvements for these simulations are summarized in Table 3.
Heterogeneous
For the remaining six simulated segments (20%), the homogeneous simulation scenario exhibits poor fit quality to the RES data. Achieving improved fit requires manipulation of more than one simulation variable without logical coupling, or abrupt material changes in multiple locations along the bed DEM. Our methodology allows such adjustments within the simulation environment, and these complex manipulations may achieve a high-quality fit to the RES reflectivity data. However, without consistent corroborating evidence for multiple material changes at the bed (such as specularity content or geographic location near the grounding line), fitting these data would simply optimize our two primary model variables independently. The resulting solutions could be non-unique or physically unrealistic. Therefore, instead of attempting to diagnose specific spatial patterns in bed heterogeneity to explain the along-track 
$R_{bed}$ variation, we assume that the bed in these locations has more complex heterogeneity than the other simulated flight segments.
The 14 km segment from D to D 
$^{\prime}$ (Figs. 1 and 2) is an example of heterogeneity. Figure 6 shows homogeneous fit quality metrics of 
$\Phi_{RS} = 0.10$, 
$\sigma_{RS}=0.48$ and 
$\chi_{RS}=0.62$, confirming a relatively poor fit between 
$R_{bed}$ and 
$R_{bed,S}$. This segment exhibits a range of nearly 
$30\,dB$ in 
$R_{bed}$ along-track, with at least two bright reflectors centered at 
$-$4.4 and 4.5 km. Adjacent to each of these reflectors are areas of dim basal reflectivity, where 
$R_{bed}$ is reduced by nearly 30 dB. In the center of the segment between 
$-$3.5 to 1.0 km along-track, variation in both 
$R_{bed}$ and 
$R_{bed,S}$ is of lower magnitude, but there is limited correlation in trends between the two series.
None of the large fluctuations in 
$R_{bed}$ are present in the homogeneous simulation scenario for the D to D 
$^{\prime}$ segment. This result implies that the bright basal reflections more likely arise from basal heterogeneity than topographic variations. Improving the fit quality between 
$R_{bed}$ and the simulated 
$R_{bed,S}$ could include the addition of water structures near the two bright reflectors, as well as adding roughness or reducing 
$\epsilon_{sub}$ in locations with dim 
$R_{bed}$ in the RES data. All of these simultaneous changes over a short distance would be speculative. Instead of attempting to fit the RES data, we diagnose this location as having complexity in the basal environment, which the homogeneous model does not capture.
Table 4 summarizes the fit quality metrics for all six simulated segments identified as heterogeneous.
Discussion
Geographic distribution of results
The simulated segments were located across the Thwaites lower catchment, with 13 (43%) east of the dividing line at 
$106.5^{\circ}$ W and 17 (57%) to the west. Simulated bed homogeneity or heterogeneity appears to contrast across this east/west division. Of the 12 segments identified as homogeneous fits, 10 were to the west. Conversely, five out of six heterogeneous segments were to the east (Fig. 2) and at least 10 km from the grounding line. We hypothesize that the disparity across lower Thwaites could be indicative of an east/west transition in the substrate materials and hydrological system.
The distribution of hydrologic fits is of particular interest in this study. Hydrologic fits at locations 15 through 21 were within a few km of the grounding zone, where seawater infiltration and grounding line retreat likely form large basal water bodies. In addition, segments 15 and 19 intersect ice shelf channels identified in Chartrand and others (Reference Chartrand, Howat, Joughin and Smith2024), and 18, 20 and 21 intersect channel locations from Hager and others (Reference Hager, Hoffman, Price and Schroeder2022).
Segments 13 and 14 are in the northeast corner of our study region (Fig. 2). These locations are the furthest from the grounding zone, lie in a location with minimal hydraulic gradient, and may be upglacier of a hypothesized hydrologic transition zone. Schroeder and others (Reference Schroeder, Blankenship and Young2013) proposed that inefficient drainage, including widespread distributed canals, dominates the upglacier hydrology. Our observations of hydrologic simulation fits to the RES reflectivity, combined with high specularity content (Fig. 1) at locations 13 and 14, are consistent with this hypothesis.
Many of the segments in western Thwaites identified as homogeneous, hydrologic, or material transition fits were beneath ice with velocity exceeding 1 km yr
$^{-1}$ (Fig. 2), and immediately downstream of the 1.7 km deep sedimentary basin identified in Borthwick and others (Reference Borthwick2025). The Borthwick and others (Reference Borthwick2025) study suggested that a deformable sediment layer extends at least to the southern limit of segment 9. Given that deformable substrate material and widespread, distributed hydrologic systems are likely beneath the fast-flowing ice of western Thwaites (Peters and others, Reference Peters2006), our results imply that this deformable sediment layer could extend to within a few km of the grounding line. MacGregor and others (Reference MacGregor2013) indicated that ice flow across western Thwaites may be controlled by topographic features, including prominent volcanoes along the western shear margin. Our simulations are also consistent with this interpretation, where model fits to bed reflectivity align with the topography and basic substrate transitions employed in our RES simulation methods.
The heterogeneity inferred across much of eastern Thwaites could be explained by the opposite interpretation (Fig. 2). Thwaites is believed to have fewer prominent topographic features regulating ice flow on its eastern side (MacGregor and others, Reference MacGregor2013), downglacier from the Schroeder and others (Reference Schroeder, Blankenship and Young2013) hydrologic transition. Schroeder and others (Reference Schroeder, Grima and Blankenship2016a) observed a sharp contrast in bed reflectivity near the eastern Thwaites shear margin, which they inferred as a change in conditions as a primary control. Previous seismic studies far upglacier of our RES survey have demonstrated a correlation between heterogeneous basal sediment distribution and topographic undulations on the order of 100 m (Muto and others, Reference Muto2019b; Clyne and others, Reference Clyne, Anandakrishnan, Muto, Alley and Voigt2020). Although we did not simulate such a scenario in this study, the heterogeneity deduced from our simulations may be consistent with a basal environment varying between bedrock, till, dilated sediments and water among rolling topography. Concentrated channels, which exhibit less pronounced but broader reflectivity peaks (Pierce and others, Reference Pierce2024a), could also contribute to 
$R_{bed}$ variations that do not align with topography across the area. Additional RES modeling and seismic studies in eastern Thwaites, similar to those found in Muto and others (Reference Muto2019b); Clyne and others (Reference Clyne, Anandakrishnan, Muto, Alley and Voigt2020); Borthwick and others (Reference Borthwick2025), would provide additional evidence for this hypothesis.
The spatial patterns in model agreement with the RES reflectivity data correlate with asymmetry in the RES specularity content data. We observe significantly higher mean specularity content in the east, as well as asymmetry between north/south and east/west oriented transects (Table 1 and Fig. 1). Similar asymmetry in an extensive 2004–05 survey was the primary evidence for a hypothesized hydrological transition 45–70 km upstream of the Thwaites grounding line (Schroeder and others, Reference Schroeder, Blankenship and Young2013).
West of 
$106.5^{\circ}$ W, our evidence points to broad areas of homogeneity, with distributed water bodies near the grounding zone. We speculate that fast-moving ice on the west side of the glacier significantly alters the local hydrological system through frictional heating and other dynamic processes at the bed surface.
Based on our simulation results and observed specularity content from the 2020 and 2022 surveys, we propose a revised interpretation of the Thwaites downstream hydrological system. Schroeder and others (Reference Schroeder, Blankenship and Young2013) hypothesized a hydrological system dominated by channelized, efficient drainage within 40–70 km of the grounding line. Our interpretation is consistent with this model for eastern Thwaites, however in western Thwaites, our results could imply more inefficient hydrology. As demonstrated by Kazmierczak and others (Reference Kazmierczak, Gregov, Coulon and Pattyn2024), the spatial variability in hydrological characteristics may play a critical role in modulating future ice mass loss across Thwaites.
Diagnostic capabilities of RES simulations
In RES data, it is common to observe variations in 
$R_{bed}$ exceeding 20 dB over a few km along-track. Bright reflectors are commonly interpreted as water structures, especially if they are located in hydraulically flat topography (Peters and others, Reference Peters, Blankenship and Morse2005; Young and others, Reference Young, Schroeder, Blankenship, Kempf and Quartini2016; Rutishauser and others, Reference Rutishauser2018). We have demonstrated through our simulation technique that variations of this magnitude may result from subglacial water (as in B to B 
$^{\prime}$) or topographic features (as in A to A 
$^{\prime}$). This finding is consistent with Pierce and others (Reference Pierce, Skidmore, Beem, Blankenship, Adams and Gerekos2024b), which demonstrated topographically induced 
$R_{bed}$ variations beneath Devon Ice Cap of nearly 30 dB. Inferring the basis for specific 
$R_{bed}$ fluctuations in RES data can be challenging with existing interpretation techniques. The simulation methodology described here may be a practical tool for this diagnosis.
For example, the reflector between A and A 
$^{\prime}$ could be classified as a ‘definite’ or ‘dim’ lake according to the scheme of Carter and others (Reference Carter, Blankenship, Peters, Young, Holt and Morse2007). A system of active lakes connected to this region underwent two discharge cycles between 2013 and 2018 (Smith and others, Reference Smith, Gourmelen, Huth and Joughin2017; Malczyk and others, Reference Malczyk, Gourmelen, Goldberg, Wuite and Nagler2020). However, during the later drainage event, no significant change in ice surface elevation was observed at this location, leaving uncertainty as to whether the segment A to A 
$^{\prime}$ is hydrologically active. In addition, spatial variability in effective bed pressure may be substantial, complicating the assumption that the region is hydraulically flat. Together, these factors introduce uncertainty into the interpretation of A to A 
$^{\prime}$ as a subglacial lake, despite the presence of radar-bright reflections. Our simulations suggest that the elevated radar reflectivity is likely controlled by bed topography alone. This supports the interpretation that no subglacial lake is present despite an elevated 
$R_{bed}$ of 15 dB and a corresponding increase in specularity content.
Independent manipulation of along-track substrate permittivity and roughness in our simulation approach can produce an array of modeling scenarios, with varying fit quality to the RES data. In this study, we have not tested all possible permutations, nor have we introduced a robust inverse framework or process to optimize the simulation/RES data agreement. As such, modeling scenarios demonstrating strong agreement with the RES data cannot be interpreted as conclusive. Rather, they provide one possible interpretation for the observed variations in RES reflectivity warranting further evaluation.
We deliberately confined our modeling variations to the ice–substrate interface, ignoring other known sources of power loss (e.g. 
$\left[L \right]_{dB}$ in Eqn (1)). Factors such as surface scattering and englacial effects including birefringence, attenuation variations and volumetric scattering also influence returned radar power (Matsuoka and others, Reference Matsuoka, MacGregor and Pattyn2012; Grima and others, Reference Grima, Blankenship, Young and Schroeder2014a; Reference Grima, Schroeder, Blankenship and Young2014b; MacGregor and others, Reference MacGregor2015; Schroeder and others, Reference Schroeder, Seroussi, Chu and Young2016b; Hills and others, Reference Hills, Siegfried and Schroeder2024). Attenuation and volumetric scattering losses will generally exceed surface scattering and birefringent effects across a wide range of conditions (Matsuoka and others, Reference Matsuoka, MacGregor and Pattyn2012; Schroeder and others, Reference Schroeder, Seroussi, Chu and Young2016b; Broome and Schroeder, Reference Broome and Schroeder2020), implying these complexities could be significant in locations identified as heterogeneous. However, without a robust optimization framework, adding further model complexity may only reduce diagnostic power. Thus, our heterogeneous segments have identified locations where a simple material model does not agree with observed RES data. We presume these locations require simulations with greater physical complexity, which we leave for future studies.
We also note that our current RES simulation approach does not account for the influence of cross-track topography (clutter) on basal returns—an important limitation and an area of active future interest. High-resolution bed DEMs derived from swath radar (Holschuh and others, Reference Holschuh, Christianson, Paden, Alley and Anandakrishnan2020) offer promise for incorporating these effects, though broader geographic coverage will be needed to support applications like the one presented here.
Sensitivity of RES simulations to model parameters
We performed most homogeneous simulations assuming bed permittivity 
$\epsilon_{sub} = 5$, roughness magnitude 
$\sigma_{bed}=0.2$ m and 
$\epsilon_{ice} = 3.18$. This is a relatively smooth interface with low dielectric contrast, possibly representing bedrock beneath glacier ice (Tulaczyk and Foley, Reference Tulaczyk and Foley2020). Such an interface is reasonable to envision beneath Thwaites Glacier; however, a wide range of alternative values for 
$\epsilon_{sub}$ and 
$\sigma_{bed}$ could exist (Peters and others, Reference Peters, Blankenship and Morse2005; Tulaczyk and Foley, Reference Tulaczyk and Foley2020) with unknown spatial distribution. Homogeneity does not require specific material properties, only that those properties do not change appreciably over the flight segment. Therefore, we assess the sensitivity of the diagnostic method developed here to material parameters 
$\epsilon_{sub}$ and 
$\sigma_{bed}$.
Segments 2, 5 and 7 in Fig. 2 were diagnosed as homogeneous in our initial analysis. These locations were sampled from across the study area, and each demonstrated high 
$\Phi_{RS}$, 
$\sigma_{RS}$ and 
$\chi_{RS}$ values for the original homogeneous scenario. For all three segments, we ran additional simulations in which 
$\epsilon_{sub}$ and 
$\sigma_{bed}$ were increased over the entire bed DEM. Figure 7 compares the changes in fit quality metrics 
$\Phi_{RS}$, 
$\sigma_{RS}$ and 
$\chi_{RS}$ for these simulations.
Fit quality metric dependence on (a) 
$\epsilon_{sub}$ and (b) 
$\sigma_{bed}$ for homogeneous simulations segments 2, 5 and 7.
Figure 7a shows the absolute change in fit quality metrics, compared to the initial value, when 
$\epsilon_{sub}$ was varied from 5 to 25 and 
$\sigma_{bed}$ remained constant at 0.2 m. Altering 
$\epsilon_{sub}$ over the entire model domain impacts the absolute returned power, but this change is neglected when evaluating relative reflections. Therefore, as expected, Fig. 7a demonstrates that the homogeneity diagnosis is relatively insensitive to the bed permittivity. Correlation coefficient (
$\Phi_{RS}$) and standard deviation ratio (
$\sigma_{RS}$) change by no more than 
$0.04$ over the range of bed permittivity tested. Range ratio (
$\chi_{RS}$) exhibited moderate sensitivity up to 
$0.17$. From these data, we propose that it is possible to use a broad range of dielectric constants without impacting the diagnosis of homogeneity or heterogeneity.
In Fig. 7b, 
$\sigma_{bed}$ varies from 0.2 to 1 m while 
$\epsilon_{sub}$ is held at 5. In these simulations, the correlation fit metric 
$\Phi_{RS}$ deviates by less than 
$0.10$ from the original value. This weak dependence may partly reflect the 720 m Savitzky–Golay filter applied during processing to suppress high-frequency noise. Shorter filter lengths may increase the sensitivity of 
$\Phi_{RS}$ to 
$\sigma_{bed}$. In contrast, 
$\sigma_{RS}$ and 
$\chi_{RS}$ change by 
$0.38$ and 
$0.35$, respectively. Increasing the bed roughness impacts returned energy by altering the small-scale bed geometry. The results in Fig. 7b demonstrate that the relative magnitude of 
$R_{beds,S}$ variation due to large-scale topography decreases as bed roughness increases. We infer from these data that the method for diagnosing homogeneity outlined in this study is sensitive to basal roughness.
Constraints on basal roughness over short distances are challenging to ascertain beneath Thwaites Glacier. Our homogeneous simulations impose a random Gaussian roughness 
$\sigma_{bed}=0.2$ m over correlation length 
$l_c=15$ m by default. As mentioned in the Methodology section, studies of basal roughness in Antarctica typically measure over much longer distances, leaving uncertainty on actual roughness values at length scales relevant to radiometrics.
The roughness magnitude chosen for most simulations, 
$\sigma_{bed}=0.2$ m, is relatively smooth for a real substrate. This magnitude is sufficient to reduce simulation artifacts such as Bragg resonance (Gerekos and others, Reference Gerekos, Bruzzone and Imai2020). We also speculate that this value may approximate glacially eroded or soft substrates beneath fast-moving ice, as is likely over large parts of our study area. As in Pierce and others (Reference Pierce2024a), the upper bound on basal roughness magnitude (
$\sigma_{bed}=1$ m) is derived from recently deglaciated Alpine topography measured at high resolution (Hubbard and others, Reference Hubbard, Siegert and Mccarroll2000). We used this value throughout the study to simulate isolated locations of bed roughness.
While we do not explicitly know basal roughness over short distances near 
$l_c$, we note that mean 
$\sigma_{RS}$ over all simulations with homogeneous fits was 
$0.80$. This indicates that the average magnitude of along-track 
$R_{bed,S}$ variation is similar to the actual 
$R_{bed}$ variation. From this, we speculate that the choice for 
$\sigma_{bed}=0.2$ m in this environment may be appropriate. We acknowledge that our treatment of basal roughness in this study is simplistic. The intent is to offer a plausible interpretations of the subglacial environment. Accurate diagnosis of basal roughness and its impact on RES over short length scales requires a more comprehensive analysis and represents a logical extension for this RES simulation methodology in future work.
RES simulation limitations
In this study, we simulated over 400 km of RES flight segments collected in 2020 and 2022. We achieved relatively high fit quality between 
$R_{bed}$ and 
$R_{bed,S}$ in 
$80\%$ of the simulated segments, which we use as evidence for the level of relative heterogeneity in the bed properties at each location. However, we acknowledge that the simulations represent less than 
$25\%$ of the 2020/2022 survey. The locations for many segments were explicitly chosen near the hypothesized features from Smith and others (Reference Smith, Gourmelen, Huth and Joughin2017) and Hager and others (Reference Hager, Hoffman, Price and Schroeder2022), and the RES data itself is more concentrated near the grounding zone than further upglacier. While we made efforts to sample a variety of locations across lower Thwaites Glacier, these inherent inconsistencies in spatial coverage may bias our findings. The spatial patterns identified through our simulations appear compelling. However, additional RES survey density with comprehensive simulation coverage would enable more definitive conclusions. Future RES simulator efficiencies combined with increased access to computing capacity may reduce the necessity for sampling, allow longer simulation lengths and enable development of a systematic inverse framework.
We note that the RES simulator outputs a focusable radargram similar to a real airborne RES experiment. However, computational constraints limit the SAR focusing aperture to 
$\sim350$ m, which is insufficient for simulating specularity content. In the future, we hope to address this constraint, as it will be helpful in reconciling hypotheses involving radiometric interpretation of subglacial hydrology.
Finally, we also acknowledge that our thresholds for 
$\Phi_{RS}$, 
$\sigma_{RS}$ and 
$\chi_{RS}$ were chosen to compare the relative fit quality between our 30 simulated segments. Our choices for these thresholds were somewhat subjective. The correlation coefficient threshold 
$\Phi_{RS} \geq 0.2$, for example, is objectively low for definitive characterization of the hydrological system using this method. However, all negatively correlated simulations are rejected with this threshold, ensuring along-track reflectivity trends are aligned. Most simulations labeled as ‘heterogeneous’ in this study met the 
$\Phi_{RS}$ criterion, but failed the more stringent 
$\sigma_{RS}$ and 
$\chi_{RS}$ thresholds due to multiple large reflectivity peaks and valleys not correlated with topography (Appendix D). Fitting these data with the simulator is technically possible, but requires multiple changes potentially speculative changes to the simulation domain. Rather than introduce additional uncertainties, we present the disparity in model fit from east to west as a finding for further exploration. We also note that interpreting homogeneity vs heterogeneity in other locations or with other RES instruments may require different threshold values to draw similar comparisons.
Conclusions
Our 30 simulated RES flight segments covered over 400 km across the lower Thwaites Glacier. We began with homogeneous simulations for each segment, where the basal material had uniform permittivity and roughness. After comparing these simulated outputs to real RES bed reflectivity variations, we identified 12 segments with a relatively homogeneous bed. An additional 12 simulations demonstrated a strong fit to the real along-track RES data if simple or coupled heterogeneity was introduced in the bed material over a large area of the modeled bed.
Most of the segments where a fit was achieved are west of 
$106.5^{\circ}$ W. The RES specularity content on this side of the glacier is consistently low, with limited orientation dependence. Conversely, simulated segments demonstrating poor fit between RES reflectivity variation and the homogeneous model were overwhelmingly observed in eastern Thwaites. These segments would require multiple spatial variations along the bed DEM to achieve a fit, and were diagnosed as having complex heterogeneity. This heterogeneity was geographically correlated with a significant orientation dependence in specularity content.
We propose an evolution of the hypothesis advocated in Schroeder and others (Reference Schroeder, Blankenship and Young2013), which cited asymmetry in RES specularity content as evidence of channelized hydrology across lower Thwaites. Channelized hydrology, combined with undulating topography, may be prevalent east of 
$106.5^{\circ}$ W longitude. This basal environment is associated with generally slower ice velocity. On the western side, widespread sediment over simple topography and fewer water channels may dominate, promoting faster ice flow.
The RES simulation technique presented in this work was applied more extensively than in previous efforts (Pierce and others, Reference Pierce2024a, Reference Pierce, Skidmore, Beem, Blankenship, Adams and Gerekos2024b). The method demonstrates promise as a tool for diagnosing subglacial conditions. It may be especially helpful in diagnosing reflectivity variations caused by topography vs material heterogeneity. The model will benefit from additional computational power and efficiency, enabling longer simulation segments with higher resolution DEMs. We also note that more comprehensive research into small-scale roughness and substrate permittivity would enhance the realism of our simulations.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2026.10160.
Data availability statement
Simulation parameters, surface / bed DEMs, and simulated focused radargrams referenced in this work are publicly available at: 10.5281/zenodo.19868888.
Airborne radar echo sounding data, including focused radargrams, bed and surface picks, returned echo power, and specularity content for the 2020 and 2022 Thwaites field seasons are also available in the repositories below:
2020 Airborne Radar Data: https://doi.org/10.18738/T8/KJKJUD
2022 Airborne Radar Data: https://doi.org/10.18738/T8/ITLQPY
Acknowledgements
This research was supported by the National Aeronautics and Space Administration (Award: 80NSSC20K1134), Korea Institute of Marine Science & Technology Promotion (KIMST) funded by the Ministry of Oceans and Fisheries (RS-2023-00256677; PM26020), the G Unger Vetlesen Foundation and Montana State University Office for Research and Economic Development. KIMST and Canadian Helicopters Limited (CHL) provided additional logistics, equipment and personnel supporting RES data collection over Thwaites Glacier. Computational efforts were performed on the Tempest High Performance Computing System, operated and supported by University Information Technology Research Cyberinfrastructure at Montana State University. We would like to acknowledge individual contributions from Dillon Buhl, Greg Nguyen, Anja Rutishauser, Jamey Stutz and Natalie Wolfenbarger, who were part of the 2020 and 2022 field campaigns to collect RES data over Thwaites Glacier. Choon-Ki Lee from KIMST also contributed a significant effort to the planning and logistics for the RES surveys over Thwaites.
Open access
