Basal-flow characteristics of a non-linear flow sliding frictionless over strongly undulating bedrock

Abstract

The flow field of a medium sliding without friction over a strongly undulating surface is calculated numerically. The results are used to elucidate the basal-flow characteristics of glacier flow and they are discussed with reference to known analytical solutions. Extrusion flow is found to become increasingly pronounced as the value of n, where n is a parameter in Glen’s flow law, becomes larger. For sinusoidal bedrock undulations, a flow separation occurs if the amplitude-to-wavelength ratio exceeds a critical value of about 0.28. The main flow then sets up a secondary flow circulation within the trough, and the ice participating in this circular motion theoretically never leaves it. The sliding velocity is calculated numerically as a function of the mean basal shear stress, the amplitude-to-wavelength ratio and the flow parameter n. For moderate and high slope fluctuations, the sliding velocity is significantly different from what would be expected from results based on the small-slope approximation.

Introduction

Knowledge of the characteristics of basal flow is mostly based on analytical solutions valid for Newtonian fluids and bedrocks that are sufficiently smooth so that , where a is the amplitude and k is the Wave number of the bed (Nye, 1969, 1970; Kamb, 1970; Morland, 1976a, b; Gudmundsson, 1997). The possibility cannot be excluded that the non-linearity of Glen’s flow law, and relaxing the condition , gives rise to flow patterns markedly different from those suggested by the analytical solutions.

In this paper. I use a numerical approach to calculate non-linear flow over strongly undulating sinusoidal beds. The sliding velocity u _b is calculated as a function of ε and n, and its sensitivity to changes in δ(δ: = 1/ kh) is investigated. The flow characteristics of basal glacier flow over perfectly lubricated beds for high values of ε and for a range of n values is analysed. This numerical work is an extension of a similar analytical study (Gudmundsson, 1997) which was limited to linear flow and to beds having small amplitude-to-wavelength ratios.

Notation

Theoretical Framework

Field equations, flow-law assumption and boundary conditions

The constitutive law for ice is taken to be

(1)

where are strain rates and are the deviatoric stresses defined by

(2)

is the second invariant of the deviatoric stress tensor

(3)

and δ _ij is the Kronecker delta. The field equations for the conservation of mass and momentum are

(4)

and

(5)

where v_i are the components of the velocity vector. The ice density ρ _r is assumed to be constant and the ratio of inertial forces to viscous forces to be small.

The bed profile is given by

(6)

where Z ₀ represents the z coordinate of the bed line. The coordinate system is tilted and the x axis makes an angle α with respect to horizontal. The z axis makes a right-angle with the glacier surface and the mean bed line (see Gudmundsson, 1997, Fig. 1). The amplitude-to-wavelength ratio is denoted by r(r : = α/λ) and is referred to as the (single wavelength) roughness of the bed.

The bed is considered to be sufficiently well lubricated as to exert no (or negligible) tangential forces on the overlying ice. This assumption is known as the assumption of perfect sliding.The shape of the bed line is assumed to be sinusoidal. For ε small, this represents no actual limitation, because to first order in ε the problem is linear and a Fourier decomposition can be used to obtain the flow field for arbitrary shaped bed lines (as long as is fulfilled for each of the Fourier components). Because of the intrinsic non-linearity of the boundary conditions (Gudmundsson, 1997), this is not true for moderate to high ε values. Using sinusoidal bed for high amplitude-to-wavelength ratios, however, facilitates comparison with analytical solutions and represents a convenient idealization of an undulating bed suitable for a study of the general characteristics of basal flow.

Ignoring the effect of regelation, the boundary conditions along the bed line are

(7)

and

(8)

where tan β(x) : = dz₀(x)/dx, is the local slope of the bed.

The upper surface is free but it is implicitly assumed that the ratio λ/h, where h is the mean thickness, is small enough to ensure negligible transfer of basal sliding variation to the surface. Formally, this means that the thinness parameter, δ : = 1/(kh), must be small compared to unity. Periodic boundary conditions are imposed in the horizontal direction.

