A constraint upon the basal water distribution and thermal state of the Greenland Ice Sheet from radar bed echoes

Abstract

There is widespread, but often indirect, evidence that a significant fraction of the bed beneath the Greenland Ice Sheet is thawed (at or above the pressure melting point for ice). This includes the beds of major outlet glaciers and their tributaries and a large area around the NorthGRIP borehole in the ice-sheet interior. The ice-sheet-scale distribution of basal water is, however, poorly constrained by existing observations. In principle, airborne radio-echo sounding (RES) enables the detection of basal water from bed-echo reflectivity, but unambiguous mapping is limited by uncertainty in signal attenuation within the ice. Here we introduce a new, RES diagnostic for basal water that is associated with wet–dry transitions in bed material: bed-echo reflectivity variability. This technique acts as a form of edge detector and is a sufficient, but not necessary, criteria for basal water. However, the technique has the advantage of being attenuation insensitive and suited to combined analysis of over a decade of Operation IceBridge survey data.

The basal water predictions are compared with existing analyses of the basal thermal state (frozen and thawed beds) and geothermal heat flux. In addition to the outlet glaciers, we demonstrate widespread water storage in the northern and eastern interior. Notably, we observe a quasilinear “corridor” of basal water extending from NorthGRIP to Petermann Glacier that spatially correlates with elevated heat flux predicted by a recent magnetic model. Finally, with a general aim to stimulate regional- and process-specific investigations, the basal water predictions are compared with bed topography, subglacial flow paths and ice-sheet motion. The basal water distribution, and its relationship with the thermal state, provides a new constraint for numerical models.

Introduction

Basal water beneath the Greenland Ice Sheet (GrIS) influences and is influenced by, the dynamics and thermodynamics of the overlying ice. A lubricated bed is a necessary condition for basal sliding, which can be responsible for up to about 90 % of the ice surface velocity . Constraining the spatial distribution of basal water is important, therefore, for understanding the dynamic state of the overlying ice and its sensitivity to external forcing. A reliable estimate of the presence of basal water can also be used as a boundary condition or as a constraint in numerical modelling and to evaluate model performance and is, as a consequence, an attractive objective.

The spatial distribution of basal water beneath the GrIS is known to arise from an interplay of different physical processes including surface melt (e.g. ), basal melting due to geothermal heat (e.g. ), frictional and shear heating (e.g. ), and transport processes (surface, englacial and subglacial), which redistribute water (e.g. ). There are, however, limited observational constraints on the ice-sheet-scale distribution of basal water, and the relationship with other glacial and subglacial properties is therefore largely unexplored and speculative. The lack of unambiguous information about basal water arises primarily because there are only a few existing observations of subglacial lakes . By contrast, the Antarctic Ice Sheet currently has over 400 identified subglacial lakes , some of which have been found to be actively draining and recharging. Instead, there is evidence that basal water beneath the GrIS exists in smaller pools as wet sediment and as part of channelised networks .

The basal temperature distribution of the GrIS determines where basal water can exist and requires basal temperatures at or above the pressure melting point (PMP) for ice (herein “thawed”). Direct basal temperature measurements are, however, sparse. At the majority of the interior boreholes – Camp Century, Dye 3, GRIP, GISP2 and NEEM – basal temperatures indicate frozen beds , with the thawed bed at NorthGRIP being an exception . Toward the ice-sheet margins boreholes generally indicate thawed beds . Indirect methodologies for discriminating frozen and thawed beds (ice-sheet model predictions, radiostratigraphy and surface properties) were recently combined by to produce a frozen-thawed likelihood map for beds beneath the GrIS. Key predictions were that the central ice divides and west-facing slopes generally have frozen beds, the southern and western outlet glaciers have a thawed bed, and basal thaw extends east from NorthGRIP over a large fraction of the north-eastern ice sheet.

Spatially variable geothermal heat flux (GHF) influences the basal temperature distribution , hydrology and flow features in the interior of the GrIS. Notably, the onset of the Northeast Greenland Ice Stream (NEGIS) is predicted to arise from basal melting that is attributed to locally elevated GHF , which, in turn, has recently been attributed to the path of the Iceland hotspot track . As with basal temperature, the sparsity of borehole measurements limits direct inference of GHF (which is related to the vertical gradient of basal temperature). Instead, a range of geophysical techniques including tectonic , seismic and magnetic models have been used to map GHF beneath the ice sheet. These models, however, differ greatly in the predicted spatial distribution for GHF and also in the absolute values . Due to the relationship between basal melt and GHF, basal water observations can be used to further refine and cross-validate GHF models .

In principle, airborne radio-echo sounding (RES) surveys provide the information and spatial coverage to infer the presence of basal water at the ice-sheet scale. Bed-echo reflectivity is the most commonly used diagnostic for this purpose (e.g. ) and is based on the prediction that, across a range of subglacial materials, wet glacier beds have a higher reflectivity than dry or frozen beds . However, due to uncertainty and spatial variation in radar attenuation (an exponential function of temperature ) bed-echo reflectivity is subject to spatial bias and can be ambiguous when mapped over larger regions . In order to mitigate spatial bias from radar attenuation, bed-echo scattering properties – including the specularity (a measure of the angular distribution of energy and, consequently, the smoothness of the bed at a radar-scale wavelength) and the bed-echo abruptness (a waveform parameter) – have also been employed in basal water detection. Although attenuation independent, use of some bed-echo scattering properties to discriminate basal water can be prone to ambiguity and can potentially lead to false-positive detection of smoother bedrock as water .

In this study we introduce a RES diagnostic for basal water that is specifically tuned to be suitable for an ice-sheet-scale assessment of basal thaw. The RES diagnostic, which we term “bed-echo reflectivity variability”, detects wet–dry transitions in bed material and acts as a form of edge detector. The technique is demonstrated to be insensitive to radar attenuation, thus reducing the likelihood of false-positive identification of basal water at the ice-sheet scale. Moreover, it also only requires local radiometric power calibration and thus enables different Operation IceBridge RES campaigns, using different radar system settings (e.g. peak transmit power, antenna gain), to be combined and mapped. However, a limitation of the technique is that it provides sufficient (not necessary) conditions for basal water, and the detected subset of basal water corresponds to sharp horizontal gradients in the bed dielectric.

The primary geographical focus of the study is to compare the spatial relationship between predicted basal water and up-to-date analyses for the basal thermal state and GHF beneath the GrIS. We observe new predictions for basal water predominantly in the northern and eastern ice sheets, which spatially correlate with elevated GHF recently inferred by . We then compare basal water and bed topography , which enables us to identify likely subglacial flow paths and storage locations beneath the contemporary ice sheet. Finally, with a view toward improving our understanding of the influence of basal water on ice-sheet motion, we compare the relationship between basal water and ice surface speed .

(a) Ice-penetrating radar flight tracks for different CReSIS radar systems. (b) Effective coverage for good-quality radar bed echoes (corresponding to peak power 10 dB above noise floor). (c) Summary of key GrIS landmarks: temperature boreholes, major drainage basin boundaries and major regions of fast-flow identified from ice surface speed . Abbreviations in (c) correspond to Camp Century (CC), Humboldt (Hu), Petermann (Pe), Ryder (Ry), Northeast Greenland Ice Stream (NEGIS), north-western margins (NWM), Jakobshavn Isbræ (JI), Kangerlussuaq (Ka), Helheim (He) and Ikertivaq (Ik). The projection is a polar stereographic north (70${}^{\circ }$ N, 45${}^{\circ }$ W) and is used in all future plots.

[Figure omitted. See PDF]

Methods

Radio-echo sounding data and bed-echo analysis

The airborne RES data used in this study were collected by the Center for Remote Sensing of Ice Sheets (CReSIS) over the time period 2003–2014. The data were taken using a succession of radar instruments: Advanced Coherent Radar Depth Sounder (ACORDS), Multi-Channel Radar Depth Sounder (MCRDS), Multi-Channel Coherent Radar Depth Sounder (MCoRDS), Multi-Channel Coherent Radar Depth Sounder: version 2 (MCoRDS v2), mounted on three airborne platforms: P-3B Orion (P3), DHC-6 Twin Otter (TO) DC8, Douglas DC-8 (DC8) . The flight-track coverage, subdivided by the radar system, is shown in Fig. a with the seasonal breakdown: ACORDS 2003 P3 and 2005 TO; MCRDS 2006 TO, 2007 P3, 2008 TO and 2009 TO; MCoRDs 2010 P3 and 2010 DC8; MCoRDs v2 2011 TO, 2011 P3, 2012 P3, 2013 P3, 2014 P3. More details on the track lengths and data segmentation can be found in . The vast majority of the data were collected in the months March–May.

