Detecting supraglacial debris thickness with GPR under suboptimal conditions

Abstract

The thickness of a supraglacial layer is critical to the magnitude and time frame of glacier melt. Field-based, short pulse, ground-penetrating radar (GPR) has successfully measured debris thickness during a glacier's melt season, when there is a strong return from the ice–debris interface, but profiling with GPR in the absence of a highly reflective ice interface has not been explored. We investigated the performance of 960 MHz signals over 2 km of transects on Changri Nup Glacier, Nepal, during the post-monsoon. We also performed laboratory experiments to interpret the field data and investigate electromagnetic wave propagation into dry rocky debris. Laboratory tests confirmed wave penetration into the glacier ice and suggest that the ice–debris interface return was missing in field data because of a weak dielectric contrast between solid ice and porous dry debris. We developed a new method to estimate debris thicknesses by applying a statistical approach to volumetric backscatter, and our backscatter-based calculated thickness retrievals gave reasonable agreement with debris depths measured manually in the field (10–40 cm). We conclude that, when melt season profiling is not an option, a remote system near 1 GHz could allow dry debris thickness to be estimated based on volumetric backscatter.

Introduction

Nearly 19% of glacier area across Eastern South Asia is debris-covered compared to a figure of only 7.3% on mountain glaciers excluding Antarctica and Greenland (Herreid and Pellicciotti, 2020). The share of glacier area covered in debris is greater than average in the eastern Himalaya, at 25% (Kraaijenbrink and others, 2017) in the East Nepal region defined by Kääb and others (2015) (Brun and others, 2017). Debris layers form in a glacier's ablation zone and are, consequently, highly relevant to glacier melt. Understanding the impact of climate change on all glaciers in High Mountain Asia is critical to predicting the future water supply of a highly populated region (Lutz and others, 2014; Brun and others, 2017; Rowan and others, 2017). The exact response of debris-covered glaciers to climate change differs from that of clean glaciers and is important because debris cover is projected to increase in a warming climate (Kraaijenbrink and others, 2017; Scherler and others, 2018). Measurements of debris thickness are needed because debris thickness determines the amount of melt. As shown through experiments (Östrem, 1959; Reznichenko and others, 2010) and in situ field measurements (Mattson and others, 1993; Nicholson and Benn, 2006), a debris-covered glacier's response to climate depends on debris thickness. Where less than a few centimeters, debris enhances underlying glacier melt because it has a lower albedo relative to ice and snow. Over a few centimeters thick, debris insulates.

Debris, predominantly rockfall from steep valley walls (Menzies and van der Meer, 2017; Mölg and others, 2018; Scherler and others, 2018), forms a layer in a glacier's ablation zone from direct deposits and melt-out of material deposited in the accumulation zone (Evatt and others, 2015; Rowan and others, 2015). Debris clast shapes and degrees of angularity vary widely, and abrupt mineralogy transitions may occur in a glacier's debris layer given the range of rock sources. Exposed supraglacial debris is loosely packed, highly porous, erratically bedded and variable in thickness over decimeter to meter spatial scales (e.g. Reid and others, 2012; Nicholson and Mertes, 2017).

At any location on a debris layer, there is a random distribution of sizes and a lack of sorting beyond that which occurs during supraglacial resedimentation (Lawson, 1979; Menzies and van der Meer, 2017). On the glacier scale, debris generally increases in thickness toward the terminus, as more debris accumulates and resurfaces. In High Mountain Asia, debris layers have been measured to range in thickness from a few millimeters at the up-glacier start of the debris layer to several meters at the glacier terminus (e.g. Ragettli and others, 2015; McCarthy and others, 2017; Rowan and others, 2017; Nicholson and Mertes, 2017). Individual debris clasts can range in size from fine sand to large boulders exceeding 10 m. The dependence of melt on thickness is incorporated into a host of 1-D (e.g. Nicholson and Benn, 2006; Reid and Brock, 2010; Lejeune and others, 2013; Giese and others, 2020) and 2-D (e.g. Reid and others, 2012; Fyffe and others, 2014) glacier melt models. For assessing glacier-scale melt and any of its impacts, debris thickness and its spatial variation must be known.

Historically, debris thickness has been determined by manual excavations (e.g. Zhang and others, 2011). Other approaches include geometrical scaling estimations of exposed debris (Nicholson and Benn, 2012; Nicholson and Mertes, 2017) and calculations from energy-balance models in concert with remote thermal imagery (e.g. Foster and others, 2012; Rounce and McKinney, 2014; Schauwecker and others, 2015) or surface height changes (Ragettli and others, 2015). Empirical relationships between thickness and remotely detected surface temperature have been derived (Mihalcea and others, 2006, 2008; Juen and others, 2014), and debris thickness has also been computed from elevation change and flux divergence by inverting a melt model (Rounce and others, 2018).

