Introduction

The influence of bed profile on the flow of an ice sheet, and in particular on its surface shape, was first analysed by Reference RobinRobin (1967). His approximate treatment relied on two assumptions: (a) that the bed slope relative to the horizontal, and in turn the surface slope, is small, and (b) that the mean shear stress σ _xy over the ice sheet thickness, is small compared to the mean deviatoric longitudinal stress where x, y are rectangular coordinates with x-axis horizontal. The role of each approximation was emphasized by Reference CollinsCollins (1968), in particular the use of assumption (b) to neglect the shear strain rate in comparison with the longitudinal strain rate, and hence the variation of the horizontal velocity a with height, which allow a simple expression for the mean longitudinal strain-rate in terms of surface accumulation and surface profile. Reference NyeNye (1959) argued that the temperature variation with height in the large cold ice sheets implies that the shear motion is confined largely to a thin “warmer” basal layer, supporting the approximation u=u(y), but also used the approximation in his calculation of an equilibrium profile, in contrast to assumption (b). Analyses to date using a temperatureindependent ice law, or ignoring temperature variation, conclude that is a smaller quantity than so a possible justification of the approximation (b) awaits a thermo-mechanical analysis. The strong inequality (b) appears unlikely. Robin’s approximations do not lead to an explicit equation for the surface profile when the bed profile and surface accumulation are given.

Reference NyeNye (1969) and Reference BuddBudd (1970[b]) also derive relations between the basal shear stress and mean longitudinal stress gradient by integrating the momentum equations through the (varying) thickness; Reference CollinsCollins (1968) obtained a similar result though referred to different axes. Reference BuddBudd’s (1970[b]) approximations of small mean bed inclination, small bed slope, and small surface slope recover the Reference RobinRobin (1967) terms together with a double integral of through the thickness, where the x-axis is along the mean bed line. The latter term is neglected by Reference RobinRobin’s (1967) assumption (b), but Reference BuddBudd (1970[b]) asserts that it is important for small bed wavelengths up to a few times the thickness and remarks that σ_xy is often much larger than σ’_x The shearstress variation appears to play no part in Budd’s subsequent estimate of a “mean viscosity through the thickness” in terms of the mean longitudinal stress and strain-rate. Noting his earlier arbitrary assumption about the depth variation of Reference BuddBudd (1970[a]) next introduced a stress function to satisfy the equilibrium equations for the stress perturbations about the solution for a uniform slab flowing down an inclined plane. For a constant-viscosity (Newtonian) fluid, the stress function is biharmonic, and sinusoidal solutions for an undulating bed of small amplitude can be constructed. When a non-linearly viscous ice law is considered there is no tractable equation for the stress function, but Budd proceeds by (implicitly) ignoring the perturbation in mean pressure to deduce that the stress function is harmonic; that is, by making an approximation where σ_zz is the stress (perturbation) normal to the plane of flow. There is no physical basis for such an approximation, so the validity of the resulting sinusoidal bed solutions (Reference BuddBudd 1970[a], Reference Budd1971) must be questioned. In fact, alternative calculations by Reference HutterHutter and others (1981)-abbreviated to HLS-imply that the approximation is not generally valid, and our later analysis demonstrates that it is invalid for most, if not all, practical conditions.

While the above approaches are based on assumptions that certain physical quantities are small enough to be neglected, they provide no systematic approximation of the full balance laws and boundary conditions. It is not clear that any of the results are valid approximations to the prescribed problems. In contrast, HLS and Hutter (1981)-abbreviated to H-have now investigated the effects of small bed undulations by series expansions in small (distinct) dimensionless parameters. Since we now propose an expansion procedure in a third parameter, extending the horizontal flat bed analysis of Reference MorlandMorland and Johnson (1980) -abbreviated to MJ-it is instructive to define the three approximation schemes and compare their ranges of application. Let the ice-sheet thickness and semi-length have magnitudes h ₀ and l respectively, and the bed undulation (or general profile) have amplitude a and wavelength (length scale) λ, so that (relative to a mean bed line) the mean surface slope ε and bed slope θ are given by

(1)

MJ have shown that the aspect ratio a is determined by the accumulation magnitude, the depth magnitude h _o, and the viscous properties of the ice, and that in most, if not all, practical situations s is small. We use ε as our expansion parameter and assume that the bed slope is not greater than the mean surface slope; thus

(2)

HLS consider the situation when the bed undulation amplitude is small compared with ice thickness, and in their expansion scheme assume that thus

(3)

ε_bis the expansion parameter. H introduces a longitudinal coordinate stretching parameter µ to describe the assumed “slow variation” of surface undulations, with surface wavelength subsequently related to bed wavelength in magnitude, again with bed amplitude small compared with ice thickness. µ is the expansion parameter and a proposed measure is

