Exploring canyons beneath Devon Ice Cap for sub-glacial drainage using radar and thermodynamic modeling

Abstract

Sub-glacial canyon features up to 580 m deep between flat terraces were identified beneath Devon Ice Cap during a 2023 radar echo sounding (RES) survey. The largest canyon connects a hypothesized brine network near the Devon Ice Cap summit with the marine-terminating Sverdrup outlet glacier. This canyon represents a probable drainage route for the hypothesized water system. Radar bed reflectivity is consistently 30 dB lower along the canyon floor than on the terraces, contradicting the signature expected for sub-glacial water. We compare these data with backscattering simulations to demonstrate that the reflectivity pattern may be topographically induced. Our simulated results indicated a 10 m wide canal-like water feature is unlikely along the canyon floor, but smaller features may be difficult to detect via RES. We calculated basal temperature profiles using a 2D finite difference method and found the floor may be up to 18°C warmer than the terraces. However, temperatures remain below the pressure melting point, and there is limited evidence that the canyon floor supports a connected drainage system between the DIC summit and Sverdrup Glacier. The terrain beneath Devon Ice Cap demonstrates limitations for RES. Future studies should evaluate additional correction methods near complex terrain, such as RES simulation as we demonstrate here.

1. Introduction

The presence (or absence) of liquid water at the glacier bed is fundamental to a variety of sub-glacial processes. Distributed hydrological systems reduce basal friction, significantly altering ice internal stress and accelerating flow rates (Creyts and Schoof, 2009; Beem and others, 2010; Irvine-Fynn and others, 2011). Sub-glacial channels can transport liquid water and thermal energy down glacier, enhancing grounding line retreat (Schroeder and others, 2013; Young and others, 2016; Dow and others, 2020). The sub-glacial water system can also sustain unique microbial communities, which may be responsible for enhanced biogeochemical erosion of the bedrock (Priscu and others, 1999; Skidmore and others, 2005). Further, these isolated habitats may provide rare terrestrial analogs for icy environments across the solar system with potential to support life (Priscu and others, 1999; Chyba and Phillips, 2002).

The ability to quickly image a large geographic area has made RES appealing for decades in the study of sub-glacial hydrology (e.g. Neal, 1979; Shabtaie and Bentley, 1987; Gogineni and others, 1998; Carter and others, 2007; Schroeder and others, 2013; Rutishauser and others, 2022). Observations of elevated basal relative reflectivity (≥15 dB over ambient) and high specularity content in RES data are commonly interpreted to indicate sub-glacial water (Carter and others, 2007; Schroeder and others, 2013, 2015; Young and others, 2016). These metrics coupled with hydraulic potential modeling represent the conventional method for identifying sub-glacial hydrology with RES observations. Despite its broad appeal, spatial heterogeneity of ice properties and bed material creates significant uncertainties, complicating interpretation of bed reflectivity from RES data. Synthetic aperture processing of coherent radar helps focus the beam pattern, reducing off-nadir clutter and improving along-track resolution. In complex topography, however, significant ambiguity in locating the bed can persist (Scanlan and others, 2020). Radar attenuation loss also varies with ice thickness, as well as temperature and impurity content of the ice. Common attenuation correction methods assume varying degrees of homogeneity in the ice and bed physical properties (Matsuoka and others, 2012; Schroeder and others, 2016), which are not always realistic. These ambiguities imply that while RES is a useful tool for characterizing sub-glacial hydrology, additional context such as geothermal heat flux, hydropotential modeling, and additional experimental methods (geophysical observations, direct access, and modeling) may strengthen or contradict RES interpretations (e.g. van der Veen and others, 2007; Murray and others, 2008; Lindzey and others, 2020; Killingbeck and others, 2021).

Until recently, basal temperatures beneath Devon Ice Cap (DIC), Devon Island in Nunavut, Canada supported the assumption of a frozen bed, largely devoid of an extensive hydrological system (Van Wychen and others, 2017; Rutishauser and others, 2018). A 2018 airborne radar echo sounding (RES) survey conducted by the University of Texas Institute for Geophysics (UTIG) identified several radar-bright, specular bed features within a hydraulically flat environment near the DIC summit. These features were interpreted as a system of sub-glacial lakes, and surrounding elevated radar reflectivity was interpreted as a distributed brine network (Fig. 1), which must be hypersaline to exist at DIC basal temperatures of −14°C or colder (Rutishauser and others, 2018, 2022).

Figure 1. Map of airborne RES surveys using the UTIG MARFA instrument in 2018 (gray) and 2023 (red). Location of sub-glacial lakes (light blue) and brine network (dark blue) proposed by Rutishauser and others (2018) and Rutishauser and others (2022) are near the summit of DIC.

The Rutishauser and others (2018) sub-glacial lakes were inferred primarily from RES data. Subsequent transient electromagnetics (TEM) and sonar studies conducted in the area indicate that a present-day sub-glacial lake is unlikely (Killingbeck and others, 2024), however no such counter-evidence exists for the brine network discussed in Rutishauser and others (2022). The original hypothesis envisioned a channelized network transporting brine from a distributed system near the DIC summit to the Sverdrup Glacier terminus. Similar transitions from distributed to channelized hydrology have been identified throughout the cryosphere in freshwater systems, although the mechanisms for their formation are poorly understood (Schroeder and others, 2013; Flament and others, 2014; Jamieson and others, 2016).

