Automated detection and characterization of Antarctic basal units using radar sounding data: demonstration in Institute Ice Stream, West Antarctica

Abstract

Basal units – visibly distinct englacial structures near the ice-bed interface – warrant investigation for a number of reasons. Many are of unknown composition and origin, characteristics that could provide substantial insight into subglacial processes and ice-sheet history. Their significance, moreover, is not limited to near-bed depths; these units appear to dramatically influence the flow of surrounding ice. In order to enable improved characterization of these features, we develop and apply an algorithm that allows for the automatic detection of basal units. We use a tunable layer-optimized SAR processor to distinguish these structures from the bed, isochronous englacial layers and the ice-sheet surface, presenting a conceptual framework for the use of radio-echo character in the identification of ice-sheet features. We also outline a method by which our processor could be used to place observational constraints on basal units’ configuration, composition and provenance.

Introduction

Processes occurring at and near the base of an ice sheet can have dramatic effects on the dynamics of the ice sheet as a whole. Subglacial hydrology, for instance, directly influences basal motion (Bartholomaus and others, 2008) and plays a role in glacial surface topography (Cooper and others, 2019). Alongside water storage, other basal processes like till deformation affect drag at the base, often to the extent that the modeled velocities of outlet glaciers vary markedly depending on basal parameters (Morlighem and others, 2010; Habermann and others, 2013; Gladstone and others, 2014). Accordingly, basal processes are intimately related to large-scale ice-sheet stability, especially among the marine sections of the West Antarctic Ice Sheet. Basal melt and its implications for ice-sheet mass balance are of particular significance in these regions (Joughin and Alley, 2011; Pollard and others, 2015; DeConto and Pollard, 2016). In West Antarctica and elsewhere, moreover, basal melt and ocean circulation are closely coupled, with important consequences for grounding line retreat and ice-sheet stability (Colleoni and others, 2018; Beckmann and others, 2019). Therefore, as ice sheets are subjected to destabilizing forces from a warming land surface and ocean, it is increasingly important to understand the complex dynamics of processes near the bed.

Ice-penetrating radar sounding is one of the principal tools available for the investigation of those basal processes (Dowdeswell and Evans, 2004). Airborne ice-penetrating radar sounding produces images, called radargrams, of the internal structure and base of an ice sheet. Because the contrasts in crystal fabric and composition among englacial and subglacial features affect the behavior of transmitted radio waves, radar can image those features, which include subglacial water (Siegert and others, 1996; Chu and others, 2016), englacial layers (Fujita and others, 1999; Siegert, 1999; Karlsson and others, 2012) and bed topography and lithology (Schroeder and others, 2014). Typically, the lower 10–30% of these radargrams are devoid of radar reflections, suggesting the structural homogeneity of that segment of the ice sheet and leading to its characterization as an ‘echo-free zone’ (Drewry and Meldrum, 1978; Drews and others, 2009). In recent years, however, due in part to improvements in radar systems and data processing, radar sounding has revealed large features in precisely that near-bed area; notably, basal units have been detected in both Greenland (Bell and others, 2014) and Antarctica (Bell and others, 2011; Wrona and others, 2018).

The compositions and formation mechanisms of these basal units remain uncertain, and it is likely that different, and sometimes simultaneous, mechanisms are at play depending on various physical parameters (Leysinger Vieli and others, 2018). We adopt a broad conception of what may be considered a basal unit, consistent with an array of theorized formation mechanisms that are not necessarily mutually exclusive. One hypothesis is that these units are accretion plumes, which form when meltwater refreezes to the ice-sheet base (Bell and others, 2011, 2014; Leysinger Vieli and others, 2018). In this model, the structures develop through glaciohydraulic supercooling, whereby pressure-related freezing-point depression causes meltwater moving uphill to freeze and accrete. In our analysis, we will differentiate between two modes of basal accretion. On one hand, ice may freeze to the base from below as a single coherent body, separated from overlying ice by virtue of its fabric or some other metamorphic transition; on the other hand, ice may accrete in successive layers, potentially highly deformed or fractured compared to typical isochronous layers. Other hypotheses for the units’ provenance similarly suggest that they are instances of extensive folding near the base, but these models invoke mechanisms other than basal accretion: they involve rheological contrasts within the ice column, variable slip conditions at the base (Wolovick and others, 2014; Wrona and others, 2018) or converging ice-sheet flow in the presence of anisotropy (Bons and others, 2016). Due to the existence of englacial debris in West Antarctica (Winter and others, 2019), we also examine sediment entrainment as a possible mechanism.

