Geometry

The geometry of the idealized ice divide is shown schematically in Figure 1. The height of the right-hand edge of the solution region was taken as one unit. It was placed horizontally 19 units from the divide, which turned out to be far enough that assumptions about the details of variation of velocity with depth on the right-hand boundary did not affect the solution near the divide at the left. The specific dimensions are not critical to the results which are described below in dimensionless quantities. Several surface elevation profiles were considered (Fig. 2).

Fig. 1.

Idealized cross-section of two-dimensional symmetric ice divide.

Fig. 2.

Range of ice sheet surface elevation (a) and corresponding horizontal (b) and vertical (c) velocity profiles for which calculations were made. Differences in shading and outline identify correspondence between the ranges. For simplicity of presentation, only results for divide thickness of 1150 m are shown in (c). Heavy solid curves show the reference surface elevation and corresponding velocity profiles approximately in balance with a uniform distribution of net balance on the upper surface. Light solid curve in (a) shows Vyalov’s profile fit from thickness and flux at right-hand edge.

With reference to the axes in Figure 1 (x horizontal, y vertical upward), it is assumed that the geometry and corresponding flow do not vary with z. More specifically, x and y components of velocity u and v may be non-zero but velocity in the z-direction is zero. Under these conditions the non-zero deformation rate components are:

Following Reference NyeNye (1957) an effective shear strain-rate may be defined as

Flow law assumptions

It is assumed that the components of deviatoric stress and deformation rate are related by a power law (Glen’s law). Restricting attention to non-zero components this is expressed as

where

In Equation (3b), n was chosen to be 3 and B to be constant. In a later section some consequences of variations of B resulting from non-isothermal conditions will be discussed.

An effective shear stress τ corresponding to ε may be defined as

The assumed flow law implies viscous incompressibility which for the two-dimensional flow is expressed as

Densification of snow and firn near the surface is not considered in these calculations.

Components of stress are given by

where p is the mean compressive stress and is not determined by the deformation rate and flow law.

Equations of motion

It is assumed that the flow is quasi-static, in which case,

where ρg is density times gravitational acceleration and is assumed constant.

Boundary conditions

Boundary conditions are illustrated in Figure 1. The upper surface was assumed to be stress free, the lower surface not to slip or melt, and the left-hand edge (x = 0) to be a plane of symmetry, which requires horizontal motion (u) and surface-parallel shear stress (τ_xy ) to be zero.

At the right-hand edge (x = 19l) the depth profile of horizontal velocity was chosen to have the form

where u _s(x) represents horizontal velocity along the upper surface. The shape of this profile corresponds to that for “laminar flow” of isothermal ice assuming Glen’s law with exponent n (Reference NyeNye, 1952). The boundary condition on u given by Equation (8) imposes a flux per unit width across the right-hand edge of

u _s(19l) was chosen to give q(19l) equal to 19lb, where by continuity b is the average downward velocity of the surface ice over the distance between the divide and the edge. This value of u _s(19l) is 23.75b for n = 3.

To establish a profile of vertical velocity at x = 19l, it is assumed that the profile shape given by Equation (8) also holds for x incrementally close to 19l, in which case differentiation of Equation (8) with respect to x gives

Equation (10) neglects small terms involving ∂h/∂x. These terms could be included and expressed in terms of ∂u _s/∂x if one were to assume a dependence of u _s on h. Such a dependence, being a priori unknown, was not assumed. Integration of ∂v/∂y upward from the motionless base using Equations (1), (5), and (10) gives

where

v _s(19l) was chosen to be −b to correspond to the average vertical velocity over the surface between the divide and 19l imposed by the choice of u _s(19l).

Choice of parameters

To further specify the problem it is necessary to choose l, b, B, and ρ, which set the geometrical scale, mass transport, ice viscosity, and body-force density from gravity. In this case these were chosen as: l = 1 km, b = 0.1 m a ⁻¹, B = 3.0 bar a^1/n , ρ = 0.9 Mg m⁻³. These are not intended to represent any particular ice sheet, but are chosen to fall within the range of these parameters for ice caps and ice sheets.

In order to examine how a solution for a particular choice of parameters can be scaled, consider dimensionless quantities

, , , . This non-dimensionalization employs l as a unit of length and b as a unit of velocity. Units of deformation rate (b/l) and stress (ρgl) are derived from these. If the non-dimensional quantities are introduced into the governing equations and boundary conditions (Equations (1) through (11b)), the equations are identical for different choices of l and b as long as the resulting dimensionless flow-law parameter