If confirmed, the brine network and associated channel network described in Rutishauser and others (2018) and Rutishauser and others (2022) represent a unique unexplored environment in the Earth's cryosphere. This opportunity motivated an additional airborne RES survey in Spring 2023, focused on the Sverdrup catchment (Fig. 1). In this article, we specifically consider the region between the DIC summit and two subaerial canyons confining the main tributaries of Sverdrup Glacier (Fig. 1). Hereafter, we will refer to this as the ‘transition zone,’ and the glacier flowing downstream from the transition zone as ‘Sverdrup Glacier.’ The transition zone is critical to establishing the existence of pathways between the hypothesized brine network and Sverdrup Glacier, and may help to further our understanding of how such hydrological transitions evolve.

A RES survey from the year 2000 and described in Dowdeswell and others (2004) mentions ‘steep-sided valleys’ and ‘bedrock troughs’ beneath the DIC. These features were qualitatively associated with accelerated ice flow, possibly due to thermodynamic conditions favoring basal melt (Van Wychen and others, 2012). Our RES survey conducted in Spring of 2023 details canyon features in the transition zone, consistent with the Dowdeswell and others (2004) description. We analyze one canyon reaching a depth of 580 m as a possible drainage route for the hypothesized sub-glacial water system proposed by Rutishauser and others (2022).

The investigation of any sub-glacial hydrology beneath the DIC transition zone must consider the influence of canyons on geothermal heat flux, which may affect basal temperatures and attenuate the radar signal. The potential for canyons to concentrate geothermal heat flux and drive local temperature anomalies has been discussed since the early 1900's (Lees, 1910). Lachenbruch (1968) considered heat conduction near topographic features, such as canyons, concluding that heat flux can be amplified many times the regional average along canyon walls and floors. In the sub-glacial context, this process could create a region warm enough for liquid water to exist in an otherwise cold-bedded glacier. Lachenbruch's model was applied to canyons beneath Jakobshavn Isbræ in Greenland, where elevated geothermal flux was predicted to enhance basal melt rates (van der Veen and others, 2007). More recently, Willcocks and others (2021) modeled 2D englacial temperature profiles near deep valleys using a finite difference approximation of the heat equation. Willcocks and others (2021) used a spatially invariant heat flux boundary condition, and found a more muted topographic effect on basal temperatures than predicted by Lachenbruch (1968) or van der Veen and others (2007).

We cannot test the extant sub-glacial hydrology hypothesis beneath the DIC transition zone without considering the effects of canyon morphology on thermodynamics and radar reflectometry. Here, models of both the englacial temperatures and radar reflectometry are applied to assist in the interpretation of radar data for identifying sub-glacial water. The results of the study show conditions are inconsistent with the existence of sub-glacial freshwater beneath the DIC transition zone. We cannot explicitly rule out hydraulic features less than 10 m wide on the canyon floor, especially if they are saline, as proposed in Rutishauser and others (2018).

More generally, this work identifies improvements to thermodynamic model boundary conditions and inherent complications of interpreting radar reflectometry in steep sub-glacial terrain. The radar modeling results suggest that bed geometries can reduce basal returns by influencing attenuation and reflectivity interpretations. We suggest that regions with complicated terrain would benefit from explicitly including a bed terrain correction. These ambiguities, taken together, can confound the identification of sub-glacial water by RES observations.

2. Data and methods

2.1 Radar data and methods

The 2023 DIC RES survey was conducted using UTIG's 60 MHz multifrequency airborne radar sounder with full-phase assessment (MARFA) instrument on an AS-350 B2 helicopter platform, as described in Lindzey and others (2020). This platform enabled closer across-track spacing (tighter turning radius) and along-track resolution (15 m, slower flight speed) near steep topography than the 2018 survey (22 m), which used a similar radar mounted on a DC-3 Basler airframe. Data were range compressed, corrected for geometric spreading loss and aircraft position, and focused in azimuth as described in Peters and others (2007).

Basal returns are identified in the radar data through manually determined brackets (one above and one below), encapsulating the identified basal reflector. The magnitude of the basal return P _r is taken as the sample with the maximum value between those bounds, at each vertical trace along-track.

We estimate local englacial attenuation using a regression model (Jacobel and others, 2009; Schroeder and others, 2013, 2016). Bed returned power P _r in decibels (as denoted by [ ]_dB) is assumed to have a linear relationship with ice thickness d _ice after correcting for geometric spreading loss G (Eqn (1a)). Geometric spreading is a function of aircraft height h, ice thickness d _ice, and ice refractive index η _ice (Eqn (1b)) (Peters and others, 2005; Haynes, 2020). Other parameters in Eqn (1a) such as system gain S, birefringent losses B, and bed reflectivity R _bed are assumed to be relatively homogeneous within the subset of nearby radar observations. The resulting best-fit slope between ice thickness and geometrically corrected bed power is twice the mean one-way attenuation rate <N _ice > in dB km⁻¹.

(1a)$$\left[P_r \right]_{dB} + \left[G \right]_{dB} = \left[R_{bed} \right]_{dB} - 2\, d_{ice}\, \left< N_{ice} \right> + \left[S \right]_{dB} + \left[B \right]_{dB}$$

(1b)$$G = ( 4\pi) ^{-3}\eta_{ice}^{-2}\left(h + d_{ice}/\eta_{ice} \right)^{-4}$$

If we consider the entire 2023 survey population using this model, we estimate <N _ice > =17.5 dB km⁻¹ (Fig. 2). This attenuation rate is more than 4 dB km⁻¹ lower than the previous estimate from Rutishauser and others (2022), which applied the same linear regression method to the 2018 UTIG RES survey data. Our value also aligns well with the mean attenuation rate of 17 dB km⁻¹ derived in Killingbeck and others (2024) using the alternative Arrhenius method. However, the root-mean-squared error (RMSE) between corrected bed power ([P _r]_dB + [G]_dB) and the regression fit is 9.2 dB, indicating high uncertainty in the attenuation loss.

