Basal conditions and ice dynamics inferred from radar-derived internal stratigraphy of the northeast Greenland ice stream

Abstract

We analyze the internal stratigraphy in radio-echo sounding data of the northeast Greenland ice stream to infer past and present ice dynamics. In the upper reaches of the ice stream, we propose that shear-margin steady-state folds in internal reflecting horizons (IRHs) form due to the influence of ice flow over spatially varying basal lubrication. IRHs are generally lower in the ice stream than outside, likely because of greater basal melting in the ice stream from enhanced geothermal flux and heat of sliding. Strain-rate modeling of IRHs deposited during the Holocene indicates no recent major changes in ice-stream vigor or extent in this region. Downstream of our survey, IRHs are disrupted as the ice flows into a prominent overdeepening. When combined with additional data from other studies, these data suggest that upstream portions of the ice stream are controlled by variations in basal lubrication whereas downstream portions are confined by basal topography.

Introduction

Ice streams connect ice-sheet interiors to their coastal margins, where mass loss occurs through melting or calving. The link between ice-sheet mass balance and ice-stream dynamics makes understanding these features critical to predicting future ice-sheet behavior. The vigor and extent of streaming flow depend on basal processes involving subglacial topography, sediment, hydrology and underlying geology (Alley and others, 1986, 2007; Blankenship and others, 1986; Anandakrishnan and others, 1998, 2007; Bell and others, 1998; Studinger and others, 2001). Direct measurement of basal properties on a basin-wide scale is not feasible, but isochronous radar reflectors (internal reflecting horizons, or IRHs) retain information about past flow (Ng and Conway, 2004; Siegert and others, 2004), and can also deform in response to current basal conditions (Catania and others, 2010; Christianson and others, 2013). Thus IRHs are a potential source of information about both paleo and modern ice dynamics (Whillans, 1976), and the basal environment. Here we use radio-echo sounding (RES) to examine deformation of IRHs in the northeast Greenland ice stream (NEGIS).

NEGIS is uniquely situated to affect ice-sheet mass balance due to its large catchment, great inland extent in comparison to other fast-flowing ice in Greenland (Fahnestock and others, 1993; Joughin and others, 2010) and termination in three large, marine outlet glaciers. Yet NEGIS differs from other ice streams in important respects, lacking topographically confined dendritic tributaries (Rignot and others, 2011) and gradually widening along-flow by incorporation of ice through the shear margins (Fahnestock and others, 2001a; Joughin and others, 2001). The onset of streaming flow is attributed to subglacial water production from vigorous basal melt due to a local geothermal anomaly (Fahnestock and others, 2001b), which requires a geothermal heat flux 20–30 times higher than that generally expected in the old continental rocks thought to underlie much of Greenland (Henriksen and others, 2000). The inference of high geothermal flux is consistent with the direct observation of a thawed bed at the nearby NorthGRIP (NGRIP) ice-core site (Dahl-Jensen and others, 2003; Buchardt and Dahl-Jensen, 2007) and is physically plausible based on ice-sheet modeling studies (Greve, 2005), but the amplitude of the locally high geothermal flux is uncertain.

The ability of the ice stream to transmit coastal forcing to the ice-sheet interior is even less constrained. Although NEGIS has likely been a steady-state feature of the ice sheet (Joughin and others, 2001), recent coastal acceleration and thinning within overdeepened fjord troughs raise questions about its stability (Pritchard and others, 2009; Joughin and others, 2010; Sasgen and others, 2012; Bamber and others, 2013; Khan and others, 2014). The presence of dilatant till (K. Christianson, unpublished information) and lack of a confining basal topographic trough (Joughin and others, 2001) under the central portion of the ice stream indicate that ice-sheet marginal thinning could be rapidly transmitted upstream via changes in basal hydrology. Without more comprehensive study of the transition between coastal and interior flow regimes, as well as the impact of the basal environment on ice dynamics at NEGIS, we cannot assess whether ongoing mass-loss trends (Khan and others, 2014) are a prelude to a more-dramatic, sustained retreat.

Here we use the internal stratigraphy revealed by radio-echo sounding to examine the regional setting of NEGIS and its flow history. Firstly, we use the continuity of IRHs observed in airborne RES data (Karlsson and others, 2012) to guide inferences about regional ice-flow dynamics across northeast Greenland. Secondly, we consider the usefulness of extending simple strain-rate models previously applied in similar areas (Fahnestock and others, 2001b; Dahl-Jensen and others, 2003) for determining local variations in accumulation and basal-melt rate within a high-resolution survey conducted in the upper reaches of the ice stream. Our results provide insight into the past and present stability of NEGIS, and allow us to identify areas prone to dynamic change as a result of ice-sheet marginal forcing.

Data

Airborne RES data