Previous Work

If the bed line is sufficiently slowly varying so that , the problem can be solved for n = 1 using standard perturbation techniques. This has been discussed by Gudmundsson (1997), where references to earlier work can be found.

Some knowledge of the sliding velocity, u _b, can be obtained from dimensional analysis (Lliboutry, 1968, 1987b). There are three dimensionless numbers that enter the problem: n, ε and δ. For Glen’s flow law, dimensional arguments give

(9)

where U _b, is a non-dimensional sliding velocity.

By using a variational theorem for non-Newtonian flow (Johnson, 1960), Fowler (1981) derived a sliding law valid for a non-linear medium. He concluded that the sliding velocity u _b should vary as 1/ε ⁿ⁺¹ in the limit . This has, however, been questioned by Schweizer (1989) and Schweizer and Iken (1992) who, through finite-element (FE) calculations, found u ^b to less dependent on ε than indicates.

Raymond (unpublished), in what seems to be the first numerical study of this problem, calculated the sliding velocity for a perfectly lubricated sinusoidal bed with no regelation using Glen’s flow law. His calculations were limited in a few ε values and they do not give the sliding velocity as a general function of ε, δ and n but they are in agreement with .

If u _b is indeed proportional to ε⁻⁽ⁿ⁺¹⁾, then it is convenient to define a function s(ε, δ ,n), called the sliding function, as

(10)

s can be thought of as a non-dimensional sliding velocity where the ε⁻⁽ⁿ⁺¹⁾ dependency has been accounted for. Linear perturbation theory (Nye, 1969; Kamb, 1970) gives .

The function s(ε, δ, n) can be expressed as a Taylor series with respect to the variable ε

(11)

where c _i are the Taylor coefficients. Note that s must be an even function of ak, since the sliding velocity is independent of the sign of a. Fowler (1981), Lliboutry (1987a) and Meyssonnier (1983) were able to give a numerical estimate of C ₀(0, 3). Using variational methods, Meyssonnier found

(12)

Meyssonnier (1983), in a comprehensive numerical study, compared numerically calculated sliding velocities with analytical estimates. By calculating the sliding velocity for values of r in the range 0.01 – 0.05 and δ in the range 0.25/π to l/π, he was able to give an estimate of c ₂(0, 3), which was c ₂(0, 3) ≈ 2.4.

Lliboutry (1993) later improved Meyssonnier’s estimate, again using variational methods, and obtained an upper bound on s for n = 3. His result is

(13)

Kamb (1970) used the fact that the effective stress is, to first order in ε, independent of x as a starting point for the development of a theory of sliding incorporating rheological non-linearity. Kamb’s expression for the sliding velocity is (Kamb, 1970, p. 703)

(14)

Using this expression, one finds that c ₀(0,3) = 4/e² ≈ 0.54. This as a value is not within the range given by expression (12), showing that Meyssonnier’s and Lliboutry’s results and those of Kamb differ somewhat. In the light of the assumptions that Kamb had to make, the agreement is in fact surprisingly good. Kamb’s theory not only gives the correct dependency of the sliding velocity on roughness but the numerical values are also almost correct. Kamb’s non-linear theory was the first theoretical work that predicted the l/ε⁽ⁿ⁺¹⁾ behaviour of u _b for .

Schweizer (1989), in a numerical study, came to the conclusion that the sliding velocity does not vary as ε⁻⁽ⁿ⁺¹⁾ but somewhat more slowly. His results are also not in accordance with results based on a dimensional analysis ( Equation (9)). They must therefore be, to some extent, inaccurate.

Sliding Law for a Non-Linear Medium

Theoretical Considerations

