Potential for stratigraphic folding near ice-sheet centers

Abstract

Lack of agreement between the deep portions of the Greenland Icecore Project (GRIP) and Greenland Ice Sheet Project II (GISP2) ice cores from central Greenland suggests that folds may disrupt annual layering, even near ice divides. We use a simple kinematic flow model to delineate regions where slope disturbances (“wrinkles”) introduced into the layering could overturn into recumbent folds, and where they would flatten, leaving the stratigraphic record intact. Wrinkles are likely to originate from flow disturbances caused internally by inhomogeneities and anisotropy in the ice rheological properties, rather than from residual surface structures (sastrugi), or from open folds associated with transient flow over bed topography. If wrinkles are preferentially created in anisotropic ice under divides, where the resolved shear stress in the easy-glide direction can be weak and variable, then the deep intact climate record at Dye 3 may result from its greater distance from the divide. Alternatively, the larger simple shear at Dye 3 may rapidly overturn wrinkles, so that they are not recognizable as folds. The ice-core record from Siple Dome may be intact over a greater fraction of its depth compared to the central Greenland records if its flat bedrock precludes fluctuations in the stress orientation near the divide.

Introduction

While folds, kink bands and shear planes are commonly observed in basal ice near the margins of ice sheets where bed-parallel shear is strong (e.g. Hudleston, 1976), most ice- sheet models predict that stratigraphic layers near ice divides in polar ice sheets should thin progressively over time, while maintaining their correct stratigraphic order. This is largely a consequence of the large scale and low spatial resolution of most ice-sheet models, which are not designed to address questions about local stratigraphic stability. Ice-sheet flow may appear to be steady when velocities are averaged over large scales, comparable to the ice thickness; however, simultaneous with such a large-scale “steady” state, layers in an ice sheet may still undergo distortion by small-scale transient- flow disturbances that average to zero mean on the large scale. Here, we take an initial look at conditions and processes that could enable one layer-disrupting process, folding, to operate on spatial scales that are small relative to the ice thickness. It may be surprisingly easy to develop folds in a polar ice sheet, even near an ice divide. In real ice sheets that have more complicated three-dimensional flow and bedrock effects, crystal anisotropy and generally unsteady flow at all scales, the potential for disturbance of the ice-core stratigraphic paleoclimate record is even greater, and additional processes such as shear band formation may also be active.

Inhomogeneities in rheological properties of the ice can alter stress patterns from those found in homogeneous ice with the same geometry, and the modified stresses can cause transient localized secondary flows that distort stratigraphic layering (e.g. Johnson and Fletcher, 1994). This is known as active deformation. For example, boudins and buckle folds can be caused, or amplified, by rheological contrasts between layers. Sub-horizontal layers of contrasting rheological properties exist in the Greenland ice sheet (Dahl-Jensen and others, 1997). Buckle folding of layers is unlikely near an ice divide where the ice flow is horizontally divergent; divide conditions are more conducive for boudinage (Staffelbach and others, 1988; Cunningham and Waddington, 1990). Almost all of the small folds observed in the Greenland Ice Sheet Project II (GISP2) core appear to have the same sense, suggesting that they are not simply drag folds on the edges of boudins. (The sense of shear appears consistent from meter to meter; however, the core was not oriented absolutely during collection and a few sections cannot be oriented relative to others (Alley and others, 1997).)

Alley and others (1997) showed that the observed folds at the scale of the ice-core diameter are associated with stripes of crystals with anomalously oriented c axes. These stripes appear to initiate in sub-vertical orientations, then rotate toward horizontal. This suggests that crystal c-axis fabric, through its strong influence on rheological properties, plays a role in the initiation of some folds.

In the absence of rheological inhomogeneity, stratigraphic layers are still deformed by the ongoing large-scale steady ice-sheet flow, in which they act as passive markers without influencing the stress and therefore the flow in their vicinity; this is known as passive deformation. For some strain-rate fields and layer orientations, passive deformation can convert subhorizontal layers into overturned folds (e.g. Hudleston, 1976). Layers that spend their entire existence in a steady-state flow field will not be passively folded; however, transient nonsteady flow at small scales within the ice sheet (e.g. flow disturbances induced by local inhomogeneity or anisotropy) could move a portion of a stratigraphic layer into an orientation that differs from the steady-state orientation of a layer at that position. The disturbed portion of the layer may then be subjected to subsequent passive folding, even if steady-state flow should subsequently be re-established at the small scale in its neighborhood. The creation of folds therefore depends on two factors (Waddington and others, 1995). First, there must be a transient source of disturbances or wrinkles in the otherwise steady-state stratigraphy. Second, the ice-sheet flow field in the vicinity of a wrinkle must be able to distort that wrinkle into an overturned fold.