The various radar system details and signal processing steps are described in detail in previous works . The centre-frequency of the radar systems is either 150 MHz (ACORDS and MCRDS) or 195 MHz (MCoRDs and MCoRDS v2), and, after accounting for pulse-shaping and windowing, the depth-range (vertical) resolution can vary from $\sim \mathrm{4.3}$ to 20 m in ice. The along-track (horizontal) resolution also varies between field seasons and is typically $\sim \mathrm{30}$ or 60 m. The radar-echo strength profiles (CSARP, level 1B data) employ fixed fast-time gain where the receiver gain is constant for each recorded pulse, which enables consistent interpretation of bed-echo power on a season-by-season basis. Whilst transmitted power can differ between seasons, since we consider local variability, offsets between seasons do not matter, which enables the combination of inter-seasonal data. Pre-2003 CReSIS data use manual gain control and hence these seasons are not included.

The procedure used to extract bed-echo power is similar to , , and aggregates power over bed-echo fading (i.e. performs a depth-range integral). Briefly, the procedure consists of the following steps. Firstly, CReSIS ice thickness (level 2) picks were used as an initial estimate for the depth-range bin of the peak bed-echo power. Secondly, a local retracker was used to locate the true depth-range bin for peak bed-echo power. Thirdly, the power was aggregated by a discrete summation over the bed-echo envelope (both before and after the peak). The summations were truncated when the power was 10 dB or less than the peak to ensure that the integral consisted of a dominant peak associated with a dominant reflecting facet. The chosen signal-to-noise (SNR) threshold of 10 dB in this study is less strict than the $\sim \mathrm{17}$ dB threshold used in and was required to extend the effective coverage to some interior regions in southern Greenland. Additionally, in this study we did not apply along-track averaging of level 1B data that was done previously in . Finally, again, to ensure a suitable SNR, a quality control measure was imposed such that peak bed-echo power must be 10 dB above the noise floor. This results in the effective coverage shown in Fig. b, which shows good- ($\text{SNR}>\mathrm{10}$ dB) and poor- ($\text{SNR}\le \mathrm{10}$ dB) quality bed-echoes. The poor quality bed echoes include a spatially coherent coverage gap in the southern interior, high-altitude data and some marginal regions.

The rationale for use of aggregated bed-echo power (over peak power) is that it serves to reduce bed-echo power variability due to roughness and thus better enables comparison with the specular bed-echo reflectivity values that are used to infer bulk material properties . Additionally, since roughness scattering loss is frequency dependent , aggregated power serves as a pragmatic way to best combine bed-echo power measurements from the 195 and 150 MHz radar systems. This is supported by the observed $\sim$ 1 dB greater scattering loss (estimated here by the dB difference between peak and aggregated power) for the 195 MHz systems at crossover points.

Key landmarks of the GrIS that are used in the data description – temperature boreholes, drainage basin boundaries and major fast flow regions – are shown in Fig. 1c.

Bed-echo power, attenuation correction and bed-echo reflectivity

The bulk material properties of glacier beds, including the presence of basal water, can, in principle, be inferred from the reflectivity of the bed echo . The reflectivity, $\left[R\right]$, is obtained from solving the decibel form of the bed-echo power equation $$\left[P\right]=\left[S\right]-\left[G\right]+\left[R\right]-\left[L\right]-\left[B\right],$$ where $\left[P\right]$ is the bed-echo power (in this case the aggregated bed-echo power, which mitigates for scattering losses from the ice–bed interface), $\left[S\right]$ is the system performance, $\left[G\right]$ is the geometric correction, $\left[L\right]$ is the attenuation loss in ice, and $\left[B\right]$ is the loss due to birefringence (ice fabric anisotropy), and the notation $\left[X\right]=\mathrm{10}{\log }_{\mathrm{10}}X$ is assumed . Rough surface losses from transmission through the air–ice interface (e.g. ) also influence the radar power budget and are discussed in more detail in Sect. . The geometric correction for a specular reflector can be defined by $$\left[G\right]=\mathrm{2}\left[\mathrm{2}\right(s+h/\sqrt{{\mathit{ϵ}}_{\text{ice}}}\left)\right],$$ where $s$ is the aircraft height and $h$ is the ice thickness and ${\mathit{ϵ}}_{\text{ice}}=\mathrm{3.15}$ is the relative dielectric permittivity of ice to give the geometrically corrected bed-echo power $$\left[P_{\text{g}}\right]=\left[P\right]+\left[G\right],$$ (e.g. ). For the majority of the CReSIS data used (2006 TO onward) $s$ and $h$ are known and Eq. () can be applied exactly. For the 2003 P3 and 2005 TO seasons only $h$ is known and $s=\mathrm{500}$ m is assumed (approximately the mean aircraft height). This approach is justifiable since in this study we are interested in local (length scale $\sim \mathrm{5}$ km) power variation, where $s$ varies slowly. It is assumed that variation in $\left[S\right]$ and $\left[B\right]$ is also negligible (again, approximations that are strengthened by consideration of local power variation); then Eq. () reduces to $$\left[P_{\text{g}}\right]=\left[R\right]-\left[L\right].$$ Finally, setting $\left[L\right]=\mathrm{2}<N>h$ gives $$\left[P_{\text{g}}\right]=\left[R\right]-\mathrm{2}<N>h,$$ where $<N>$ (dB km${}^{-\mathrm{1}}$) is the (one-way) depth-averaged attenuation rate .

A fundamental ambiguity in bed-echo reflectivity analysis is that there are two unknowns in Eq. (): $<N>$ and $\left[R\right]$. Two approaches are typically used to determine $<N>$: (i) forward modelling using estimates of attenuation as a function of englacial temperature (e.g. ) and (ii) empirical determination using the linear regression of bed-echo power and ice thickness (e.g. ). Attenuation follows an Arrhenius (exponential) relationship with temperature and a linear dependence on the concentration of ionic impurities: primarily hydrogen (${\mathrm{H}}^{+}$) but also chlorine (${\mathrm{Cl}}^{-}$), and ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) . On an ice-sheet scale, the uncertainty when forward modelling $<N>$ is so high that it can be prohibitively challenging to accurately calibrate $\left[R\right]$ . This is due to both uncertainty in ice-sheet model temperature fields, the ionic concentrations, and the tuning of the parameters in the Arrhenius equation . Empirical determination of $<N>$ using bed-echo power is also subject to sources of potential bias. In particular, the regression methods can be ill-posed when there is rapid spatial variation in attenuation , or when there is a thickness-correlated distribution in bed-echo reflectivity .

We will later demonstrate that, unlike absolute values of $\left[R\right]$, local variability in bed-echo reflectivity is highly insensitive to modelled values of $<N>$ (Sect. ). However, despite acting as a very weak constraint, an initial estimate for ice-sheet-scale variation in $<N>$ is still necessary to calculate reflectivity variability. The estimate for $<N>$ relies on previous work by and uses the M07 Arrhenius equation , the Greenland Ice Sheet Model (GISM) temperature field from as updated in , depth-averaged ionic concentrations from the GRIP ice core , and the Greenland ice thickness data set in . The temperature field derives from a full 3-D thermomechanical simulation over several glacial-interglacial cycles and resolves the flow on a model resolution of 5 km, which causes some smoothing of the temperature field in narrow outlet glaciers near to the coast. In the calculation of $<N>$ the temperature field was firstly interpolated to a 1 km resolution, then vertically scaled using the 1 km representation of the ice thickness. The geothermal heat flux in GISM was initially taken from but further adjusted with Gaussian functions around the deep ice core sites to match observed basal temperatures. Vertical temperature profiles are within 1–2 ${}^{\circ }$C when compared to available in situ measurements.

Using this model framework, it is that predicted that $<N>$ varies by a factor $\sim \mathrm{5}$ over the GrIS, ranging from $\sim \mathrm{6}$ dB km${}^{-\mathrm{1}}$ in the colder northern interior to $\sim \mathrm{30}$ dB km${}^{-\mathrm{1}}$ toward the warmer south-western margins (refer to Fig. 5a in for a spatial plot).

Calculating bed-echo power and reflectivity variability