(4)

While the relations (2) allow a a≈h ₀, and hence a treatment of a typical isostatic bed profile, the inequalities (3) and (4) exclude this situation. In all cases the permitted bed slope θ is not greater than the small expansion parameter, necessarily much smaller in the H scheme (4), but the θ bound in scheme (2) is the pre-determined aspect ratio ε. In general,

(5)

so that θ≫ε is permitted in schemes(3) and (4) when the amplitude and wavelength are sufficiently “large” and “small” respectively. For such undulating bed profiles θ=ε _s with 0 < s <1, and our scheme must be modified by coordinate stretching with the parameter θ. In this case the lead-order surface profile is necessarily flat and not influenced by the accumulation distribution, as is clear from the MJ analysis. A global solution will require a distinct expansion procedure for the margins with appropriate matching, and possibly matching between basal layer and outer flow. The HLS analysis in fact supposes that accumulation has no lead-order effect and that the lead-order surface profile is flat, which implies ε_b ≫ ε Thus , which is the "short wavelength-large slope" undulation not treated by the present analysis.

HLS is a central-zone analysis, not completed to the margins. The first iteration on the flatsurface solution is expressed as an additive superposition of the effects of bed profile with zero accumulation and the effects of retaining a flat bed but with surface accumulation; this determines the surface slope to order ε_b. In the first problem a sinusoidal bed is adopted, and in both problems the surface perturbation is assumed to be sinusoidal. The mean profile with slope ε determined by a smooth accumulation distribution is lost, but an oscillatory accumulation is assumed to determine a perturbation slope of order ≫ ε. The main results of HSL are the prediction of surface effects due to the bed undulations. The solutions are necessarily numerical, but one conclusion is that the transfer of bed undulations to the surface depends strongly on the sliding law. While amplitude transfer is maximal for wavelengths of three to five times the thickness for small bed inclination χ and a particular range of sliding velocity for the power law of sliding, this is not a general feature as predicted by Budd. The HSL perturbation scheme presented fails when χ is too small ( < c. 0.01); that is when or bed drop over the ice-sheet length is not greater than the ice thickness, which must be the useful range for large ice sheets. Larger χ will have application to steep, thin, glaciers.

The H analysis is immediately restricted to wavelengths λ≫h ₀, and for there is a further restriction λ≪l(μ≫ε),; while the "large slope" case θ≫ε requires λ≪(a/h ₀)l The range of permitted wavelengths, amplitudes, and slopes is therefore very limited. For illustration a sinusoidal bed profile is adopted and it is assumed that the surface profile is the superposition of a known equilibrium profile for a flat bed and a small sinusoidal form of prescribed amplitude. Reference NyeNye’s (1959) equilibrium profile is adopted, which meets the slope magnitude bound only in a central zone, but, by the inequality (4), several undulation wavelengths span a distance much greater than d. Further, when there are many undulations over the ice-sheet span. λ≫l, by (4) μ≫ε and the lead-order central surface profile is again flat. The illustration is therefore of doubtful value. Flat-bed calculations were also performed with the given equilibrium profile for inclination χ over the range 10^-3 to 0.2 to evaluate the second-order stresses. However, the lead order equations explicitly assume χ= O(1), so only the range 10^-1 to 0.2 is appropriate.

We now present a global analysis incorporating the mean slope is determined from the accumulation, extending the MJ analysis to bed topography with slope θ bounded by ε and to arbitrary mean bed inclination χ. For steep glaciers with χ=O(1) the validity of the perturbation scheme imposes a restriction on the thickness when a general viscous ice law is adopted, but no restriction with the conventional models with a single response coefficient. The lead-order solution reduces to an ordinary differential equation for the ice thickness, showing that bed topography is fully transferred to the surface, in contrast with the associated flat-bed solution. Our analysis is not valid for “short wavelength-large slope” undulations, but these results suggest that maximum transfer of bed topography to the surface is probably found in steep. thin glaciers. When η=O(1) the gravity term in the longitudinal momentum equation dominates the pressure gradient and provides the balance with the shear stress gradient. For large ice sheets with the balance requires an extension of the MJ analysis with a horizontal flat bed, and the lead-order solution gives an ordinary differential equation for the surface profile in which the bed profile enters in distinct ways. However, the bed profile is significant to lead order only’ if its amplitude is not negligible compared to ice thickness, which is the case for beds formed (in part) by isostasy. If the bed amplitude is of order r: compared with ice thickness, then the contributions to the surface profile and stresses are necessarily much smaller than our lead-order solution, and we have not constructed the next terms of the series. They become useful only when surface slopes of order a compared to the mean slope are relevant to measurements and predictions.

