2.1. Formulation

Consider a vertical cross section of firn, taken in the ice-flow direction, x. The firn moves in this direction at velocity u(x, t) due to the underlying ice motion and receives surface accumulation at rate a(x, t), where t is time (see Fig. 2). Let z(x, t) denote the depth of an isochronal layer or ‘isochrone’ below the surface. This coordinate ignores variations in surface elevation so that plotting z against x shows the isochrone under a flat horizon, as in Figure 1a. Our analysis of layer geometry considers also the slope of isochrones, the hinges occurring on them, the loci of these hinges (called hinge lines) and the migration velocity associated with hinge lines. Figure 2 introduces these terms, which are explained in more detail as the theory develops.

Fig. 2. Key terms and symbols in our mathematical model. Accumulation rate a (upper plot) and ice-flow velocity u are external forcings generating the isochrone pattern (lower plot). An isochrone has depth profile z(x) and local slope s ( = ∂z/∂x), with s defined to be positive and negative, respectively, where the layer plunges and rises along flow. Hinges are positions where s = 0 and can be troughs or crests, as indicated by the solid dots and open circles on the deepest layer. The locus of hinges from different layers tracing the axes of a fold forms a hinge line. Examples of two types of hinge lines are shown: a trough line (bold curve) and a crest line (dashed curve). As we descend on a hinge line (with depth z), the hinge migrates in the sense that its horizontal position xvaries; the corresponding layer age t increases. dx/dt is the local hinge migration velocity, and (dx = dz)⁻¹ is the local dip of the hinge line. In our theory, Equations (20) and (23) predict these quantities.

As each firn column translates, it receives surface accumulation but stretches or compresses laterally due to gradients of u. If we assume the rate of submergence at any depth to equal the surface burial rate a (expressed as a velocity), then a layer in the column deepens at the rate dz/dt (≡∂z/∂t + u∂z/∂x) = a – z∂u/∂x; thus an isochrone evolves according to the partial differential equation

This kinematic description assumes plane strain (zero flux in the y-direction and infinitely wide flow) and incompressibility for the firn. If necessary, the firn density can be used to convert the units of a to kgm⁻² a⁻¹or to mw.e.a⁻¹.

Throughout this paper, u and a are taken as functions of x only, not of time t. In the resulting steady flow, isochrones of different ages form a stationary pattern, in the sense that each isochrone deforms into the shape of successively older (and deeper) isochrones as it evolves in time. It is therefore possible to construct the entire pattern of isochrones by tracking the evolution of a single isochrone that lies initially at the surface. In this context, the variable t takes the meaning of the age of the isochrones.

Equation (1) is simplistic because it ignores firn densification, which is clearly important in the depth range of interest. While this process does not fundamentally cause isochrones to undulate (we shall see a justification shortly), it affects their depth. Compaction of the firn under its own weight occurs at a rate dependent on rheology, temperature and accumulation rate, resulting in a firn density ρ that increases nonlinearly with depth (e.g. Reference Herron and LangwayHerron and Langway, 1980; Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others, 2010). This causes material at depth to experience a submergence rate less than a. As the isochrones deepen, their shape retains a history of differential submergence caused by spatial variations in a, but densification reduces their mean separation.

In order to account for firn densification in tracking isochrones, several authors (Reference Arcone, Spikes and HamiltonArcone and others, 2005; Reference EisenEisen, 2008; Reference Woodward and KingWoodward and King, 2009) have posited a depth-dependent submergence rate, scaled to the local accumulation rate by the density–depth profile ρ(z), which is assumed invariant across the radargram. Specifically, Sorge’s law (Reference BaderBader, 1954) then gives the submergence rate as ρ ₀ a/ρ(z), where ρ ₀ is the firn surface density. Assuming steady forcings as before, we use these assumptions to formulate the mass-conservation equation

of which Equation (1) is now the special case ρ(z) ≡ ρ ₀.

Although Equation (2) is a more complicated model of isochrone evolution than Equation (1), it can be put in the same form as the latter by the change of variable

which yields

The densification description adopted by us here (and by the aforementioned authors) thus merely introduces a depth correction that is uniform horizontally, implying that densification does not cause isochrones to undulate. However, in Section 4 we discuss caveats concerning this description that lead us to rethink its validity. The issues concern the assumption that ρ is invariant with x, which affects not only Equation (2) but also the practice of compiling radargrams, where it is common to use wave-speed data from a few horizontal positions (obtained via common-midpoint surveys, as in Fig. 1a, or via ρ in firn cores) to convert radar time traces to depth at all horizontal positions.

We proceed with Equation (4) and first show that it can be simplified to a canonical form for any spatially non-uniform flow velocity u(x)>0. Define the ‘transformed depth’ variable Z(x, t) = u(x)f(x, t)/u ₀, where u ₀ = u(x)=0) is a constant reference velocity, here taken at the up-flow end of our domain. Then Equation (4) becomes

and changing the variable x to X (‘transformed distance’) with

yields the canonical kinematic wave equation

Some radar surveys have been made where u increases along flow, for example in the onset zones of ice streams, where we may suppose a linear velocity model

in this case Z = (1 + kx)f and Equation (6) leads to

and

