Recumbent folding of divide arches in response to unsteady ice-divide migration

Abstract

Arches in stratigraphic layers directly under a flow divide (Raymond bumps) are predicted by models of steady ice-sheet flow, and have been observed in several ice domes. Here, we model the evolution of these layers when a formerly stationary divide migrates rapidly to a new position, then again becomes stationary, leaving the arched layers in a flank position. As they are then carried downstream with the flow, these abandoned arches can develop into recumbent folds. These folds can occur over a wide range of divide migration speeds. The shearing flow that produces these recumbent folds also distributes the folded layers over a wide distance downstream from the original divide location. If the divide offset is abrupt, ‘pre-cores’, or material lines comprising core-relative isochrones, can be used to quickly identify portions of an abandoned Raymond bump that would be overturned at any future ice-core site downstream. If, as appears to be the case in Greenland, the divide is never stable long enough to produce a mature arch, folds of this type would not occur. The most likely place to find such folds might be the flank of an ice ridge bounded by unsteady ice streams.

1. Introduction

Over large areas in the polar ice sheets, internal layers have been detected by ice-penetrating radar. Radar energy is reflected at contrasts in dielectric properties. In the upper few hundred meters, density variations are a leading cause. At greater depths, ice-fabric contrasts and electrical-conductivity contrasts cause most or all reflections. The internal layers are widely interpreted to be isochrones, because conductivity and density contrasts usually indicate climatic boundaries. Although fabric contrasts usually develop after deposition, they too can often be tied to climatic boundaries. Layers are often readily observed in the upper 50–90% of the ice depth, where the strain associated with flow is relatively uniform. At greater depths, layers are more difficult to detect due to geometrical power loss and attenuation, and possibly because strain can be inhomogeneous in the deepest parts of the ice sheets.

Arches in stratigraphic layers directly under an ice divide have been predicted by models of steady-state ice-sheet flow. Models of ice sheets using power-law rheology produce distinctive velocity patterns within a few ice thicknesses of a divide that is frozen to its bed (Raymond, 1983). Since the effective viscosity of a power-law material increases with decreasing deviatoric stress, ice deep under such a divide tends to be more viscous than elsewhere. Calculations of isochrones near ice-sheet divides, using thermomechanical finite-element models based on such a non-linear rheology, produce a distinctive arching of stratigraphy in a region approximately two ice thicknesses wide (Hvidberg, 1993, 1996; Hvidberg and others, 1997). This arch has been called a Raymond bump (Hindmarsh, 1996; Vaughan and others, 1999).

Arched layers under divides have also been observed with ice-penetrating radar in several ice sheets (Nereson and others, 1998; Conway and others, 1999; Vaughan and others, 1999). For example, Nereson (1998) and Nereson and others (1998) estimated the rate of divide movement from the inclination of the axis of arched layers observed in Siple Dome, West Antarctica. They also examined the possibility that the arch resulted from reduced accumulation due to wind scouring over the divide. While their work favored the rheological explanation, they could not rule out low accumulation as the cause. In this paper, we start from arches produced by rheological contrasts, but our approach would also apply to arches caused by differential accumulation. Divide movement has been modeled (Anandakrishan and others, 1994; Marshall and Cuffey, 2000). This movement has been suggested as a reason why no arch is observed at the Greenland summit (Hempel and Thyssen, 1992; Hindmarsh, 1996; Hvidberg and others, 1997).

Changing geometry of an ice sheet was proposed by Hudleston (1976) as a cause for folds observed at the margin of Barnes Ice Cap, Nunavut, Canada. Advance or retreat of the ice-sheet margin could change the orientations of streamlines relative to bed-parallel foliation planes. Under these altered flow conditions, the foliation could be sheared, creating recumbent folds some distance further downstream. A similar analysis can be applied to layers higher above the bed and near the ice-sheet center. Waddington and others (2001) viewed stratigraphic folding near ice-sheet centers as a two-stage process. Initially the flat stratigraphy is disturbed locally by small-scale dynamic processes, forming gentle open folds. Then the kinematics of the large-scale flow shears these folds, overturning appropriately oriented limbs. In Jacobson and Waddington (2004) (see also Jacobson, 2001) we extended this idea by actually following the finite kinematic deformation of fold limbs as they moved along particle paths. Using a kinematic model of steady ice-sheet flow, we showed that, under some circumstances, a gentle open disturbance of the stratigraphic layers could be deformed into an overturned fold that would disturb the order of the layers in an ice core located downstream. We introduced the concept of pre-cores, which are material lines that would become vertical at the core site at a specific time in the future. These core-relative isochrones allowed us to determine the minimum slopes that limbs of open folds at any location would need to have in order to develop into recumbent folds at a later time at a particular observation site (e.g. in a future ice core).

