Ground-penetrating radar imaging reveals glacier’s drainage network in 3D

. Hydrological systems of glaciers have a direct impact on the glacier dynamics. Since the 1950’s, geophysical studies have provided insights into these hydrological systems. Unfortunately, such studies were predominantly conducted using 2D acquisitions along a few proﬁles, thus failing to provide spatially unaliased 3D images of englacial and subglacial water pathways. The latter has likely resulted in ﬂawed constraints for the hydrological modelling of glacier drainage networks. Here, we present 5 3D ground-penetrating radar (GPR) results that provide high-resolution 3D images of an alpine glacier’s drainage network. Our results conﬁrms a long-standing englacial hydrology theory stating that englacial conduits ﬂow around glacial overdeepenings rather than directly over the overdeepening. Furthermore, these results also show exciting new opportunities for high-resolution 3D GPR studies of glaciers.

geometry and its temporal variations, thus hampering a deeper understanding of these seasonal variations (Hart et al., 2015;Church et al., 2020).
The glacier's hydrological system can be probed using a variety of methods. Direct observations have been made from boreholes (Fountain et al., 2005), tracer testing (Nienow et al., 1996), speleology (Gulley, 2009;Gulley et al., 2009;Temminghoff et al., 2019) and geophysical measurements. The latter include predominantly active (Peters et al., 2008;Zechmann et al., 2018;Church et al., 2019) and passive (Podolskiy and Walter, 2016;Lindner et al., 2019;Nanni et al., 2020) seismic 30 measurements or ground-penetrating radar (GPR) measurements (Moorman and Michel, 2000;Stuart et al., 2003;Irvine-Fynn et al., 2006;Harper et al., 2010;Baelum and Benn, 2011;Hansen et al., 2020). Most GPR studies, published so far, relied on two-dimensional (2D) data, where GPR measurements were acquired along profiles. 2D data sets are typically unable to image complex subsurface structures, such as glacier drainage networks, and the resulting interpretations may thus be inconclusive. However, spatially unaliased 3D GPR data sets (i.e., datasets with a data point spacing smaller than the dominant wavelength of the GPR signals (Sheriff and Geldart, 1995;Grasmueck et al., 2005)) are uncommon in glaciological applications. This is unfortunate, because 3D GPR provides subsurface images that can be viewed from arbitrary directions, thus allowing for an unequivocal interpretation. Such an approach would be particularly useful for glacier drainage networks, and should be 40 feasible because of the strong reflections caused by the very pronounced electrical impedance contrasts at ice/water interfaces (Reynolds, 1997).
In this study we demonstrate the feasibility of and opportunities offered by 3D GPR imaging on glaciers. Therefore, we performed a 3D GPR survey over a temperate alpine glacier to delineate the drainage network and provide much needed hydrological observations to confirm long-standing glacier hydraulic theory. Furthermore, such GPR surveys are able to provide and C-C' are shown in Figure 2 and the crossing points of the GPR profiles are represented by a yellow and red dot. Co-ordinates for all plots are local swiss co-ordinates LV03. Orthophoto was provided by Swiss Federal Office of Topology: Reproduced by permission of swisstopo (JA100120), ©2020 swisstopo (JD100042) 3 Methods

55
To detect and characterise the drainage network located at the glacier's tongue, we acquired a high-resolution 25 MHz 3D GPR survey. This survey was motivated by the results of earlier 2D surveys including active seismic data (Church et al., 2019), repeated GPR measurements on a coarse grid of 2D lines (Church et al., 2020) and a drilling campaign in 2018, which provided initial imaging of a potential drainage network. Between 15 th July 2020 and 23 rd July 2020, GPR data were acquired over the area expected to harbour the englacial conduit network (Fig. 1b).

60
The survey covered an area of 140,000 m 2 , within which the ice thickness varied between 25 and 110 m and where the glacier bed forms a distinct overdeepening (Fig. 1b). The common-offset GPR data were collected using a Sensor & Software PulseEkko ® system with an antenna separation of 4 m and carried by hand at approximately 1 m above the glacier ice surface.
The sampling rate of the GPR system was 1 GHz, giving a time resolution of 1 ns, and thereby providing a vertical spatial resolution in temperate ice of 0.168 m. The use of a large sampling frequency allows small topographical changes from trace to trace of <0.2 m to be observed (King, 2020). A GPR stacking of 4 improved the signal-to-noise ratio and allowed the GPR data to be acquired with average walking speed of approximately 0.4 m per GPR trace. For all GPR lines, a high precision global navigation satellite system (GNSS) continuously recorded the x, y, z coordinates of the centre point between the transmitting and receiving antennas every second. The average accuracy of the GNSS during GPR acquisition was 0.008 m.
The survey area was covered with 281 profiles, resulting in a total profile length of approximately 85 km. The 3D GPR data 70 were collected perpendicular to the ice flow direction and profiles were interspaced 2 m. The line spacing was chosen, such that the diffractions and reflections within the ice body would not become aliased for our antenna frequency of 25 MHz. To ensure data were consistent across the duration of the survey, two GPR lines were always repeated from the previous day's acquisition. Furthermore, six orthogonal profiles were collected along the glacier flow.