is the same. This guarantees complete geometrical and dynamic similarity of solutions. These calculations consider one case

from which solutions for equivalent combinations of l, b, B, and ρ can be derived. The actual range for thin ice caps with high accumulation and large ice sheets with low accumulation runs from about 0.3 to 14 × 10⁻³. A second case of substantially different scale is examined below.

Calculation method

The velocity distribution was calculated using quadrilateral finite elements for incompressible two-dimensional flow. The velocity within an element is determined in terms of the velocity at the four corner nodes based on subdivision along diagonals into four constant strain-rate triangles with velocity at the virtual center node chosen to give constant dilatation throughout the element. Pressure is assumed constant in an element. A solution region is divided up into elements and nodal velocities and element pressures are found to satisfy the equations of motion, incompressibility, and the boundary conditions cast in variational form. Non-linearity arising from the flow law is treated using a Newton–Raphson iterative solution technique. More information about the numerical method is given by Reference RaymondRaymond (unpublished) and Reference HookeHooke and others (1979).

Calculations were done using several grids in order to gain some idea of the grid dependence of the solution. Grids of 6 × 20 and 11 × 20 were used. First they were applied to the full 19l length of the solution region. In order to obtain a closer longitudinal spacing of nodes, the same grids were used over the distance 0 to 9.5l using results from the full profile to define boundary conditions on the right-hand edge at 9.5l. Solutions for surface velocity from the different grids show agreement to two significant figures or better.

Choice of surface-elevation profile

A number of surface profiles between the divide and 19l were considered (Fig. 2).

The heavy solid line in Figure 2a represents the surface topographic profile chosen for principal attention in the following analysis. It will be called the reference surface profile. It is represented by

where

, and . The parameters , and represent dimensionless divide thickness, limiting surface slope approaching the divide, and slope gradient of the surface. This particular surface-elevation profile was singled out because it gives a relatively uniform vertical velocity v _s along the surface (Fig. 2c), which would be approximately in local balance with a uniform distribution of annual net accumulation of b. No attempt was made to find a profile which would give a perfectly uniform distribution of v _s.

One interesting feature of the profile given by Equation (12) is a non-zero slope as the divide is approached. Apparently any rounding of a divide crest happens on a scale smaller than one ice depth, which is the resolution limit of the numerical calculation grid. From the calculation results, this seems essential to match uniform balance distribution. If the slope were flatter (or zero) on a scale of one ice depth, then the outward flow from the divide and the corresponding downward motion there would be too small (Fig. 2b and c). This feature of a rather sharp divide crest seems to be characteristic of real ice divides (see e.g. Reference Paterson and ColbeckPaterson, 1980, fig. 6; Reference BuddBudd, 1969, fig. 6.4; Reference WeertmanWeertman, 1961, fig. 2), although they must be rounded on some small scale.

Although the profile given by Equation (12) was selected from a number of more or less arbitrary trials, it turns out to be represented well by one derived theoretically assuming ice transport is related to thickness and surface slope through laminar flow theory (Reference VyalovVyalov, 1958) (see Fig. 2a).

Other choices of profiles in the range shown in Figure 2a gave velocity distributions which would require rather bizarre distributions of accumulation rate to produce steady geometry (Fig. 2c). Even though these other profiles are unrealistic, the results from them are of interest. Below they are examined to test the sensitivity of the pattern of variation of deformation with depth near the divide to differences in geometry and surface velocity pattern.

Velocity versus depth

The finite-element solution assumes that velocity varies linearly within subtriangles of the quadrilateral elements. Figure 3 shows computed profiles as the computation method actually sees them. Profiles of velocity versus depth are normalized to allow comparison of profile shape at different longitudinal positions x. The shapes do not depend on x except within a zone near the divide having a width of several ice depths. Profiles for horizontal velocity u (Fig. 3a) approximate that predicted by laminar flow theory except near the divide, where shearing is more uniformly distributed over depth.

Fig. 3.

Normalized velocity vesus normalized depth for various distances from the divide. Distances are in units of l (approximately ice depth). Heavy lines show finite-element results from the reference surface-elevation profile (solid line of Fig. 2a). Bars in (b) indicate the range found at the divide and positions x/l > 4 for all surface profiles shown in Figure 2a. Light lines show approximate analytical descriptions of depth profiles at and distant from the divide. Dotted light lines show depth profiles for flow of a linear fluid in the reference geometry with viscosity chosen to give nearly uniform vertical velocity at the upper surface. For n = 1 the profile shapes are not sensitive to distance from the divide.

