Stratigraphic variation is a well-known effect of spatial variations in accumulation rates (Reference BillingsBillings, 1971). By variation, we mean alterations in horizontal stratigraphy caused by depositional changes and not by deformation processes. Correlations between spatial variations in snowfall accumulation rates and surface slope are well known from Antarctic studies (Reference Black and BuddBlack and Budd, 1964; Reference Gow and RowlandGow and Rowland, 1965; Reference WhillansWhillans, 1975a, Reference Whillansb; Reference Mosley-ThompsonMosley-Thompson and others, 1995; Reference Van der Veen, Mosley-Thompson, Gow and MarkVan der Veen and others, 1999; Reference Vaughan, Corr, Doake and WaddingtonVaughan and others, 1999; Reference HamiltonHamilton, 2004; Reference Spikes, Hamilton, Arcone, Kaspari and MayewskiSpikes and others, 2004). Therefore, it is reasonable to hypothesize that this process, repeated year after year and in concert with ice flow, may be a significant cause of the extensive stratigraphic variation we have recorded in long-distance short-pulse radar profiles of firn (Reference Arcone, Koppenjan and LeeArcone, 2002; Reference Arcone, Spikes, Hamilton and MayewskiArcone and others, 2004, in Reference Arcone, Spikes and Hamiltonpress), as well as in profiles recorded in ice by others (Reference Welch and JacobelWelch and Jacobel, 2003). Such features include the progressive steepening of fold limbs, the unexpected changes in fold-hinge loci direction, and the failure of hinges to move horizontally at the ice-flow rate as they are buried.
Ideally, cores for climate studies should be retrieved where accumulation rates gradually change over hundreds of kilometres so that each core represents a significant area and an undisturbed temporal record. However, almost any spatial sampling area will have minor topography. Although the relative or even absolute surface slopes are small across many regions of Antarctica, they are important. Reference Black and BuddBlack and Budd (1964) used stake measurements every 1.6 km over a 200 km up-slope transect at about 30° to the mean wind direction, and found that deposition generally increased on windward slopes and decreased on leeward slopes. Their slopes varied between −0.6° and 2.3°. Reference WhillansWhillans (1975b) used 50 measurements over 200 km and concluded (Reference WhillansWhillans, 1975a) that deposition could be sensitive to slope changes as small as 0.1°. This spatial bias could then pose a trade-off on any leeward slope between the possible loss of annual layers and the possible gain of many more layers for a given core depth, while any windward slope might offer a trade-off between more accumulation and fewer annual layers. A deeper core would then add these effects because of ice flow. Given this sensitivity of accumulation rates to such small values and changes in slope, radar is then a necessary tool both for tracking the continuity and measuring the depths of horizons to determine accumulation rates, and for obtaining profiles over the very long distances needed to recognize the impact of spatial variability upon the stratification caused by climate. Our purpose is to provide an understanding of the genesis of the spatial variability in relation to topography, wind and ice velocities, so that its presence can be recognized and for which its overprint on temporal variability caused by climate change could, at least partially, be quantified.
Our objectives are to measure and correlate the accumulation rate and topographic slope more exactly than previous studies, and to quantify the dips and accumulation rates within folds that occur in a radar profile of firn stratigraphy. We then use the accumulation rate variability to construct a model to reproduce these fold features. We chose a long profile with good alignment between ice, wind and profile directions, measured the accumulation rates from a recent horizon whose time delay could be calibrated in terms of depth from ice-core density profiles, and derived slopes from the global positioning system (GPS) elevation data. At several locations we then delineated and measured the migration of fold-hinge axes and compared these distances with those expected from the ice speed. At two locations we measured the variation of time-averaged accumulation rates and horizon dips with depth, and compared these measurements with the results of our model. Much of the stratigraphic variation in horizon depth compares with the local elevation changes. Therefore, all profiles are presented, and stratification measurements were made with respect to a normalized surface (as appears in a recorded profile) to provide a uniform reference, remove the distortion added by the surface slope and to facilitate recognition of stratigraphic variation and comparisons with model results.
The profile we discuss was recorded along a 177 km transect in West Antarctica during 1999, the first of the four years of the International Trans-Antarctic Scientific Expedition (ITASE) (Reference MayewskiMayewski, 2003). A major objective of this research program was to obtain widely spaced ice cores to study recent changes in climate. This and all other transects were predetermined, both for core sampling efficiency and for safety reasons, so they do not always line up exactly with (or are orthogonal to) the ice-flow direction. The profile gives the best-known alignment with ice flow and wind direction of all that were recorded during the ITASE program. Our ice-flow vectors were determined over a 1 year period using repeated GPS measurements. Wind vectors were recorded by automatic weather stations (AWSs) operating since 1980 and 1997 respectively at Byrd and Swithinbank stations (located later). The transect included many hills and valleys and was aligned between 27° and 37° of the wind direction and 25° of the ice direction. The variation of ice speed along the part of the transect we analyze precluded the possibility of stratigraphic deformation by compression. In addition, the firn is not old enough to have accumulated much strain history. In obvious contrast with the studies of Reference Black and BuddBlack and Budd (1964) and Reference WhillansWhillans (1975a), we used a GPS to measure elevations and to calculate surface slope more accurately.
We recorded the firn stratigraphy with a 400 MHz short-pulse radar whose pulse duration allowed us to vertically resolve horizons to about 40 cm separation and provided sensitivity to layers possibly thinner than 10 cm (Reference Arcone, Spikes, Hamilton and MayewskiArcone and others, 2004, in Reference Arcone, Spikes and Hamiltonpress). This degree of resolution enabled us to follow particular, dated horizons along which we could track the average historical accumulation rates. Without actually following individual horizons, several innovative studies (Reference Fujita and MaeFujita and Mae, 1994; Reference FujitaFujita and others, 1999; Reference Seigert and FujitaSeigert and Fujita, 2001; Reference Eisen, Wilhelms, Nixdorf and MillerEisen and others, 2003) have shown that higher-frequency radar horizons in firn must result from density variations. Reference Spikes, Hamilton, Arcone, Kaspari and MayewskiSpikes and others (2004) show that these horizons are, indeed, isochrones, as has long been assumed. Reference Arcone, Spikes, Hamilton and MayewskiArcone and others (2004, in Reference Arcone, Spikes and Hamiltonpress) argue that the general variation of horizon amplitude with depth, their −25 to −50 dB reflection coefficients and the consistent phase structure of their wavelets are consistent with the stratification processes caused by the post-depositional snow metamorphosis and sublimation that form couplets of hoar and ice (Reference AlleyAlley, 1988).
Despite some similarities in stratigraphic structure, this stratification process, the lack of erosional unconformities, the advective movement of the firn and its general alignment with wind direction distinguish the genesis of our stratigraphy from that of other sedimentary strata deposited over variable topography, such as seen in seismic profiles of deep-ocean contourites (Reference Michels, Rogenhagen and KuhnMichels and others, 2001; Reference VianaViana, 2001).
2. Description and Location of Transect
The location of our ground-penetrating radar transect and the core sites along it are shown in Figure 1. We recorded the profile on 11–12 December 1999. The transect is superimposed on a RADARSAT image, which has a texture that appears to be related to the topography. The imagery and the location data of Reference Price, Bindschadler, Hulbe and BlankenshipPrice and others (2002) show that the transect crossed tributary D1 of Bindschadler Ice Stream (BIS), starting at approximately 50 km and ending at about 84 km. Our base station for this and all other ITASE transects (Reference MayewskiMayewski, 2003) was Byrd Surface Camp (BSC). The transect extended 177 km southwest from BSC and dropped 560 m in elevation to the Swithinbank AWS (core site 99-2). Beginning at about 70 km from BSC, the transect traversed long, rolling hills 30–60 m in relief and spaced on the order of 10 km. The strike of the first three of the largest hills appears to be oriented about 20° from perpendicular to the transect direction. The orientation of all other hills cannot be discerned from the image. The snow surface was generally smooth, with no large sastrugi encountered during travel. RADARSAT imagery indicated no crevasses, nor were any encountered or detected with a separate radar we used for crevasse detection (Reference Arcone, Delaney and NoonanArcone and Delaney, 2000).
Our GPS data show that the ice flow is parallel with our transect at BSC, 20° to the west of our transect at 99-1, and 25° to the west at 99-2. The mean wind direction is 30° to the east of our transect at BSC (measured daily since 1980; data courtesy of the Antarctic Meteorological Research Center, University of Wisconsin-Madison) and 37° to the east at 99-2 (measured daily since 1997). Our repeated GPS surveys of stakes provide ice speeds of 11 m a−1 at BSC and 48 m a−1 at 99-1 (87 km) and 99-2 (177 km). The velocity data of Price and others (2002), whose GPS grid overlaps the first 100 km of our transect, show that ice speed reaches about 40 m a−1 by 53 km and about 49 m a−1 at about 65–80 km. It then decreases to about 30 m a−1 by 95 km.
3. Radar Methods
We used a Geophysical Survey Systems Inc. (GSSI) model SIR 10 B control unit and the GSSI 1.2 W peak power 400 MHz model 5103 antenna transducer unit. The 4.1 ns duration of our mainly 3/2-cycle pulse gives about 0.41 m layer resolution in firn of refractive index n = 1.5 (density, ρ = 600 kg m−3). We pushed the 5103 unit ahead of a Tucker Sno-Cat® (Fig. 2); in subsequent years we dragged it in a sled (Reference Arcone, Spikes, Hamilton and MayewskiArcone and others, 2004, in press). We recorded differential GPS position every 30 s, corresponding to about 60 m at typical traverse speeds, and used the data to normalize our profile distances between 1 km segments.
We recorded our radar data at 16-bit samples per trace and with time-range gain. Thus our maximum dynamic range in a profile is 96 dB, but the useful range considering noise was probably <90 dB. We recorded at a time range of 600 ns and at 2048 samples per trace. We used a running stack of 32 and a trace acquisition rate of 24 s−1 so that the effective recording rate was about 1 trace per 2.7 m at the traverse speed of approximately 2 ms−1. Profiles recorded in later years extended the time range to 1500 ns (Reference Arcone, Spikes, Hamilton and MayewskiArcone and others 2004, in press), but at a greatly compromised acquisition rate of 1 trace per 15 m, thus sacrificing some data quality.
3.2. Data processing
We alleviated direct-current gain offsets and electronic noise with a wide-band filter, removed constant time-delay nearsurface clutter with a horizontal filter, removed the timevariable gain used during recording, and then applied an inverse range-dependent gain function to compensate for geometric spherical beam spreading losses. The progressive accumulation of two-way electric field transmission losses through each interface in firn is theoretically insignificant.
We used the echo delay formula,
where c = 0.3 m ns−1, to transform echo time delay, t, measured in nanoseconds (ns), into depth, d, in meters. We used an empirically determined relationship (Reference Kovacs, Gow and MoreyKovacs and others, 1996), verified with earlier datasets (Reference CummingCumming, 1952), to transform density, ρ, into dielectric permittivity, e. Equation (1) was applied meter by meter to the densities obtained from our core profiles to produce a slightly nonlinear depth scale that varied insignificantly between core sites. The quantity
Our spatial processing of the profiles consisted of additional stacking (beyond that used during recording), which greatly improved the signal-to-noise quality, and distance normalization using our GPS data to provide an equal number of traces per kilometer. As discussed earlier, we do not correct for elevation, so that the variations in horizon depths relative to the surface are clear. With elevation corrections, it would be difficult to see these depth variations along the stratigraphic dips because they compare with the variations in elevation. The extremely shallow stratigraphic dips precluded any need to migrate our profiles. However, the stacking weakened reflections along steeper dips, as will be discussed.
3.3. Antenna considerations and reflection dynamics
The complexities of the antenna directivity are not important, because of the gentle slope of the layers. We calculate that our reflections originated within an approximately coneshaped volume from beneath the antenna, given the gentle dip of the stratigraphy and the width of the first Fresnel zone. for our situation, the Fresnel zone is the circular area over which distance from the antenna does not vary by more than about λ/8 (0.06 m in firn of ε = 2.6) and has a diameter where d is depth beneath the antenna and λ is the in situ wavelength. At an average ε = 2.6, D is only about 5.2 m at 56 m depth. The recording rate and the Fresnel zone determine the theoretical distance over which the reflections are integrated. This distance varied from 4.3 m at a depth of 5 m to 7.8 m at a depth of 56 m.
4.1. General profile and accumulation rates
We show a 350:1 horizontally compressed version of the entire profile in Figure 3. An important variation feature that we discuss later is the general increase in the dip of the limbs with depth within any one fold. This feature is strongest at 70–130 and 150–177 km, which are the regions of higher ice velocity. A second important feature discussed later is the general failure of the loci of a stratigraphic fold to monotonically migrate down-ice (to the right in Fig. 3) with depth, at the rates expected from the ice speeds. A third feature is the significant change in depth of some horizons over long distances. The horizon marked by vertical up-pointing arrows varies from >60 m depth near 30 and 40 km, where it dips below the record, to 30 m depth at 177 km. This horizon is dated to the year 1734 (personal communication from D. Dixon, 2003). A second horizon, marked by smaller arrows and dated to 1875, shows a similar variable-depth character.
Figure 4 shows the calibration data we used to convert the time delay of any particular horizon in Figure 3 into a cumulative water equivalent (w.e.) value in m a−1. The time delay was then converted to depth using the meter-bymeter procedure discussed in section 3.2, and then to years from the dating. The BSC calibration in Figure 4a was used for the computation of the near-surface accumulation rates in Figure 5a, while all calibrations and an interpolation (Fig. 4b) were used for the deeper rates of Figure 5b.
The accumulation rates plotted in Figure 5a are 8 year w.e. averages for 1992–99 (start of 1992 to the end of 1999, at the surface). We picked the time delays at every kilometer along the leading edge of a horizon half-cycle that we dated to the start of 1992 using the ice-core data from core site BSC (2.5 m depth; Reference Kreutz, Mayewski, Meeker, Twickler and WhitlowKreutz and others, 2000) and from core site 99-1 at 87 km (4.2 m depth; personal communication from D. Dixon, 2003). We assume that the near-surface density profiles are consistent along the entire transect because cores from 0 km (BSC) and 177 km (99-2) show similar near-surface rates of densification as a function of depth (Fig. 4a). Although density was not obtained for the first 4 m at core site 99-1, the 1992 horizon at BSC does track to about 4.2 m depth at this location.
We chose this shallow horizon so that each accumulation rate would correspond with the topography at that location, whereas deeper horizons include accumulation from some distance up-glacier. The horizon appears to bifurcate at two locations, but where it does, the correct continuation was guided by both the amplitude and general dip of neighboring horizons. Our accuracy in measuring time delay is estimated from the width of a wavelet half-cycle (about 1.3 ns; Fig. 2). Therefore, the error in our accumulation rates varies from a maximum of 6.8% for the minimum time delay of 14.7 ns (at 118 km) to a minimum of 2.5% for the maximum delay of 49.3 ns (77 km). The estimated average error caused by timing inaccuracies is then 4.7% for the entire profile. The accuracy in the dating is about 0.1 year, which increases the 8 year w.e. average error to about 6%.
We calculated along-profile slopes from the GPS elevation profile. The slope at each km represents an average over ±1 km from that point, and we did not smooth the data. The maxima in the plot are not necessarily the absolute maxima because our direction of travel relative to the wind was off by as much as 37°. Slope maxima indicate windward slopes facing up-glacier (to the left) to the northeast, and minima indicate leeward slopes, facing down-glacier (to the right) to the southwest. On average, the transect dipped 0.18° to the southwest. Consequently, leeward slopes are always negative, but some maximum slopes, on windward sides of hills, may still be negative.
The correlation between slope and accumulation rate is strong, regardless of slope value. In some cases, even small perturbations in slope, such as at 106 and 114–123 km, are in phase with corresponding perturbations in accumulation rate. At a few locations, such as at 12 and 136 km, they are out of phase by 180°. At other locations, such as at 40–50 and about 76.3 km, they are out of phase by about 1 km, which may result partially from the 2 km method of computing the slope.
Figure 5b shows plots of average accumulation rates calculated for 1875–1999 and 1734–1999. The rates are based on density values that occur at depths which correspond with time delays of 200–580 ns. We used the density profile at BSC for 0–22 km, an interpolation between the BSC and 99-1 profiles for 23–65 km, the profile at 99-1 for 66–173 km, and the 99-2 profile for 174–177 km. Core site 99-2 appears unique because the density reached 0.82 kg m−3 at only 38 m depth. The radar horizons distinctly change in amplitude below this depth (discussed later), and the elevation drops into a basin after about 173 km.
In Figure 5b there is considerable variation for the two horizons. The flat tops within the 1734 horizon represent sections where this horizon dipped below the profile depth. Between 5 and 70 km the 1875–1999 rates vary from about 0.12 to 0.19 m w.e. a−1, and from 1992 to 1999 they have a similar range of about 0.15–0.21 m w.e. a−1. The six prominent maxima for 1892–1999 between 70 and 125 km in Figure 5a also occur in Figure 5b at almost the same locations, but the 1875–1999 range of about 0.085–0.185 m w.e. a−1 is lower than the 1992–99 range of 0.08–0.28 m w.e. a−1. However, from 125 to 177 km the 1875–1999 rates vary from about 0.06 to 0.15, the 1734–1999 rates vary from 0.06 to 0.13 m w.e. a−1, and the 1992–99 rates vary from about 0.09 to 0.135 m w.e. a−1. Therefore, although the tracking of maxima and minima is generally comparable for all horizons to at least 300 years, the larger, long-term ranges of variation between 70 and 125 km appear to diminish with depth. This is reasonable because any rate determined at a single location is an average over an increasingly longer distance upstream as depth increases at that location, and so the variability should be reduced as the depth increases. This smoothing of long-term accumulation rates will be verified in section 5.5 where we compare model and profile rates vs depth.
4.2. Details: 0–45 km
The profile segment presented in Figure 6 contains relatively gentle topography, with slopes varying between about 0.1° and −0.5° (Fig. 5a). We measured the ice speed at BSC at 11 m a−1. According to Reference Price, Bindschadler, Hulbe and BlankenshipPrice and others (2002), the speed at 45 km is about 30 m a−1. Therefore, the maximum down-ice movement at the bottom of this section (about 300 years) should be ~3 km at the left side and ~8 km near the right side of the profile, considering that the transect direction is no more than about 20° to that of the flow.
The full-depth profile in Figure 6 shows several folds, some of which appear to migrate down-ice. However, any characteristic feature of the folds, such as the loci of hinges (indicated by the dashed lines), shows either no migration or no more than half of what is expected. for example, the hinges starting at 32 and 35–36 km migrate only about 3 km by the bottom of the profile. As shown later by our modeling, this lack of expected movement is caused either by relatively isolated accumulation anomalies or by the effect of ice movement upon the periodicity within the accumulation variation. In comparison with the next 45 km, the periodicity in the accumulation rate is weak, and so the anomalous behavior may dominate in this section.
4.3. Details: 45–90 km
From 45 to 70 km the topography continues to be relatively gentle (Fig. 7). Over this section, slopes vary from about 0.15° to −0.25°. The transect crosses into BIS at about 50 km. After about 67 km the elevation drops steadily and then fluctuates more strongly from 70 to 90 km where the slope varies between about ±0.8°. The hills and valleys continue into the next section and have topographic wavelengths varying from about 8 to 11 km, with an average near 9.6 km.
After 70 km the topography becomes more variable. At 75 km the most negative leeward surface slope of about −0.8° is encountered, along with a sharp accumulation minimum at about 0.12 m w.e.a−1. At 77 km, just slightly beyond a valley bottom, the slope passes through 0°, elevation begins to increase, and the 8 year accumulation rate reaches 0.28 m a−1 (Fig. 5a). At this location the accumulation peak lags behind the peak of the slope by about 2 km, but all subsequent maxima track almost perfectly with slope maxima, as do the rate minima. The accumulation rate remains strong on this complex rise, then drops to a minimum at the minimum in the slope situated at 82 km, and then rises to a maximum at 88 km, just slightly beyond the slope maximum (Fig. 5a).
The ice speed in this section varies from about 30 m a−1 at 45 km, to at least 48 m a−1 from 65 to 87 km. These ice speeds should give a range of displacement for any variation feature of 8–13 km over 300 years, considering that the transect is angled to the ice-flow direction by about 20°. However, although the variation in Figure 6 becomes relatively strong after 70 km, there is no appreciable or consistent migration to the southwest with depth for any of the variation features. for example, the profile detail shows that the synclinal fold-hinge locus starting at about 77 km is vertically stable to about 10 m depth. It then moves about 2 km down-glacier to about 27 m depth, and then is stable to 56 m. Within this fold, the horizon that starts at ~10 m depth at 73 km descends to ~25 m depth by 78 km.
An unusual, lens-like feature occurs at about 43 m depth between 85 and 90 km. In this vicinity the loci of fold hinges noticeably change direction with depth, a feature which occurs often in the rest of the profile. Along any vertical section, the dips of the fold limbs are generally always positive or always negative, with both cases including values of 0° relative to the surface. This feature might appear to represent a single anomaly in accumulation rate or an unconformity. However, our modeling (section 5) will show that it results from the continual superposition of horizons whose spatially varying accumulation rates are progressively out of phase with each other because of ice motion.
4.4. Details: 90–135 km
Within this segment (Fig. 8), surface slopes vary from −0.9° to 0.5° between 90 and 115 km and then vary from about −0.6° to 0.2° between 115 and 135 km. Both the detail in the middle profile and Figure 5a show how well the accumulation rate tracks with the slope until about 128 km. A rare anticorrelation between surface slope and accumulation rate is centered at about 135 km.
According to Reference Price, Bindschadler, Hulbe and BlankenshipPrice and others (2002), ice speed at about 92 km is 40 m a−1, and drops to about 25 m a−1 by about 100 km. This latter speed could be expected to give a fold-hinge locus displacement of about 7 km near the bottom of the profile, again considering the 20° misalignment between transect and flow directions. However, the consistent migration with depth of the fold-hinge locus to the southwest at this distance is only 2.5 km. Although the velocity gradient is about −2 m a−1 between 92 and 100 km, any compressional folding caused by this gradient should be an insignificant cause of these direction changes, in view of the strong fluctuations in accumulation rate along this section, and is unlikely because the firn is too young to have developed much strain history. At 125 km, where the ice speed may drop to 20 m a−1, the locus displacement is only 3 km instead of the expected 5–6 km. More common are loci of fold hinges whose directions vary about the vertical near 94, 111, 120 and 129 km.
The radar failed to record stratigraphy in a few places. The blank section below ~35 m depth near 97 km is artificial. At this location the magnitudes of the stratigraphic dips exceed about 0.4° and were sufficient to cause destructive interference between consecutively stacked traces. Reference Arcone, Spikes, Hamilton and MayewskiArcone and others (2004) give quantitative analysis as to why this is the case. However, weak reflections below ~44 m depth at 108–112, 115–121 and 131–135 km are associated with dips of less value. These zones of weak reflections continue near the bottom of the next section (Fig. 9;135–146 km) and are most conspicuous near its end where our core density profile revealed that the firn/ice transition was reached at only 38 m depth (discussed below), and where the accumulation rate was low. Therefore, we speculate that the edges of these zones also represent transitions to ice of density greater than about 820 kg m−3, below which density contrasts would be greatly diminished. The arching horizons above them suggest they were also generated in areas of low accumulation rates.
4.5. Details: 135–177 km
The surface slopes in this segment (Fig. 9) vary from about 0.2° to −0.9°. As can be seen in the detail profile or from Figure 5a, the accumulation variations track very well with the slope. We cannot estimate ice speed throughout this section. However, the horizon depth variations between 135 and 155 km are far less severe than between 70 and 135 km. Based on our modeling presented next, this is consistent with the decreased variability in accumulation rates over this section, and with an ice speed that is probably about 20 m a−1. After 155 km the degree of the variation increases, and so does the range of the accumulation rates. At some point the ice speed does also, and we later estimate that by about 160 km it is about 25 m a−1.
The righthand side of this segment ends at core site 99-2, which was located in a basin where we measured ice speed at 48 m a−1. The ice core shows the firn/ice transition at a relatively shallow 38 m depth (unpublished data), and the location itself has a very low accumulation rate (Fig. 5a). We traced the age of the horizon at this depth back to core site 99-1 and found it is about 300 years old. Therefore, approximately 14 km of displacement after 300 years might be expected for the locus of any fold feature here. However, as with the previous sections, there is no consistent migration down-ice of any variation feature, and the few loci that we traced for the fold hinges are primarily vertical.
5. Interpretation from Modeling
5.1. Model parameters
We seek to explain the unusual profile features noted above with a simple model. We take the glacier radar surface to be flat, as viewed in a radar profile, and determine the thickness, y, of an annual deposition with a cosinusoidal accumulation rate such that
The coefficient A is the amplitude of the variable part of the accumulation rate and C is a constant, background rate. We chose A = 0.125 m a−1 and C = 0.225 m a−1 for snow at density ρ = 400 kg m−3. These parameters assure there is always horizon separation and that the rate will always be positive because it varies between 0.1 and 0.35 m a−1. The quantity x is distance, v is ice speed (such as 10 m a−1), a is time in years, and λ is the wavelength of topographic variation, which we assume to be 10 km (peak to peak) to simplify comparisons. The actual average hill spacing near 99-1 is about 9.6 km.
We make no assumptions regarding the stability or shape of the surface topography. We seek only the effects of a systematically displaced series of sinusoidally varying thickness functions upon a radar record. We do not correct our model profile for elevation, because this would make the effects we seek less clear, and difficult to compare with the radar profile. Consequently the introduction of any distance phase lag between the thickness function and the surface slope would be meaningless because we only seek the effects of a variable accumulation rate and ice movement at one point in time. Migration of the surface topography (Reference Black and BuddBlack and Budd, 1964; Reference WhillansWhillans, 1975a) would complicate a time evolution picture of the stratigraphy, but whether this truly happens is debatable (Reference WhillansWhillans, 1975a).
The first accumulation surface is measured relative to the instantaneous topographic surface and is followed by succeeding surfaces at depths that are characterized by greater ages, a, and greater displacement down-ice, va. The depth, y n, to the nth surface at any distance x is the addition of all the layer thicknesses above it, and is expressed mathematically as
This summation produces some interesting features that will be seen in the examples to follow.
The initial density of 400 kgm−3 for all layers is a good approximation for that of the first 10 m of snow, which is the thickness of all our modeled layers. We then multiplied the depths, y n, by a densification (or compaction) function, which simulates the densification at site 99-1, to find a new depth , and applied this function to all of our models. The mathematical form of this empirical function is
where y n and are negative numbers. The expression in braces was derived by integration of an empirical exponential density function,
to derive the total mass at a depth y. Although this compaction consideration provides the approximately correct depths of 40–50 m for the 300 year surface, it marginally affects the qualitative results discussed below. Our densification function is continuous and provides no abrupt density contrasts. Therefore, the individual horizons we plot could represent very thin density discontinuities, such as those caused by the millimeter-thick ice crusts and associated hoar layers that occur so frequently in ice cores (Reference GowGow, 1968; Reference Arcone, Spikes, Hamilton and MayewskiArcone and others, 2004).
5.2. Changes in uniform ice speed
Figure 10a–c show cases where the ice speed is steady. In Figure 10a the speed is relatively low at 10 m a−1. The accumulation rate differences appear to propagate down-ice uniformly with depth, the dips of the folds become steeper with depth, and the loci of the fold-hinge apices appear to form straight lines. However, after 300 years the deepest hinge axis has moved only 1.5 km down-ice instead of the 3 km expected from the velocity, and reproduces the unexpectedly reduced migration of fold hinges down-ice seen in the profiles. This reduced migration is explained mathematically by either a simple trigonometric argument or an integration of Equation (3). Trigonometrically, the addition of two or more sinusoids, each with a progressively increasing incremental change in phase, ϕ, produces a new function with half the phase shift. By trigonometric identities
Figure 10b shows the case where ice speed is 30 m a−1. The down-ice migration of the last surface is 4.5 km instead of the expected 9 km. The straight line shows that a locus connecting the fold hinges would deviate from linearity, so that the fold feature is beginning to deform. In addition, the variation in depth of the lowest horizon is far less than in Figure 10a, so that the glacial speed relative to the accumulation wavelength smooths the variation.
Figure 10c shows the case of v = 50 m a−1. The variable direction of the fold-hinge loci is now significant. This ice speed gives one wavelength (10 km) of ice travel after 200 years. The 200 year surface is the horizontal (relative to the surface) horizon at about 29 m depth. Along any vertical section near and across this horizon, the hinge orientation periodically changes from anticlinal above to a down-ice displaced, synclinal below. As shown previously, changes in direction of the fold-hinge loci occur at 80–90, 103–135 and 157–167 km. We model the one between 80 and 90 km below.
5.3. Anomalous changes in ice speed
Figure 11a presents another possible way in which the fold-hinge loci could change direction between relatively up- and relatively down-glacier. In this example, the speed varies with age, starting out at 50 m a−1 for the most recent 100 years, then 10 m a−1 for the second, and finally 100 m a−1 for the third most recent 100 years. These changes translate into additional phase lags in Equations (2) and (3). for example, at year 210, the total translation results from 100 years of motion at 50 m a−1, plus 100 years of motion at 10 m a−1, and 10 years of motion at 100 m a−1. Rather than produce more down-ice translation at the bottom, the pattern appears to produce up-ice translation. for example, the loci of the anticlinal fold hinges that start at 5 and 15 km form a curved line that first slopes down-glacier, then tends toward vertical at about 25 m depth, then starts to slope back up-glacier at about 29 m depth, and then tends back towards vertical at 32 m depth.
The 10 m a−1 speed model of Figure 10a predicts the largest depth fluctuations of the deeper horizons in Figure 3, but the model is unrealistic because such slow speeds occurred only in the first quarter of the profile where topographic and accumulation variations were relatively small. The more realistic model of Figure 10c predicts depth variations of about 7.5 m for any horizon older than 120 years. However, between 70 and 130 km the horizon dated to 1875 varies as much as 14 m, and that to 1734 much more. Although the variable-speed model of Figure 11a produces up to 10 m change in horizon depth after 300 years along with variable directions in fold-hinge loci, and many speed variations are possible, any significant change in speed should cause significant distortion throughout the entire profile in Figure 3, which is not the case. In addition, there is no other evidence that such changes have ever occurred in such a short time-span. Therefore, we think that the situation depicted in Figure 10c provides a more reasonable explanation for loci direction changes because it assumes a single glacial speed.
5.4. Anomalous change in accumulation rate
We have argued that a periodically varying accumulation rate precludes tracing the burial history of an ice particle along the loci of any particular feature, such as a fold hinge. Figure 11b shows that the correct burial trajectory can be traced from an accumulation anomaly localized only in space and superimposed on a uniform deposition rate, as might occur year after year on an isolated windward, monocline surface. A nearly identical model was developed and used by Reference Leonard, Bell, Studinger and TremblayLeonard and others (2004) to explain an accumulation anomaly seen in the Vostok (East Antarctica) ice core. As with the previous models, the surface is normalized to horizontal and we assume each successively deeper layer had the same anomalous accumulation distribution when it was at the surface. The anomalous peak amplitude of 0.25 m a−1 is superimposed on a constant rate of 0.35 m a−1, all for snow at 400 kgm−1 density. The anomaly is represented by a displacement y a, and has a Gaussian distribution expressed as
where v =50 m a−1. This distribution is seen most clearly in the first accumulation horizon. Each horizon represents 10 years of accumulation. The depths of the successive 10 year surfaces are summed as in Equation (3). The anomalous deposition is centered at 4.5 km, as expressed by the delay of 4.5 km in the exponential. The velocity then centers the anomaly at 5 km distance after 10 years.
The maximum displacement of any feature within the 300 year surface might be expected to be 15 km after 300 years. However, only the leading edge of the anomalies (delineated by the dashed line in Fig. 11b) propagates at the ice speed. In the lowest surface it has moved 15 km down-ice. Other features have spread with depth but are displaced far less. Therefore, trajectories that indicate the correct ice speed could only be recovered from the leading edge of local anomalies that are well isolated from any other variations in accumulation. The spread of the variation at depth is reasonable because the displacements of the accumulation disturbances in the recent, upper layers translate through to the lower ones.
Radar profiles of down-glacier spreading of stratigraphic variability with depth caused by a single anomaly are given by Reference Vaughan, Corr, Doake and WaddingtonVaughan and others (1999, fig. 2c) and by Reference Leonard, Bell, Studinger and TremblayLeonard and others (2004). Vaughan and others show that an anomaly of about 0.5 km width at the surface becomes almost 1.3 km wide at about 96 m depth. Our modeling suggests their interpretation of ice speed from the locus connecting the crests of the synclinal hinges is approximately correct since it does not differ significantly from a locus joining the down-glacier, leading edge of the disturbance, while the up-glacier edge appears vertical.
Figure 11c shows a similar, Gaussian-shaped anomaly that is localized in time (120–130 years BP), as well as in space (centered at 12 km distance), and is superimposed upon a periodic accumulation rate and a steady ice flow at 30 m a−1. Mathematically, this anomalous thickness, y a, is expressed as
The anomaly doubles the accumulation rate to a peak of 0.7 m a−1 during a 10 year period. As in Figure 11b, the extra accumulation delays the reflections from all subsequent layers. Consequently the anomaly distorts the stratigraphy beneath it and creates a vertical series of secondary fold hinges. A similar effect would result if the anomaly were expressed as an increase in density, which would also supply extra time delay.
There is no indication in our data or modeling of any feature that consistently tilts down-glacier, starting from the bottom of the glacier, as suggested by Reference Whillans and JohnsenWhillans and Johnsen (1983) and Reference Jacobel, Gades, Gottschling, Hodge and WrightJacobel and others (1993). In these cases, they surmised that the cause of the folding occurred far upstream, where basal boundary conditions changed during the transition from inland ice to ice-stream flow.
5.5. Comparisons with profile features
We recreate two features that appear to be central expressions of the stratigraphic variation process. They are (1) steepening fold dips associated with asymptotically vertical hinge loci; and (2) periodically changing fold dips associated with fold hinges whose loci directions waver about vertical. We first attempt to reproduce the slopes and the accumulation rates associated with the fold limbs located near 77 km (within BIS) in Figures 3 and 7, with the isolated accumulation anomaly model of Figure 11b. We chose this feature because the change from a dipping to a vertical direction of the fold-hinge locus at about 25 m depth appears similar to that of our model, and because the ice speed is nearly 50 m a−1 across this section. The synclinal hinges begin near the start of the hilly area crossing BIS and there is little topographic variation up-glacier from this section. The section is not yet in the midst of where the strong, periodically varying accumulation rate occurs. Consequently, this feature should have some aspects of an anomaly localized in space but not in time, because the hillside is a likely permanent feature in response to sub-ice topography.
Figure 12a and b compare the model with the feature. The vertical dashed lines are where the comparisons between accumulation rates (Fig. 12c) and slopes (Fig. 12d) were measured. In Figure 12b the vertical line is placed at 76.5 km where the anomaly appears to originate, and about 260 years of accumulation are considered. Similarly, the vertical line in the model is placed where the anomalous accumulation is centered at the surface. The model uses a densification function that simulates the density profile at 99-1 (at 87 km), and the maximum age is 300 years. As mentioned earlier, in the real profile the horizon at 10 m depth at about 73 km drops to about 25 m depth at 78 km. In the model, the 12 m horizon on the left side drops to about 23 m depth.
The yearly-averaged w.e. accumulation rates in Figure 12c are referenced to the surface (end of year 1999). We obtained horizon dates by tracking the dated horizons from site 99-1 at 87 km. The accumulation rate function (given in Fig. 12a) was chosen to generate the measured rate of 0.28 m a−1 for the 8 year old surface at 77 km (Fig. 5a). Qualitatively, the model variations track those of the profile. Quantitatively, the model gives similar peak values, as it should, and a similar asymptotic rate of about 0.12 that is consistent with the calculations at this distance in Figure 5b. The rates decrease from a peak value because lower-accumulation-rate snow has moved beneath the peak of the accumulation. Figure 12d shows that the model closely simulates the relatively rapid decrease in slope, from 0° (the normalized surface) to about 0.3° at 15–16 m depth, the knee in the slope variation, and then the near-steady value of about −0.22° to the bottom of the record.
Figure 13a and b compare the model in Figure 10c with a segment of the profile located near core site 99-1 at 87 km. We chose this feature because the change from anticlinal folding relatively up-glacier, to synclinal relatively down-glacier at about 4 km before the topographic peak, occurs in both model and profile, and because the ice speed across this section is close to 50 m a−1. The model uses a sinusoidal accumulation rate amplitude of 0.125 m a−1 superimposed upon a steady rate of 0.30 m a−1. A horizontal (relative to the surface) horizon occurs at about 38.3 m depth (time delay of 386 ns) in the profile. This depth dates to about 1813 (personal communication from D. Dixon, 2003), or an age of 186 years. As discussed above, this horizon should occur at about 200 years BP (1999) given the glacial speed of about 50 m a−1 and the accumulation wavelength variation of about 10 km.
The slight discrepancy in years between the model and profile horizontal stratum is explained by the slightly less hill spacing of the profile data. The age, a n, in years of the horizontal line is found from the simple formula
At a wavelength of 10 km and speed of 50 m a−1, a n = 200 years. If the wavelength is reduced to the actual value of 9.6 km, then a n = 192 years, which improves the agreement with the data, but at a speed of 48 m a−1 it is still 200 years. We do not think that this computation should consider the vector component of the speed in the direction of the transect (46 m a−1), because the depth of this horizon would then steadily increase to infinity as the transect changed direction.
Figure 13c shows good agreement between the model and the profile yearly-averaged w.e. accumulation rates, and peak and asymptotic values at 86.5 km. The asymptotic value of about 0.11–0.12 agrees with the calculations for this distance in Figure 5b. The w.e. values decrease from the peak because lower-accumulation-rate snow has moved down-glacier and beneath the peak of the accumulation. Figure 13d shows good agreement between fold slopes, which vary sinusoidally and are always negative.
A second possible example of this phenomenon appears in the boxed area in Figure 9 between about 158 and 165 km. It occurs at about 35 m depth at the midpoint of the fold-hinge locus that changes direction. This segment was recorded within a series of four hills whose average peak-to-peak spacing is 4.9 km. From 160.5 to 162 km the folds change from anticlinal to synclinal going down-glacier, with an intervening horizon at about 35 m depth. This horizon can be traced back to core site 99-1 where it coincides with the same 200 year horizon discussed in Figure 13. Using Equation (9), the ice speed at this location should then be about 25 m a−1. Although ice speed increases to 48 m a−1 by another 16 km, this is not an unreasonable value given the decrease in ice speed seen in the data of Reference Price, Bindschadler, Hulbe and BlankenshipPrice and others (2002) after 110 km.
The strong correlation between slope and accumulation rate occurs for almost all slopes along our transect and confirms the earlier work of Reference Black and BuddBlack and Budd (1964) and Reference WhillansWhillans (1975b), but on a more continuous scale. Our calculations show no threshold slope value for this process. The large ranges of accumulation rates for all horizons, and the similar periodicity and placement of maxima and minima provide evidence that at depth in this region of West Antarctica the ice motion requires >300 years of deposition to provide long-term, yearly-averaged accumulation rates that are significantly smoother with distance than those of the near-surface. Certainly, greater depths imply less spatial variability in long-term averages.
The consequence of this process is strong stratigraphic variability when viewed on a length scale of tens of kilometers. We conclude that the primary causes of the variability are the effect of ice movement upon periodically variable accumulation rates, and local anomalies in accumulation. We also conclude that tracing the loci of variation features is not a correct procedure for deriving ice speeds unless the accumulation anomalies are isolated.
Our model provides a possible tool to assist the interpretation of the depositional history of an ice core and of a region in general. Features such as large dips and vertically aligned fold hinges should indicate strong surface accumulation anomalies, while wavering fold hinges should indicate spatially periodic accumulation rates. Improvements to the model should result from input accumulation rates that were recorded along transects more aligned with the wind direction, from better continuity in the GPS velocity data, and from modeling parameters that allow for a variable topographic wavelength and ice speed. The transect directions relative to those of the wind and the ice flow would affect the appearance of the stratigraphic variability in the profile. If we had crossed perpendicularly to the wind and ice-flow directions, then accumulation might have correlated better with elevation change, rather than with slope (Reference SpikesSpikes, 2003; Reference Spikes, Hamilton, Arcone, Kaspari and MayewskiSpikes and others, 2004). In our case, we were within 25° of the ice-flow direction and so the ice speed along the transect was within 10% of the maximum value.
The larger fluctuations of horizon depth are primarily caused by anomalous accumulation rates. Such fluctuations within the firn regime often appear to grow with continued burial into the englacial regime, and should superimpose upon the deformations that result from convergent flow over subglacial topography. Unfortunately, the resulting fold dips may be too steep to track with the slow recording rate of present high-frequency radar systems when using high rates of trace stacking needed to suppress noise, and must await developments in faster digitization rates of real-time signals. Since these density horizons are isochronous, they are also characterized by the chemically caused conductivity changes sensed by lower-frequency and deeper-probing radars. Ironically, because of their lower frequency, they can record traces at a far faster rate, and so, at the cost of horizon resolution, might be able to track variable stratigraphic horizons more continuously at depth.
US National Science Foundation awards 9814589 and 00880355 (to S.A.A.) supported this work. We thank our program manager J. Palais for her support, our fellow ITASE members P. Mayewski and D. Dixon for supplying ice-core data, and J. Friddell, E. Williams and M. Johnston for graphics and processing support. We also thank R. Bindschadler and an anonymous reviewer for their helpful comments and suggestions.