The restriction on bed slope in Equation (2) is strictly θ=0(1). so formally our lead-order solution could be applied to the case θ≃ε^1/2 with a/h ₀≃ε^1/2 when λ≃h ₀ to describe a “short wavelength-large slope” situation and bed-profile effects. However, any regular perturbation scheme assumes that coefficients in subsequent terms remain of order unity. Because the third derivative of the lead-order stream function with respect to the scaled longitudinal coordinate arises in a higher-order longitudinal momentum balance, MJ propose that it must be bounded in the margin, which then implies it is everywhere of order unity in the analysis with a horizontal flat bed. Here we show that the stream function depends on the bed profile, and with the same proposal impose the stronger restriction

(6)

For example, if we allow a magnitude ε^-1/2 on the right-hand side. then with a/h ₀≃ε^1/2 for a lead order effect, which still implies the restriction λ≫h ₀ of the H scheme. That is, under the restriction (6), we cannot treat undulations with a length a few times the ice thickness.

Following MJ we assume a power law of sliding with coefficient depending on the ice thickness, equivalent to dependence on overburden pressure. The requirement of a small slope at the margins again imposes a restriction on the coefficient as ice thickness approaches zero. The accumulation/ablation at the free surface is a prescribed function of elevation (above some horizontal). For a linear sliding law on an inclined bed it is found that the slope at the upper margin is uniquely determined if there is ablation at that elevation, consistent with the flat-bed analysis, but there are two possible slopes (equilibrium profiles) if there is accumulation there. However, global solutions with small slope up to both margins do not appear to exist when the ice sheet is not symmetric about the ice divide. A general non-linearly viscous, incompressible f uid ice law for constant temperature is adopted (Reference GlenGlen, 1958; Reference MorlandMorland, 1979) with the assumption that dependence on the square of the strain-rate tensor is not greater than that on the strain rate. The latter alone is incorporated in the conventional model which we also adopt for illustration. The analysis treats exactly the apparent singularity at zero stress associated with Glen’s law with exponent n > 1. For all inclination χ the result

(7)

holds, contrary to Robin’s assumption. A corresponding analysis incorporating the significant temperature-dependent rate factor (Reference MorlandMorland, 1979) is now required to determine whether the strong inequality (7) is weakened (or inverted) by the typical variation of temperature with depth in a cold ice sheet. With the usual ice law the approximation (7) implies that the longitudinal strain-rate is much smaller than the shear strain-rate, but this would not follow if the term ignored in the general relation was in fact dominant. One-dimensional experiments cannot distinguish the terms (Reference MorlandMorland, 1979).

The Steady Ice Sheet

The dimensionless variables introduced by MJ are used from the outset. All length scales have unit h₀, an ice thickness magnitude. Figure 1 shows the ice sheet with mean bed line inclined at angle χ(–π/2<χ<π/2) to the horizontal; X, Y are rectangular coordinates along and normal to the mean bed line in the plane of motion, Y=H(X) is the free surface, Y=F(X) is the given bed profile. For numerical convenience the origin is chosen to be the left-hand margin, upper margin if χ > 0, with the mean bed line defined to pass through the origin; it does not, in general. pass through the right-hand margin. The dimensionless stress ∑ and pressure P have unit pgho where p is the ice density. Longitudinal and normal velocities U, V, vertical accumulation Q per unit horizontal cross-section, and normal basal drainage B, have unit q_m which is a magnitude of maximum accumulation/ablation/drainage. Immediately H, P, V, Q, B do not exceed order unity in magnitude. Let Z be the vertical elevation above the origin, then it is supposed that Q = Q(Z) and B = B(Z) are prescribed. The net accumulation ablation over the surface plus drainage over the base must vanish for a steady ice sheet.

Fig.1. The ice sheet and bed profiles.

Introducing a stream function ψ by

(8)

satisfies mass balance identically. The longitudinal and normal momentum equations for slow viscous flow with negligible inertia terms are

(9)

With the zero traction and accumulation conditions on the surface Y= H(X) are

(10)

(11)

(12)

Basal drainage and a power law of sliding between the tangential velocity Ub and shear traction on are given by

(13)

(14)

It is assumed that basal slip and shear traction are in the same direction, and that, m≽1. Dependence of the positive coefficient Λ on ice thickness is shown to be equivalent to dependence on normal bed pressure. Λ magnitude of A will be assumed in the analysis so that the sliding law influences the motion; perfect slip and non-slip are the limits Λ→0 and Λ→0∞. The margin conditions are

(15)

where X_M is the unknown ice sheet span. The zero global mass flux provides the additional condition to determine the unknown span.

A general incompressible viscous fluid ice law may be written (Reference GlenGlen, 1958; Reference MorlandMorland, 1979)

(16)