It is instructive to see how it can be shown convincingly, albeit not rigorously, that as and for the sliding function s(ε, δ, n) approaches asymptotically a function which is only dependent upon n. A somewhat similar argument has previously been given by Lliboutry (1968, 1987b) and Kamb (1970). The argument, given below, suggests a certain form of the sliding law but does not prove that it must have this form. The proof of this result has been given by Fowler (1986, 1987).

The idea is to express the pressure variation along the bed as a function of the amplitude, the wavelength, the sliding velocity and the viscosity, and incorporating the effects of the non-linearity of the flow law by using an effective viscosity. The main underlying assumption is that

.

To this end, I start by writing a plausible expression for the pressure distribution, given by

(15)

where

(16)

C ₀ is some unknown constant, p _a, is the atmospheric pressure and η _eff is an effective viscosity that will be defined below. In order to see the plausibility of this expression, note that δp must vanish as . That δp(x) depends linearly on ε can be thought to be the result of a Taylor expansion with respect to ε where only the first term is retained. The product ηu_bk gives the right dimension for the pressure. It would also have been possible to use the basal shear stress T _b to get the right dimensions but then η would not have been a part of the expression. Note also that k cannot be substituted by l/a because, if must vanish (another way of seeing this is that if must change sign). It is also clear that the variation of δp with x must be of the form cos kx if ε is sufficiently small.

To get an estimate of the effective viscosity η ^eff, consider a column of ice which extends from the bed up to some distance l where the disturbance due to the bed undulations has disappeared. As the ice moves a distance Δx, in the time Δt = Δx/u_b, it is stretched/compressed to the length l + Δz, where Δz = Δxdz₀/dx. The (average) strain rate is given by

(17)

where the brackets are used to indicate that these are averaged strain rates over a given vertical distance l, that is

(18)

For the average vertical strain rate, one gets

(19)

where use was made of the fact that l must be proportional to 1/k, since the only parameters entering the problem that have the dimension length are a and l/k, and that using a would again be incorrect since the strain rates must change sign as . (This expression could also have been obtained by calculating the average of the vertical strain rates over the whole glacier using results from the first-order perturbation theory.) A similar argument for the shear rates lead to

(20)

The effective strain rate έ is

(21)

leading to

(22)

Glen’s flow law can be written as

(23)

where

(24)

Inserting Equations (22) and (24) into Equation (16) gives

(25)

Force equilibrium requires

(26)

where γ is a path along the sine curve, C(γ) is the path length, σ is the stress tensor and is a unit normal vector

(27)

A change of varibles dx = ds (1 + (dz₀/dx)²))^−1/2 gives for the drag (the x component of Equation (26))

(28)

where tan β = dz₀/dx. The last term on the righthand side is of second order in ε and the exact boundary condition on z ₀, Equation (8), shows that is then also of second order. Up to first order in ε, the drag can therefore be calculated according to

(29)

giving

(30)

where c ₀ is the same unknown constant introduced in Equation (16). After integrating and re-arranging terms, one gets

(31)

Expression (31) brings out the asymptotic 1/εⁿ⁺¹ behaviour of the sliding velocity as and shows how strongly it depends on the bedrock roughness.

Equation (31) shows that s(ε, δ, n) as a function of n will be equal to , where č is some unknown number. Since , it follows that

, leading to

(32)

It follows from the first-order perturbation solution for the pressure that c ₀ = 2 for a linear Newtonian medium (Nye, 1969, 1970; Kamb, 1970). Assuming that the effect of a non-linear flow law on the pressure variation along the bed can be described by a corresponding change in the effective viscosity, c ₀ in Equation(16) will remain the same for different values of n. The estimate of the effective viscosity given above may, however, not accurately describe the effect that a change of n has on n _eff. The proper integration length l in Equation (18) will, for example, depend on n if the deformation is concentrated in a different way near the bed for non-linear than for linear behaviour.

Numerical Calculations of Flow Over a Sinusoidal Bed

Solution procedure

