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

. There is widespread, but often indirect, evidence that a signiﬁcant 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 reﬂectivity, but unambiguous mapping is limited by uncertainty in

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 wetdry 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 Peter-mann 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.
T. M. Jordan et al.: Greenland basal water from radar bed echoes 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. van de Wal et al., 2008), basal melting due to geothermal heat (e.g. Dahl-Jensen et al., 2003;Rogozhina et al., 2016), frictional and shear heating (e.g. van der Veen, 2013), and transport processes (surface, englacial and subglacial), which redistribute water (e.g. Rennermalm et al., 2013;Chu, 2014). 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 (Palmer et al., 2013(Palmer et al., , 2015Howat et al., 2015;Willis et al., 2015). 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 (Chu et al., 2016) as wet sediment (Christianson et al., 2014) and as part of channelised networks (Livingstone et al., 2017).
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 (Weertman, 1968;Gundestrup and Hansen, 1984;Dahl-Jensen et al., 1998;Cuffey et al., 1995;MacGregor et al., 2016), with the thawed bed at North-GRIP being an exception (Andersen et al., 2004). Toward the ice-sheet margins boreholes generally indicate thawed beds (MacGregor et al., 2016). Indirect methodologies for discriminating frozen and thawed beds (ice-sheet model predictions, radiostratigraphy and surface properties) were recently combined by MacGregor et al. (2016) to produce a frozenthawed 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 (Dahl-Jensen et al., 2003;Greve, 2005;Rogozhina et al., 2016), hydrology (Rogozhina et al., 2016) and flow features (Fahnestock et al., 2001) 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 (Fahnestock et al., 2001), which, in turn, has recently been attributed to the path of the Iceland hotspot track (Rogozhina et al., 2016;Martos et al., 2018). As with basal temperature, the sparsity of borehole measurements limits direct inference of GHF (which is related to the vertical gradient of basal tem-perature). Instead, a range of geophysical techniques including tectonic (Pollack et al., 1993), seismic (Shapiro and Ritzwoller, 2004) and magnetic (Fox Maule et al., 2009;Martos et al., 2018) 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 (Rogozhina et al., 2012). Due to the relationship between basal melt and GHF, basal water observations can be used to further refine and cross-validate GHF models (Siegert and Dowdeswell, 1996;Schroeder et al., 2014;Rogozhina et al., 2016).
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. Peters et al., 2005;Jacobel et al., 2009;Oswald and Gogineni, 2008) 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 (Bogorodsky et al., 1983;Martinez et al., 2001;Peters et al., 2005). However, due to uncertainty and spatial variation in radar attenuation (an exponential function of temperature Corr et al., 1993) bed-echo reflectivity is subject to spatial bias and can be ambiguous when mapped over larger regions (Matsuoka, 2011;Mac-Gregor et al., 2012;Jordan et al., 2016). 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) Young et al., 2016) and the bed-echo abruptness (a waveform parameter) Gogineni, 2008, 2012) -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-sheetscale 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 upto-date analyses for the basal thermal state (MacGregor et al., 2016) and GHF (Shapiro and Ritzwoller, 2004;Fox Maule et al., 2009;Martos et al., 2018) 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 Martos et al. (2018). We then compare basal water and bed topography (Morlighem et al., 2017), 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 (Joughin et al., 2010(Joughin et al., , 2016.

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) (Paden, 2015). The flighttrack coverage, subdivided by the radar system, is shown in Fig. 1a with the seasonal breakdown: ACORDS 2003 P3 and2005 TO;MCRDS 2006P3, 2008MCoRDs 2010 P3 and2010 DC8;MCoRDs v2 2011TO, 2011P3, 2012P3, 2013P3, 2014. More details on the track lengths and data segmentation can be found in Mac-Gregor et al. (2015a). 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 Gogineni et al., 2014;Paden, 2015). 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 ∼ 4.3 to 20 m in ice. The along-track (horizontal) resolution also varies between field seasons and is typically ∼ 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-byseason 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 Oswald and Gogineni (2008), Jordan et al. (2016), Jordan et al. (2017) 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 ∼ 17 dB threshold used in Jordan et al. (2016Jordan et al. ( , 2017 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 Jordan et al. (2016Jordan et al. ( , 2017. 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. 1b, which shows good-(SNR > 10 dB) and poor-(SNR ≤ 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 (Oswald and Gogineni, 2008). Additionally, since roughness scattering loss is frequency dependent (MacGregor et al., 2013), 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 ∼ 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 (Bogorodsky et al., 1983;Peters et al., 2005). The reflectivity, [R], is obtained from solving the decibel form of the bed-echo power equation 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 (Zwally et al., 2012) and major regions of fast-flow identified from ice surface speed (Joughin et al., 2010(Joughin et al., , 2016. Abbreviations in (c) correspond to Camp Century (CC), Humboldt (Hu), Petermann (Pe), Ryder (Ry), Northeast Greenland Ice Stream (NEGIS), north-western margins (NWM), Jakobshavn Isbrae (JI), Kangerlussuaq (Ka), Helheim (He) and Ikertivaq (Ik). The projection is a polar stereographic north (70 • N, 45 • W) and is used in all future plots.
where [P ] is the bed-echo power (in this case the aggregated bed-echo power, which mitigates for scattering losses from the ice-bed interface), [S] is the system performance, [G] is the geometric correction, [L] is the attenuation loss in ice, and [B] is the loss due to birefringence (ice fabric anisotropy), and the notation [X] = 10log 10 X is assumed (Matsuoka et al., 2010). Rough surface losses from transmission through the air-ice interface (e.g. Grima et al., 2014;Schroeder et al., 2016a) also influence the radar power budget and are discussed in more detail in Sect. 4.6. The geometric correction for a specular reflector can be defined by where s is the aircraft height and h is the ice thickness and ǫ ice = 3.15 is the relative dielectric permittivity of ice to give the geometrically corrected bed-echo power Finally, setting [L] = 2 < N > h gives where < N > (dB km −1 ) is the (one-way) depth-averaged attenuation rate (Matsuoka et al., 2010). A fundamental ambiguity in bed-echo reflectivity analysis is that there are two unknowns in Eq. (5): < N > and [R]. Two approaches are typically used to determine < N >: (i) forward modelling using estimates of attenuation as a function of englacial temperature (e.g. MacGregor et al., 2007;Matsuoka et al., 2012a;Chu et al., 2016) and (ii) empirical determination using the linear regression of bed-echo power and ice thickness (e.g. Jacobel et al., 2009;Schroeder et al., 2016b). Attenuation follows an Arrhenius (exponential) relationship with temperature and a linear dependence on the concentration of ionic impurities: primarily hydrogen (H + ) but also chlorine (Cl − ), and ammonium (NH + 4 ) (Corr et al., 1993;Wolff et al., 1997;MacGregor et al., 2007MacGregor et al., , 2015b. On an ice-sheet scale, the uncertainty when forward modelling < N > is so high that it can be prohibitively challenging to accurately calibrate [R] (Matsuoka, 2011;Mac-Gregor et al., 2015b;Jordan et al., 2016). 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 (MacGregor et al., 2007(MacGregor et al., , 2015b. 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 (Matsuoka et al., 2012a), or when there is a thickness-correlated distribution in bed-echo reflectivity (Jordan et al., 2016).
We will later demonstrate that, unlike absolute values of [R], local variability in bed-echo reflectivity is highly insensitive to modelled values of < N > (Sect. 2.6). 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 Jordan et al. (2016) and uses the M07 Arrhenius equation MacGregor et al. (2007), the Greenland Ice Sheet Model (GISM) temperature field from Huybrechts (1996) as updated in Goelzer et al. (2013), depth-averaged ionic concentrations from the GRIP ice core (MacGregor et al., 2015b), and the Greenland ice thickness data set in Bamber et al. (2013a). The temperature field derives from a full 3-D thermomechanical simulation over several glacialinterglacial 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 Bamber et al. (2013a) ice thickness. The geothermal heat flux in GISM was initially taken from Shapiro and Ritzwoller (2004) but further adjusted with Gaussian functions around the deep ice core sites to match observed basal temperatures. Vertical temperature profiles are within 1-2 • C when compared to available in situ measurements.
Using this model framework, it is that predicted that < N > varies by a factor ∼ 5 over the GrIS, ranging from ∼ 6 dB km −1 in the colder northern interior to ∼ 30 dB km −1 toward the warmer south-western margins (refer to Fig. 5a in Jordan et al. (2016) 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 σ [P g ] and σ [R] . The numerical calculation is analogous to how topographic roughness (the rms height) is calculated from bed elevation profiles (Shepard et al., 2001), 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 MacGregor et al. (2016), 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. 2.5) -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 nonstandard approach to bed-echo data analysis, we now take a closer look at the statistical properties. The formula for σ [P g ] follows from the variance of Eq. (4) and is given by where σ [L] is the standard deviation in attenuation loss, and σ [R], [L] is the covariance of bed-echo reflectivity and attenuation loss. Using [L] = 2 < N > h and assuming < N > can be approximated as constant (justifiable at the 5 km length scale that is considered); then Eq. (6) becomes where σ h is the standard deviation of ice thickness and σ [R],h is the covariance of bed-echo reflectivity and ice thickness. Equations (6) and (7) consider the variability of logtransformed 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 (Bogorodsky et al., 1983;Peters et al., 2005). This includes prior applications that have considered the distribution and standard deviation of dB reflection amplitudes in the context of water detection (e.g. Wolovick et al., 2013;MacGregor et al., 2013). In regions where σ [R],[L] ≈ 0 (bed-echo reflectivity has negligible covariance with attenuation loss), Eqs. (6) and (7) are approximated by and the loss component of σ [P g ] is solely modulated by σ [L] which is proportional to the product < N > σ h . Whilst an approximation (in certain circumstances the second and third terms on the right-hand side of Eqs. (6) and (7) can be of comparable magnitude), this scenario provides an intuitive way to understand the interrelationship between σ [P g ] , σ [R] and σ [L] . Two example profiles for and h are shown in Fig. 2. Figure 2a is an example from the interior of the ice sheet where σ [L] is relatively low and h is thick (∼ 2.8 km). Subsequently the profiles for σ [P g ] and σ [R] are very similar in appearance, with the most notable difference at a distance of ∼ 360 km where there is higher σ [P g ] due to the subglacial trough. The peaks in σ [R] are later related to wet-dry bed material transitions in Sect. 2.5. It is also important to note that σ [R] can be greater than σ [P g ] (e.g. at a distance of ∼ 342 km), which is an effect that can be explained by the covariance between attenuation loss/ice thickness and bed-echo reflectivity in Eqs. (6) and (7). Figure 2b is a representative example from toward the ice-sheet margins, where σ [L] is higher due to more rapid variation in h (more complex bed topography) and higher values of < N > (warmer ice). In this case, σ [P g ] is noticeably greater than σ [R] , with the differences largely attributable to higher σ [L] as anticipated by Eq. (8). The examples in Fig. 2 highlight that, at this length scale, higher values of σ [R] arise primarily due to a large single transition in [R]. However, there are also instances where multiple transitions result in local variability peaks (e.g. at distance ∼ 380-385 km along Fig. 2a).
When calculating σ [P g ] , σ [R] and σ [L] , bed-echo coverage gaps within a 5 km bin (see Fig. 1b) 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. 2b, where, aligned with the deep subglacial trough at distance ∼ 1324 km, there are along-track gaps in σ [P g ] , σ [R] and σ [L] .
It is important to clarify the difference between the use of bed-echo power and reflectivity variability in this study from previous radioglaciology studies (Neal, 1982;Peters et al., 2005;Carter et al., 2007;Grima et al., 2014). 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. 2.1). We are therefore able to focus on power variability that is a result of along-track changes in the bed dielectric.  Fig. 2b). 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 (×8 the resolution of the ice-sheet-scale plots).

Distributions for bed-echo power and reflectivity variability
Spatial distributions for the variability measures: are shown in Fig. 3a, b, c along with ice thickness (Morlighem et al., 2017) in Fig. 3d. In general, σ [P g ] 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 σ [P g ] that is attributable to σ [L] , results in a more uniform ice-sheet-scale distribution of σ These examples serve to further illustrate the spatial interrelationship between σ [P g ] , σ [R] and σ [L] in a typical interior region with lower σ [L] and a typical marginal region with higher σ [L] . Its is clear that the interior example has very similar spatial distributions for σ [P g ] and σ [R] , whereas the marginal example has higher σ [P g ] associated with the higher σ [L] that   Fig. 3c). Later in the study σ [R] > 6 dB is used as the threshold criteria for diagnosing the presence of basal water. 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 σ [R] .
The corresponding frequency distributions for σ [P g ] and σ [R] are shown in Fig. 4. Both demonstrate a strong positiveskew, with a long-tail extending to higher values. The mean and standard deviation for σ [P g ] is greater than σ [R] . This is consistent with the commonly observed result that making an attenuation correction to [P g ] acts to reduce the overall decibel range for [R] (e.g. Oswald and Gogineni, 2008;Schroeder et al., 2016b), hence more closely resembling the predicted dB range for bed materials (Peters et al., 2005).
A quantification and discussion of crossover statistics for σ [R] is given in Appendix A.

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 where ǫ ∼ ice and ǫ ∼ bed are the complex dielectric permittivity of the glacier ice and bed material respectively. The rela- Table 1. Dielectric and reflective properties of subglacial materials based on a compilation of past values by Bogorodsky et al. (1983), Martinez et al. (2001) and Peters et al. (2005). 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 [R]. tive (real) part of the permittivity, ǫ bed , is the primary control on [R]. 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 Bogorodsky et al. (1983), Martinez et al. (2001) and Peters et al. (2005). 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 (Martinez et al., 2001;Peters et al., 2005). In general, lower values of ǫ bed and [R] occur for dry or frozen bed materials (approximately ǫ bed < 7 and [R] < −14 dB), whilst higher values occur for wet bed materials (approximately ǫ bed > 7 and [R] > −14 dB). Dielectric mixing between bed materials can also occur at the length scale of the Fresnel zone (∼ 100 m), which results in further averaging of the observed reflectivity (Berry, 1975;Peters et al., 2005). 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 σ [R] . A simple two-state dielectric model, does, however, enable us to physically motivate the water diagnostic in terms of dielectric properties (Fig. 5). The model assumes that the along-track sample window is comprised of two different bed materials: the dry background bed material with permittivity ǫ dry and reflectivity [R dry ], and the wet material with permittivity ǫ wet and reflectivity [R wet ]. 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

Bed material
where [R] = [R wet ] − [R dry ] is the reflectivity difference between wet and dry beds. Equation (11) is derived by con-

Reflectivity variability,
[R] (dB) sidering 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 Fig. 5c, and shows that, for fixed [R], σ [R] is maximised when f = 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. Jacobel et al., 2009;Chu et al., 2016). In this study, a similar thresholding approach is applied to the distribution of σ [R] (Fig. 4b). The threshold choice for basal water (σ [R] > 6 dB) corresponds to the region greater than the σ [R] = 6 dB contour in Fig. 5c and requires a minimum wet-dry reflectiv-ity difference of [R] > 12 dB. In general, [R] > 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 = 0.5, [R] = 17 dB) has σ [R] = 8.5 dB. The contours in Fig. 5c demonstrate that small perturbations to even mixing (f = 0.5) produce similar σ [R] , and hence that water detection is insensitive to the discretisation of the along-track sample window (Sect. 2.3). Overall, the threshold choice (σ [R] > 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. 2.3, the 6 dB variability threshold applies to the distribution of log-transformed reflectivity.
The bed-echo power aggregation in Sect. 2.1 partially mitigates for roughness-induced scattering loss and the alongtrack power variability associated with this. Supporting evidence is that the crossover analysis for σ [R] (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 σ [R] value increasing.
Table 1 also indicates that, in exceptional circumstances, [R] > 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. 1b.

Basal water distribution and robustness to attenuation model bias
The initial basal water predictions (σ [R] > 6 dB, pre-sensitivity analysis) are shown in Fig. 6 (red, blue and green data), and correspond to ∼ 3.5 % of bins containing predominantly good-quality bed echoes (Fig. 1b). A full geographic analysis of the spatial distribution is performed in Sect. 3. To demonstrate the robustness of the predictions, we performed a sensitivity analysis with respect to the modelled attenuation correction < N > (Sect. 2. 2) The analysis considered a series of increasingly large (uniform, multiplicative) perturbations to < N > and then tested whether σ [R] > 6 dB also held for the perturbed model. Examples of persistent water predictions for ±20 % (red and green data) and ±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 ∼ 2.6 % and 2.1 % of the along-track bins for ±20 % and ±50 % respectively). The sensitivity analysis tests the robustness of the water predictions to a number of different physical scenarios. uniform perturbation considered in Fig. 6 was proposed by MacGregor et al. (2015b) 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 (Jordan et al., 2016), ±20 % is a reasonable estimate for temperature-related uncertainty. Subsequently, in the comparison with other data sets in Sect. 3 the subset of red and green points in Fig. 6 is used. Inherent bias in the Arrhenius equation parameters could be significantly higher than temperature uncertainty (MacGregor et al., 2015b). However, since the spatial structure for the basal water distribution under the ±50 % perturbation is largely preserved, this is unlikely to significantly alter the conclusions that are drawn. Additionally, whilst σ [R] > 6 dB is used as the threshold, Fig. 3c demonstrates that a less conservative threshold (σ [R] > 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 σ [R] with respect to uncertainty and model bias in < N > (particularly compared with bed-echo reflectivity, [R]). An analogous sensitivity analysis by Jordan et al. (2016) demonstrated that systematic over-and underestimates in < N > lead to pronounced ice-thickness-correlated biases in the distribution for [R] in northern Greenland (Fig. B1 in Jordan et al., 2016).

Results
The basal water distribution is now compared with existing analyses for the basal thermal state ( Fig. 6. 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. 1b) 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
In Fig. 7a the basal water predictions are underlain by the basal thermal state synthesis (frozen/thawed likelihood) map by MacGregor et al. (2016). 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 bedecho data and is therefore independent of our basal water predictions.
The reflectivity variability water diagnostic enables a positive discrimination of basal thaw, since σ [R] > 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 (σ [R] < 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.
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. 1c. 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. 7b-f. The most obvious disagreement is the quasilinear "corridor" of basal water in the north-central ice sheet (Fig. 7d). This feature  Martos et al. (2018) 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).
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. 7b), in the far north (Fig. 7c), to the east of GRIP (Fig. 7e), and around Kangerlussuaq (Fig. 7f). 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 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. van der Veen, 2013). 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 an-ticipated to be ∼ 55-70 mW m −2 in the interior of the ice sheet (Dahl-Jensen et al., 2003;Greve, 2005;Buchardt and Dahl-Jensen, 2007). In the water-GHF comparison we therefore define elevated GHF (i.e. likely to produce basal melt) as > 60 mW m −2 . This definition is also informed by the lower range of values (37-50 mW m −2 ) that are typically associated with non-altered ancient continental crust (Artemieva, 2006;Rogozhina et al., 2016). In Fig. 8 the basal water predictions are underlain by three different GHF models: the seismic model by Shapiro and Ritzwoller (2004) and two models derived from magnetic anomalies by Fox Maule et al. (2009) and . The GHF model by Shapiro and Ritzwoller (2004) is based on the correlation between a 3-D tomographic model of the crust and mantle temperature. The GHF models by Fox Maule et al. (2009) andMartos et al. (2018) are based on a thermal model of the lithosphere with the lower boundary defined by the Curie depth, which is determined from magnetic anomalies. Martos et al. (2017) further describes this approach and the additional spectral processing method used to produce Fig. 8c. An older tectonic GHF model by Pollack et al. (1993) is not considered and a spatial plot for the GrIS can be found in Rogozhina et al. (2012)  profiles and thermomechanical model inversions (Weertman, 1968;Dahl-Jensen et al., 1998, 2003Greve, 2005;Buchardt and Dahl-Jensen, 2007;Petrunin et al., 2013) are provided by Rezvanbehbahani et al. (2017), Martos et al. (2018) and demonstrate general consistency between Fig. 8c and local estimates at GRIP, NEEM, NorthGRIP and Camp Century. Local estimates of GHF at Dye 3 (∼ 20-25 mW m −2 ) are significantly lower than all three GHF models. In interpreting Fig. 8, we limit the comparison to the ice sheet interior where the spatial correlation between GHF and basal water should be strongest. The model by Shapiro and Ritzwoller (2004) , Fig. 8a, predicts low GHF (< 60 mW m −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 Fox Maule et al. (2009) , Fig. 8b, 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 Martos et al. (2018) , Fig. 8c, exhibits strong overall spatial correla-tion 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. 8 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 Martos et al. (2018) is a consequence of the spectral method (window size and overlap), which has an effective resolution of ∼ 75 km.

Comparison between basal water distribution, bed topography and subglacial flow paths
In Fig. 9 the basal water predictions are underlain by the most recent Greenland bed topography digital elevation model cussion 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 (Shreve, 1972), which was calculated using the bed elevation and ice thickness surfaces at a grid cell resolution of 600 m (derived from Morlighem et al., 2017). The flow-routing algorithm was implemented in Ar-cGIS using the inbuilt flow accumulation tool on a hydraulic potential surface that had been processed for sink removal using the method of Wang and Li (2006). Likely hydrological flow paths were identified by excluding flow paths where 50 or fewer neighbouring cells cumulatively contribute to a given location. Figure 9 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. 9b and north-western margins in Fig. 9c. Figure 9b also demonstrates that basal water is present along sections of the "mega-canyon" feature iden-tified by Bamber et al. (2013b) -for example, north-west of the intersection (80 • N, 50 • 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. 9d).

Comparison between basal water distribution and ice surface speed
In Fig. 10 the basal water predictions are underlain by ice surface speed (Joughin et al., 2016), 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 Joughin et al. (2010). 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. 10b).
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. 10c), and the Helheim catchment (Fig. 10d).
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.

Basal water, basal thermal state and temperature
The basal water distribution in this study and the basal thermal state synthesis by MacGregor et al. (2016) 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 MacGregor et al., 2016) 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 MacGregor et al., 2016) all occur where the ice surface speed is low. This includes the north-central ice divide (Fig. 7d), the region east of GRIP (Fig. 7e), and the region west of Kangerlussuaq (Fig. 7f). The regions of agreement and disagreement are, perhaps, unsurprising, since three of the four methods employed by MacGregor et al. (2016) -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 MacGregor et al. (2016) 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. Shapiro and Ritzwoller, 2004;Fox Maule et al., 2009).
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 highreflectivity 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: ∼ 10 km south of GRIP and ∼ 7 km north-east of Camp Century (Fig. 7). At GRIP this is less surprising, since the basal temperature is 6 degrees below the PMP (Dahl-Jensen et al., 1998;MacGregor et al., 2016) and GHF is likely to be elevated in this region (Fig. 8). 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 (Weertman, 1968;MacGregor et al., 2016). 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 (Beem et al., 2017). 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 −13 • C) (Murray et al., 2012) and at the Devon Ice Cap in the Canadian Arctic (Rutishauser et al., 2018).

