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

is the second invariant of the deviatoric stress tensor

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

and

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

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 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) : = 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.

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 ₀ 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

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 + (dz₀/dx)²))^−1/2 gives for the drag (the x component of Equation (26))

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

giving

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

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

It follows from the first-order perturbation solution for the pressure that c ₀ = 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 ₀ 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.

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).

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 Reference Nye,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 Reference Fowler,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)

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 Reference Meyssonnier,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 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.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.

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 (Reference Gudmundsson,Gudmundsson, 1997). With n increasing but ε and δ fixed, moves progressively closer to the bed. This is in agreement with other numerical calculations (Reference Raymond,Raymond, unpublished).

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

should only exist for ε < 1/2 (Reference Gudmundsson,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 (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, 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 Reference Gudmundsson,Gudmundsson (1994a).