Figure 2. Attenuation rate determined by linear regression fit of ice thickness vs bed returned power (corrected for geometric spreading). The full RES survey is shown in blue, while the ‘Near Canyon’ data (orange/red) is the subset of points located within 3 km of Canyons A and B from Figs. 6a, b.

A geographically constrained footprint may minimize heterogeneity in bed and ice properties, improving the fit between Eqn (1a) and the RES data (Schroeder and others, 2016). For this study, we also re-calculated the attenuation regression model using only a subset of the data within 3 km of the canyons identified in the transition zone (see Fig. 6). For this data subset, which we refer to as ‘near-canyon,’ our estimate of <N _ice > = 33.6 dB km⁻¹, with a RMSE of 10.7 dB (Fig. 2). This alternative calculation is much higher than previous attenuation estimates and does not reduce uncertainty. Therefore, we use the original estimate of 17.5 dB km⁻¹ in all subsequent calculations unless otherwise noted.

2.2 RES simulation method

We used the RES simulator from Gerekos and others (2018), to create modeled radargrams with topography similar to canyons identified beneath the DIC. This method, applied to study sub-glacial hydrology as described in Pierce and others (2024), enables us to quantify the theoretical RES response to changes in assumed bed conditions. The simulator ingests a digital elevation model (DEM) with resolution l _f = 5 m, derived from the radar observations of bed and surface elevation along a flight line (Fig. 3a). Each facet within the bed and surface DEM are assigned dielectric constants $\epsilon$ for ice, rock, or water based on the location and hydrological feature we wish to model (Table 1).

Figure 3. (a) 3D view of RES simulation DEM and flight path over canyon topography at DEV3_PER0a_Y86a (label ‘B’ in Figs. 6a, b). Note that simulation increments y _i along the flight path are not drawn to scale. Actual simulated increments in y are v _a/PRF = 1 m (Table 1) (b) Cartoon of ray tracing in a 2D cross-section of a simulation. Given our simulation radius R = 300 m and l _f = 5 m, returned power is estimated by tracing over 10⁴ rays at each position increment y _i.

Table 1. Summary of parameters as described in Pierce and others (2024) used in RES simulations of DEV3_PER0a_Y86a

^aNative PRF for MARFA radar is 6400 Hz. Simulations are run at 30 Hz to replicate effective PRF after coherent summation during data processing.

^bFujita and others (2000), Tulaczyk and others (2000) and Peters and others (2005).

^cWe chose a value near the low end of the range reported for saturated bedrock (Tulaczyk and others, 2000) to maximize contrast between rock and ice.

^dA value between seawater and freshwater, as reported in Peters and others (2005) was used.

The simulator estimates the returned power profile P _r for along-track location y _i by superimposing the returned power from many individual rays directed at the ice surface within a simulation radius R. The energy's path and field strength are traced through the ice surface, then reflected from the bed (Fig. 3b). By compiling the data for each discrete location y _i along a flight path, we build a simulated radargram to approximate returned power from the topographic and material model applied. Comparable to the actual radar data, simulation outputs are range compressed and focused (Pierce and others, 2024).

Attenuation loss rate is controlled in the model through the complex ice dielectric constant, $\epsilon _{ice} = \epsilon _{ice}'( 1 + i\, \tan \delta )$. The loss tangent ($\tan \delta$) is related to attenuation rate via Eqn (2), where λ = 5 m is the radar center wavelength, and $\eta _{ice} = \sqrt {\epsilon _{ice}'}$ is the ice refractive index. We chose $\epsilon _{ice} = 3.18 ( 1 + .0018i)$, where $\tan \delta$ was selected to emulate the calculated 1-way attenuation rate <N _ice > −17.5 dB km⁻¹, and $\epsilon _{ice}'$ is well documented in the literature (Fujita and others, 2000; Peters and others, 2005). As with the real RES data, we can convert between raw power P _r and R _bed using Eqn (1a). Unlike the real RES data, <N _ice > was explicitly defined through the loss tangent, meaning there is no uncertainty in modeled attenuation loss.

(2)$$\left< N_{ice} \right> \,= \,{20\pi \eta_{ice}\log_{10}e\over \lambda}\, \tan \delta$$

We modeled two scenarios near the intersection of Canyon B and DEV3_PER0a_Y86a, as outlined below.

1. 1. Frozen bed: all bed DEM facets have dielectric constant $\epsilon _{rock}$.

2. 2. 10 m canal: a strip of facets 10 m wide has dielectric constant $\epsilon _{H_2O}$. All other facets have dielectric constant $\epsilon _{rock}$, as in the frozen bed simulation.

There are a wide range of alternative possible water structures at the Canyon B floor, including Röthlisberger channels, Nye channels, or saturated tills. We have chosen a continuous, flat canal because its geometric profile and dielectric contrast with the surrounding bed material will produce the greatest reflectivity response. We speculate that any continuous channels through the transition zone are likely to have a footprint smaller than 10 m, due to the limited hypothesized source of liquid water draining from the DIC summit, and minimal basal melt from the generally frozen bed. Pierce and others (2024) also demonstrated that the reflectivity response curve for a simulated flat canal is asymptotic above 10–20 m. Therefore, our simulation scenarios represent end-members on the continuum of possible hydrologic structures. Simulation results will diagnose the radar's overall sensitivity to water within this topographic environment, but may not conclusively determine the precise water structure without further modeling.