In the rest of this paper, we analyse the problem in the transformed coordinates and use lower-case symbols for them (write x for X, z for Z and a for A) for convenience, with the understanding that reduction to the canonical form has been made and that we are referring to transformed layer geometries unless we indicate otherwise. Then Equation (7) becomes

It is this simple wave equation that ultimately governs the spatial structure of the layer folds. Notice if a(x) and u ₀ are multiplied by the same factor, t can always be rescaled to recover the original equation; this means that different forcings with the same ratio a(x)/u ₀ generate the same isochrone pattern. Such non-uniqueness does not alter our findings on the pattern-formation mechanism below, but matters in inversions for the forcings from the layers (we shall see an example of such inversion in Section 3.3).

2.2. Propagation of information

Wave Equation (10) can be solved exactly by the method of characteristics (Reference Carrier and PearsonCarrier and Pearson, 1988) to give

along the characteristic line or ‘characteristic’

Here d/dt is the total derivative, and Equations (11) and (12) track the positions of material points (x, z) parametrically. Material at position x on an isochrone aged t had the initial position x ₀ at the surface when t = 0. As it moves down-glacier at velocity u ₀, it deepens at a rate determined by the local accumulation rate (the use of transformed coordinates already accounts for the effects of densification and longitudinal flow extension/compression). Its z-value thus travels and varies along the characteristic. Figure 3 shows this idea on the x–t plane (distance–age domain), where isochrones are lines of constant t. Note that in this paper, t refers always to time or age, and not to the two-way travel time in the geophysical measurement of radar traces.

Fig. 3. The x–t plane and the characteristic line (or ‘characteristic’) along which information travels. Depth z and slope s of isochrones, both functions of distance x and age t, may be regarded as 3-D surfaces over the plane.

Now, integrating Equation (11) with x taken from Equation (12) and with the initial condition z = 0 at t = 0, yields the explicit solution

(we have used the substitution ζ = x ₀ + u ₀ τ), which calculates each material particle’s depth from the total overlying accumulation that it accrued as it travelled the horizontal distance u ₀ t to its current position. While z evaluated with Equation (13) at a given time describes an isochrone’s shape, z(x, t) can be portrayed generally as a three-dimensional (3-D) surface over the x–t plane. Isochrones on a radargram are projections of fixed-time samples of this surface onto the x–z plane.

A key insight in this paper is that the slope of isochrones also obeys a kinematic wave equation. Let s = ∂z/∂x be the isochrone slope (Fig. 2). Since z is referenced to the surface and z and x are transformed depth and distance, s is not the true slope against the horizontal, but related to it. Differentiating Equation (10) with respect to x gives

where b is the spatial gradient of the surface accumulation rate. Both this wave equation and Equation (10) apply at positions x where a is differentiable (we assume this to be the case throughout our analysis), including where a changes smoothly from accumulation (a > 0) to ablation (a < 0). However, note that a < 0 causes unconformity in the stratigraphy by removing parts of isochrones and making them discontinuous.

Solving Equation (14) with the method of characteristics yields

which describes how isochrone slope varies along the same characteristics (as in Equation (12)) in response to the overlying accumulation rate gradient. The corresponding integral solution satisfying s = 0 at t = 0 (the surface isochrone is ‘flat’) is

Hence we may also regard s(x, t) as a continuous surface over the x–t plane – we call this the ‘slope surface’ (Fig. 4) – with fixed-time samples of it representing the slope variation of individual isochrones. And since slope information propagates on the x–t plane, we may consider how s varies across the radargram (x–z plane). These ideas help us understand the spatial structure of layer folding and are explored further in Section 2.4. The solutions above may in principle be found in the untransformed coordinates, but with much more algebra.

Fig. 4. The function s(x, t) as a 3-D ‘slope surface’. The axes are distance x, isochrone age t and isochrone slope s. Intersection of the slope surface with the horizontal plane defines a system of hinge lines (dashed).

2.3. Retrieving surface accumulation rate

Before turning to the folds, we use Equation (13) to establish two ways of extracting the accumulation rate pattern from isochrones that we build upon later. As with Equation (13), the following results rest on the steady-flow assumption and we are working in the transformed coordinates.

Method I is called ‘shift-differencing’. Consider an isochrone aged t and a slightly deeper, older isochrone aged t + δt. Shift these isochrones down-glacier and up-glacier, respectively, by the distance u ₀ δt/2, then subtract their depths. Equation (13) shows that

Therefore, a at position x can be found by subtracting the depth of the younger isochrone at x − u ₀ δt/2 from the depth of the older isochrone at x + u ₀ δt/2, and dividing the result by δtt. This method of estimating the accumulation rate pattern a(x) avoids recasting Equation (4) as a formal inverse problem (e.g. Reference EisenEisen, 2008) and has an error of order δt ² caused by truncation of higher-order terms in Equation (17).

Method II involves subtracting the depths of the two isochrones without shifting them and also yields information about a. In this case

The new differencing retrieves the accumulation pattern at a distance u ₀ t up-glacier, with an error of order δt. If t = 0, the upper isochrone is simply the firn surface, and Equation (18) shows that a(x) can be inferred from a shallow isochrone (e.g. Reference Gray, King, Conway, Smith and NgGray and others, 2005; Reference Woodward and KingWoodward and King, 2009).

