1. Introduction
The distribution of stress and velocities in glaciers and ice sheets has been treated previously by several authors. The effect of longitudinal strain on the state of stress at the base of a glacier or ice sheet was studied by Reference RobinRobin (1967) and Reference BuddBudd (1969) has found its theoretical basis in the articles of Reference CollinsCollins (1968), Reference NyeNye (1969), and Reference BuddBudd (1970[a]). In these articles attention is focused on the significance of the longitudinal stresses set up by the flow of ice over protuberances on the bed. The idea is to improve the basal shear-stress formula, which, to zeroth order, is independent of the material behaviour of ice and only involves the local glacier depth and the local inclination of its surface.
Steady-state stress and surface variations in infinite ice slabs with small undulations of their beds have equally been treated, notably by Reference Yosida and YoshidaYosida (1964) and Reference BuddBudd (1970[b]) who also gives further information on the pertinent literature. Yosida considered Newtonian viscous flow with no slip of a medium of uniform thickness down a uniform slope with small harmonic undulations superimposed on it. Budd generalized this solution by allowing perfect slip along the bed and by using an approximate solution technique to account for a non-linear viscous behaviour. His method gave answers to the question of how small bedrock irregularities are transferred to the surface topography and in what way they influence the stress distribution. In particular Reference BuddBudd (1970[b]) was able to give precise statements regarding the significance of the longitudinal stresses set up by the flow. He could explain under what circumstances the basal shear stress fluctuates in sympathy with the surface slope.
Budd’s careful analysis requires the knowledge of a mean velocity across the depth of the ice slab. More precisely, in his work the velocity boundary conditions are determined by the shape of the base contour, the mean down-slope velocity, and the equation of continuity formulated as an integrated mass balance. This equation contains the accumulation rate, so that, in principle, the effect of the accumulation rate on the topography could be determined. Moreover, for a time-dependent accumulation rate, one would expect the surface undulations to be a direct result of the accumulation rate rather than the bedrock topography. The corresponding time-dependence of the surface topography will, under certain circumstances, result in surface waves travelling down the glacier akin to those treated by Reference NyeNye (1960, Reference Nye1963[a], Reference Nye[b]) in his kinematic-wave theory.
It is possible, if one so desires, to challenge the above work by, first, extending the analysis to include steady-state and transient accumulation rates and, secondly, to improve it by replacing the velocity boundary condition by a kinematic condition at the top surface. For it is quite clear that accumulation and ablation should enter the problem through a condition regarding the time-evolution of the glacier surface rather than an integrated mass balance. Moreover, the boundary condition at the ice-rock interface treated by Yosida as a no-slip condition and by Budd as a perfect slip could be generalized to account for velocity dependent friction. In other words, the determination of stress and velocities in glaciers and ice sheets and the propagation of surface waves should be unified and be derivable from a single general concept.
There exists a rational approach which delivers this unification, and explicit calculations can relatively easily be performed for an infinitely long slab with almost parallel sides. The method roughly proceeds as follows: As a first step, the general governing equations consisting of the balance laws of mass and momentum, the constitutive relations (Glen’s flow law in Nye’s three-dimensional extension), and the boundary conditions are derived. The latter include a statement about the sliding mechanism at the base, they express the stress tensor at the top surface in terms of the atmospheric conditions and involve a statement about the rate of change with time of the surface geometry in terms of the accumulation rate. These equations are perfectly general and no assumptions need to be invoked regarding the geometry of the ice sheet. The solution of a boundary-value problem is rather complex at this generality, however, due to two circumstances, first, Glen’s law is non-linear, and, secondly, the geometry of the top surface is part of the solution of the boundary-value problem to be solved. The second step in attacking the complex problem is thus a restriction to simple geometries and simple flow configurations, or at least to patterns close to such situations. A realistic model for glaciers and ice sheets is the parallel-sided infinite slab of which the flow pattern and stress distribution was determined by Reference NyeNye (1953, Reference Nye1957). Problems, to which its steady-state behaviour can be considered as a first approximation can, among others, be formulated as follows:
-
(i) What is the influence of small bedrock undulations on the stress and velocity-distribution and on the surface topography?
-
(ii) How does steady-state position-dependent accumulation rate influence the stress distribution and the surface geometry and, can this effect be separated from that of the bedrock undulations?
-
(iii) What are the characteristic features of surface waves? In other words, does the kinematic wave theory suffice for the description of these waves, or does it have to be extended?
-
(iv) Would such an extended wave theory perhaps give rise to a possible explanation of seasonal and surge-type waves?
Of these questions only question (i) has been sufficiently treated; no systematic study of the effect of steady accumulation rate is known to us, and the transient response is uninvestigated except for some remarks made by Reference BuddBudd (1970[b]). Moreover, since Budd’s solution is based on an approximate integration technique, and because he does not seem to treat the influence of the inclination of the ice slab and of the sliding mechanism systematically, it is justified to reinvestigate item (i) at least from this restricted point of view.
The aim here is an attempt to provide answers to the above-stated questions. In the present paper we shall derive the basic governing equations, but shall deal in detail only with the steady-state response. The transient behaviour and, in particular, surface wave problems are reserved to another paper (Reference HutterHutter, 1980) as is a direct approach to the longitudinal strain effects on basal shear (Reference HutterHutter, 1981). Hence, only items (i) and (ii) above will be treated here. Attention will exclusively be focused on the ncarly-parallel-sided slab. This restriction is tantamount to the neglect of effects of accumulation on the geometry of the ice slab, which therefore limits consideration to length scales which are not asymptotically large as compared to glacier thicknesses. The results of the present paper, however, are meaningful on length scales over which ice thicknesses do not change considerably.
The rationale in achieving an analysis of the nearly-parallel-sided slab is, first to establish field equations and boundary conditions for stress and velocity in a suitably non-dimensional form, secondly to replace the boundary-value problem as stated for the true geometry by an approximate one of the strictly parallel-sided slab and, thirdly to solve this approximate boundary-value problem by a regular perturbation approach. The non-dimensionalizadon is different from the usual ones in liquid-film theories; contrary to the latter, it is valid for all inclinations γ of the slab, including γ = 0. In the application of the perturbation procedure we shall, however, restrict considerations to values of γ which are bounded away from zero, but nevertheless indicate how calculations should proceed for small γ or γ = 0.
2. Basic governing equations
The model of a glacier or ice sheet is depicted in Figure 1. It shows a vertical cut through an infinitely wide glacier; the upper surface and bed are thus cylindrical. In the mean these surfaces are assumed to be plane over most parts of the glacier. Apart from the glacier head and the region at the snout, the bottom and top surfaces are thus approximately parallel, or for a distance that is large compared with the ice thickness may be regarded as approximately parallel. The Cartesian coordinate system is taken such that the x-axis is parallel to the mean direction of the top or bottom surface and the local depth is measured perpendicular to the x-axis. Both top and bottom surface are inclined with respect to the x-direction so that, in general α and β do not vanish. Yet, we shall assume that the deviation of the glacier surface from its mean γ = D and of the bottom topography from the line γ = 0 are small. Denoting the distances of these surfaces from the x-axis byγ s and γ b respectively, we thus may assume |γ b| and |γ s–D| to be small in comparison with the glacier thickness.
Fig. 1. Nearly parallel-sided slab. Definition of geometry.
We restrict considerations to plane strain and presume the ice to be isotropic, but generally we shall not strictly assume it to obey Nye’s extension of Glen’s flow law. The reason is the singular behaviour of Glen’s flow law at small stress deviators, which leads in certain problems to mathematical singularities, which can be removed by a realistically modified constitutive relationship. As for Glen’s flow law this constitutive relationship is an incompressible non-linear viscous fluid model. We do not presume that temperature effects are negligible, but allow for a variation of temperatures perpendicular to the main flow. The glacier may thus be cold or temperate. For most of the analysis the latter will be assumed.
2.1. On an extension of Glen’s flow law
Glen’s power law carries the disadvantage that stress deviators grow infinitely fast at small stretchings. In some of the numerical solutions to be presented below this turns out to lead to singular behaviour. Rather than deal with these singularities mathematically, we shall, when necessary, replace Glen’s flow law by another one that does not show this singular behaviour. The idea is not new and has independently been presented by Reference HutterHutter (1979) and Reference ThompsonThompson (1979). To introduce this more general constitutive relation, let d be the stretching tensor and t’ the stress deviator. The class of constitutive relationships given by
in which A and n are constants and B is a function of the second stress deviator invariant t 11 2 = ½t’ : t’, embraces Glen’s flow law as well as its extension
The constant k introduces a nearly linear Newtonian behaviour at small stretchings and removes the singularity mentioned above. For brevity we shall call Equation 1 with the specification (2) the generalized Glen flow law. A and n are known from the usual law and k can be determined from creep experiments at low strain-rate. Glen’s flow law is obtained for k ≡ 0.Footnote *
Below we shall introduce a non-dimensionalized version of the function B, denoted by
For the generalized Glen flow law this function reads
To find a numerical value for
In the subsequent calculations we shall therefore choose n = (2, 3) and correspondingly
In cold ice the temperature varies very nearly with a distinct coordinate, say y; such a dependence can be approximately taken into account by assuming A to be y- dependent. For such a case the function
A(y) is obtained, if in the Arrhenius law
Before proceeding a word of caution seems to be in order. In what follows, perturbation solutions for small bedrock undulations with amplitude є will be sought, and these perturbation solutions will be constructed considering є → 0, keeping
2.2. Field equations
In order to simplify the subsequent analysis, the governing field equations will be written in dimensionless form. To this end a characteristic length D and a characteristic time τ will be introduced. Stresses are then non-dimensionalized with the reference stress ρgD and velocities with V = D/τ. With these scales and with the plane-strain assumption, the balance laws of mass and momentum and the constitutive relations assume the form
where an explicit dependence of
is the second invariant of the stress deviator. Moreover, ℱ and
whose values depend on the selection of D and τ (or U). Choosing for D the mean glacier thickness guarantees that the dimensionless stresses vary between —1 and 1. It is then straightforward to see that ℱ is very small (< 10−6) for even unrealistically large values of the velocity U. This justifies the neglect of the local and convective acceleration terms and allows for a relatively free choice of the characteristic velocity U or time τ. Three different choices are directly suggested by the problem, namely,
-
(i) to simply assign a specific value to
(= 1);
-
(ii) to choose U according to realistic longitudinal surface velocities;
-
(iii) to scale velocities with the accumulation rate.
The first choice is mathematically convenient, the second and third are physically suggested. For case (ii) emphasis is laid on a proper treatment of velocity profiles and accumulation is regarded as a secondary effect. In the third case, accumulation is regarded as important. Realistic values for surface velocities and accumulation rate are, perhaps, 100 m a−1 and 1 m a−1 respectively. If we then choose D = 100, a = 1 m/a, n = 3, A = 10−19 dyn −n cm2n a−1 the orders of magnitude shown in Table 1 are obtained. If one scales equations with the choice (ii) one guarantees that the dimensionless velocity u and the stresses (σ x , σ y ) are of order one. The dimensional quantities are then obtained from the non-dimensional ones according to
and realistic values for these are
Table I. Orders of magnitudes for variouschoices of the characteristic velocity U or time τ
In what follows Equations (3), which form a set of five partial differential equations for five unknowns, will be solved neglecting inertial terms ℱ → 0. The emerging system then contains only one dimensionless number, namely
Frequently it is advantageous to work with the stress deviator rather than with stresses. Denoting the former by primes these are related by
Introducing the dimensionless pressure p as an unknown, the momentum equations read in this case
When written in terms of the stress deviators the field equations (3) comprise six equations for six unknowns, the new unknown being the pressure.
Finally, we mention that there are still other possibilities for the non-dimensionalizations. For instance, one may non-dimensionalize stress using a reference stress ρgD and may introduce a dimensionless pressure p according to
Denoting dimensionless quantities by overhead bars, it can be shown that the field equations (2) become in this case
where
This is the non-dimensionalization used in the theory of viscous liquid films. The system of Equations 9 has the advantage of being better adjusted to the fact that the flow is gravity driven. Certain spurious singularities which occur for small γ with Equations 3 can be shown to be removable when using Equations 9 yet Equations 9 are definitely improper for γ = 0, whereas Equations 3 remain valid. The subtleties shall be pointed out at the appropriate places.
2.3. Boundary conditions
To complete the formulation of the problem, the field equations must be complemented by boundary conditions. At the base the latter will be expressed in either one of the following two ways. At first, we assume perfect slip with no cavity formation and thus may write
Here ū is assumed to be a prescribed sliding velocity, which as we shall see must be independent of time and space, see below. It is not the same as the tangential speed, but is proportional to it. For ū = 0 Equation (11) includes the no-slip condition.
There is another possibility of formulating the boundary conditions at the base of a glacier without cavity formation. In this case one sets
The first of these is the tangency condition. The second, on the other hand expresses a velocity-dependent friction law akin to that of Reference WeertmanWeertman (1957). For
At the top surface we require that stress is continuous and that the conditions of a kinematic surface are satisfied. As regards the latter, it should be pointed out that the free surface is not a material surface, since due to accumulation and ablation particles are added or removed. If F: = Y s(x, t)—y = 0 is the equation of the surface we must thus have
where a(x, t) is the dimensionless accumulation rate, which is positive as an accumulation rate and negative as an ablation.
The boundary conditions of stress on the other hand read
Here α is the surface inclination so that
The minus sign stems from the definition of positive angles α; p is the dimensionless atmospheric pressure, which, in general, is a function of x and t. We shall neglect this variation and assume p to be constant.
With the inclusion of the boundary conditions at the base and at the top surface a well-set boundary-value problem has been obtained, in which γ s (x, t) must also be considered as unknown. The equation for its determination is Equation (13).
Before proceeding further it is worthwhile to give an order of magnitude for the accumulation rate. For this purpose note that
With a dim = 10 m a−1, D and τ as before, a value of a ≤ 10−3 to 10−4 is obtained. Hence the dimensionless accumulation rate is very small, and to zeroth order it may safely be neglected. Clearly, these values are hinged to the choice of the characteristic time; it precludes analysis of accumulation effects on glacier geometry. Moreover, the friction constant in the first of Equations (12) is related to its dimensional counterpart by
and is an appropriately chosen constant.
2.4. Decoupling of the governing equations into steady-state and transient problems
One of the major questions raised in this article is the same as that asked and already answered to a high degree of satisfaction by Reference BuddBudd (1968, Reference Budd1969, Reference Budd1970[a], Reference Budd[b]) namely to evaluate the influence of the bedrock undulations on the surface topography. Another is to estimate the effect of the accumulation as a function of position or time or both. This second problem will lead to surface waves travelling down the glacier, akin to the kinematic waves studied by Reference NyeNye (1960, Reference Nye1963 [a], Reference Nye[b]).
Below, the separation of the total fields into steady-state and time-dependent fields will be demonstrated although ultimately in this paper only the steady-state problem will be treated, because this paper will serve as a reference for others (Reference HutterHutter, 1980, Reference Hutter1981).
To separate the time-dependent from the time-independent problem, let
be the decomposition of the total fields into steady-state (circumflexes) fields and quantities which are deviations from this state (tildes). By introducing the separations 18 in the field equations and boundary conditions, boundary-value problems for the steady-state and the purely transient part of the problem can be obtained. This decomposition being straightforward, we shall give below the resulting equations only.
(a) Steady-state problem
Formally, the steady-state equations are obtained from Equations (3), (11), (12), (13), and (14) by merely dropping local time derivatives and using circumflexes. This essentially amounts to formal changes in the kinematic boundary condition at the top surface, which now becomes
and in appropriate changes of the momentum equations, if R x and R y should be kept.
For brevity the equations will not be repeated here.
(b) Transient problem
For the solution of the transient problem the solution to the steady problem is assumed to be known. In the transient problem we then only search for the deviations from the total fields. This process of separation is straightforward in principle, but somewhat messy if no resort is made to a simplifying assumption. We shall do this and postulate that the steady-state stresses are much larger than the corrections due to a possible time-dependent effect. Based on this assumption the constitutive relations in the fourth and fifth of Equations 3 can be linearized. This linearization in the stresses cannot be over-emphasized. In particular, it does not imply that ũ and ṽ need be small. A linearization in these quantities might indeed be inappropriate.
The basic governing equations are:
in which
are a steady-state and a transient second stress-deviator invariant. On the other hand, the boundary conditions become:
At the base:
The formulae on the left apply when the boundary condition is of kinematic nature, those on the right when a velocity-dependent sliding condition is used;
At the top surface:
for the kinematic condition
and for the boundary condition of stress
In these equations α is the total angle of surface inclination,
Equations 24 complete the formulation of the problem. In the following sections a systematic approach will be shown which allows us to attack various problems occurring in the velocity and stress distribution of glaciers and large ice sheets.
3. Steady-state solution for the nearly parallel-sided slab
The equations of the steady-state problem will now be subject to an approximate solution scheme. The key idea in this is the assumption that the glacier surface and the bottom topography do not deviate considerably from the infinite or semi-infinite parallel-sided slab. On this presumption an iterative scheme will be developed which allows us, first, to calculate the influence of the bottom topography on the surface geometry and, secondly, to estimate the influence of a steady accumulation or ablation rate on the stress and velocity distribution as well as on surface topography. Some of these questions have already been answered by Reference BuddBudd (1970[a], Reference Budd[b]), however his solution is approximate, and here we aim at a justification or disproof of his approach.
3.1. Steady-state behaviour of an infinite slab on a bed with small undulations
(a) Perturbation scheme
The fundamental assumption on which the following perturbation approach is based is that amplitudes of the base are small in comparison to the mean glacier thickness. Consequently, we may write
In this equation Λ(x) and its derivative Λ′(x) are O(1). The representation 25 suggests searching for a perturbation solution in the form
in which
If the representations 26 are substituted into the field equations and boundary conditions of the steady-state problem and if, in the respective equations, terms of equal powers in є are collected, a hierarchy of boundary value problems is obtained. For instance, the zeroth- and first-order field equations are as follows:
-
(i) zeroth-order field equations
(27) -
(ii) first-order field equations
(28)
Higher-order equations could also easily be deduced. They are of little interest, however.
The derivation of the approximate boundary conditions is based on the idea that the condition valid on the boundary surface is replaced by an appropriate condition on the mean bed and mean surface respectively. To explain the procedure, consider Equation (11), in which ū and y B(x)must be regarded as known functions of x; û and û, on the other hand are functions of the form
Developing these functions in terms of Taylor series expansions of the second variable and substituting the resulting series into the equations results in new forms of the boundary conditions in which all functions are evaluated at y = 0, rather than at the actual bed. Hence we may say that the simple relations (11) on the complex geometry yB(x)are replaced by formally more complex equations on a much simpler geometry, namely y = 0. The complexity of the new equations is not really increased as compared to Equations (11), however, because the small parameter є will enter and allow us to make use of the expansions (26). If this is done the kinematic boundary conditions at the base read, to first order,
In much the same way the conditions (12) can also be handled. To first order this gives
Again, higher-order equations are easy to obtain. In what follows we shall either use Equations (29) or (30).
Next, the boundary conditions at the glacier surface are transformed to appropriate conditions at a mean surface, y mean =
In view of Equations 26 this can be transformed into a hierarchy of kinematic boundary conditions. Because a is of the order of 10−2 to 10−4 it is also justified to write
so that A is of order 1 at most. With the aid of Equations (32) and 26, the above boundary conditions to first order thus become
Next, the boundary conditions of stress, Equations (14), in which all quantities carry a circumflex, must be expressed as a condition on the mean surface. To this end, notice that all quantities on the left are evaluated at the surface, while p, on the right is a given function of x. Proceeding as before we may replace the first of Equations (14) to zeroth order by
and to first order by
Summarizing, Equations (27), the first of Equations (29) or (30) and of (33), and Equations (34) comprise together the zeroth-order boundary-value problem. The solution to this problem is assumed to be readily available. Knowing it, the solution to the first-order boundary-value problem can then be constructed using Equations (28), the second of Equations (29) or (30) and of (33), and (35). This boundary-value problem is linear.
(b) Zeroth-order solution
Consider a doubly infinite strip of thickness
in which the boundary conditions (34) have already been incorporated. The solution in Equation (36) is independent of the constitutive response. Notice also that
In view of the third of Equations 27, ū must be a constant.
The solution in Equation 37 is based on a purely kinematical boundary condition. If, at the base use is made of the first of Equations (30) instead, ū in Equation (37) must be replaced by
Hence, to zeroth order, the two types of boundary conditions at the base are equivalent.
(c) First-order solution
To construct a solution of the first-order boundary-value problem, it is advantageous to introduce stress and stream functions. Accordingly, we set
thereby satisfying the first three of Equations 28 identically. The remaining equations then transform into the system
Equations 41 are a second-order system of linear partial differential equations for Φ, and ψ with y-dependent coefficients. For a Newtonian fluid, f(y) = ½, g(y) = 2, so the coefficients are constant in this case. Note that the function
The system (41) must be complemented by appropriate boundary conditions. These are at the base: if Equations (29) are used
and if Equations 30 are used instead,
At the top surface the second of Equations (33) and Equations (35) must hold. When written in terms of stream and stress functions, these conditions read
where
Before we attempt to find a solution to the above problem we would like to emphasize that the difference in the boundary conditions at the base, Equations (43) and (44), is now explicit and, perhaps, also visible after integration. Note further, that the surface topography enters the formulation and must be determined along with the determination of the entire problem. The explicit occurrence of the accumulation rate complicates this problem slightly. It becomes easier if one sets
This choice transforms Equation (45) into
in which
Before we proceed, we indicate how the general steady-state problem with steady-state accumulation rate can be solved. To this end, it is not hard to see that the influence of the bedrock undulations and that of a steady accumulation rate can be fully separated. Indeed let
(d) Approximate integration scheme
The differential equations (41) are not easy to integrate, in general; it is thus appropriate to reserve a separate section for this integration. Here we would like to demonstrate, how the first-order equations can be integrated approximately. Such an approximate scheme was introduced by Reference BuddBudd (1970[a]) without any mention of its approximate character. Because of the prominence of his article and our conclusions in this article as to the falsity of the approximate integration scheme, it will be briefly outlined here. To this end we write Equations (28) in terms of deviator stresses
From the continuity equation and the constitutive relations involving
whereas the constitutive relation for
where ψ’ is the stream function of this approximate problem (primes will indicate that the perturbation pressures
The differential equations which must be solved in this approximate approach are Equations (50) and (51). As compared to Equations (41) they carry the advantage of being separated insofar as the general solution of the Laplace equation for the stress function can be constructed in advance and may, in a second step, then be used as a forcing function in the wave equation (51) for the stream function ψ’. Of course, the suitability of the approximate boundary-value problem in which
3.2. Harmonic perturbations from uniform flow for zero accumulation rate
In this section we shall solve the boundary-value problem for the special case that A(x) = 0. We assume that Λ(x) is given as a Fourier series,
It suffices only to look at a single Frequency, and so we will simply write
We also expect the surface topography to be sinusoidal, but shifted with respect to Equation (53), so that
The functions
will be called the filter function and the phase shift, respectively. The former determines that fraction of the protuberance amplitude that is, at a certain frequency, carried over to the surface. It could also be interpreted as the ratio of the surface to the basal amplitude of the inclination. According to the observations and the calculations of Reference BuddBudd (1970[a]) this function should have a maximum for wavelengths which are between 3 to 4 times the glacier thickness. The phase shift, on the other hand, is close to π/2.
Because the Equations (50) and (51) will be shown to lead to false results and the approximate solution will only be used to demonstrate this point, the derivation of the solution will not be presented here. Interested readers may consult the first author (Hutter). Here we only list the relationships for h 1 and h 2 which will allow the evaluation of the filter and the phase-shift functions when Glen’s flow law is used. One obtains
in which
V p(y) and V m(y) are the functions
where
With these formulas the evaluation of the filter function ℱ becomes a routine matter. We shall postpone the presentation of numerical results until the next chapter.
(b) Exact solution for a Navier–Stokes fluid
Unfortunately, an exact integration of Equations (41) is hardly possible for general non-linear material behaviour. This case, although being of a somewhat academic nature, is nevertheless important because it allows a check of the numerical scheme needed for the non-linear fluid. Since for this case f(y) = ½ and g(y) = 2 the coefficients in Equation (41) are constant. With
the emerging differential equations for the coefficient functions assume the form
The boundary conditions (43) and the second and third of (45) on the other hand become: at y = 0,
at y = 1,
The form of Equations (65) suggests searching for exponential solutions. Indeed, the reader may check himself that
is the most general solution of Equations (65). The eight unknown coefficients are determinable from the eight boundary conditions (66) and (67). The corresponding calculations are very tedious, even though they are straightforward, so that it does not seem to be appropriate to present them explicitly. After lengthy manipulations one obtains:
and
where
and
With the stress and stream functions being determined as functions of the surface and bottom amplitudes h 1, h 2 and B = 1, the former two follow from Equations (48) by mere substitution. If this is done, we obtain
The filter function ℱ and the phase lag angle φ are therefore given by
This completes the construction of the exact solution.
(c) Numerical solution for non-linear rheology
For the harmonic bed undulations (Equation (53)) the stream and stress functions may be assumed to have the form of Equations (64). When substituted into Equations (41), the following system of ordinary differential equations for the unknown functions Φ1, Φ2, ψ1 and ψ2 are obtained:
These must satisfy the boundary conditions (43) and (48), which with the aid of Equations (54) and (64) may be written
and
Introducing the vector
it can readily be shown that Equations (76) correspond to the linear vector differential equation
where A is given by
The boundary conditions (77) and (78), on the other hand, can be written as follows: At y = 0
and at y = 1:
To solve the boundary-value problem (80), (81), and (82), it is best to integrate it numerically with, say, a Runge–Kutta scheme. Because of the linearity of the problem the principle of superposition can be used. Accordingly, we solve Equation (80) with the following five different initial conditions:
For each of these initial conditions the values of g 1,…, g 4 in Equations (82) can be evaluated; they are
The correct solution is obtained if
which is a linear system of four equations for the unknowns X j. Once they are determined we may integrate Equation (80) once more with the initial condition
to obtain the final solution for the stress and stream functions. First-order stresses and velocities are then obtained by simple differentiations. A program was written solving the boundary-value problem in the indicated manner. Solutions of it can be obtained for various different physical situations, some of which we shall discuss below.
3.3. Effect of steady accumulation rate
In the last subsection accumulation was set to zero, but bottom undulations were present. Here we shall assume a non-vanishing accumulation rate, but shall set all bottom undulations aside. If both are small their effects can be superimposed, as we have already seen above.
It was mentioned earlier that the zeroth-order solution constructed in the Equations (36) and (37) did not account for geometric effects of steady-state accumulation rate. This limits considerations to length scales which are not asymptotically large as compared to glacier thicknesses. The effect of accumulation rate on length scales somewhat larger than glacier thickness may, however, still be analysed.
To solve the boundary-value problem described by Equations (41), (43), and (48) we shall set Λ ≡ 0 and, further assume that
Stress and stream functions may then be given as shown in (64) and
In these equations f is defined in Equation (79), A(y) in (80b) and g 1, …, g 4 in (82). The construction of numerical solutions to the above boundary-value problem follows essentially the approach of Section 3.2(c) and will not be repeated here.
For a Navier–Stokes fluid, an analytical solution to the boundary value problem (84) can be constructed, which, in approach is very similar to that in Section 3.2(b). These solutions were used as a further check of the numerical integration scheme. For brevity they will not be presented here.
3.4. Influence of differences in boundary conditions at the base
All the foregoing calculations are based on the kinematic boundary conditions (29). We have seen that, to zeroth order, this condition is equivalent to (30); but when first-order equations are considered, (29) and (30) differ from each other. It is interesting to investigate in what respect the difference in boundary conditions results in differences in the state of first-order stresses and velocities.
The boundary-value problem that must be solved is now Equations (41), (44), and (48). For vanishing accumulation rate and harmonic solutions (54) and (64) the system (80) evolves as do the boundary conditions (82), but Equations (81) must be replaced by
where
It is then straight forward to show that the differential equations (80) may be written as
with
This system must be subject to the boundary conditions
and
The boundary-value problem (87) to (90) agrees formally with that of Section 3.2(c). The technique for numerical solution therefore also agrees with the one presented there.
4. Results
Several interesting features can be obtained from a numerical exploitation of the theoretical results derived above. These will be discussed now in due order. We begin with the effect of bottom undulations, continue to discuss the influence of accumulation rate and the effects of sliding at the bottom, and end by discussing the effect of temperature variations on the longitudinal strain effects in cold ice. In the following, accumulation is neglected except where explicitly stated.
4.1. Transfer of bottom protuberances to the surface
We begin with the presentation of the filter function ℱ and the phase lag angle φ. Analytical expressions were obtained for a Newtonian fluid and for the case that no simplifications were introduced in the stress representations. It was shown, moreover, that analytical solutions for the first-order stresses and velocities could also be obtained when first-order pressure corrections were discarded. Figures 2 and 3 show the corresponding results. In Figure 2a the filter function ℱ is plotted against a dimensionless wavelength λ = 2π/w for the case when there is no slip along the basal surface. It is clearly seen that the transfer of bottom undulations depends strongly on the mean inclination γ of the ice slope. Generally, ℱ grows with growing λ, and from the graph it appears that it does so in a monotonic way. Consequently, small-wavelength protuberances are filtered out, but large wavelengths are transferred to the surface, even though in an attenuated form. This result is surprising as it is Figure 2 entirely different from Reference BuddBudd’s (1970[a]) prediction. His results look qualitatively like those in Figure 3 in which the filter function ℱ, as obtained from the approximate solution of Section 3.2(a) shows a clear maximum at a preferred wavelength λ. The similarity of this curve with that of Budd stems from the fact that both Budd and we neglect the perturbation pressure in the calculations leading to Figure 3. Quantitative differences, on the other hand, emanate from differences in the boundary conditions at the free surface. Incidentally, the behaviour of ℱ at wavelengths which are below λ = 0.7 is somewhat awkward, because ℱ does not monotonically decrease with decreasing λ but shows an intermediate maximum at λ ≈ 0.45. This indicates that the approximate solution must fail at small wavelengths. Since the least approximations have been made in the calculations leading to Figure 2 we must reject both Budd’s and our approximate solutions. We have also shown in Figure 2b the phase lag angle φ. Except for small wavelengths it is nearly independent of wavelength and depends on the inclination γ. For very small γ this phase lag is close to —π/2.
Fig.2. (a) Fitter function ℱ plotted against dimensionless wavelength λ = 2π/ω for a Navier-Stokes liquid and parameterized for various mean inclinations γ. No slip at the bottom. (b) Phase lag angle φ for the same situation as in (a).
Fig.3. (a) Fitter function ℱ plotted against dimensionless wavelength λ = 2π/ω for a Navier–Stokes liquid and parameterized for various mean inclinations γ. Approximate solution as derived in Section 3.2a. No slip at the bottom. The spike at small values of λ indicates a failure of the approximate integration procedure at small wavelengths.
Results have also been calculated for an ice slope sliding down its bed. As far as ℱ and φ are concerned these are virtually indistinguishable from those obtained with the no-slip condition, hence no further results need be presented.
There is still the possibility of a dominant transfer of bottom undulations to the surface at a certain distinct wavelength when non-linear constitutive behaviour is considered. That this possibility can also be excluded is demonstrated in Figure 4. In Figures 4a and 4b the filter function ℱ is displayed as a function of λ = 2π/ω as obtained by use of the numerical integration scheme described above for a generalized Glen flow law with exponent n = 2 and 3 (the value of
Fig.4. Filter function ℱ plotted against λ = 2π/ω for various mean inclinations γ and for a generalized Glen flow law with n = 2, (a), and n = 3, (b). The value of
Phase lag angles could also be shown for the above cases, yet differences from Figure 2b are insignificant so corresponding graphs are not shown.
Finally, the influence of the
Fig.5. Filter function ℱ plotted for γ = 0.1 against λ = 2π/ω parameterized for various values of
4.2. Basal stresses
The effect of longitudinal strain can be estimated by evaluating the stresses tangential and perpendicular to the boundary at the base. It is easily seen that
where
Θ and Σ depend on y, λ = 2π/ω, p, γ, and ū (which through Equation (38) may be expressed in terms of
These functions are indicators for the significance of the effect of longitudinal strains upon basal shear and for the limitations of the mathematical approach applied herein. Within errors O(ϵ 2), the second of Equations (91) gives the pressure normal to the base. Thus, we must necessarily have (Θ0, Σ0) = O(1). The functions Θ0 and ϕ T are plotted for no slip and for various inclination angles γ, in Figures 6 and 7 against λ. It is clearly seen that perturbation amplitudes for the shear stresses grow with decreasing wavelength. They are generally larger the smaller γ. Of course, this property is directly traceable to the first of Equations (92). The singularity associated with the limit γ → 0 is, however, spurious as explained previously already. At large λ, Θ0 grows with decreasing γ, but at small wavelengths, where the shear stresses become large and a singularity seems to develop, it is the reverse. The singularity is not a numerical peculiarity, because it has also been obtained for n = 1 using the analytical solution represented by Equations (64) to (74). It follows, therefore, that the perturbation solution must break down for wavelengths λ < c. 1. The reader should be aware of this limitation of our perturbation scheme; that it develops is not surprising, however. For if at fixed amplitude undulation wavelengths become smaller and smaller, their higher derivatives will correspondingly increase and eventually no longer satisfy the assumptions of small perturbation theory. Of course, the validity also depends on the value of ϵ.
Fig. 6. First order amplitude of shear stress Θ0 plotted against λ and for various values of the inclination angle γ. The exponent in the generalized Glen flow law is n = 2, and
Fig. 7. (a) Same as Figure 6, but for n = 3 and
For ū = 0 phase lag angles ϕ τ are shown in Figure 7b. They are all nearly independent of λ and close to zero. Similar results have also been obtained for n = (1, 2) for which reason the corresponding graphs will be omitted.
Besides shear, the other significant stress component is the stress normal to the surface. To order ϵ2 this is given by σ y . For no slip and n = (2, 3) its first-order amplitudes Σ0 (see Equations (93)) are plotted in Figure 8 and 9a. Generally, a weak dependence of Σ0 on λ can be observed, and, except for γ ˃c. 0.2 and λ < 3, Σ0 increases with increasing λ. Only for γ ˃ c. 0.2 and λ < 3 is a growth of Σ0 with decreasing λ observed. This seems to be a property of the non-linear material behaviour, because it did not occur for n = 1. Particularly interesting in Figure 8 and 9a is the fact for all wavelengths λ > 2, and for all inclination angles Σ < 1. Only at small wavelengths can the onset of a singularity be observed. Hence, the sign of σ y will, for ϵ < 1 and for most wavelengths, be that of
Fig. 8. First-order amplitude of normal stress Σ0 plotted against λ and for various values of the inclination angle γ. The exponent in the generalized Glen flow law is n = 2 and
Fig. 9. (a) First-order amplitude of normal stress ∑0 plotted against λ and for various values of the inclination angle γ. The exponent in the generalized Glen flow law is n = 3 and
4.3. Surface velocities
It is interesting to see whether bed undulations could also be seen at the surface when velocities are looked at. Zeroth-order velocities can be evaluated from Equation (37) by setting y = 1; they will not be plotted here. First-order velocities can easily be calculated from the stream function ψ by mere differentiation. This allows the representation
where the index s signifies evaluation at the surface and u s is the zeroth-order surface velocity given in Equation (37) (for y = 1). Clearly, first-order velocities will depend on the bottom boundary condition. Here, we present results for no slip and postpone corresponding results for sliding to Section 4.5.
Results for n = 2 and n = 3 do not differ significantly, so that only those for n = 3 will be presented. In Figure 10 the first-order longitudinal velocity amplitude U 0 and the corresponding phase lag angle ϕ u are plotted; generally, the amplitudes U 0 grow with growing wavelength; they are larger, the smaller γ. For small wavelengths they tend to zero. Since U 0 = O(1) for small γ, it is the bottom protuberance amplitude ϵ which indicates whether first-order velocity corrections should be accounted for. The phase lag angle ϕ u undulates, being positive at moderate wavelengths and negative for larger ones. For small γ, ϕu is nearly zero, so that bottom perturbations and longitudinal surface velocities are nearly in phase in this case. For steep glaciers this, however, is no longer so.
Fig.10. (a) First-order amplitude of longitudinal velocity U0 as a function of wavelength λ and plotted for various inclination angles γ. The exponent in Glen’s flow law is n = 3 and
The first-order surface velocity component in the y direction is plotted in Figure 11. Accordingly there exists a distinct wavelength λ (which depends on mean inclination) for Figure 10 which V 0 is maximized; moreover, V 0 is O(1) for large γ and becomes negligibly small for small γ. The phase lag angle ϕv is shown in Figure 11b. For large λ it is nearly independent of wavelength, but shows a clear dependence on γ. For small γ, bottom undulations and surface velocity components are nearly out of phase.
Fig.11. (a) First-order amplitude of normal surface velocity V0 as a function of λ and γ. The constants in Glen’s flow law are n = 3 and
4.4. Effect of steady accumulation rate
We showed in Section 3.8 how steady accumulation rate can be accounted for and also indicated that explicit solutions for a Navier–Stokes fluid were constructed. Calculations were performed for a sinusoidal variation of the accumulation rate. This could be regarded as one term of a general harmonic analysis. We chose
Dependent upon the value of Re A, wavelength (which is of the order 100–102) and the surface velocity, Re H may be as large as O(1). Calculations for H = 1 thus lead to realistic results. We refrain from presenting all pertinent details, but shear-stress and normal-stress amplitudes Θ0 and Σ0 may be presented. They are displayed in Figure 12 and demonstrate that the effect of accumulation rate on basal stresses is of the same order as that due to bottom protuberances (compare Figure 8 and 9a). Neglect of accumulation rate in the evaluation of longitudinal strain effects is therefore not in general permissible. Incidentally, this same inference can be drawn also from the longitudinal velocity amplitudes, which for H = 1 are even larger than those for bottom undulations.
Fig.12. (a) First-order amplitude of shear stress Θ0 for sinusoidal steady-state accumulation rate,
4.5. The effect of sliding at the bottom
All previous numerical results have been presented using the no-slip condition at the bottom. The general theory has, however, been developed for two different sliding conditions, namely that
-
(i) a sliding velocity ū is prescribed and
-
(ii) a Weertman-type sliding law,
, is used.
There are at least two good reasons for studying both these cases. First, zeroth-order solutions are not different if we simply relate ū and
In the calculations described below, ū has been varied, taking values between 10−5 and 10−3. (If τ = 5 d, D = 500 m, ū = 10−3 then the sliding velocity is 1 m/d. ū = 10−3 may thus be regarded as a reasonable non-surging upper bound.) For each inclination angle γ and for m = 2, the corresponding value of
Fig.13. Plot of ū =
Before we present the results, it should be mentioned that, as was the case for no slip, the perturbation scheme applied will fail for small wavelengths. For no slip, the lower bound was roughly λ = 1; with sliding it increases and is as large as λ ≈ 5 when ū = 10−3. The reason for this behaviour is a rapid growth of
4.5.1. Transfer of bottom topography to surface
Filter functions for the transfer of bottom undulations to the surface were plotted for both the sliding conditions mentioned above and for five different values of the sliding velocity ū in the range 10−5 ≤ ū ≤ 5 × 10−3. For case (i) (ū prescribed at the bottom) filter functions do not change in character, when compared with those for no slip. Transfer of bottom undulations to the surface is very small for small wavelengths and increases with increasing λ. In particular, no distinct wavelength exists for which transfer would be maximized. Moreover, with increasing ū the dependence of the filter function on γ becomes less and less. Corroboration for this is provided by Figure 14, in which at most the onset of an intermediate maximum of the filter function can be observed, and this only for specific values of γ. For small γ (<0.01) and for λ < 2 these results must not be taken for granted, since first-order stresses become large in this case and perturbation solutions must fail.
Fig.14. Filter function for the transfer of bottom undulations. Dimensionless sliding velocity at the bottom is ū = 5 × 10−4 and the parameters in the generalized Glen flow law are n = 3 and
The corresponding behaviour for the Weertman boundary condition is quite different. At small ū (= 10−5, 10−4) filter functions show the same characteristic behaviour as they do for no slip. Maximum transfer occurs at very large wavelengths. At ū = (5×10−4, 10−3) and for small mean inclinations γ ≤ 0.05 a clear maximum exists between 3 < λ < 5, which is akin to the result obtained by Budd in an analysis we claimed to be invalid, see Figure 15. Our results seem to give the correct answer to the original problem tackled by Budd, but unlike Budd, who claims that a distinct wavelength for optimal transfer of undulations exists independent of the sliding law, we conclude that in order for this optimum to exist, sliding velocities must be right; if they are too low no clear maximum develops; if they are too high (ū = 5× 10−3, for instance, not graphically displayed here) such a maximum has disappeared again. At ū = 5 × 10−4 it exists for γ = (0.01, 0.05).Footnote * With D = 1000 m and τ = 1 d the corresponding basal sliding velocity would be 0.5 m/d. It thus appears that the cause of the predominant transfer frequency is Weertman-type sliding paired with appropriate mean inclinations. This must be so at least for glaciers which are temperate throughout. Whether this behaviour persists also for cold ice sheets has yet to be seen.
Fig.15. (a) Filter function ℱ plotted against wavelength λ for various inclination angles γ. The dimensionless sliding velocity is ū = 5 × 10−4 and the parameters in Glen’s flow law are n = 3 and
4.5.2. Basal stresses
Differences in sliding conditions are also detectable, when first-order stresses are looked at. For the same range of sliding velocities and for both boundary conditions at the base the first-order basal stresses were evaluated. The results may, perhaps, be summarized as follows: For boundary condition (i) first-order shear stresses react critically to sliding velocities, provided that mean inclinations are small. For larger values of γ, the behaviour is very similar to that for no slip. Figure 16 may serve as evidence for this. It shows the first-order shear stress amplitude and the corresponding phase angle for the case that ū = 5 × 10−4. When compared with Figure 7, it is seen that both amplitude and phase angle depend more significantly on γ than is the case for no slip. The perturbation solution becomes questionable for γ = 0.01. This becomes even more pronounced for larger values of ū.
Fig.16. (a) First-order amplitude of shear stress Θ0 as a function of λ and γ. Note that results for γ = 0.01 are outside the range of the plot. Parameters in the generalized Glen flow law have the values n = 3,
It is quite different for the Weertman-type boundary condition (ii)! Here sliding seems to “stabilize” the perturbation scheme, as first-order corrections decrease in size with increasing sliding velocity. As an example we show the first-order shear stresses for ū = 5 × 10−3 (Fig. 17). For wavelengths λ ˃ 1 shear-stress amplitudes are all O(i) irrespective of the mean inclination angle. Corresponding phase angles are very small. For smaller sliding velocities the behavior is similar with less pronounced independence of the shear stresses on γ (compare with Figure 7) and with an increasingly singular behaviour at small wavelengths, λ < 2 say.
Fig.17. (a) First-order shear stress amplitude Θ0 as a function of λ and γ. Notice the regular behaviour for all λ > 1 and all γ. Parameters for Glen’s flow law are n = 3,
As far as first-order normal stresses σ y are concerned, results are reasonably insensitive to differences in boundary conditions ((i) versus (ii)), but the dependence on sliding conditions is pronounced. As was shown in Figure 9 for no slip, first-order normal stresses depend critically on γ; when γ is large and λ is small, Σ0 becomes large. With growing sliding velocity ū, the dependence on γ becomes smaller and smaller, but the singular behaviour at moderately small values of λ is also carried over to small inclinations γ. As an example, we show in Figure 18 the first-order normal stress amplitude Σ0 for the Weertman-type boundary condition and a sliding velocity ū = 5 × 10−3. As can clearly be seen, the perturbation scheme must fail at moderately small wavelengths at all indicated mean inclination angles. At smaller sliding velocities the behaviour lies between those shown in Figures 9a and 18.
Fig.18. First-order amplitude of normal stress Σ0 as a function of λ and γ. Parameters in the generalized Glen flow law are n = 3 and
We conclude by stating that first-order basal shear stresses must depend critically on the sliding law, in particular at small angles γ. Of the two sliding laws, that due to Weertman (ii) leads generally to smaller first-order stresses at comparable sliding velocities. It implies that the perturbation analysis shown here remains valid for a wider range of wavelengths (down to λ = 1) and inclination angles γ than is the case for the other boundary condition (i). Of course, condition (ii) is more reasonable than (i), but this is no proof of the suitability of the former. Our analysis rather indicates that first-order shear stresses in a temperate ice sheet may depend critically on the basal boundary condition. As long as experts debate the superiority of one over the other, it seems therefore illusory to estimate longitudinal strain effects in temperate glaciers.
4.5.3. First-order velocities
It is now no longer surprising that differences in the boundary conditions manifest themselves also in the first-order longitudinal velocities. For the boundary conditions (i) the velocity amplitude U 0 depends strongly on the value of ū; for the Weertman-type condition (ii) it does not. As an example, the graphs of Figure 19 may serve. They show the first-order longitudinal velocity amplitude U 0 for a value of ū = 5 × 10−3. In Figure 19a the results for boundary condition (i) are displayed. Only plots for γ = (0.4, 0.2) are shown as those for smaller γ lie outside the range of the ordinate. This behaviour is quite unlike that for boundary condition (ii)! Here, U 0 remains O(1) or smaller for all inclination angles 0.01 ≤ γ ≤ 0.4. There is a smooth transition in this case from the behaviour for no slip (see Fig. 10) to that of Figure 19b, where sliding velocities increase with no indication of unusual behaviour. The distribution of U 0 for boundary condition (i) is completely different. At certain sliding velocities it behaves qualitatively as shown in Figure 19a, at others more like that in Figure 19b.
Fig.19. (a) First-order longitudinal velocity amplitude U0 as a function of λ and γ. The parameters in the generalized Glen flow law are n = 3 and
Interestingly, the vertical surface velocity is insensitive to the form of the boundary condition at the bottom. But sliding does affect the size and the wavelength behaviour of V 0 and also its phase. For corroboration the reader may compare Figures 11a and 20, the latter being valid for a Weertman-type boundary condition and a sliding velocity u = 5 × 10−3.
Fig.20. Amplitude of normal velocity V0 as a function of λ and γ. Parameters in the generalized Glen flow law are n = 3 and
Many more results could be shown, but space limitations prevent us from doing so.
4.6. Temperature variations and the longitudinal strain effects in cold ice
All the foregoing calculations were performed for ice sheets with spatially independent material properties. For temperate glaciers such conditions are satisfied to a sufficient degree of accuracy, but in cold ice sheets such as Greenland and Antarctica, material properties change drastically with depth, because of the appreciable dependence of the Arrhenius factor on temperature. It is interesting to investigate to what extent the above analysis remains correct when temperature variations perpendicular to the main flow direction are taken into account.
First, results for an ice sheet which is cold throughout and therefore adheres to the bed have been obtained by Hutter and Spring (1979). They demonstrate that transfer of bottom undulations, first-order basal shear and normal stresses, and first-order surface velocities, are all attenuated by vertical temperature variations.
If the ice sheet reaches the pressure-melting point at the base, sliding is possible; the inferences drawn by Hutter and Spring may then change. Our analysis remains valid when material properties change with depth. All that is needed is an adjustment of the form of the function
where Q is the activation energy, k the Boltzmann constant, T 0 the surface temperature, and T(y) the temperature distribution with depth. Calculations were performed for
in which θ 1 = 258.8 K, θ 2 = 1.05, L* = 3.70, and erf is the error-function. This choice corresponds to an ice sheet 1400 m thick whose bottom boundary is at the melting point (272.12 K) and of which the upper surface is at the temperature 258.8 K; needless to say, Equation (97) is a solution of the heat conduction equation.
In the zeroth-order solution, temperature variations manifest themselves in the longitudinal velocity distribution; the first of Equations (37) remains valid as an expression for
Calculations for the first-order problem were performed for sliding at the bottom according to both sliding laws (i) and (ii), and ū was given the values 10−4 and 10−3. The results may, perhaps be summarized as follows:
-
(1) As regards the transfer of bottom undulations, filter functions are reduced in general by the vertical temperature variation. This is uniformly so for large wavelengths, but not necessarily so for smaller ones. Moreover, the effect of maximum transfer of bottom protuberances to the surface, which has been observed for some angles (see Figure 14), may be amplified by temperature variations. Corroboration is given by Figure 21, which shows the filter function for ū = 10−3 and Weertman-type sliding (maxima for sliding law (i) are less pronounced).
-
(2) First-order shear-stress amplitudes are also reduced by the inclusion of temperature variations. This is uniformly so for all investigated sliding velocities and for either sliding law. As an example, we show in Figure 22 the shear-stress amplitudes Θ0 as obtained when the temperature distribution given by Equation (97) and a Weertman-type boundary condition were used. The dimensionless sliding velocity was ū = 10−4. For comparison, results are also shown for a temperate glacier and an ice sheet whose material properties are kept constant and which adheres to the bed. Figure 22 is sufficient proof of the importance of both sliding and temperature variation.
-
(3) Equally interesting are results regarding first-order normal stresses Σ0. Whereas sliding (as opposed to no slip) has not led to a reduction in the first-order normal stresses (see Figs 9a and 18), the inclusion of temperature variations brings a reduction, especially at wavelengths 1 < λ < 5, at large mean inclinations γ, and at small sliding velocities. For no slip it seems to be largest. Figure 23 may serve as an example, it shows results for ū = 10−4 and the temperature distribution of Equation (97), and a comparison for constant temperature.
-
(4) First-order velocities U 0 and V 0 are also changed considerably by the inclusion of vertical temperature variations. For γ ˃ 0.1 these changes seem to be insignificant, but for smaller values of the mean inclination γ the changes are substantial. Figure 24 may serve as an example.
Fig.21. Filter function ℱ for a cold ice sheet plotted against λ for various values of γ. The temperature distribution was according to Equation (97) and the parameters in Glen’s flow law are n = 3,
Fig.22. First-order amplitude of shear stress as a function of λ and γ. The temperature distribution was according to Equation (97) and the parameters in Glen’s flow law were n = 3 and
Fig.23. First-order amplitude of normal stress as a function of λ and γ. Shown are results for ū = 10−4 and a temperature distribution according to Equation (97). For comparison results are also shown for a constant temperature distribution. The parameters in Glen’s flow law are n = 3 and
Fig.24. First-order amplitude of longitudinal velocity as a function of λ and γ. Results are shown for ū = 10−3 and a temperature distribution according to Equation (97). For comparison results are also shown for a constant temperature distribution. The parameters in Glen’s flow law are n = 3,
5. Concluding remarks
In this article, first-order stresses and deformations in an ice slope have been determined. The fundamental governing equations have been formulated on the basis of the continuum mechanics of slow viscous flow of non-linear fluids. In the solutions, efficient use has been made of perturbation approximations, first to separate steady-state from time-dependent problems, and, secondly, by assuming that bottom undulations deviate only slightly from a plane. Using the amplitude of the bottom undulations essentially as the perturbation parameter, first-order corrections of the velocities and stresses could then be determined from a linear boundary-value problem. A detailed analysis of the steady-state flow of a nearly parallel-sided ice slope has been given, and emphasis laid on the transfer of bottom undulations to the surface, on longitudinal strain effects on basal stresses and surface velocities, and on the effect of spatial variations of these. Results can be summarized as follows.
Contrary to previous claims, there is no universally preferred wavelength at which transfer of bottom protuberances to the surface is optimal. Generally, this transfer increases with increasing wavelength, yet at very distinct mean inclination angles and distinct values of the sliding velocity, filter functions do show a preferred wavelength at which transfer of bottom undulations to the surface is maximal. These wavelengths lie between 3 < λ < 5, the maxima appear to be more pronounced when vertical temperature variations are taken into account than when constant material properties are assumed throughout. (We do not know whether mean inclinations and sliding velocities in the instances when preferred wavelengths were observed in the field were such that a preferred wavelength transition to the surface does exist.) The differences between our theoretical findings and those of Budd are important as they indicate that the wavelength spectra of steady-state surface undulations may not possess a distinct wavelength the cause of which could uniquely be related to bottom protuberances. If sliding and inclination are right, they do, otherwise one might speculate that x-dependent zeroth-order solutions, of the kind which must occur below an ice fall, might provide the answer to it (see Reference WilliamsWilliams, 1979).
As far as basal stresses are concerned, the analysis has indicated that their order of magnitude strongly depends on the spatial temperature distribution and on the type of slidinglaw. Longitudinal strain effects on basal stresses are roughly twice as large in temperate ice as they are in cold ice. Moreover, for a Weertman-type sliding law they are smaller than for a sliding law in which sliding velocities are prescribed. This result is particularly important, as it indicates that the first-order stress distribution depends critically on the sliding law. Therefore as long as experts debate the proper form of the sliding law it is illusory to estimate longitudinal strain effects on basal stress.
Differences due to basal boundary conditions are even more pronounced in the surface velocities than they are in the basal stresses. For Weertman-type sliding, first-order surface velocities are generally smaller than when sliding velocities are prescribed.
Analysis of the accumulation-rate effects on the basal stress distribution showed further that it is in general not justified to neglect their effect.
The above results are all subject to the validity of the applied perturbation scheme. Results indicate that this scheme breaks down at small wavelengths; the lower bound of wavelength to which solutions remain valid depends on the inclination of the ice slab and the sliding law. Typically, wavelengths must not be smaller than the ice thickness and for small inclinations of the ice slab may have to be as large as three times the ice thickness.
As a consequence of these facts, calculations of the effects of longitudinal strain on basal shear as performed by Reference BuddBudd (1971) are questionable, and conclusions from them must most likely be invalid. We could, if we so desired, repeat the Budd’s analysis with our improved approach, but shall not do it here for the following reason. Bottom and top surfaces were assumed to deviate only slightly from straight parallel lines and protuberance amplitudes were assumed to be small compared with the glacier thickness. Such assumptions are somewhat unrealistic; top and bottom surface are rarely parallel, but their change along the axis is usually slowly varying. Stress variations with ice flow over undulations is thus better analysed using this less stringent hypothesis. Hutter (1981) presents the pertinent details.
The analysis also showed that the temperature variation, which manifests itself in a position-dependent stress-strain relationship, may be taken into account without further difficulties. Contrary to the usual understanding, however, the ice slope need not be divided into a series of strips of uniform properties and, in particular, no restriction to a linear viscous law is required. Material properties may vary continuously according to independently-taken temperature records. Owing to space limitations, results have only been sketched in this regard. A detailed analysis for no slip has been given by Hutter and Spring (1979), and results for sliding at the bed have only been demonstrated for a very limited range of applications. They showed that under no circumstances is it allowable to disregard temperature variations in cold ice, and in particular, it is shown that filter functions and basal stresses are reduced by temperature variations.
The above discussion is restricted to steady-state problems, but the theory was presented initially also for time-dependent problems. In particular, it would be interesting to see whether a kinematic wave theory could be derived from our governing equations. Formally this is the case; in fact Hutter (1980) gives the corresponding analysis. Explicitly he shows that the governing equations of surface waves travelling down an inclined slope can be derived from the above set of equations by merely invoking a long wavelength. To a first approximation, kinematic wave theory does evolve from the corresponding analysis, but calculations indicate it to be oversimplified. For details the reader is referred to Reference HutterHutter (1980).
Still further extensions are possible, however. The influence of global curvature effects is one of these, the behaviour at small wavelengths would be another. Future investigations will have to deal with these as well as different sliding laws and cavity formation.
Acknowledgements
The analytical part of this paper was performed jointly by the first and second authors. The finite-difference program and the execution of the numerical computations were performed by the third author. The layout and final draft of the manuscript are due to the first author, who, of course, is also responsible for all errors which may remain.
The work of U. Spring was financially supported by the Swiss Nationalfond für Wissenschaftliche Forschung.