where D is the strain-rate tensor, A(T) a temperature dependent rate factor, and σ₀, D₀ are stress- and strain-rate units, chosen to normalize the response coefficients φ₁, φ₂. Here we adopt a temperature-independent model by taking a constant value for A corresponding to a temperature 250 K when A(273 K)= 1. Thus

(17)

The commonly adopted model sets χ₂=0,χ₁=χ₁(I₂), and inverts the relation (16); thus

(18)

We present the analysis for the general law (16). but adopt the common model for illustration. The reciprocal relations in the third and fourth of Equations (18) allow direct use of Glen’s law for temperature 273 K is

(19)

where n depends on the stress range chosen for the approximation. The singularity in 01and J₂ →0 represents an infinite viscosity at zero stress, and an alternative bounded approximation over the shear stress range 0→ 10⁵ N m^-2 is the Reference Colbeck and EvansColbeck and Evans(1973) polynomial law

(20)

ϕ₁ cannot be determined analytically, but is bounded, and ω is used directly in the solutions.

The MJ analysis for a horizontal flat bed shows that solutions for laws (19) and (20) are not significantly different. However, we again present the analysis to allow a singularity of the form (19) by setting

(21)

so that remains bounded. For Glen’s law ϕ₁ is constant (order unity), and the bounde ϕ₁case is recovered by setting α= 0(n- 1). We also assume that the stress term is not of greater magnitude than that of the ϕ₁ term, and that ϕ₂ is bounded; thus Hence the non-zero stress components for plane flow are

(22)

where

(23)

and for most, if not all, practical situations (MJ),

(24)

It is immediately clear that changes in ∑_ZZ are of the same magnitude as changes in when P is not negligible compared to the extra viscous stresses, and inequalities (24) suggest that P will always make a lead-order contribution even when velocity gradients are large, so Reference BuddBudd’s (1970[a], Reference Budd1971) harmonic stress-function approximation appears to have no practical validity.

Perturbation Method

If P, ψ, and appropriate gradients do not exceed unity in magnitude, then a series expansion in v gives the lead-order solution from Equations (9) and (11):

(25)

The X=0(1) solution is a static ice reservoir with horizontal surface resting against a steep embankment, and the X=0(l) solution is a slab of uniform depth on a bed with small mean inclination. Neither case is a leading approximation to a global solution with upper and lower margins. To avoid the surface result of the second option in Equation (25), it is necessary to scalethe longitudinal coordinate X by a small factor ε. so that the shear-stress gradient in the longitudinal momentum balance enters to lead order.

Following MJ, set

(26)

where ε, γ, ψ and all relevant derivatives with respect to ξ are of order unity (or less), together with

(27)

Assume β is order unity (or less)-strictly εβ=0(l). Thus the mean surface slope β has magnitude ε and the bed slope β is restricted to order β so that f′(ξ) does not exceed unity in magnitude. For later validity we assume also that f′′(ξ) and f′′′(ξ) do not exceed order unity.Footnote ^* The stream function scaling is given by the accumulation balance. The momentum equations become

(28)

and the corresponding stresses are

(29)

where

(30)

Suppose now that , then the leadingorderre lations are

(31)

The response coefficients are evaluated at and the Φ₂ contributions to

are retained for later estimations. Immediately the longitudinal deviatoric stress is smaller th an the shear stress ξ_xy ⁰ by a factor ε. Boundary conditions (10) to (14) become

(32)

(33)

where

(34)

is assumed to be of order unity for the influence of the sliding law. The distinct momentumbalances for X of order unity and X of order ε will be presented separately.

Finite Bed Inclination

For X = O(1) the momentum balances (28) with stress expressions (31) have terms of order

(35)

respectively. By the second set of magnitudes, and the relations (21), (23),

(36)

Then the fir st set requires that the balance

(37)

and the inequality (36) implies the h₀, Ø₂ restriction

(38)

This limits h₀ to the order of 10 m if Ø₂ = O(1), but the small-slope solution is valid for thicker glaciers as Ø₂ decreases, with no limit for the usual model Ø₂ Ξ 0. The Ø₂ contribution to ∑_xx ⁰,∑_yy ⁰ is of order unity when the equality in (38) holds. The mean surface slope or aspect ratio ε:= vⁿ is smaller than ε=v^n/(n+1) obtained for the horizontal flat bed (MJ) and for x=O(ε). Set

(39)

and allow the equality in (38) to hold. The lead-order solution of the normal balance given by thesecond of Equations (28) subject to the surface condition given by the second of Equations (32)is

(40)

and of the longitudinal balance given by the first of Equations (28) subject to surface condition given by the first of Equations (32) is

(41)