Both methods here rely on an isochrone separation small enough to limit truncation errors, but not so small that it becomes dominated by measurement uncertainty in the layer depths. Both methods require the age of isochrones (or at least their age difference) to be known, for example from a dated firn core at the site. As this is not always available, in Section 3.2 we develop a method that could be used without dated layers. Finally, while methods I and II will not work on isochrones that have been made discontinuous by ablation (or wherever z is undefined), they can be used on continuous isochrones that bracket ablation unconformities, as long as the differentiability condition following Equation (14) is satisfied. In this case, shift-differencing the isochrones will make them intersect and retrieve negative as well as positive values of a.

2.4. Structural architecture of layer undulations

Recurrent properties of folded isochrones in the firn have been recognized in numerous flow-aligned radargrams from West Antarctica. It is noticed that the hinge axes of layer folds dip at angles that reflect a migration velocity less than the ice-flow velocity, and that these axes can dip up- as well as down-glacier (Reference King, Morse, Alley, Anandakrishnan, Blankenship and SmithKing and others, 2002; Reference Arcone, Spikes and HamiltonArcone and others, 2005; Reference Gray, King, Conway, Smith and NgGray and others, 2005). By choosing suitable forms of u(x) and a(x), these authors have used numerical models of isochrone evolution to recreate layer patterns and simulate the observed dips. Reference Arcone, Spikes and HamiltonArcone and others (2005), in particular, posed sinuisoidal and Gaussian functions for a(x) to study the controls on the dip. Some radargrams also exhibit folds whose amplitudes increase with depth (Reference Vaughan, Corr, Doake and WaddingtonVaughan and others, 1999, Reference Vaughan, Anderson, King, Mann, Mobbs and Ladkin2004; Reference King, Morse, Alley, Anandakrishnan, Blankenship and SmithKing and others, 2002; Reference Arcone, Spikes and HamiltonArcone and others, 2005; Reference Gray, King, Conway, Smith and NgGray and others, 2005), a property that has been observed on deep radar layers in ice streams (Reference Ng and ConwayNg and Conway, 2004). Examples of the structures described here can be found in Figure 1a.

Given that we know the folds form from the cumulative effect of flow and submergence, what mechanisms govern their structure? And can we understand them without simulation? We advance a theory of folding that invokes the slope surface.

3.1. Synthetic radargram

We first use the forward model to create a radargram for which everything is known. In this example, we assume the accumulation rate pattern in Figure 6a, and we ignore firn densification and assume a constant ice-flow velocity u = u ₀ = 40m a⁻¹, which allows us to work directly in the transformed distance and depth coordinates. With these steady forcings, Equation (10) is solved numerically by finite difference (with upwinding and explicit time-stepping) to evolve an isochrone from the surface, thus generating a sequence of isochrones with an age increment of 0.025 years. Edge effects were avoided by assuming the 10 km horizontal domain to be periodic; we halted the simulation at t = 150 years before this caused the isochrone pattern to repeat at depth. The isochrones were subsampled at 2.5 year intervals to form the radargram, which is shown in Figure 6b.

Fig. 6. Analysis of a synthetic radargram. (a) Prescribed accumulation rate pattern. (b) Synthetic radargram made from model-simulated isochrones with a constant age increment of 2.5 years. (c) Isochrone slope map on the x–t plane. (d) Isochrone slope map on the x–z plane. In (b–d), solid and dashed curves trace the hinges of fold troughs and fold crests, respectively.

As expected, down-pointing folds develop in the shallow subsurface under major peaks in the surface accumulation (see Fig. 6a and b). Deeper down, however, the undulating pattern becomes less predictable and shows a level of complexity like that shown in Figure 5.

To study the fold pattern, we traced hinge lines across the radargram by linking the stationary points of isochrones found by differentiating their depth (which also yields their slope). Trough and crest lines, distinguished in Figure 6b by different line types, display a topology confirming our predictions. Despite their irregular shape they form loops, dividing the radargram into jigsaw areas of opposite isochrone slopes; the loops are closed unless they meet the surface. Hinge transitions of the kinds described in Section 2.4.1 are abundant. The upper/lower end of each loop locates where a pair of crest and trough lines emerges/ vanishes; such pairing causes hinge lines to switch between crests and troughs in the x-direction. At hinge transitions, the trough and crest lines join smoothly with zero gradient and do not form cusps. We know from our theory that this is because the underlying slope surface is smooth.

Different portions of the hinge lines in Figure 6b show a variety of dips, some pointing up-flow but more often down-flow (cf. Arcone and others, Reference Arcone, Spikes and Hamilton2005). Two dips are common. One of them is nearly vertical. The other dip has a down-glacier slope of ∼0.01, as shown by several loops that swing towards the lower right. This dip seems to track the mean direction of advection coupled with burial, since material moves down-glacier by 6 km in the time it takes to reach the depth (∼60 m) of the oldest isochrone (t = 150 years). Both dips here are associated with the zones of steep isochrone slope discussed below.