The solution procedure applied to solve the system of differential equations defined by Equations (1), (4) and (5), subjected to the boundary conditions (7) and (8), was the method of finite elements (FE). Commercially available FE, software, the general-purpose program MARC, was used for the calculations (MARC Analysis Research Corporation, 1992). The technical details of the numerical calculations, such as mess generation, implementation of the boundary conditions and post-processing, have been explained by Gudmundsson (1994a).

Verification of Numerical Results

A comprehensive testing of the correctness of the numerical results obtained with MARC was done. A description of the testing procedure and the results obtained is found in Gud-mundsson (1994b).

Flow in parabolic channels—shape factors

As a part of the testing procedure, shape factors for parabolic channels having various aspects ratios, and for n ranging from 1 to 5, were calculated and compared to previously published results. The shape factor f is here defined as being related to the flow velocity U ₀ along the centre line according to

(33)

Since the number of calculated cases is larger and the estimated numerical errors smaller than in previous numerical studies of this problem (Nye, 1965; Echelmeyer, 1983), the results of the FE calculations are given in Table 1, as they may be of general interest.

Sliding Velocity

The sliding velocity was determined by calculating numerically the sliding function s(ε, δ, n) for n = 1 to n = 5, for r = 0.001 to r = 1.0, and for

in the range from 0.025 to 1.0. Some of the results are shown in Figure 1, which depicts In s as function of In ε for n = [1,2,3,4,5], for δ = 0.05/2π.

Fig. 1. In s as a function of In ε and n for ç = 0.05. Every symbol represents the result of one calculation. The constant slope of the curves for

shows that

in that limit.

The sliding velocity as a function of n and the amplitude-to-wavelength ratio

If

for

(see Equation (31) and (9)), then S must be independent of ε for ε small. The numerical results depicted in Fig. 1 are in an agreement with this theoretical result. The slope of the In U _b : n curves is given in Table 2. The calculated slope is −(1.017 + 0.986 n) or within 2% of the theoretical slope of −(1 + n).

Table 1. Velocities, discharge and shape factors for parabolic channels. W is the aspect ratio, defined as the ratio of the channel’s half-width ( b_p) to its depth (a_p). U_o is the velocity at the centre Line normalized by 2 a_pAT_b ⁿ and Q is the discharge normalized by 2 Aa³T_b ⁿ. The shape factor was calculated according to f = ((n + 1) U_O)l /ll. The numbers in parentheses are from Nye (1965). Calculated values are estimated to deviate less than 0.22 %from exact ones

Table 2. Linear regression coefficients for In U_b = a_n + b_n In E for the range 0< E < 0.125 and 8 = 0.05;27r. The caLcu Lated numbers Jor bn are in agreement with the theoreticaL estimate bn = - (n + 1) at the Limit E -t 0

Figure 1 shows that the range of validity of Equation (31) depends strongly on n and that this range decreases with increasing n. This is in contrast to statements made by Fowler (1981, p.675), who argued that it was only εⁿ⁺¹ that had to be small for Equation (31) to be valid. That would mean that the range of ε values, for which Equation (31) is accurate, should increase with increasing n and not decrease as the numerical results show. As an example, for n = 5, ε must be smaller than about e⁻² ≈ 0.135 (or r < 0.0215) for Equation (31) to be an approximation valid within 20%. Although a value less than 0.02 for r is not unreasonably small, this fact considerably reduces the usefulness of Equation (31). For n = 1, on the other hand, it is not until ε = 1.0 is reached that s(ε, 0.05/2π, 1) differs more than 20% from the limit.

Another interesting finding from Figure 1 is that, for ε and δ fixed, In s changes by a constant amount for every unit change in n, i.e. In (s(ε, δ, n)) — In (s(ε, δ, n + 1)) = : In Δ (ε, δ). Hence, . This is in accordance with Equation (32), which was shown to be correct for . That the function Δ(ε, δ) does not depend on n is of some practical value, since one does not therefore have to calculate s for every possible value of n. An interpolation or extrapolation based on several different n values will suffice.