In both simulations, we superimpose random isotropic Gaussian roughness on the bed DEM, with σ _bed = 0.2 m or 0.35 m over a correlation length of 15 m. In part, this is added to reduce the potential for simulation artifacts such as Bragg resonance (Gerekos and others, 2020). This is also done to provide an approximation of divergent radar energy scattering due to basal roughness. A more detailed discussion of basal roughness, including rationale for the choices of σ _bed and l _c are found in Appendix C.

Because the modeled attenuation is known, variation in R _bed must result from changes in topography or bed material. Therefore, comparing simulated results to actual RES data helps to constrain the possible sub-glacial environmental conditions that may produce a given along-track reflectivity profile. It represents an alternative, and perhaps more selective approach to simply equating sub-glacial water with high values of R _bed.

All other relevant parameters used in the RES simulation setup are tabulated in Table 1. These values were chosen to mimic the geometry and radiometric properties of a helicopter-based UTIG MARFA radar experiment.

2.3 2D thermodynamic modeling method

We take a 2D modeling approach similar to Willcocks and others (2021) in order to estimate thermodynamic conditions at the bottom of the canyon features beneath the transition zone. However, our approach explores two alternative basal heat flux boundary conditions near the 580 m canyon beneath DIC. The model constrains the likely basal thermal regime, which we consider in conjunction with the RES evidence to evaluate the possibility of a present-day drainage network in the Sverdrup transition zone.

We assume the englacial temperature (T) is described by the steady-state heat equation with advection as presented in Paterson (1994) (Eqn (3a)). The form of Eqn (3a) assumes the ice has constant thermal diffusivity k, defined as the ratio between thermal conductivity K and the product of specific heat capacity c _p and density ρ (3b). Thermal advection occurs due to ice motion with ice velocity vector ${\bar {\bf w}}$.

(3a)$$k\nabla^2\, T - {\bar{\bf w}} \cdot \nabla T = 0$$

(3b)$$k = {K\over \rho \, c_p}$$

As in Willcocks and others (2021), we create a simplified 2D model representing the topographic cross-section desired. The ice surface and bed geometry are approximated with a series of plane slopes, as shown in Fig. 4 for a canyon feature identified in survey flight DEV3_PER0a_Y86a. The simplified canyon has depth d, bottom width c _b, nominal terrace ice thickness h, and canyon walls sloped at angle θ. Figure 4 superimposes bed and surface elevation profiles derived from radar transect DEV3_PER0a_Y86a to demonstrate the similarity between actual canyon topography and a modeled geometry with d = 530 m, h = 170 m, $\theta = 45^\circ$, and c _b = 380 m.

Figure 4. Representative model geometry is superimposed over actual canyon cross-section derived from radar bed and surface picks for DEV3_PER0a_Y86a (label ‘B’ in Figs. 6a, b). Similar figures demonstrating modeled and actual canyon topography at all modeled radar transects are found in Appendix B, Figs. 18–27.

The idealized surface and bed geometry is discretized into a 2D matrix with resolution Δx = Δz = 10 m. Equation (4a) is the central finite difference approximation for Eqn (3a), where i and j are indices in the horizontal (x) and vertical (z) directions. Each pixel above the canyon is assigned ice thermal diffusivity k from Table 2. Vertical velocity is described by (4b) above the ice/rock boundary. The annual ice accumulation rate, b = 0.17 ± 0.05 m yr⁻¹ was based on the 1963–2000 mass balance performed in Mair and others (2005) across four shallow core sites within the upper Sverdrup catchment. In order to confine the model domain to two dimensions, we assume advection only in the z direction and neglect ice motion in x and y (Paterson, 1994).

(4a)$$\eqalign{T_{i, j} & = {1\over 4}\left[T_{i-1, j} + T_{i + 1, j} + T_{i, j-1} \left(1 + {w_z \Delta x\over 2k} \right)\right.\cr & \quad\left. + T_{i, j + 1} \left(1 - {w_z \Delta x\over 2k} \right)\right]}$$

(4b)$$w_z = -bz_{i, j}/h$$

Table 2. Summary of parameters used in 2D heat transfer model

^aPaterson (1994).

^bMair and others (2005).

^cRutishauser and others (2018).

The 2D thermodynamic model iteratively solves Eqn (4a). The model converges at iteration m, when the criterion in Eqn (5) is met.

(5)$$\sum_i\sum_j{{\left(T_{i, j}^m - T_{i, j}^{m-1}\right)^2}} \leq0.001$$

The boundary conditions for the model domain are:

1. 1. Constant ice surface temperature T _s. We assume $T_s = -23^\circ {\rm C}$ at the DIC summit (1900 m), and increases 4.4°C per km of elevation loss (Gardner and others, 2009; Rutishauser and others, 2018).

2. 2. No horizontal heat flux at the left and right side of the domain (dT/dx = 0).

3. 3. One of the following basal heat flux boundary conditions:

   1. (a) Constant: Vertical heat flux applied the ice/bedrock interface is spatially constant and equal to the regional average geothermal heat flux, Q _z(x) = 65 mW m⁻¹ K⁻¹ (Rutishauser and others, 2018).

   2. (b) Lachenbruch: Q _z(x) = q _an(x) ⋅ 65 mW m⁻¹ K⁻¹ varies along the the x direction at the ice/bedrock interface as defined in Lachenbruch (1968). The solution for q _an(x) is tabulated in Lachenbruch (1968) for θ = 15 − 90°.