We did not present a comprehensive model of the transient dynamic processes that could generate the initial disturbances. For some approaches to that problem, see Azuma and Goto-Azuma (1996), Castelnau and others (1998) and Thorsteinsson and Waddington (2002). In this paper, we evaluate the possibility that a divide arch might be a seed for this kinematic folding, if the divide subsequently moves. Such a move could put the abandoned arch into a flow regime dominated by bed-parallel shear, causing it to eventually overturn.

Here we modify the kinematic flowband model from Jacobson and Waddington (2004) to replicate the velocities calculated by a finite-element model that produces divide arches in steady-state ice sheets. This allows us to compute conventional stratigraphic isochrones (including Raymond bumps) near ice divides. By further modifying our kinematic model to handle an evolving ice-sheet surface, we can track the evolution of stratigraphic layers, as a previously stationary divide moves smoothly for a limited distance. This allows us to explore the impact of divide migration velocity on the likelihood of creating recumbent folds, when an abandoned Raymond bump was the source of initial layer disturbances.

In addition, we can use the kinematic model to find precores. These pre-cores can provide a relatively simple tool to predict which arched layers in an abandoned Raymond bump will become recumbent folds, and which will not, when a divide moves. The major controlling factors are (1) duration of steady divide behavior prior to an abrupt divide move, (2) divide-offset distance during the move, and (3) location of the ice core of interest.

We do not attempt to model folding that could be observed in a core drilled at the summit of an ice sheet. While the seed disturbance for our folding is created under a divide, the folding itself occurs some distance downstream. However, it is conceivable that the divide could later migrate back over a region containing layers previously folded by this mechanism.

2. Geometry and Velocity Model

The coordinate directions in our flow model are x, horizontal aligned with the flow, and z, vertical. The corresponding velocity components are u and w. Although this is not essential to our approach, for simplicity we assume plane strain, with a constant flowband width, so that strain rates perpendicular to the flowband are zero. The surface and bed profiles, S(x) and B(x), are chosen to approximate the geometry of the north flank of Siple Dome (Nereson, 1998). Ice thickness is given by h(x) = S(x) – B(x). With a uniform net accumulation, b = 0.1 ma^–1, and divide thickness H = 1000 m, the model characteristic time is T = H/b = 10 000 years. Ice flux Q(x) is defined as the integral over the ice thickness of the horizontal velocity. The flux, Q(x), is zero at the position of the divide, x _div. Our flowband geometry and coordinate system are illustrated in Figure 1. The horizontal velocity is calculated from

(1)

Fig. 1. Isochrones for a steady-state ice sheet. The heavier lines mark intervals of 1 T. The flowband bed B(x) (close to z/H = 0 at the divide) and surface S(x) (about z/H = 1 at the divide) approximate the north flank of Siple Dome. The divide thickness, H, is approximately 1000 m, and T is approximately 10 000 years. Note the arches under the divide (x/H = 0).

The horizontal-velocity shape function û(x, ž) is derived from a dynamic flow model. Except near divides, horizontal-velocity profiles at different locations tend to have similar shapes; as a result, û(x, ž) is a function primarily of ž = (z – B)/(S – B), the normalized height above the bed. As Equation (1) shows, ù(x, ž) is normalized such that its integral over the ice thickness is equal to unity. The depth-averaged velocity ū(x), which provides dimensions and scale for the horizontal velocity u(x, z), is derived from standard assumptions about the ice flux and the flowband thickness (Reeh, 1988).

