Airborne RES data

We use airborne radar data collected in May 1999 over northeast Greenland (Reference GogineniGogineni and others, 2001) to calculate IRH continuity (Reference Karlsson, Rippin, Bingham and VaughanKarlsson 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 (Reference Karlsson, Rippin, Bingham and VaughanKarlsson and others, 2012). These radar data are geolocated with dual-frequency GPS data, and provide thickness accuracy of ±10m for >90% of the survey (Reference GogineniGogineni and others, 2001). Typical trace spacing is 150200 m (Reference GogineniGogineni 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 (Reference VallelongaVallelonga 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 (Reference Joughin, Smith, Howat, Scambos and ScambosJoughin 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 (Reference Welch and JacobelWelch and Jacobel, 2003; Reference Welch, Jacobel and ArconeWelch and others, 2009). Each trace is a stack of 2000 radar pulses. Typical raw trace spacing is - 8 m. Processing followed Reference Christianson, Jacobel, Horgan, Anandakrishnan and AlleyChristianson 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 (Reference ChenChen, 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 (Reference Gades, Raymond, Conway and JacobelGades 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 (Reference LimLim, 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.

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 (Reference Matsuoka, Morse and RaymondMatsuoka and others, 2010). An IRH formed in this manner represents a single moment in time, i.e. each IRH is isochronal (Reference WhillansWhillans, 1976; Reference SiegertSiegert, 1999; Reference Kanagaratnam, Gogineni, Gundestrup and GundestrupKanagaratnam 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 (Reference WhillansWhillans, 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 Reference Karlsson, Rippin, Vaughan and CorrKarlsson and others (2009, Reference Karlsson, Rippin, Bingham and Vaughan2012) and Reference Sime, Karlsson, Paden and GogineniSime 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

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 (Reference Karlsson, Rippin, Bingham and VaughanKarlsson 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 (Reference WeertmanWeertman, 1976; Reference WhillansWhillans, 1976; Reference Whillans and JohnsenWhillans and Johnsen, 1983; Reference Waddington, Neumann, Koutnik, Marshall and MorseWaddington and others, 2007; Reference Koutnik and WaddingtonKoutnik 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 (Reference Clarke, Lhomme and MarshallClarke and others, 2005; Reference Lhomme, Clarke and MarshallLhomme and others, 2005; Reference Waddington, Neumann, Koutnik, Marshall and MorseWaddington and others, 2007; Reference Koutnik and WaddingtonKoutnik 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. Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969; Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock and others, 2001b,Reference Fahnestock, Abdalati, Luo and Luoc; Reference Dahl-Jensen, Gundestrup, Gogineni and MillerDahl-Jensen and others, 2003; Reference Waddington, Neumann, Koutnik, Marshall and MorseWaddington and others, 2007; Reference ChristiansonChristianson and others, 2013). Here, following Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock and others (2001b), we apply a 1-D approach to our high-resolution survey of NEGIS, and obtain useful but necessarily limited results.

Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 (Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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) (Reference NyeNye, 1951, Reference Nye1957, Reference Nye1959, Reference Nye1963). This is a ‘plug-flow’ model, in which the horizontal velocity and vertical strain rate are depth-independent. Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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, Reference Dansgaard and JohnsenDansgaard 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 Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 (Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969; Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 (Reference Fahnestock, Abdalati, Luo and LuoFahnestock and others, 2001c):

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 Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 (Reference Joughin, Smith, Howat, Scambos and ScambosJoughin 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 (Reference RasmussenRasmussen and others, 2006) at NGRIP and traced to our field site in airborne RES data (Reference GogineniGogineni and others, 2001) (Figs 1 and 2).

Fig. 2.

(a) Continuity index of northeast Greenland. Ice speed contours are from Reference Joughin, Smith, Howat, Scambos and ScambosJoughin and others (2010). Ice flow is from lower left to upper right. Positions of NGRIP (NorthGRIP Members, 2004) and NEGIS (Reference VallelongaVallelonga 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 (Reference BamberBamber and others, 2013).

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

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 (Reference Joughin, Smith, Howat, Scambos and ScambosJoughin 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, Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 Reference Holschuh, Christianson and ChristiansonHolschuh 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 Reference WeertmanWeertman (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; Reference Whillans and JohnsenWhillans and Johnsen, 1983; Reference ChristiansonChristianson 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 (Reference Holschuh, Christianson and ChristiansonHolschuh 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; Reference Jacobel, Gades, Gottschling, Hodge and WrightJacobel and others, 1993), basal conditions (Reference WhillansWhillans, 1976; ; Reference Whillans and JohnsenWhillans and Johnsen, 1983; Reference Catania, Hulbe and ConwayCatania and others, 2010; Reference ChristiansonChristianson and others, 2013), or accumulation rate coupled with ice-flow effects (Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969; Reference WeertmanWeertman, 1976; Reference WhillansWhillans, 1976; Reference Whillans and JohnsenWhillans 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 (Reference SiegertVallelonga 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 Reference Joughin, Fahnestock, MacAyeal, Bamber and GogineniJoughin 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 Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock and others (2001b).

Strain-rate model interpretation

For each point on a flowline (10 m posting) in our application, our method (which follows Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 (Reference BalesBales and others, 2009; Reference EttemaEttema and others, 2009; Reference BurgessBurgess and others, 2010; Reference VallelongaVallelonga 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² (Reference Dahl-Jensen, Gundestrup, Gogineni and MillerDahl-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 Reference Fahnestock, Abdalati, Joughin, Brozena and BrozenaFahnestock 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 (Reference Eisen, Hamann, Kipfstuhl, Steinhage and SteinhageEisen 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°) (Reference Joughin, Smith, Howat, Scambos and ScambosJoughin 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 (Reference Fahnestock, Abdalati, Luo and LuoFahnestock 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 (Reference Waddington, Neumann, Koutnik, Marshall and MorseWaddington 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 (Reference Dahl-Jensen, Gundestrup, Gogineni and MillerDahl-Jensen and others, 2003; Reference Clarke, Lhomme and MarshallClarke and others, 2005; Reference Lhomme, Clarke and MarshallLhomme and others, 2005; Reference Waddington, Neumann, Koutnik, Marshall and MorseWaddington and others, 2007; Reference Koutnik and WaddingtonKoutnik 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.