Figure 5 compares the basal heat flux boundary conditions 3a and 3b for the simulated geometry in Fig. 4. The Lachenbruch boundary condition predicts large variations in heat flux between the canyon walls, and is most applicable if the thermal diffusivity contrast at the interface is high. As thermal diffusivity contrast between bedrock and ice approaches zero, basal heat flux will approach the Constant boundary condition. Therefore, although we cannot know the precise variation in basal heat flux Q _z(x) along the canyon, comparing temperature profiles using both boundary condition 3a and 3b provides a useful constraint on temperature uncertainty due to the topography beneath DIC.

Figure 5. Geothermal heat flux boundary conditions for DEV3_PER0a_Y86a. The Lachenbruch boundary condition exhibits large variations in heat flux along the x direction, which is dependent upon the canyon geometry.

For all simulations, we did not incorporate latent heat into the thermodynamic model. Therefore, similar to Willcocks and others (2021), any resulting temperatures in excess of the ice pressure melting point represent melt potential, and a likely presence of liquid water. We also neglect any frictional heating in this model. We expect basal drag to be negligible across most of the DIC transition zone, with the possible exception of the location nearest Sverdrup Glacier at DEV3_PER0a_Y87a (Van Wychen and others, 2017).

3. Results

3.1 RES results

The 2023 survey included 10 transects running approximately perpendicular to ice flow across the transition zone (Figs. 6a, b). The flight closest to Sverdrup Glacier (DEV_PER0a_Y88a) was characterized by locations of thin ice (<100 m) with significant clutter from nearby canyon walls, and consequently is excluded from our analysis. The remaining 9 flights each demonstrate broad areas of elevated bed reflectivity between narrow regions with dim R _bed (Fig. 6a). The pattern in specularity content (Fig. 6b) co-varies with relative bed reflectivity.

Figure 6. (a) Relative reflectivity and (b) specularity content from the 2023 RES survey in the transition zone between the DIC summit and Sverdrup Glacier. The locations of A and B on DEV3_PER0a_Y86a are shown. The gray diamonds trace the canyon's locations from the hypothesized brine network to Sverdrup Glacier. (c) 1D focused radargram for DEV3_PER0a_Y86a. Points A and B are the canyon features mapped on panels a and b above. (d) Relative reflectivity and specularity content for DEV3_PER0a_Y86a. (e) and (f) zoomed in images of Canyons A and B from the radargram in (c).

The 1D focused radargram for one of these flights, DEV3_PER0a_Y86a, shows the existence of at least two narrow canyons lying between broad, flat terraces under the DIC (Fig. 6c). The canyons are labeled A and B in Fig. 6, and have estimated depths of 370 m and 530 m along the DEV3_PER0a_Y86a flight path, respectively. A and B are located directly up-glacier from Sverdrup's main tributaries (Figs. 6a, b).

The DEV3_PER0a_Y86a topographic lows at A and B are strongly correlated with minima in both R _bed and specularity content (Figs. 6c, d). The contrast in R _bed between the terraces and the canyon floor is 30–35 dB. Each of the parallel survey flights across the transition zone exhibit a similar relationship between topography, reflectivity, and specularity content (see Appendix A, Figs. 9–17). Figures 6a, b plot the adjacent topographic lows for each flight line, and Canyon B's path can be traced from the Rutishauser and others (2022) hypothesized brine network location to Sverdrup Glacier. Near the proposed brine network, Canyon B has a depth of ~100 m, while ice 410 m thick covers the terraces. As you follow Canyon B toward the Sverdrup Glacier, radar shows the ice thinning to around 100 m, while the canyon deepens to 580 m (Table 3).

Table 3. Model parameters used for each intersection between survey flight line and Canyon B

^aDistance measured from nearest edge of proposed brine network depicted in Fig. 6 along Canyon B path.

^bIce surface elevation at location where survey line intersects Canyon B.

^cExample survey line referenced in all figures.

If any sub-glacial water exists, it would follow the hydraulic potential gradient from the brine network to Sverdrup Glacier (Rutishauser and others, 2022). Sub-glacial hydraulic potential (ψ) is defined in Eqn (6) as a function of surface elevation z _s and ice thickness d _ice, where ρ _ice and $\rho _{H_2O}$ are density of ice and water, respectively. In the transition zone the surface topography is relatively smooth, sloping toward Sverdrup Glacier. By contrast, the steep canyon topography results in large ice thickness changes. Therefore the canyons exert strong control on the local hydraulic potential, which will drive any water toward the canyon floor. The hydropotential analysis in Rutishauser and others (2022) corroborates a likely drainage pathway from the hypothesized brine network to Sverdrup Glacier corresponding to Canyon B's location.

(6)$$\psi = d_{ice}\left({\rho_{ice}\over \rho_{H_2O}} - 1 \right) + z_s$$

Elevated bed reflectivity (≥15 dB above surroundings) and/or specularity content in RES data is frequently correlated with sub-glacial water (Peters and others, 2005; Schroeder and others, 2013, 2015; Young and others, 2016). However, the area beneath the DIC transition zone demonstrates the opposite radiometric signature. If there was sub-glacial water flowing out of the hypothesized brine-network, the canyons are the most likely locations for subglacial water, yet, they exhibit the lowest lowest reflectivity and specularity content.