In this study we use simple standard deviation measures to quantify the variability of bed-echo power and reflectivity along the flight tracks, denoted by ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ and ${\mathit{\sigma }}_{\left[R\right]}$. The numerical calculation is analogous to how topographic roughness (the rms height) is calculated from bed elevation profiles , and assumes an along-track window of length 5 km evaluated at 1 km intervals (taking the standard deviation of all the points within the window). The choice of horizontal length scale was chosen for consistency with the basal thermal state mask by , which was deemed an appropriate scale for integration of radar data with thermomechanical models at the ice-sheet scale. At this 5 km length scale, high values of reflectivity variability – a diagnostic for wet–dry transitions (also, see Sect. ) – act as a form of edge detector. The rationale for using this approach (rather than a conventional edge detector) is that it does not impose that a single along-track transition is present, which is desirable when aiming to also detect multiple smaller water bodies relative to the window size.

Since this application of along-track variability is a non-standard approach to bed-echo data analysis, we now take a closer look at the statistical properties. The formula for ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ follows from the variance of Eq. () and is given by $${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}=\sqrt{{\mathit{\sigma }}_{\left[R\right]}^{\mathrm{2}}+{\mathit{\sigma }}_{\left[L\right]}^{\mathrm{2}}-\mathrm{2}{\mathit{\sigma }}_{\left[R\right],\left[L\right]}},$$ where ${\mathit{\sigma }}_{\left[L\right]}$ is the standard deviation in attenuation loss, and ${\mathit{\sigma }}_{\left[R\right],\left[L\right]}$ is the covariance of bed-echo reflectivity and attenuation loss. Using $\left[L\right]=\mathrm{2}<N>h$ and assuming $<N>$ can be approximated as constant (justifiable at the 5 km length scale that is considered); then Eq. () becomes $${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}=\sqrt{{\mathit{\sigma }}_{\left[R\right]}^{\mathrm{2}}+\mathrm{4}<N{>}^{\mathrm{2}}{\mathit{\sigma }}_h^{\mathrm{2}}-\mathrm{4}<N>{\mathit{\sigma }}_{\left[R\right],h}},$$ where ${\mathit{\sigma }}_h$ is the standard deviation of ice thickness and ${\mathit{\sigma }}_{\left[R\right],h}$ is the covariance of bed-echo reflectivity and ice thickness. Equations () and () consider the variability of log-transformed variables (i.e. the dB or additive form of the radar power equation). This differs from the variability of the linear variables (i.e. the multiplicative form of the radar power equation) expressed in log space. The first advantage to our use of log-transformed variables is that we can better separate the variability contributions into separate components. The second advantage is that we can make a clearer connection to the dB reflection amplitudes that are typically used in radioglaciology . This includes prior applications that have considered the distribution and standard deviation of dB reflection amplitudes in the context of water detection (e.g. ).