75
All GPR common offset data were processed using a combination of an in-house MATLAB-based toolbox (GPRglaz) (Grab et al., 2018) and SeisSpace ProMAX 3-D. The processing was based upon a typical 3D seismic data processing workflow.
Initially, the GPR data were assigned to their corresponding GNSS data. Since the GNSS data were recorded every second and GPR data were recorded approximately every 0.3 seconds, the GNSS data were linearly interpolated to provide the same temporal resolution as the GPR data. The data were then corrected for time zero, and a Butterworth bandpass filter was applied 80 in order to suppress any noise outside the GPR frequency band and thus to increase the signal-to-noise ratio. Overlapping GPR data between different acquisition days were used to investigate whether a data matching filter was required. However, no amplitude matching was required, due to the fact that the GPR equipment used produced stable and repeatable GPR data.
Subsequently, the data were 3D binned using a master grid of 2 m between GPR profiles (inline spacing) and 0.5 m between GPR data points within the profile (crossline spacing). 3D binning consists of assigning each GPR trace to the closest bin 85 centre. As a result of the variable walking speed and walking around crevasses during acquisition, a number of bins had more than a single GPR trace assigned (known as over-fold), whereas other bins (i.e. bin located within a crevasse) were empty. GPR data points were removed from bins that had more than one GPR data point and therefore, this resulted in the bins containing either a single GPR data point or no GPR data point. Such a processing step is required in order not to leave an amplitude imprint on the data during the interpolation and regularisation stages. The data were interpolated and regularised, such that 90 each bin within the survey boundary had an interpolated GPR trace at the centre of the bin.
A 3D Kirchhoff migration algorithm re-positioned the reflected and diffracted signals back to their correct subsurface location. The Kirchhoff migration algorithm was performed using EM wave propagation velocity of ice (0.167 m ns -1 ) and furthermore, the algorithm corrected for amplitude losses from geometrical spreading. Prior to interpretation the data were corrected for attenuation using a Q compensation (Irving and Knight, 2003), a second Butterworth bandpass filter was applied to further improve the signal-to-noise ratio, and the data were stretched from the time to depth domain using a constant velocity of 0.167 m ns -1 .
The 3D interpretation was performed in dGB Earth Sciences OpendTect. The ice-bed interface was manually picked, linearly interpolated, smoothed and constrained to within the survey area. Similarly, the drainage network was picked at a similar horizontal spacing, linearly interpolated and smoothed. The spatial extent of the drainage network was determined from the 100 GPR depth slices (Fig. 3) and also by extracting the root-mean squared amplitude between the surface and the ice-bed interface.
Finally, the reflected GPR amplitudes were extracted from both the ice-bed interface and the drainage network by using a 2 m window (±1 m) centred around the feature.

105
Displaying 3D models adequately is generally a non-trivial task. Below, we discuss the GPR results using a variety of vertical and horizontal cross sections. In our view, such data sets are represented best in form of movies showing scans along different directions. We therefore highly encourage the readers to check the digital supplement.
Selected 2D profiles of the 3D GPR data cube are shown in Figure 2. A water-filled englacial conduit can be identified as a continuous specular strong reflector (Fig. 2a, b, c). The ice-bed interface is identifiable as a weak reflection (Fig 2), indicating 110 that subglacial water is not present. However, in isolated areas, the ice-bed interface has been identified as strong ice-bed reflections (Fig. 2c) and thereby indicating the presence of subglacial water.
The lateral extent of the englacial and subglacial network can be characterised by analysing horizontal slices of the 3D GPR data cube. A slice at 2216 m.a.s.l shows a strong meandering reflection in the northern part of the survey (Fig 3a). The strong reflection is traceable on the eastern edge of the survey before fading out towards the southern edge. At 2213 m.a.s.l (Fig 3b), 115 there is a continuation of the strong reflection, but it becomes more diffusive in the central part of the survey area. At 2210 m.a.s.l (Fig. 3c), we observe another strong meandering reflection that heads southwards towards the terminus of the glacier.
There is an approximately 6 m topographic difference between the drainage network in the survey's northern edge and the southern edge indicating that the imaged drainage network has a shallow global inclination along the flow (<1 degree).