Specularity content attempts to quantify the diffuse or specular quality of radar reflections from a target by comparing returned power over two different apertures (Schroeder and others, 2015). A surface oriented normal to the radar signal's travel path with low roughness compared to the radar wavelength will exhibit high specularity content. The surface of a sub-glacial lake often meets this description. However, in the DIC transition zone topography, the terrace features are broad and relatively flat over our processing aperture lengths (700 m and 2 km). By contrast, the width of Canyon B from brink to brink is less than 2 km, with large topographic variation in between. Therefore, we expect higher specularity content on the terraces than in the canyons, even if sub-glacial water is present.

The strong positive connection between bed reflectivity and topography has implications to the interpretation of sub-glacial water within the transition zone topography. Conventional RES analysis might altogether exclude the Canyon B floor from consideration as having existing hydrology. However, we observe a 10 dB local reflectivity peak coincident with the Canyon B floor (Fig. 6d). This could be a reasonable increase for a canal-like structure in the canyon. Within the context of a >30 dB reduction in bed reflectivity across Canyon B  such a radiometric signature will appear dim. Therefore, we consider the RES simulation results for interpreting the along-track R _bed profile and assessing the potential sensitivity of radiometrics to liquid water in this sub-glacial environment.

3.2 RES simulations of Canyon B topography

The simulated data in Fig. 7 indicate bed reflectivity reductions up to 30 dB from the terraces to Canyon B's floor in each scenario, irrespective of the roughness value σ _bed. The pattern mimics the trend in the actual RES data for DEV3_PER0a_Y86a in the vicinity of Canyon B. Both simulations exhibit a small local peak in R _bed at the Canyon B floor. In the frozen bed scenario, this peak is 10 dB above the minimum R _bed value (Fig. 7d). The results are virtually indistinguishable when σ _bed increases from 0.2 m to 0.35 m.

Figure 7. Simulated relative bed reflectivity results along an 8 km segment of DEV3_PER0a_Y86a, centered on Canyon B, (a) with a 10 m flat canal at the floor and (b) with a frozen bed and no canal. (c) and (d) zoomed view of R _bed response near the bottom of Canyon B.

The presence of a 10 m canal increases the magnitude of this small peak at the Canyon B floor to 16.2 dB when σ _bed = 0.2 m (Fig. 7c). The peak increases to 20.5 dB when roughness increases to σ _bed = 0.35 m. This indicates a 6.2–10.5 dB response resulting from the canal vs the frozen bed. Similar simulations from Pierce and others (2024) of a 10 m canal in flat topography resulted in a 14–20 dB response, suggesting that the canyon topography attenuates the reflectivity response. The power returned from the canal also remains at least 10–15 dB below bed reflectivity along the terraces. This result implies sub-glacial water in the canyon, even a specular canal of appreciable size, would appear dim on a radargram. Thus, the radar signature of such a canal could be easily overlooked in RES interpretation.

The magnitude of the reflectivity peak at the bottom of Canyon B closely matches the frozen bed simulation (Fig. 7d). Our simulations also demonstrate that RES has low sensitivity to liquid water in this environment. The simulation scenarios represent the end members for a spectrum of possible hydrological structures in this location. RES is unlikely to determine the presence of alternative structures with lower reflectivity signatures such as Röthlisberger channels, Nye channels, or hydrated tills (Pierce and others, 2024). Therefore, the simulation results imply the RES data is inconclusive regarding the presence of liquid water. We move on to consider local thermodynamic conditions at the Canyon B floor, to assess the potential for the presence of water, and thus the potential for such alternative water bodies.

3.3 2D thermodynamic model results

Figure 8 shows modeled bed temperatures along the 14 km path of Canyon B, using the parameters in Table 3. The Canyon B floor temperature increases from −10.0 ± 1.0°C near the brine network to −3.0 ± 1.0°C closest to Sverdrup Glacier, assuming constant basal heat flux. Canyon B floor temperatures are 1.3–4.5°C higher when the Lachenbruch heat flux boundary condition is applied, reaching a maximum of 1.4 ± 1.1°C near Sverdrup Glacier. Along the terraces, ice thins over the same path, corresponding to a decrease in bed temperatures from −12.6°C to −16.9°C. Bed temperatures along the terrace are in strong agreement for both heat flux boundary conditions.

Figure 8. (a) Modeled 2D bed temperature profile locations in the Sverdrup transition zone. Colors indicate temperature of the Canyon B floor and terraces (constant boundary condition). The grayscale indicates ice thickness, derived from ArcticDEM (Porter and others, 2018) and bed elevation data provided in Rutishauser and others (2018). Location of Rutishauser and others (2022) hypothesized brine network is shown in blue. (b) Comparison of modeled bed temperatures for the Constant and Lachenbruch boundary conditions along Canyon B. The points represent the positions where the flight lines intersect Canyon B and bed temperatures were modeled. Shaded areas represent the modeled temperature range due to uncertainty in accumulation rate b.

We have not included latent heat into our model, therefore any temperatures above the ice pressure melting point indicate freshwater melt potential (Fig. 8b). Saline hydrology could exist at lower temperatures, with freezing point depression dependent on salt concentration (Badgeley and others, 2017; Rutishauser and others, 2022).

4. Discussion

4.1 Existence of liquid hydrology

Our 2D thermodynamic model results indicate a continuous, saline drainage network is thermodynamically possible in the Sverdrup transition zone. Temperatures between −7.7°C and −11.3°C closest to the hypothesized brine network would sustain liquid at salinity lower than previously proposed in Rutishauser and others (2022). However, the findings from Killingbeck and others (2024) indicate a nearby present-day sub-glacial lake (Rutishauser and others, 2018) is unlikely. By extension, it is not unreasonable to question the existence of the brine network proposed in Rutishauser and others (2022). Combined with the radar evidence presented in this paper, brine drainage channels of sufficient size for RES detection are improbable in the DIC transition zone.