We examine the folds more generally by considering isochrone slope s, which varies with x and z on the radargram, but, as explained in Section 2, may be conceived to vary with xand t. Figure 6c and d plot the functions s(x, t) and s(x, z) in colour. These ‘slope maps’ are easy to compile for this radargram because we know the age and depth of its isochrones. Features on s(x, z) are distorted versions of those on s(x, t), but can be superposed directly on the stratigraphy in Figure 6b.

We focus on s(x, t) (Fig. 6c) because it shows the essence of slope information travel. Intriguingly, we see two distinct directions of constant-slope migration, shown by the vertical zones of red and diagonal zones of blue, both descending from major accumulation anomalies (cf. Fig. 6a). The red is aligned with peaks in a(x); the blue emanates from the lee of these peaks and inclines at ∼40ma−¹|(u ₀). As these colours respectively signify strongly positive and strongly negative slopes, their pattern means that steep slopes on the descending flanks of layer folds migrate at velocity v≈0, whereas steep slopes on the ascending flanks of layer folds migrate at v≈u ₀.

This pattern is explained by slope evolution along the characteristics (defined in Section 2.2), which descend across the x–t plane at velocity u ₀. Equation (15) shows that on this plane, s depends on its initial value on its characteristic (at t = 0) and on its alteration by b(x) along the characteristic. The observed pattern stems from the contest between these inheritance and modification factors, as may be understood by considering two end-member scenarios. First, slopes that are shallow initially (of either sign) can be changed easily by a strong accumulation anomaly on their course, to adopt its sign at its x-position. Thus a major peak or valley in a(x) causes isochrones beneath it to plunge or rise steeply, respectively, forming a vertical zone of high slopes (e.g. red in Fig. 6b). Second, an initially steep slope (of either sign) inherited from a strong accumulation anomaly will stay steep unless strong opposite anomalies lie on its characteristic to negate it. In this case, a set of steep fold limbs drifts down-flow with little slope change to form a diagonal zone; an example in Figure 6b is the rising limbs between x = 2 km, z = 10m and x = 6 km, z = 45m.

These scenarios together explain why each accumulation anomaly creates a pair of vertical and diagonal zones. Moreover, since our wave equation is linear, where these zones intersect they disrupt each other by partial cancellation of layer slopes (Fig. 6c). The current analysis identifies such spatial arrangement, which we summarize in Figure 7, as the hallmark of layer patterns formed under steady forcings. We learn that although the layer folds in different steady-flow radargrams may look distinct, they obey the same general rules of patterning, with the main difference being the location and size of accumulation anomalies. We also learn the usefulness of the slope map s(x, t) as a tool for elucidating the structure of these folds.

Fig. 7. Schematic of the criss-cross arrangement of distinct vertical and diagonal zones on the x–t plane of a steady-flow radargram, as caused by anomalies in accumulation rate a(x) (cf. Figs 6c and 10g). Red and blue zones signify where isochrones plunge and rise, respectively. Constructive superposition of layer slopes occurs where same-coloured zones intersect, reinforcing each other’s slope contributions. Destructive superposition of layer slopes occurs where different-coloured zones intersect and offset each other’s contributions.

3.2. Extracting a (x) and layer ages through optimization

Compiling the map s(x, t) for any radargram requires the age of its layers to be known. Here, making further use of the synthetic radargram, we devise a method (based on method I in Section 2.3) capable of retrieving the accumulation rate a(x) and layer ages simultaneously from steady-flow radargrams.

Given the depths of two adjacent isochrones, z ₁(x) (shallower) and z ₂(x) (deeper), we can use Equation (17) to estimate a(x) by means of three steps:

These steps involve ‘sh ift-differenci ng’ the layers (by the shift distance Δ), inferring their age difference δt from Δ, and finding the accumulation rate by dividing the depth difference δz by δt. As before, x and z here are transformed coordinates.

In scenarios where either δt or u ₀ is unknown, the correct shift distance Δ to use in step (1) is unknown. However, given multiple pairs of isochrones, we can choose a combination of Δ’s for these pairs to optimize the agreement between the a(x)’s reconstructed from them, thus also determining the C’s. In this method, a key assumption is that the isochrone pattern is stationary, generated by steady forcings.

In an experiment, we tested this method by using all 61 layers of the synthetic radargram (Fig. 6b) and using the knowledge that their age difference is uniform so that the same Δ can be used in step (1) to shift-difference all consecutive layer pairs. For each choice of Δ, this step yielded 60 difference profiles (the δz’s). We then calculated the variance of these δz’s at each x-position and used the spatial mean of this variance as a single measure of mismatch between the difference profiles. Figure 8a shows that the mismatch is minimized (agreement between the difference profiles maximized) when Δ = 100 m. With our additional knowledge of u ₀ = 40ma−¹ in the synthetic problem, this optimal shift distance correctly recovers δt = 2.5years in step (2). The corresponding estimates of a(x), found from step (3) and plotted in grey in Figure 8b, match the actual accumulation pattern in Figure 6a closely, despite a small scatter due to truncation error.