Basal units affect both ice-sheet stratigraphy and ice flow, making them worthwhile topics of investigation. Where large-scale deformation of englacial layers is implicated as the formation mechanism (Bons and others, 2016), their relation to ice-sheet stratigraphy is direct. However, both observational and modeling studies have demonstrated that even isolated, accreted basal units can significantly alter the stratigraphy of the surrounding ice column (Bell and others, 2014; Leysinger Vieli and others, 2018), making ice cores difficult to interpret. Of particular relevance to the discussion of basal processes, accreted basal units can also affect ice flow through deformation and softening during the melt and refreezing processes (Bell and others, 2014); these units, therefore, are closely linked to basal dynamics and ice-sheet stability.

Our investigation of Antarctic basal units is twofold: first we present a novel technique for the automatic detection of these units in radargrams, then we use the algorithm output to theorize about the formation and composition of a specific basal unit in Institute Ice Stream, West Antarctica (Fig. 1a). Using the output of our automated detector, we consider the feasibility of several possible formation mechanisms. We introduce this algorithm and discussion as a conceptual framework, serving as a launching point for further echo-character-based analysis of these structures in the future.

Fig. 1. West Antarctic basal structures selected for analysis: (a) Radargram of a segment of survey transect C31a. The most prominent units visible in this radargram are labeled, hereafter referred to as structure 1 and structure 2. These structures are viewed in cross-section, with structure 1 ~ 300 m in height and structure 2 between 100 and 150 m thick. Ross and others (2019) have discussed both of these features in detail. (b) Flight lines (transects) of the British Antarctic Survey's study of Institute and Möller ice streams, West Antarctica. Transect C31a is highlighted in blue, tracing a path from 80°S, 81°W to 83°S, 75°W. A yellow circle denotes the approximate intersection point of the two structures with transect C31a, near 82°S, 77°W. (Both structures stretch over several transects parallel to C31a, and Ross and others (2019) provide a complete mapping.) The grounding line, based on data from the National Snow & Ice Data Center (Bindschadler, 2011), is traced in purple and shading is based on BEDMAP2 surface elevation (Fretwell and others, 2013).

Methods and Data

Processing method

Our approach relies on a novel form of synthetic aperture radar (SAR) processing applied to a set of raw data from the British Antarctic Survey's investigation of the Institute and Möller ice streams, shown in Figure 1 and described in more detail below. Typically, SAR processing allows for finer along-track resolution and an increased signal-to-noise ratio (SNR) through coherent summation of spatially adjacent radar returns (Peters and others, 2007). The disadvantage of this technique, however, is that spatially adjacent pulses are not always in phase. When this is the case, as it is for sloping specular reflectors, coherent summation results in destructive interference and renders those features weaker or invisible on the resulting radargram (Holschuh and others, 2014). This phenomenon is commonly seen wherever englacial layers are steeply sloping; images of such regions are frequently characterized by patches of noise where destructive interference has muted the signal. Castelletti and others (2019) have addressed this problem through a technique called Layer Optimized SAR (LOSAR) processing. LOSAR prevents destructive interference by correcting for the differences in phase between adjacent pulses before summing them. For every pixel, the processor iterates over a range of phase shifts and selects the one that optimizes the SNR of the layer return upon coherent summation. Because the slope of the layer dictates the degree to which adjacent pulses will be out of phase with one another, the optimal phase shift is a direct function of layer slope.

Phase shift response function

As part of its operation, the LOSAR algorithm calculates SNR values as a function of phase shift for all pixels in a radargram. The relationship between SNR and phase shift forms the basis of our automatic detection approach and is hereafter denoted as the ‘phase shift response function’. Plotting the phase shift response function for pixels in specific ice-sheet features demonstrates that each responds differently to applied phase shifts. We illustrate this in Fig. 2 using pixels selected from the bed, an englacial layer and the basal unit marked ‘1’ in Fig. 1a.

Fig. 2. LOSAR processed echo power as a function of applied phase shift for three features: bedrock (green), englacial layer (black), and structure 1 (blue). Thin lines plot the relationships for individual pixels and thick lines plot the average over six pixels for each feature.