The temperature modeling does demonstrate significant temperature deviations of 2–18°C between the floor and terraces along Canyon B. Although modeled basal temperatures are elevated throughout the canyon, temperatures do not exceed the ice pressure melting point along most of the path, precluding a freshwater network. At approximately 14 km from the proposed brine network, near DEV3_PER0a_Y86a, the pressure melting point could be reached if the Lachenbruch boundary condition is assumed, but remain 3°C below freezing if the heat flux is constant.

The temperature models are generally consistent with the RES data and simulation results. The bed reflectivity patterns observed in the survey flights we examined imply either a frozen Canyon B floor, or one that contains small, disconnected liquid water features. We have also demonstrated through both simulated and real RES data that returned energy is reduced by up to 30 dB in the canyons. This low returned energy in the canyon appears to reduce the radar's sensitivity to liquid water in our simulations, making small-scale hydrological features challenging to interpret. It is important to note, however, that we simulated a 10 m channel in a simulation environment with a 5 m resolution. There is some risk a feature spanning only two facets in the backscattering model will be less prominent than in reality.

4.2 Thermodynamic modeling and sub-glacial hydrology

Both boundary conditions in our 2D thermodynamic model produced generally congruent results for bed temperatures along the Canyon B terraces and floor. The difference in mean terrace temperature between the two boundary conditions was less than 0.5°C along the entire canyon length. At the canyon floor, both boundary conditions demonstrate a similar trend, although the Lachenbruch boundary condition predicts warmer temperatures due to the large variation in heat flux along the topography (Fig. 5).

In previous work, the Lachenbruch (1968) heat flux anomaly has been used as evidence for elevated bed temperatures within sub-glacial canyon features (van der Veen and others, 2007). Lachenbruch provides a solution for heat flux along a simplified 2D topographic surface, assuming the domain is bounded in the z direction at that surface. The derivation does not consider thermal contact with another material above the topographic surface, such as ice, which may insulate or enhance energy conduction. The Lachenbruch boundary condition may be a reasonable approximation when thermal diffusivity contrast at the interface is large and positive. This is true in the case of sub-aerial topography ($k_{air} \gtrapprox 10\, k_{rock}$). However, the bedrock underlying the transition zone is likely comprised of layered gneiss, sandstone, and limestone (Okulitch, 1991). These materials have a mean thermal diffusivity of 1.1 mm² s⁻¹ (Seipold and Huenges, 1998; Waples and Waples, 2004; Tiwari and others, 2020), matching the value for ice in Table 2. Although no other study has thoroughly evaluated sub-glacial heat flux near topography, we expect the Lachenbruch boundary condition overestimates heat flux because it ignores the overlying ice. Therefore, we contend the constant heat flux case represents basal temperatures more accurately along Canyon B.

Our 2D thermodynamic model also ignores advective heat transfer in both the x (left to right) and y (into the page) directions (see Fig. 4). This allowed us to simplify Eqn (3a) for a 2D use case which considers only the topographic cross-section. Ignoring advection in the x direction is almost certainly reasonable. The coordinate system is approximately oriented with ice flow along the y-axis, indicating ice velocity and temperature gradients in the x direction are much lower than in y or z.

Ignoring advection in the y direction may be less straightforward. Assumed surface temperature gradients along Canyon B's length are two orders of magnitude lower than the largest vertical temperature gradients calculated in our 2D thermodynamic models (~10⁻⁴ vs ~10⁻²°C m⁻¹). Conversely, maximum local ice surface velocity of 20 − −40 m yr⁻¹ is two orders of magnitude higher than the accumulation rate used in our model (Van Wychen and others, 2017). Therefore, the advection term from Eqn (3a) (${\bar {\bf w}} \cdot \nabla T$) in a 3-dimensional model could be of similar magnitude along the y and z axes. For our scenario along Canyon B, this implies our 2D solution likely predicts marginally higher temperatures than an equivalent 3D model. This reinforces our primary conclusion that connected drainage is unlikely beneath Canyon B unless it is saline. More broadly, attempts to replicate the thermodynamic modeling in this paper should carefully consider the validity of a 2D vs 3D model, especially in regions of fast-flowing ice where frictional heating and lateral advection may be relevant.

4.3 Impact of topography on RES interpretation

The canyons beneath DIC's relatively thin ice demonstrate significant potential limitations for conventional RES analysis techniques. First and foremost, we have demonstrated that topography can produce significant variations in along-track relative bed reflectivity. Elevated basal reflectivity >15 dB is a common criterion used to infer sub-glacial water from RES results (Peters and others, 2005; Young and others, 2016; Rutishauser and others, 2018). Our interpretation of the 2023 DIC RES survey data demonstrates that more nuance is necessary. Conventional RES analysis would exclude the canyons altogether due to their dim reflections, while the terraces might be considered as potential locations for sub-glacial water due to their elevated R _bed of 30 dB or more. However, the terraces are the least likely location for sub-glacial water in this environment due to local hydropotential gradients. RES simulations may be a useful tool in combination with conventional techniques for resolving the likely structure, extent, and shape of sub-glacial water. This is especially true near topographic features where high variation in apparent basal reflectivity is expected.