Fig. 8. Retrieval of accumulation rate from the synthetic radargram layers in Figure 6b. Results of two methods are presented. (a) Mismatch between δz profiles derived from shift-differencing all consecutive layers by the same shift distance Δ, plotted against the value of Δ. The point of minimum mismatch identifies the optimal shift. (b) Retrieved accumulation rate patterns: grey curves show the results of the method in (a), the black curve the result of the method in (c) and (d). (c) Layer pairs chosen to illustrate the method in Figure 9. (d) Mismatch (log-10 scale) across the parameter space of shift distances Δ₁ and Δ₂. The combination of Δ₁ and Δ₂ with minimum mismatch yields the black curve in (b).

In practice, many radargrams are not furnished with local firn cores to constrain the age of their stratigraphy (e.g. via annual-layer counting); when ‘picking’ layers from them there is no way to ensure uniform age difference between the picks. But we can circumvent this requirement because steps (2) and (3) in Equation (24) can be combined to give

Consequently, when shift-differencing multiple layer pairs to find a(x), we can minimize the mismatch between profiles of δz/Δ rather than of δz. The optimization now involves choosing as many Δ’s as there are layer pairs, not just one Δ.

We tested this method on the synthetic radargram taking only the two pairs of isochrones shown in Figure 8c. Figure 9 illustrates the shift-differencing procedure on the x–t plane. Two shift distances Δ₁ and Δ₂ were chosen to difference the pairs to yield two estimated profiles of δz/Δ, whose variance was summarized as a mismatch in the same way as before. The best combination, Δ₁ = 395m and Δ₂ = 201m, was found by searching for minimum mismatch over the parameter space (Fig. 8d). The corresponding reconstructed a(x) is the black curve in Figure 8b. It is smoother than the original accumulation forcing, indicating that truncation error associated with the separation between the input layers causes low-pass filtering. Nevertheless, with only four isochrones as input, this reconstruction recovers the peaks and troughs of the forcing remarkably well.

Fig. 9. The procedure (on the x–t plane) of shift-differencing two pairs of isochrones whose age differences are unknown and unequal

Besides its ability to retrieve accumulation rates, this method has three useful properties that are worth noting:

1. 1. The optimal shift distance Δ for each pair of layers divided by the flow velocity u ₀ gives their age difference (step (2) in Equation (24)).

2. 2. Even if u ₀ is unknown, the method can reconstruct the normalized accumulation rate pattern, a(x)/u ₀ (see Equation (25)).

3. 3. Since the method assumes steady forcings, layer patterns formed by a and u that varied temporally as well as spatially (so that a/u ₀ is not invariant) will cause fundamental mismatch between the δz/Δ profiles from different layer pairs that cannot be reconciled by any shift combination. In this sense, the method is able to discern the existence of unsteadiness, although not quantify its history.

3.3. Radargram from Bindschadler Ice Stream

Our real case study uses the radargram in Figure 1, obtained by ground-based radar on a traverse at the onset zone of Bindschadler Ice Stream in the 2001/02 austral summer. Along the 52.7 km radar line, surface velocity increases from 59 to 111m a⁻¹. Radar traces were recorded by a 100MHz pulseEKKO radar at 0.8 ns time interval and every 5m along flow. Conversion of two-way travel time to depth used the wave speed in firn measured by a common-midpoint survey at km 33, which shows that this speed decays smoothly from ≈0.23mns^–1 at the surface to a constant value of ≈0.168mns^–1 below a depth of 50m, consistent with an expected increase of firn density with depth. Migration was found to cause negligible change to the depth and slope of isochrones and so was ignored. For details of radar processing, see Woodward and King (Reference Woodward and King2009).

The modern ice-flow dynamics of this region have been characterized by field and remote-sensing observations and numerical modelling (e.g. Joughin and Tulaczyk, Reference Joughin and Tulaczyk2002; Price and others, Reference Price, Bindschadler, Hulbe and Blankenship2002), but little is known about their history over centuries or longer. For Bindschadler Ice Stream, interpretation of its deep radar layers (Siegert and others, Reference Siegert, Payne and Joughin2003) and of features downstream of it on the Ross Ice Shelf (Hulbe and Fahnestock, Reference Hulbe and Fahnestock2007) does not suggest that its flow has undergone major fluctuations. However, deep radar-layer folds upstream of the onset zone have been interpreted for changes in flow direction there ∼1.5 ka ago (Siegert and others, Reference Siegert2004).

We picked the depth profiles of 17 clearly visible layers in Figure 1a and smoothed them by low-pass Fourier filtering to suppress noise for the later slope calculations. In view of the flow acceleration along the radar line, we first rendered the profiles in the transformed coordinates xand z(Section 2.1). We used the linear velocity model in Equation (8) with u ₀ = 59ma⁻¹and k = 0.0167km⁻¹; x is then related to true distance via Equation (9a). For the firn-densification correction in Equation (3), we assumed the density profile ρ(z) = ρ_i − (ρ_i − ρ ₀)e^−cz with ρ _i = 917 kgm^–3, ρ ₀ = 400 kg m^–3 and c = 1/35m⁻¹ (Arcone and others, Reference Arcone, Spikes and Hamilton2005) derived from International Trans-Antarctic Scientific Expedition (ITASE) Core 99-1, the nearest firn core (≈25 km upstream) from our line (Fig. 1b). Figure 10b shows the distance conversion and Figure 10a plots the layers after transformation, which squashes them laterally and depresses them towards the right-hand side of the domain.