Passive ice flow can have two competing effects on preexisting wrinkles. Simple shear parallel to the layers tends to overturn disturbances, generating folds (Fig. 1a). However, pure shear tends to stretch out and remove wrinkles (Fig. 1b). Disturbances to “perfect” steady-state layering are certain to exist or to be created by as-yet unspecified means in real and imperfect ice sheets, even near their centers; we focus here on deriving stability fields for those wrinkles. We identify regions in which the wrinkles tend to overturn into recumbent folds, and regions in which the wrinkles tend to be stretched and flattened. We also estimate upper limits to (a) the time that it takes for each unstable wrinkle to fold, and (b) the distance that the ice travels during the folding process.

Fig. 1. Effect of ice flow on disturbances in stratigraphic layering. (a) Simple shear overturns wrinkles, while (b) pure shear flattens wrinkles.

In this simple initial approach to near-divide folding, we use analytical models to approximate folding and flow. Numerical flow models that accurately calculate the finite strains during deformation (Jacobson, 2001) will ultimately give a more detailed picture of folding in more geometrically realistic ice sheets with more complex flow patterns; however, our analytical approach allows us to more easily identify key processes and parameters.

An Overturning Criterion

We take x and y as horizontal coordinates, with x in the direction of flow and z vertical. The corresponding velocity components are (u, v, w); we address motions (u, w) in the x, z plane.

Rotation of disturbed layers

In order to understand the behavior of wrinkles, we first express how a layer segment changes its slope with time. In Figure 2, AB is a segment of an isochrone. In a reference frame moving with the ice at point A, we can express the velocity of point B relative to A in first-order Taylor series around point A. This expansion is exact if the velocity gradients are homogeneous over this small region. We can also express the evolution of this segment (δx, δz) of the isochrone by

(1)

Using the quotient rule and Equations (1), the slope m = δz/δx of a passive marker evolves according to

(2)

where an over-dot indicates a material time derivative. When ice flows toward positive x, the first term rotates non-horizontal layer segments clockwise by horizontal simple shear, the second term flattens non-horizontal segments by longitudinal extension and vertical compression, and the third term rotates non-vertical layer segments.