We estimated the attenuation rate of 17.5 dB km⁻¹ using a linear regression fit of ice thickness vs geometrically corrected bed returned power ([P _r]_dB + [G]_dB from Eqn (1a)). Figure 2 shows this regression model in blue. Our attenuation estimate is 4.3 dB km⁻¹ lower than the rate of 21.8 dB km⁻¹ in Rutishauser and others (2022). Killingbeck and others (2024) indicated that the high attenuation rate in Rutishauser and others (2018), which was derived using the same linear regression technique, may have contributed to mis-identifying sub-glacial lakes near the DIC summit. They suggest the Arrhenius method (Matsuoka and others, 2012) may be more reliable for evaluating englacial attenuation rates, and estimated average attenuation over DIC ~17 dB km⁻¹ (Killingbeck and others, 2024).

Although our attenuation estimate is similar to the average value in Killingbeck and others (2024), our depth-averaged temperatures along the terraces from our 2D models ranged between −17°C to −18°C. This is colder than the ice near DIC's summit, and would translate to an expected attenuation rate for pure ice around ~14 dB km⁻¹ using an Arrhenius method (Matsuoka and others, 2012). Therefore, it is possible we have still overestimated attenuation rate in this study. Further, attenuation uncertainty exceeded 9 dB across the 2023 survey.

In Fig. 2 we also re-calculate a ‘near-canyon’ regression using a subset of the data within 3 km of the Canyon A and B locations shown in Figs. 6a, b. In this data subset, the apparent 1-way attenuation rate nearly doubles to 33.6 dB km⁻¹. Both the full RES survey and near-canyon subset exhibit similar mean bed returned power where ice thickness is less than 450 m. However, mean returned bed power in the near-canyon data is 12 dB less than the full RES survey for ice >450 m. In the near-canyon data, ice exceeding 450 m is only found between the canyon walls, while thinner ice (<450 m) is generally on the adjacent terraces. This is a significant contrast with the full RES data set, which includes locations near the DIC summit where thicker ice covers topography with smaller magnitude slopes. This bias toward low returned power near the canyons inherently affects attenuation rates when derived via linear regression. The resulting regression fit may overestimate the survey-wide attenuation rate and increase uncertainty. This can result in overestimated bed reflectivity and misinterpretation of RES data (Killingbeck and others, 2024).

The linear regression method in Equation (1a) assumes relative homogeneity in ice material properties and bed reflectivity (Eqn (1a)) to estimate attenuation. Ice material properties, including attenuation rate, are known to increase with temperature (Matsuoka and others, 2012). As we have shown through modeling, the canyon ice is likely much warmer compared to other locations beneath DIC. However, our RES simulations demonstrated a similar reduction in bed returned power despite bed reflectivity and attenuation rate values that were known and constant. This leads us to the conclusion that heterogeneity in ice temperature and the bed environment were not the primary mechanism for the reduced bed returns in the canyons.

We therefore speculate that the canyon geometry itself is the cause of the reduced echo power, and the effect manifests itself as attenuation uncertainty. The linear regression model assumes that geometric spreading of the radar signal [G]_dB is known explicitly as a function of ice thickness, surface elevation, and radar positioning (Eqn (1b)). The derivation of Eqn (1a) assumes minimal variation in geometric spreading over the radar's pulse-limited footprint (Peters and others, 2005; Haynes and others, 2018). For thick ice covering a gently sloping bed, this is often reasonable. However, we speculate that near the steep, narrow canyon walls in the DIC transition zone, d _ice variation within the pulse-limited footprint exceeds a threshold value for this assumption to remain valid. Further research should assess the need for topographic corrections at appropriate slope angle and ice thickness thresholds. Additionally, alternative attenuation correction methods which account for topographic heterogeneity may benefit radar data processing.

5. Conclusion

A 2023 RES survey of the DIC enabled detailed analysis of canyons connecting the summit area with Sverdrup Glacier. The canyons reach depths of nearly 600 m, and represent the most likely route for draining any sub-glacial brines as hypothesized in Rutishauser and others (2022). However, specularity content and bed reflectivity values are low in the canyons, contradicting the pattern expected in the presence of hydrological features. Through the use of RES simulations, we demonstrated that radar returns from sub-glacial water within the observed topography may appear dim compared to the surrounding bed, making RES a challenging tool for sub-glacial water diagnosis in this environment. Our simulated results indicated a large, specular canal (10 m wide) beneath the canyon is unlikely, however we cannot preclude the existence of smaller, less continuous water bodies based solely on RES reflectivity and specularity content.

We estimated basal temperatures beneath the deepest canyon identified using a finite difference approximation of the heat equation in two dimensions. The results indicate that temperatures in the canyons may be elevated up to 18°C above basal temperatures on the terraces. These thermodynamic conditions may enable a saline hydrological system to exist, as discussed in Rutishauser and others (2022). Conditions exceeding the pressure melting point are unlikely along most of the canyon's floor, implying sub-glacial freshwater cannot exist.

Our work also considers some limitations of conventional RES analysis applied to sub-glacial hydrology in complex terrain. Traditional interpretation of radar data equates high reflectivity and specularity content with bodies of water. We have demonstrated that this is overly simplistic in complex terrain. Comparing actual RES bed reflectivity to simulated returns, as in Pierce and others (2024), may offer a complimentary method for identifying and characterizing sub-glacial water with greater specificity. Further, common methods for estimating englacial attenuation with regression analysis may be flawed where the bed environment exhibits significant heterogeneity. Resulting attenuation estimates may exhibit abnormally high uncertainty, provoking flawed interpretations of basal conditions. Future radar analysis may consider corrections for basal topography, placing greater emphasis on returned power instead of relative bed reflectivity, or alternative attenuation estimation techniques (Matsuoka and others, 2012).