Fig. 10. Model analysis of the BAS Line 11 radargram in Figure 1a. (a) 17 picked isochrones (including the surface) in the transformed coordinates x and z. (b) Conversion between transformed distance x and true horizontal distance based on Equation (9a). (c) Initial estimates of the normalized accumulation rate from each of the 16 layer pairs (see step 1 of the method in Section 3.3). (e) Final estimates of the normalized accumulation rate (step 2 of the method) and (d) the corresponding optimized shift distances for the layer pairs. In (c) and (e), grey curves show estimates found from individual layer pairs, black curve shows the mean of these estimates, and dashed curve shows their standard deviation. Dotted curve in (e) shows a(x)/u ₀ found by shift-differencing the shallowest subsurface layer against the surface layer according to the method in Equation (17); it overlaps the grey curve corresponding to this pair of layers. (f) Inferred age of each layer vs the mean of its true depth. The age–depth profile from ITASE Core 99-1 is included for comparison. (g) Map of isochrone slopes (in transformed coordinates) on the x–t plane compiled using the layer ages in (f). (h) Map of isochrone slopes on the untransformed radargram domain. Also shown are hinge loci of fold troughs (solid curves) and fold crests (dashed curves) traced from Figure 1a.

Of interest are the hinge-slope structures of the layer pattern and their links with the accumulation rate pattern. Unlike in the synthetic study, however, here a(x) is not known, nor are layer ages needed for compiling the slope map s(x, t). We estimated these by the inversion method developed in Section 3.2, assuming the picked layers to be isochronal and their pattern stationary.

This problem is more complicated than the example in Figure 8c because we need to choose 16 shift distances (one Δ for each pair of consecutive layers) to maximize the agreement between the δz/Δ’s found from shift-differencing the respective layer pairs. We determined the Δ’s in two steps: by finding rough guesses and refining them. In the first step, the method in Figure 9 was applied to two pairs of consecutive layers at a time, but to all 15 consecutive pairs of layer pairs. For the top- and bottommost layer pairs this yielded the Δ guess directly; for the other layer pairs this yielded two estimates of Δ, which we averaged to form the guess. In Figure 10c, the grey curves plot the normalized accumulation profiles a(x)/u ₀ evaluated by Equation (25) with these Δ guesses. Although they look shifted vertically from each other, their peaks and valleys show correlation.

In the second step, we adjusted the Δ guesses to optimize the agreement between the δz/Δ profiles found from shiftdifferencing all 16 layer pairs together; accordingly we redefined the measure of mismatch based on all these profiles. The optimization used a Monte Carlo scheme with many trials (Press and others, Reference Press, Teukolsky, Vetterling and Flannery1992) and the Δ guesses as starting point. In each trial, we chose a pair of consecutive layers at random and improved its Δ by minimizing the mismatch, keeping the other Δ’s fixed. This procedure brought different estimates of a(x)/u ₀ into focus without the need for exhaustive search on the 16-dimensional space of the Δ’s. The mismatch decayed rapidly at first and stabilized at :≈7 × 10⁻⁴ m²after about 100 trials. Figure 10d shows the optimal shift distances. As expected, these distances are larger for layer pairs that lie further apart. Since the edge of this radargram caused the δz/Δ profiles to be shorter than the domain, in both steps above we evaluated their variance (and a(x)/u ₀) only where they exist and overlap.

In Figure 10e, the grey curves show the normalized accumulation rate profiles a(x)/u ₀ calculated with the optimal shifts. Their broad agreement means that different layer pairs tell similar stories about past accumulation; their scatter is worse than in the synthetic study but remains small compared with their mean. The scatter could arise from many factors. One factor is numerical errors from the layer picking, conversion of radar time to depth, and the truncation approximation in Equation (17). A second factor is processes unaccounted by our model, such as temporal changes in u and a, lateral flow convergence or divergence, and spatial variability in densification. A third factor is geometries that render the radargram imperfect, including out-of-plane reflections and offset of the radar line from the ice-flow direction.

In some studies, the accumulation rate distribution was estimated from the depth and age of a shallow layer (e.g. Gray and others, Reference Gray, King, Conway, Smith and Ng2005; Woodward and King, Reference Woodward and King2009). To see how well this method performs, in Figure 10e we plot the normalized accumulation rate found by shift-differencing our uppermost subsurface layer and the surface layer (see the dotted curve). The method in Equation (17) was used, with u ₀ δt equated to the optimal shift distance of the layer pair (this effectively constrains the age of the lower layer). As our two-step optimization is also based on Equation (17) and the same shift, the two methods yield identical results of a(x)/u ₀ (the dotted curve overlaps the grey curve for the layer pair). The shallow-layer result here resembles our mean optimized estimate for a(x)/u ₀ (black curve in Fig. 10e) due to the limited scatter discussed above. However, when used on radargrams whose forcings have varied substantially, the shallow-layer method will not give valid accumulation estimates for times preceding the age of the lower layer.

We next converted all the optimal shift distances Δ_i (i = 1–16) into age differences and summed these to estimate the age of the picked layers. We used the equation