We use airborne radar data collected in May 1999 over northeast Greenland (Gogineni and others, 2001) to calculate IRH continuity (Karlsson and others, 2012). Although additional data exist, we limit our analysis to 1999 for three reasons. Firstly, the May 1999 field data provide spatially comprehensive coverage of this area (Fig. 1a). Secondly, these data were acquired using a P-3B turbo-prop aircraft, which flies lower and slower than the DC8 used in more-recent missions, allowing better resolution of IRHs. Thirdly, we wish to minimize the number of radar systems used in our analysis, to avoid changes in IRH continuity index arising from changes in instrumentation rather than changes in the ice sheet (Karlsson and others, 2012). These radar data are geolocated with dual-frequency GPS data, and provide thickness accuracy of ±10m for >90% of the survey (Gogineni and others, 2001). Typical trace spacing is 150200 m (Gogineni and others, 2001).

Fig. 1. Location of NEGIS radar survey. (a) Location of ground-based radar survey (white box) and aerial radar profiles (brown lines) at NEGIS. Inset of Greenland indicates area of panels (a) (black box) and (b) (black dot). Airborne radar profiles BB'-DD' are shown in Figure 3. Positions of NGRIP (NorthGRIP Members, 2004) and NEGIS (Vallelonga and others, 2014) ice cores are annotated (black dots). (b) Zoom of white box in (a), including NEGIS radar profiles (brown lines). Synthetic radar profiles along flowlines aa'–cc' (dotted orange lines) are shown in Figure 5. Black dot indicates the location of the NEGIS ice core. Ice flow is from lower left to upper right. Velocities are derived from interferometric synthetic aperture radar (InSAR) data (Joughin and others, 2010). Note that the northwestern shear margin steps outward just upstream of our survey, visible most prominently in the 50ma^-1 velocity contour. Projection here and in all figures is polar stereographic with a central meridian at 45° W and latitude of true scale at 70° N. Elevation is relative to World Geodetic System 1984 ellipsoid.

Ground-based surveys

During summer 2012, we collected -350 line kilometers of ground-based RES and dual-frequency GPS data -140 km downstream of the onset of streaming flow at NEGIS (Fig. 1). In our study region, the ice stream is -25 km wide, and reaches a peak velocity of 65 ma-¹. This field site was selected as the farthest downstream location safe for ground surveys; widespread surface crevassing appears down-flow as ice speed increases. Radar data were generally acquired along profiles orthogonal to flow, crossing both shear margins and extending an additional 5-10km into the slow-moving ice flanking the ice stream (Fig. 1b).

Radar data were collected using a monopulse system with a center frequency of - 3 MHz (Welch and Jacobel, 2003; Welch and others, 2009). Each trace is a stack of 2000 radar pulses. Typical raw trace spacing is - 8 m. Processing followed Christianson and others (2012), and included bandpass filtering, correction for antenna separation, interpolation to a constant trace spacing (8 m), and time-wavenumber migration (assuming constant radar-wave velocity of 169 m µs^-1). Dual-frequency GPS data were collected concurrently with RES data, and processed kinematically using differential carrier-phase positioning (Chen, 1998). Position uncertainties for each trace are ±0.05 m in the horizontal dimensions and ±0.1 m in the vertical dimension. Processed RES data have a vertical precision of - 3 m. In all profiles, we clearly imaged the basal interface and numerous continuous IRHs to within 50 m of the bed except in areas of steeply dipping IRH folds in the ice-stream shear margins (Fig. 4, further below). IRHs were automatically digitized by an algorithm that traced three wavelet half-phases centered on the local amplitude maxima chosen at any point along a profile by an observer (Gades and others, 2000). Surface, IRH and bed-elevation datasets assembled from the RES and GPS data were gridded to 100 m posting using an inverse-distance weighted nearest-neighbor surfacing algorithm and used to calculate glaciostatic hydropotential. Interpolation artifacts were removed using adaptive noise removal, in which the noise is estimated as the mean local variance in a 0.09 km² window (Lim, 1990). The final surfaces are not forced to match the available data, to allow accurate noise estimates in data-sparse regions. The effects of these choices vary between datasets. For ice thickness, gridded and directly measured data agree to within ±10 m. However, for IRHs in the regions of short-wavelength features (i.e. the steeply pitching IRHs of greatest interest to our study), the gridded data underestimate the true IRH fold amplitude by up to 50 m. We draw our interpretations in light of these limitations.

Methods

Internal reflecting horizons