Recent studies have employed ground-based ground-penetrating radar (GPR) on supraglacial and englacial debris to determine thickness (e.g. Mackay and others (2014) in Antarctica; Wu and Liu (2012) in China's Tien Shan; Petersen and others (2019) in Wyoming, USA). McCarthy and others (2017) used 225, 450, 900 and 1200 MHz GPR to observe continuous horizons from debris–ice interfaces in 16 of 29 profiles covering ~ 600 m on Lirung Glacier, Nepal. They characterized the radar debris events as high scatter relative to the ice reflections. Nicholson and Mertes (2017) used 200 and 600 MHz dual frequency GPR to validate thickness calculations on Ngozumpa Glacier, Nepal, and Nicholson and others (2018) used the same GPR to collect ~ 3.5 km of radar profiles on Ngozumpa. The ice surface was distinct enough in the majority of the profiles to be picked manually, but a clear debris–ice interface was not always present.

On a larger spatial scale, Huang and others (2017) used 1.27 GHz satellite Synthetic Aperture Radar to isolate volume scattering power caused by debris cover on Koxkar Glacier, Tianshan Mountains, China, which they inverted for debris thickness by decomposing signals into various glacier targets. Huang and others (2017)'s results agreed well with field measurements of the thicknesses of Koxkar Glacier's debris layer, which is primarily composed of gravels and coarse sand with 12% porosity.

Our aim was to determine if a GPR system operating at a bandwidth relevant to remote airborne systems that radiate near 1 GHz can be used to determine glacial debris thickness. Airborne platforms such as drones and low-flying planes achieve extensive spatial coverage but remain relatively unexplored for measuring supraglacial debris. Our study took place on Changri Nup Glacier in the Nepal Himalaya and differs from previous ones because (1) our antennas were elevated slightly (about one free space wavelength); (2) we made direct simultaneous observations of debris thickness and lithology; (3) we collected data in the post-monsoon and did not detect a debris–ice interface with our GPR and (4) we performed laboratory simulations to aid and support our interpretations.

The elevation of our antennas makes our conclusions relevant to a future debris-sensing approach on airborne platforms. The transferability of our findings is not straightforward, as an airborne platform has inherent complications we did not face, such as off-nadir surface returns (i.e. ‘surface clutter’, Holt and others, 2006; Choudhary and others, 2016) interfering with desired on-nadir returns. We do not provide specific recommendations for overcoming challenges on airborne platforms but rather simply assert that our findings may be relevant for developing ~1 GHz airborne systems that can detect the base of supraglacial debris layers.

The fact that we did not detect an ice–debris interface with our ground-based GPR differs from other studies, most likely due to survey timing but also perhaps because of differences in debris lithology, porosity or moisture content. Previous studies show successful 1 GHz GPR signal penetration of dry rounded rocks in the laboratory (Liu and others, 2013) and embedded within wet glacial till near 160 MHz (Arcone and others, 2014). The successful penetration achieved in these studies suggests that our signals passed through the debris–ice interface without significant reflection, and we set up experimental tests to investigate this behavior.

With the initial observation of no distinct, detected debris–ice boundary, we sought a technique to measure debris thickness from GPR data lacking an interface delineation. Backscatter from angular clasts is significant at 960 MHz; accordingly, we hypothesized that volumetric backscatter itself may indicate debris depth. No previous study has directly addressed volumetric scatter of larger-clast debris like that characterizing Changri Nup. To this end, we compared the depth of volumetric backscatter from 960 MHz GPR data with manual ground-truth measurements of debris thickness. We justify our indirect backscatter approach with experimental studies and propose a methodology for determining debris thickness from GPR data that, for any number of reasons, do not show a clear debris bottom.

Study area and methods

Field site: Changri Nup Glacier

The ~ 4 km long Changri Nup Glacier (27.987° N, 86.785° E) in the Nepal Himalaya (Fig. 1) flows from above 5700 m a.s.l. southeast to its terminus at 5240 m a.s.l. (Vincent and others, 2016). It terminates on land short of the Khumbu Glacier in the Mount Everest/Sagarmatha region (Fig. 1 inset). Its lower reaches are covered by a debris layer ~2.3 km long by 0.7 km wide (Vincent and others, 2016).

Fig. 1. GPR profiles A–E taken on Changri Nup Glacier with a 960 MHz antenna unit. The colorbar indicates debris depth determined by volumetric backscatter in GPR profiles. The star indicates the likely source of the debris ridge down the midline of the glacier's accumulation zone. Note the large ice sail intersecting transect B.

Debris on Changri Nup matches the surrounding bedrock in its lithology: granite, pegmatite, gneiss, pelite, calc-silicate and amphibolite, with minerals including K-feldspar, quartz, biotite, muscovite, hornblende, plagioclase and sillimanite (Searle and others, 2003). Changri Nup's ablation zone has a nearly continuous debris mantle punctuated by ice cliffs, surface ponds and a large ice sail (visible in Fig. 1, Evatt and others, 2017). The debris surface contains angular clasts of widely varying sizes (Figs 2a, b). The debris layer generally shows inverse grading (Fig. 2c), with finer particles buried deeper and closer to the debris–ice interface as a result of supraglacial resedimentation (Lawson, 1979; Menzies and van der Meer, 2017). However, larger, monolithic clasts are found nestled in or superimposed on the sandy boulder-gravel (Hambrey and others, 2008) (Fig. 2d). Large rocks insulate underlying ice, causing differential ablation (Fig. 2e). Exposed ice cliffs contain englacial material (Fig. 2f) that was likely incorporated in the accumulation zone.

Fig. 2. Pictures of the felsic granitic and metamorphic debris layer on Changri Nup, with scales marked. (a, b) Clasts are angular and range in size from fine sand to boulders ~ 10 m in diameter, and at any location there is a random distribution of sizes. (c) The debris shows inverse grading, although (d) larger clasts are nestled in or superimposed on finer debris in places. Exposed ice reveals (e) differential ablation and (f) melt-out of englacial debris.

GPR and other field measurements

We used Geophysical Survey Systems, Inc.'s SIR 3000 control unit and Model 3101D ‘900’ MHz antenna unit (pulse dominant wavelength 960 MHz). We collected five GPR profiles on three across-glacier (A, B, C) and two along-glacier (D, E) transects shown in Fig. 1 from 14 to 27 November 2015. This timing, after the monsoon and before the winter cyclones, likely minimized the presence of wet debris. The profiles were spatially distributed over the debris layer of Changri Nup, although along-glacier profiling was limited in favor of collecting across-glacier profiles, especially those on full-width transects A and B. We located transect end points and recorded our paths using a handheld Garmin76 GPS, which has a manufacturer stated accuracy of ±15 m.

We recorded 1024 samples at 16-bit resolution in a 50 ns time range for each trace (although parts of D and E were collected with 100 ns range, apparent in Fig. A2). Logistics of profiling over a very rough surface required recording traces every 10 cm in single point mode; subsequent 10 × stacking provided a 1 m horizontal resolution. As shown in Figure S2, the 960 MHz antenna unit was raised by 27 cm (transect B) or 19 cm (other transects) to place the ground surface in the antennas’ far-field, mitigating complicating effects of a highly variable topography and also attempting (unsuccessfully) to separate the surface reflection from the direct coupling (DC), which is the direct pulse from the transmitter to receiver. One benefit of not raising the antenna unit enough to completely separate the debris surface return from the DC was that closer placement enhanced the chances of detecting a signal from the debris–ice interface. The height of the antenna unit produced a far field spherical wave, the curvature of which approximated a plane when intersecting the surface within the footprint of the antenna directivity. Hereafter, the debris–ice interface is referred to as ‘the interface’; the distance to the interface is equivalent to the debris thickness.

During data collection, the antenna unit operator was always adjacent to the antennas (Fig. S2), but subsequent tests confirmed no interference from body reflections. The face of the antenna unit was kept in approximate alignment with the surface. Although the surface was not flat, most deviations were within ±30° from horizontal. Such deviations and associated changes in polarization relative to transect direction did not impact results because the antenna unit has a 70° 3-dB two-way beamwidth in air regardless of polarization (Arcone and others, 1986).

We made detailed observations of ground debris, including 353 measurements of thickness, by digging to the interface or measuring the height of blocks situated on ice. Sizes varied widely, and debris clasts were generally very jagged and angular. The majority of clasts were less than one 960 MHz wavelength in air (31 cm), while an estimated one-half to two-thirds of debris volume was comprised of rock fragments less than the wavelength in debris (18 cm, calculated with the dielectric properties derived in the rock box experiments). In total, 94% of ground-truth measurements via excavation reached a distinct glacier surface. In the remaining 6% of measurements (indicated by ‘minimum measured z’ in Figures 5, 6, A1 and A2), the debris was frozen, and excavation was prohibitively difficult beyond the recorded depth.

We were unable to determine the debris dielectric constant in the field using a common midpoint survey because we used a radar system where transmitter and receiver were inseparable (see box in Fig. S2), and so we determined the dielectric constant for time-depth conversion subsequently. Post-processing steps in RADAN 7.0 software included resetting time zero to minimize dead time, stacking sequential traces and bandpass filtering to reduce noise, horizontal filtering to remove the clutter of antenna unit reflections, linear range gain to compensate signal loss caused by spherical beam spread, spiking deconvolution to reduce signal duration, and Hilbert magnitude transform to capture returns’ positive-sign energy envelopes (Yilmaz, 1987; Arcone, 1996).

Performing a Hilbert transform on radar data provides a better representation of the total reflected energy than is given by the several wave cycles of the radar wavelet. The Hilbert transform uses a spatial envelope to convert the radar wave time series into magnitude (GSSI, 2011–2017), which is more sensitive to dielectric differences than signed amplitude (i.e. voltage). The Hilbert transform improved interpretability (Zheng and others, 2016) by amplifying the difference between debris and ice, a subtle difference in our raw data.

In our subsequent analysis, we derived a threshold parameter to relate depth-integrated raw energy to debris thickness. Using our ground-truth measurements, we calculated the percent of integrated scattered energy that corresponded to debris depth. We selected this method after first noting the absence of a clear and persistent reflector and then observing that the radargrams had bright and dark areas but no associated, straightforward backscatter drop-off. Integrating under a positive-valued curve is a way to use backscatter as a metric even when scans vary.

By employing the method of leave-one-out cross validation (LOOCV; Arlot and Celisse, 2010), we simultaneously calculated a backscatter drop-off threshold (τ) and assessed the statistics of applying that calculated threshold to the entire dataset (i.e. LOOCV provided τ and uncertainty). LOOCV allows calibration of a model estimator (here, τ) with n observations (here, our physical debris depth measurements). In each of n iterations, the LOOCV model is trained with all but one debris depth measurements, and the depth estimates calculated with the resulting τ are validated against the observation that was left out of the τ calculation. Each ground measurement point, then, gives a separate threshold τ; each τ is compared with the measurement that was left out of its calculation to give its RMS difference.

Rock box experiments

To aid analysis and interpretation of the Changri Nup data, we recorded profiles over a constructed rock box, analogous to Changri Nup Glacier in terms of the clast mineralogy, clast shape and sub-debris interface reflectivity. The reflectivity of an interface between two materials has a strength dependent on the dielectric contrast between the two materials. The reflection coefficient (Γ) at an interface can be calculated from

(1)$$\Gamma = {n_{X} - n_{Y}\over n_{X} + n_{Y}} ,\; $$

where n ($= \sqrt {\it \epsilon }$) is the refractive index and X and Y are the two materials in contact (Rees, 2001). The quantity ε is a material's relative dielectric permittivity, for which in this manuscript we use the terms permittivity and dielectric interchangeably.

In the constructed rock box, freshly broken rocks sourced from a quarry local to Dartmouth College were angular and faceted quartzo-feldspathic gray gneiss, similar to the debris layer on Changri Nup Glacier. The rocks were placed in a trough (Fig. 3) in three clast size groups, small (1–4 cm), medium (5–9 cm) and large (10–20 cm), over 4 cm thick dry pine boards (ε = 1.9–2.0) above densely packed pine shavings (ε = 1.4). The intention was to construct an interface whose reflectivity (Γ = 0.11 for pine-debris) was on the order of (i.e. very low) but greater than that on the glacier to investigate wave behavior at a weakly reflective interface. A greater reflectivity in our constructed experiment setup than on the glacier allowed us to observe processes that had occurred during data collection on Changri Nup but which were not detected (in other words, enhancing them for the purpose of our investigation).

Fig. 3. (a) Annotated photograph of trough with rocks of three clast sizes resting on 4 cm thick dry pine boards over bales of pine shavings. The rock box experiment materials were selected to simulate the reflection magnitude of ice underlying dry gravel debris. Aluminum foil placed beneath the rocks in some of the tests showed the location of the debris bottom and, thus, aided determination of wave velocity. The rocks protruded above the 28.6 cm high board edge by 0 cm (small, b), 1 cm (medium, c), 5 cm (large, d).

We used the rock box to simulate the field conditions of dry angular clasts above a low dielectric substrate (the setup described above) and then above perfectly reflecting aluminum foil to (a) prove rock penetration, (b) examine the 960 MHz GPR's direct coupling duration and (c) measure time delays through the clasts in order to determine debris permittivity. We profiled with both a 960 MHz GPR (the field frequency) and a 2.6 GHz one. We used 2.6 GHz signals to explore how GPR signals interact with clast sizes and debris thicknesses greater than a wavelength. The 960 MHz and 2.6 GHz signals have dominant wavelengths of ~18 and ~ 6.7 cm in debris, respectively. We profiled with both antenna units at a height of ~ 50 cm, which allowed us to separate any surface reflection from the DC, and the strong reflection given by the foil allowed straightforward ε calibration from the time delay through the three rock sections. Using velocity $\equiv c/\sqrt {\epsilon }$, where c is the speed of light, we converted signal return time into depth. The size of the setup limited our ability to vary clast depth. Nevertheless, we were able to vary frequency, rock size and bottom reflectivity to investigate the dielectric constant and test the important factor of debris–ice reflectivity on an analog to the weak debris–ice interface encountered on Changri Nup.

Results and analysis

Changri Nup Glacier

The diffractions and reflections in all of the collected data resemble those in Fig. 4, an arbitrarily chosen 75 m section of transect B. The depth range of 4 m corresponds to a time range of 45 ns for a dielectric constant of 3. The radargram shows raw data, with the only post-processing being a wide background removal filter that removed the direct coupling but was wide enough (301 traces) to preclude removal of any coherent subsurface reflection horizons. In the vertically exaggerated profiles of Figures 5 and 6, raw data are dominated by columns of narrow diffractions, the widths of which show that they are caused by narrow targets.

Fig. 4. A sample 75 m profile along transect B, showing raw data, without stacking or Hilbert transformation. The only post-processing applied is a wide (N = 301) background removal filter that reduces receiver artifacts but leaves diffractions and reflections unaffected. A depth of 0 m is the debris surface. Figure S1 shows details extracted from this profile, highlighting many faint events that are hardly visible in the 75-m-long radargram.

Fig. 5. Top 1 m of the Hilbert-transformed and 10 × stacked GPR data, GPR-based thickness retrievals (red), and ground-truth measurements ($\bullet$ and $\triangle$) along profiles over three across-glacier transects on Changri Nup Glacier. Uncertainty (yellow) is placed above and below a smoothed debris retrieval for ease of interpretation. All three run from climber's right to climber's left across the glacier; the 25 m gap in A indicates a corrupted file, and the 120 m gap in B is colocated with the prominent ice sail. The y-axis depth scale was calculated using ε = 3, and 0 m is the debris surface. Note that the antenna unit was 8 cm higher over transect B. z is thickness, and minimum measured z is above frozen debris that prevented digging to the ice surface.

Fig. 6. Top 1.5 m of the Hilbert-transformed GPR data, GPR-based thickness retrievals (red) and ground-truth measurements ($\bullet$ and $\triangle$) along profiles over two along-glacier transects on Changri Nup Glacier, starting at their down-glacier ends. These profiles are much shorter than those in Fig. 5 and, for visual interpretation, are not stacked here. As in Fig. 5, z is thickness and a depth of 0 m is the debris surface.

In neither the uncompressed nor stacked data are there hyperbolic diffractions sufficiently wide to be interpreted accurately for dielectric permittivity (Yilmaz, 1987; Arcone, 1996), which would not be a simple exercise because the antenna unit was elevated. Figures 5 and 6 show the near-surface of the Hilbert-transformed across-glacier and along-glacier profiles, respectively. All five profiles show the following:

1. (1) absence of clear interface reflections

2. (2) energy concentrated in the near-surface followed by a significant signal decrease (see Fig. S3 for detailed depiction)

3. (3) debris surface reflections partially interfered with by the DC

4. (4) backscatter that is primarily volumetric.

Backscatter originates mostly in the near-surface (Figs 5 and 6; see Figs A1 and A2 for all collected data). Performing a Hilbert transform on the data, using an envelope function to convert the data to a positive instantaneous magnitude, gives the smeared and streaked character to the backscatter in the radargrams. The received signals are dominated by volumetric backscatter from the debris; subsurface returns arrive immediately after the surface reflection, which is weak and partly masked by the DC.

The backscattered returns vary significantly in strength and delay, even over short distances. The absence of a clear debris–ice transition in the radar data was what necessitated investigating whether the radar waves did, in fact, penetrate to and through the interface or were instead attenuated by scattering.

Rock box experiments

Using the experimental set up (Fig. 3), we collected GPR data with and without aluminum foil under the debris because aluminum foil is a very strong reflector. The strong bottom reflection with foil definitively showed that the 960 MHz pulse significantly penetrated all three clast sizes – with minor to no frequency loss or waveform dispersion (Fig. 7) after a round trip in excess of 57 cm (~ 3 in situ wavelengths) – and that more penetration was possible. Experimental tests showed that signal loss could not explain the lack of interface in our field data because energy clearly penetrated the glacial debris.

Fig. 7. (a) A radar profile collected with a raised 960 MHz antenna unit over the rock box in Fig. 3 that had three partitioned clast size sections, under which aluminum foil provided a strong bottom reflection (BR). Total debris depths were 28.6, 29.6 and 33.6 cm over the small debris, medium debris and large debris sections, respectively, and the DC and faint surface reflection (SR) are also marked. (b) The trace marked by the black line in (a), with the same labels as in (a).

The radargram in Fig. 7a shows a surface reflection that is smooth over all sections despite the rough surface, and the waveforms of the distinct backscattered events in the trace replicate the transmitted waveform (Fig. 7b). Locations of surface and bottom reflections were determined precisely from the known trough dimensions.

Raising the antenna 60 cm over the rock box yielded separation of the surface reflection from the DC (see Fig. S3 for more details on how measuring the DC here helped interpret field data). Separating the DC and surface reflection permitted calculation of energy loss with depth from a single scan. The dampening of the signal shown in Fig. 7b can be extrapolated to show that by 6.2 m depth the 960 MHz signal will be completely attenuated. This provides a maximum possible thickness of debris theoretically observable with 960 MHz. Hypothetically, with our method, 6.2 m would give τ = 100%. In practice, a lower τ is necessary because our method makes use of backscatter drop-off. A likely maximum depth for our method, then, would be close to 4.6 m, which corresponds to a τ = 75%.

Unlike the 960 MHz signal, which penetrated all clast sizes, the 2.6 GHz signals underwent scattering loss that precluded any visible returned energy from the foil layer except beneath the small clasts (Fig. S9). The 2.6 GHz signal's difficulty penetrating medium and large clasts was likely because the medium debris compared in size to an in situ wavelength of 6.7 cm and the large clasts exceeded a wavelength in size. This demonstrates the importance of frequency selection; when clast size is comparable to or larger than a radar wavelength, the radar may not penetrate the debris and receive returns for deriving a threshold τ. It follows from the 2.6 GHz wave behavior that the 960 MHz signal would face similar difficulty penetrating debris clast sizes ≥18 cm. Fortunately, most clasts encountered on Changri Nup were <18 cm in maximum dimension (the larger sizes seen in Fig. 2 did not make up the majority of the debris). Field observations of this were corroborated by a comparison of spectra from the rock box and field data. Average reflected spectra from the rock box varied from ~ 770 MHz for the large clasts to 1 GHz for the small ones. Reflected spectra from the field data ranged from 950 to 1000 MHz (Fig. S4), suggesting that the majority of clasts encountered on Changri Nup were smaller than the large ones in the rock box.

By measuring the relative permittivity of mineralogically similar, dry porous debris in the rock box, we found that a lack of dielectric contrast between glacier ice (ε = 3.18) and debris likely accounted for the absence of interface detection on Changri Nup. The relative permittivity of the debris, regardless of clast size, was ~3 (Table 1). The reflection from the debris–ice interface on Changri Nup Glacier, like the reflection from the experimental debris–substrate interface, was coherent and weak because, at ε = 3, the average permittivity of the rock–air–debris matrix is comparable to that of ice. The real part of the permittivity of pure ice at temperatures −18 to 0 °C with 60 MHz–10 GHz electromagnetic waves has long been known to vary insignificantly from 3.18, and the imaginary part is insignificant (Auty and Cole, 1952; Cumming, 1952; Mätzler, 1996; Arcone, 2002). Because the reflectivity of an interface is a function of the difference between the ε values of the two materials (Eqn (1)), this similarity (ε = 3 vs ε = 3.18) suggests that the reflectivity of the debris–ice interface was very weak. The ε similarity, therefore, provides a reasonable explanation for the absence of a continuous ice horizon in our field data.

Table 1. Average refractive index (n) and dielectric constant (ε = n ²) for six randomly chosen measurements in each size classification of the experiment debris

Thickness retrieval: Changri Nup

We could not identify any interface returns in the traces we collected on Changri Nup. Instead, to locate the interface on each trace, we used a proportion (i.e. threshold τ) of the integrated area under the gain-corrected, Hilbert-transformed curve. We used 313 debris thicknesses manually measured on the across-glacier transects to inform a calculation of τ via LOOCV (Arlot and Celisse, 2010), iteratively calculating and evaluating the τ threshold statistic for each measurement. To give the thresholds in Table 2, we averaged the m lowest-RMSE thresholds, where m is 10% of the number of ground-truth measurements. In other words, we determined the integral whose upper bound, on average, best matched the location of the debris–ice interface and used a single threshold computed from the lowest 10% RMSEs. Note that we used only the three long transects (A, B and C) for this computation but applied the resulting threshold to transects D and E, as well.

Table 2. Ground-truth manual debris thickness measurements and GPR-based thickness retrievals for each transect

Given the surface heterogeneity of debris and the fact that debris thicknesses do not follow a normal distribution, the median and interquartile (IQ) range statistics represent the measurement and retrieval datasets better than the mean (μ) and standard deviation (σ); their consistency with the calculated retrievals supports our LOOCV-based threshold. A different threshold is given for transect B because the antenna unit was elevated 8 cm more during data acquisition than it was for the other four transects. A threshold trained using transects A and C was applied to A, C, D and E. The uncertainty ranges are determined by the mean RMS difference in measured – retrieval depths. (1) ‘Total’ includes measurements to frozen debris, but only the measurements to the interface were used to train the threshold values and compute the statistics in column 5. There were gaps in transects A and B (as explained in the caption to Fig. 5) such that data were not collected along the entire length.

Much of Changri Nup's debris is sourced from a large rocky spur generating many rockfalls (star in Fig. 1). We detected debris from this single, dramatic spur buried in the accumulation zone (GPR profiles not shown) during the 2015 field season. Debris incorporated into Changri Nup's accumulation zone subsequently reemerges on the surface as a thick ridge in the glacier's midline to the west of the salient ice sail intersecting transect B in Fig. 1. Profiles D and E traverse this ridge and were excluded from the LOOCV because they are not directly comparable to A–C. Given their location on this ridge, transects D and E cover high concentrations of englacial debris and a potentially gradational ice–debris interface where active englacial melt-out is occurring (see Fig. 2f, larger in Fig. S21, which shows part of transect E). In contrast, the characteristics of transects A–C represent most of the ablation zone on debris-covered glaciers; they cover glacier ice with comparatively few englacial debris inclusions and a distinct debris boundary (Fig. 2d, e).

Training the threshold with LOOCV resulted in a τ of 42% (and 53% for transect B); that is, the depth by which 42% (53%) of the integrated energy over 1024 samples had been scattered matched most closely with the ground-truth points. Over transect B, the antenna unit was elevated 27 cm rather than 19 cm during data collection. Separate thresholds are appropriate because this elevation difference rendered the signals not directly comparable, as gain was applied at different depths. The threshold for B is greater because the deeper returns experience higher post-processing gain, which was applied starting at the antenna. The upper integration bound τ is shown in black on the top 1/4 of an example Hilbert-transformed trace in Fig. 8 (and of the entire trace in the inset). The threshold's numerical value is a function of the number of samples and is larger when fewer than 1024 samples are taken per scan. The red curves in Figures 5 and 6 are thickness retrievals calculated with each transect's corresponding τ. After calculating the thickness retrieval depth relative to the antenna unit, we subtracted the antenna unit height to get interface depth from the surface. The small number of non-physical negative retrievals (Table 2, column 6) were left out of Figures 1, 5 and 6. The yellow curves in Figures 5 and 6 indicate the uncertainty range, which is calculated from the mean RMSE of the LOOCV (12.6 cm or 9.7 cm in Table 2, column 7); the figures show an uncertainty smoothed using a 50 m moving average for visual interpretation.

Fig. 8. An example Hilbert-transformed trace from 298 m along transect B (measured from transect start on climber's right) with the threshold applied. The upper bound of integration, which approximates the bottom of the debris layer, is shown as a black line. Only 256 samples are shown for clarity, but all 1024 are shown in the figure inset.

Discussion

Absence of interface

Our experiments revealed a debris ε of 3, which does not give a dielectric contrast with glacier ice (ε = 3.18). Two factors determine the dielectric constant of the rocky clast debris layer: rock mineralogy and porosity of the air–debris matrix. A reflection is expected at an abrupt change between materials with different dielectric constants. One possible explanation for not detecting an interface on Changri Nup is that the debris to glacier ice transition is a progressive one, with debris clasts incorporated into the ice at the ice surface. Field observations along transects A–C rejected this explanation, and bottom rock geometry is similar to surface geometry, which provides a smooth reflection (see Fig. 7).

In the experiments, we observed strong returns from metal foil placed at the debris base. Therefore, we saw evidence that the waves did penetrate the rocky debris and reach the interface, suggesting that it was the lack of dielectric contrast – rather than loss of signal – that accounted for the absence of a detectable interface in our Changri Nup data. For completeness, we note that signal could have been attenuated significantly if the debris had been wet, but this was not the case as we profiled in the driest time of year with mean temperatures well below freezing.

McCarthy and others (2017) detected a debris–ice interface on 16 of 29 profiles, meaning that nearly half of their recorded profiles lacked an interface. They attributed a missing interface to high scatter or to debris that was too thick or thin. Eleven of their profiles showing an interface were recorded with frequencies lower than 900 MHz, which would decrease scatter but would not be able to measure thin layers. Nicholson and Mertes (2017) and Nicholson and others (2018) profiled Ngozumpa Glacier at 200 and 600 MHz and identified a debris–ice interface throughout their data. None of these studies give a detailed description of debris mineralogy or porosity, and it is certainly plausible that variations from ours in either or both would account for a dielectric contrast at the interface. Another key difference between our study and the aforementioned is the season of data collection. We collected our GPR profiles in the dry post-monsoon season, before winter, whereas McCarthy and others (2017), Nicholson and Mertes (2017) and Nicholson and others (2018) profiled in the spring. In the ablation season, they likely encountered a wet interface, a ‘saturated layer’ of water on ice, below mostly dry debris. A wet interface gives a much stronger reflection than a dry one (Gades and others, 2012).

Ground-truth

The thickness retrievals do not connect the ground-truth points exactly for two main reasons. First, field measurements of debris cover are inherently difficult, and porosity (ϕ) varies markedly over the 2 km of profiles. Second, we used a single value for the dielectric constant to calculate the depth scales on all five transects. Given debris’ heterogeneous porosity and mineral composition, a single value is an average, and some deviations are to be expected. Radar waves are transmitted in a finite, effective, two-way transmit–receive beamwidth (70°); consequently, received signals represent an average of conditions proximal to the antennas. In contrast, the ground-truth depth measurements are point measurements that are not necessarily representative of their surrounding areas. Furthermore, the training points in the LOOCV bias the thickness retrieval. Due to practicality and field constraints, measurements of large rocks were common. We used the experimentally derived ε = 3 for depth-time conversion of thickness retrievals, but solid blocks have much higher dielectric constants (field measurements; Hubbard and others, 1997). Because a wave travels more slowly in a material with a greater ε, we would not expect the thickness retrievals to match depth measurements of blocks exactly.

Backscatter

The scatter dominating our records – both from the field and from the rock box – is volumetric. The rock box results confirmed that there is a surface reflection (Fig. 7). Thus, later returns are indeed from within the debris and are not wide angle surface scatter. Multiple volumetric scattering is not the scattering in our data; it is clear the events are from single scattering because time delays give reasonable and consistent values for dielectric constant in the rock box (Table 1), and the backscattered waveforms in both the field and rock box replicate the transmitted pulse.

Volumetric point backscatter (i.e. simple hyperbolic diffractions) from jagged debris accounts for the general appearance of the field profiles before stacking and Hilbert transform. The shapes of the surface debris are the same as the shapes of underlying debris, although sizes are generally larger on the surface (Fig. 2). Streaks in the transformed profiles indicate energy returned from multiple depths at a single transect location. These streaks are likely due to the combination of diffractions from the edges of angular debris and the Hilbert enveloping of each pulse.

Some returns may have been caused by ice-embedded debris, which could explain events that occur deeper than the thickness retrievals in Fig. 5 and, more frequently, in Fig. 6. Englacial debris, which is more common on transects D and E that run along the debris ridge described in the ‘Thickness retrieval: Changri Nup’ subsection, would affect the calculation of a threshold τ; partly for that reason, measurements taken along transects D and E were excluded from the LOOCV. Any reflectors below the interface would lead to an underestimation of τ and, if that τ were applied to places with less englacial debris, an underestimation of debris thickness, as well. We applied the τ derived on transects A and C to the ridge transects as a way to circumvent a threshold accuracy issue; although applying A and C's threshold to D and E could have caused an overestimation of debris thickness depending on the amount of englacial scattering, this did not happen. The values of thickness retrievals in Fig. 6 are not consistently and significantly greater than measurements. Because the majority of a debris-covered ablation zone has an abrupt and distinct ice interface under which englacial debris is not highly concentrated, the complications posed by englacial scatterers are unlikely, and it is recommended to collect calibration measurements outside the small areas where englacial melt-out of highly-concentrated englacial debris occurs.

Although deeper returns and resonance features exist in our data, most of the signal's near-surface return energy is volume scattering from within the debris layer. Huang and others (2017), too, interpret debris thickness from the volume scattering power of electromagnetic waves (in their case, polarimetric L-band 1.27 GHz synthetic aperture radar). They assumed that volume scattering in the ice is negligible and mathematically derived debris thickness from a coherency matrix and target decomposition. Their approach using a broadband pulse worked only for debris thicknesses <50 cm, whereas our approach does not have such a thickness limitation within the range of likely debris thicknesses.

Leave-one-out cross validation

The Hilbert-transformed curves show that scattering is high in the field debris. The LOOCV provided a way to find the interface based on where the scattering decreases. We used this statistics-based approach of integrating the energy down the waveform because the decrease in scattering is not a sharp transition (i.e. edge-detection and transition-point-detection do not work). The LOOCV suggested that the time by which 42% (53% for transect B) of the integrated energy has been scattered corresponds to the depth of the debris layer. The other 58% (47%) of the 1024 samples used for our data collection comes from the area of the trace that is close to the noise floor (below black line in Fig. 8 inset); therefore, in spite of its large percentage in terms of integrated power, it is inconsequential.

Although the mean retrieval depth is generally less than the mean depth of field measurements, the measured averages are well within the retrieval error range for transect B and for A, C, D and E when considered together. When considering transects A, C, D and E separately, measured depths are still approximately within the uncertainty. The largest difference is for along-glacier transect E, where the average measured thickness exceeds the retrieval's upper error bound. This is unsurprising given the dearth of ground-truth points and measurement bias toward thick, exposed blocks.

The error bounds are large because the variation in measured depths is large; after all, debris thickness varies horizontally on the centimeter scale. Calculating τ separately for transects A and C rather than together yields 41.72% and 41.42%, respectively; this confirms the appropriateness of applying 42% to all radar profiles collected with the 19 cm elevated antenna units.

As evident in Table 2, the standard deviation of debris thickness measurements exceeds the mean in two of five transects. The measurement statistics reflect the highly variable thickness of the supraglacial debris layer; while no retrieval standard deviations exceed the means, thickness retrievals also exhibit variability reflective of the nature of debris.

Uncertainty

The RMSE of the difference between colocated measured depths and retrieval depths using τ = 42% (53% for transect B) gives an uncertainty in calculated thickness retrievals of ±9.7 (12.6) cm. The uncertainty is large because the debris cover is highly heterogeneous, with large variations in thickness over short distances. Applying a single threshold glacier-wide on Changri Nup (τ = 42% for a 960 MHz antenna unit elevated 19 cm with 1024 samples per scan) is appropriate, especially since small variations in τ over the glacier surface do not change the retrieval significantly.

The finding of ε = 3, informed by our rock box experiments, is reasonable for the mineralogy of Changri Nup's debris layer and the observed air content across all five transects, particularly since the permittivity value remains constant over all three clast sizes in the rock box. However, although the laboratory debris was chosen to emulate Changri Nup's in mineralogy and clast angularity, the two inevitably had slightly different electrical properties. The value of the dielectric constant affects the threshold derived through the LOOCV, but any existing, small difference in the dielectric constant was compensated for by training the threshold exclusively with measurements of Changri Nup's debris thickness.

Porosity varies over the surface of a debris layer, and variations in porosity translate to uncertainty. In areas of lower porosity (including large blocks), ε would be greater and the velocity less, which would change the depth scale calculated from return times. Conversely, in more porous areas, the depth scale would expand relative to what is shown with ε = 3. Using a porosity of ϕ = 50% (based on field observations; see photographs in the Supplementary material), the complex refractive index method (Eqn (2); Arcone and others, 2016) gives a refractive index for the solid rock (n _sr) in a debris rock–air matrix of 2.46 and, thus, a relative permittivity ε_sr of 6.07:

(2)$$( 1-\phi_{\rm debris}) n_{\rm sr} + \phi_{\rm debris}n_{\rm air} = n_{\rm debris}$$

where n _air = ε^1/2 = 1 and $\phi _{\rm debris} = 50\percnt$. Equation (2) is rearranged to give

(3)$${n_{\rm debris}-\phi_{\rm debris}n_{\rm air}\over ( 1-\phi_{\rm debris}) } = n_{\rm sr}$$

A debris sample that electromagnetic waves pass through with a 1.4 ns two-way travel time is 25 cm thick with $\phi _{\rm debris} = 50\percnt$ but 21 and 30 cm with $\phi _{\rm debris} = 70\percnt$ and $\phi _{\rm debris} = 30\percnt$, respectively. Such adjustments to the depth scale caused by porosity differences would have an effect on the thresholds in Table 2 but a very minor one. For example, reducing ϕ_debris from 50 to 30% increased the threshold by only 1%.

A value of 6.07 for ε_sr is within the range of most solid feldspars (Hubbard and others, 1997). A change in rocktype to limestone (with ε_{calcite-dolomite} = 9), the plausible end member for higher values of dielectric constant in the Khumbu Himalaya, would give $\epsilon _{\rm debris} = ( 0.5\times \sqrt ( 9) + 0.5\times ( 1) ) ^{2} = 4$ with $\phi _{\rm debris} = 50\percnt$ or 5.7 with $\phi _{\rm debris} = 30\percnt$. Both give a greater dielectric contrast with glacier ice. A reported dielectric constant of 4.03 on the debris cover of Koxkar Glacier (Huang and others, 2017) supports the fact that debris covers vary in their dielectric properties due to variations in porosity, mineralogy and/or moisture across High Mountain Asia.

Glacier-scale considerations

Given debris input distribution and mechanisms of debris transport, debris thickness increases with distance down-glacier, from the clean ice–debris transition to the terminus (e.g. Kirkbride, 2000). Both our measurement means (Table 2, column 5) and thickness retrieval means (Table 2, column 7) show deviation from this glacier-wide trend. Local-scale variability exists everywhere in the debris cover (Nicholson and others, 2018), resulting from not only surface features such as ice cliffs, rivers and ponds but also the location of debris sources. We expect thickness to be greater in the glacier flowlines passing from or through rockfall-prone areas and to be less in the flowlines of large clean areas of the accumulation zone.

In future study, it would be useful to compare debris thickness distribution on different glacier flowlines. A first-order look at the lateral variability along our profiles reveals that A and B do, indeed, have thicker debris along the ridge sourced from the eroding rock spur mentioned in the ‘Thickness retrieval: Changri Nup’ subsection and marked in Fig. 1. Evaluating debris thickness by glacier flowline fits into an important area of future investigation: the spatial distribution of debris thickness across a debris cover. We recommend that future studies incorporate regularly spaced manual depth measurements, which could help inform robust methodologies for interpolating and extrapolating debris thickness values glacier-wide. Due to the logistical difficulty of operating a GPR on the debris surface (see Figs S11–S21), future studies could benefit from distributed ground measurements coupled with radar deployed on drones or low-flying planes, provided that clutter and other airborne operation challenges could be addressed (e.g. source discrimination following Benham and Dowdeswell, 2003; Holt and others, 2006).

Conclusions

This study presents a method to calculate debris thicknesses along five transects of Changri Nup Glacier based upon the depth decay of the volumetric backscatter that dominates our recorded profiles. We pursued this approach because we did not detect a reflection from the ice surface, something we ascribed to the absence of dielectric contrast between glacier ice and Changri Nup debris at the time we measured it. Dielectric contrast is affected by debris mineralogy, porosity and moisture, and our elevated, high-frequency GPR could have detected an interface where dielectric contrast existed. As confirmed by experiments, electromagnetic waves certainly penetrated the debris into the underlying glacier ice.

We explored how GPR signals generally interact with a debris cover of dry, angular, felsic debris: we showed 960 MHz penetration through irregular rocky debris, volumetric backscatter and weak surface scatter despite debris sizes comparable to an in situ wavelength, and a permittivity (ε = 3) relatively constant over a range of clast sizes.

Although a detected debris–ice interface is the most straightforward way to measure debris thickness with GPR, our approach offers an alternative when insufficient dielectric contrast exists. Detection of an interface is most likely when radar measurements are coincident with glacier surface melt; therefore, we recommend that debris GPR campaigns take place in the ablation season (pre-monsoon in High Mountain Asia). When not possible, the approach presented herein offers a method to arrive at debris thickness: calculate a depth measurement-informed threshold under Hilbert-transformed traces.

Our findings also hold potential for studying Martian lobate debris aprons, which soundings from the Shallow Radar on the Mars Reconnaissance Orbiter suggested are comprised mostly of water ice (Plaut and others, 2009) and, thus, are likely to be debris-covered glaciers (Holt and others, 2008). Much of this debris also has ε = 3 (Petersen and others, 2018), and its thickness is poorly constrained due to absence of clear ice–debris interface returns (Baker and Carter, 2019). In combination with point measurements of debris depth, perhaps taken by robot, our method offers a means to calculate debris thicknesses regionally.

The 42% threshold that matches debris thicknesses on Changri Nup Glacier is not totally transferable to other glaciers but, nevertheless, may indicate debris thickness on layers with mineralogy and porosity similar to Changri Nup's. Future study could assess the transferability of the specific threshold to other glaciers. Future research into the variability in dielectric across a debris layer is necessary to quantify the impact of assuming a constant ε. The dielectric constant can vary most significantly (a) when limestone clasts are present in the debris layer or (b) when debris is wet.

We explored the possibility of using GPR to measure supraglacial debris thickness when a dielectric contrast is absent at the debris–ice interface. The volumetric backscatter threshold, calibrated and validated with field measurements of debris thickness, provides a method. Our study lends support to scaling debris surveys by using GPR on airborne platforms because we found that the frequency deployable on such systems (~ 1 GHz) can resolve debris thickness, even in the absence of a melting ice surface.