so ∑_xy ⁰ is of order unity and is everywhere directed down the mean bed line. Hence σ_xy is of thesame order as the overburden pressure. Thus U_b always has positive component down the meanbed line, so ξ= sgn X, and there is no ice divide. Hence (∑_xx - ∑_yy )⁰ is of order ε in contrast toorder ε² for the horizontal flat bed (and small inclination bed), but is still a factor ε; smaller than ∑_xy ⁰ Using the definitions given in Equations (30), (23), (21), and (31), i₀ = the left -handside of relation (41) becomes . To proceed analytically we need an inversionin terms of J₂ at J₃ =

. We will assume the third and fourth of relations (18)with ξ₂=0, and adopt the form (19) or (20) for further analysis. Then the relation (41), since ε/δ=(10 m/h ₀) ⁿ , gives

(42)

or equivalently

(43)

where for Glen’s law, Equation (19),

(44)

and for Colbeck and Evans’ law, Equation (20),

(45)

Now Equation (44) defines g(t) of order unity since the argument of g in Equation (43) is oforder unity, independent of the magnitude of h₀, while Equation (45) implies g(t) ξ I if h₀;> 10 m.However, h₀ ξ 10 m corresponds to a stress unit of 105 N m ^-2, which is therefore the magnitude of both p and σ_Xy, and is also the limit of the polynomial law of Equation (20). Thus Equation(45) is not applicable to deviatoric stresses much greater than 10⁵ N ^m-2; that is for h₀ ;> 10 m. An alternative bounded viscosity law is therefore required to describe the response over a greaterstress range corresponding to larger ho when the bed inclination is finite. Thus any realrestriction on ho depends on the magnitude of Ø₂ through the bound (38). We use Equation (45)with h₀ = 10 m for illustration. Defining

(46)

integration of Equation (43) yields

(47)

The basal drainage and sliding laws given by the first two of Equations (33) give

(48)

where

(49)

The vertical elevation of both base and surface is approximated by Ȥ to lead order since themean bed drop dominates over thickness variation when x= O(1). Finally, the accumulationcondition given by the third of Equations (32) gives a first-order differential equation for thethickness η(ξ) of the form

(50)

subject to end conditions

(51)

The second condition determines the unknown range ξ_M, but in fact this is pre-determined bythe condition of zero mass flux

(52)

which is compatible with the perfect differential form provided that 1 + m( 1 - s) > 0 where .A stronger restriction is guaranteed bybounding e₁(ξ) necessary for bounded velocities. Assuming that Q^* does not vanish at bothmargins and is sufficiently smooth, the perfect differenti al shows that and hence we require s = 1. Then

(53)