Ice-penetrating radar waves reflect off interfaces (IRHs) that mark dielectric contrasts in the ice column. Dielectric contrast generally results from changes in ice chemistry due to variations in ionic content of atmospheric fallout (Matsuoka and others, 2010). An IRH formed in this manner represents a single moment in time, i.e. each IRH is isochronal (Whillans, 1976; Siegert, 1999; Kanagaratnam and others, 2001). Post-depositional deformation produces folds in this isochronal stratigraphy. Here we use the spatial pattern of IRH folds to assess ice-stream spatio-temporal stability (Whillans, 1976). After establishing that folds in upstream portions of the ice stream are stable and likely form in situ, we use a suite of one-dimensional (1-D) strain-rate models to infer basal conditions consistent with the observed IRH deformation patterns.

Radar continuity index

Under ideal conditions, IRHs observed in a single radar sounding record the history of fallout events. However, the number and strength of IRHs can be reduced by any process that disrupts them (e.g. folding, extensive crevassing, and melting) or that reduces returned radar power from them (e.g. surface crevassing, increased radar-wave attenuation due to warmer ice, or large IRH slopes that scatter energy). Continuous IRHs demonstrate that such processes are weak or absent, placing upper limits on estimates of spatial and temporal changes in ice flow. Tracing IRHs across large regions is difficult, and does not easily yield a quantitative measure of continuity. Instead, we follow Karlsson and others (2009, 2012) and Sime and others (2014) and use the IRH continuity index, a single, readily calculated value for each sounding that increases with the number and strength of IRHs detected. Discontinuity in IRH number or strength will produce easily detected spatial changes in .

is defined as the mean of the absolute value of the derivative of the power of each radar trace within a time window (n ²–n ₁) chosen to eliminate effects from the surface and the basal interface, or

(1)

where P is the reflected power (dB) for each time sample i, Ar is the one-way range (depth if calculated over the entire ice column) and N is the number of samples. Since we are interested primarily in changes in, and not its absolute value, Ar is simply a scaling factor, and is set to 1 for convenience (Karlsson and others, 2012). We calculate using only airborne RES profiles, and average IRH continuity over 20 traces (~3–4km width dependent on the trace spacing), which generally corresponds to 1-3 ice thicknesses.

Strain-rate modeling

Internal layers record the history of ice flow, and of mass input and loss through snow accumulation and basal melting (Weertman, 1976; Whillans, 1976; Whillans and Johnsen, 1983; Waddington and others, 2007; Koutnik and Wadding-ton, 2012). Fully simulating the complete stratigraphic history of an ice sheet is frequently an underdetermined problem, but requires at least a fully coupled higher-order model of ice and tracer dynamics constrained by a range of ice-core, geophysical and climate data (Clarke and others, 2005; Lhomme and others, 2005; Waddington and others, 2007; Koutnik and Waddington, 2012). Recognizing the great computational and observational difficulties of solving this problem fully, a wealth of information has been extracted from simplified one- and two-dimensional treatments of ice flow (e.g. Dansgaard and Johnsen, 1969; Fahnestock and others, 2001b,c; Dahl-Jensen and others, 2003; Waddington and others, 2007; Christianson and others, 2013). Here, following Fahnestock and others (2001b), we apply a 1-D approach to our high-resolution survey of NEGIS, and obtain useful but necessarily limited results.

Fahnestock and others (2001b) used 1-D vertical strain-rate models to estimate spatial variations in time-averaged long-term (~8 thousand years ago (ka)) accumulation and basal-melt rate across north-central Greenland during the Holocene. Modeled age-depth relationships were then compared with those derived from ice cores to reveal spatial changes in accumulation and basal-melt rate (Fahnestock and others (2001b)). Physically, an increase in basal melting drops all layers above equally; an increase in accumulation rate and the resulting changes in ice flow have a depth-dependent effect on the depths of IRHs, limited by the restriction that ice cannot flow through the bed. Since the effects of surface accumulation and basal melt on the relative spacing of layers differ, this modeling approach provides at least the sign of accumulation rate and basal-melt rate anomalies, and some guidance on magnitudes, without fully constraining time evolution of either.

1-D strain-rate model

Following Fahnestock and others (2001b), the simplest ice-sheet strain-rate model directly applies the definition of strain rate to the entire ice column: vertical strain rate is equal to the change in length (accumulation rate) per unit length (ice thickness) (Nye, 1951, 1957, 1959, 1963). This is a ‘plug-flow’ model, in which the horizontal velocity and vertical strain rate are depth-independent. Fahnestock and others (2001b) modified this model to include basal melt (termed the Nye-melt model). If the ice sheet is frozen to the bed, the assumption of depth-independent vertical strain in the Nye and Nye-melt models is invalid. For these areas, Dansgaard and Johnsen (1969) introduced a model in which the vertical strain rate is constant from the surface to some height z = h above the ice-sheet bed, and from there decreases linearly to zero at the ice-sheet base. Following the terminology of Fahnestock and others (2001b), we refer to this section of the ice column (between z = 0 and z = h) as the ‘shear layer’. Further details of these models can be found in the literature (Dansgaard and Johnsen, 1969; Fahnestock and others (2001b)). In our analysis, we minimized the squared difference f(λ_h, b)² between calculated ages and ice-core ages of IRHs using nonlinear ordinary least-squares fitting (Fahnestock and others, 2001c):

