A Glacier Flow Model Incorporating Longitudinal Deviatoric Stresses

Abstract

A model for glacier flow is developed which incorporates longitudinal deviatoric stress contributions to the field equations. The underlying assumptions may be applied to develop models for various situations. Here, they are developed for steady-state and non-steady-state sliding glaciers in plane flow. The models reduce to a proper generalization of plane-flow pseudo-hydrostatic theory if longitudinal deviatoric stresses are neglected in comparison to basal shear stresses. Solution of this simpler reduced model allows an estimate to be made of the magnitude of the longitudinal deviatoric stress to test if it is negligible or, more generally, investigate under what conditions it can be neglected. The steady-state model predicts that longitudinal deviatoric stresses are negligible for very arbitrary non-uniform sliding-law distributions provided that the following conditions exist: the region must be distant from an ice divide or terminus and subject to normal (not extreme) accumulation or ablation. On the other hand, examples are produced where, under non-steady conditions, longitudinal deviatoric stresses are important and even dominant.

Résumé

Le modèle d’écoulement de glacier développé introduit les effets de la composante longitudinale du déviateur des contraintes dans les équations. Les hypothèses sous-jacentes peuvent être appliquées au développement de modèles pour des cas de figures variées. Ici, elles sont envisagées pour des glaciers glissant en écoulement plan, à l’état stationnaire et non stationnaire. Le modèle se ramène à une généralisation adéquate de la théorie pseudo-hydrostatique en écoulement plan si la contrainte longitudinale déviatrice est négligée devant la contrainte de cisaillement basale. Une solution de ce modèle plus simple permet une estimation de l’amplitude de la contrainte déviatrice longitudinale en tant que test de son importance, ou de façon plus générale, en tant que recherche des conditions qui la rendent négligeable. Le modèle d’état stationnaire conduit à des contraintes longitudinales négligeables pour de très arbitraires lois non uniformes de glissement sous réserve de l’existence des conditions suivantes: les régions doivent être éloignées d’une diffluence ou d’un front et soumises à des accumulations ou ablations normales (extrema exclus). D’autres part des exemples produits montrent que dans des conditions non stationnaires des contraintes longitudinales déviatrices ont un effet important et même dominant.

Zusammenfassung

Es wird ein Modell für den Gletsherfluss entwickelt, das längsgerichtete, ablenkende Druckkomponenten in die Feldgleichungen einführt. Die getroffenen Annahmen können zur Aufstellung von Modellen für verschiedene Situationen verwendet werden; hier werden Modelle für den stationären und nicht-stationären ebenen Fluss gleitender Gletscher dargestellt. Die Modelle vereinfachen sich zu einer geeigneten Verallgemeinerung der pseudohydrostatischen Theorie für ebenen Fluss, wenn ablenkende Normaldrucke gegenüber der Scherspannung am Untergrund vernachlässigt werden. Die Lösung dieses einfacheren Modells erlaubt eine Beurteilung der Frage, ob das Ausmass des längsgerichteten, ablenkenden Druckes vernachlässigbar ist, oder allgemeiner, eine Untersuchung, unter welchen Bedingungen dieser vernachlässigbar werden kann. Das stationäre Modell lässt erkennen, dass bei konstantem Fluss ablenkende Normaldrucke gegenüber der Scherspannung am Untergrund für weitgehend willkürliche Verteilungen im nicht-gleichformigen Gleitgesetz vernachlässigt werden dürfen, wenn folgende Bedingungen erfüllt sind: Das Gebiet muss entfernt von einer Eisscheide oder einem Gletseherende liegen und darf nur normale (nicht aussergewöhnliche) Akkumulation oder Ablation aufweisen. Andrerseits werden Beispiele angeführt, wo unter nichtstationären Bedingungen die ablenkenden Längsdrucke wichtig und sogar bestimmend sind.

Introduction

