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 (Reference Nye,Nye, 1969, Reference Nye,1970; Reference Kamb,Kamb, 1970; Reference Morland,Morland, 1976a, Reference Morland,b; Reference Gudmundsson,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 (Reference Gudmundsson,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
where 
are strain rates and 
are the deviatoric stresses defined by
and δ ij is the Kronecker delta. The field equations for the conservation of mass and momentum are
and
where vi 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
where Z 0 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 Reference Gudmundsson,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 (Reference Gudmundsson,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
and
where tan β(x) : = dz0(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 Reference Gudmundsson,Gudmundsson (1997), where references to earlier work can be found.
Some knowledge of the sliding velocity, u b, can be obtained from dimensional analysis (Reference Lliboutry,Lliboutry, 1968, Reference Lliboutry,1987b). There are three dimensionless numbers that enter the problem: n, ε and δ. For Glen’s flow law, dimensional arguments give
where U b, is a non-dimensional sliding velocity.
By using a variational theorem for non-Newtonian flow (Reference Johnson,Johnson, 1960), Reference Fowler,Fowler (1981) derived a sliding law valid for a non-linear medium. He concluded that the sliding velocity u
b should vary as 1/ε
n+1 in the limit 
. This has, however, been questioned by Reference Schweizer,Schweizer (1989) and Reference Schweizer, and IkenSchweizer and Iken (1992) who, through finite-element (FE) calculations, found u
b to less dependent on ε than 
indicates.
Reference Raymond,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 ε−(n+1), then it is convenient to define a function s(ε, δ ,n), called the sliding function, as
s can be thought of as a non-dimensional sliding velocity where the ε−(n+1) dependency has been accounted for. Linear perturbation theory (Reference Nye,Nye, 1969; Reference Kamb,Kamb, 1970) gives 
.
The function s(ε, δ, n) can be expressed as a Taylor series with respect to the variable ε
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. Reference Fowler,Fowler (1981), Reference Lliboutry,Lliboutry (1987a) and Reference Meyssonnier,Meyssonnier (1983) were able to give a numerical estimate of C 0(0, 3). Using variational methods, Meyssonnier found
Reference Meyssonnier,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 2(0, 3), which was c 2(0, 3) ≈ 2.4.
Reference Lliboutry.Lliboutry (1993) later improved Meyssonnier’s estimate, again using variational methods, and obtained an upper bound on s for n = 3. His result is
Reference Kamb,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 (Reference Kamb,Kamb, 1970, p. 703)
Using this expression, one finds that c
0(0,3) = 4/e2 ≈ 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/ε(n+1) behaviour of u
b for 
.
Reference Schweizer,Schweizer (1989), in a numerical study, came to the conclusion that the sliding velocity does not vary as ε−(n+1) 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 Reference Lliboutry,Lliboutry (1968, Reference Lliboutry,1987b) and Reference Kamb,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 Reference Fowler,Fowler (1986, Reference Fowler,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
where
C
0 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 ηubk 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/ub, it is stretched/compressed to the length l + Δz, where Δz = Δxdz0/dx. The (average) strain rate is given by
where the brackets are used to indicate that these are averaged strain rates over a given vertical distance l, that is
For the average vertical strain rate, one gets
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 
The effective strain rate έ is
leading to
Glen’s flow law can be written as
where
Inserting Equations (22) and (24) into Equation (16) gives
Force equilibrium requires
where γ is a path along the sine curve, C(γ) is the path length, σ is the stress tensor and 
is a unit normal vector
A change of varibles dx = ds (1 + (dz0/dx)2))−1/2 gives for the drag (the x component of Equation (26))
where tan β = dz0/dx. The last term on the righthand side is of second order in ε and the exact boundary condition on z
0, Equation (8), shows that 
is then also of second order. Up to first order in ε, the drag can therefore be calculated according to
giving
where c 0 is the same unknown constant introduced in Equation (16). After integrating and re-arranging terms, one gets
Expression (31) brings out the asymptotic 1/εn+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
It follows from the first-order perturbation solution for the pressure that c 0 = 2 for a linear Newtonian medium (Reference Nye,Nye, 1969, Reference Nye,1970; Reference Kamb,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 0 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 Reference Gudmundsson,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 0 along the centre line according to
Since the number of calculated cases is larger and the estimated numerical errors smaller than in previous numerical studies of this problem (Reference Nye,Nye, 1965; Reference Echelmeyer,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
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
The sliding velocity as a function of n and the amplitude-to-wavelength ratio
If
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 ( bp) to its depth (ap). Uo is the velocity at the centre Line normalized by 2 apATb n and Q is the discharge normalized by 2 Aa3Tb n. The shape factor was calculated according to f = ((n + 1) UO)l /ll. The numbers in parentheses are from Reference Nye,Nye (1965). Calculated values are estimated to deviate less than 0.22 %from exact ones
Table 2. Linear regression coefficients for In Ub = an + bn 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 Reference Fowler,Fowler (1981, p.675), who argued that it was only εn+1 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−2 ≈ 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)
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
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 ε
2. 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
and 
are obtained through a least-squares approximation.
Table 3. Taylor coefficients of the sliding function s(ε, δ, n) = ∑i=0 ∞ C2×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 ç2 for δ = 0.05/2π and ç < π/2 (r < 0.25). The symbols represent calculated values and the lines are least-squares approximations using
With increasing ε, the sliding-velocity variation changes markedly. The numerical results show that for large ε values (1n ε > 1),
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 0 for n = 3 from Table 3 agree with the estimate (12) from Reference Meyssonnier,Meyssonnier (1983) based on theoretical arguments. The value of c 2(0, 3) = 1.163 ± 2% from Table 3 is, however, not in accordance with his estimate for c 2(0, 3), which was c 2(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 Reference Lliboutry.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. (Reference Lliboutry.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 Reference Gudmundsson,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
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(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 
)
one finds
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 Reference Gudmundsson,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
Extrusion Flow above the Peak of the Sinusoid
Figure 5 shows the percentage increase ofthe velocity
Fig. 5. The velocity increase of
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
Extrusion Flow above the Trough of the Sinusoid
As a measure ofthe magnitude of the extrusion flow, the ratio
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
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
Fig. 7. Vertical position of
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 (Reference Gudmundsson,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, VX(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.minz(u(3π/2, Z))/ub
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
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 Reference Gudmundsson,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 (Reference Michael, and O’Neill.Michael and O’Neill, 1977; Reference Hasimoio, and SanoHasimoio and Sano, 1980; Reference Sherman,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 ° (Reference Moffat,Moffatt, 1964).
Reference Pozrikidis,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 Reference Pozrikidis,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 (Reference Sherman,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
Table 4.
Above the trough of the sine curve (
Regions of local extrusion flow are associated with both these stationary points (
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 (Reference Sherman,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
Acknowledgements
This research was supported by Swiss Science Foundation grant 20-29619.90. I am grateful to A. Iken, K. Hutter. C. F. Raymond and D. Bahr for thorough comments, which helped to improve an earlier version of this manuscript.