In Jacobson and Waddington (2004), the horizontal-velocity shape function, û(x, ž), was derived with the shallow-ice approximation for isothermal ice (Hutter, 1983; Paterson, 1994, p. 262), which yields a shape function that is independent of x. Here, we use results from a finite-element flowband model (Nereson and others, 1998) to get polynomial approximations to û, yielding at the divide and at 10H downstream. At intermediate locations, we use a weighted average of these two profiles,

(2)

The weighting function, η(x), varies with distance from the divide, and is chosen to optimize the fit between horizontal-velocity profile shapes in this kinematic model and in the finite-element model (see, e.g., Nereson and Waddington, 2002). The flank and divide shape-function profiles are compared in Figure 2a. The flank shape function, has greater curvature at depth than the shape function from the shallow-ice approximation, because the finite-element model accounts for the polythermal nature of the ice, with warmer, softer ice near the bed. The divide shape function on the other hand, has greater curvature near the surface when compared to because the full strain-rate tensor is used when evaluating the ice viscosity.

Fig. 2. (a) Horizontal-velocity shape functions for divide û_d (solid line), flank û_f (dashed line) and isothermal shallow ice û _SIA (dot-dashed line). (b) Corresponding vertical-velocity shape functions, ŵ(ž) defined in Equation (4).

The vertical velocity is derived from the horizontal velocity by incompressibility (δ_x u + δ_zw = 0). Velocity is zero at the frozen bed, yielding

(3)

where Equations (1) and (2) provide u(x, z). In our model, we evaluate all terms that arise in Equation (3). However, we can gain some insight into the formation of Raymond bumps by examining the largest term contributing to w(x, z) when surface and bed slopes are small. Under those conditions, it is straightforward to show that w(x, z) from Equation (3) is dominated by a term where ḣ(x) is the rate of thickening, and where is defined as the vertical integral of

(4)

If the ice sheet is also steady (ḣ(x) = 0), then w(x, z) is dominated by –ḃ(x) ž(x, ž) (Reeh, 1988). Since ž(x, ž) expresses the leading contribution to the shape of the actual vertical velocity, and is the only contribution when bed and surface slopes are zero and û(x, zž) is invariant with x, w(x, ž) is often called a vertical-velocity shape function.

Figure 2b shows w(x, ž) at a divide and at 10H from a divide. The flank shape function ž _f(ž) is nearly linear, except close to the bed, while the divide shape function ŵ _d(Ž) is generally smaller at any given height ž, and more curved (almost quadratic) in the lower two-thirds of the ice. This means that w(x, z) decreases faster with depth under the divide than on the flank. Slower downward velocities directly under a divide produce the arch in the isochrones.

Particle paths are calculated numerically using the pair of differential equations for position, and ż = w(x, z), where over-dots indicate time derivatives.

For more details on this flowband model see Jacobson and Waddington (2004) and appendix A of Jacobson (2001).

3. Folding with Abrupt Divide Offsets

Consider a divide that stays at x = 5 H long enough to develop a mature Raymond bump, and then moves rapidly to x = 0H. The original divide-arch stratigraphy is then 5H away from the divide, and subject to a strain rate that differs significantly from that which formed the arch. How will the stratigraphy of this abandoned Raymond bump evolve, now that it is in a flank position and being carried further downstream?

To explore this situation, we first used our kinematic model to compute isochrones under a steady ice divide, similar to those in Figure 1, using accumulation and temperature patterns characteristic of Siple Dome. However, in this new case, the divide was located at x = 5H. Then, by shifting the steady ice-sheet surface profile, we abruptly shifted the divide by 5 ice thicknesses, to x = 0H. Using the steady-state flow field associated with this freshly offset divide, we then numerically tracked a large number of particles that started on each arched layer. Figure 3 shows the evolution of two arches. The upper layer does not overturn by 25H, though it could overturn further downstream. The leading limb of the lower arch becomes vertical by x = 12H. With this modeling approach, the angles of both fold limbs can be found at any location.

Fig. 3. Evolution of two isochrones in an abandoned divide arch after the divide moves abruptly from x/H 5 to x/H = 0. Ages at the time that the divide moves are 1:6T (upper) and 3T (lower). The layers are plotted subsequently at intervals of 0:5T.