A number of attempts have been made to incorporate the effects of longitudinal deviatoric stresses in models of glacier flow. The deviatoric stress and Cauchy stress are related by , and so have common shear components. In plane flow, in Oxy (Fig. 1) with Ox parallel to the bed inclined at an angle α to the horizontal, τ_xx is termed the longitudinal deviatoric stress, and τ_b is the shear traction τ_xy = σ_xy on the bed. Weertman (1961) was concerned with the situation near an ice divide where the longitudinal deviatoric stress |τ_xx| may be larger than basal shear stress τ_b. Robin (1967), Collins (1968), and Budd (1971) were concerned with the effect of basal undulations upon surface profile; both Robin and Collins assumed that |τ_xx| dominated τ_b. Robin showed by calculation that to the south of Camp Century, North Greenland, |τ_xx| lies between 0.9 and 1.5 bars with τ_b less than 0.8 bars. Nye (1969) considered the effect of τ_xx in the equilibrium equations. These earlier papers included longitudinal deviatoric stresses in an ad hoc manner so that even for the case of plane flow there does not exist a general model which incorporates longitudinal deviatoric stresses.

Fig. 1. A section of a glacier or ice sheet sliding over a flat inclined bed. Positive directions are shown for accumulation b(x), velocity components (u, v), and stress components.

There are other situations not treated in the above papers for which a more general model is required. The common assumptions are that |τ_xx/τ_b|≫ 1 for floating ice shelves whereas |τ_xx/τ_b|≪ 1 for grounded ice sheets. It is clear that a transition region exists near the grounding line where both terms are important. As another example, Shoemaker (1981) found it necessary to invoke the assumption |τ_xx/τ_b.|>>1 in studying transient creep slump in glacier reservoirs. In general, non-steady-state conditions might be expected to produce values of |τ_xx/τ_b| at least of order unity, for example the values cited by Robin above; however there has been no analytical evaluation of such situations. Steady-state problems might produce values of |τ_XX/τ_b| of order unity under conditions of high accumulation or ablation or when basal slope or sliding conditions change rapidly.

The conventional pseudo-hydrosttic theory arising when τ_xx is neglected in comparison with τ_b, so termed by Robin (1967), is shown by Morland and Johnson (1980, 1982) to be the lead-order approximation on the global scale for steady plane flow over a slowly undulating bed. The surface slope relative to the mean bed line provides a small parameter for a perturbation series solution.

We use the abbreviation p-h for this approximate theory. Since the full plane-flow equations do not allow significant analytic progress, there is need for a less restrictive approximation which incorporates effects of the deviatoric stress τ_xx. This paper presents a simplified model which incorporates the deviatoric stress τ_xx in conjunction with Glen’s flow law for the viscous shear response of ice and the incompressibility approximation. The model satisfies the field equations in an average sense, and reduces to p-h theory as a limiting case. It is developed here only for isothermal plane flow of a sliding glacier over a plane inclined bed, but could be extended to a slowly undulating bed. A shape-factor approximation for side drag (Nye, 1965) could be included, and effects of a prescribed temperature profile can be incorporated.

It should be remarked that application of the model does require extensive and delicate numerical calculation even though only a single independent variable is involved.

Steady-State Formulation

In rectangular Cartesian axes Ox_i (i = 1,2,3), with the equivalences x₁, = x, x₂= y, in the flow plane Oxy (Fig. 1), Glen’s flow law for the shear response and the incompressibility approximation are given by

(1)

where ė_ij are the strain-rate components and τ is the second invariant of the deviatoric stress given by

(2)

Henceforth, we adopt n = 3, although the model can be developed for any arbitrary positive integral value of n. The factor A is significantly temperature dependent, particularly near melting. Let (u, v) denote the velocity components in Oxy. We assume that |∂u/∂y| ≫ |∂v/∂x| and henceforth neglect ∂v/∂x. In the p-h theory longitudinal gradients are much smaller than gradients through the thickness and |v| ≪ |u| in the mean (Morland and Johnson, 1980), and non-negligible τ_xx should not annul both these strong inequalities. This assumption does not apply to floating shelves. Hence for plane flow, where ė₃₃ =ė₂₃ ė₁₃ = 0_, τ₃₃ = τ₂₃ = τ₁₃ = 0, τ₁₁ = – τ₂₂, and Equations (1) reduce to

(3)

The average down-slope velocity U(x) is defined by

(4)

where u_b(x) is the basal sliding velocity and h(x) is the glacier thickness.

