Folding in strongly anisotropic layers near ice-sheet centers

Abstract

Conditions for passive folding near ice-sheet centers are derived treating the ice as an anisotropic viscous medium. Vertical uniaxial compression at a dome, or pure shear stress near a ridge divide, both tend to stretch and flatten folds, while horizontal simple shear deformation tends to overturn folds. Overturned folding is likely in a given vicinity in a steady-state flow field, if the initial slope of a layer disturbance exceeds the ratio of compressive/extensive deformation to shear deformation. Analytical equations for particle tracks in steady state allow us to model the evolution of layers with initial slope disturbances (``wrinkles’’). the effects of anisotropy are explored using an analytical solution for the strain rate as a function of a vertically symmetric c-axis orientation distribution, called a cone fabric. Stronger anisotropy (small cone angle) makes the material softer in horizontal shear, and facilitates folding, i.e. elements with smaller slopes can be overturned for the same stress. the relation between anisotropy and folding is complicated by the fact that, for a range of cone angles, the material is also softer in compression, which opposes folding. Simulating a layer with spatially variable tilt of cone symmetry axes, and accounting for fabric development, demonstrates that variations in the fabric cause localized flow variations that could create the initial perturbations.

1. Introduction

Our knowledge of past climate conditions comes largely from ice cores and other sedimentary records. Interpretation of these records relies on the assumption that the stratigraphic layering has remained intact. for marine and lacustrine sedimentary cores, bioturbation and tectonic displacements are major concerns. for ice-core records, layer stratigraphy is usually assumed to have remained intact except possibly very close to the bed. However, the loss of correlation of δ¹⁸O (Alley and others, 1995) and other climatic indicators below 2750 m depth (some 270 m above the bed) in the Greenland Icecore Project (GRIP) and Greenland Ice Sheet Project 2 (GISP2) ice-core records, and visual observations of small-scale folding in the bottom third of the GISP2 core, point to the possibility of layer disturbance at considerable heights above the bed, even near an ice divide (Alley and others, 1997). Layer disturbances due to varying flow over bed topography are generally expected to have wavelengths comparable to the bed topography and to be found within one wavelength of the bed. Therefore the small features seen in the Greenland cores are unlikely to be directly caused by bedrock topography. Folding of surface features such as sastrugi is unlikely, due to the longitudinal stretching of layers near ice-sheet centers.

Figure 1 shows a simple description of the deformation states near ice-sheet centers, with longitudinal stretching and horizontal shearing competing to flatten out or to overturn any ``wrinkles’’ in the layering. This is very different from the conditions near the bed or at the margins, where horizontal compression is common (Hudleston, 1976; Hambrey and Lawson, 2000).

Fig. 1 Evolution of layer disturbances in steady-state flow fields typical of near-divide environments. (a) Pure shear deformation tends to flatten disturbances, while (b) simple shear can overturn the disturbance to cause recumbent folding.

Layer disturbances such as wrinkles (undulations, or open folds in annual layers) can lead to overturned or recumbent folds that have the potential to severely compromise the chronological interpretation of ice cores. In an ice sheet that has an overturned fold, an ice core could conceivably go through the same layer three times. the origin of layer disturbances has not yet been fully resolved. Jacobsen (2001) shows that the special arches found in layers under steady ice divides (Raymond,1983) can act as a disturbance source for overturned folds if the divide then moves. More generally, transient flow may leave layers in disturbed orientations near bedrock bumps (Hudleston,1976). Almost all of the various potential surface sources (sastrugi, accumulation variations) of layer disturbances are unlikely sources of wrinkles near ice divides, due to the strong stretching of annual layers near ice-sheet centers (Waddington and others, 2001). Variations in rheology and anisotropy are the most likely sources of layer disturbances within an ice sheet (Dahl-Jensen and others, 1997). Several authors have pointed out that asymmetry of the fabric might cause flow inhomogeneities that could lead to layer disturbances (Azuma and Goto-Azuma,1996; Castelnau and others, 1998).

Here we derive a criterion for folding of layer disturbances injected in an otherwise steady-state flow field in vertical uniaxial compression or pure shear stress, combined with a simple shear stress state. Assuming homogeneous strain-rates equations for particle tracks allows us to model the evolution of layers. Waddington and others (2001) followed a similar approach to assess folding potential for isotropic ice. Here the effects of anisotropy on the folding criterion are explored using the anisotropic flow law formulated by Thorsteinsson (2000, 2001, 2002).