The relationships shown in Fig. 2 are consistent with the known radar sounding character of their sources. Layers (black curves in Fig. 2) are highly specular reflectors (Schroeder and others, 2015) and are therefore expected to have a narrowly optimized phase shift range. Because the signal is specularly reflected in one direction, there is one phase shift that significantly increases the SNR of the return; this behavior is exemplified by the small peak width of the phase shift response function for the layer. As mentioned previously, that optimal phase shift is a function of layer slope (Castelletti and others, 2019). The bed (green curves in Fig. 2) is a far more diffuse reflector; it produces returns of significant magnitude across a broad range of directions and is therefore enhanced by a wide range of phase shifts. This is a physical basis for the wide plateau of optimal phase shifts observed in the graph. The basal unit in Institute Ice Stream, with its smaller plateau, appears to behave somewhere between the diffuse scattering bed and specular layers. Phase shift response functions for pixels from structure 1 are the blue curves in Fig. 2. We leverage these feature-dependent differences in phase shift response functions during our automated classification approach as well as during our analysis of the features themselves.

Classification algorithm

These distinctive phase shift response functions allow us to develop a technique that classifies feature types. After eliminating the surface and the bed, our approach assigns model phase shift response functions to the two remaining feature types: isochronous internal layers and potential basal units. We then generate the phase shift response functions for every pixel, examining its unique graph and matching it to the model function that it most resembles – thereby designating it as part of a basal unit or part of a layer.

In order to remove the surface and bed echoes, we begin by separating the radargram into an upper and lower half. We select the brightest pixels per range line in each; these can reasonably be expected to be part of the ice surface and bed, respectively. (This is possible after the removal of the brightest return from the flight line). We then perform a linear interpolation between the selected bright pixels in the upper half to pick the surface and a cubic interpolation between bright pixels in the lower half to pick the bed (Fig. 3). After interpolating, we apply a simple region-growing algorithm to properly remove the entire diffuse reflection from the bed. With the surface and bed removed, the approach focuses on the interior of the ice sheet.

Fig. 3. Ice surface and bed. The ice-sheet surface is traced using a linear interpolation (shown in yellow) and the bed is traced using a cubic interpolation (shown in red). Both interpolations were performed using brightness threshold picks, which can be seen scattered around the respective curves.

In the absence of the surface and bed, any remaining features are assumed to be either layers or potential basal units. We distinguish between these using their characteristic phase shift response functions. The first differentiation test we use involves the conversion of each phase shift response function to a probability density function, plotting the frequency of power values achieved across the entire range of phase shifts applied. These probability density functions will look different for different features, likely as a result of the differences in ‘tailedness’ among phase shift response functions for various features; this tailedness is similar to abruptness in character (Oswald and Gogineni, 2008; Young and others, 2016), but applied to the azimuth response, like specularity. We apply a Gaussian fit to the model probability density functions for a layer and a basal unit, produced by averaging over a collection of pixels from each feature. Finally, we conduct the comparison, generating this power distribution function from every pixel and applying a Gaussian fit. This allows us to compare the standard deviation of the Gaussian from every pixel to the standard deviations of the model layer and basal unit Gaussians. Pixels are thereby sorted on a binary basis, into either the layer group or the basal unit group.

We end by conducting a short secondary test, re-evaluating only those pixels that have been classified as part of a basal unit. Specifically, we reclassify as layers any pixels whose phase shift response functions are not centered near zero. We consistently observe that the phase shift response functions for pixels in the basal unit are centered around zero (as in Fig. 2) likely because of the presence of many small specular segments of varying slopes (explained in depth in the discussion). Of course, the phase shift response functions for a layer can be centered at zero, too; this occurs when the layer is flat and no phase shift is necessary to correct for the slope. This implies that possessing a phase shift response function centered at zero is necessary but not sufficient to classify a pixel as a basal unit; for this reason, it is possible that pixels in flat layers could fail the second test (examples of which can be seen in Fig. 4). Figure 4 does show, however, reasonable success in characterizing flat layers in general, suggesting that most of these pixels have already been correctly classified as layers by the preliminary test and are not subject to the second test.

Fig. 4. Output of automatic feature classification algorithm: (a) The original radargram of the Institute Ice Stream basal unit, the algorithm's input. (b) Classification of pixels by the algorithm. Noise, ice surface and bed are eliminated first, and remaining pixels are sorted into basal unit or layer classes afterwards.

Raw data

The radar sounding data we analyze were gathered in an airborne geophysical investigation of the Institute and Möller ice streams (Fig. 1b), which are located in the Weddell Sea sector of the West Antarctic ice sheet and whose drainage accounts for 20% of the ice sheet by area (Rose and others, 2015). The campaign was conducted by the British Antarctic Survey during the 2010–2011 field season and the radar sounding data were collected by the Polarimetric-radar Airborne Science INstrument (PASIN) (Siegert and others, 2016; Jeofry and others, 2018).

