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, Reference Creyts and Schoof2009; Beem and others, Reference Beem, Jezek and Van Der Veen2010; Irvine-Fynn and others, Reference Irvine-Fynn, Hodson, Moorman, Vatne and Hubbard2011). Sub-glacial channels can transport liquid water and thermal energy down glacier, enhancing grounding line retreat (Schroeder and others, Reference Schroeder, Blankenship and Young2013; Young and others, Reference Young, Schroeder, Blankenship, Kempf and Quartini2016; Dow and others, Reference Dow2020). 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, Reference Priscu1999; Skidmore and others, Reference Skidmore, Anderson, Sharp, Foght and Lanoil2005). Further, these isolated habitats may provide rare terrestrial analogs for icy environments across the solar system with potential to support life (Priscu and others, Reference Priscu1999; Chyba and Phillips, Reference Chyba and Phillips2002).
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, Reference Neal1979; Shabtaie and Bentley, Reference Shabtaie and Bentley1987; Gogineni and others, Reference Gogineni, Chuah, Allen, Jezek and Moore1998; Carter and others, Reference Carter2007; Schroeder and others, Reference Schroeder, Blankenship and Young2013; Rutishauser and others, Reference Rutishauser2022). 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, Reference Carter2007; Schroeder and others, Reference Schroeder, Blankenship and Young2013, Reference Schroeder, Blankenship, Raney and Grima2015; Young and others, Reference Young, Schroeder, Blankenship, Kempf and Quartini2016). 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, Reference Scanlan, Rutishauser, Young and Blankenship2020). 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, Reference Matsuoka, MacGregor and Pattyn2012; Schroeder and others, Reference Schroeder, Seroussi, Chu and Young2016), 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, Reference van der Veen, Leftwich, Frese, Csatho and Li2007; Murray and others, Reference Murray, Corr, Forieri and Smith2008; Lindzey and others, Reference Lindzey2020; Killingbeck and others, Reference Killingbeck, Dow and Unsworth2021).
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, Reference Van Wychen2017; Rutishauser and others, Reference Rutishauser2018). 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, Reference Rutishauser2018, Reference Rutishauser2022).
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 (Reference Rutishauser2018) and Rutishauser and others (Reference Rutishauser2022) are near the summit of DIC.
The Rutishauser and others (Reference Rutishauser2018) 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, Reference Killingbeck2024), however no such counter-evidence exists for the brine network discussed in Rutishauser and others (Reference Rutishauser2022). 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, Reference Schroeder, Blankenship and Young2013; Flament and others, Reference Flament, Berthier and Rémy2014; Jamieson and others, Reference Jamieson2016).
If confirmed, the brine network and associated channel network described in Rutishauser and others (Reference Rutishauser2018) and Rutishauser and others (Reference Rutishauser2022) 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 (Reference Dowdeswell, Benham, Gorman, Burgess and Sharp2004) 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, Reference Van Wychen2012). Our RES survey conducted in Spring of 2023 details canyon features in the transition zone, consistent with the Dowdeswell and others (Reference Dowdeswell, Benham, Gorman, Burgess and Sharp2004) 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 (Reference Rutishauser2022).
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, Reference Lees1910). Lachenbruch (Reference Lachenbruch1968) 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, Reference van der Veen, Leftwich, Frese, Csatho and Li2007). More recently, Willcocks and others (Reference Willcocks, Hasterok and Jennings2021) modeled 2D englacial temperature profiles near deep valleys using a finite difference approximation of the heat equation. Willcocks and others (Reference Willcocks, Hasterok and Jennings2021) used a spatially invariant heat flux boundary condition, and found a more muted topographic effect on basal temperatures than predicted by Lachenbruch (Reference Lachenbruch1968) or van der Veen and others (Reference van der Veen, Leftwich, Frese, Csatho and Li2007).
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 (Reference Rutishauser2018).
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 (Reference Lindzey2020). 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 (Reference Peters2007).
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, Reference Jacobel, Welch, Osterhouse, Pettersson and Macgregor2009; Schroeder and others, Reference Schroeder, Blankenship and Young2013, Reference Schroeder, Seroussi, Chu and Young2016). 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, Reference Peters, Blankenship and Morse2005; Haynes, Reference Haynes2020). 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−1.
If we consider the entire 2023 survey population using this model, we estimate <N ice > =17.5 dB km−1 (Fig. 2). This attenuation rate is more than 4 dB km−1 lower than the previous estimate from Rutishauser and others (Reference Rutishauser2022), 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−1 derived in Killingbeck and others (Reference Killingbeck2024) 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.
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, Reference Schroeder, Seroussi, Chu and Young2016). 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−1, 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−1 in all subsequent calculations unless otherwise noted.
2.2 RES simulation method
We used the RES simulator from Gerekos and others (Reference Gerekos2018), 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 (Reference Pierce2024), 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).
(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 104 rays at each position increment y i.
Summary of parameters as described in Pierce and others (Reference Pierce2024) 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 (Reference Fujita, Matsuoka, Ishida, Matsuoka and Mae2000), Tulaczyk and others (Reference Tulaczyk, Kamb and Engelhardt2000) and Peters and others (Reference Peters, Blankenship and Morse2005).
cWe chose a value near the low end of the range reported for saturated bedrock (Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000) to maximize contrast between rock and ice.
dA value between seawater and freshwater, as reported in Peters and others (Reference Peters, Blankenship and Morse2005) 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, Reference Pierce2024).
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−1, and $\epsilon _{ice}'$
is well documented in the literature (Fujita and others, Reference Fujita, Matsuoka, Ishida, Matsuoka and Mae2000; Peters and others, Reference Peters, Blankenship and Morse2005). 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.
We modeled two scenarios near the intersection of Canyon B and DEV3_PER0a_Y86a, as outlined below.
1. Frozen bed: all bed DEM facets have dielectric constant $\epsilon _{rock}$
.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 (Reference Pierce2024) 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, Reference Gerekos, Bruzzone and Imai2020). 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 (Reference Willcocks, Hasterok and Jennings2021) 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 (Reference Paterson1994) (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}}$
.
As in Willcocks and others (Reference Willcocks, Hasterok and Jennings2021), 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.
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−1 was based on the 1963–2000 mass balance performed in Mair and others (Reference Mair, Burgess and Sharp2005) 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, Reference Paterson1994).
Summary of parameters used in 2D heat transfer model
aPaterson (Reference Paterson1994).
bMair and others (Reference Mair, Burgess and Sharp2005).
cRutishauser and others (Reference Rutishauser2018).
The 2D thermodynamic model iteratively solves Eqn (4a). The model converges at iteration m, when the criterion in Eqn (5) is met.
The boundary conditions for the model domain are:
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, Reference Gardner2009; Rutishauser and others, Reference Rutishauser2018).2. No horizontal heat flux at the left and right side of the domain (dT/dx = 0).
3. One of the following basal heat flux boundary conditions:
(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−1 K−1 (Rutishauser and others, Reference Rutishauser2018).
(b) Lachenbruch: Q z(x) = q an(x) ⋅ 65 mW m−1 K−1 varies along the the x direction at the ice/bedrock interface as defined in Lachenbruch (Reference Lachenbruch1968). The solution for q an(x) is tabulated in Lachenbruch (Reference Lachenbruch1968) 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.
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 (Reference Willcocks, Hasterok and Jennings2021), 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, Reference Van Wychen2017).
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.
(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 (Reference Rutishauser2022) 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).
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, Reference Rutishauser2022). 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 (Reference Rutishauser2022) corroborates a likely drainage pathway from the hypothesized brine network to Sverdrup Glacier corresponding to Canyon B's location.
Elevated bed reflectivity (≥15 dB above surroundings) and/or specularity content in RES data is frequently correlated with sub-glacial water (Peters and others, Reference Peters, Blankenship and Morse2005; Schroeder and others, Reference Schroeder, Blankenship and Young2013, Reference Schroeder, Blankenship, Raney and Grima2015; Young and others, Reference Young, Schroeder, Blankenship, Kempf and Quartini2016). 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, Reference Schroeder, Blankenship, Raney and Grima2015). 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.
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 (Reference Pierce2024) 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, Reference Pierce2024). 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.
(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, Reference Porter2018) and bed elevation data provided in Rutishauser and others (Reference Rutishauser2018). Location of Rutishauser and others (Reference Rutishauser2022) 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, Reference Badgeley2017; Rutishauser and others, Reference Rutishauser2022).
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 (Reference Rutishauser2022). However, the findings from Killingbeck and others (Reference Killingbeck2024) indicate a nearby present-day sub-glacial lake (Rutishauser and others, Reference Rutishauser2018) is unlikely. By extension, it is not unreasonable to question the existence of the brine network proposed in Rutishauser and others (Reference Rutishauser2022). 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 (Reference Lachenbruch1968) heat flux anomaly has been used as evidence for elevated bed temperatures within sub-glacial canyon features (van der Veen and others, Reference van der Veen, Leftwich, Frese, Csatho and Li2007). 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, Reference Okulitch1991). These materials have a mean thermal diffusivity of 1.1 mm2 s−1 (Seipold and Huenges, Reference Seipold and Huenges1998; Waples and Waples, Reference Waples and Waples2004; Tiwari and others, Reference Tiwari, Boháč, Dieška and Götzl2020), 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−4 vs ~10−2°C m−1). Conversely, maximum local ice surface velocity of 20 − −40 m yr−1 is two orders of magnitude higher than the accumulation rate used in our model (Van Wychen and others, Reference Van Wychen2017). 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, Reference Peters, Blankenship and Morse2005; Young and others, Reference Young, Schroeder, Blankenship, Kempf and Quartini2016; Rutishauser and others, Reference Rutishauser2018). 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−1 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−1 lower than the rate of 21.8 dB km−1 in Rutishauser and others (Reference Rutishauser2022). Killingbeck and others (Reference Killingbeck2024) indicated that the high attenuation rate in Rutishauser and others (Reference Rutishauser2018), 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, Reference Matsuoka, MacGregor and Pattyn2012) may be more reliable for evaluating englacial attenuation rates, and estimated average attenuation over DIC ~17 dB km−1 (Killingbeck and others, Reference Killingbeck2024).
Although our attenuation estimate is similar to the average value in Killingbeck and others (Reference Killingbeck2024), 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−1 using an Arrhenius method (Matsuoka and others, Reference Matsuoka, MacGregor and Pattyn2012). 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−1. 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, Reference Killingbeck2024).
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, Reference Matsuoka, MacGregor and Pattyn2012). 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, Reference Peters, Blankenship and Morse2005; Haynes and others, Reference Haynes, Chapin and Schroeder2018). 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 (Reference Rutishauser2022). 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 (Reference Rutishauser2022). 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 (Reference Pierce2024), 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, Reference Matsuoka, MacGregor and Pattyn2012).
Data
Radar data cited in this work are available at https://doi.org/10.18738/T8/ZN15UF. RES simulation and modeled temperature outputs can be accessed at https://doi.org/10.5281/zenodo.10934128.
Acknowledgements
We thank the Nunavut Research Institute and the communities of Grise Fiord and Resolute Bay for granting us permission to conduct research on Devon Island. This research was supported by the National Aeronautics and Space Administration (Award: 80NSSC20K1134), the G Unger Vetlesen Foundation, and Montana State University Office for Research and Economic Development. Computational efforts were performed on the Tempest High Performance Computing System, operated and supported by University Information Technology Research Cyberinfrastructure at Montana State University. Canadian Helicopters Limited (CHL) and Natural Resource Canada's Polar Continental Shelf Program (PCSP) provided logistics, equipment, and personnel supporting RES data collection over the Devon Ice Cap. We would like to acknowledge individual contributions from Dillon Buhl, Corey Skender, Clint Baumgartner, and Joshua Brett. Their expertise and dedication to the project ensured a safe and successful 2023 field survey of the DIC. Finally, we would like to thank the two reviewers for their constructive feedback, which resulted in significant improvements to the final paper.
Appendix A
1D focused radargrams, aligned with bed reflectivity and specularity content, for all survey lines running approximately perpendicular to Sverdrup Glacier ice flow. Canyon A and B are shown as appropriate. Flight paths for the radargrams are mapped in Figure 6a. Flight DEV3_PER0a_Y79a is the flight closest to the hypothesized brine network, and each transect is numbered in increasing order approaching Sverdrup Glacier (Figs. 9–17).
Radargram, reflectivity, and specularity content for DEV3_PER0a_Y79a.
Radargram, reflectivity, and specularity content for DEV3_PER0a_Y80a.
Radargram, reflectivity, and specularity content for DEV3_PER0a_Y81a.
Radargram, reflectivity, and specularity content for DEV3_PER0a_Y82a.
Radargram, reflectivity, and specularity content for DEV3_PER0a_Y83a.
Radargram, reflectivity, and specularity content for DEV3_PER0a_Y84a.
Radargram, reflectivity, and specularity content for DEV3_PER0a_Y85a.
Radargram, reflectivity, and specularity content for DEV3_PER0a_Y86a.
Radargram, reflectivity, and specularity content for DEV3_PER0a_Y87a.
Appendix B
Modeled Canyon B geometry for each flow-perpendicular survey line. The actual canyon cross-section from the radar data is superimposed to demonstrate fit. The locations for each canyon cross-section are mapped in Fig. 6a. Flight DEV3_PER0a_Y79a is the flight closest to the hypothesized brine network, and each transect is numbered in increasing order approaching Sverdrup Glacier (Figs. 18–26).
Canyon B model geometry for DEV3_PER0a_Y79a.
Canyon B model geometry for DEV3_PER0a_Y80a.
Canyon B model geometry for DEV3_PER0a_Y81a.
Canyon B model geometry for DEV3_PER0a_Y82a.
Canyon B model geometry for DEV3_PER0a_Y83a.
Canyon B model geometry for DEV3_PER0a_Y84a.
Canyon B model geometry for DEV3_PER0a_Y85a.
Canyon B model geometry for DEV3_PER0a_Y86a.
Canyon B model geometry for DEV3_PER0a_Y87a.
Appendix C. Basal roughness approximation
Basal roughness is commonly discussed in the literature over longer length scales (hundreds of meters to km) which are easily captured by RES (Bingham and Siegert, Reference Bingham and Siegert2009; Hoffman and others, Reference Hoffman2022). When creating the bed DEM for our simulations, we incorporate roughness over distances significantly longer than the radar's along-track resolution of ~15 m explicitly, because this surface is created from the RES bed picks. However, roughness over distances near the radar wavelength, which contributes to divergent scattering of the radar signal, is not easily captured through RES. In order to approximate roughness over these shorter distances, we use a correlation length l c = 15 m several times our DEM resolution (l f = 5 m) as an approximation of small-scale roughness.
To identify the appropriate roughness magnitude, we downloaded ATLAS/ICESat-2 geolocated photons which transect a small area of exposed canyon and terrace features immediately west of DIC (Neumann and others, Reference Neumann2023). We assume this adjacent sub-aerial topography is representative of the sub-glacial topography in the transition zone.
We filter the data for observations with a minimum of ‘medium’ signal confidence as described in the user documentation (Neumann and others, Reference Neumann2023), resulting in mean along-track resolution of 9–16 cm. We sampled the top of three terrace features similar to the terrain near Canyon B, over a distance of 500–700 m. A Savitzky-Golay filter with window length 15 m was applied to each sample data set to approximate a ‘smooth’ surface at our chosen l c (S SG(x)). We then define surface roughness as the standard deviation of the difference between measured photon height S(x) and S SG(x) in Eqn (C1), where N is the number of photons in the data set.
Measured σ for our three sampled locations ranged from 0.22 m to 0.35 m. We assume the average value of our basal roughness beneath DIC in the transition zone falls within this range. Simulations with the Canyon B bed topography were run with σ bed = 0.2 m and 0.35 m (Fig. 27).
(a) Location of ATLAS/ICESat-2 photons (C-C’) with Landsat-09 imagery from July 31, 2023 as background (USGS, 2023). Canyons, brine network, and Sverdrup Glacier are labeled for reference. (b) Photon height and location of three sampled terraces. (c)–(e) show photon height S(x) and smoothed height S SG(x) for each sample, with resulting calculated roughness σ.
Open access