Basal water and geothermal heat flux
The comparison between basal water and the different GHF models in Fig. 8 (Shapiro and Ritzwoller, 2004;Fox Maule et al., 2009;Martos et al., 2018) demonstrates greatest consistency with the distribution by Martos et al. (2018). 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 Rezvanbehbahani et al. (2017) 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 (Huybrechts, 1996). Other studies include investigating the sensitivity of ice-sheet dynamics to the thermal boundary condition (Seroussi et al., 2013) or basal lubrication (Shannon et al., 2013), and thermal models of the underlying lithosphere (Rogozhina et al., 2016).
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 (Rogozhina et al., 2016;Martos et al., 2018). The magnetic GHF map in Fig. 8c, 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. 8c from Petermann Glacier to NorthGRIP and follows a path previously anticipated by Forsyth et al. (1986). 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 (Cooper et al., 2016) and Humboldt catchments (Livingstone et al., 2017), along with the prominent megacanyon feature which extends from the NorthGRIP region in the south to the Petermann Glacier in the north (Bamber et al., 2013b). The comparison between the predicted flow paths and basal water in Fig. 9 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. 9b 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. Schroeder et al., 2014) and use of a coupled hydrological ice flow model (e.g. Le Brocq et al., 2009).

Basal water and ice-sheet motion
Both observational (e.g. Moon et al., 2014;Tedstone et al., 2013) and theoretical studies (e.g. Creyts and Schoof, 2009;Schoof, 2010) point toward a complex spatio-temporal relationship between basal water and ice surface speed in fastflowing 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 reemphasise 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 (Dawes, 2009;Henriksen, 2008) and in the NEGIS drainage basin (Christianson et al., 2014). 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 Layberry and Bamber (2001) 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 Gogineni (2008, 2012), 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 Gogineni (2008, 2012) will generally not be able to discriminate water in many outlet glaciers and tributaries including Petermann and the northwestern margins . This is because these regions tend to exhibit a diffuse scattering signature (associated with fine-scale roughness), whereas the method proposed by Gogineni (2008, 2012) 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 Rippin (2013) and Jordan et al. (2017) 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 Bons et al., 2016;Dow et al., 2018). The features have been attributed to a supercooling freeze-on process , stick-slip mechanisms,  and the rheological anisotropy of ice (Bons et al., 2016). 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. 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 < 5 km in horizontal extent (Palmer et al., 2013(Palmer et al., , 2015Howat et al., 2015;Willis et al., 2015) 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 Chu et al. (2018). 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. Peters et al., 2005;Jacobel et al., 2009;Schroeder et al., 2013;Young et al., 2016) 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 bedreturned power. Notably, we did not account for (i) power modulation due to birefringent propagation (ice fabric anisotropy, e.g. Matsuoka et al., 2012b), (ii) rough surface transmission (e.g. Grima et al., 2014;Schroeder et al., 2016a) 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 (Martin et al., 2009;Drews et al., 2012). However, the zoomed-in panels in these regions (e.g. Figs. 7b,d,9d) 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. MacGregor et al., 2016) and surface melt (e.g. van de Wal et al., 2008). The degree of surface-induced power variability depends on both the surface permittivity and the roughness regime (relative to the radar wavelength) (Schroeder et al., 2016a). 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 MacGregor et al. (2016). In turn, the central results in this study -the new regions of basal water and thaw identified in Fig. 7 -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 (MacGregor et al., 2016), three different GHF model distributions (Shapiro and Ritzwoller, 2004;Fox Maule et al., 2009;Martos et al., 2018), bed topography (Morlighem et al., 2017) and predicted subglacial flow paths and ice surface speed (Joughin et al., 2010(Joughin et al., , 2016. Central to the methods in the study was the use of bedecho 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 Mac-Gregor et al. (2016). 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 Martos et al. (2018).
The comparison with bed topography (Morlighem et al., 2017) 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.
Data availability. 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 CRe-SIS at https://data.cresis.ku.edu/data/rds/ (last access: March 2018) and are documented in Paden (2015). The profile in Fig. 2a is from data segment 2012050804 from the 2012 P3 season and the profile in Fig. 2b is from data segment 2014051601 from the 2014 P3 season. The Greenland basal thermal state synthesis (MacGregor et al., 2016), ice thickness and topography data sets (BedMachine V3) (Morlighem et al., 2017), and ice surface speed (Joughin et al., 2016), are archived by NSIDC at https://doi.org/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)   Appendix A: Crossover analysis for bed-echo reflectivity variability To assess the internal consistency of the bed-echo reflectivity variability data, Fig. 3c, we carried out a crossover analysis at flight-track intersections (defined as bin centres separated by < 2.5 km in correspondence with 5 km window size). In this analysis we decomposed the data by the radar system categories in Fig. 1a using MCoRDS v2 (the most spatially extensive and recent measurements) as a benchmark. We also discounted adjacent postings from the same flight track.
The crossover statistics are shown in Table A1 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 σ [R] is underestimated) and is logically consistent with σ [R] 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 σ [R] in Table A1 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 bedecho power/reflectivity where the Fresnel zone (length scale ∼ 100 m) defines a spatial overlap. As a form of edge detector, the purpose of the σ [R] 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 alongtrack data in Fig. 2 shows that σ [R] can rapidly fluctuate at a 5 km length scale. We therefore point toward the high degree of spatial structure for σ [R] in Fig. 3c and the water predictions in Figs. 7-10 as evidence for the robustness of the approach.