Example flight-track profiles for bed-echo power and variability, $\left[P_{\text{g}}\right]$ and ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$, bed-echo reflectivity and variability, $\left[R\right]$ and ${\mathit{\sigma }}_{\left[R\right]}$, attenuation loss and variability, $\left[L\right]$, ${\mathit{\sigma }}_{\left[L\right]}$ and ice thickness, $h$. Example (a), left panel, is from the north-central interior of the ice sheet and (b), right panel, is from the north-western margins (locations both shown in Fig. ). The threshold that is later used for water detection, ${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB, is indicated by the dashed pink line and applies to the right axis. The values for $\left[R\right]$ are relative with zero mean. The variability measures are all calculated at a 5 km length scale with 1 km posting.

[Figure omitted. See PDF]

In regions where ${\mathit{\sigma }}_{\left[R\right],\left[L\right]}\approx \mathrm{0}$ (bed-echo reflectivity has negligible covariance with attenuation loss), Eqs. () and () are approximated by $$\begin{array}{}& {\displaystyle }{\mathit{\sigma }}_{\left[P_{\text{g}}\right]}& {\displaystyle }\approx & {\displaystyle }\sqrt{{\mathit{\sigma }}_{\left[R\right]}^{\mathrm{2}}+{\mathit{\sigma }}_{\left[L\right]}^{\mathrm{2}}},& {\displaystyle }& {\displaystyle }\approx & {\displaystyle }\sqrt{{\mathit{\sigma }}_{\left[R\right]}^{\mathrm{2}}+\mathrm{4}<N{>}^{\mathrm{2}}{\mathit{\sigma }}_h^{\mathrm{2}}},\end{array}$$ and the loss component of ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ is solely modulated by ${\mathit{\sigma }}_{\left[L\right]}$ which is proportional to the product $<N>{\mathit{\sigma }}_h$. Whilst an approximation (in certain circumstances the second and third terms on the right-hand side of Eqs. () and () can be of comparable magnitude), this scenario provides an intuitive way to understand the interrelationship between ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$, ${\mathit{\sigma }}_{\left[R\right]}$ and ${\mathit{\sigma }}_{\left[L\right]}$.

Two example profiles for $\left[P_{\text{g}}\right]$, ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$, $\left[R\right]$, ${\mathit{\sigma }}_{\left[R\right]}$, $\left[L\right]$, ${\mathit{\sigma }}_{\left[L\right]}$ and $h$ are shown in Fig. . Figure a is an example from the interior of the ice sheet where ${\mathit{\sigma }}_{\left[L\right]}$ is relatively low and $h$ is thick ($\sim \mathrm{2.8}$ km). Subsequently the profiles for ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ and ${\mathit{\sigma }}_{\left[R\right]}$ are very similar in appearance, with the most notable difference at a distance of $\sim \mathrm{360}$ km where there is higher ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ due to the subglacial trough. The peaks in ${\mathit{\sigma }}_{\left[R\right]}$ are later related to wet–dry bed material transitions in Sect. . It is also important to note that ${\mathit{\sigma }}_{\left[R\right]}$ can be greater than ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ (e.g. at a distance of $\sim \mathrm{342}$ km), which is an effect that can be explained by the covariance between attenuation loss/ice thickness and bed-echo reflectivity in Eqs. () and (). Figure b is a representative example from toward the ice-sheet margins, where ${\mathit{\sigma }}_{\left[L\right]}$ is higher due to more rapid variation in $h$ (more complex bed topography) and higher values of $<N>$ (warmer ice). In this case, ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ is noticeably greater than ${\mathit{\sigma }}_{\left[R\right]}$, with the differences largely attributable to higher ${\mathit{\sigma }}_{\left[L\right]}$ as anticipated by Eq. (). The examples in Fig.  highlight that, at this length scale, higher values of ${\mathit{\sigma }}_{\left[R\right]}$ arise primarily due to a large single transition in $\left[R\right]$. However, there are also instances where multiple transitions result in local variability peaks (e.g. at distance $\sim \mathrm{380}$–385 km along Fig. a).

When calculating ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$, ${\mathit{\sigma }}_{\left[R\right]}$ and ${\mathit{\sigma }}_{\left[L\right]}$, bed-echo coverage gaps within a 5 km bin (see Fig. b) were accounted for by neglecting bins where less than half the data corresponded to good bed echoes. The effects of this filtering step are demonstrated in Fig. b, where, aligned with the deep subglacial trough at distance $\sim \mathrm{1324}$ km, there are along-track gaps in ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$, ${\mathit{\sigma }}_{\left[R\right]}$ and ${\mathit{\sigma }}_{\left[L\right]}$.

It is important to clarify the difference between the use of bed-echo power and reflectivity variability in this study from previous radioglaciology studies . These studies focused on the statistical distribution of the peak echo power as a result of phase modulation by interfacial roughness. By contrast, in this study we suppress roughness effects by performing a depth integral for power over the echo envelope (Sect. ). We are therefore able to focus on power variability that is a result of along-track changes in the bed dielectric.

Spatial distributions for (a) bed-echo power variability ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$, (b) attenuation loss variability ${\mathit{\sigma }}_{\left[L\right]}$, (c) bed-echo reflectivity variability ${\mathit{\sigma }}_{\left[R\right]}$ and (d) ice thickness, $h$. Zoomed-in panels with flight-track data at true buffer size (5 km) are shown for the north-central ice sheet (pink bounding box, containing profile in Fig. a) and north-western margins (black bounding box, containing profile in Fig. b). The profiles are indicated in bold green in the ice thickness zoomed-in panels. In panels (a–c) higher variability data is stacked on top of lower variability data, which acts to emphasise higher variability. The zoomed-in panels all have the same scale ($\times \mathrm{8}$ the resolution of the ice-sheet-scale plots).

[Figure omitted. See PDF]

Distributions for bed-echo power and reflectivity variability

Spatial distributions for the variability measures: ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$, ${\mathit{\sigma }}_{\left[L\right]}$, ${\mathit{\sigma }}_{\left[R\right]}$ are shown in Fig. a, b, c along with ice thickness in Fig. d. In general, ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ has a strong ice thickness dependence and increases toward the margins where ice is thinner. The attenuation correction, which primarily acts to reduce the component of ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ that is attributable to ${\mathit{\sigma }}_{\left[L\right]}$, results in a more uniform ice-sheet-scale distribution of ${\mathit{\sigma }}_{\left[R\right]}$ than ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$. Notably, there are localised patches of higher ${\mathit{\sigma }}_{\left[R\right]}$ present in both marginal and interior regions (which are later attributed to the presence of basal water). The ice-sheet-scale trends in ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ and ${\mathit{\sigma }}_{\left[L\right]}=\mathrm{2}<N>{\mathit{\sigma }}_h$ may be related to spatial variation in $<N>$ and bed roughness (which correlates with ${\mathit{\sigma }}_h$).

Frequency distributions for (a) bed-echo power variability, ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ (corresponding to Fig. a); (b) bed-echo reflectivity variability, ${\mathit{\sigma }}_{\left[R\right]}$ (corresponding to Fig. c). Later in the study ${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB is used as the threshold criteria for diagnosing the presence of basal water.

[Figure omitted. See PDF]

The two zoomed-in regions in Fig.  include the flight-track profiles in Fig. . (north-central ice sheet, pink bounding box; north-western margins, black bounding box). These examples serve to further illustrate the spatial interrelationship between ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$, ${\mathit{\sigma }}_{\left[R\right]}$ and ${\mathit{\sigma }}_{\left[L\right]}$ in a typical interior region with lower ${\mathit{\sigma }}_{\left[L\right]}$ and a typical marginal region with higher ${\mathit{\sigma }}_{\left[L\right]}$. Its is clear that the interior example has very similar spatial distributions for ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ and ${\mathit{\sigma }}_{\left[R\right]}$, whereas the marginal example has higher ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ associated with the higher ${\mathit{\sigma }}_{\left[L\right]}$ that occurs in the subglacial troughs and more complex topography toward the edge of the ice sheet. The marginal example also demonstrates that the power variability associated with the subglacial troughs is largely removed for ${\mathit{\sigma }}_{\left[R\right]}$.

The corresponding frequency distributions for ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ and ${\mathit{\sigma }}_{\left[R\right]}$ are shown in Fig. . Both demonstrate a strong positive-skew, with a long-tail extending to higher values. The mean and standard deviation for ${\mathit{\sigma }}_{\left[P_{\text{g}}\right]}$ is greater than ${\mathit{\sigma }}_{\left[R\right]}$. This is consistent with the commonly observed result that making an attenuation correction to $\left[P_{\text{g}}\right]$ acts to reduce the overall decibel range for $\left[R\right]$ (e.g. ), hence more closely resembling the predicted dB range for bed materials .

A quantification and discussion of crossover statistics for ${\mathit{\sigma }}_{\left[R\right]}$ is given in Appendix A.

Dielectric and reflective properties of subglacial materials based on a compilation of past values by , and . The bulk values take into account typical ranges of saturation and porosity for the dielectric mixing of water and ice with the background material. The relative dielectric permittivity of ice is 3.15, which means that dry (just the background dielectric) or frozen material (a mixture of the background dielectric with ice) produces a similar range for $\left[R\right]$.

Bed material	Permittivity $\mathit{ϵ}$	Reflectivity $\left[R\right]$ (dB)
Groundwater	80	$-\mathrm{2}$
Wet till	10 to 30	$-\mathrm{11}$ to $-\mathrm{6}$
Wet sandstone	5 to 10	$-\mathrm{19}$ to $-\mathrm{11}$
Dry/frozen granite	5	$-\mathrm{19}$
Dry/frozen limestone	4 to 7	$-\mathrm{26}$ to $-\mathrm{14}$
Dry/frozen till	2 to 6	negligible to $-\mathrm{19}$
Dry/frozen sandstone	2 to 3	$-\mathrm{37}$ to $-\mathrm{16}$

Interpretation of reflectivity variability as a sufficient diagnostic for basal water

Radar bed-echo reflectivity depends on the dielectric contrast between glacier ice and bed material. For a specular, nadir reflection, the Fresnel power reflection coefficient is given by $$\left[R\right]=\mathrm{10}{\log }_{\mathrm{10}}{\left|{\displaystyle \frac{\sqrt{{\mathit{ϵ}}_{\text{bed}}^{\sim }}-\sqrt{{\mathit{ϵ}}_{\text{ice}}^{\sim }}}{\sqrt{{\mathit{ϵ}}_{\text{bed}}^{\sim }}+\sqrt{{\mathit{ϵ}}_{\text{ice}}^{\sim }}}}\right|}^{\mathrm{2}},$$ where ${\mathit{ϵ}}_{\text{ice}}^{\sim }$ and ${\mathit{ϵ}}_{\text{bed}}^{\sim }$ are the complex dielectric permittivity of the glacier ice and bed material respectively. The relative (real) part of the permittivity, ${\mathit{ϵ}}_{\text{bed}}$, is the primary control on $\left[R\right]$. A summary of dielectric and reflective properties of glacier bed materials at typical ice-penetrating radar frequencies is given in Table 1 and is collated from , and . The permittivity and reflectivity range for each material arises due to sub-wavelength dielectric mixing between either ice or water and the bed material, and takes into account typical saturation and porosity values . In general, lower values of ${\mathit{ϵ}}_{\text{bed}}$ and $\left[R\right]$ occur for dry or frozen bed materials (approximately ${\mathit{ϵ}}_{\text{bed}}<\mathrm{7}$ and $\left[R\right]<-\mathrm{14}$ dB), whilst higher values occur for wet bed materials (approximately ${\mathit{ϵ}}_{\text{bed}}>\mathrm{7}$ and $\left[R\right]>-\mathrm{14}$ dB). Dielectric mixing between bed materials can also occur at the length scale of the Fresnel zone ($\sim \mathrm{100}$ m), which results in further averaging of the observed reflectivity .

Interpretation of bed-echo reflectivity variability, ${\mathit{\sigma }}_{\left[R\right]}$, as a diagnostic for basal water. (a, b) Schematics of the two-state dielectric model for single and multiple along-track transitions in dry to wet bed material. Both scenarios are identically parameterised by the wet$/$dry mixing ratio $f$ (visually, the fraction of blue to yellow) and wet–dry reflectivity difference, $\mathrm{\Delta }\left[R\right]=\left[R_{\text{wet}}\right]-\left[R_{\text{dry}}\right]$. (c) Phase-space plot for ${\mathit{\sigma }}_{\left[R\right]}$ as a function of $f$ and $\mathrm{\Delta }\left[R\right]$. ${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB is used as a threshold for positive discrimination of basal water (corresponding to green, red and yellow regions).

[Figure omitted. See PDF]

Due to the range of possible bed materials at the ice-sheet scale it is not possible to formulate a unique dielectric model for diagnosing water from ${\mathit{\sigma }}_{\left[R\right]}$. A simple two-state dielectric model, does, however, enable us to physically motivate the water diagnostic in terms of dielectric properties (Fig. ). The model assumes that the along-track sample window is comprised of two different bed materials: the dry background bed material with permittivity ${\mathit{ϵ}}_{\text{dry}}$ and reflectivity $\left[R_{\text{dry}}\right]$, and the wet material with permittivity ${\mathit{ϵ}}_{\text{wet}}$ and reflectivity $\left[R_{\text{wet}}\right]$. For simplicity, it is assumed that each along-track measurement is in one of the wet or dry states, with the wet–dry mixing ratio parameterised by $f$. In this formulation, a single body of wet material or multiple smaller bodies of wet material have the same formula for the reflectivity variability given by $${\mathit{\sigma }}_{\left[R\right]}=\mathrm{\Delta }\left[R\right]\sqrt{f^{\mathrm{2}}(\mathrm{1}-f)+(\mathrm{1}-f{)}^{\mathrm{2}}f},$$ where $\mathrm{\Delta }\left[R\right]=\left[R_{\text{wet}}\right]-\left[R_{\text{dry}}\right]$ is the reflectivity difference between wet and dry beds. Equation () is derived by considering the weighted variance for two discrete random variables and does not account for non-linear variations due to variable scattering coherence. A phase-space plot for ${\mathit{\sigma }}_{\left[R\right]}(f,\mathrm{\Delta }[R\left]\right)$ is shown in Fig. c, and shows that, for fixed $\mathrm{\Delta }\left[R\right]$, ${\mathit{\sigma }}_{\left[R\right]}$ is maximised when $f=\mathrm{0.5}$ (i.e. an even mixing of wet and dry materials).

Past diagnosis of basal water typically associates the upper tail of the reflectivity distribution with water, prescribing a threshold above which the bed is interpreted as wet (e.g. ). In this study, a similar thresholding approach is applied to the distribution of ${\mathit{\sigma }}_{\left[R\right]}$ (Fig. b). The threshold choice for basal water (${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB) corresponds to the region greater than the ${\mathit{\sigma }}_{\left[R\right]}=\mathrm{6}$ dB contour in Fig. c and requires a minimum wet–dry reflectivity difference of $\mathrm{\Delta }\left[R\right]>\mathrm{12}$ dB. In general, $\mathrm{\Delta }\left[R\right]>\mathrm{12}$ dB is only possible for a mixture of wet and dry (or frozen) bed materials (Table 1). For example, an even mixing of groundwater and dry granite ($f=\mathrm{0.5}$, $\mathrm{\Delta }\left[R\right]=\mathrm{17}$ dB) has ${\mathit{\sigma }}_{\left[R\right]}=\mathrm{8.5}$ dB. The contours in Fig. c demonstrate that small perturbations to even mixing ($f\ne \mathrm{0.5}$) produce similar ${\mathit{\sigma }}_{\left[R\right]}$, and hence that water detection is insensitive to the discretisation of the along-track sample window (Sect. ). Overall, the threshold choice (${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB) is fairly conservative and is deliberately intended to reduce false-positive detection of basal water (at the expense of reduced overall detection). Finally, as discussed in Sect. , the 6 dB variability threshold applies to the distribution of log-transformed reflectivity.

The bed-echo power aggregation in Sect.  partially mitigates for roughness-induced scattering loss and the along-track power variability associated with this. Supporting evidence is that the crossover analysis for ${\mathit{\sigma }}_{\left[R\right]}$ (see Appendix A) demonstrates there to be no significant bias for the 150 MHz radar systems (ACORDS and MCRDS) versus MCORDS v2 (the 195 MHz radar system used as a benchmark). Additionally, whilst small-scale roughness transitions (transitions from specular to diffuse scattering) will often correlate with wet–dry transitions, this scenario will act to amplify the useful signal component with the ${\mathit{\sigma }}_{\left[R\right]}$ value increasing.

Table 1 also indicates that, in exceptional circumstances, $\mathrm{\Delta }\left[R\right]>\mathrm{12}$ dB could be generated in frozen/dry regions that partially contain sandstone or till that is close to matching the permittivity of ice. However, if present, these regions are likely to have indistinct bed echoes and will not be included in the effective coverage in Fig. b.

Basal water distribution and robustness to perturbations in the attenuation rate estimate, $<N>$. The original predictions (${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB) are represented by all three colours. Persistent water predictions (${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB for $\pm \mathrm{20}$ % and $\pm \mathrm{50}$ % perturbations to $<N>$) are indicated by the subset of green and red points, and the subset of red points respectively. The subset of red and green points is used in the rest of the paper.

[Figure omitted. See PDF]

Basal water distribution and robustness to attenuation model bias

The initial basal water predictions (${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB, pre- sensitivity analysis) are shown in Fig.  (red, blue and green data), and correspond to $\sim \mathrm{3.5}$ % of bins containing predominantly good-quality bed echoes (Fig. b). A full geographic analysis of the spatial distribution is performed in Sect. . To demonstrate the robustness of the predictions, we performed a sensitivity analysis with respect to the modelled attenuation correction $<N>$ (Sect. ) The analysis considered a series of increasingly large (uniform, multiplicative) perturbations to $<N>$ and then tested whether ${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB also held for the perturbed model. Examples of persistent water predictions for $\pm \mathrm{20}$ % (red and green data) and $\pm \mathrm{50}$ % (red data) perturbations to $<N>$ are indicated. As the perturbation size increases this results in a slight decrease in the overall percentage of water predictions (corresponding to $\sim \mathrm{2.6}$ % and 2.1 % of the along-track bins for $\pm \mathrm{20}$ % and $\pm \mathrm{50}$ % respectively).

The sensitivity analysis tests the robustness of the water predictions to a number of different physical scenarios. Firstly, there is inherent bias in the Arrhenius equation parameters. For example, an empirical correction similar to the uniform perturbation considered in Fig.  was proposed by to model unaccounted frequency dependence in the electrical conductivity. Secondly, there is bias in the model temperature field ($<N>$ is approximately equivalent to depth-averaged temperature). Thirdly, the bias is due to assumed ionic concentration values. It is hard to formally quantify the possible range of these uncertainties but, based on solution variability for $<N>$ using ice-sheet model temperature fields , $\pm \mathrm{20}$ % is a reasonable estimate for temperature-related uncertainty. Subsequently, in the comparison with other data sets in Sect.  the subset of red and green points in Fig.  is used. Inherent bias in the Arrhenius equation parameters could be significantly higher than temperature uncertainty . However, since the spatial structure for the basal water distribution under the $\pm \mathrm{50}$ % perturbation is largely preserved, this is unlikely to significantly alter the conclusions that are drawn. Additionally, whilst ${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB is used as the threshold, Fig. c demonstrates that a less conservative threshold (${\mathit{\sigma }}_{\left[R\right]}>\mathrm{4}$ dB) retains the majority of contiguous regions with no water predicted (e.g. the interior of southern Greenland).

It is important to emphasise the robustness of ${\mathit{\sigma }}_{\left[R\right]}$ with respect to uncertainty and model bias in $<N>$ (particularly compared with bed-echo reflectivity, $\left[R\right]$). An analogous sensitivity analysis by demonstrated that systematic over- and underestimates in $<N>$ lead to pronounced ice-thickness-correlated biases in the distribution for $\left[R\right]$ in northern Greenland (Fig. B1 in ).

Results

The basal water distribution is now compared with existing analyses for the basal thermal state (Sect. ), geothermal heat flux (GHF) (Sect. ), bed topography and subglacial flow paths (Sect. ) and ice surface speed (Sect. ). The basal water predictions in the results are always indicated by red circles and correspond to the set of red and green points in the sensitivity analysis, Fig. . In regional zoomed-in panels the circles are fixed to be 5 km in diameter (a true representation of the along-track window size and the effective resolution of the radar method). In ice-sheet-scale plots the buffer size of the water predictions are increased for visualisation purposes. The radar flight tracks represent where there are good bed echoes (Fig. b) and hence indicate the effective coverage.

In interpreting the maps it is important to emphasise that the basal water predictions in this study correspond to a subset of flight-track data where basal water is present. Specifically, they correspond to where there are rapid spatial transitions and gradients in the bed dielectric (i.e. small, finite, water bodies or the boundaries of larger water bodies). The predictions therefore act as a sufficient constraint on the distribution of basal water rather than being a fully comprehensive flight-track map for the water extent. Additionally, since the vast majority of the radar measurements were collected before the onset of summer surface melt, to a first approximation the basal water predictions correspond to the winter storage configuration.

Comparison between basal water distribution and basal thermal state synthesis by . (a) Ice-sheet scale. Major regions of disagreement (water in likely frozen regions) are highlighted in the zoomed-in panels. (b) Camp Century. (c) Far north. (d) North-central ice sheet. (e) East of GRIP. (f) Around Kangerlussuaq. The zoomed-in panels all have the same scale ($\times \mathrm{5}$ the resolution of a).

[Figure omitted. See PDF]

Comparison between basal water distribution and geothermal heat flux (GHF) models. (a) Seismic GHF model by . (b) Magnetic GHF model by using satellite data. (c) Magnetic GHF model by derived from spectral methods using airborne data. The colour bar scale is the same in all panels and is truncated to emphasise the spatial variation in panel (c).

[Figure omitted. See PDF]

Comparison between basal water distribution and basal thermal state synthesis

In Fig. a the basal water predictions are underlain by the basal thermal state synthesis (frozen/thawed likelihood) map by . The synthesis employed four independent methods: (i) assessment of thermomechanical model temperature fields, (ii) basal melting inferred from radiostratigraphy, (iii) basal motion inferred from surface velocity, and (iv) basal motion inferred from surface texture. The four methods were then equally weighted, leading to a likelihood map for frozen beds, thawed beds and uncertain regions. Importantly, the prediction did not utilise radar bed-echo data and is therefore independent of our basal water predictions.

The reflectivity variability water diagnostic enables a positive discrimination of basal thaw, since ${\mathit{\sigma }}_{\left[R\right]}>\mathrm{6}$ dB is deemed a sufficient (but not necessary) criteria for basal water. Positive discrimination of frozen regions is not, however, possible. This is because low-reflectivity variability (${\mathit{\sigma }}_{\left[R\right]}<\mathrm{6}$ dB) could correspond to many different scenarios: a frozen region, a drier region at or above the PMP or a wet region that is smoothly varying in bed-echo reflectivity. Since basal water enables a positive discrimination of thaw, red circles in likely thawed (pink) regions indicate agreement and red circles in likely frozen (blue) regions indicate disagreement with the basal thermal state synthesis. Absence of basal water in likely frozen regions is an indicator of general consistency between the two methods.

Comparison between basal water distribution, bed topography and major subglacial flow paths (blue lines). (a) Ice-sheet scale. (b) Petermann catchment. (c) North-western margins. (d) North-central ice sheet. To improve clarity the radar flight tracks are removed from (a) and a hillshade is applied to the bed topography. The zoomed-in panels are all have the same scale ($\times \mathrm{4}$ the resolution of a).

[Figure omitted. See PDF]

There is general agreement (water in predicted thawed regions) for the beds of major outlet glaciers and their upstream regions. This includes Helheim, Kangerlussuaq, Jakobshavn and the other fast-flowing regions identified in Fig. c. There is also general agreement between basal water and the extent of predicted thaw in the NEGIS drainage basin. Major regions of disagreement (water in predicted frozen regions) are highlighted in the zoomed-in panels, Fig. b–f. The most obvious disagreement is the quasilinear “corridor” of basal water in the north-central ice sheet (Fig. d). This feature tracks close to the central ice divides and extends from the NorthGRIP region in the south toward Petermann Glacier in the north. There are also noticeable areas of disagreement to the north and east of the Camp Century borehole (Fig. b), in the far north (Fig. c), to the east of GRIP (Fig. e), and around Kangerlussuaq (Fig. f). There is also an absence of water in many predicted frozen regions, indicating consistency. This includes parts of the southern interior, north of the NEGIS drainage basin, and the majority of the interior region between the Camp Century and NEEM boreholes.

Comparison between basal water distribution and ice surface speed (logarithmic-scale). (a) Ice-sheet scale. (b) Humboldt, Petermann and Ryder. (c) Kangerlussuaq and region to south. (d) Helheim (north of panel) and Ikertivaq (south-west of panel). The zoomed-in panels are all have the same scale ($\times \mathrm{4}$ the resolution of a).

[Figure omitted. See PDF]

Comparison between basal water distribution and geothermal heat flux models

The basal temperature of glacier ice is governed by surface temperature, GHF, strain heating from internal deformation, frictional heating, and diffusive and advective heat transport (e.g. ). In the interior of the ice sheet, close to the ice divides, GHF, vertical advection, and diffusion are the dominant processes which influence basal temperature. In this scenario, the thermodynamic (temperature) equation can be approximated by the classical Robin model which predicts that basal melting occurs when GHF is above a certain threshold. For typical values of ice thickness and surface accumulation rate, which control the rate of vertical heat advection, the minimum GHF forcing for melt is anticipated to be $\sim \mathrm{55}$–70 mW m${}^{-\mathrm{2}}$ in the interior of the ice sheet . In the water–GHF comparison we therefore define elevated GHF (i.e. likely to produce basal melt) as $>\mathrm{60}$ mW m${}^{-\mathrm{2}}$. This definition is also informed by the lower range of values (37–50 mW m${}^{-\mathrm{2}}$) that are typically associated with non-altered ancient continental crust .

In Fig.  the basal water predictions are underlain by three different GHF models: the seismic model by and two models derived from magnetic anomalies by and . The GHF model by is based on the correlation between a 3-D tomographic model of the crust and mantle temperature. The GHF models by and are based on a thermal model of the lithosphere with the lower boundary defined by the Curie depth, which is determined from magnetic anomalies. further describes this approach and the additional spectral processing method used to produce Fig. c. An older tectonic GHF model by is not considered and a spatial plot for the GrIS can be found in along with a discussion of the caveats of the different types of model. A summary of local GHF estimates using borehole temperature profiles and thermomechanical model inversions are provided by , and demonstrate general consistency between Fig. c and local estimates at GRIP, NEEM, NorthGRIP and Camp Century. Local estimates of GHF at Dye 3 ($\sim \mathrm{20}$–25 mW m${}^{-\mathrm{2}}$) are significantly lower than all three GHF models.

In interpreting Fig. , we limit the comparison to the ice sheet interior where the spatial correlation between GHF and basal water should be strongest. The model by , Fig. a, predicts low GHF ($<\mathrm{60}$ mW m${}^{-\mathrm{2}}$) over the vast majority of the central and northern interior. There is therefore no correlation between elevated GHF and basal water. The model by , Fig. b, predicts elevated GHF around GRIP and the southern and eastern boundaries of the NEGIS basin, and basal water is also present in this region. There is, however, no correlation between elevated GHF and basal water along the ice divides north of NorthGRIP. The model by , Fig. c, exhibits strong overall spatial correlation between basal water and elevated GHF in the interior of the northern ice sheet. Notably, there is a striking correlation between elevated GHF and the quasilinear corridor of basal water that extends from NorthGRIP toward Petermann Glacier. All three models predict regions of elevated GHF in the southern interior including the Dye 3 region. However, there is only isolated radar evidence for basal water.

In the comparison between the flight-track water predictions and GHF distributions in Fig.  it is important to bear in mind that the GHF distributions are evaluated at a lower spatial resolution. For example, the resolution of the GHF distribution by is a consequence of the spectral method (window size and overlap), which has an effective resolution of $\sim \mathrm{75}$ km.

Comparison between basal water distribution, bed topography and subglacial flow paths

In Fig.  the basal water predictions are underlain by the most recent Greenland bed topography digital elevation model (DEM) . To motivate further discussion about water storage locations and hydrological connectivity, a predicted subglacial flow path network is also included. The network structure is governed by gradients in the hydraulic pressure potential , which was calculated using the bed elevation and ice thickness surfaces at a grid cell resolution of 600 m (derived from ). The flow-routing algorithm was implemented in ArcGIS using the inbuilt flow accumulation tool on a hydraulic potential surface that had been processed for sink removal using the method of . Likely hydrological flow paths were identified by excluding flow paths where 50 or fewer neighbouring cells cumulatively contribute to a given location.

Figure  demonstrates that the vast majority of the basal water predictions are well aligned with predicted subglacial flow paths. This alignment is most visually pronounced toward the margins and zoomed-in panels are shown for the Petermann catchment in Fig. b and north-western margins in Fig. c. Figure b also demonstrates that basal water is present along sections of the “mega-canyon” feature identified by – for example, north-west of the intersection (80${}^{\circ }$ N, 50${}^{\circ }$ W). In the interior of the ice sheet, where the horizontal gradients in ice thickness are small, local minima in the hydraulic potential surface should correlate with topographic depressions. The water storage locations in the interior generally conform to this behaviour (Fig. d).

Comparison between basal water distribution and ice surface speed

In Fig.  the basal water predictions are underlain by ice surface speed , which is based on a temporal average from 1 December 1995 to 31 October 2015. The ice surface speed is determined using interferometric synthetic aperture radar (InSAR) as described in . Whilst there is a complex overall relationship between basal water and ice velocity, there are some clear spatial patterns. Notably, in the topographically less constrained northern and western outlet glaciers, basal water is often concentrated in the fast-flow onset regions and tributaries whilst it is absent from the main trunks. This behaviour is particularly evident for the Petermann Glacier catchment (Fig. b). In the topographically more constrained south-eastern outlet glaciers, there is widespread evidence for basal water storage in both the fast-flowing glacial troughs and upstream regions. This includes both the Kangerlussuaq catchment and the tight network of subglacial troughs to the south (Fig. c), and the Helheim catchment (Fig. d).

In the interior of the ice sheet basal water is predicted near to the head of the NEGIS ice stream. However, basal water is also predicted in some of the slowest-flowing regions of the ice sheet interior, notably, close to the central ice divides between NorthGRIP and Petermann and north-east of GRIP.

Discussion

Basal water, basal thermal state and temperature

The basal water distribution in this study and the basal thermal state synthesis by represent two independent approaches to predict where the bed beneath the GrIS is thawed. There is greatest agreement (basal water in likely thawed regions identified by ) toward the ice margins where ice surface speed is generally higher. The most noticeable regions of disagreement (basal water in likely frozen regions identified by ) all occur where the ice surface speed is low. This includes the north-central ice divide (Fig. d), the region east of GRIP (Fig. e), and the region west of Kangerlussuaq (Fig. f). The regions of agreement and disagreement are, perhaps, unsurprising, since three of the four methods employed by – ice velocity, surface texture and radiostratigraphy – associate basal thaw with present (or past) ice sheet motion. In general, a thawed bed is a necessary (but not sufficient) condition for appreciable basal motion, and there is likely to be a subset of thawed regions where basal motion is negligible. This subset naturally incorporates water/thaw near the ice divides (since driving stress is low) and in the eastern ice sheet (since ice flow is topographically constrained).

Another key difference between the water predictions in this study and the thaw predictions by is that their study employed techniques better tuned to identify continuous regions of basal thaw, whereas the basal water predictions are localised. This provides another means to reconcile regions of disagreement, since in some instances the basal water predictions may correspond to localised patches above the PMP in an otherwise frozen region. A final explanation for the discrepancies is that the model temperature fields included in the basal thermal state synthesis were often tuned around knowledge of GHF at the time (i.e. ).

There is no evidence for basal water at the location of the temperature boreholes, which, based on the resolution of our method, corresponds to within 5 km. Since high-reflectivity variability is not necessary for thaw, this is consistent with both frozen and thawed borehole temperatures. Water is, however, observed fairly close to two frozen boreholes: $\sim \mathrm{10}$ km south of GRIP and $\sim \mathrm{7}$ km north-east of Camp Century (Fig. ). At GRIP this is less surprising, since the basal temperature is 6 degrees below the PMP and GHF is likely to be elevated in this region (Fig. ). The basal water predictions near Camp Century are more surprising, since in the late 1960s basal temperatures were measured to be 11.8 degrees below the PMP . One possible explanation, which was recently invoked to explain the presence of a lake less than 10 km from the South Pole (where the bed is frozen), is that the basal water is yet to reach thermal equilibrium . Another possible explanation is that the presence of hypersaline water could result in a depression of the PMP. This situation arises at Lake Vida in East Antarctica (where liquid water exists at $-\mathrm{13}$ ${}^{\circ }$C) and at the Devon Ice Cap in the Canadian Arctic .

Basal water and geothermal heat flux

The comparison between basal water and the different GHF models in Fig.  demonstrates greatest consistency with the distribution by . Notably there is a pronounced spatial correlation between elevated GHF and the new predictions of basal thaw in the northern ice sheet. A recent machine-learning-derived map for GHF beneath Greenland by is also consistent with there being extensive basal thaw in this region. However, establishing definitive attribution of regions of the basal melt to GHF forcing (rather than frictional and strain heating, low advection from colder ice above, and surface meltwater storage) will require integration with thermomechanical ice-sheet models. The basal water predictions could also be used as a constraint in a wide variety of other numerical modelling contexts. For example, experiments with 3-D models to reconstruct the full ice temperature history over the last glacial cycle(s) can constrain the minimum GHF required to produce basal melting at the predicted basal water locations . Other studies include investigating the sensitivity of ice-sheet dynamics to the thermal boundary condition or basal lubrication , and thermal models of the underlying lithosphere .

Recent analyses imply that much of the spatial variation in GHF beneath the northern GrIS can be explained by Greenland's passage over the Iceland mantle plume between roughly 35 and 80 million years ago . The magnetic GHF map in Fig. c, alongside gravity data (Bouger anomalies), was recently used by to infer the most likely passage of the hotspot track. The most likely predicted path (corresponding to going forwards in geological time) follows the quasilinear region of elevated GHF in Fig. c from Petermann Glacier to NorthGRIP and follows a path previously anticipated by . The spatial correlation between elevated GHF and the quasilinear basal water corridor provides an additional source of evidence for the predicted path.

Basal water, bed topography and subglacial flow paths

There is growing evidence that much of the present day subglacial flow path network beneath the GrIS is palaeofluvial in origin. This includes the dendritic flow path networks in the Jakobshavn and Humboldt catchments , along with the prominent mega-canyon feature which extends from the NorthGRIP region in the south to the Petermann Glacier in the north . The comparison between the predicted flow paths and basal water in Fig.  enables a revised assessment of the hydrological flow paths that are likely to be utilised in the contemporary ice sheet. For example, the flow routing analysis demonstrates that basal water originating in the Petermann catchment is likely to route through sections of the canyon toward the ice-sheet margins. Fig. b supports this hypothesis, since the there is evidence for basal water along the majority of the canyon. However, it is important to stress that more rigorously assessing hydrological connectivity will require incorporation of DEM uncertainty when performing the flow routing (e.g. ) and use of a coupled hydrological ice flow model (e.g. ).

Basal water and ice-sheet motion

Both observational (e.g. ) and theoretical studies (e.g. ) point toward a complex spatio-temporal relationship between basal water and ice surface speed in fast-flowing regions of the ice sheet. This ultimately depends on the details of how the subglacial drainage system responds to surface meltwater. It is therefore essential to re-emphasise that the basal water predictions generally correspond to the winter storage (pre-surface melt) configuration and may, therefore, be of limited utility in understanding spatio-temporal patterns related to ice dynamics.

In addition to basal water and temperature, spatial variation in the underlying geology and lithology of the GrIS (notably, presence or absence of deformable sediment) will also influence ice-sheet motion. It is widely anticipated that much of the interior of the ice sheet is underlain by hard pre-Cambrian rocks, with more limited sedimentary deposits toward the margins and in the NEGIS drainage basin . It is therefore entirely plausible that much of the basal water predicted in the interior lies on a hard undeformable bed (particularly in the context of the igneous rock that would be associated with the geological remnants of the Iceland hotspot track) and therefore experiences little motion due to bed deformation.

Comparison with past RES analyses of basal water and disrupted radiostratigraphy in Greenland

Despite acknowledged calibration issues due to both variable radar system performance and spatial variation in attenuation, the bed-echo reflectivity analysis of 1990s PARCA RES data by anticipated some of the water predictions in this study. This includes prior basal water predictions in the NEGIS onset region and the upstream areas of the Kangerlussuaq, Petermann and Humboldt glaciers.

There is mixed agreement between the basal water predictions in this study and those from , who performed joint bed-echo reflectivity and scattering analysis of the 1990s PARCA data. In general, better agreement with our results occurs in smoother topographic regions in the ice-sheet interior, such as close to the NorthGRIP borehole. Since the effects of spatial bias due to attenuation uncertainty are lower in the interior of the ice sheet, this is where bed-echo reflectivity as a water diagnostic should be more robust. Additionally, the water detection method proposed by will generally not be able to discriminate water in many outlet glaciers and tributaries including Petermann and the north-western margins . This is because these regions tend to exhibit a diffuse scattering signature (associated with fine-scale roughness), whereas the method proposed by is specifically tuned to detect water bodies that exhibit a spatially continuous (near-) specular scattering signature. By contrast, comparison between the water predictions in this study and the radar-derived bed roughness maps in and demonstrate a lack of modulation by bed roughness, with basal water present in rougher marginal regions and a generally smoother ice-sheet interior.

Basal units of disrupted radiostratigraphy are widely present in Greenland RES data sets . The features have been attributed to a supercooling freeze-on process , stick-slip mechanisms, and the rheological anisotropy of ice . In the fast-flow initiation regions of some outlet glaciers (e.g. Petermann and the northern tributaries of NEGIS), these basal units are closely aligned with the basal water predictions.

Limitations of bed-echo reflectivity variability as a RES technique to detect basal water

Bed-echo reflectivity variability provides a practical way to automate the detection of a subset of basal water with high confidence at the ice-sheet scale. In particular, as a form of edge detector, the technique is well-tuned to detect finite water bodies with sharp horizontal gradients in water content. These attributes are thought likely to be common to basal water in the (likely hard-bedded) interior of the ice sheet. It is, however, important to note that the approach will fail to identify basal water with a homogeneous dielectric and reflective character. This includes the centre of large subglacial lakes (based on the resolution of our method lakes greater than 5 km in horizontal extent) and regions of more uniformly saturated subglacial till. Since all identified subglacial lakes in Greenland are $<\mathrm{5}$ km in horizontal extent we believe that the former scenario is likely to be rare. However, extensive regions of saturated till that evade detection are likely to be present, particularly beneath larger outlet glaciers. This interpretation is supported by comparison with the bed-echo reflectivity and basal water maps of the Petermann catchment in . Specifically, there is a good agreement between the two water maps in the interior and fast-flow initiation region, but this study fails to predict the basal water (likely to be wet sediment) in the main trunk of the outlet glacier. Therefore, if we are focusing on the catchment-scale subglacial hydrology of Greenland outlet glaciers or other glaciologically similar regions of Antarctica, a suite of existing RES techniques to detect and characterise basal water (e.g. ) is better suited. Finally, it is important to add that part of the subglacial water budget is likely to comprise groundwater , which would be practically undetectable by RES.

As is generally the case in RES analysis, certain simplifications were made in this study when interpreting bed-returned power. Notably, we did not account for (i) power modulation due to birefringent propagation (ice fabric anisotropy, e.g. ), (ii) rough surface transmission (e.g. ) or (iii) scattering by near-surface water. The first of these mechanisms could potentially influence power variability near the ice divides, since abrupt fabric transitions can be present in these regions . However, the zoomed-in panels in these regions (e.g. Figs. b, d, d) indicate that the water predictions occur over a multiple range of flight-track orientations relative to the ice divide (which should correlate with orientation of the principle dielectric axes). We therefore can discount ice fabric having a dominant influence on the results. The second and third of these mechanisms will influence power variability primarily in faster-flowing regions toward the margins, as these regions have higher surface roughness (e.g. ) and surface melt (e.g. ). The degree of surface-induced power variability depends on both the surface permittivity and the roughness regime (relative to the radar wavelength) . It is, however, important to note that we are likely to get the majority of false positives (elevated power variability due to surface modulation) in regions where the bed is predicted to be thawed with high certainty by . In turn, the central results in this study – the new regions of basal water and thaw identified in Fig.  – are likely to be largely unaffected by surface modulation.

Summary and conclusions

This study placed a spatially comprehensive observational constraint on the basal water distribution beneath the GrIS and, hence, regions of the bed at or above the PMP of ice. The distribution of basal water is influenced by, and has influence on, multiple ice sheet and subglacial properties and processes. Subsequently, with a focus on ice-sheet-scale behaviour, we performed an exploratory comparison with related data sets for the GrIS. This included an up-to-date synthesis for the basal thermal state , three different GHF model distributions , bed topography and predicted subglacial flow paths and ice surface speed .

Central to the methods in the study was the use of bed-echo reflectivity variability (rather than bed-echo reflectivity) as a RES diagnostic for basal water. Our use of this diagnostic (a form of edge detector) was motivated by its insensitivity to radar attenuation at the ice-sheet scale, and the pragmatic advantages when performing data combination for multiple RES field campaigns. The reflectivity variability diagnostic is, however, only able to detect wet to dry (or wet to frozen) transitions in bed material and is a sufficient (but not necessary) condition for basal water. It will therefore need to be combined with other information to fully map the extent of basal water and classify basal water bodies.

There was much agreement between the basal water distribution and the thawed marginal regions predicted by . However, we identified regions of basal water and thaw in the interior of the ice sheet that were previously classified as likely to be frozen. The most extensive new region of predicted thaw is a quasilinear corridor feature which extends from NorthGRIP in the south to Petermann in the north. This feature, and the majority of basal water in the northern interior, spatially correlate with elevated GHF inferred from magnetic data by .

The comparison with bed topography and predicted flow paths demonstrated good overall agreement between the basal water storage locations and the geometric constraints imposed by the hydrological pressure potential. However, many of the basal water predictions in the ice-sheet interior occur where ice surface speed (and hence basal motion) is negligible. One plausible explanation is that much of the interior lies on a hard and undeformable bed. Future investigation of basal control on GrIS dynamics should integrate information about basal water and the basal thermal state with better constraints on bed lithology and geology.

The basal water radar data set is available at https://doi.pangaea.de/10.1594/PANGAEA.893097. Additional data generated in this study are available on request from Thomas Jordan. The CSARP level 1B RES data are available from CReSIS at https://data.cresis.ku.edu/data/rds/ (last access: March 2018) and are documented in . The profile in Fig. a is from data segment 2012050804 from the 2012 P3 season and the profile in Fig. b is from data segment 2014051601 from the 2014 P3 season. The Greenland basal thermal state synthesis , ice thickness and topography data sets (BedMachine V3) , and ice surface speed , are archived by NSIDC at 10.5067/R4MWDWWUWQF9, https://nsidc.org/data/idbmg4 (last access: March 2018) and https://nsidc.org/data/NSIDC-0670/versions/1 (last access: March 2018) respectively. The GHF maps by Martos et al. (2018), Fox Maule et al. (2009) and Shapiro and Ritzwoller (2004) are available at https://doi.pangaea.de/10.1594/PANGAEA.892973 (last access: March 2018), http://www.dmi.dk/dmi/dkc09-09.pdf (last access: March 2018) and https://doi.pangaea.de/10.1594/PANGAEA.734145 (last access: March 2018).

Crossover analysis for bed-echo reflectivity variability

To assess the internal consistency of the bed-echo reflectivity variability data, Fig. c, we carried out a crossover analysis at flight-track intersections (defined as bin centres separated by $<\mathrm{2.5}$ km in correspondence with 5 km window size). In this analysis we decomposed the data by the radar system categories in Fig. a using MCoRDS v2 (the most spatially extensive and recent measurements) as a benchmark. We also discounted adjacent postings from the same flight track.

Crossover statistics for ${\mathit{\sigma }}_{\left[R\right]}$ for the different radar instrument classes in Fig. a.

Data comparison	Mean (dB)	SD (dB)
MCoRDs v2 – MCoRDs v2	N/A	1.13
MCoRDs – MCoRDs v2	$-\mathrm{0.43}$	1.18
MCRDS – MCoRDs v2	$-\mathrm{0.07}$	1.69
ACORDS – MCoRDs v2	0.00	1.24

The crossover statistics are shown in Table  and demonstrate standard deviations that range from 1.13 dB (MCoRDS v2 – MCoRDS v2) to 1.69 dB (MCoRDS v2 – MCRDS). The relatively high standard deviation for MCoRDS v2 – MCRDS is likely to arise because the MCRDS coverage is almost exclusively in outlet glacier regions where power variability due to surface roughness and crevassing will be higher. There is no significant bias due to radar centre frequency (MCoRDS v2 and MCoRDS are for 195 MHz, whereas MCRDS and ACORDS are for 150 MHz). However, MCoRDS does have a small negative bias. Since the underlying mechanism for this is unclear we do not empirically correct the data. We note that this is a conservative approach (since ${\mathit{\sigma }}_{\left[R\right]}$ is underestimated) and is logically consistent with ${\mathit{\sigma }}_{\left[R\right]}$ being used as a sufficient condition for basal water (i.e. basal water may be present that is not above the prescribed variability threshold).

The crossover standard deviation values for ${\mathit{\sigma }}_{\left[R\right]}$ in Table  should not be interpreted as standard errors (and are likely significant overestimates) as the flight-track windows do not necessarily sample the same region of the glacier bed. This contrasts with performing crossover analysis of bed-echo power/reflectivity where the Fresnel zone (length scale $\sim \mathrm{100}$ m) defines a spatial overlap. As a form of edge detector, the purpose of the ${\mathit{\sigma }}_{\left[R\right]}$ metric is to identify a signal attributable to basal water within a 5 km region (rather than the coarse grain of an average value). Additionally, the along-track data in Fig.  shows that ${\mathit{\sigma }}_{\left[R\right]}$ can rapidly fluctuate at a 5 km length scale. We therefore point toward the high degree of spatial structure for ${\mathit{\sigma }}_{\left[R\right]}$ in Fig. c and the water predictions in Figs. – as evidence for the robustness of the approach.

The study was initiated by JLB and co-advised by MJS as part of the Basal Properties of Greenland Project. TMJ wrote the paper, analysed the radar data, and, with advice from DMS, developed the water detection method. CNW and MAC developed the GIS mapping environment and modelled the subglacial flow paths. JDP and DMS both advised on radar data processing. YMM contributed the new heat flux distribution and advised on the interpretation of the heat flux models. PH contributed the GISM temperature field and advised on ice-sheet thermodynamics. All authors commented on the draft of the paper and contributed to the interpretation of the science.

JLB is an advisory editor of The Cryosphere. The authors declare that they have no conflict of interests.

Acknowledgements

The project was primarily supported by UK NERC grant NE/M000869/1 as part of the Basal Properties of Greenland project. Thomas M. Jordan was also supported by an EU Horizons 2020 grant 747336-BRISRES-H2020-MSCA-IF-2016. Michael A. Cooper was supported by the UK NERC grant NE/L002434/1 as part of the NERC Great Western Four + (GW4+) Doctoral Training Partnership. Dustin M. Schroeder was supported by a grant from the NASA Cryospheric Sciences Program. We acknowledge the use of data products from CReSIS generated with support from NASA grant NNX16AH54G. We would like to thank Reinhard Drews and an anonymous reviewer for their thoughtful and constructive comments, and Kenichi Matsuoka for being our review editor. We would also like to thank Winnie Chu, Stanford University, for her helpful discussions. Edited by: Kenny Matsuoka Reviewed by: Reinhard Drews and one anonymous referee