In Figure 2, the logarithm of the sliding function S is shown as a function of n for a few roughness values r and a fixed corrugation value of ç = 0.005. It is seen that a straight line always results but with a slope dependent on ε. For ε = 0.01 π (the smallest value of ε in the figure)

(35)

with a standard deviation of 0.010. Equation (32) predicted this power-law behaviour but with a slope of . Note that relation (35) gives In

for n = 1, although s was defined in such a way that In s(0, 0, 1) = 0. This small discrepancy can be attributed to numerical errors and the finite values of δ and ε.

Fig. 2. In s as a function of n for several roughness values r and ç = 0.05.

For δ = 0.05/2π, S was approximated using the first six terms in Equation (11) for the range . The resulting curves are seen in Figure 3, which shows s as a function of ε ². The values of the Taylor coefficient are given in Table 3. If the numerical value of s is needed for some other n, s should be calculated using the numerical values for the first six terms in Equation (11) from Table 3 and interpolated or extrapolated using the fact that

where and are obtained through a least-squares approximation.

Table 3. Taylor coefficients of the sliding function s(ε, δ, n) = ∑_i=0 ^∞ C_2×i(δ,n) ε^2×i for δ = O.05/2π. These values can be used for δ = 0 with less than 2% error

Fig. 3. s as a function of ç² for δ = 0.05/2π and ç < π/2 (r < 0.25). The symbols represent calculated values and the lines are least-squares approximations using

. Table 3 gives the values for C _{2 × i}.

With increasing ε, the sliding-velocity variation changes markedly. The numerical results show that for large ε values (1n ε > 1),

, where r depends on n (see Fig. 4). The best straight-line approximation to γ is given by γ = −1.11 – 0.228n with a standard deviation of 0.010, i.e. for .

Fig. 4. 1n U ^b as a function of 1n ε for n = 1 to n = 5 and ç = 0.05 for large ε values (ε > 2.7). The straight lines show the best linear approximations through calculated values given by the symbols.

Comparison with Previous Estimates of the Taylor Coefficients of the Sliding Function

The values of c ₀ for n = 3 from Table 3 agree with the estimate (12) from Meyssonnier (1983) based on theoretical arguments. The value of c ₂(0, 3) = 1.163 ± 2% from Table 3 is, however, not in accordance with his estimate for c ₂(0, 3), which was c ₂(0, 3) ≈ 2.4. This deviation could be caused by the limited number of calculations done by Meyssonnier, numerical errors or both. Despite this difference, there seems to be no reason to doubt the general correctness of Meyssonnier’s findings.

The numerical results are also in agreement with Lliboutry’s (1993) theoretical upper estimate for s(0, 0, 3) (see Equation (13)) and his tabulated values for S. His upper estimate is, however, up to several times greater than the calculated values. (Lliboutry (1993) used the symbol V for the sliding function.)

The Dependency of the Sliding Velocity on Glacier Thickness

The sliding function s depends on the ratio of the glacier thickness to the wavelength of the bedrock undulations, as well as on ε and n. In Gudmundsson (1994b), the dependency of u _b on δ is discussed and it is shown that most of the variation of s with δ can be described by writing

(36)

Changing the value of δ used in Fig. 1 by a factor of 2 has almost no effect on the results shown in that figure. A relative change in ε has a much larger effect on the value of s than a corresponding change in δ. Nevertheless, by calculating the slope of the s(ε) curve for different δ values, it was found that the small but finite value δ is responsible for the 0.64% deviation of c ₀(0.05/2π, 1) from unity seen in Table 3.

The ratio of the internal deformation velocity to the sliding velocity

If one writes the surface velocity u _s as a sum of the deformation velocity u _d and the sliding velocity (which in general is only approximately correct for )

(37)

one finds

(38)