Results

Figure 4 shows the output of this automatic feature classification algorithm, with the ice-sheet feature indicated by pixel color. As is visible from Figure 4b, the classifier has processed raw ice-penetrating radar data and produced a workable map of a basal unit in this region of Institute Ice Stream, including its complex morphology. (We remark on structures 1 and 2 separately in the discussion.) There are, to be sure, areas that seem prone to inaccurate categorization, as demonstrated by the scattered yellow pixels above structure 2 and the band of yellow near the surface. With regard to the pixels above structure 2, it is important to note that these mid-depth layers have also undergone significant deformation, and so may bear some resemblance to structure 2 in their echo character without having achieved the same level of stark stratigraphic contrast. Because these pixels are located primarily within locally flat segments, however, we believe they are likely true misclassifications, possibly slipping through the second test because their phase shift response functions are centered at zero.

The classification of ice near the surface, however, may be more grounded in physical properties; while this misclassification provides an opportunity for improvement if the algorithm is to be used purely to find basal units, it seems also to indicate that there is some real difference in the shallow subsurface, as has been previously shown (Grima and others, 2014a). This suggests another possible use for the processor as an indicator of heterogeneity in the ice sheet that is not directly identifiable from the radargram itself.

We emphasize that what we provide here is for now a conceptual framework, demonstrating the feasibility of phase shift response as a classification mechanism. We have tested it on a limited dataset, which we present as a case study. Future work could be aimed at eliminating empirical thresholds and improving the model functions used, making the processor location-general. In other words, phase shift response could serve as a starting point for future image processing or machine learning algorithms. As we have remarked, phase shift response is grounded in physical characteristics of the source feature; this implies that the fundamental technique should theoretically be applicable regardless of ice-sheet conditions.

Discussion

Derivation of the scattering function

Here, we focus specifically on the potential basal unit in Institute Ice Stream, which has been noted elsewhere by Bingham and others (2015), Winter and others (2015) and Ross and others (2019). The radar echo character of an ice-sheet component – that is, the manner in which it interacts with transmitted radar pulses – is an expression of its wavelength-scale geometry (Campbell, 2007; Grima and others, 2014b). Therefore, the echo character of the basal unit allows us to place constraints on its structure and formation. One graphical representation of a feature's echo character is its scattering function, which describes how returned signal power depends on the radar look angle (Ulaby and others, 2014).

The following derivation indicates that there is in fact a one-to-one relation between the phase shift response function and the scattering function of a feature. As detailed above, the phase shift response function describes, for a given pixel, the dependence of signal power on phase shift. Consistent with the demonstration by Holschuh and others (2014) that power loss during SAR processing depends on layer slope, Castelletti and others (2019) have shown that there exists a direct relationship between phase shift and optimal layer slope:

(1)$$\Delta\varphi = -{4\pi f\Delta x n_{\rm i}\sin{\theta}\over c}$$

where Δφ is the phase shift between adjacent range lines that optimizes SNR, f is central frequency, Δx is distance between range lines, n _i is ice's refractive index, c is the speed of light, and θ is the local slope of the layer, expressed as an angle to the horizontal.

The slope of a layer can be considered analogous to look angle; changing the layer slope with respect to an instantaneous aircraft position produces the same effect as changing the aircraft position with respect to a fixed layer. This conversion of phase shift to incidence angle provides a means of converting between phase shift response function and scattering function. Finally, because the relationship between layer slope and phase shift is monotonic, those two functions will have the same shape, scaled to different axes. Phase shift response function is therefore a good proxy for the scattering function of an observed englacial feature.

Echo character of Institute Ice Stream structures

Having established a one-to-one relationship between a feature's phase shift response function and its scattering function, we can analyze the basal unit's phase shift response function as a representation of its echo character. We individually address four possible compositions for this basal structure: typical isochronous layers, entrained sediment, an isolated body accreted to the base and heavily deformed or fractured layers.

Typical isochronous layers

As is visible in Figure 2, the phase shift response function (and, by the previous derivation, the scattering function) of a regular isochronous layer is distinctive. As would be expected for a specular reflector, a layer's scattering function possesses one clear peak, as it reflects at precisely one angle. Structures 1 and 2, as demonstrated by Figure 5, exhibit scattering functions that are easily distinguishable from that of an archetypal isochronous layer. Both structures have scattering functions with wide plateaus, indicating some disruption in the consistent specularity that characterizes overlying layers. Furthermore, these plateaus are visible not only in the average scattering function for the feature, but in the scattering functions of individual pixels, too, shown by the thinner lines in Figure 5. This seems to indicate some fundamental departure from typical layer specularity in the geometry of the basal unit.