Finally, we consider the origin of disturbances to the layering. Fast changes in stress state, which might arise as ice flows over bumps or as a divide migrates, can leave the fabric in a different symmetry to that which it would have developed in the new stress state. We will refer to fabric that has symmetry different from the fabric that would have developed under steady conditions as ``unadjusted fabric’’. Smaller-scale heterogeneity, such as stripes (Alley and others, 1997), is also a potential source for layer disturbances at all depths. the deformation of a material with an unadjusted fabric can cause layer disturbances. We examine the effects of unadjusted fabric on the deformation, using both an analytical formulation for a given fabric (Thorsteinsson,2000, 2001), and a model where the fabric evolves with the deformation (2000, 2002).

2. Overturned Folds

2.1. Particle paths and deformation

For a homogeneous three-dimensional deformation field in vertical uniaxial compression or pure shear combined with simple shear, we can solve for the position of a particle at any given time using the coupled first-order differential equations (Ramberg,1975)

(1)

where and is the velocity gradient; u_i is the velocity component in direction where z is vertical).

The particle path, in steady state, is then given by the solution to Equation (1)

(2)

where x₀ = x(t =0) and z₀ = z(t = 0). This solution assumes that the displacement at the origin x_i = 0 is zero, and that the L_ij’s are constants in time and space. Because ice in an ice sheet tends to move into regions of increasing shear, this assumption is not strictly correct. However, Waddington and others (2001) show that the solution (Equation (2)) is still useful, and Jacobsen (2001) shows that finite-strain calculations can illuminate layer behavior where the simple solution (Equation (2)) is not applicable. We use Equation (2) to derive a criterion for folding below. All the shear deformation takes place in the plane, and we can formulate our equations in two dimensions, noting that variations in the y direction are allowed.

2.2. Non-dimensional numbers

It is convenient to define several non-dimensional numbers, in order to simplify the equations. First we define the slope m

(3)

where A and B are two points on a layer, connected by a straight line (see Fig. 1).

We define a shear number, S, as

(4)

The shear number is a measure of how easy it is to deform the material in horizontal simple shear, relative to vertical pure shear or uniaxial compression.When S is large, shear dominates over compression and folding is likely. It is also convenient to define a stress ratio, as

(5)

where τ is the shear stress and σ is the vertical compressive stress.

Using local strain rates to define a characteristic time-scale

(6)

we non-dimensionalize time as T = t/T. We also use L, a characteristic length scale for the layer disturbance, to non-dimensionalize spatial variables, i.e. X = x/L and Z = z/L.

2.3. Passive folding

To derive a criterion for folding, we consider two particles, A and B, as shown in Figure 1. Particle A at time t = 0 is at x =0; z = A, x_A(t =0; x₀ =0; z₀ = A), and particle B is at x_B(0;L,0). for overturned folds to occur, particle A must reach an x position greater than that of particle B at some later time, x_A(t;0,A) > x_B(t;L,0). Waddington and others (2001) used a similar folding criterion. An equivalent approach is to consider when a line with an initial slope m₀ = A/L reaches vertical, and consequently, since the velocity gradient is uniform, overturns. By letting t → ∞ we find the smallest initial slope that can reach vertical. Using Equation (2) to solve for x_A(t; 0,A) – x_B(t;L, 0) >0, where t → ∞, we get

or using the non-dimensional notation

(7)

The time it takes particle A to overtake particle B (x_A(t; 0,A) – x_B(t; L, 0) = 0) is

or in non-dimensional form,

(8)

As discussed by Waddington and others (2001), the shear number S is not uniform along particle paths as we assume here; layers tend to move toward regions of increasing simple shear, and the overturning time (Equation (8)) tends to be an upper limit. By calculating the finite strain along particle paths, Jacobsen (2001) assessed the error introduced by the uniform strain-rate assumption. the time error is small, and the assumption a good one, except in regions where S_m0 ≈1, and in regions strongly affected by bedrock topography.

2.4. Anisotropic strain rates