Clearly, abandoned Raymond bumps can be seed disturbances for subsequent recumbent folds in ice cores, if the Raymond bump is sufficiently well developed, if a divide moves sufficiently far and sufficiently rapidly, and if the future ice-core site is sufficiently far downstream. Next, we describe how pre-cores can be used to explore the boundaries of the parameter space that allows the formation of these recumbent folds, without the need to track all possible arches through all possible divide migration histories.

4. Potential for Folding Assessed with Pre-Cores

Conventional isochrones are contour lines of age, showing the time since the ice accumulated at the surface of the ice sheet. These isochrones are also evolving material lines marking former surfaces of the ice sheet. Figure 1 shows these isochrones for our model. Note the arches under the divide.

As discussed in Jacobson and Waddington (2004), we can also calculate isochrones relative to other reference lines or surfaces, such as the vertical material line at a possible core site. During its evolution, this material line traces out curves that can be called core-relative isochrones, or pre-cores. Each pre-core shows where ice sampled in the core would have been found at a different earlier time. We can use our kinematic model to track ice in a hypothetical future core backward in time to find its pre-cores;pre-cores for a core at x = 20H are plotted in Figure 4. Figure 4 also shows corerelative isochrones downstream from the core site. These could be called post-cores; however, since post-cores are inherently less interesting to ice-core glaciologists, we will not consider them further.

Fig. 4. Pre-cores (dashed) and post-cores for a core at x/H = 20 (vertical line). The heavier lines mark intervals of 1 T. Selected particle paths are drawn as dotted lines. The three lower panels illustrate the change in shape of a disturbance at three points, (p), (q) and (r), along a particle path. In this example, the overturning limb of the fold at each point is aligned with the corresponding pre-core. θ is the pre-core slope angle (measured relative to the –x axis).

The lower set of panels in this figure illustrates how precores help us to assess the potential of observing recumbent folds in a core. For example, the fold at point (q) (the core) has a near-vertical righthand limb, which should stand out as a recumbent fold-in-the-making if observed in a core sample. At approximately 0.4T earlier, this particular structure would have been a symmetric open fold at point (p), with a righthand limb angle of 12°. This is the same slope as the pre-core at (p). At each point along a particle path, the pre-cores show the orientation that a layer must have in order to be vertical when it reaches this particular core site.

To evaluate whether an initially open fold in the stratigraphy could become a recumbent fold at a particular core location further downstream, we can superimpose that open fold on the pre-cores for that particular core location. Any portion of the open fold in an actual layer that is steeper than the pre-core at that point, but less than vertical (90°), will have rotated past vertical by the time it reaches the core location. Conversely, any actual layer that is less steep than the pre-core at the same location will not be overturned at the core site.

Now we revisit our example in Figure 3, in which a divide stayed at x = 5H long enough to develop a mature Raymond bump, and then moved rapidly to 0H. In Figure 5 we show the arched portion of the isochrones (the abandoned Raymond bump) left at x = 5H, together with pre-cores in the vicinity, taken from Figure 4. The pre-cores are based on the new flow pattern when x _div = 0H. The region where the isochrones in the arch are steeper than the pre-cores is outlined in gray. Any portion of an isochrone inside this region will rotate through vertical as it moves downstream to the core site at 20H.

Fig. 5. Offset (5H) surface-relative divide-arch isochrones (solid lines) (repeated from Fig. 1) overlaid on pre-cores (dashed lines) for a core at 20H. This is a portion of the flowband shown in Figure 4. The heavy line near z/H = 0 is the bed. The gray band outlines the area in which the isochrones are steeper than the pre-cores, indicating that these layers will overturn before they reach the core.

We can also confirm the effectiveness of our pre-cores predictions. Figure 6 again shows the evolution of the two arched layers in Figure 3. We have now added the pre-cores for the core at x = 20H, together with a heavy gray line marking the region in which layer slopes exceed pre-core slopes, i.e. in which layers are predicted to overturn before reaching the core site. The upper layer falls outside the region outlined in gray. As predicted, this fold limb does not overturn by 20H. The deeper layer falls within the overturning zone predicted by pre-cores. Our layer-tracking method confirms the simple pre-cores prediction; this arch forms a recumbent fold in a core at x = 20H. We can also see at a glance exactly which other parts of this abandoned Raymond bump would be overturned in a future core at x = 20H.

