Introduction
Ice streams and outlet glaciers account for the majority of ice loss from the Antarcticice sheet (e.g. Bamber and others, Reference Bamber, Vaughan and Joughin2000; Rignot and others, Reference Rignot2019), but observational data are limited by a short observational period and inherent difficulty in viewing the subsurface. The primary technique to efficiently probe the glacial subsurface at large spatial scales is airborne radio-echo sounding (e.g. Evans and others, Reference Evans, Dowdeswell and Cofaigh2004; Chu and others, Reference Chu2016; Jordan and others, Reference Jordan2017), which leverages the transparency of ice to specific electromagnetic frequencies. Since the 1960s (Evans and Smith, Reference Evans and Smith1970), airborne radar has been used to image large swaths of the Greenland and Antarctic ice sheets. These airborne surveys record the depth-varying return of electromagnetic waves reflected by englacial reflective layers, interpreted as isochrones (i.e. layers of snow deposited at the same time; e.g. Eisen and others, Reference Eisen2008; Grima and others, Reference Grima, Blankenship, Young and Schroeder2014; Cavitte and others, Reference Cavitte2018). Specifically, transmitted radio waves scatter due to variations in dielectric properties between isochronal layers, generating an image of the ice column.
Isochrones have long provided insights into ice-sheet history, due to the extended timescale that they record (e.g. Hudleston, Reference Hudleston2015; Macgregor and others, Reference Macgregor2016). These layers are passively advected with flow, recording cumulative deformation of the ice mass (e.g. Ng and Conway, Reference Ng and Conway2004; Siegert and others, Reference Siegert, Ross, Corr, Kingslake and Hindmarsh2013; Bingham and others, Reference Bingham2015; Winter and others, Reference Winter2015; Bons and others, Reference Bons2016; Young and others, Reference Young2018). Isochrones have been interpreted as tracers at ice divides to confirm theoretical ice deformation mechanisms in a relatively simple deformational regime (Raymond, Reference Raymond1983; Vaughan and others, Reference Vaughan, Corr, Doake and Waddington1999). Similarly, in instances where radargrams exactly follow ice flowlines, paleo-reconstructions and steady-state imprints of isochrone deformation have been demonstrated (Wolovick and others, Reference Wolovick, Creyts, Buck and Bell2014; Gudlaugsson and others, Reference Gudlaugsson, Humbert, Kleiner, Kohler and Andreassen2016; Wolovick and Creyts, Reference Wolovick and Creyts2016; Holschuh and others, Reference Holschuh, Parizek, Alley and Anandakrishnan2017). However, large deformations in the ice column make it challenging to reliably infer the past flow history from radar sounding interpretations.
A major issue with interpreting isochrones as a proxy for dynamics is the limited coverage of radar transects, compared to the deformational regime in many areas of the ice sheets. At an ice divide where flow is pseudo-2-D, out-of-plane deformation can be neglected and the interpretation is relatively straightforward (e.g. Raymond, Reference Raymond1996; Nereson and others, Reference Nereson, Raymond, Waddington and Jacobel1998; Conway and Rasmussen, Reference Conway and Rasmussen2009). However, in regions where the flow is complex or the radar transects are not aligned with the flow axis, out-of-plane deformation can confound interpretation (e.g. Vieli and others, Reference Vieli, Hindmarsh and Siegert2007; Siegert and others, Reference Siegert, Ross, Corr, Kingslake and Hindmarsh2013; Bons and others, Reference Bons2016). A salient example of complex flow are Antarctic ice streams (i.e. regions of fast ice flow), where our understanding of ice flow limits our understanding of overall ice-sheet mass loss. The influence of basal resistance on isochrone deformation has been shown along flowlines (Wolovick and others, Reference Wolovick, Creyts, Buck and Bell2014; Gudlaugsson and others, Reference Gudlaugsson, Humbert, Kleiner, Kohler and Andreassen2016; Wolovick and Creyts, Reference Wolovick and Creyts2016; Holschuh and others, Reference Holschuh, Parizek, Alley and Anandakrishnan2017), but a method has not been developed to compare radar profiles of arbitrary orientation. Many studies interpret isochrone deformation as paleo-ice flow without methods to quantify the conclusions (Ng and Conway, Reference Ng and Conway2004; Christianson and others, Reference Christianson2013; Siegert and others, Reference Siegert, Ross, Corr, Kingslake and Hindmarsh2013; Keisling and others, Reference Keisling2014; Bingham and others, Reference Bingham2015; Winter and others, Reference Winter2015; Bons and others, Reference Bons2016). Alternatively, theoretical studies have developed isochrone advection as tracers in large-scale ice-sheet models (Born, Reference Born2017; Clarke and others, Reference Clarke, Lhomme and Marshall2005; Hindmarsh and others, Reference Hindmarsh, Leysinger Vieli, Raymond and Gudmundsson2006). An alternative methodology is necessary to confirm to observational interpretations in regions of complex flow. Moreover, the limited coverage of radar surveys motivates resampling of an existing survey along a profile of interest (e.g. a flowline), or in locations where other field data have been collected.
A benchmark setting of particular interest is flow over and around a basal sticky spot. In regions where ice flow is facilitated by a weak, failing bed, sticky spots are areas with anomalously high basal resistance. They were first identified by the warping and generation of streaklines on the ice surface, and subsequently studied as specific examples of basal conditions (e.g. Stokes and others, Reference Stokes, Clark, Lian and Tulaczyk2007; Sergienko and Hulbe, Reference Sergienko and Hulbe2011). Observationally, sticky spots have been imaged from the surface through seismic surveys (e.g.Winberry and others, Reference Winberry, Anandakrishnan, Wiens, Alley and Christianson2011; Luthra and others, Reference Luthra, Anandakrishnan, Winberry, Alley and Holschuh2016, Reference Luthra2017). These anomalous patches provide much higher basal resistance than the surrounding region and therefore are thought to resist the majority of ice flow (Alley, Reference Alley1993; Anandakrishnan and Alley, Reference Anandakrishnan and Alley1994; MacAyeal and others, Reference MacAyeal, Bindschadler and Scambos1995; Joughin and others, Reference Joughin, Tulaczyk, Bindschadler and Price2002; Price and others, Reference Price, Bindschadler, Hulbe and Blankenship2002; Joughin and others, Reference Joughin, MacAyeal and Tulaczyk2004), compensating for up to 13% of the driving stress for some ice streams (Alley, Reference Alley1993). The locations and role of sticky spots are therefore an important consideration in characterizing mass loss from the ice sheets.
Characterizing the flow over a sticky spot is confounded by the relatively short timescale of the observational record. Present-day velocity maps provide a snapshot of the extent and magnitude of fast ice flow in amazing detail (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011, Reference Rignot, Mouginot and Scheuchl2017). Streamlines are the paths of passively advected parcels in the present-day flow field. Alternatively, features called streaklines, or ‘flow stripes’, represent short wavelength surface roughness associated with streaming flow (Bindschadler and Vornberger, Reference Bindschadler and Vornberger1998). Streaklines are thought to record past ice flow history over an irregular bed (Fahnestock and others, Reference Fahnestock, Scambos, Bindschadler and Kvaran2000). The connection between streamlines, streaklines, and internal stratigraphy remains unclear, with studies demonstrating that isochrone folds lack correspondence to surface streaklines (Campbell and others, Reference Campbell, Jacobel, Welch and Pettersson2008).
In this paper, we consider the flow and deformation of isochrones around the central sticky spot of Whillans Ice Stream, shown in Figure 1. This sticky spot is the initiation point of stick-slip events that facilitate the large-scale motion of the ice stream (e.g. Winberry and others, Reference Winberry, Anandakrishnan, Wiens, Alley and Christianson2011; Lipovsky and Dunham, Reference Lipovsky and Dunham2017; Barcheck and others, Reference Barcheck, Tulaczyk, Schwartz, Walter and Winberry2018). Active and passive seismic surveys have been deployed to characterize the location of this stick-slip patch (e.g. Winberry and others, Reference Winberry, Anandakrishnan, Wiens, Alley and Christianson2011; Barcheck and others, Reference Barcheck, Tulaczyk, Schwartz, Walter and Winberry2018) and determine the properties of the bed (e.g. Luthra and others, Reference Luthra, Anandakrishnan, Winberry, Alley and Holschuh2016). However, the timescale over which the sticky spot has been active and the size of the sticky spot remain unknown. Luthra and others (Reference Luthra, Anandakrishnan, Winberry, Alley and Holschuh2016) found that the sticky spot is likely due to a local topographic high in the bedrock, diverting subglacial meltwater from this region. We consider the along-flow deformation of internal layers over the Whillans central sticky spot and develop a method to generate synthetic radargrams along flow. Specifically, we compare radar features between upstream and downstream radar profiles over the Whillans central sticky spot to place observational constraints on flow history. By considering englacial layers, deformation can be inferred throughout the ice column, compared to only observing streaklines at the surface (Fahnestock and others, Reference Fahnestock, Scambos, Bindschadler and Kvaran2000; Hulbe and Fahnestock, Reference Hulbe and Fahnestock2007; Glasser and Gudmundsson, Reference Glasser and Gudmundsson2012; Roberts and others, Reference Roberts, Warner and Treverrow2013).
We demonstrate that the deformation of internal layers suggests flow diversion around the sticky spot, following surface streaklines. The length scale over which this deformation occurs suggests this basal sticky spot has been active for centuries. Based on our synthetic radargrams, we identify distinct features that propagate undisturbed for kilometers along flow. This confirms that a well-lubricated bed supports Whillans Ice Stream and basal slip, u b, dominates the velocity profile. These results are consistent with the use of a Shelfy-Stream model for ice stream trunks (MacAyeal, Reference MacAyeal1989). From the streakline displacements and idealized modeling, we estimate the size of the sticky spot to cover a larger area than the stick-slip patch determined by active and passive surveying. This analysis underscores the complexity of isochrone interpretation, with isochrone rotation and translation occurring across radargrams.
Methods
We develop a methodology to construct synthetic radargrams in regions of fast-flowing ice. This method employs certain assumptions in this implementation, and we refer to it as the Synthetic Advected Radargram Algorithm.
Model assumptions
Regions of fast flowing ice are often facilitated by a weak underlying bed (e.g. Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000). Basal strength inversions for many ice streams indicate that the bed supports nearly none of the driving stress (Joughin and others, Reference Joughin, Tulaczyk, Bindschadler and Price2002, Reference Joughin, MacAyeal and Tulaczyk2004; Cuffey and Patterson, Reference Cuffey and Patterson2010). Therefore, we can consider the analytical solution of ice speed with depth, z,
with basal velocity, u b, englacial deformation, u e, flow-law prefactor, A, flow-law power, n, basal shear stress, τb and ice column height, H (Cuffey and Patterson, Reference Cuffey and Patterson2010). Basal velocity in regions with failing, weak, unconsolidated till is a function of basal yield strength, τc, normal stress, N and pore pressure, p (Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000). For regions of fast flow, basal shear stress is very small (τb → 0) and basal velocity dominates englacial deformation in the ice column (u b ≫ u e). Despite the local influence of sticky spots, fast-flowing ice streams are generally considered to follow the physics of a ‘Shelfy-Stream’ model (MacAyeal, Reference MacAyeal1989; Joughin and others, Reference Joughin, MacAyeal and Tulaczyk2004). If we consider the sticky spot to be a transient feature that slows down most of the ice column as a unit, basal velocity dominates and we can assume that englacial deformation is negligible (u b ≫ u e). We therefore argue that our radargrams are passively transported along flow, uniform with depth. While this assumption may not be valid in all settings, this assumption will be tested in this analysis.
Methodology
To generate synthetic radargrams, we develop a technique to sample an advected radargram along an alternative acquisition line, using the assumption of negligible ice column deformation. As shown in Figure 2a, this method requires one radargram upstream (Profile A), a downstream profile of interest (Profile B) and a hypothesized map-view flow field, u(x, y) (variations with depth, z, are neglected). Profile A can then be passively transported along the flow field by a given spatial step size, $\Delta x_\parallel$, and sampled with bilinear interpolation where Profile A intersects Profile B, shown as magenta dots. The spatial coordinates of Profile A, x, are advected using Euler integration,
with sufficiently small step size, $\Delta x_\parallel$, less than the along track resolution of the Profile A radargram. The sampling of the return power of Profile A onto the synthetic radargram of Profile B is computed with bilinear interpolation of the four nearest values of Profile A. By sampling and storing the interpolated values of Profile A along Profile B, we return a synthetic radargram that would be generated by the hypothesized flow field, shown as the dashed magenta line in Figure 2a. A link to a well-commented version of the code is included in the Acknowledgements for ease of use and modification.
While this technique works conceptually, in practice the along-track resolution of the radargrams is prohibitive for this technique to be computationally feasible over multi-kilometer length scales. Therefore, we introduce efficiency improvements to make the synthetic radargram generation tractable for real applications. The most costly component of this technique is interpolation from the advected radargram to the synthetic radargram at each spatial step. By reducing the number of points to search for the nearest neighbors, we can improve the runtime dramatically. We achieve this by considering only points in Profile A that are a distance, Δx ⊥, or closer to Profile B at a given spatial step (shown as the blue region in Fig. 2a), thereby reducing the search and interpolation times dramatically. To pre-compute the number of computational steps performed during the advection of Profile A, a search along the advected path of Profile A is performed to output the first and last positions where Profile A intersects Profile B. Finally, this algorithm does not require information from prior spatial steps to compute the interpolation at a given spatial step (i.e. an ‘embarrassingly parallel’ algorithm). Therefore, we implement this method in parallel by splitting the interpolation of advected profiles across processors. Each segment of the synthetic radargram is generated separately, and the entire synthetic radargram is concatenated as a final processing step.
Field setting and data
To contrast the behavior of ice at depth to surface observables, we consider the information provided by both velocity streamlines and surface streaklines (Glasser and Gudmundsson, Reference Glasser and Gudmundsson2012). Velocity streamlines (Fig. 1) are derived by tracking the path of a parcel of ice through satellite-derived surface velocity estimates (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011). These streamlines are generated by advecting particles along the modern surface velocity field, again using the Euler integration (Eqn (3)). In the trunk of Whillans Ice Stream, the plug-like flow profile enables a uniform velocity field, leading to parallel, uniform streamlines (Fig. 2b). In contrast, streaklines are an indicator of past flow history, being generated by basal asperities in the paleo-flow field and maintained as surface topographic expressions (Fahnestock and others, Reference Fahnestock, Scambos, Bindschadler and Kvaran2000; Hulbe and Fahnestock, Reference Hulbe and Fahnestock2007; Glasser and Gudmundsson, Reference Glasser and Gudmundsson2012; Roberts and others, Reference Roberts, Warner and Treverrow2013). Streaklines can be easily seen in MODIS imagery of the region (Scambos and others, Reference Scambos, Haran, Fahnestock, Painter and Bohlander2007), and wrap around the sticky spot (Fig. 2c).
In this analysis, we consider two airborne radar profiles collected by the Center for Remote Sensing of Ice Sheets during 2013–2014 field season using the Multi-Channel Coherent Radar Depth Sounder (MCoRDS) radar system (Shi and others, Reference Shi2010). We consider only the high-gain channel of the MCoRDS L1B data product for imaging the deep ice layers at depth. The radar data were collected with the MCoRDS 2 radar system using a transmit power of 1500 W and bandwidth of 180–210 MHz (Shi and others, Reference Shi2010). Specifically, we consider two radar profiles collected over Whillans Ice Stream: one upstream and one downstream of the sticky spot (Fig. 1). The structure of the radargrams suggest a complex flow history imprinted into the isochrones (Fig. 3). The upstream profile (Fig. 3a) contains many regions with steeply dipping layers, characterized by the linear features, which lack clear reflections (Holschuh and others, Reference Holschuh, Christianson and Anandakrishnan2014). Ice advected from the suture zone where Whillans and Mercer Ice Streams merge lacks any clear internal reflectors. This is likely a result of the highly deformed ice distorting the isochrones to such an extent that they become destroyed, and indistinguishable in radar returns. The downstream profile (Fig. 3b) exhibits many complex flow features including large anticlines, synclines, steeply dipping layers and destroyed layers from the Whillans–Mercer suture zone. The interpretation of these features is difficult, and verification of an interpretation is even more challenging. Applying our methodology to these compelling radargrams allows for a quantitative comparison of disparate flow hypotheses over the Whillans central sticky spot.
Results
We apply the Synthetic Advected Radargram Algorithm to sample the advection of Profile A along Profile B, across the central sticky spot of Whillans Ice Stream. First, we consider the case where Profile A follows present-day streamlines, with the synthetic radargram shown in Figure 3c. Comparing the synthetic radargram to the actual radargram across Profile B, we find that the synthetic radargram does not reproduce many of the complex isochrone structural features. It matches the location and width of the suture zone, but the steeply dipping layers occur in the wrong position. Moreover, it does not reproduce the anticline and syncline features of the actual radargram. Second, we consider the case where Profile A is advected downstream along surface streakline paths, with the synthetic radargram shown in Figure 3d. Comparing this synthetic radargram to the actual radargram, shows that many of the complex structural features are reproduced in great detail. This deformational hypothesis recovers the location and widths of the anticline, syncline, steeply dipping layers and suture zone.
Analysis
To better understand these disparate flow hypotheses, we conduct an idealized modeling study of flow over a sticky spot in the Ice Sheet System Model (ISSM) (Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012). The model is an idealized (rectangular) ice stream, described by the Shelfy-Stream Approximation (e.g. Morland, Reference Morland1987; MacAyeal, Reference MacAyeal1989; Bueler and Brown, Reference Bueler and Brown2009), with a central, circular sticky spot of anomalously high friction coefficient (Fig. 4a). Tuning the ratio of the friction coefficient and driving stress reproduces velocities similar to the trunk of Whillans Ice Stream (Fig. 4b). Despite a pronounced deceleration of ice flow over the sticky spot, the ice streamlines are not noticeably skewed by the sticky spot, and rather stay parallel to the driving stress. Alternatively, we can advect particles through the velocity field in the same way that englacial ice would experience advection, analogous to radargram advection (Fig. 4c). The advected particles experience a distinct rotation and wrapping around the sticky spot, similar to the streakline flow hypothesis. This idealized case demonstrates that these two flow hypotheses are not conflicting, instead streaklines store a more nuanced view of ice flow compared to streamlines.
Therefore, this analysis suggests that a paleo-flow field similar to surface streaklines persists at depth. It also demonstrates that the assumption that negligible deformation exists in this part of the ice column is valid. In this region, isochrone features persist over tens of kilometers in the ice stream trunk.
Discussion
In regions of fast flow, synthetic radargrams can reproduce complex features in the ice column, and supply observational constraints for paleo-flow hypotheses. This analysis demonstrates that isochrone advection can provide additional information of flow history at depth, in conjunction with surface streaklines (Fahnestock and others, Reference Fahnestock, Scambos, Bindschadler and Kvaran2000; Hulbe and Fahnestock, Reference Hulbe and Fahnestock2007; Glasser and Gudmundsson, Reference Glasser and Gudmundsson2012; Roberts and others, Reference Roberts, Warner and Treverrow2013). In contrast to the results of Campbell and others (Reference Campbell, Jacobel, Welch and Pettersson2008), these results suggest that isochronal deformation corresponds to the paths defined by surface streaklines. This may be a matter of scale, as this study considers only a portion of the ice stream trunk, not the entire length. This discrepancy cautions the use of this methodology for longer spatial scales. Additionally, the advected and re-sampled basal reflector is not physical, and should not be interpreted as real.
Our analysis also underscores the limitations of isochrone interpretation in regions of fast flow. In regions where radargrams are not aligned with ice flowlines, out-of-plane deformation adds complexity to the interpretation of flow history. Even in this study area, where flow around a sticky spot is relatively straightforward, the observed features (Fig. 3b) are unintuitive without a method to validate interpretation. Synthetic radargrams are therefore an enabling tool to provide testing of flow history hypotheses in regions with multiple radargrams along flow.
This analysis illustrates how large anticline and syncline features can be generated by the out-of-plane advection of features across a radargram (Fig. 5). The steeply dipping layers in Profile A are advected across the plane of Profile B, generating the distinct anticline and syncline features. Figure 5a illustrates that when an upstream feature is advected downstream in a simple (linear) deformational setting, the angle, or obliquity, with which the radar transect crosses flow results in an intuitive expansion of the radar feature. However, in a slightly more complex (harmonic) deformational setting, Figure 5b, the obliquity of the radar transect results in a distinctly different image of a downstream advected radar feature.
The method we have developed will hold for arbitrary 3-D flow fields, u(x, y, z). This includes situations where the bed is not weak and there is internal deformation in the ice column ($u_{\rm b} \not \gg u_{\rm e}$), as in Eqns (1–2). Using a 3-D flow field, the evolution equation described by Eqn (3), still holds for a 3-D implementation. Using depth-varying velocity, the radargram will distort as it is advected through the sampling profile. Sampling of the advected domain may be confounded by the domain passing out of the sampling region. In practice, a suitable 3-D flow field must be generated, such as with a higher-order ice-sheet model (Pattyn, Reference Pattyn2003; Pattyn and others, Reference Pattyn2008). The output of ice-sheet models could be validated by using their velocity fields to generate synthetic radargrams and comparing to observations. As with the 2-D implementation, the bed geometry of the synthetic radargrams is not physical, and not interpretable.
For this application, further quantification of uncertainty would be beneficial to provide a metric to validate models. One possibility is to use methods developed in image processing to quantify the similarity of the synthetic radargram to the observed radargram (e.g. Chopra and others, ). This might provide a more quantitative comparison of flow hypotheses. Additionally, a quantitative similarity metric could be used an objective function in an inverse model to tune an unknown flow field in a region of interest. The inverse approach would allow a model to tune the velocity field in a 3-D sense across the entire ice column, in contrast to hypotheses informed solely by surface observations. Quantitative metrics of image similarity are still in a nascent state of research, therefore this is beyond the current scope of this methodology.
Another application of this method, separate from hypothesis testing, assumes that the flow field in a region is well constrained. In this case, the input flow field is the known quantity, and existing radar surveys can be resampled along profiles of interest. This technique could resample radar surveys along flowlines, improve the coverage of existing surveysor resample radar profiles where other field data were collected to provide additional, direct context for field experiments. This technique would need to be applied carefully, as it might only be appropriate at a local scale. However, this study demonstrates that in regions underlain by flat, weak beds, radargrams can advect in a meaningful way over tens of kilometers. Applying this methodology in regions with more basal heterogeneity would require further analysis. Promising features that could be resampled across surveys are image-obstructing diffractors (Catania and Neumann, Reference Catania and Neumann2010) or steeply dipping layers (Holschuh and others, Reference Holschuh, Christianson and Anandakrishnan2014).
Conclusions
We developed a method to resample radargrams in regions of fast flow, and provided some insights into isochrone deformation at depth. Specifically, our technique allowed us to compare seemingly disparate flow hypotheses across the central sticky spot of Whillans Ice Stream. Our analysis provided a direct visual comparison between the hypothesized flow regimes for this region, by comparing synthetic radargrams to actual radargrams acquired in the region. We found that flow along surface streaklines persists at depth, and provided a better description of past flow history, compared to present-day surface velocity streamlines. Overall, our study provides a quantitative method to test interpretation of isochrone structure, and a potential avenue for resampling existing surveys in regions with well-constrained flow.
Acknowledgments
This research was partially supported by an NSF CAREER award and the George Thompson Postdoctoral Fellowship at Stanford University. The authors thank CReSIS for the radar data analyzed here (Paden and others, Reference Paden, Li, Leuschen, Rodriguez-Morales and Hale2014). The authors also thank anonymous reviewers for comments that greatly improved the manuscript. This code developed in this paper is available open-source under the GNU General Public License, version 3 on the SIGMA research group (SImulations of Geophysical Multi-phAse flows) webpage (https://pangea.stanford.edu/researchgroups/sigma/sigmagitlab). The code is hosted by Stanford University's Center for Computational Earth and Environmental Science (CEES) through the GitLab repository manager.