The corresponding profiles of vertical velocity v (Fig. 3b) show downward velocity, which is significantly less at all depths than expected from a linear depth variation and approaches a parabolic depth variation at the divide.

Deformation rate and deviatoric stress versus depth

The finite-element calculation assumes deformation rate and stress are constant in subtriangles of the quadrilaterial elements; therefore, the solutions for these quantities are step functions along any profile through the body. Results are here plotted as smooth curves through points evaluated at the mid-points of segments traversing elements. Figure 4a and b shows profiles over depth at the divide and 10l from the divide. This latter position was chosen to be both remote from the divide so to be characteristic of a flank position (Fig. 3) and also remote from the right-hand edge, so to be free from details about the assumed boundary condition there.

Fig. 4.

Deformation rate and stress versus depth. Heavy lines show finite-element results from solid-line surface elevation profile of Figure 2a. In (a) and (b) light lines show possible analytical descriptions discussed in text. Equations (13) define ε⋆ and τ⋆. Numbers indicate distance from divide in units of l, which is approximately ice depth.

The results for deformation rate are shown in units of

which is close to the vertical strain-rate at the divide surface. Deviatoric stress is shown in corresponding units of

defined by the flow law represented in Equations (3) and (4).

At the divide the vertical strain-rate d_yy = −d_xx varies approximately linearly with depth (Fig. 4a) with a consequent fairly simple distribution of deviatoric stress which decreases toward the bed (Fig. 4b).

Remote from the divide, d_yy = −d_xx tends to be nearly constant near the surface with the principal decrease to zero at the bed occurring in the lower half of the depth (Fig. 4a). As a result

and τ are also roughly constant near the surface. At greater depth they increase strongly under the influence of d_xy and τ_xy which are the principal contribution near the bed. A simple feature is that τ_xy increases linearly with depth. In response to τ_xy , d_xy increases with depth. A linear increase of d_xy near the surface is similar to what would be expected for a linear fluid of constant viscosity. In this case it occurs in response to a nearly constant effective viscosity near the surface controlled by the nearly depth-independent vertical strain-rate there as explained by Reference NyeNye (1957).

At the divide τ_xy = 0 by symmetry, and distant from the divide it varies linearly with depth. Figure 4c shows the evolution of the depth profile of τ_xy as the divide is approached. Within 2l of the divide, there is a distinct deviation from a linear variation versus depth.

Approximate analytical description of velocity, deformation rate, and deviatoric stress at the divide.

At the divide, the vertical variation of v is close to being parabolic (Fig. 3b), and may be represented approximately as

Symmetry at the divide and Equations (14), (1), (2), (5), and (13a) indicate the relevant deformation-rate components would be

These are plotted in Figure 4a as light lines where they differ perceptibly from the numerical results. The deviations are small.

From symmetry at the divide and Equations (1) and (15a), it follows that

to first order in x. Therefore Equation (14) implies that the depth profile of u approaches being linear as the divide is approached (x → 0), which seems compatible with Figure 3a.

Deviatoric stress components corresponding to Equations (15) may be calculated from Equations (3), (4), and (13) to find

These do not differ perceptibly from the numerical results on the plotting scale of Figure 4b.

Although Equations (13) to (17) appear to be reasonable approximations, it is important to recognize that they apply only in the limit x → 0 and even then are not exact. This is apparent from the free upper-surface boundary condition which requires τ_xy , d_xy , and ∂u/∂y to be small at y = h when ∂v/∂x and ∂h/∂x are small. In fact Figure 3b shows that v(0, y)/v _s(0) determined numerically is larger than given by Equation (14) near the upper surface; correspondingly the profile of u(x, y)/u _s(x) would maintain some upward concavity near the surface right up to the divide. At the base, an expansion of the stress to first order in x shows that Equations (14) and (16) lead to an unrealistic stress singularity at u = 0. To avoid this ∂u(x, 0)/∂y must be zero to first order in x. This would imply the profile u(y) would become concave downward near the base as the divide is approached. This tendency seems to show up in the numerical results one half depth from the divide. Apparently the precise profile of u versus y as x → 0 has an inflection.

Approximate analytical description of velocity, deformation rate, and deviatoric stress distant from divide