Note that Δ_b can be quite large. For r = 0.1, ç = 0.05 and n = 3, one finds, for example, using the values from Table 3, that Δ_b = 0.58. Changing r to r = 0.05 gives Δ_b = 0.99. Sliding over a hard bed without bed separation can be significant.

Basal-Flow Characteristics

As explained by Gudmundsson (1997), second-order perturbation theory for a linear medium predicts that the horizontal velocity component will have a local maximum () above the peak of the sinusoid as long as ε < 0.138, and a local minimum () above the trough of the sinusoid as long as ε < 1/2. Regions of extrusion flow develop above and below

. The region of extrusion flow above

extends up to a point called

. In the following, use of scaled coordinates, (X, Z) := k(x, z), and scaled velocities (V_X, V_z): = (v_x,v_z)/u_b will be made.

Extrusion Flow above the Peak of the Sinusoid

Figure 5 shows the percentage increase ofthe velocity

, compared to the horizontal velocity at the bed line at X = π/2, as a function of ε for n = 1 to n = 5. Symbols stand for numerically calculated values. In agreement with the analytical prediction for n = 1, no

points were found for ε > 0.138. The long dashes show the theoretical asymptotic slope for

. A good agreement between numerical and analytical results is found.

Fig. 5. The velocity increase of

with respect to the velocity at the bed at X = π/2 for ç = 0.05.

Extrusion flow is not limited to a linear medium. It is also found for n > 1. Indeed, as n increases, the range of ε values, for which extrusion flow is found, becomes progressively larger. Extrusion flow is therefore enhanced by the non-linearity of Glen’s flow law. For n = 3, extrusion flow is, for example, found for ε values of up to approximately ε = 0.2 (or r ≈ 0.03) and the velocity

is about 10% larger than the local basal sliding velocity directly below

. The vertical displacement of

for n = 1 as ε varies follows closely the theoretical prediction (Gudmundsson, 1997). With n increasing but ε and δ fixed,

moves progressively closer to the bed. This is in agreement with other numerical calculations (Raymond, unpublished).

Extrusion Flow above the Trough of the Sinusoid

As a measure ofthe magnitude of the extrusion flow, the ratio

(39)

giving the minimum of the absolute horizontal velocity at X = 3π/2 as a percentage of the velocity at the bed, can be used. This ratio is shown in Figure 6 for ç = 0.05 as a function of ε for a few different values of n.

The theoretical prediction for n = 1 was that

should only exist for ε < 1/2 (Gudmundsson, 1997). However, for n = 1, a minimum was found for ε up to 1.17 and not only up to ε = 1/2. The most direct explanation for this deviation from the theoretically predicted range is that the theoretical value is based on a perturbation analysis which is only valid for

. Since the numerical calculations show

to exist even for ε > 1, it is clear that the perturbation approach could never have given the correct answer. The asymptotic change of the velocity decrease as

, shown as a solid line, is however reproduced. There is therefore a good agreement between theory and numerics at ε values where an agreement can be expected.

For all calculated values of n, the velocity decrease (shown in Figure 6 as a negative velocity increase) becomes larger as ε increases from zero, reaches a maximum and decreases again. There is always some ε value above which no extrusion flow is found. Similarly to the situation above the peak of the sinusoidal, extrusion flow above the trough becomes progressively larger in magnitude and exists up to higher ε values as n increases. Comparison of Figure 6 with Figure 5 shows that extrusion-flow behaviour is more dominant above the trough than above the riegel. Above a trough, a 30 – 40% decrease in horizontal velocity with height over a distance of approximately 1/k is possible and could, for example, cause a considerable inversion of a borehole inclination.

Fig. 6. Relative decrease of

with respect to the velocity at the bed.

is found for ε = 0 at Z = 1 and it moves towards the bed with increasing ε (see Fig. 7). The dashed-dotted line in Figure 7 denotes Z = −ε, which is the vertical position of the bed line. As long as the symbols remain above the dashed-dotted line,

exists.

remains more or less at a constant height of approximately 1 (or at λ/2π in dimensional units) above the bed for all values of ε.