Ice obeys a power law with a stress exponent n ≈3; most other minerals of geophysical interest also haven > 1. the analytical anisotropic flow law that we use here is based on the homogeneous stress assumption (often referred to as the Sachs (1928) model). the derivation of the velocity gradient, L_ij, as a function of cone angle can be found in 2000, 2001. Cone angle α is a half-apex angle of a cone within which the c axes of all the crystals are assumed to be uniformly distributed. Vertically symmetric fabrics are widespread in polar ice sheets (Gow and others, 1997; Thorsteinsson and others, 1997). the Appendix gives the equations for the non-zero components of the velocity gradient in uniaxial compression combined with simple shear (UC&SS) stress state, and in combined pure and simple shear (PS&SS) stress state, as functions of cone angle. the stress states considered are

(9)

The anisotropy breaks the one-to-one correspondence between stress and strain-rate components that exists for isotropic materials. We can examine how the anisotropy changes the deformation rate, relative to the isotropic response, by plotting the ratio of anisotropic deformation rate to the deformation rate of a corresponding isotropic material subjected to the same stress state.

Figure 2 shows the normalized vertical deformation rate, as a function of cone angle α and the stress ratio in combined uniaxial compression (σ) and simple shear stress (τ) state. the normalized deformation rate can be viewed as an enhancement factor for this particular deformation-rate component. the peak enhancement in vertical strain increases with increas- ing shear stress, and the maximum (for a given value of ) moves from α ≃ 60˚ for to α ≃ 35˚ for The vertical strain rate varies greatly as a function of cone angle. Figure 3 shows the component. the shearing rate L₁₃ depends less on the relative magnitude of shear stress to compressive stress when When the enhancement varies greatly with both and α.

Fig. 2 The normalized vertical deformation rate, as a function of stress ratio and cone angle α, in Note how the compressive deformation rate has a maximum near 60˚ when and near 35˚ for larger shear stress.

Fig. 3 The normalized shear deformation rate, as a function of the stress ratio and cone angle α, in

We now examine the dependence of the shear number S on the anisotropy and stress. Figure 4 shows the shear num- ber S as a function of cone angle and stress ratio in combined uniaxial compression and simple shear (UC&SS) stress state. the folding criterion (Equation (7)) shows that the slope, m, must be larger than S–¹ for overturned folding to be likely. Figure 5 shows the shear number S as a function of cone angle and shear stress in combined pure and simple shear (PS&SS) stress state. Folding is more difficult in PS&SS than UC&SS since, with equal layer thinning, the extension in the x direction (and therefore layer flattening) is stronger in PS&SS. for isotropic material obeying a constitutive equation of the form (Glen’s flow law (Glen, 1958)), where σ_e is the effective stress, it can be shown that in UC&SS and in PS& SS.

Fig. 4 Shear number S as a function of cone angle and stress ratio in UC&SS stress state. the slope m₀ of a layer disturbance must be greater than S–¹ for folding to develop.

Fig. 5 Shear number S as a function of cone angle and stress ratio in stress state. the slope m₀ must be greater than S–¹ for folding to develop.

The evolution of layers is calculated using Equation (2). Figure 6 shows the evolution of the same initial wrinkle in isotropic and anisotropic layers subjected to the same stress field. the initial amplitude is A =1, the initial base length is L = 10, and the initial slope of a straight-line segment between points A and B is m₀ = 0.1. the stress state is PS& SS, with Under these conditions, point A will never overtake point B in the isotropic medium (Fig. 6a); overturned folding is unlikely to occur. In the anisotropic medium, on the other hand, folding occurs rapidly (T _f ≈ 0.35) under the same conditions (Fig. 6b).

Fig. 6 Evolution of a wrinkle in PS&SS stress state. the ratio of shear to compressive stress is(a) Isotropic material, α = 90˚. the line segment between A and B has a slope m₀ = 1/10, and the top point A will never overtake the bottom point B (see Fig. 1). (b) Strongly anisotropic material, with a cone angle α = 20˚. In this case the wiggle overturns rapidly. the layer profiles are shown at equal time intervals, and at the same time in both panels.

3. Origin of Layer Disturbances