Figure 3a and b shows that the vertical profiles of u and v have shapes approximating those given by Equations (8) and (11), when x ≧ 4l. These may then be used to estimate deformation rate and deviatoric stress.

The x-differentiated form of Equation (8) given by Equations (10) and (11b) indicate

where

Differentiation of Equation (8) with respect to y and the assumption that ∂v/∂x is negligible gives

where

The corresponding effective strain-rate is

From this is seen that ε _s is the value of ε at the upper surface (y = h) where ε is dominated by d_xx =−d_yy , and ε _b is the value of ε at the bed y = 0 where ε is dominated by d_xy = d_yx .

The deviatoric stress components and effective stress are easily calculated from Equations (18a), (18b), and (18c) and the flow law (3a), (3b), and (4).

The resulting deformation rate and deviatoric stress components are plotted in Figure 4a and b as light lines where they are distinguishable from the numerical results. The major deviation occurs near the surface where d_xy from Equation (18b) is distinctly less than calculated numerically. This is apparent also in Figure 3a where it is seen that the numerically determined velocity profiles show shearing more evenly distributed over depth than for the laminar flow profile on which Equation (18a) is based. This discrepancy exists because the laminar flow profile neglects near-surface softening by the contribution of vertical strain-rate to the effective viscosity (Reference NyeNye, 1957).

Longitudinal variation of surface velocity and basal shear stress

The longitudinal variation of u _s (Fig. 2b) and τ _b (Fig. 4c) are shown in Figure 5a and b for comparison with predictions from standard “laminar-flow” theory (Reference NyeNye, 1952). Distant from the divide the agreement is good, but within 4l distinct deviations are apparent. This is consistent with expectations from consideration of the depth profile shapes for u and τ _xy considered previously (Figs. 3 and 4).

Fig. 5.

Surface velocity and basal shear stress versus distance from divide. Heavy lines show finite-element results calculated for the reference surface elevation profile of Figure 2a. Light line shows corresponding predictions from “laminar-flow” theory.

Of particular interest is the rather sharp drop in basal shear stress in the near-divide zone. At an intuitive level it may be understood as the result of two constraints: first, basal shear stress must go to zero at the divide by symmetry; second, it must maintain a high value to allow significant flow to take place because of the non-linear dependence of flow velocity on basal shear stress.

Dependence of velocity at the divide on geometry, scale, and mass balance

A number of different symmetric surface elevation profiles were examined (Fig. 2a). Some of these yielded rather complex distributions of surface velocity (Fig. 2b and c). For all of these, the shape of the depth profile of v at the divide is closely the same as indicated by the range bars in Figure 3b.

To test the consequences of a more extreme variation in divide geometry and also to consider a different scale, calculations were made for the geometry shown in Figure 6. This geometry approximates a flow-line profile across the divide of the Devon Island ice cap in the vicinity of the Devon Island core-hole sites (Reference PatersonPaterson, 1976), where ice depth is about 300 m and balance rate is about 0.2 m a⁻¹. The flow-line profile is based on surface and bed elevation maps from radio echo-sounding (personal communication from W. S. B. Paterson).

Fig. 6.

Geometry and finite-element grid simulating a flow-line profile of the Devon Island ice cap. Curves show stream lines determined from finite-element calculated velocity field. Temperature was allowed to vary with depth and temperature-dependent flow law was used as explained in text. Spacing of vertical grid lines is 400 m. The finite-element grid extended over a 6.8 km length of the profile approximately centered on the surface divide. Edges of the calculation region are omitted in the figure. Vertical exaggeration × 2.

For the calculation, B in the flow law (Equations (3)) was chosen in several ways. First, it was chosen to be constant B = 2.06 bar a^1/n so that the solution matched the surface horizontal velocity measured at the core-hole sites. In the divide region of the Devon Island ice cap, the temperature varies from −23 °C near the surface to about −18 °C at the bed. To account for this, B was chosen to vary with temperature as

with B _r = 2.06 bar a^1/n , T _r = 253 K (−20 °C), Q = 59 kJ mol⁻¹ (14 × 10³ kcal mol⁻¹), and R is the universal gas constant, which with the measured variation of temperature with depth gives results distinguishable but not dramatically different from the case of B constant. Figure 6 shows the flow pattern for the non-isothermal case.