where m is the number of the layer and t_m its age. Figure 10f plots the estimated age of each layer against its original (untransformed) mean depth on the radargram. This age–depth relationship gives older ages than found by layer counting in ITASE Core 99-1 (Steig and others, Reference Steig2005), implying lower accumulation rates in our study area than at the core site tens of kilometres away. Indeed, our mean untransformed value of a(x) is 0.273ma−¹ or 0.109mw.e. a−¹ (Fig. 11a, top plot), which is 20% less than the mean accumulation rate 0.136mw.e. a−¹ derived for Core 99-1 over the period 1713–2000 (Dixon and others, Reference Dixon, Mayewski, Kaspari, Sneed and Handley2004). Our age–depth relationship also shows an upward bend caused by firn densification, but we emphasize that it represents a spatial mean only. The local age–depth relationship must vary across the radargram because the radar layers undulate.

Fig. 11. Comparison of radar layers simulated by the forward model (curves on the lower plots) with those in the recorded radargram in Figure 1a (grey background in the lower plots). Upper plots show the accumulation rate forcings used in two different simulations: (a) the mean reconstructed a(x) from Figure 10e; (b) the a(x)’s reconstructed from individual layer pairs from Figure 10e. The accumulation rate scales are shown in both the units of velocity and water-equivalent units.

The slope maps s(x, t) and s(x, z) of this radargram (Fig. 10g and h) share strong similarities with those of the synthetic case study. The former map (Fig. 10g) was compiled using the layer ages estimated above and uses transformed distance as its horizontal axis. The latter map (Fig. 10h) uses true depth and true distance as axes. Both maps were made by interpolating the slopes of the 17 layers and without picking additional layers. As in the synthetic study, s(x, t) here shows the criss-cross arrangement of vertical red and diagonal blue zones, descending respectively from the peaks and lee sides of accumulation anomalies. The blue zones in Figure 10g are straightened from those in Figure 10h and incline at a velocity near u ₀. Compared with Figure 6c and d, both maps show more fine-scale features, some of which presumably stem from the variability that causes the minor differences between the retrieved profiles of a(x)/u ₀ in Figure 10e.

These results imply that the layer pattern on this radargram can be explained by time-invariant forcings through our proposed folding mechanism. Agreement between a(x)/u ₀ found from the different layer pairs reflects limited historical variations in this ratio and is consistent with the steady-flow assumption behind our analysis. If we use Equation (10) to forward-model the pattern with u ₀ and the mean reconstructed a(x) as inputs, and show the results in true coordinates, then isochrones simulated at the age of the picked layers are found to mimic the latter well in terms of the general location and structure of key folds (Fig. 11a). Improvement in matching the depth of isochrones is found if we use the a(x)’s reconstructed for the times spanned by successive layer pairs (Fig. 11b).

Inferences from these findings about the ice stream’s history are non-unique because our model precludes temporal variations in a and u that might have occurred. Specifically, three different interpretations are possible:

1. 1. Ice flow and accumulation conditions in the onset zone have been steady over a time going back to the age of the deepest picked isochrone (386 years in Fig. 10f). This interpretation implies centurial-scale ice-stream stability.

2. 2. Both a and u changed over time, but co-varied in a way to maintain approximately the same a/u ₀ ratio, producing an apparently stationary isochrone pattern.

3. 3. The accumulation rate pattern has been stable but has been migrating up- or down-glacier, producing an apparently stationary isochrone pattern. This is possible if we suppose an accumulation forcing of the form a(x − wt), where w is the migration speed. Putting this forcing in our model of isochrone evolution in Equation (10), we see that the latter can be rewritten as

   (27)

   

   with x ^* = wt being a moving coordinate. As this equation has the same form as the original (albeit with advection velocity u ₀ − w), layer patterns generated by it will appear to have been caused by steady forcings when observed at a given time, but are actually translating at speed wwith respect to the stationary frame of reference. Migrating surface accumulation patterns not only could, but do, exist, as have been inferred for various dune-like features in East Antarctica (see review by Eisen and others, Reference Eisen2008).

In both the second and third interpretations, the layer ages derived from Equation (26) and used to form the age–depth relationship and s(x, t) (Fig. 10f and g) would be incorrect.

Here we favour the first interpretation for the following reasons. The second interpretation requires the accumulation rate to be coupled to the ice flow, most likely via interaction between atmospheric processes and surface topography. While such coupling is possible, for example via snow redistribution by katabatic wind (e.g. Black and Budd, Reference Black and Budd1964; King and others, Reference King, Anderson, Vaughan, Mann, Mobbs and Vosper2004), we cannot envisage how it causes concerted changes where variations in the ice-flow velocity u modulate the mean rate of accumulation and its anomalies by the same factor. The first interpretation needs a mechanism to fix a(x) in space, and we find evidence for this and against the third interpretation. Figure 12 shows convincing correlation between the retrieved accumulation pattern a(x) and the (along-flow) surface slope of the ice stream and also between the surface slope and the ice-stream bed topography. We interpret these correlations to mean that bed topographic undulations are inducing surface topography via ice-flow dynamics (this process is well recognized and studied in the literature; e.g. Gudmundsson, Reference Gudmundsson2003); in turn, the surface topography creates anomalies in a, enhancing it on up-glacier-facing slopes and reducing it on down-glacierfacing slopes. Information given by Arcone and others (Reference Arcone, Spikes and Hamilton2005) for their figure 5, which shows a similar correlation between slope and accumulation on the ITASE traverse near the Core 99-1 site (Fig. 1b), suggests that our up-facing and down-facing slopes are windward and leeward, respectively.