Fig. 5. LOSAR processed echo power as a function of applied phase shift for the two sections of the basal unit, as numbered in Figure 1: structure 1 (blue) and structure 2 (green). Thin lines plot the relationships for individual pixels and thick lines plot the average over six pixels for each feature.

Entrained sediment

Given the presence of englacial debris in West Antarctica (Winter and others, 2019), we also consider the possibility that the basal unit represents some pocket of entrained sediment. The cumulative effect of a collection of point scatterers is to produce similar reflections regardless of the incidence angle (Aglyamov and others, 2017). On average, therefore, returned power should be approximately constant over the entire range of angles, producing a flat scattering function. The scattering functions for both structures have clear maxima, as shown in Figure 5, so we find that sediment is an unlikely composition for either of these specific features.

Single accreted body

Basal freeze-on has typically been expected to produce units containing both the refrozen ice itself as well as deformed internal layers (Bell and others, 2014; Leysinger Vieli and others, 2018). We first consider the boundary case, in which the basal unit is a single coherent body of ice accreted at the base, distinguishable from overlying ice by virtue of a difference in crystalline fabric or some other metamorphic transition. The basis of radio echo sounding, however, is the existence of dielectric contrasts within the ice column (Siegert, 1999); as a result, we would expect echoes from the surface but not the interior of a homogeneous body, because that body would not contain the requisite contrasts to produce a reflection. We do, however, see echoes from within both basal structures; a collection of the individual phase shift response functions in Figure 5 is produced by pixels in the structures’ interiors. Consequently, these structures seem unlikely to have formed by the accretion of a single coherent body onto the base of the ice column.

Deformation and folding

A number of processes have been shown to produce large-scale deformation of internal layers near the base. These include mechanisms that do not invoke basal freeze-on and therefore would suggest that the basal structures are entirely composed of previously linear layers that have since been heavily deformed. Such mechanisms include moving sticky and slippery patches (Wolovick and others, 2014) and convergent flow coupled with anisotropy (Bons and others, 2016).

Of course, basal freeze-on generates large-scale folding as well, as has been shown through both modeling studies (Leysinger Vieli and others, 2018) and radar observations (Bell and others, 2014). In these models, basal structures might be composed of a combination of accreted ice and heavily deformed layers.

Either way, the scattering functions shown in Figure 5 seem indicative of the presence of folded and deformed layers in some capacity. As stated above, and further exemplified in Figure 2, the basal unit scattering functions are not characteristic of perfectly specular or perfectly diffuse reflectors. Instead, they appear to represent an echo character that is in some sense intermediate, as might occur in the presence of deformed layers. In particular, the functions could be explained by some process that ‘shortens’ the specular interface. Major deformation could produce that effective shortening by folding layers to such an extent that what was once a long interface would become a collection of specular segments of varying slope. This mechanism seems particularly reasonable for structure 2, the distance of which from the base makes accretion of refrozen meltwater unlikely; it likely began as typical englacial layers that were deformed, or possibly fragmented, to such a degree that they no longer resemble other layers in their echo character.

As mentioned, such massive folding can be a product of a number of dynamic processes. We believe this offers another application of our algorithm and of the phase shift response function in particular. The phase shift response function, as demonstrated, provides a way to approximate a feature's scattering function. Process-specific models, like those conducted by Leysinger Vieli and others (2018), may produce process-specific scattering signatures. Future work could compare the scattering functions we expect from models of various processes to the actual scattering functions that this algorithm allows us to estimate. Different processes, for example, may produce scattering functions with different depth dependencies, allowing us to use the phase shift response function to distinguish between various folding and deformation mechanisms.

Conclusion

With their implications for basal processes and ice-sheet morphology, basal units remain a subject worthy of investigation. Here, we provide an algorithm aimed at facilitating the analysis of their formation mechanisms and compositions. The framework we present here is, in theory, applicable to any coherent radar sounding data. Applied at scale, it could eliminate the need to manually identify and analyze basal units in large radar sounding datasets. As a case study, we applied this technique to a prominent basal unit in Institute Ice Stream. In particular, the algorithm allows us to approximate the features’ scattering functions and echo character. In this case, the structures seem to be the result of large-scale deformation and folding, the cause of which remains to be investigated. That investigation might take the form of a comparison between observed and modeled scattering functions; moreover, the same approach could be used to investigate any number of hypotheses as to a unit's formation mechanism and composition.