Near ice-sheet margins, horizontal compression tends to fold even apparently insignificant layer disturbances. Near a divide, however, horizontal extension (Fig. 1) prevents most disturbances from evolving into recumbent folds. Therefore special conditions are required to create layer disturbances that will fold. Unsteady ice-divide migration might provide layer disturbances on a scale comparable to the ice thickness (Jacobsen,2001), just as unsteady flow near margins can cause folds (Hudleston,1976). There are also several possible sources for small-scale layer disturbances near ice divides. Sastrugi or accumulation variations at the surface seem unlikely sources for the small-scale disruption of stratigraphy observed at GRIP and GISP2 (Waddington and others, 2001). Variations in rheological properties seem to be the most likely source of small-scale layer disturbances near ice divides. Differences in stiffness between horizontal layers in horizontal extension can cause boudinage (Smith, 1977; Staffelbach and others, 1988; Cunningham and Waddington,1990), while variations in the anisotropy may lead to layer disruption in several ways. Inhomogeneous anisotropy causes inhomogeneous deformation, which after some time leaves wrinkles. Thus, if we can create wrinkles on a previously flat layer we have a mechanism to make folds. the problem then reduces to how to generate spatially heterogeneous fabric in flat layers.

3.1. Tilted cones

If the symmetry axis of the cone is non-vertical the deformationwill be a combination of pure and simple shear, even if the imposed stress state is only vertically directed pure shear stress (Azuma and Goto-Azuma,1996; Thorsteinsson, 2000). In contrast, vertically applied pure shear stress can only flatten and smooth disturbances in horizontally lying isotropic layers. to explore the deformation of a layer with variably tilted symmetry axes, we prescribe a sinusoidal variation of the tilt (to a maximum of ±40˚) of the cone symmetry axis (see Fig. 7a). This choice of spatial function for the tilt is completely arbitrary; the important point is that the tilt changes with X. We then apply vertical pure shear stress on the layer. Figure 7b shows the initially flat region within a layer, or a region, with a 20˚ cone whose symmetry axis tilt varies sinusoidally. Figure 7c shows the shape of the region after the applied deformation. the dashed box shows the shape of the same region if all the 20˚ cones are vertical; it has been subjected to ~2% vertical strain. Note that there is about 106 exaggeration for the vertical scale. the sense of shear is in the direction of the tilt (of the cone), creating the asymmetry in the resulting shape of the initial region. the ratio of maximum to minimum vertical strain rate, for a 40˚ and 0˚ tilt, respectively, is ~3.5. the troughs are displaced away from the center (X = 25), and the limbs thicken around the center. the spatially variable anisotropy has created wrinkles where none would be created if the ice were isotropic. Moreover, the wrinkles were created rapidly; note the small strain on the region with vertical cones (dashed box).

Fig. 7 (a) A layer with 20˚ cones that have a sinusoidal variation of the tilt angle, with a maximum tilt of ±40˚. (b) An initially flat region within a layer that has the fabric shown above, and (c) the same region after ~2% deformation in a pure shear stress field. the dashed box shows how the region responds if the cones are not tilted.

3.2. Evolution of tilted fabric

So far we have considered only cases in which the fabric is frozen (not evolving). In reality the fabric evolves during deformation, because the c axes of the ice crystals rotate due to intra-crystalline slip and constraints imposed by neighboring crystals (Wenk and Christie, 1991). We use the fabric evolution model described in Thorsteinsson (2000, in press) to calculate the c-axes rotation and the velocity gradient as a function of time. Since the deformation rate changes continually during the deformation, we retain dimensional time in these calculations. the calculations were done at –10˚C and a compressive stress of 20 kPa. Other flow-law and model parameters can be found in Thorsteinsson and others (1999) and 2000, 2001. However, our primary focus is on the differences attributable to cone tilts.

When the symmetry axis of the cone is tilted away from vertical, in vertical uniaxial compression, the strain rates are very different from those arising from the deformation of a vertically symmetric fabric. Figure 8 shows layer thickness as a function of time for layers with an initially 40˚ cone tilted by 0˚, by 20˚ and by 40˚ from vertical. the layer thickness h as a function of time is

Fig. 8 Evolution of layer thickness, in uniaxial compression, as a function of time for different tilts of the original cone fabric. the three lines show the layer thickness for the initial condition of a 40˚ cone tilted by: 0˚ (top), 20˚ (middle) and 40˚ (bottom). the fabric then evolves freely in vertical uniaxial compression stress (using Thorsteinsson (2002) model, with mild-nearest-neighbor interaction and σ =–0.2 bar).

(10)