Fig. 12. Correlation between the retrieved mean accumulation-rate distribution a(x) and the ice stream’s cross-sectional geometry along the radar line. (a) Surface and bed-elevation profiles of the ice stream. (b) Surface slope in the ice-flow direction (solid curve) and the bed topography (dashed curve; plotted in inverted scale to highlight its relationship with surface slope). (c) The retrieved mean accumulation rate distribution a(x) derived from Figure 10e. Comparison of (b) and (c) shows high accumulation rates on up-glacier-facing slopes and low accumulation rates on down-glacier-facing slopes, if we define these slope types by using the mean slope as reference.

These mechanisms imply a basal origin and stationarity for the anomalies on a(x) and flow stability in the onset zone. While this interpretation is specific to this case study and may not hold for other areas of the ice sheet, we note that the effect of basal topography on accumulation rates has been recognized independently by other radar studies (e.g. Rotschky and others, Reference Rotschky, Eisen, Wilhelms, Nixdorf and Oerter2004). Furthermore, while the flow-stability interpretation may not seem sensational, it adds to our knowledge of Bindschadler Ice Stream. One way to test it is to compare our layer-age estimates in Figure 10f with the layer-counted ages in a firn core drilled in our study area, or with the layer-counted ages in ITASE Core 99-1 after further radar profiling to link the isochrones between our study area and the core site.

3.4. Hinge and slope migration

We use the layer patterns in the synthetic and the real radargrams to evaluate the hinge and slope migration rates predicted by theory, given by Equations (23) and (20), respectively:

Working in transformed coordinates, we tested each equation by a separate procedure. For each radargram we first compiled, at various x-positions, values of a(x) and b(x) as well as of s and ∂s/∂x on its picked layers (the 61 layers of the synthetic radargram and the 17 picked layers in the real radargram). The quantity u ₀ and the age difference between successive layers are known.

Equation (23) predicts the local dip of hinge lines, which have been traced on both radargrams (Figs 6b and 10h). For each hinge-line segment between successive layers, we used the equation with the compiled data to calculate the values of dx/dz at its two end points, then compared their mean with dx/dz measured directly for the segment from the transformed radargram. We made this comparison for all traceable hinge-line segments, and defined their dip angles as tan−¹ (dx/dz ÷ 100) (so that 0° points downward) rather than measuring both up- and down-glacier dips against the horizon, which causes a sign discontinuity across the vertical. The factor of 100 here puts the angle more in line with the vertically exaggerated perspective of the radargrams shown in this paper.

Figure 13 plots the results, distinguishing those for crest and trough lines by different symbols. The strong correlation in Figure 13a for the synthetic radargram reveals the predictive power of Equation (23). By contrast, Figure 13b shows a weak, almost non-existent correlation for the real radargram. We attribute this to unsteadiness in its forcings (Section 3.3) and the large separations between picked layers, which introduce significant errors in the dx/dz estimates in the comparison. We experimented with resampling the synthetic radargram at greater age interval and found that the scatter in Figure 13a is also caused by such errors, because it increased when we used fewer, more widely spaced isochrones.

Fig. 13. Testing hinge-line Equation (23) with data from: (a) the synthetic radargram in Figure 6b and (b) picked layers on the radargram used in our real case study. Vertical axes plot hinge-line angles measured from the transformed radargrams. Horizontal axes plot hinge-line angles predicted by Equation (23). See Section 3.4 for our definition of these angles.

Equation (20) predicts the migration velocity v of a constant isochrone slope as we cross from one layer to a neighbouring (deeper) layer. In other words, starting from a given point on the first layer where its slope s is known and using the age difference between the layers, this equation can be used to predict the ‘target’ point on the deeper layer where the same slope should occur. We tested Equation (20) by comparing the layer slope measured at the target point with its expected slope, for multiple positions on the layers on the real and synthetic radargrams. Figure 14a and b show the results. These plots exclude test positions where low isochrone curvature is expected to cause high migration velocities in Equation (20) and, together with the separation between layers, large prediction errors. Agreement between measured and expected slopes is strong for the synthetic radargram (Fig. 14a) and much weaker for the real radargram (Fig. 14b) (due again to the errors mentioned for the last test), but, unlike in Figure 13b, a noticeable correlation can be discerned here for the real radargram.

Fig. 14. Testing slope-migration Equation (20) with data from: (a) the synthetic radargram in Figure 6b and (b) picked layers on the radargram used in our real case study. Horizontal axes plot isochrone slope s measured on layers of the transformed radargrams. Vertical axes plot isochrone slope s measured at positions on deeper layers adjacent to the original layers, positions predicted by Equation (20). In both panels, the grey points and black points, respectively, exclude test positions where absolute value of the isochrone curvature |∂s/∂x| is <5% and <20% of the maximum |∂s/∂x| found on the radargram.