Here we will not consider how well calculated results can be fit to the full range of measured flow and deformation data on the Devon Island ice cap, which are complicated by compaction of firn at the surface and transverse strain-rates. Attention is focused on the theoretical depth variation of v at the flow divide. This particular divide geometry shows considerable complexity in comparison to the isothermal symmetric flat-bedded case considered in the previous section. Foremost are the uneven bed, a shift of the surface topographic divide from the bed divide, and a shift of the flow divide from the surface topographic divide showing some features suggested by Reference Haefeli and KingeryHaefeli (1963). Also the scale is different; the parameter

(in comparison to for the reference symmetric divide). Even with these major differences in divide shape and scale, the shape of the calculated depth profile of v at the flow divide is consistent with Figure 3b, with results for both isothermal and non-isothermal variations of B falling in the range bars shown.

Dependence of velocity distant from the divide on geometry

The shapes of depth profiles of u and v away from the various divides shown in Figure 2a were fairly consistent and held to the patterns shown in Figure 3a and b. However, the deviations were larger than at the divide, as shown by the range bars in Figure 3b.

Deviations are expected from two sources: first, a wide range of longitudinal strain-rates, which will result in some differences in the depth distributions of effective viscosity (Reference NyeNye, 1957), and second, a wide range of longitudinal strain-rate gradient and associated longitudinal stress gradient (Reference BuddBudd, 1968), which will result in differences in stress distribution and may affect the depth variation of the shear strain-rate parallel to the surface. Some effect from the first of these has already been noted. The second appears to be the more significant, because the largest deviations of profile shape occurred in cases where |∂² u/∂x ²| was largest. In this regard, the results show that when ∂² u/∂x ² > 0, the depth profile of u is modified such that shear becomes more evenly distributed over depth and the depth profile of v becomes more concave downward, that is, such that vertical motion is more confined to near the upper surface. Opposite effects on profile shapes are found for ∂² u/∂x ² < 0. These are interesting systematic effects of longitudinal stress gradients, which are simply noted here but will not be considered in detail. Based on the fact that ∂² u/∂x ² = 0 at the flow divide by symmetry and the predominant effect on velocity depth profile shape from ∂² u/∂x ², we also see a partial explanation of why the profile shapes at divides are so uniform.

Effect of ice rheology

To test what role the non-linear flow law of ice plays in the special pattern of deformation near the divide, calculations were done for n = 1, which corresponds to a linear Newtonian fluid.

When the viscosity, surface topographic profile, and right-hand boundary condition are matched, so that v and ∂u/∂x are nearly uniform along the upper surface, then the variation of u over depth shows a parabolic depth variation consistent with a Poiseuille flow profile. Furthermore, this shape persists up to the divide at least within the longitudinal resolution of the numerical calculations. This x-independent profile shape for u and the corresponding one for v are shown in Figure 3.

When viscosity, surface topographic profile, and right-hand boundary conditions are such that large variations in v and ∂u/∂x occur along the surface, then deviations from these Poiseuille profile shapes over depth are perceptible and show the same systematics depending on the sign of ∂² u/∂x ² as in the non-linear case. This shows that the changes in profile shape associated with ∂² u/∂x ² being non-zero arise principally from a redistribution of the depth variation of shear stress τ_xy associated with longitudinal stress gradients and only secondarily from rheological non-linearity. A simple explanation will be considered separately in a paper in preparation.

Width of the divide zone

Since there is no distinguishable near-divide zone for n = 1, a special stress distribution pattern imposed by the divide geometry alone is not the principal factor. In contrast to the case for n = 1, in which viscosity is constant, when n = 3 the effective viscosity varies with depth. At the divide effective viscosity increases downward. Distant from the divide, it increases upward. These facts are seen from the distribution of ε in each case (Figure 4a and b) and the flow law (Equation 3b). The principal cause of the near-divide zone must be this difference in effective viscosity arising from the non-linearity. The difference arises from general constraints. At the divide, ε is small and correspondingly viscosity is large at depth because ∂u/∂y = 0 by symmetry and ∂u/∂x → 0 at the bed because of its frozen condition. Away from the divide, ∂u/∂y becomes large at depth yielding relatively low viscosity there.