Fig. 2. Velocity of point B relative to point A, fora layer segment of slope m = tan (θ = δz/δx. The leading edge of an upward wrinkle is (momentarily) stable if the relative velocity (heavy vector), found by summing four contributions in Equation (1), also has slope m. Grey arrows indicate magnitude and sign of vector components.

At every point (x, z) in a steady ice sheet, there is an isochrone whose slope can be represented by m ₀(x, z), which is independent of time. However, from the perspective of a lump of ice embedded in a particular layer, the slope of that material layer can change with time as the ice moves, because the slope must always match the spatially non-uniform (but temporally steady) slope m ₀(x, z) at each point along a steady streamline. Figure 3a and b show a net mass-balance pattern and a corresponding steady-state ice- sheet profile, resembling a typical transect across the Greenland ice sheet. However, the horizontal motion is modeled entirely by basal sliding, which is calculated from the local ice thickness and surface slope.

Fig. 3. (a) Mass-balance pattern that produces steady-state ice sheet (b) which moves entirely by sliding (i.e. no recumbent folds possible). Layer segments (thick black bars) rotate to match steady isochrones (dashed grey lines) as they move along particle paths ( solid lines). When the model includes steep gradients in accumulation rate (c) orbed topography (d), steady- state layer segments experience even larger rotations as they move along streamlines. The grey bars show the history of a disturbed-layer segment. Note the different horizontal scales in (b), (c) and (d).

The thickness H _* and the accumulation rate at the divide have been used as characteristic values for length and velocity to non-dimensionalize Figure 3. The corresponding characteristic time is , and non-dimensional variables will be indicated by tildes. For both central Greenland, where H _* ∼3000 m and (Bolzan and Strobel, 1994), and Siple Dome, Antarctica, where H _* ∼1000 m and (Nereson and others, 1996), the characteristic deformation time-scale is τ _* ∼10⁴ years.

Gradients in ice thickness and mass balance must exist between divide and terminus of even the simplest steady ice sheet. These gradients can change the slope m ₀ of a layer (indicated by thick black bars) as it moves along a streamline. In Figure 3c, accumulation is increased substantially in the narrow zone shown in the inset. This spatially non-uniform accumulation also causes stronger, more localized rotation of layers as they move through the steady-state layer pattern which is now more structured. In Figure 3d, a particular layer segment (indicated by black bars) must rotate substantially to continually match the steady layer pattern (dashed lines) as it moves over a bedrock obstacle. Since there is no horizontal shear in this model, no folding is expected. In an ice sheet with horizontal shear, a disturbed layer that is steepening may not be an incipient fold, if the steady isochrones also steepen along the streamline, while a layer that is flattening could actually be in the process of overturning into a recumbent fold, if a steady-state layer moving along the same path would flatten even faster. For this reason, we define the momentarily stable, or critical, wrinkle at any position as the disturbed-layer segment that is neither flattening nor steepening relative to a steady-state layer following the same streamline.

Letting θ = tan⁻¹ (−m) be the angle of a layer segment relative to horizontal (Fig. 2), application of the chain rule shows that

(3)

To find the orientations m _crit of layer segments (if any) that rotate at the same rate as the steady-state layers at the same position, first eliminate the rate of change of slope W from Equations (2) and (3) to get

(4)

Forming the corresponding equation for the steady isochrone slope m ₀, and subtracting it from Equation (4) leaves

(5)

which has a trivial solution, m _crit = m ₀, and a single nontrivial solution

(6)

showing that there can be at most one other angle at which a disturbed layer would rotate at the same rate as the steady isochrone. Since has only two zero crossings for all m, it must take the opposite sign when the slope m of a disturbed layer lies between m ₀ and m _crit, compared to outside this interval. A vertical layer, which we expect to lie outside this interval, experiences shear that causes it to overturn, so any layer with a slope between m ₀ and m _crit must rotate back toward the steady-state layer orientation m ₀.

The steady isochrone slopes m ₀(x, z) are small near divides. For example, in the Greenland ice sheet near the summit, |m ₀| ≤ 0.1 (Jacobel and Hodge, 1995), and by Equations (2) and (3). Expressing the difference in slope between a disturbed layer and the local steady isochrones by

(7)

the difference in slope between the critical wrinkles m _crit(x, z) and the steady isochrone m ₀(x, z) is

(8)

We use Δm _crit(x, z) to compare wrinkles at different places in the ice sheet.

Shear number

When we introduce a dimensionless shear number

(9)

the slope difference between critical wrinkles and steady isochrones is

(10)

measures the relative roles of horizontal simple shear and pure shear, which are the two dominant deformation patterns in ice sheets. is similar to other deformation indices used in structural geology (e.g. kinematic vorticity number (Means and others, 1980) or Index of Simple Shear (Hudleston and Hooke, 1980)).

Limitations

The stability criterion (Equation (10)) should be used with care. Although these critical wrinkles momentarily define the boundary between steepening and flattening layer segments, they tend to be advected into regions of greater horizontal shear, where they too may fold. In this sense, they overestimate the regions where wrinkles are stable.

In addition, Equation (10) does not infallibly identify incipient folding in the presence of complicated strain inhomogeneity. Steady layers change orientation along streamlines in the absence of accumulation or bed irregularities; however, the vertical exaggeration in Figure 3b is close to 100: 1, and slopes are also exaggerated by this factor. The rotation also occurs slowly, since the time taken for ice to traverse a streamline is comparable to τ _* or longer. Both m ₀(x, z) and are very small, and will have negligible impact on estimated layer stability. Equation (10) should give adequate stability estimates for an ice sheet with a smooth bed and mass-balance distribution. The vertical exaggeration in Figure 3 c is closer to 10:1, suggesting that changes in m ₀(x, z) and along streamlines as a result of accumulation irregularities are likely to be noticeable; however, Equation (10) accounts for these changes.

There is no vertical exaggeration in Figure 3d, and changes in m ₀(x, z) and along streamlines as a result of flow over bedrock topography can be substantial. In this case, even measuring slopes relative to the steady isochrone pattern can be inadequate. Figure 3d shows the spatially varying steady-state isochrone slope pattern (black bars) encountered by a parcel of ice containing a disturbed layer (grey bars) as it follows a streamline up and over a bedrock bump. The angle between the steady isochrones and the disturbed layer closes as the ice ascends the bump, and opens again as the flow descends the lee side. However, this increasing angle (which is what our Δm _crit measures) is not an incipient fold; since the horizontal velocity in this example is entirely taken up by basal sliding, i.e. δ_zu = 0, folding of disturbed layers by horizontal simple shear is not possible. The vertical extension that widens the angle in the lee of a bump merely recovers vertical compression that the ice experienced as it climbed the stoss face. Longitudinal compression alone cannot overturn a layer, and when the bedrock is rough we cannot conclusively identify wrinkles that will fold just by accounting for rotations relative to the steady isochrone pattern. Rather than modify our Δm _crit criterion further to account for bedrock topography, it will be more practical in future to analyze the total strain along each particle path (Jacobson, 2001). Because of this limitation, we restrict our analysis to cases where bedrock relief is small.

Ice-Sheet Flow Model

Now we use Equation (10) to determine where wrinkles are likely to overturn or flatten in an ice sheet. To construct analytical models of wrinkle behavior, we must capture the most robust characteristics of the large-scale flow in real ice sheets.

A modeling strategy

We describe the large-scale ice-sheet flow near steady ice- sheet centers by drawing on previous kinematic analyses by Nye (1959) and Reeh (1988, 1989). We model ice deformation using Glen’s flow law (Glen, 1958), with a stress exponent of n = 3 (e.g. Paterson, 1994, p.91). This relationship assumes that polar ice is homogeneous and isotropic. Other, more recent constitutive relations for ice incorporate effects of crystal anisotropy (e.g. Azuma, 1994) and grain-size (e.g. Goldsby and Kohlstedt, 1997; Guffey and others, 2000). Small-scale inhomogeneities in fabric and grain-size may be the prime initiators of the secondary flows that generate new wrinkles on stratigraphic layers by active deformation. However, Glen’s flow law adequately reproduces the large- scale features of the steady-state flow field that we invoke to investigate passive folding of wrinkles that might be produced on that steady layer pattern by small-scale transient processes. Rather than solve momentum-conservation equations to find the flow pattern for a complicated ice- sheet geometry (e.g. Raymond, 1983; Schott and others, 1992; Hvidberg, 1996), we use a shape-function approach (e.g. Reeh, 1988, 1989) to calculate velocity. Shape functions representing the depth variation of the velocity components are usually derived from solutions of the momentum-con- servation equations for a relatively simple geometry. The magnitude of the steady velocity field is derived from mass conservation; the ice flux through a gate of width W(x) across a flow-band is equal to the integrated accumulation rate upstream in that flow-band.

A sandwich model

Since thickness and accumulation rates generally have low spatial gradients near ice-sheet centers, we use uniform thickness H(x, y) = H_* and uniform accumulation rate Because we describe the flow by kinematic shape functions, we do not need a surface slope to drive ice flow. Reeh (1989) called this a sandwich model.

Recalling that non-dimensional variables are indicated by tildes, we approximate the large-scale steady flow field near the center of polar ice sheets, using the well- known parallel-sided slab approximation for steady kinematic flow of isothermal ice frozen to a flat bed. The horizontal velocity takes the non-dimensional form (e.g. Raymond, 1983)

(11)

where n is the exponent in Glen’s flow law (Glen, 1958). A describes how horizontal strain rate is partitioned between longitudinal and transverse strain rates. For an ice ridge in plane strain, Λ = 1, while for the corresponding circular ice dome, Λ = 1/2 due to lateral spreading; this result was noted by Nye (1959) and by Reeh (1989). The vertical velocity field is

(12)

which is the same for ridges or domes. In this simple model, all isochrones in steady state are flat (m ₀ = 0). Substituting velocity gradients derived from Equations (11) and (12) into Equation (9) yields the shear number

(13)

which increases with distance from the ice divide, and with depth .

Stability boundaries

The contour , along which the critical wrinkle slope (from Equation (10)) is a constant, is readily found from Equation (13). Figure 4a shows these critical stability boundaries for a range of wrinkle slopes Δm _crit under a divide ridge (right) and under a circular dome (left). Above each contour (parameterized by a particular value of Δm _crit ), wrinkles with that Δm are flattening (Fig. 1b), while below the line, wrinkles with that Am are overturning (Fig. 1a). Under the ice divide at , any transient wrinkle in the stratigraphy tends to stretch and flatten. Away from the ice divide, however, wrinkles can be overturned, and wrinkles with progressively smaller slopes become unstable with increasing depth and distance from the divide; this is a consequence of the increasing dominance of bed-parallel simple shear strain rates over longitudinal stretching strain rates.

Fig. 4. (a) Stability fields for various wrinkle slope disturbances Δm _crit in a simple sandwich ice-sheet flow model. ( b, c ) Upper bounds to dimensionless overturn times in units of the characteristic time ,for steep wrinkles (Δm = 1.0) and for wrinkles with lower slope (Δm = 0), respectively. Grey trajectories place upper bounds on the displacements of these same wrinkles during the folding process. In all panels, right side shows results under a ridge divide (plane strain), while left side is result under a circular ice dome.

Overturn times

Suppose that at time , a wrinkle with slope is introduced at some point where the shear number is . If the wrinkle slope m _p (= Δm _p in the sandwich model) at P exceeds the critical slope , we approximate the time T required to turn this wrinkle into a recumbent fold by the time required for point B in Figure 2 to overtake point A. This is a reasonable criterion, because as the leading fold limb AB on an upward wrinkle reaches vertical, the folding process accelerates; the horizontal simple shear and the normal strain rate begin to work in tandem rather than in opposition. In Equation (2), the second term changes sign as m changes sign, but the first term does not. Therefore, wrinkles spend relatively little time with near-vertical edges; once they reach vertical, they quickly overturn into recumbent folds.

If the strain-rate field is homogeneous in the region traversed by the wrinkle during the time required to fold, then the Equations (1) (with ) are readily solved for

(14)

(15)

The time at which gives the overturning time

(16)

Since , the height of a wrinkle (Equation (15)) decreases even as its leading edge steepens. At the time when the fold overturns, the wrinkle height is reduced by the factor

(17)

We now estimate the time required to overturn a wrinkle with slope m _p in the sandwich model. Substituting velocity gradients derived from Equations (11) and (12) into Equation (16) leads to an equation for the contour line joining all points at which a wrinkle introduced with slope m _p would take time to overturn,

(18)

Figure 4b and c show these contours of non-dimensional overturning time for wrinkles with slopes of Δm = 1.0 and Δm = 0.1, respectively. Within 10H _* of an ice divide, recumbent folds can form from wrinkles with slope Δm = 1 in times as short as 0.03τ _*, and even low-angle wrinkles (Δm = 0.1) can overturn in as little as 0.25τ _*; these times correspond to a few hundred years and a few thousand years, respectively, at both Siple Dome and the Greenland summit.

Wrinkles of a given initial slope overturn more slowly under a dome than under a ridge because, in a dome, only half of the simple shear in the horizontal plane occurs in the direction of ice flow.

Because ice moves from P to positions where the shear number tends to exceed , folding times given by Equation (16) tend to overestimate the time required to fold. The estimates are more reliable for shorter overturn times because the ice is less likely to move into a significantly different strain-rate regime.

Displacement during folding

An unstable wrinkle moves along its streamline during the time required to overturn. We can estimate this displacement for a wrinkle that had slope m _p at time at position where the velocity is . First, we define ; this is the depth below point P at which the vertical velocity would reach zero (i.e. an effective “bed”) if the vertical strain rate were actually homogeneous as we assume. We also define ; this is the distance upstream from point P at which the horizontal velocity would vanish (i.e. an effective “divide”) if the horizontal strain rate were actually homogeneous. Then,

(19)

The wrinkle advects toward the effective bed at an exponentially decreasing rate set by . It is also swept downstream at an exponentially growing rate (set by ) that scales with distance from the effective divide and height above the effective bed. The trajectories (grey lines) in Figure 4b and c show the estimated displacements (Equation (19)) during the dimensionless time required to fold a selection of wrinkles of initial slopes Δm = 1.0 and Δm = 0.1 under circular domes and ridges. Glearly, wrinkles of modest slope can be overturned close to their points of origin.

To derive Equation (19), we assumed that the velocity gradients were homogeneous throughout the region traversed during the folding process. This assumption becomes inappropriate when the folding process takes times that are comparable to τ _* or longer, since is a characteristic time to move a significant distance through the ice sheet. Where the homogeneous strain-rate assumption clearly breaks down, we show the trajectories by dotted lines. The assumption causes Equation (19) to overestimate the distance traveled prior to folding; wrinkles can fold even closer to their origin points than we show in Figure 4b and c. We could include the strain-rate gradient terms in this analysis. However, the complexity of an analytical solution could obscure further insights into characteristic scales and times, and the solution would still be approximate. To progress further, we recommend use of numerical models incorporating finite strains and ice velocities that can take account of site-specific bed and surface topography (Jacobson, 2001).

Sources of Wrinkles

Sastrugi and snowdrifts

We now consider processes that could produce wrinkles. Surface features such as sastrugi or drifts are unlikely to be folded near the ice divide, because their slopes are reduced both by firn compaction and by longitudinal straining in the upper layers of the ice sheet. For example, a snowdrift might have a slope of Δm ≈ 1.0 at the surface, in snow with a density of 300 kg m⁻³ . After vertical compaction into ice of density ≈900 kg m , the slope would be reduced to approximately Δm = 0.3. It would be further flattened by both vertical compression and longitudinal extension associated with ice flow. In steady state, the vertical velocity (Equation (12)) also gives the relative thickness of layers at each depth; at depth , a layer retains a fraction of its initial thickness at the surface. Gonservation of mass implies that it must stretch horizontally by the factor . Under an ice divide in plane strain, slopes at depth are reduced by the factor relative to their initial value at the surface. Under a dome, only half of the horizontal spreading is directed in the flow direction, so slopes at depth are reduced by . According to Equation (12), this factor is 0.14 for a ridge and 0.23 for a dome at 50% depth; the snowdrift (with slope Δm = 0.3 after compaction into ice) would have slopes Δm = 0.04 and 0.07, respectively. At 80% of full depth, the corresponding slopes are Δm = 0.002 and 0.007. A snowdrift subjected to this severe flattening would be unlikely to catch up with or cross a fold stability boundary (Fig. 4a) until it was near the bed, and remote from the ice divide.

Bedrock topography

Standing waves form in layers flowing over undulating bedrock (e.g. Fig. 3d). Their amplitude decays with height above the bed, with a scale length comparable to their wavelength (Whillans and Johnsen, 1983). If the direction of flow should alter, the existing standing waves might not match the steady layer pattern for the new flow regime. Layer “disturbances” Δm would result from differences between the old and the new steady isochrone patterns. The maximum slope disturbances that can be generated in this way are smaller than the slope of the bedrock itself, and can be found only within a wavelength of the bed. In central Greenland, bedrock undulations on scales of a few kilometers have slopes on the order of m = 0.1 or less (Hempel and Thyssen, 1993; Jacobel and Hodge, 1995). A change in flow direction associated with divide migration at the last deglacial transition (e.g. Marshall and Guffey, 2000) could have created long- wavelength slope disturbances that may affect the lowermost 10% of the ice sheet today. At Siple Dome, bedrock slopes are an order of magnitude smaller, and so are less likely to generate wrinkles during transient flow.

Bedrock features with steeper slopes probably have wavelengths too short to disturb layering as high as 0.1H_* above the bed, where disturbances occur in the central Greenland cores. Bedrock standing-wave patterns also cannot form the very short wavelength folds seen in the ice cores (Alley and others, 1997), or generate the initial steep disturbances required to cause folding relatively high above the bed near the divide.

Internal wrinkle generation

It is difficult for surface or bed processes to create adequate slope disturbances on layers at mid-depths in the ice. However, plate 2 of Alley and others (1997) shows small folds and fabric disturbances that were observed at mid-depths in the GISP2 ice core. Inclined layers that must form parts of disturbances with scale sizes that exceed the core diameter were observed at mid-depths in both Greenland Icecore Project (GRIP) and GISP2 cores (Alley and others, 1995). Most wrinkles at mid-depths are probably created in situ by inhomogeneous deformation associated with spatial variability in rheological properties (viscosity and anisotropy) on a variety of spatial and temporal scales.

Temperature strongly affects ice viscosity; between −10° and −50°C, viscosity changes by >2 orders of magnitude (Paterson, 1994, p. 97). However, temperature varies smoothly and slowly in polar ice sheets. While it affects deformation patterns at the H_* scale, significant temperature-induced flow variability does not occur at the 10⁻² to 10¹ m scales of the wrinkles that we consider.

Soluble impurities may soften polar ice, and variations in impurity content could create local hard and soft layers, possibly through their effect on crystal size (e.g. Cuffey and others, 2000). When a horizontal layer of a more viscous material is embedded between layers of less viscous material, and then subjected to horizontal extension, the more viscous layer can act as a stress guide. It can develop necks that evolve into boudin structures resembling link sausages, i.e. a regular pinch-and-swell wave on its boundaries. Smith (1977) described the growth of boudins by a perturbation analysis. Layers of contrasting viscosity exist in central Greenland (Dahl-Jensen and others, 1997), and existence of boudins near the Greenland ice divide was postulated by Staffelbach and others (1988) and by Cunningham and Waddington (1990). Small boudins were seen in the GISP2 ice core. When boudins form, their wave-like surfaces could provide wrinkles to be folded by subsequent passive deformation.

Rheological contrasts may occur in patches as well as in layers. As ice flow proceeds near an ice divide, a harder lump will be compressed more slowly in the vertical direction than the material to either side, and the stratigraphic layers above and below a harder lump will drape around it. Even when the large-scale steady flow is two-dimensional, the wrinkles generated around three-dimensional inhomogeneities are three-dimensional; they can be folded by subsequent flow in any direction.

The fabric stripes and associated small folds reported by Alley and others (1997) appear to result from local variations in c-axis orientation, which is probably the primary control on ice deformation at small scales. Ice crystals shear much more readily along their basal plane than along any other plane, and polar ice can develop a strong crystal-fabric anisotropy with depth (Alley, 1992) in which the c axes (normal to the basal planes) tend to cluster vertically. This makes the ice soft for horizontal shear and hard for vertical compression and horizontal extension. Vertical variations in the strength of this c-axis clustering (Gow and others, 1997; Thorsteinsson and others, 1997) and associated rheological properties (Dahl-Jensen and Gundestrup, 1987; Dahl-Jensen and others, 1997; Thorsteinsson and others, 1999) have been observed, particularly across large climatic-change boundaries. Under tension, boudins may also develop on layers with stronger anisotropic fabric. Lateral variations in the strength of c-axis clustering would make ice inhomogeneous for vertical compression; ongoing deformation could result in wrinkles in the stratigraphic layering at the corresponding scales. At the scale of a few ice crystals, random variations of crystal fabric could create lumps that are harder or softer to vertical compression than their surroundings. This may be a factor in the origin of the fabric stripes and small folds reported in the GISP2 core by Alley and others (1997). At larger scales, random variations in the number density or area of the fabric stripes themselves could create harder or softer zones. Under conditions where recrystallization becomes a factor (Alley, 1992), some patches of crystals may recrystallize before neighboring patches; this could lead to spatially inhomogeneous crystal fabric that would in turn make some regions harder or softer to vertical compression than neighboring regions. Calculating the response of these zones to the applied stress would require a full anisotropic flow model (e.g. Azuma, 1994; Thorsteinsson and others, 1999); nevertheless, effects on stratigraphy should qualitatively resemble effects in an isotropic viscous fluid containing lumps of differing viscosities.

Ice with the strong vertical c-axis fabric found under both ice-divide and flank locations can deform most easily in horizontal simple shear. Within approximately one ice thickness on either side of an ice divide, the dominant stresses are vertical compression and horizontal extension, and there is very little resolved shear stress in the direction of easy deformation. As a result, small changes in stress orientation can cause large changes in the deformation rate and pattern. Azuma and Goto-Azuma (1996) and Castelnau and others (1998) suggested that small stress variations associated with flow over sloping bedrock may subject ice in the divide region to stronger anisotropic-flow instabilities, compared to ice on the flanks. On the flanks the dominant stress is horizontal shear, which is oriented approximately in the direction of easy glide for the existing fabric. There, small fluctuations in the stress orientation exert only a small influence on the deformation rate. In our terminology, the divide region may generate more frequent or steeper wrinkles than flank regions.

Application to Ice-Core Sites

The Greenland summit

Comparison of the GRIP and GISP2 ice cores (Grootes and others, 1993; Taylor and others, 1993; Alley and others, 1995) showed that the stratigraphic records agreed down to approximately 90% of the ice depth, but below that level the stratigraphy was disturbed in one or both cores. The GISP2 ice core is situated approximately 10H _* from the ice divide, which is currently less than one ice thickness to the east of GRIP (Hvidberg and others, 1997). The flow pattern is intermediate between a plane-strain ridge and a circular dome. Figure 4a shows that at GISP2, wrinkles with slopes of Δm ≈ 0.5 are unstable over more than half the depth, and slopes as small as 0.05 are unstable in the lowermost 500 m of the core. At the ice divide near GRIP, this simple theory predicts that all wrinkles should be flattened and eliminated by strain; however, the ice divide has not always been near GRIP (e.g. Anandakrishnan and others, 1994; Marshall and Guffey, 2000) . If the divide was previously as little as 2H _* = 6 km from GRIP, wrinkles of slopes as gentle as Δm = 0.1 would have been unstable and could have overturned in the lowermost 10% of the ice sheet, where the stratigraphy now differs between the two cores. If the disturbances seen in the GRIP core (Alley and others, 1995) result from folding as we have described, then the GRIP site cannot always have been on or near the ice divide over the time-scales required (Fig. 4b and c) to overturn layer disturbances.

In addition, if flow near the divide is subjected to stronger anisotropic instability (Gastelnau and others, 1998) then the number density of wrinkles generated may be larger near the divide, making the divide region more vulnerable to fold disturbances than our simple model predicts. Unfortunately, at GISP2 and GRIP, the disturbed ice is nearly 10⁵ years old, and its origin location relative to an ice divide is unknown, so we cannot test whether the disturbed ice originated in a divide zone (Marshall and Guffey, 2000).

Our analysis of wrinkle folding depends only on the slope, and not on the amplitude of the wrinkle. However, in general, there may be a spatial-scale hierarchy for wrinkle generation (Alley and others, 1997). At the GRIP and GISP2 core sites, a few small folds such as those seen by Alley and others (1995) at the scale of 10⁻² to 10⁻¹ m appear to be created first, in increasing numbers over time (or depth). The errors that these small folds introduce into the paleoclimate stratigraphy occur at vertical scales too short to be readily detectable. Larger folds, if they occur, appear to arise only later, after larger total strains have been experienced, or after larger shear numbers have been encountered, or after enough small wrinkles, folds or stripes have been created such that variations in their number density create rheological inhomogeneity at increasingly larger scales. Where the growing length scale of vertical disruption exceeds the depth and temporal resolution of the cores, the GRIP and GISP2 cores lose coherence.

Dye 3

The Dye 3 ice core (Dansgaard and others, 1985) was recovered from a flank position approximately 15 ice thicknesses from the summit of the Greenland South Dome. The upper 95% of its length appeared to be stratigraphically intact at the temporal resolution at which comparisons with the GRIP and GISP2 cores are possible (Johnsen and others, 1992). A combination of factors could contribute to this greater fraction of intact record at Dye 3 than at GISP2 or GRIP.

Since the Dye 3 core site is farther (in units of H _*) from the divide than the GISP2 site, a larger fraction of the Dye 3 core contains ice originating on the flanks. Gastelnau and others (1998) suggested that flank ice may have been subjected to fewer anisotropic-flow instabilities than ice from a divide region. The South Dome ice divide may have existed at or near its current location, making Dye 3 a remote flank site, since late Wisconsinan time when the deepest undisturbed ice was deposited. An ice-flow reconstruction (Reeh and others, 1985) suggested that ice originating 1 km from the divide (in the divide zone) at approximately 16 ka BP is now 180 m above the bed at Dye 3. Existence of large-scale disturbances only in ice deeper than this is not inconsistent with creation of wrinkles predominantly in the ice-divide zone. Alternatively, a relatively high shear number S near the bed at Dye 3 in the Wisconsin ice (which is soft for bed- parallel shear) could also have created disruptive folds from wrinkles with relatively large vertical extent but low initial slope after that basal ice left the divide.

As ice moves farther from the ice divide, the shear number increases, and the time to overturn a wrinkle decreases (Fig. 4). In deep ice far from a divide, unstable wrinkles overturn rapidly. Once such a wrinkle becomes a recumbent fold with flattened limbs, it may no longer be recognizable as a strati- graphic defect. If wrinkles are created at a relatively uniform rate throughout the ice sheet, then as the shear number increases, the number of wrinkles still visibly unfolded (see, e.g., Alley and others, 1997) decreases due to increased efficiency of the folding mechanism.

If the Dye 3 core contains a relatively high number of fold defects that can no longer be recognized as defects, then the fact that the climate record is still intact at the level of temporal smoothing at which it can be matched to the GRIP record (Johnsen and others, 1992) implies that there is an upper limit to the vertical scale of folds at Dye 3, and this limit is less than the resolution of the paleoclimate record.

Siple Dome

A 1000 m ice core was recovered from Siple Dome in January 1998. The present accumulation rate is approximately 10 cm a ice equivalent (Nereson and others, 2000), and the characteristic time for flow is τ _* ≈ 10⁴ a. A 100 kyr climate record can be expected (Nereson and others, 1996). The core quality was relatively poor in the lower sections, and it is not yet clear whether there are visible disruptions or folds in the core. Preliminary inspection did not reveal strong evidence of disruptions or folds, but work is ongoing. Since the bed is flat, our analysis of fold stability is not restricted by bed topography. Based on radar internal layering, however, Nereson and others (1998) concluded that the ice divide has been migrating northward for at least several thousand years at ≈0.3 m a⁻¹. The ice-core site, which is currently 0.5H _* from the ice divide, was probably directly under the divide several thousand years ago. Although the West Antarctic ice sheet is thought to have undergone more dynamic change than the Greenland ice sheet during the last deglaciation, during the past 15 kyr, Siple Dome was probably never much thicker than it is today (unpublished model results based on the observed depth–age). Siple Dome is located between Ice Streams G and D, which have probably been persistent low-profile fast-flow features in West Antarctica; this suggests that the ice-divide position may have been more stable at Siple Dome through the last deglaciation than at Summit Greenland.

The stability boundaries suggested by Figure 4 are probably still approximately correct; however, the boundaries and overturn times could be asymmetric as a result of different flow conditions experienced by ice on opposite sides of the moving divide. An ultimately more satisfying approach would be to follow the time-dependent finite strain of layers as they move in the transient-flow pattern associated with the migrating divide.

Because the bedrock is so flat at Siple Dome, the interaction of variable stress orientation (dictated by bed slope) with anisotropic fabric (Gastelnau and others, 1998) should be less effective at producing wrinkles in the divide zone at Siple Dome in comparison with central Greenland. If so, the Siple Dome core should be disturbed less by folds than ice in central Greenland.

Conclusions

In our simple model of folding near ice-sheet centers, two scales of ice-sheet flow are important. We envision an ice sheet whose flow is steady when averaged over large volumes (on scales comparable to the ice thickness); however, the flow can be unsteady on small scales (e.g. 10⁻² to 10¹ m), at which tran- sient-flow disturbances (which average to zero on the large scale) can arise due to spatial variations in rheological properties. These transient-flow disturbances can create wrinkles on the otherwise steady-state stratigraphy. We approximate the large-scale steady flow using a simple steady-state kinematic ice-flow model based on shape functions. This large- scale steady flow either passively folds, or else flattens and removes, the small wrinkles.

Steeper wrinkles fold more readily than low-angle wrinkles. Wrinkles of a given slope tend to flatten if introduced at shallow depths or near an ice divide, but tend to fold at progressively faster rates at greater depths and greater distances from a divide. This suggests that, if wrinkles are generated uniformly throughout an ice sheet, then a steady ice divide is the best place for an ice core. However, if anisotropic-flow instability near divides generates steeper wrinkles, or a higher number density of wrinkles there, then there may be an optimal distance from an ice divide where the potential for folding disturbances in a particular depth or age section of an ice core can be minimized.

Although our model is simple, it predicts the existence of overturned folds near ice-sheet centers and ice-core sites. Real ice sheets, which have additional complications such as bedrock topography, extensive variations in ice fabric and texture, and transient flow at all spatial scales, can only be more vulnerable to folding, and to other layer-disrupting processes such as shear band formation. Further progress in understanding the folding potential in polar ice sheets will be achieved by investigating the active and potentially anisotropic deformation processes that initiate wrinkles and shears (e.g. Thorsteinsson and Waddington, in press), and the large finite strains that can overturn layers in inhomo- geneous and time-varying deformation fields.