and hence Q_M ^* (net accumulation and drainage) is positive at an upper margin and negative at alower margin; the roles of ξ;=0 and ξ = ξ_M are reversed when X changes sign. This analyticbehaviour ensures that derivatives of ψ₀ up to the third, which arise in the next-order momentumbalances, are bounded provided that e^''' ₁ is bounded; that is, g^"(t) is bounded as 1→0. For Glen’slaw, Equation (44), this requires n=1 or n ξ 2. Recall thatJ" and f” were assumed to be oforder unity, and hence η^” and η^“” are also of order unity.

Since the thickness η(ξ) is determined independently of the bed profile f(ξ), the latter issimply superposed on η(ξ) to determine the surface profile η(ξ) that is, there is full transmissionof bed profile to surface profile relative to the corresponding flat -bed solution (not relative to aslab of uniform depth). As noted earlier, the present solution does not apply to “shortwavelength-large slope” undulations. The velocity distributions involve both f(ξ) and η(ξ)independently. The basal shear traction is simply

(54)

By the asymptotic results (53), the marginal slopes satisfy the explicit relation

(55)

We now note the main features for a linear sliding law (m = 1) with coefficient and accumulation linear in elevation,

(56)

where q₁ = 0(1). This form supposes that Q^* varies no more than of order unity over the beddrop. Immediately, for all ice laws g(t),

(57)

Small Bed Inclination

Now consider χ= εχ₀, χ₀ = 0(1), which is a realistic situation for large ice sheets. Of the momentum term magnitudes displayed for χ= 0(1) in the sequence (35), only the final sin χ term of the longitudinal balance changes, from order unity to order ε. A lead -order shear stress gradient then implies that

(58)

and is identical to the estimate of aspect ratio for the horizontal flat bed (MJ). The restriction (36)

(59)

The factor multiplying εφ₂, independent of h_o, is the term θ^{(2α+1) (2–2α)} of the MJ analysis, and for all practical values does not exceed order unity. Thus with the assumption |ϕ ₂⩽1|, the φ₂ ≤ 1 terms in both longitudinal and normal momentum balances are smaller than tile lead-order terms, ԑ and unity respecti vely, by a factor ԑ at least, and P=P₀. The lead-order normal balance with surface conditions given by the second of Equation (32) again gives Equation (40), but now the pressure gradient of order ԑ enters the longitudinal balance, so with the surface condition given by the fir st of Equation (32),

(60)

Thus ∑ = O(ε) while from Equations (31) and (58), . Here cos χ = 1, sin χ = εχ₀. Thus ∑⁰ _xy=0 at position(s) that is, where thesurface is horizontal. Clearly U_b = 0 at ξ=ξ_d from the sliding law, and we show later that U=0 at ξ = ξ_d which is therefore an ice-divide–there is no flow across ξ=ξ_d.

The inverse relations (42) and (45) are again obtained with sin X replaced by X_o - γo and (h₀/ 10 m) replaced by (ε_ho/10 m) = θ^1/(n+1), since here The deviatoric stress limit in (45) now implies that, which only requires. Hence

(61)

Substituting the first of Equations (60) and the first and third of Equations (61) in the sliding law(33)2 and drainage law given by the second and first of Equation (33) yields

(62)

and hence a longitudinal velocity

(63)

Since

For X₀ - Y₀ = or < 0, confirming that ξ = ξ_d is an ice divide.

Finally, the accumulation condition, the third of Equation (32), gives the differential equation

(64)

(65)

These results are analogous to equations (69) and (70) of the MJ horizontal flat -bed analysis, with n - f replacing, but Y₀ - X₀, and not Y₀ - β, replacing Yo. Thus the second-orderdifferential equation (64) for n(E) cannot be expressed simply as an equation for the thickness η= η-f, and so the bed profile is not simply superposed on a corresponding flat -bed surface profile.

For bounded e₁(ξ) at the margins it is again necessary that , which is sufficient to ensure a zero mass-flux integral of the differential equation (64). Since the sliding law is independent of bed inclination and profile solutions, the differential equation mustgive bounded slopes γ_M for any (non-zero) χ₀ - γ_M at a margin. γ_M = χ₀ at both margins is the special case of horizontal surface at both margins. In general the net accumulation Q-B = Q_M ^* at a margin (ξ=0 or ξ=ξ_M is not zero, and the sliding law is independent of the accumulation, so the leading-order balance requires (η–f)^m(1–s)(γ,₀γ–β to be of order unity at the margins, and hence s = 1 for bounded γ_M’ Then

(66)

Thus a horizontal margin surface γ_M = χ₀, or surface tangential to the bed γ_M = β_M, imply Q^* _M = 0, and vanishing Q^* _M implies one of these situations. In general

(67)

and γ_M =β_M > or > 0 at the left- and right-hand margins respectively.

Given the end conditions η-f=0 and η’ evaluated at ξ=0 (left-hand margin say) by Equation (66), the differential quation can be integrated numerically until η-1 aξain vanishes, determining ξ = e;M. We should note that for a given set of prescribed values of Q_M, P_M, X₀, the algebraic equation (66) may have no real solution Y/M, may have a unique solution, or may have more than one solution. This can be illustrated for the linear sliding law m = I when

(68)

For ablation, Q< 0 at the margins, and the above inequalities imply that 2_YM > X₀ + B_M and YM > X₀ at a left-hand margin, and that 2γ_M<χ₀ + β_M and γ_M<χ₀ at a right-hand margin. Hence the positive root and negative root in Equation (68) apply at the left- and right-hand margins respectively; that is, the surface slopes above the horizontal from either margin. For accumulation, a solution requires that at a lefthand margin and χ₀<γ _M <β _M at a right-hand margin; that is, the bed and surface both slope below the horizontal from either margin. Under these conditions both roots give compatible solutions, and hence two possible equilibrium profiles. We have not been able to construct a global analysis of the differential equation to confirm or reject either root. In practice, for this small mean bed inclination, we would expect ablation/drainage to occur at margins, and hence a unique margin slope when in= 1.

The above criteria were concerned only with bed and surface slopes relative to the horizontal, independent of the value or sign of the mean bed inclination Xo. Clearly, for X = εX₀ = 0(ε), a horizontal mean bed line may be chosen and the differential equation (64) referred to coordinates . To lead order (cos X= 1, sin X=ξX₀), ξ =ξ, and Equation (64) is recovered for il(k) since f=f is the profile superposed on the mean bed line.

It remains to consider possible singular behaviour at a margin n -f’z 0 and an ice divide contributes to the shear stress gradient in the next order momentum balance, following MJ we propose that it must remain of order unity. From Equations (61) and (62). the anticipated restriction

is necessary. As γ₀→χ₀ at an ice divide, η→η_d≠f, and Since a given ice law and given sliding law must allow valid solution for any sensible Q and B near the divide, we can consider a variation with determine any required restrictions on the laws. This explanation was omitted from the MJ discussion. The essential behaviour of the differential equation (64) as is

(69)

where c is a positive constant of order unity. Thus

(70)

and . That is, there must be accumulation at a convex ice divide and ablation at a concave ice divide.

The dominant terms of

are

(71)

Thus , we find and (of order unity) at an ice divide where requires l < 1 ; that is, l = 1. Direct differentiation of the separate dominant terms shows that bounded or 11 > 3, as in the MJ analysis. Since at a margin then conditions (i) and (ii) ensure that 17’ = YM, 1]", 1]'" are all bounded (of order unity), and in turn is of order unity.

Once n(ξ) is determined from the differential Equation (64), the lead-order stresses are given by Equation (40) and Equation (60), and velocities by Equation (61), Equation (62).

It is instructive to re-express the lead-order Equation (64) in the physical variables which were the starting point of the MJ analysis, namely

at a margin, where h is the surface height normal to the bed, z is the height above a horizontal, λ(h - J) is the sliding coefficient with λ ~ λ; (h - f) at a margin, with bed profile I referred to the physical length scale, and q is the accumulation with value - qo (ablation) at the margin. Thus, for given q(z ), the differential equation and starting slope are independent of the height scale ho, and therefore the profile is independent of the magnitude ho. That is, ther’e is a unique profile and centre height, small or large, governed by a length scale defined by the variation of q with height, though the aspect ratio

and hence the span 2/, depend explicitly on the magnitude. Our dimensionless illustrations are computed by prescribing the normalized Q(Z) on a length scale ho, which, if q(z) is given, implies that a solution with centre height ho exists. Since ho is determined by q(z), solutions cannot exist for arbitrary Q(Z) and in particular the variation must allow a balance of accumulation and ablation.

Illustrations and Conclusions

The following illustrations all use the Colbeck and Evans’ law, Equation (45), with the factor h₀/10 m replaced by (εh₀/10 m) = θ¹² = 0.3 for direct comparison with the MJ results for ahorizontal flat bed. The latter show that results for Glen’s law, Equation (44), are notsignificantly different. A linear sliding law (m = 1) with coefficient is used. The majority ofsolutions are computed for accumulation linear in height and zero basal drainage:

(72)

where Q₁ = 1 + Q₀, as in MJ, to set Q( 1) = 1. Q₀ >0 represents ablation at Z =0 and in turn aunique slope at a margin at Z = 0. Effects of bed inclination and bed topography are presentedseparately.

(i) Inclined flat bed

Choosing X₀ = 0, so that Z = Yand

(73)

describes a flat bed inclined at angle ξX₁ (β = - X₁) to the horizontal. A margin at is the upper margin if X₁ > 0, the lower margin if X₁ < 0. We set Z = Y =0 at E = 0, so all profiles are described with respect to common horizontal and vertical axes OξY. Solutions for X₁ in the range (horizontal bed) to 5 for various Qo have been constructed, integrating the profile given by Equation (64) from the upper margin ξ = 0. Further computations have been made with Q₁ ≠ 1 + Q₀ and also for quadratic and exponential Q(Z). In all cases with ablation at the upper margin, Qo > 0, the surface first rises into an accumulation zone, but finally approaches a horizontal, not returning to the (receding) bed at a finite span. Thus the solution approaches an ice reservoir with horizontal surface over which there is uniform accumulation, balanced by decreasing horizontal outflow over increasing depth, illustrated in Figure 2. For Q₀ < 0 in the distribution given by Equation (72), with accumulation at the upper margin, with either of the negative starting slopes given by Equation (68), the surface slope increases to zero (a concave ice divide) where the surface has dropped to an ablation level. The surface then appears to remain horizontal, not returning to the bed, also illustrated in Figure 2.

To avoid the reservoir solution, the computations were repeated for - X₁ in the range to 5,so that the starting point E =0 is the lower margin, at which there is ablation (Qo >0). In all casestreated, the surface rises to an ice divide at an accumulation level at which mass balance is(Fig.2.)attained, and then descends towards the rising bed, but while Yo remains order unity, it is foundthat - η" increases dramatically as η - J approaches zero. This violates the assumption of orderonederivatives, since it" enters the differential equation (64), and implies that a global solutionwith small slope at the upper margin is not possible. As noted above, the solution with prescribedsmall slope at an ablating upper margin approaches an unbounded reservoir. For Xl = 0(horizontal bed), the profile is symmetric with order unity it" everywhere, but as -Xl is increased, the zone of large - it" at the upper margin increases. The solutions appear reasonable close to themargin for - Xl ξ 0.1, but with increasing - Xl a contradictory situation of no ablation beyondthe ice divide is reached e.g. - Xl = 0.5, Q(~m ξ 1.3) = 0.3 > O. In such situations M is close to Figure 3 shows a comparison between Xl = 0 and Xl = - 0.5 with Qo = 1, Ql = 2, with the zoneof large - it" indicated by the dashed lin e. This situation was repeated for all the alternative Q(Z).

Fig.2. Ice reservoir solutions.

Since mass balance is attained at the ice divide, zero net flux must also be attained in the upper zone beyond the divide. Here, the rising bed and Q increasing with elevation requires the surface to steepen to reduce the accumulation at the upper levels and increase the ablation zone, provided that the rise of the bed does not eliminate the ablation zone too quickly. A validsolution up to the upper margin will therefore require an order-one slope (at least) in some marginal zone, so the full equ ations must be solved there and matched to the small-slope solution. However, since the small-slope solution is uniquely determined by lower margin starting conditions, it is valid until - it becomes large, and so is not influenced by a passive boundary-layer solution necessary near the upper margin. Thus the small-slope solution provides (Fig.3. and Fig.4.) initial conditions at an appropriate point ξ to continue the full equations through to an upper margin which is part of the solution.

Fig.3. Surface profile for horizontal and inclined flat beds with Q₀ = 1, Q₁ = 2.

Fig.4. Surface profile for Q =-1 +2Z_Y0.

A modified form of Q,

(74)

was used to demonstrate that a valid global solution can be obtained. Since accumulation is required at a convex ice divide, Yo = 0, the form of Equation (74) implies that an ice divide can not occur. Thus γ₀> 0; everywhere, and the distinct mass balances in lower and upper regions which caused previous difficulties, do not arise. Mass balance now determines a span ξm = 2QoIQI xf. Figure 4 shows a valid small slope solution for Q₀ = 1, Q₁ = 2, X₁ = 1.

(ii) Topography on a horizontal mean bed

For a flat horizontal bed, f = 0, X₀ = 0, with Q = Q(Z) the ice sheet is symmetric about the ice divide. Recall that the ice-divide height η_d and position ξ_d are determined uniquely by the secondorder differential Equation (64) with prescribed margin η(0)=0 and margin slope η'(O) =YM’ In the symmetric case the same centre height is reached by integrating from either margin.

Fig.5. Surface perturbation Δ induced by bed undulations.

Suppose the bed is perturbed at one side of the divide only, then the height of the ice divide reached by integrating from the margin at that side is perturbed from the ice-divide height attained by integrating over the flat -bed side. Provided that the perturbation in 1Jd is of order e, continuity of the lead order n is satisfied, and it is anticipated that the bed perturbation only influences the next term in the series for n. However, if the bed topography induces a lead-order effect, then (in non-symmetric cases) the unique ice-divide heights attained by integrating from the two margins with prescribed margin slopes will in general differ. Thus we can anticipate that the solution of the small-slope Equation (64) starting from either margin will again lead to large curvature near the far margin in an attempt to preserve mass balance. This occurs in all the computed solutions, with n" becoming large through negative values in some cases as n - f → 0 analogous to the inclined -bed results and also large through positive values when n - f remains positive and the profile opens up again. In a class of single-hump perturbations it was found that

(75)

Solutions for periodic bed undulations of the form

(76)

starting from a margin at E = 0 have been computed for various amplitudes f₀, phase shifts f₂, and integers f₁ so that the wavelength is an integer multiple of the semi-length E_M/2 = 1.54 for the flat bed profile. Since the ma rgin is not, in general, in phase with E=0, large curvature occurs there and there is a small marginal zone in which the small-slope solution fails. The choice

, Imples representing two and four full waves over the flat bed span. Figure 5 shows in the distribution given by Equation (72), but as with all the other examples, there is no clear indication of the transmission of undulating topography to the surface. The large negative value of L) at the far margin, indicated by the dashed continuation of the curve is a consequence of the profile reappearing as described earlier. We can only expect to exhibit a correlation by anal ysing the many-wave bed form, which requires an analysis valid without the restriction (6),

and is necessarily an t:-order effect since

for wavelength A approximately equal to ho. We are now developing an appropriate perturbation scheme.

Finally, the isostasy relation for a bed of highly viscous fluid (Reference WeertmanWeertman 1973):

(77)

where ρ_i, ρ_b are the ice and bed densities, is investigated with Weertman’s value ρ_i/ρ_b =1. The differential equation (64) now determines n (or f ) directly, and n is clearly symmetric about the ice divide so no far-margin problems arise. Figure 6 shows a comparison between the flat-bed solution, surface n. and surface n, and ice thickness n - f given by Equation (77), with Q₀ = 1 In Equation (72).

Fig.6. Comparison of the flat-bed projile η_f the isostalic profile η, and ice thicklless η-f.