Fig. 6. Evolution of two divide-arch isochrones following abrupt divide offset from x = 5H to x = 0H (as in Fig. 3) is compared with pre-cores folding prediction. The dotted lines show pre-cores for a core at 20, and thick gray line outlines zone in which overturning was predicted by pre-core analysis. Layer ages at the time of the divide move are 1,6T (upper) and 3T (lower). The layers are plotted subsequently at 0.5T intervals. The 1.6T isochrone is above the zone of potential folding for a core at x/H = 20, and does not show signs of folding until well past 20, while the 3T isochrone is overturning by x/H = 12, confirming our pre-core analysis.

The distance that the divide moves affects the nature of the shear flow in which an abandoned Raymond bump finds itself. In Figure 7, we show the overturning regions for divide offsets of 1, 2 and 8H in addition to 5H. In all cases, the core location is at x = 20H, i.e. 20H from the final divide location. As before, a gray band outlines the region in which the isochrones are steeper than the pre-cores and therefore would be overturned later in the core. For the smaller divide moves, only deeper, older ice can be overturned. For the larger divide moves, shallower isochrones can be overturned, and larger volumes of ice can be affected. For offsets greater than a few H, in which the arch is left completely in the flank-flow regime following the move, the affected volume does not significantly increase in size.

Fig. 7. Thick gray curves outline volumes in which layers will be overturned in a core at x/H = 20, for divide offsets of x/H = 1,2, 5 and 8. Dashed curves show contours of pre-core slope in degrees. Thin gray and black curves show contours of negative and positive layer slopes (–2°, –4°, –6°, 2°, 4°, 6°), for arches left at x/H = 5 after a divide offset of 5 km (see Fig. 5). Thick gray curves are constructed by joining locations (open circles) where slope of precores is equal to slope of layers forming the divide arch. This procedure selects layers that will have a slope of 90° (i.e. will be vertical) at core site.

Figure 7 also illustrates the construction procedure whereby we identify these gray regions of overturning for the core at x = 20H. The most straightforward way to locate points at which two sets of curves (pre-cores and actual layers) are locally parallel is to find locations at which the slopes of both sets of curves are equal. Dashed curves show contours of the slope angle of the pre-cores in Figure 5 for a core at x = 20H. These pre-core angle contours are roughly horizontal in the regions in which overturning can occur, with the gentlest angles at the bottom (θ = 0°). The thin solid curves show contours of slope angles of the isochrones forming the arch in Figure 5. The roughly elliptical negative and positive isochrone slope contours highlight the left and right shoulders of the abandoned divide arch that originally formed at x = 5H. The region in which the arch will be overturned when the ice reaches x = 20H is outlined by the points of intersection, where these two sets of contours take equal values. For example, the points of intersection of the 2°, 4° and 6° contours are marked with open circles in Figure 7. At these points, the pre-cores are parallel to the isochrones.

Finally, similar maximum slopes of approximately 7° are found over a range of arched isochrones in Figure 5, and the shallowest part of each of the four overturning regions in Figure 7 coincides with pre-cores having slopes of 6–7º (see dashed curves). Because the overturning boundaries occur where the isochrone slopes are equal to the pre-core slopes, we see that the folding region in this example always appears to include the steepest parts of some arches, regardless of the divide offset. However, when the divide offset is less than a few ice thicknesses, folding is restricted to the steep sections of the deeper layers.

The shape and extent of these overturnable regions can be sensitive to the ice-velocity shape-function profiles used in the flow model, and to the choice of a core site where the layers will be sampled. Nevertheless, the existence of zones in which layers can be overturned is a robust feature of abrupt divide offsets.

5. Folding with Gradual Divide Movement

Although the pre-cores technique can give a straightforward and quantitative assessment of potential folding, it is predicated on essentially instantaneous divide motion. If the divide movement is slow enough, the stratigraphy may adjust to the new flow pattern without overturning. To explore this effect, we now assess how folding possibilities change when an ice divide moves from its initial location to its final location at a range of speeds. For these tests, we track layers with our kinematic flow model. Whereas Nereson and others (1998) and Nereson and Waddington (2002) found the shape of layers with steady divide motion, here we will include the transients associated with (i) starting from a stationary divide at time T ₁, and (ii) stopping the divide migration at a later time, T ₂.