We retain the assumptions that the gradient of each stress component in the longitudinal direction is much smaller than its gradient through the thickness, and that the surface slope h′(x) is small. If τ_xx = 0(τ_xy), then by Equations (3) the longitudinal velocity u must have comparable gradients in both directions, and then by incompressibi1ity ∂v/∂y = 0(∂u/∂y). These comparisons do not follow from the global coordinate and variable scalings of Morland and Johnson (1980, 1982), but are expected to arise only on a local scale where more rapid longitudinal changes are initiated. With a surface slope magnitude as a small parameter, the free surface conditions correct to lead and first order are given by (Morland and Johnson, 1982)

(5)

Inertia terms are negligible so momentum balance is given by the equilibrium equations

(6)

In the p-h theory, |τ_xx| ≪ |τ_xy| integrating the equilibrium equations (6) subject to surface conditions (5), with an x-derivative yielding a lower order of magnitude than a y-derivative, gives the lead-order stress distribution (Morland and Johnson, 1982)

(7)

which show that τ_XY = 0(h′σ_yy) if α = 0(h′), but τ_xy = 0(σ_yy) if α = 0(1).

Now allow τ_XX = 0(τ_xy). Since ∂τ_xx/∂x is still a lower order of magnitude than ∂τ_Xy/∂y, the lead-order expression given by the second of Equations (7) for τ_xy is still valid, and τ_xx only influences the first-order correction. Integrating the second of Equations (6) subject to the surface given by the first of Equations (5) to include first-order terms requires only a lead order τ_xy in ∂τ_xy/∂x, significant when α = 0(1). Then

(8)

which is again a p-h result, not influenced by τ_XX, and the sin α term is only a first-order contribution when α = 0(1). We cannot obtain an expression for τ_xx independent of the flow solution, but a lead-order relation for the mean longitudinal deviatoric stress is given by integrating the longitudinal balance given by the first of Equations (6) over the thickness, using the second of Equations (5) and the first of Equations (7). For τ_xx = 0(τ_b), with τ_b ≈ _ρgh sin α if a = 0(1) and τ_b = 0(pghh′) if α = 0(h′), the sin α term in Equation (8) does not contribute to the lead-order terms in either case. Define

(9)

then, to 0(h′τ_b),

(10)

The basal shear traction τ_b must be related to the sliding velocity u_b by a prescribed sliding law. Following Nye (1959), we adopt, for illustration,

(11)

with variable coefficient λ(x) to allow for a variation of basal conditions. An explicit exponent 1/3 is introduced for subsequent calculations. Finally, integrating the incompressibility relation ė_kk = 0 through the thickness gives the mass-balance result for steady conditions,

(12)

where b is the net accumulation (surface accumulation less any basal drainage), that is, the volume of ice added per unit time per unit horizontal cross-section.

Reduction of the system of differential equations in x and y to a problem in one independent variable x requires integration of the constitutive Equations (3) through the thickness. In p-h theory where τ_xx is absent to zero and first order, the zero-order integration using the linear τ_xy expression in Equations (7) is simple provided that A is assumed constant; that is, for temperature-independent ice or isothermal conditions. With typical temperature profiles in a cold glacier there will be significant variation of A with depth, particularly in the regions near melting, and use of the average value for the integrand factor is not then a convincing approximation. Given an empirical temperature profile, a realistic variation of A with y can be prescribed. For the present illustration we assume A is constant, which is appropriate to a near-temperate glacier. It is not possible, however, to prescribe the actual variation of τ_xx with y, since this requires solution of the equations in x and y, so to evaluate the mean velocity integrals in Equation (4) we replace τ_xx in the velocity gradients in Equations (3) by its mean value t_xx. This is exact only if τ_xx is independent of y. There is no evidence, however, to suggest an alternative realistic variation of τ_xx which could be used to evaluate these integrals. Since the model delivers only mean values, the approximation cannot be tested for consistency, so this assumption of the y variation of τ_xx must be noted as the basic weakness of the approximation.

Now integration of the velocity gradient in the first of Equations (3) gives u_s – u_b, where u_s = u(x, h) is the longitudinal velocity at the surface, and the repeated integral on the right-hand side of Equation (4) gives U(x). Differentiation of the integral in the middle of Equation (4) yields an integral of the velocity gradient in the second of Equations (3). The results are

(13)