Fig. 7. Vertical position of

for n = 1 to n = 5 and for ç = 0.05.

Flow within the Trough of the Sinusoid

The velocity within an overdeepening as a fraction of the sliding velocity is often of interest. It is, for example, sometimes important to know at what roughness values ice within an overdeepening effectively remains there without taking part in the overall glacier motion. What exactly is meant by saying that the ice does not move depends, of course, on what part of the overdeepening — at the bed or only “close” to the bed — one is referring to but it turns out that this is, at least for the case of a perfectly lubricated bed, relatively unimportant.

The ratio of the local basal velocity at the base of the trough ((X,Z) = (3π/2, −ε)) to the sliding velocity u _b is one possible measure of the magnitude of the flow within the trough in relation to the mean flow along the bed line. Figure 8 shows that, as ε increases, this ratio at first becomes larger, reaches a maximum and then decreases. The somewhat surprising increase is a consequence of the extrusive nature of the flow. The maximum increase is larger for non-linear than for linear flow.

Other useful measures of the magnitude of ice movement within a trough would be the ratios of both: (1) the basal velocity at the base of the trough ((X,Z) = (3π/2, −ε), and (2) the minimum velocity at X = 3π/2 for some Z, to the basal velocity at the top of the riegel ((X,Z) = (π/2, ε)). Depicting these ratios gives essentially the same information as does Figure 8 (Gudmundsson, 1994b). For ε greater than about 1.5, there is almost no ice movement through the trough.

Fig. 8. Velocity (Χ, Ζ) = (3π/2, −ε) as a fraction of the sliding velocity (ç = 0.05).

Flow Separation

The minimum of V _x along the vertical line X = 3π/2 is shown in Figure 9. For ε large enough to exclude extrusion flow, the minimum is found at the bed (Ζ = −ε). Note that for ε > 1.8, V_X(X = 3π/2, Z = −ε) is negative, i.e. the medium in the lowest part of the trough flows in the opposite direction to the main flow. This is a clear indication of a flow separation.

Fig. 9. The ratio of the minimum of the horizontal velocity above the trough of the sine wave to the sliding velocity, as a function of ε, i.e.min_z(u(3π/2, Z))/u_b

An example of a flow separation is given in Figure 10. It is an enlarged part of Figure 11, which shows the flow above and within an overdeepening for λ = 50 m, a = 20 m, h = 200 m, n = 1 and

. Figure 11 again only shows a part of the whole configuration, which had the dimensions 200 m × 200 m. The velocities have the dimension m a⁻¹. Although the velocities within the overdeepening are small compared to the velocity at the riegel, they are a significant fraction of the sliding velocity, which for this particular case is 0.50 m a⁻¹. Figure 10 shows how the main flow induces a secondary flow circulation in a clockwise direction. A separation line is formed (shown as long dashes in Figure 10), separating the main flow from the induced flow. The ice below the separation line will theoretically circulate there for ever, never leaving the trough.

Fig. 10. Detailed view of a recirculation within a trough of a sinusoid. Only a part of the FE model is shown.

Fig. 11. Recirculation pattern within a through of a sinusoid showing a flow separation. The direction of the main flow is from left to right. The vectors indicate the direction of the flow at each FE node.

Frequency Doubling

Due to the non-linearity of the boundary conditions at the bed line, the bedrock perturbation will in general cause flow perturbations having higher harmonics than the bedrock undulation itself. A calculation of the first four harmonics of the flow field along the bed line showed the relative amplitude of higher harmonics to increase strongly with ε and the frequency doubling to be more pronounced for non-linear flow. This issue has been discussed in more detail by Gudmundsson (1994a).

Discussion

Possibly the most striking result is the onset of a flow separation within the trough of the sinusoid when the number ε exceeds a critical value. Review of the literature has not revealed any other examples of flow separation for gravity-driven flow nor for flow over a perfectly slippery boundary. The phenomenon of corner eddies in Stokes flow is nevertheless widespread (Michael and O’Neill, 1977; Hasimoio and Sano, 1980; Sherman, 1990, p.258 – 65). Corner flow, as an example, driven by circumferential motion with no-slip boundary conditions, is known to form so-called Moffatt corner eddies if the angle of the corner is less than about 146.3 ° (Moffatt, 1964).

Pozrikidis (1987) did numerical calculations of shear-driven creeping flow of Newtonian material in a channel constricted by a plain wall and a sinusoidal wall. He concluded that, for every channel width, there is a critical amplitude-to-wavelength ratio for flow separation. For wide channels, this ratio corresponds to ε ≈ 1 (see Pozrikidis, 1987, Fig. 7). This Value for the critical number is somewhat smaller than obtained here for gravity-driven flow. One should not expect perfect agreement, since the driving mechanism of the flow is not the same and because the local properties of Stokes flow depend in general strongly on the global structure of the flow (Sherman, 1990). Note that, due to the limited spatial resolution of the FE mesh, the possibility that the critical value for flow separation is somewhat less than 1.8 cannot be ruled out.

Glaciologists often tend to think about basal flow of glacier ice in terms of the simple analytical solution for a plane slab flowing down an inclined plane. The presence of even small basal undulation leads, however, to a different picture of basal flow than the plane-slab solution suggests. The horizontal velocity can increase with depth and even at only moderate amplitude-to-wavelength ratios of 0.28 a flow reversal takes place. It should therefore not come as a surprise if the stratigraphy of basal ice becomes strongly disrupted by glacial flow. An ε value of 1.8 can hardly be considered to be unrealistically large, and it must, in general, be expected that glacier beds have regions where the ε values are this big.

Where an overdeepening is found, with a depth-to-width ratio corresponding to an ε value of about 1.8, the ice will most probably move directly over the trough and ice within will be stagnant. A possible candidate for this type of flow regime is the spectacular overdeepening found at Konkordiaplatz, Aletschgletscher, Switzerland, which is about 1000 m long and 400 m deep.

Conclusions

With the use of both analytical and numerical methods, the flow characteristics and the sliding velocity of a highly viscous medium flowing under the influence of gravity over a perfectly lubricated sinusoidal bed have been analysed.

Directly above the peak of the sine curve (kx = π/2), a local maximum of the horizontal velocity component develops if is listed as a function of n for in Table 4 for n = 1 to n = 5. The value of increases with increasing n, showing that the non-linearity of the flow law makes the range of ε values for which

exists larger.

Table 4.

as functions of n for δ ﹤﹤ 1. The value of

is based on an analytical solution (Gudmundsson, 1997). All other values are based on numerical calculations

Above the trough of the sine curve (

), a local minimum of the horizontal velocity component develops if

. The calculated values of

for

are also given in Table 4.

Regions of local extrusion flow are associated with both these stationary points (

). The non-linearity of the flow law increases the range ε values for which

is found, as the table shows. The velocity decrease with respect to the velocity at the bed also becomes more pronounced with increasing n (see Fig. 6).

Flow separation occurs for ε > 1.8 for at least 1 ≤ n ≤ 5, for both perfectly lubricated beds and for no-slip boundary conditions. Based on the large number of analytical, experimental and computational demonstrations of its existence, flow separation for no-slip boundary conditions is known to be a universal feature of laminar flow in corners (Sherman, 1990, p. 265).

The sliding velocity as a function of n, ε and δ has been calculated. The numerical results agree with Equation (31). For finite values of ε, the sliding velocity can be determined by using the Taylor coefficients listed in Table 3.

In u _b depends linearly on n for all values of ε and δ as Equation (32) suggests. Using this fact, it is possible to calculate s as a function of n for all values of n by using Table 3 to calculate the sliding function for several n values and then interpolate or extrapolate the results assuming 1n

for some

and

.