To look at the effect of divide migration speed on the potential for folding, we make the divide position, x _div(t), a function of time, and redefine the surface profile, S(x – x _div), and the flux, Q(x – x _div), to be functions of the distance from the divide. Then we can calculate particle paths for an ice sheet with changing geometry. The flux vanishes at the divide, and the surface height at the divide is a maximum, with zero gradient. The blending function, η(x – x _div) in Equation (2), for the divide and flank velocity profiles is also a function of this distance from the moving divide.

To simplify calculations while moving the divide, we assume that the bed is flat (B(x) = 0), that the accumulation rate (ḃ) is constant and uniform, and that the surface assumes the simple analytical profile derived from the shallow-ice approximation (see, e.g., Jacobson and Waddington, 2004). Such a profile is an approximation even when the divide is fixed, and more so when divide movement introduces some degree of asymmetry. However, as argued in Nereson and Waddington (2002), slopes around the divide are small, and the asymmetry associated with the dynamical processes causing the divide to migrate does not have a significant impact on the shape of arches in the isochrones. Their kinematic model of divide migration incorporated an even simpler flat-surface profile.

Our divide movement function has the form

(5)

The divide arch isochrones can be generated and folded in a seamless manner by starting a set of particles at the surface between x = 4.75H and 5.25H at time t = 0. At a later time their subsequent positions define an isochrone. From t = 0 to t = T ₁, the layer develops the standard steady Raymond bump centered at x _div = 5H. At t = T ₁ the divide position starts migrating at speed V _div, and at T ₂ it stops. We continue to track the particles for t > T ₂.

Figure 8 shows the evolution of this layer for four different migration speeds, V _div. In all cases, the divide starts to move at T ₁ = 3T. The shape of the divide arch when the move ends (T ₂) is highlighted. The current divide position, x _div(t), is marked on each isochrone and under it on the bed with an open diamond symbol. In Figure 8a, the divide moves very fast (the T ₁and T ₂ isochrones are indistinguishable), and the layer overturns at x ≈ 12H, similar to the 3T isochrone in Figure 6.

Fig. 8. Folding of an isochrone when a divide moves to a new position at various rates. The divide is at x/H = 5 from t = 0T to 3T to build an arch. It then moves at speed V _div to x/H = 0. The isochrone at the end of this move (∆t elapsed time) is highlighted. Isochrones are plotted at 3.0, 3.5, 4.0, 4.5, 5.0, 6.0 and 8.07) V _div is expressed in units of H/T = ḃ = 0.1 ma^–1. The divide position is marked with a on the isochrones and on the x axis. The dotted lines are paths for three of the initial points (the two ends and the center). The bed at z/H = 0 is flat.

In Figure 8b, the divide migration speed is reduced by a factor of 1/50, yet the overturning edge of the former Raymond bump becomes vertical only a few ice thicknesses further downstream (≈14H). At this point it is also 0.5T older and slightly deeper than the overturning isochrone in Figure 8a. This trend continues in Figure 8c. There is, in effect, a trade-off between shearing the arch sooner (with a high V _div) and moving it deeper before shearing (with a lower V _div). This trade-off also means that it is difficult to evaluate the effect of the migration speed without also varying the depth and age of the particular layer that we consider in the divide arch.

In Figure 8d, the divide motion is slow enough that the sample layer does not overturn until close to 20H. The moving divide also leaves a raised ‘plateau’ in the layering behind it. After the divide stops at 0H (at 5T), it starts to produce a new divide arch (see the raised segment on the lowest isochrone, of age 8T).

In these examples, the divide moves only to 0H. If it continued to move in the same direction (to the left), the former divide arch would be sheared even faster. The same qualitative picture would hold, but the arches would be found even further downstream from the current divide and they could be even more overturned due to the larger bed-parallel shear strain. Conversely if the divide later moves back (to the right), the shear experienced by the arch would be less. But unless the divide returns significantly beyond its original position, the arch will remain on the flank and will continue to overturn.

6. Where the Folds might be Found