To examine the width of the divide zone, one is led to compare longitudinal strain-rate ∂u/∂x to horizontal shearing rate ∂u/∂y. Averaged over depth, these are: 〈|∂u/∂x|〉 = |v _s/h|, 〈∂u/∂y〉 = u _s/h. Thus ∂u/∂x is larger than ∂u/∂y on average only when |v _s| > u _s. Assuming steady-state geometry such that v _s is adjusted to a uniform mass balance and horizontal velocity averaged over depth is proportional to u _s with a constant factor k, volume conservation implies −v _s x = ku _s h. For a laminar flow profile shape, k is (n + 1)/(n + 2), and for the shapes close to the divide it approaches 1/2, so it is not precisely constant. However, for the present argument, it is sufficient that it is nearly constant and close to one. That implies |v _s | > u _s when x < h. From this argument, the width of the near-divide zone is scaled to one ice depth independent of divide thickness h and divide balance. The one ice depth represents a scaling of a transition zone and it is not inconsistent with discernible effects out to several ice depths as occurs for n = 3 (Figs 3, 4, and 5).

Effect of non-isothermal conditions

The example based on the Devon Island ice cap geometry showed that if the temperature variation over depth is not large (e.g. 5 deg), then the results derived above can still be reasonably applied. However, with larger temperature variations over depth they become invalid. Table I shows examples of the effect on vertical velocity at the divide for the reference surface profile caused by several assumed temperature distributions. For this purpose, temperature versus depth was parametrized as

where T _s is surface temperature and G is geothermal gradient. This profile form approximates expected theoretical and measured shapes of temperature profiles. For a typical geothermal gradient of 0.025 deg m⁻¹, the temperature difference over a 10³ m thick ice sheet would be about 16 deg and this gives results which fall distinctly out of the range of variation caused by geometrical and kinematic effects considered previously.

Table I.

Comparison of normalized vertical velocity v(o, y)/v(o, h) for various temperature profiles given by Equation (20). Ice sheet surface profile is given by Equation (12) with l = 10³ m. Results are based on 6 × 20 grid

Summary

The salient features of the flow pattern for the reference isothermal symmetric ice sheet are:

• i. depth profile shapes at the divide such that v is approximately parabolic and u approaches being linear;

• ii. depth profile shapes distant from the divide such that u is approximately consistent with predictions for laminar flow of isothermal ice with a corresponding v derivable from incompressibility.

• iii. the transition from near-divide to divide-distant patterns occurring over several ice depths.

Based on calculations from the different geometries: (i) holds relatively independent of divide geometry and corresponding surface motion pattern; deviations from (ii) exist for geometries such that surface longitudinal strain-rate gradient is large. Theoretically, the width of the near-divide zone relative to thickness is insensitive to thickness and surface velocity. By comparison with calculations for a linear fluid (n = 1), which shows no similar near-divide zone, (i) depends on the assumed non-linearity of the ice rheology and specifically that n = 3. The principal factor preventing general applicability of (i) for ice divides is variation of temperature with depth.

Age versus depth

Figure 7 shows age versus depth. At the divide, particles move vertically downward and age versus depth is found from

Fig. 7.

Age versus depth for the reference surface-elevation profile at divide (x = 0) and at x = 19l determined from integration along particle paths (solid lines) and from local vertical velocity (light line). Unit of time T⋆ is |v _s (0)|/h(0). which is the time that would be requiredfor a particle deposited at the divide to reach the bed if it were to maintain its initial downward velocity.

An approximation can be calculated from Equation 14 to find

when it is assumed that h and v _s are independent of time.

Away from the divide, age versus depth is properly calculated by integration along particle paths, which allows for the initial elevation of a particle and changes in vertical velocity pattern during its horizontal displacement. Particles near the bottom originated near the divide and experienced much of their vertical displacement near the divide. As a result, age at a given depth near the bottom is greater than would be calculated from Equation (21) using the local variation of v with depth (Fig. 7).

Relative advantages of divide and flank positions for coring

Ice divides have a number of advantages for coring: the oldest ice will be thickest there (Fig. 7); the profile shape of vertical velocity with depth and the corresponding pattern of vertical straining should remain unchanged with changes in thickness and mass balance (Fig. 3b); there is no horizontal displacement resulting in a complicated history of vertical straining associated with motion over irregular bed topography. Against these advantages must be balanced possible drawbacks arising from the narrowness of the divide-like pattern. Shifts of divide position of only several ice depths could result in a complex history of vertical displacement as a particle experienced alternately divide-like and flank-like patterns of vertical movement. Based on the results found here, ice divide positions appear to be best for the deepest ice, but flank positions may be more simple for shallow and intermediate-depth ice.