Spatial distribution of drainage network 120
The 3D GPR survey imaged an active meltwater drainage network within the survey boundary. It comprises both, an englacial and subglacial drainage network. The entire drainage network was identified from the GPR data, and reflected amplitudes from the drainage network were extracted (Fig 4a). The conduit network can be delineated by areas of high amplitude (red in Fig.   4a). Furthermore, the drainage network can be split into four separate components labelled in Fig. 4b: (A) A meandering well-defined englacial conduit spanning the overdeepen oriented perpendicular to the glacier flow,    (C) The englacial conduit in B flows into a subglacial drainage system, the subglacial drainage network consists of a single main conduit (Fig 4a), that has a sinusoidal nature, (D) The subglacial conduit in C encounters a basin and flows into another englacial conduit towards the terminus of the 130 glacier.
Given the glacier ice flow direction (N-S) and the ice-thickness distribution, the conduit is expected to flow from north to south. The high-resolution results allow the width of the drainage network to be examined. The mean width of the northern sinusoidal englacial conduit, which flows across the overdeepening (Fig. 4b Section A), is 8 m. As the network flows southwards around the overdeepening (Fig. 4b Section B) the width increases to 11 m. Furthermore, the mean width of the subglacial 135 drainage conduit (Fig. 4b Section C) is 17 m, and finally at the southern edge of the survey site, the mean width of the englacial drainage network is 25 m (Fig. 4b Section D). The thickness of the conduit in section A has previously been investigated in Church et al. (2020), and it was found to be at the limit of the 25 MHz GPR vertical resolution at 0.4 m, when assuming the conduit is water-filled. The conduit thickness in sections B, C and D are also at the limit of the vertical resolution as only a single reflection is visible (Fig. 2c). If the conduit thickness was beyond the vertical resolution, two separate englacial reflec-140 tions would be visible representing the top and bottom of the conduit. Consequently, the channels throughout the study area are thinner than 0.4 m, and therefore, their shape is significantly smaller in the vertical direction than lateral.

Basal reflected amplitude
The amplitude of the ice-bed interface was extracted from the basal horizon as well. This provided insights into the local basal conditions. In the southern region of the survey site, the ice-bed interface reflection has identical amplitude and spatial 145 distribution to the drainage network ( Fig. 5a and c white arrows), thereby indicating this area is positioned along the main drainage network identified in Figure 4. In the northern region of the survey site, instead, there are numerous isolated high amplitudes basal reflections (Fig. 5a), mostly situated within localised flat areas (Fig. 5 b red arrows) indicating a pooling of water. In addition, ubiquitous areas of high amplitude basal reflections are present along the southern region of the survey, indicating the presence of subglacial water (Fig. 5a and c pink arrows). In comparison to the isolated patches in the northern 150 region, these high amplitude basal patches in the southern region appear to be connected to each other indicating the possibility of an additional connected subglacial drainage system.

Comparison of 2D and 3D GPR processing
2D GPR surveys along single profiles are the current approach in today's glaciological applications, although such data sets have imaging limitations. Besides being unable to provide 3D subsurface images, 2D GPR techniques assume all reflections 155 originate from the vertical plane of the acquisition. Complex englacial structure or basal geometry can result in a reflection originating from outside of the acquisition plane, in turn resulting in distortions to the final processed GPR image. Figure 6 shows an example of such a distortion. These distortions are particularly severe for complex geometry of alpine glaciers 3D GPR acquisition and processing are able to mitigate these distortions. The Rhonegletscher GPR data cube is the product of a 3D processing workflow and with the employment of 3D migration over conventional 2D migration, the distortions from 160 out-of-plane reflections are removed and an improvement in lateral resolution is gained. A 3D migration effectively collapses the Fresnel zone in both inline and crossline directions, thereby reducing the lateral resolution to the 3D bin size (2 m x 0.5 m).
This lateral resolution improved leads to improvements in subsurface imaging, as two closely laterally-separated reflectors are able to be imaged as two individual reflectors.
A comparison between GPR data processed with two different workflows (2D GPR workflow ( Figure 7a) and 3D GPR 165 workflow (Fig. 7b)) highlights the imaging differences on both the englacial conduit network and the ice-bed interface. Generally, both workflows produce similar subsurface images however, there are subtle differences that indicate a more unambiguous interpretation with the 3D GPR workflow. The ice-bed interface in the 3D GPR data cube has increased reflector continuity in comparison to the 2D workflow ( Fig. 7 brown arrows). Furthermore, the englacial conduit imaged using a 3D GPR workflow has fewer artefacts and is absent of events that are incorrectly intersecting the englacial conduits ( Fig. 7 blue arrows).