It is evident from Figure 6 that arches in different isochrones are displaced relative to each other as well as folded by the shearing ice flow. The deepest isochrones experience substantial shear, but little downstream motion because the bed is frozen. By the fourth snapshot (2T after the move), the upper layer (age 1.6T at the time of the move) is approximately 6H further downstream relative to the lower layer (age of 3T at the time of the move). Figure 9a shows the evolution of a column of arched isochrones. The divide was stationary at x = 5H for at least 9T (the age of the deepest isochrone at T ₁), i.e. all these layers had previously known only steady-state conditions. The subsequent divide motion followed the same pattern as in Figure 8b, with the move from 5H to 0H occurring in ∆t = 0.5T. The initial snapshot of the arches was taken at T ₁, the onset of divide motion, and the second snapshot 2T later.

Fig. 9. In each panel, a set of arched isochrones (on the left) ranged in age from 1.5T to 9T at t = T ₁. As in Figure 6, the divide then moved from 5H to 0H in time ∆t = 0.5T, then stopped. The same layers are shown (on the right) at time t = T ₁ + 2T. (a) The initial isochrones (left) formed under a steady-state divide. The deep arches are well developed, and after 2T their overturned limbs span 1H. However, if a divide is subsequently highly mobile, the steady-state initial configuration may be improbable. In (b), the divide was initially stationary at x/H = 5 for only 2T. Layers that were younger than 2T at t = T ₁ had experienced only steady-state conditions, and so show identical deformation histories to corresponding layers in (a). However, arches in the deeper layers were not as well developed as in (a), and as a result their subsequent recumbent folds, while still present, are less pronounced.

Because this arching and shearing affects a wide span of stratigraphic layers, the folds are eventually found over a substantial distance downstream from the originating divide. In Figure 9a, a core drilled anywhere from 4H to 17H downstream from the divide at 5H would sample a folded layer at some depth. Thus if a moving divide were to produce folding, quite a number of potential core sites could be affected by folds at later times, due to this mechanism alone. However, the overturned layers would be restricted to different depth intervals at different core sites at any particular time.

We should, however, keep in mind our assumptions about the lateral strain of the flowing ice. Our model with plane strain assumes that there is no variation orthogonal to the modeled flow. In three dimensions, such a divide arch has the form of a ridge with a large lateral extent. On the other hand, the arches under an axi-symmetric dome occur in a restricted circular region a few ice thicknesses in diameter (Hvidberg, 1996). When the divide moves away, leaving the arches in a flank position, this region spreads out laterally as well as downstream, reducing the tendency for recumbent folding (Waddington and others, 2001). We have not modeled folding under such a dome. However, we expect folds to be less common than in the corresponding ridge case. Sampling folds in the dome case would also be more difficult, since the abandoned arch would be confined to a narrow region downstream (in the new flow pattern) from the original divide position.

7. Role of Stationary Times

Of course, if the divide has moved throughout the past 2T, it may be unrealistic to assume that it had previously been in a steady state, or even stationary for at least 9T. Realistically, we should also consider ice divides that move at least as often as they are stationary. For example, Figure 9b explores the case where the divide was stationary at 5H for only 2T. Isochrones that were older than 2T at the first snapshot, i.e. at the time that divide motion began, spent the first part of their existence in flank flow, and therefore had less time to develop into arches.

To produce these deeper isochrones we used a hybrid velocity field. For the flux and surface-elevation calculations the divide stayed at 5H, but the velocity-profile blending function, η(x) (Equation (2)), used x _div = 0H during the flank-flow period. This made it easy to select the starting points for particle paths that would be under the 5H divide at the start of the 2T stationary period.

Before the divide starts to move, the arch amplitudes of the deeper layers in Figure 9b are reduced compared to their steady-state counterparts in Figure 9a. Because the slopes of the deeper arches are gentler in Figure 9b, the severity of folding is reduced in the second snapshot. However, this reduction is moderated by the fact that these lower-amplitude arches are also, on average, deeper in the ice and subject to greater shear than layers of comparable ages in Figure 9a.

The length of time that the divide is stationary at x = 5H prior to moving clearly modifies the development and subsequent folding of arches in the deeper layers. We have compared arch development for a range of divide-stationarity times between 3 T and 1 T. As a function of depth, the arch angles do not vary strongly until the stationary time drops to 1 T. On the other hand, as a function of isochrone age, the arch angles vary significantly for all these times, because the depth of the arch of a particular age also varies with divide residence time, due to the differences in downward-velocity pattern at a flank as compared to a divide. As we found when varying the divide movement speed, it is difficult to change just one variable in the complex of factors that affect the folding potential of an offset divide arch.