when the lead-order expansion in the second of Equations (7) is written τ_Xy = τ_b(l – y/h). Note that the second and third together with the first of Equations (13) are two independent relations between the five variables U, u_b, τ_b, τ_xx, h, since they result from the two independent constitutive relations in the first and second of Equations (3). These, together with the sliding relation (11), mass-balance Equation (12), and the longitudinal equilibrium equation (10), complete the system of five equations, which will now be considered in various situations.

It is clear from the construction that having τ_XX of order τ_xy does not change the lead-order p-h Stress distributions, but influences the velocity field given by Equations (3). If τ_xx is neglected in comparison with τ_xy (one order of magnitude smaller in p-h theory, then the first of Equations (3) is the lead order p-h relation, while the second of Equations (3), and hence the third of Equations (13), are superfluous to the lead-order solution. The required system is Equations (10), (11), (12), and the second of Equations (13), for four variables U, u_b, τb, h, though the p-h equations can be solved directly for both x and y variation, not just in an average sense (Morland and Johnson, 1980, 1982). Thus the present model is consistent with p-h theory when τ_XX is negligile, but allows the magnitude of τ_xx, as measured by the mean stress t_xx, to be estimated.

Dimensionless Variables

We are concerned primarily with effects induced by changes over longitudinal scales of order the thickness, so we introduce dimensionless coordinates based on a thickness h(0). The origin, in the region of interest, is supposed not to be at a divide, so that u_b(0) ≠ 0, τ_b(0) ≠ 0, and U(0) is a non-zero measure of longitudinal velocity. Define

(14)

Thus ξ, H, ū_b, T_b, and T_xx, when τ_XX = 0(τ_xy),_ are order-unity variables, but away from a divide b is small. λ(ξ) describes a prescribed variation of sliding conditions over a thickness length scale. Typical smooth variations of accumulation b over this length scale will not induce longitudinal deviatoric stresses, so we adopt uniform b for our examples (which would not be realistic over the scale of the span); thus b is a constant dimensionless parameter measuring the accumulation magnitude relative to the mean longitudinal velocity. The other natural dimensionless parameters are α and

(15)

which measure the ratios of the overburden pressure and mean longitudinal deviatoric stress to the basal shear stress, and basal sliding velocity to mean longitudinal velocity, respectively, at the chosen origin. Away from a margin, q will normally be a large parameter. Note that 0 < ū₀ < 1 since the longitudinal velocity increases towards the surface.

We can now eliminate U and u_b from the steady-state equations in dimensionless variables to obtain three equations for T_b, T_xx, and H. The constants Ah(0) τ³ _D(0) and λ(0) are eliminated in favour of the parameters ū₀ and T₀ by evaluating the second of Equations (13), and Equation (11) at ξ = 0. Integrating the mass-balance Equation (12) and eliminating U between the second of Equations (13) and Equation (11) gives an algebraic relation

(16)

Eliminating U’ between the third of Equations (13) and Equation (12) with u_s given by the first of Equations (13), and then U by the integral of Equation (12), gives the slope expression

(17)

Finally, the longitudinal balance Equation (10) becomes

(18)

where h′ can be eliminated by Equation (17). A steady-state solution will not exist for all choices of parameters defining conditions at ξ = 0. For example T₀<0 and b > 0 implies positive accumulation and compressive flow near ξ = 0 which would require H to increase with time. Such situations are reflected by rapid growth or instability of the numerical solution of Equations (16) – Equations (18). No numerical instabilities have arisen with sensible choices of the parameters.

The restriction to steady state lies solely in the mass balance equation (12), which incorporates the surface accumulation condition for h stationary in time. The constitutive Equations (3), equilibrium Equations (6), and sliding relation (11), apply equally to non-steady flow. Thus the system of four Equations (1), Equations (11), and second and third of (13) with the first of Equations (13) may be applied to non-steady conditions, but cannot determine the five variables U, u_b, τ_xx, and h, nor any actual time variation of the solution. However, if a profile h(x) is prescribed for some fixed time, these four equations determine the instantaneous distributions of U, u_b, τ_b, and t_xx, and in particular the magnitude of t_xx relative to τ_b. Such estimates naturally depend on the choice of h(x) and the parameters defining conditions at x= 0. An application which starts from a steady-state profile determined with a given basal sliding condition, and investigates the new velocity and stress distributions induced by a change of the sliding function λ(ξ), is presented later. Since the time scale of change between steady states will be many hundreds of years or more, there will be a finite period following the change of λ(ξ) during which H(ξ) remains approximately at the initial steady profile.

Eliminating and λ(0) by the second of Equations (13) and Equation (11) as before, and u_b by Equation (11), and comparing U given by the second of Equations (13) and the integral of the third of Equations (13), gives the integral equation

(19)

Integrating Equation (18) gives a second integral equation

(20)

The two integral Equations (19) and Equations (20) for T_b and T_XX are preferred for numerical calculation. First-order differential relations can be obtained by comparing the third of Equations (13) with the derivative of the second of Equations (13) and rewriting Equation (10) as an expression for .

If we set T_xx ≡ 0, then Equation (18) reduces to the p-h result for T_b, given by y= 0 in the second of Equations (7), which involves both H and h′. Then Equation (16) is a first-order ordinary differential equation for H(ξ) in a region of constant accumulation, subject to H(0) = 1, which is a special case of the equation for arbitrary b(x) (Morland and Johnson, 1982). That is, Equation (16) determines the profile predicted by p-h theory. Given a valid small surface-slope profile, the full perturbation solution (Morland and Johnson, 1980, 1982) now automatically predicts a longitudinal stress τ_xx one order of magnitude smaller than τ_xy for a bounded viscosity law; there are singularities at the free surface if a power law with exponent greater than unity is used. With the present model, which retains the integrated stress T_xx in comparison with T_b, we have Equation (17) which determines T_xx(ξ) in terms of H(ξ). Thus, given the p-h profile H(ξ), we can use Equation (17) to test the underlying approximation |T_xx|≪T_b, and confirm, or reject, in this mean sense, the applicability of the p-h theory. Solution of the complete set of Equations (16) – (18) for steady state is necessary when T_xx is not everywhere negligible in comparison with T_b, and will show in particular how the p-h profile is modified. The same approximations also provide an estimate of T_xx in the non-steady-state case. However, in the non-steady-state applications of the model illustrated later, T_xx does not remain negligible and p-h theory is inappropriate.

Steady-State Illustrations

To obtain a 10% change from the p-h results, the constitutive equations (3) and slope expression (17) suggest that (T_xx/T_b)² should exceed 0.1 and 0.3 respectively in the region of interest, say T_xx/T_b ⩾ 0.2. We will first take T₀ = 0 as the initial condition for T_xx, at ξ = 0, to avoid imposing a significant longitudinal deviatoric stress artificially, and investigate the possible build-up of T_xx over a region 0<ξ<10 for different values of ūb and q. Bed inclination α will not have any significant effect on T_xx, so for simplicity we take α = 0. Two values of the ratio of the initial overburden pressure to basal shear stress are considered; q= 400 which corresponds to ice-sheet thickness 1.6 km and basal shear stress 0.36 x 10⁵N m⁻², and q = 25 which corresponds to thickness 100 m at the same basal shear stress. A range 1/10 to 1/1.001 for the ratio ū₀ of initial basal velocity to mean longitudinal velocity is considered.

First we eliminate possible influence of accumulation (small b) and variation of sliding coefficient λ(ξ) by taking b = 0, λ = 1. Table I gives a summary of the solution of the initial-value problem defined by Equations (16) – Equations (18). In all cases the arbitrarily selected initial condition T₀ = 0 causes T_xx to increase rapidly in a narrow boundary layer of dimensionless width H^⋆. (The values shown in Table I were chosen subjectively, but the layer is rather prominently delineated if q = 400, rather less so if q = 25.) Outside the boundary layer T_xx/T_b enters a plateau region where it increases slowly at an almost constant but increasing slope. Average values of T_xx/T_b outside the boundary layer are shown in Table I as well as average values of the slope (T_xx/T_b)′. The results show that H^⋆ increases rapidly as ū₀ → 1 or as q → 0. This is to be expected since internal adjustment to an arbitrary initial condition, in this case T₀ = 0, can be accomplished only through internal deformation which vanishes as ū₀ → 1 or q → 0. The same reasoning explains the increase in (T_xx/T_b)′ with increasing u_b(0)/U(0) or decreasing q. Table I shows that plateau values of T_xx/T_b in excess of 0.2 are possible, but only under conditions of very low sliding resistance (ū₀ → 1), very thin ice, or some combination of these two limiting conditions. Certainly, near the terminus both of these limiting situations exist. On the other hand, near an ice divide where q would be large, it is not clear what magnitude ūo = u_b(0)/U(0) would approach.

Table I. Dimensionless Width H^⋆ of Boundary Layer. Average Plateau Value of T_xx/T_b Outside Boundary Layer Over Interval (0.10) and Slope (T_xx/T_b)′ in Plateau Region. All For α = b = T₀ = 0, λ = 1

Near an ice divide (U(0) → 0), b becomes large, and the present dimensionless formulation is not satisfactory. It can be seen that in the slope expression (17), the terms involving b become dominant if |T_XX| is order unity or less. Calculations were made corresponding to finite non-zero b and it was verified that values of T_xx/T_b of order unity can easily arise under steady-state conditions near ice divides. This was the situation which Weertman (1961) studied.

As earlier remarked, these conclusions concerning the order of magnitude of deviatoric stresses near glacier termini and ice divides were expected. The conclusion that |T_xX/T_b| ≪ 1 in regions removed from ice divides and termini, provided there is an “appreciable” sliding resistance and |h′| ≪ 1, is in agreement with the global perturbation analysis (Morland and Johnson, 1980, 1982).

The next question is whether a rapid change in λ(ξ) can produce values of T_xx/T_b of order unity. The situation corresponding to a central reservoir with a constant reduction of sliding resistance λ was examined for a wide range of other parameter values. The results are in line with Table I. That is, the effect of an elevated step in λ is to produce a transition boundary layer with a wide of order H^⋆ where variables change rapidly and continuously in adjusting to plateau values of T_xx/T_b similar to those given in Table I. Note that the boundary layers also show up clearly in Figures 2 and 3 and in the corresponding Figure 4 for non-steady-state cases to be discussed. Thus, the model predicts that under steady-state conditions and for fixed q and ū_o, the order of magnitude of the maximum values of |T_XX/Tb| in a region of length of order the thickness is determined primarily by the value of λ in the region rather than by the distribution of λ (ξ) outside the region, except when the boundary layer spreads over that distance. Moreover, discontinuities in sliding resistance do not by themselves produce large longitudinal deviatoric stresses.

Fig. 2. Illustrations of the effect of q and T₀ upon the dimensionless variables H, T_xx, and T_b in steady-state flow. The T_xx curves start from (0,−1) and the T_b curves start from (0,1). Two cases shown correspond to α = 0, b = 0, λ = 1, T₀ = −1, ū_o = 0.5, with (i) q = 400 and (ii) q = 25.

Fig. 3. The effect of a sudden reduction in sliding coefficient in a central region upon dimensionless variables T_b and flux F. Height H and the pseudo-hydrostatic basal shear stress T_ph result from a steady-state solution with λ ≡ 1, q = 400 ū₀ = 0.2, α = b = 0. T_b and F result from steady-state H input to non-steady-state equations with λ = 1/5 on 0.25 < ξ < 0.5. Arrows indicate a near-jump discontinuity.

Fig. 4. The effect of a sudden reduction in sliding coefficients from unity in a central region 0.25 < ξ < 0.5 upon T_xx/T_b. Cases shown are q = 100, λ = 1/8 (---), q = 400, λ = 1/5 (—), q = 100, λ = 1/2 (⃜), all for u₀ = 0.2.

The effect of starting from an order unity T_xx is illustrated by Figure 2 for the situation of a uniform sliding coefficient λ = 1 on (0,8) with T₀ = −1, b = α = 0. The glacier profile corresponding to q = 400 is very flat and T_xx is close to zero except in a very narrow boundary layer where it increases rapidly from T₀ = −1 to a small positive value. Similarly, T_b increases rapidly to a constant level 1.13. T₀ affects this plateau value of T_b however; if T₀ = 0, for example, the plateau value is very close to unity. The case when q = 25 corresponds to a position near the terminus. With h′ ≈ −0.06 the assumption |h′| << 1 begins to be in question and if integration is carried out sufficiently far this assumption fails entirely. The magnitude of H′ near a terminus is largely a result of the sliding law and initial parameters, mainly b. The wide boundary layer is evident. While T_xx increases, T_b also increases so that the profile is affected primarily by Tb rather than by a build-up of T_xx.

Non-Steady-State Illustration

In this section we present results for a glacier which is initially in steady state and which suddenly, by which we mean on a time-scale short compared to the characteristic creep slump time of the central region (Shoemaker, 1981), has its sliding resistance lowered in a central sub-region of dimensionless length Ξ. For convenience, we takes Ξ to be of order unity, except for extreme cases where H^⋆ is also of order unity; nothing is lost by such a choice. A reduced sliding resistance could result from water injection onto the base of the central region. The physical situation is closely analogous to the phenomenon of creep slump in central galcier reservoirs (Shoemaker, 1981).

To model this non-steady-state problem the steady-state equations were first solved using a uniform sliding law λ ≡ 1. The resulting H(ξ) was then substituted into the non-steady-state equations (19) and equations (20) with λ ≡ const. < 1 in the central region and unity outside. It is supposed that H(ξ) remains approximately steady for a period which is small in comparison with the slump time but sufficiently long for stress and velocity distributions to adjust to the new basal conditions. Figure 3 illustrates an example where α = b = T₀ = 0, q = 400, ū = 0.2, and λ = 1/5 on the central region 0.25 < ξ < 0.5. Here F represents a dimensionless flux HU/U(0), which is uniformly unity in the steady state since b = 0, and T_ph is a dimensionless form of the p-h basal shear stress given by the second of Equations (7) with y = 0 and α = 0:

(21)

Except for a transition zone adjoining the origin where it increases rapidly from zero, T_ph is almost unity. Note that from Equation (17) with b = 0, h′(0) = 0 results from the artificially imposed initial condition T₀ = 0. Because steady-state conditions apply on 0 < ξ < 0.25, F is unity there but decreases to the right. Because of the reduced sliding resistance in the central region, T_b, which is unity to four figures on 0 < ξ. < 0.24, changes almost discontinuously to another plateau value, constant to three figures, on 0.25 < ξ < 0.5 and then at ξ = 0.5 to a final plateau where it decreases slowly from unity. The behaviour of T_xx may be deduced from the appropriate |T_XX/T_bl curve of Figure 4. It is noted that T_xx changes sign within the boundary layer adjoining ξ = 0.25; that is, the flow changes from extensional to compressive.

The |T_XX/T_b| curves in Figure 4 are very nearly linear in the central region outside the boundary layer; the slopes of these lines are labelled in Tables II, III, and IV as T’ for various examples. Calculations also show that in jumping to a region of increases sliding resistances λ, the thickness of the transition boundary layer is negligible as compared with the thickness corresponding to a decrease of λ by the same factor. We find that the plateau values of |T_XX/T_b| to the right of the central region also increase with to the right of the central region and can be of order one or larger. From this, it can be seen that for a boundary-value problem where, in general, the entire glacier would be in a non-steady state, |T_XX/Tb| could be of order unity throughout, except near a zero of T_xx.

Table II. Effects of ū₀ on Deviatoric Stress in Non-Steady-State Situation. Values of T_xx/T_b at Different ξ for λ = 1/5 in Central Region, Ξ = 0.25, All For q = 400. T′ Values are (T_xx/T_b)′ (Nearly Constant) in Central Region Outside Boundary Layer.

Table III. Effect of q Upon Deviatoric Stresses, All Cases for λ = 1/2 in Central Region and ū₀ = 1/3. T′ Value for q = 25 is Not Meaningful Because of Boundary Layer

Table IV. Effects of λ Upon Deviatoric Stresses. All Cases for q = 200 and ū₀ = 1/2

The effects of parameter variations are shown in Tables II, III, and IV and are summarized below. In all cases the central region of reduced sliding resistance extends from 0.25 to 0.50; Ξ = 0.25.

• (i) Provided only that the length Ξ is sufficiently large, values of T_xx/Tb > 0.3 are encountered for practical ranges of parameters q, u_b(0)/U(0), and the jump in λ. ^*

• (ii) Table II shows the effect of ū₀ on (T_XX/T_b). With values given at ξ = 0.25, 0.5, and 0.75, as well as the average value of (T_xx/T_b)′ on 0.25 < ξ < 0.50 (shown as T′) outside the boundary layer, a sketch of each case can be made. In all cases λ = 1/5 and q = 400. T′ values increase with increasing ū₀, reflecting the fact that adjustment to the reduced sliding resistance by shear stress is less efficient at the higher ū₀ values so that higher values of T_xx are induced.

• (iii) Table III illustrates that a change in q does affect T′ For all cases λ = 1/2 and ū₀ = 1/3. (For q = 25, because H* ≈ 0.6 from Table I, the interval (0.25, 0.75) lies entirely within the boundary layer.) This result may be interpreted as follows: for fixed ū₀, λ, and h(0), if T_b(0) is changed there will be no change in T′ in adjusting to the jump in sliding resistance. (Note that H[ξ] changes as q varies.)

• (iv) Table IV shows that, as expected, decreasing λ produces rapidly increasing T_xx/Tb values. All cases are for q = 200 and ū₀ = 1/2.

Discussion

The model consists primarily of a set of assumptions which may be applied to produce equivalent models in situations not addressed here, e.g. those involving radial flow, non-sliding cold glaciers, other sliding laws, or other viscous constitutive relations such as a polynomial relation which exhibits finite viscosity at zero stress.

The plane-flow model developed here is a generalization of p-h theory designed to treat mean flow features when |T_XX/Tb| is of order unity, and reduces to p-h theory when |T_XX/Tb everywhere. It does, however, still incorporate an estimate of the integrated longitudinal deviatoric stress T_xx if the profile is determined by the simpler p-h relations based on the approximation |T_XX| < T_b, which can therefore be used to test the applicability of a p-h solution. If the estimated T_xx distribution satisfies |T_XX| T_B throughout the region of interest, then the p-h solution is valid, provided of course that the fundamental assumption of small surface slope, | h′ (ξ) |<< 1, is confirmed. When the estimated T_xx is not negligible in comparison with T_b, then the present model can be used to determine stresses, velocities, and profile which are compatible with a significant integrated longitudinal deviatoric stress T_XX.

It has been shown that under a wide range of conditions the p-h theory is valid for steady-state glaciers away from a margin or divide. However, examples have demonstrated that the p-h theory is not applicable to a particular class of non-steady problems in which |T_XX/T_b) becomes of order unity. It is conjectured that longitudinal deviatoric stresses are more generally important in transient problems, so that the estimation procedure for T_xx described above should always be applied before accepting p-h results.