where H is the initial thickness of the layer. Since it is easiest to deform crystals oriented close to 45˚, the cone tilted by 40˚ has the highest vertical strain rate initially. the fabric-evolution calculations were carried out to 0.25 equivalent strain. the cone tilted by 40˚ accumulated 0.23 vertical strain at t ~ 8 kyr, the cone tilted by 20˚ accumulated 0.17 vertical strain at 10 kyr and the cone tilted by 0˚ accumulated 0.13 vertical strain at 10 kyr. If there is no fabric evolution (frozen fabric) then the cone tilted by 40˚ accumulated 0.21 strain at t ~ 8 kyr, the cone tilted by 20˚ accumulated 0.24 strain at 10 kyr and the cone tilted by 0˚ accumulated 0.21 strain at 10 kyr. the c axes rotate toward vertical, and, as Figure 9 shows, in the process they do not preserve the cone fabric, in the sense that the c-axis distribution tightens and migrates toward vertical. In addition, the boundary enclosing all of the axes can become elongated, and the c-axis density nonuniform.This demonstrates that a region of tilted cones will evolve to look very similar to adjacent regions with vertical fabric. the rate of layer thinning depends on how many crystals are close to, and pass through, the soft orientation, which is 45˚ in this case. Therefore the total strain for the cone tilted by 40˚ is greater when fabric evolution is active than for frozen fabric, but smaller when the tilt is 0˚.

Fig. 9 Fabric evolution in uniaxial compression for 40˚ cones. the fabric is shown on Schmidt plots. Initial fabric is shown in the top row, and fabric after 8 kyr in the lower row. the cone is initially (a) vertical, (b) tilted by 20˚, and (c) tilted by 40˚.

Tilted cones are therefore transient features capable of generating layer disturbances.

3.3. Stripes

In a study of the GISP2 ice core, Alley and others, 1997 found that the expected vertical c-axis fabric was interrupted by planes of grains, with the c axes oriented approximately in the dip direction of the planes. These features were called stripes because of their appearance when intersected by vertical thin sections. the stripes are typically one to very few grains in thickness. If they are orientated close to vertical, their length is a few grains, but if they are orientated further away from vertical, they tend to be longer. Alley and others, 1997 hypothesized that stripes form through organized polygonization, such that after polygonization of the center grain the stress on neighboring grains will tend to rotate them so that they will be parallel to the polygonized grain. the stripe would then grow through internal spinning of grains at the end of a stripe, i.e. shear on the plane of the stripe would force the grain at the end into the same orientation.

From our exploration of the tilted cone fabric above, we know that non-vertical alignment of crystals will complicate the deformation, and potentially give rise to layer disturbances. Stripes have a tilted symmetry axis, relative to the bulk ice, and are therefore a potential source of small-scale folding. Moreover, if the density of stripes is spatially variable, increasing with depth for instance, then stripes might affect the flow on a larger scale.

4. Summary and Conclusions

Knowledge of the processes responsible for small-scale folding near ice-sheet centers is important in order to assess the integrity of ice-core records. Understanding the origin of layer disturbances near ice-sheet centers is a challenging problem. Several authors have suggested that anisotropy may play a significant role in that process (Azuma and Goto-Azuma,1996; Castelnau and others, 1998). We simulated the effects of tilted cones on the shape of a horizontal layer and found that tilted cones can indeed produce distorted layers in stress states that would not produce any wrinkles with isotropic ice. Accounting for fabric evolution does not change that conclusion.

The presence of rheological inhomogeneity and anisotropy changes the distribution of stress. In general, if the ice does not fracture or tear, then stress will be increased where the ice is relatively stiffer with respect to the applied stress. Similarly, where ice is made softer, stress will decrease. By taking the view that we can specify the stress pattern, we cannot incorporate this stress redistribution in our analysis. Therefore our results should be interpreted as a qualitative or semi-quantitative guide to the effects of anisotropy on folding.

We have shown that if regions of unadjusted fabric, i.e. fabric with different symmetry than expected in the current deformation field, exist in ice sheets, that fabric is a likely source of subsequent layer disturbances. the remaining question is then about the source of the anomalous fabric. Evidence from ice cores, such as stripes, shows that local heterogeneity can form, but the concentration or distribution of such features is not well known.

Under steady-state stress conditions, the vertically oriented anisotropy that tends to develop near ice domes can greatly reduce the slope of disturbances needed to generate overturned folds. Flow instabilities associated with an asymmetric fabric in initially flat-lying layers could be a source of layer disturbances that evolve into overturned folds, even in the dominantly extensile flow regime near ice divides.