8. Folding Prospects at Various Domes

The discrepancy between the Greenland Icecore Project (GRIP) and Greenland Ice Sheet Project 2 (GISP2) ice cores brought the issue of possible folding in ice sheets to prominence (Grootes and others, 1993; Alley and others, 1995). However, this divide may not be a good candidate for the folding mechanism that we analyze in this paper, if there has been too much divide movement to develop the initial arch. Radar profiles of the current summit show little sign of arches, raising questions about its long-term stability (Hempel and Thyssen, 1992; Hindmarsh, 1996; Hvidberg and others, 1997). The position of this divide is sensitive to the dynamics of the whole ice sheet (Anandakrishan and others, 1994; Marshall and Cuffey, 2000). Modeling of the whole Greenland ice sheet over the last glacial cycle has suggested a bimodal pattern in the divide location, with a separation of about 70 km between dominant glacial and interglacial positions (Cuffey and Marshall, 2000; Marshall and Cuffey, 2000). Movement within their 20 km model grid size is less well constrained, but the glacial position appears to be more stable than the current interglacial one. If the glacial divide stayed within 6 km (2H) for intervals of 20 kyr then there could now be folded Raymond bumps in the ice sheet, but if the glacial divide was always on the move, then such folds are unlikely.

The clearest examples of Raymond bumps in radar profiles are seen at the West Antarctic ridges such as Siple Dome (Nereson and Raymond, 1996; Nereson and others, 1998), Roosevelt Island (Conway and others, 1999) and Fletcher Promontory (Vaughan and others, 1999). These are highly elongated domes, with ridge-like flow (plane strain) near their summits. Their boundaries are strongly controlled by bed topography and the surrounding ice streams and ice shelves.

Fletcher Promontory has a well-developed divide arch that is attributed to divide rheology. The tilt of the arch crest suggests a divide movement rate of 2ḃ (Vaughan and others, 1999). The Raymond bump at Roosevelt Island is (perhaps) the best-developed arch observed to date; however, this arch is probably accentuated by the thinning of a much thicker ice sheet over the past several thousand years (Conway and others, 1999). The divide position is constrained by a bedrock platform and the surrounding ice shelf, so there is little likelihood of significant divide movement in the past or future.

The arches at Siple Dome are not as prominent, possibly because the current divide movement rate is on the order of 5ḃ (Nereson and Raymond, 1996; Nereson and others, 1998). However, this dome is bounded by active and relict ice streams (Scambos and Nereson, 1995), so it is possible that this divide location has experienced the kind of stop– start motion that our model requires.

9. Conclusion

Arches in stratigraphic layers under a divide form a potential source of recumbent folds on the flanks of an ice dome. The arches must go through three stages: initial development, divide offset and subsequent folding. The divide must be fixed long enough for the arches to develop, with times on the order of several T = H = ḃ. Then, there must be a period of relatively rapid divide migration, followed by a period in which the abandoned Raymond bump remains on the flank. The divide migration puts the arch into a new flow regime, in which it is no longer an equilibrium feature. During the subsequent period of divide stability, the arches can overturn. The folds are found only downstream from the original stable-divide position, i.e. on the flank in the new geometry. Such a fold would be found under a divide only if the divide subsequently moved back past its original position. Divide movement is a possible mechanism for producing folds in ice-sheet stratigraphy, but it is not a promising explanation for folds at or near current divides.

Both tools that we introduced are useful for assessing recumbent folding of abandoned Raymond bumps. When a divide migration episode is rapid, pre-cores offer a simple and effective tool to identify all layer segments in a Raymond bump that would be overturned in a future core at any specified downstream site by the post-migration steady flow field, without the need to track all those layers through time. However, tracking particular layers with a transient flow model shows how folds evolve and where recumbent folds first develop. The layer-tracking approach also allows us to identify overturning layers even when a divide migrates relatively slowly, such that the layers undergo significant displacement and strain before the post-migration flow pattern is established.