(2)

where t_m is modeled IRH age, λ _h is accumulation rate, b is either shear-layer thickness or basal-melt rate, depending on the model used, d(i) is depth of IRH i, t(i) is the age for IRH i, a is the age uncertainty of IRH i, and n is the number of dated IRHs. Here we set a = 1 to weight all layers equally, as we assumed the depth of the IRH of age i to be absolute and did not attempt to translate the age-depth uncertainty between different radar data (airborne (150 MHz) vs ground (3 MHz)).

As noted by Fahnestock and others (2001b), discrepancy between the modeled and ice-core derived age-depth relationship can demonstrate that the incorrect model has been applied. In particular, an inferred negative shear-layer thickness implies that the strain within the ice column is insufficient to accommodate the accumulation history in steady state; this is equivalent to finding annual layers thicker and deeper than the model can achieve. Thus, if the best-fit solution for shear-layer thickness was negative in the Dansgaard-Johnsen model, following Fahnestock and others (2001b), we applied the Nye-melt model, which can move layers down faster and with greater thickness in response to basal melt. If the Dansgaard-Johnsen model achieved a positive shear-layer thickness, the Nye-melt model was not applied, and basal melt was assumed negligible.

Strain-rate modeling along flowlines

We applied these strain-rate models to ice-parcel flowlines that traverse our high-resolution ground-based survey. Flowlines passing through chosen points on the farthest downstream across-flow line in our survey (profile AA' in Fig. 1b) were calculated by following the negative gradient of the present-day velocity field (Joughin and others, 2010). IRH depths along these flowlines were extracted for seven bright, continuous IRHs above the Bølling/Allerød transition at NGRIP (NorthGRIP Members, 2004). These IRHs were dated using the GICC05a timescale (Rasmussen and others, 2006) at NGRIP and traced to our field site in airborne RES data (Gogineni and others, 2001) (Figs 1 and 2).

Fig. 2. (a) Continuity index of northeast Greenland. Ice speed contours are from Joughin and others (2010). Ice flow is from lower left to upper right. Positions of NGRIP (NorthGRIP Members, 2004) and NEGIS (Vallelonga and others, 2014) ice cores are plotted (white dots). Detail boxes (b-d) include the 50ma^-1 velocity contour, airborne radar profiles BB', CC' and DD' (shown in Fig. 3), and bed topography (Bamber and others, 2013).

We fit only IRHs from the Holocene (the last - 1 0 ka) to ensure a temporally stable accumulation rate (Alley and others, 1993; Taylor and others, 1997), and we necessarily assume constant interior ice-sheet geometry during this period.

Results

IRH continuity

The spatial pattern of the radar continuity index provides insight into the history of ice flow in northeast Greenland, which in turn informs hypotheses about future changes. Four key observations are:

1. 1. In upstream regions (Fig. 2b; profile BB' in Fig. 3; profile AA in Fig. 4), continuity is low only in the ice-stream shear margins, with high continuity both outside and inside the ice stream.

2. 2. In downstream regions (Fig. 2d), continuity drops abruptly as ice flows over basal obstacles at the head of a marine overdeepening, and remains low down-flow.

3. 3. In central regions (Fig. 2c; profiles CC' and DD' in Fig. 2 and 3), IRH continuity drops as the ice crosses shear margins or flows over prominent basal topography, but generally recovers down-flow.

4. 4. The highest IRH continuity in all of northeast Greenland is in the upstream regions of NEGIS, including NGRIP and our survey (Fig. 2b).

Fig. 3. IRH continuity index and radargrams of (a) profile BB', (b) profile CC' and (c) profile DD' in Figures 1 and 2. In all panels, ice flow is into the page. Continuous IRHs are found at lower elevation inside the ice stream relative to outside. Gaps in continuity index of profiles BB' and DD' are due to RES profile splicing.

In upstream portions of the ice stream (Fig. 2b), low IRH continuity is restricted to the shear margins. Individual profiles show prominent folding only in these areas (profiles AA' and BB' in Figs 3 and 4), which reduces returned power and thus . The high continuity on both sides of the shear margins and within the ice stream suggests long-term stability, as discussed next.

Fig. 4. Ice-penetrating radar profile along AA' (Fig. 1b). (a) Surface elevation showing the central ice-stream depression and two shear-margin troughs. (b) Radargram showing internal stratigraphy of profile AA'. Seven roughly equally spaced continuous reflectors between the orange dotted lines were used as input to strain-rate models. The orange dotted lines represent reflectors with ages 3.4 ka (upper reflector) and 10.0 ka (lower reflector). Color bar shows absolute value of longitudinal strain rate calculated from InSAR data (Joughin and others, 2010). Note that longitudinal strain-rate amplitude highs are co-located with complex IRH folds at both shear margins. Vertical exaggeration is ~15x. Ice flow is into the page.

We note that if the ice stream did not exist until very recently (-centuries or less), then none of the ice inside the ice stream flowed through the shear margins, and the continuous IRHs inside the ice stream were inherited from a time before the shear margins formed. However, it would be very surprising for an ice stream to suddenly form at a time of no major external forcing ~10ka into the Holocene. Furthermore, Fahnestock and others (2001b) presented evidence of enhanced melting under the ice stream (also see below), and argued that the ice sheet was in long-term balance with this melt, inconsistent with a recent change. Hence, we suggest that the high-continuity IRHs in the ice stream flowed through the low-continuity shear margins.

Our detailed ground-based survey allows us to assess the reason for the low IRH continuity within the shear margins. As discussed by Holschuh and others (2014), the IRHs here show reduced returned radar power because they dip steeply as the limbs of folds (Figs 3 and 4). As first noted by Weertman (1976), basal ice flows over an area of higher basal shear stress (a ‘sticky’ bed), causing local along-flow compression and vertical extension that moves the near-basal IRHs upward; this is reversed as the ice then encounters a ‘slippery’ bed. For this site, we argue that the downward increase in fold amplitude, and the relation of the fold positions to inferred variations in basal lubrication, indicate that the shear margin includes bands of ‘sticky’ and ‘slippery’ bed.

Horizontal flow velocities are fast enough relative to the mean vertical flow from basal melting and ice strain that for steady-state ice flow over a relatively short horizontal distance, such as crossing a shear margin, the path followed by a particle in the ice will approximate the shape of observed IRHs (Whillans, 1976; Whillans and Johnsen, 1983; Christianson and others, 2013). Analogous to an airplane rising and falling on standing lee waves over a mountain range and then flattening away over the plains, the upward and downward deflections of particle paths caused by the bands of varying basal friction produce the observed IRH folds between regions of nearly horizontal IRHs both up-flow and down-flow.

In the downstream region of NEGIS, drops abruptly as ice flows over a prominent basal headscarp into the marine overdeepening (Fig. 2d). Visual inspection of RES profiles shows that IRHs are generally absent farther down-flow. Initially, returned radar power is reduced due to incoherent backscatter off surface crevasses and steeply dipping IRHs. However, in order to enter the overdeepening, ice must flow over -500 m amplitude undulations in basal topography (Fig. 2d). Radar continuity may not recover down-flow because IRHs are permanently tilted past the radar imaging threshold (Holschuh and others, 2014). Alternatively, due to high strain rates as the ice flows over these features, shear zones or recumbent folds may have developed in the ice that permanently disrupt IRH continuity. The absence of surface crevassing, generally smooth bed topography, and lower strain rates in the overdeepening indicate that IRHs could flatten here, increasing . We do not observe this (Fig. 2d). Other factors that can reduce, such as an increase in radar-wave attenuation due to warmer ice at depth, would increase more gradually than the observed drop in, and would vary with depth. Thus, we argue that IRH continuity is permanently disrupted primarily due to high strain as the ice enters the marine overdeepening, although we cannot rule out contributions from changes in the ice surface character caused by melting or crevassing, especially in the ablation zone.

In the central portion of the ice stream (Fig. 2c; profiles CC and DD' in Fig. 3), there is a region where low continuity appears coincident with both high surface strain rates and flow over bedrock obstacles but recovers down-flow. This suggests that low continuity up-flow is attributable to IRH steady-state folds within the ice column, similar to those we observe within our ground survey (Fig. 2b; profile AA' in Fig. 4). As such, our ground survey is likely illustrative of regional conditions.

Inferred surface-accumulation and basal-melt rate

In general, IRHs deflect in response to changes in bed topography (Robin and Millar, 1982; Jacobel and others, 1993), basal conditions (Whillans, 1976; ; Whillans and Johnsen, 1983; Catania and others, 2010; Christianson and others, 2013), or accumulation rate coupled with ice-flow effects (Dansgaard and Johnsen, 1969; Weertman, 1976; Whillans, 1976; Whillans and Johnsen, 1983). At NEGIS, the depth of the 10ka IRH varies by - 1 km along a 40 km transect across-flow in comparison to ~400 and ~20 m variations in bed and surface elevation, respectively (Fig. 4), indicating different conditions along different flowlines. In addition, major fold trains occur in the shear margins independent of topographic undulations. As discussed above, these fold trains generally do not persist as ice enters streaming flow (Fig. 2), indicating that they are probably an in situ steady-state response to variations in basal conditions. As we do not know the details of changing basal ice motion in response to changing basal lubrication or rheology, we take the simplest approach possible to test this hypothesis, i.e. we apply the 1-D strain-rate model to the age-depth relation of IRHs for a representative set of along-flow profiles to infer changes in surface accumulation rate and basal melting: outside the ice stream (profile aa'), in the central trunk (profile bb') and in the southeastern shear margin (profile cc') (Fig. 1). We first present the results, and then consider limitations, identify robust outcomes and discuss future modeling directions. We note in advance that the high-frequency variations along flowlines are indicative but not quantitatively accurate.

Outside NEGIS

Profile aa' runs parallel to the northwest shear margin -10 km outside NEGIS (Fig. 1). The average inferred accumulation and basal-melt rates along this flowline are 0.15 and 0.05 ma^-1, respectively. IRH geometry generally mirrors hydropotential (Fig. 5a). A broad minimum in inferred basal-melt rate of 0.00ma^-1 occurs between 6 and 10 km (coincident with the lowest hydropotential and likely low basal friction) (Fig. 5d). The inferred accumulation rate increases from 0.12 ma^- 1 to 0.16ma^_1 within a broad surface depression (4-12 km in Fig. 5d), consistent with topographic snowdrift trapping.

Fig. 5. Best-fit strain-rate model solutions along flowlines aa'–cc' shown in Figure 1. (a–c) Synthetic radar profiles interpolated along flowlines (a) aa', (b) bb' and (c) cc ' . Color bars are hydropotential, calculated from the same high-resolution ground-based survey. Ice flow is from left to right. (d–f) Inferred accumulation (AR) and basal-melt rate (BMR) from strain-rate modeling. We interpret the patterns shown here, rather than the absolute values, in the text. Note that all ordinate axes have identical scales. Color bars are surface elevation.

Central trunk

Profile bb' runs through the center of NEGIS, and the oldest IRHs likely originate from near the ice stream's onset (Fig. 1). The average inferred accumulation and basal-melt rates are 0.13 and 0.11 ma^-1, respectively. The highest inferred surface accumulation rate (-0.15 ma^-1) occurs in the surface depression near 3 km (Fig. 5e). A shallow core collected at NEGIS (Vallelonga and others, 2014) near profile bb' shows a 0.03 ma^{- 1} decrease in accumulation rate over the last 400 years, consistent with the past position of the core site as it advected through the topographically influenced accumulation-rate field we inferred from Holocene IRHs.

Shear margins

Profile cc' lies within a region of slightly enhanced flow that, if better defined and faster-moving, might be identified as a southern tributary of NEGIS (see Joughin and others, 2001). This profile obliquely crosses the southeastern shear margin (Fig. 1). IRHs are steeply folded (32-37km in Fig. 4 and 2-6 km in Fig. 5c), and the ice compresses and stretches, complicating inferences of accumulation and basal-melt rate (Fig. 5f), but the average inferred accumulation and basal-melt rates are 0.13 and 0.09 ma-¹, respectively. The highest inferred accumulation rate occurs within the shear-margin surface depression (4-12 km in Fig. 5f), again consistent with snowdrift trapping, as suggested by Fahnestock and others (2001b).

Strain-rate model interpretation

For each point on a flowline (10 m posting) in our application, our method (which follows Fahnestock and others (2001b) converts the entire integrated history of basal melting, surface accumulation, and layer thinning from ice flow into an accumulation rate and basal-melt rate, which are thus long-term averages with some poorly defined weighting. Values of these rates at the upstream ends of our profiles thus sample conditions even further upstream. Variations in these values along a profile cannot be interpreted directly to learn local rates quantitatively, but the signs of shifts in accumulation or basal-melt rate should track local anomalies in these quantities affecting the ice above.

Our results are consistent with those of Fahnestock and others (2001b)) in finding relatively similar average accumulation rates for all three flowlines, which agree with accumulation rates for the region estimated using other techniques (Bales and others, 2009; Ettema and others, 2009; Burgess and others, 2010; Vallelonga and others, 2014). In turn, our findings agree with those of Fahnestock and others (2001b) that there have been no strong trends in accumulation rate during the Holocene. The fine-scale structure indicates enhanced accumulation in surface depressions, including those that flank the ice stream, as expected.

The average inferred basal-melt rate outside the ice stream is 0.05 ma-¹, which is significantly lower than the average basal-melt rate inside the ice stream (0.11 ma-¹) and in the line crossing the southeastern margin (0.09ma^-1), but all are quite elevated. The heat generation from sliding plus deformation (basal shear stress times surface velocity) is ~70mWm^-2 outside the ice stream, and -300 mWm-² within the ice stream, with only some of that heat generated at the bed, and with smaller values upstream along the flowlines. Conduction into the ice from the bed at NGRIP is ~55 mWmr² (Dahl-Jensen and others, 2003), and values in our study region are probably not hugely different from that. Thus, although basal shear heating provides an important component, the basal melt rates we infer require locally elevated geothermal flux of ~500 mWm-² outside the ice stream, and ~1000mWmr² within the ice stream and similar for the flowline crossing the southeastern margin. Our elevated values partially reflect conditions upstream closer to the geothermal source of the ice stream, and confirm the high regional geothermal flux suggested by Fahnestock and others (2001b).

Strain-rate model limitations

Uncertainty in both radar and InSAR data must be considered when assessing the strain-rate model results along flowlines. Older dated IRHs may have large age uncertainties due to miscorrelation of IRHs recorded by different radar systems and uncertainty in tying ice-core depth to radar depth (Eisen and others, 2007), although we believe the IRHs used here are correctly identified. The flowlines are also subject to uncertainty in ice speed (-3%) and direction (up to 10°) (Joughin and others, 2010). We conducted sensitivity analyses to determine the effect of these uncertainties on our solutions. The best-fit solution for IRH ages and depths does not vary for a variety of assumed initial conditions, and shows no sign of numerical instabilities or local minima. Solutions are relatively insensitive to changes in basal-melt rate, which agrees with other studies (Fahnestock and others, 2001c). We estimate formal numerical uncertainties of ±0.01 ma^-1 and ±0.05 ma^-1 in accumulation and basal-melt rate, respectively. However, as the local layer approximation (Waddington and others, 2007) generally is invalid in this area, our results should not be interpreted quantitatively, so directly assigning these errors to our results is misleading. Nevertheless, as discussed earlier, there are robust differences between flowlines.

More complex strain-rate models, incorporating either multiple dimensions or additional fitting parameters (Dahl-Jensen and others, 2003; Clarke and others, 2005; Lhomme and others, 2005; Waddington and others, 2007; Koutnik and Waddington, 2012), would provide more-accurate results if they were supplied with proper input data. Frequently these models include upstream corrections and assume relatively homogeneous basal conditions. Additionally, they often incorporate knowledge of ice-column temperature, ice flux, scaling relationships between past ice speed and accumulation rate, and spatio-temporal variation in accumulation rate from other sources (namely ice cores). Unfortunately, these data do not exist for NEGIS, nor would simple scaling relationships hold in a region of streaming flow. One-dimensional strain-rate models including additional parameters (e.g. basal sliding fraction or shear-layer shape functions) might also be applied, but may not greatly improve results, because accumulation and basal-melt rate are the strongest controls on layer height and thickness. Furthermore, none of these methods incorporates basal hydrology, lithology, and side shear in the ice column, which likely have strong effects at NEGIS. Rather than adding additional parameters that we cannot adequately constrain using existing data, we chose to use the simplest (Dansgaard-Johnsen and Nye-melt) models to provide initial guidance for future field and modeling studies.

Discussion

Folds in isochronal IRHs are formed due to ice deflection over bed obstacles (Robin and Millar, 1982; Jacobel and others, 1993) or ice deformation in response to variable basal conditions (Weertman, 1976; Whillans, 1976; Whillans and Johnsen, 1983; Catania and others, 2010; Christianson and others, 2013). Ice advection and/or relaxation can then modify folds so that their location and morphology are inconsistent with the current ice-flow configuration (Jacobel and others, 1996; Ng and Conway, 2004). In this study, we used radar internal stratigraphy to assess current basal conditions and their influence on streaming flow in northeast Greenland. First, we calculated IRH continuity index to locate areas of deformation. Then we inferred generation mechanisms from IRH fold location plus fold morphology, and assessed the strengths and limitations of using strain-rate inversions to further constrain the physical processes underlying IRH deformation.

We used IRH continuity index (Karlsson and others, 2012) of airborne RES data to map IRH disruption in northeast Greenland (Figs 2 and 3). In downstream portions of NEGIS, IRHs are disrupted as ice flows over basal obstacles and enters an overdeepening (Fig. 2d). In contrast, upstream regions (Fig. 2b) have uniformly high continuity except within the shear margins. Consideration of the full suite of data shows that this pattern is very unlikely to have formed either by advection of low-continuity ice to its current position or by very recent ice-stream formation in the current position, and instead that the pattern was produced by nearly steady flow. Firstly it is highly unlikely that an actively advecting feature resulting in low continuity would serendipitously occur everywhere just at the shear margins, including matching the step-out in ice-stream width within our survey (Fig. 1 and 2a and b). There is no evidence of low continuity outside the shear margins, another indication that this signal has not advected from elsewhere.

Secondly, it is also highly unlikely that these low-continuity regions in the shear margins formed due to IRH deformation following initiation of the ice stream within the last - 0 . 1 to 1 ka in a region previously lacking streaming flow, such that the continuous IRHs within the modern ice stream never flowed through a shear margin. In addition to the inherent implausibility of this behavior during the relatively stable Holocene, this hypothesis can also be rejected based on the pattern of IRH deflection. Very rapid melting under the newly formed ice stream would be required to explain the lower elevation of IRHs inside the ice stream, but, in such a scenario, all IRHs would descend almost equally due to very rapid basal melting. In reality, the deeper IRHs within the ice stream have been deflected more (Figs 3 and 4), consistent with steady-state basal melting under the ice stream, likely as a result of basal shear heating. We note that this conclusion is the same as that reached by Fahnestock and others (2001a); our detailed survey confirms their results and extends them to finer spatial scale.

The age-depth relationship of IRHs in our survey indicates basal melt is more rapid within the ice stream than outside and that basal melt has been persistent during the Holocene. Fahnestock and others (2001b) hypothesized that enhanced lubrication due to basal melt is likely the primary cause of streaming flow. This suggests that feedbacks between ice dynamics and basal hydrology influence the extent and vigor of streaming flow along the upstream reaches of NEGIS. Our inferred patterns of accumulation and basal-melt rate support this inference. Importantly, this work highlights that both ice dynamics and geothermal flux result in heightened basal melt beneath the ice stream, and the extent to which these processes are interlinked in this system should be explored further. Although we do not place much reliance on the quantitative details of our results (given the difficulty of observing deep layers in the shear margins and the complexity of variable flow over the bands of strongly contrasting basal shear stress), the inferred pattern of higher vs lower melting in the shear margins is consistent with physical expectations.

Our strain-rate model results have two robust features. Firstly, basal melt is more than twice as high within the ice stream as outside on the north side due to enhanced geothermal flux as well as basal shear heating under the region of streaming flow, and basal melt has been persistent during the Holocene. Secondly, inferred accumulation rates are generally similar to other modern estimates (Fahnestock and others (2001b); Bales and others, 2009; Burgess and others, 2010; Vallelonga and others, 2014), with local positive anomalies occurring in surface lows, likely as a result of snowdrift trapping (Joughin and others, 2001; Fahnestock and others (2001c)).

In general, the pattern of accumulation and basal-melt rate is consistent with current ice-stream geometry and shear-margin location, indicating that major changes in ice-stream location or vigor are unlikely to have occurred during the Holocene upstream of our survey, consistent with the results of Fahnestock and others (2001b). Although the simple strain-rate models applied here map all observed deformation into two parameters (accumulation and basal-melt rate), they allow us to infer shear-margin fold generation via reasonable mechanisms, and our results are consistent with direct examination of the basal interface in RES and seismic data (K. Christianson, unpublished information), strengthening the case that basal hydrological feedbacks contribute to the present geometry of NEGIS. For flowlines transiting a shear margin, at least a 2-D treatment is necessary to fully capture the deflection of IRHs in response to changing basal lubrication (Weertman, 1976; Christianson and others, 2013), as well as changing surface and mass balance (Dahl-Jensen and others, 2003). Therefore, we strongly endorse the use of higher-dimension, more-complex modeling. However, full characterization of the ice stream will likely require a careful consideration of basal hydrology, lithology and a spatially variable basal flow law, processes that are not included in most ice-sheet models. Additionally, the current sparse coverage by comprehensive geophysical data limits application of these more complex models to this region.

Conclusions

As the only Greenlandic ice stream that extends deep into the interior of the ice sheet, NEGIS may be uniquely capable of influencing ice-sheet mass balance via changes in ice dynamics (Joughin and others, 2001; Khan and others, 2014). Our study of radar stratigraphy in NEGIS shows that in upstream portions of the ice stream IRHs deform in response to variations in bed conditions as the ice enters streaming flow. The deformation of ice consistent with spatial variation in basal shear stress, and the lack of basal topographic confinement, suggest control of ice flow by basal-water routing. Although our study indicates that the ice stream has been relatively stable during the Holocene, if our interpretation is correct, sufficient thinning could alter subglacial hydropotential with possible consequences for the ice-stream catchment. Farther downstream, the ice-stream geometry is constrained by basal topography. The ability of ocean-marginal forcing to propagate across the transition zone between these two regimes, where both basal topography and shear-margin location limit and steer fast flow, is unclear. Thus, further field and modeling investigations of upstream and central portions of the ice stream are warranted in order to understand the role of NEGIS in the past, present and future mass balance of the ice sheet.