Geometry of drainage network
This is the first time that a glacier's drainage network is imaged in 3D with GPR data, thus providing an unprecedented, highresolution information of the geometry from such a system. In our case, it has a meandering nature throughout the survey area, and the drainage network's width increases towards the terminus of the glacier. Moreover, it consists of a single dominant 175 conduit that alternately flows englacial and subglacial. Such a drainage network is known as an efficient channelised network.
Studies from both cold-ice (Chandler et al., 2013) and temperate-ice (Nienow et al., 1996(Nienow et al., , 1998Mair et al., 2002)  have shown that early in the melt season, the glacier's drainage network is slow and inefficient. Typically, it evolves into a fast channelised drainage network just before the peak of the glacier's discharge. Since the peak discharge for Rhonegletscher is typically mid-August (GLAMOS, 2017) and thus occurred one month after the acquisition, the drainage network is expected 180 to be in a channelised configuration. The 3D GPR imaging results confirm that the network consists of a single, wide conduit, that alternately flows through englacial and subglacial channels, known as Röthlisberger channels (Röthlisberger, 1972).
The theoretical shape of englacial drainage conduits is circular, however the drainage network imaged on Rhonegletscher is up to a maximum of 60 times wider than its thickness. Such an observation contradicts the theory of circular conduit crosssectional shape proposed independently by Shreve and Röthlisberger (Röthlisberger, 1972;Shreve, 1972) but is in line with 185 the further development by Hooke et al. (1990). Lliboutry (1983) hypothesised that when water encounters an overdeepening, the water flows along the glacier's flank as so-called gradient conduits, therefore avoiding the deepest part of the overdeepening (Cook and Swift, 2012). According to Lliboutry's theory, these conduits should be located at the same altitude as the lowest point of the riegel that produces the overdeepening. The 3D GPR data suggest that in the case of Rhonegletscher, the flow paths indeed route meltwater around the overdeepening rather than across it ( Fig. 4b Labelled: B). Similarly, the elevation of the englacial conduit that is located around the overdeepening coincides with the elevation of the riegel and also the proglacial lake level. These observations support the long-standing theory, but was never verified by field observations. Subglacial and englacial water flows as a response of changing hydraulic potentials. This hydraulic potential can be estimated by assuming spatially uniform flotation fraction as described in Flowers and Clarke (1999) and within the imaged 195 Rhonegletscher's drainage network the water flow routing followed the gradient of the hydraulic potential and not along a single hydraulic potential contour.

Water accumulation in temperate glaciers
In addition to the detection of the primary drainage network, the 3D GPR data provided evidence of subglacial water accumulation. GPR-based detection of large amounts of subglacial water, such as subglacial lakes is well established (Ridley et al., 200 1993;Siegert et al., 2005;Palmer et al., 2013), but their spatial extent is often unclear as a result of the limitations of 2D GPR surveys. This is in contrast to our 3D GPR data set, which can be used to delineate such extents. In addition, 3D GPR has the potential to identify smaller subglacial water accumulations, such as expected to occur within water-filled cavities. Subglacial cavities can form, when the sliding ice uncouples from the glacier bed as a result of either rapid glacier sliding or pronounced bed roughness (Nye, 1970). Two types of subglacial cavity system are generally distinguished -isolated cavities and linked cavities -and both cavity systems alter the glacier's dynamics (Lliboutry, 1976(Lliboutry, , 1979Hoffman et al., 2016;Rada and Schoof, 2018). The high amplitude reflections along the basal interface (Fig 5a) represent water accumulations and appear to be isolated from each other, thereby indicating that they are potentially isolated water-filled cavities forming an inefficient drainage network. It is interesting to note both an inefficient drainage network and an efficient network can coexist in overdeepenings (Hooke et al., 1990;Rada and Schoof, 2018). Although our data provide an instantaneous image of these systems, repeated 3D 210 GPR surveys could also yield insights into their temporal dynamics.

Future of 3D GPR within glaciology
A 3D approach as presented within this contribution is feasible and highly beneficial over 2D GPR for detecting and quantifying dimension of a glacier's hydrological network. For large-scale investigations in Greenland and Antarctica, it will be more challenging to conduct 3D GPR surveys as a result of their spatial distribution and therefore, 2D GPR acquisition will likely 215 continue to prevail. However, future radar surveys could be complimented with the use of 3D GPR acquisition, in order to reduce the ambiguity of interpretations in places of interest.
The major limiting factors of such future 3D GPR surveys are the rate of acquisition and the accessibility of the field site.
Both of these concerns can be addressed with drone technology. Drone technology is often used in cryosphere research (Gaffey and Bhardwaj, 2020) however GPR-based drone surveys are currently limited to landmine detection (Colorado et al., 2017;