Because the model employs averaging of various quantities over the thickness, so that the elliptic partial differential equations are approximated by ordinary differential equations, it cannot be applied to study the y-distribution of velocity or stress resulting, for example, from flow over or around an obstacle. Similarly, it is not expected that the description by the mdoel of the effect of abrupt changes in, for example, sliding resistance, will be accurate near the discontinutity. For example, the boundary layers exhibited here in the ordinary differential-equation solutions occur down-stream of a discontinuity, while an elliptic partial differential equation would produce a boundary layer on both sides of a discontinuity. On the other hand, the model should give a reasonable description of the global effects of such discontinuities over span lengths of order the thickness.

All of the results were obtained for short lengths not exceeding ten times the thickness at the starting place. It was decided not to obtain complete steady-state glacier profiles for several reasons. First, the phenomena investigated here do not require consideration of more than a short length. Second, it can be concluded from results presented that for reasonable parameter values the condition |T_XX/T_b| <<1 is satisfied everywhere except at the terminus or near an ice divide. Finally, the study of profiles near a terminus is a complex situation which requires more detailed treatment. We remark that parameters can be chosen such that the condition |h′|<< 1 is satisfied everywhere, while p-h theory requires a sliding law such that T_b 0 at the terminus to maintain bounded slope.

It is believed that greater than six-figure accuracy for T_b and H and at least two-figure accuracy for T_xx was obtained for all cases. The arithmetic operations involved in computing T_xx produce a loss in accuracy. In addition to the usual tests for accuracy, the steady-state and non-steady-state systems must produce identical results if a steady-state profile is used as input to the non-steady-state equations, allowing a cross-check between two distinct numerical problems.