Steady flow of a viscous ice stream across a no-slip/free-slip transition at the bed

Abstract

An exact, analytic description of the steady, downhill flow of a viscous ice stream across a no-slip to free-slip transition in basal boundary condition suggests ways in which such transitions might be detected by observations of ice-stream surface properties. We find that the best expression of this transition is in the free-surface profile, which dips over the point of transition and becomes horizontal far downstream. The streamwise velocity at the surface shows a gradual change across the transition, and this is in disagreement with previous studies which suggest a marked change. Basal shear stress and basal pressure exhibit singularities at the point of transition. As concluded in previous study, the former singularity is a likely point of strong basal erosion. The latter singularity is problematical, however, because it violates the assumptions which make the exact solution possible. It is not clear how this problem can be overcome without appealing to thermodynamics, non-Newtonian material properties, cavitation or wake separation.

1. Introduction

We examine the flow of an idealized ice stream across a no-slip to free-slip transition in basal boundary condition. Our goal is to determine the expression of this transition on the velocity and elevation of the upper ice surface; and, in so doing, to examine ways in which observations of the ice surface might diagnose the basal conditions. Our interest is motivated by the abrupt sliding velocity gradients observed on some ice streams, which drain Antarctica and Greenland, and on several mountain glaciers. These abrupt gradients can stem from a variety of causes generally associated with conditions at the bed. Grounding lines, for example, present a situation where the basal traction of an ice stream changes abruptly. Other examples include transitions between frozen and wet basal conditions, and possible gradients in the properties of deformable subglacial sediments. An understanding of the surface expressions of these basal conditions may help in the interpretation of ice-stream surface properties, which are becoming more easily observed as a result of satellite technology (e.g. Bindschadler and Scambos, 1991). Other motivations for our study include (1) the question of how stress concentrations, which may affect bedrock erosion and pressure-melting temperatures, might be generated by basal sliding transitions (Hutter and Olunloyo, 1981), and (2) the question of how the age profile of ice cores might be disturbed by the flow effects of basal sliding transitions (Weertman, 1976).

In brief, our problem concerns the steady two-dimensional stream of ice, treated as an incompressible fluid of constant viscosity, down an inclined plane of infinite length. At a single point on this plane, a switch from a no-slip to a free-slip basal boundary condition is assumed. Upstream of this point, constant thickness and steady Poiseuille flow (Batchelor, 1967, p. 183) are assumed. Far downstream, steady plug flow consistent with the free-slip condition at the bed is assumed, but the thickness is undetermined.

While this flow geometry and the viscous rheology are largely unrealistic, they allow an exact analytic solution of the governing equations. This solution is valuable because it may capture some of the essential features of the more complicated phenomena encountered in Nature. As easily appreciated after examining the mathematical discussion to follow, relaxation of these strict simplifications may render the problem intractable by analytic means alone.

Our problem is similar to that examined by Hutter and Olunloyo (1980) who were interested in the transition between melted and frozen beds of mountain glaciers. The essential result of their study is that shear and normal stresses at the bed are concentrated at the transition, and are in fact singular. This concentration may account for the efficient erosion of the glacial substratum suggested by observations (Hutter and Olunloyo, 1981). Our problem differs from that posed by Hutter and Olunloyo in that we allow for a free surface which might reveal a signature of the changes in basal conditions. Our study also differs in certain details of the method used to solve the governing equations. Hutter and Olunloyo (1980) used the Wiener-Hopf method to solve the equations, but made an approximation to simplify the analysis. We dispense with this approximation by following the mathematical analysis of Richardson (1970) concerning a low Reynolds number jet. (We are indebted to Professor K. Hutter for bringing this reference to our attention.)

To highlight physical insights which can be appreciated without resorting to lengthy, and somewhat specialized, mathematical development, we determine first several a priori results from the governing equations. Mechanical energy considerations permit us to identify the algebraic nature of the singularities of the stream function, vorticity and pressure which we encounter later in the exact solution. Of particular interest in light of previous research is that singularities in any variable can only occur at the point on the basal plane where the no-slip/free-slip transition occurs. We additionally establish by means of the maximum principle for harmonic functions that the vorticity of the flow is everywhere positive. As a consequence, recirculation, which might upset the chronological sequence of ice strata, is forbidden.

2. Mathematical Model

Following Hutter and Olunloyo (1980), we consider the simplest problem which embodies an abrupt change in the basal friction of an ice stream. To that effect, we consider the flow of a layer of viscous fluid overlying a boundary which changes abruptly from no-slip to stress-free. The upper surface of the layer is assumed stress-free, and can adopt a shape consistent with the conditions at the base. To simplify the problem to the utmost, we restrict our attention to the two-dimensional, steady-state case. Furthermore, we shall assume that the fluid is incompressible and homogeneous. The driving mechanism will be provided by gravity; in other words, the layer of fluid flows down an inclined plane.

If we pick the x-axis to lie along the bottom boundary and the z-axis normal to the latter, then the steady Navier-Stokes equations read (Batchelor, 1967, p. 147).

(2.1)

In the above equation, u and w are the velocity components in the x and z directions, p is the pressure, ρ is the constant density of the fluid, μ is its constant viscosity, g is the gravitational constant and a is the angle between the basal boundary and the horizontal. Subscripts represent partial differentiation with respect to the subscripted variable.

As mentioned earlier, we assume that one part of the bottom boundary, say x < 0, is rigid and the other part, say x > 0, is stress-free. As a result, the boundary conditions along the bottom are:

. (2.2)

At the upper surface, which has the following representation

(2.3)

we require the stresses to be continuous, i.e. to have the same values as those in the air which we assume to be at rest. Thus, if Π is the constant pressure in the air at the surface interface, then

. (2.4)

The formulation of the mathematical model is completed with the addition of the kinematic boundary condition which makes the assumption that there is no mass exchange with the atmosphere at the free surface, viz.

. (2.5)

Note that, if the bottom boundary were rigid throughout, the flow would be a simple Poiseuille flow. We shall assume that this is the flow which exists far upstream.

In spite of all the idealizations made, the above problem is too complicated for an analytical treatment. We shall therefore simplify the formulation still further by assuming that the slope α is small. This will enable us to do a perturbation analysis in α, which will yield a simpler, linearized version of the problem.

3. Small Slope Approximation

Series expansion of variables

Since u,w,p and h are functions of α, for the small α case, we look for a solution of the above problem as a Taylor series in a, viz.

(3.1)

where superscripts designate the order of the various terms in the expansion. The velocity fields have no zeroth-order terms because there is no motion when the slope is nil. One must remember, however, that even when the slope is nil there is a pressure field and the free surface has a specific form; hence there is a need for zeroth-order terms in the series expansion for p and h.

We substitute these series representations into Equations (2.1)–(2.5) in order to obtain a sequence of boundary-value problems for the fields of various orders. We start with the boundary conditions (2.4) which require the greatest care. They read

(3.2)

The typical field, say ϕ ⁽ⁿ⁾ entering in the above expressions is evaluated at ; thus, it can be Taylor expanded again in powers of α. In other words,

(3.3)

Therefore, condition (3.2) is equivalent to

(3.4)

Zeroth-order Problem

We are now in a position to consider separately the zeroth-and first-order problems. The zeroth-order problem is simple because it does not involve the flow:

(3.5)

with

. (3.6)

The solution to Equations (3.5) and (3.6), which can be found by inspection, implies that the upper surface is parallel to the bottom and the pressure is hydrostatic

. (3.7)

First-order Problem

To determine the flow, we must consider the first-order problem, which reads

(3.8)

with

. (3.9)

The kinematic boundary condition (2.5) requires the vertical velocity to vanish at the upper surface

. (3.10)

Finally, the bottom boundary conditions (2.2) imply that

. (3.11)

Dimensionless Variables

The above development completes the formulation of the first-order problem on which we shall now focus. We rewrite this problem in a slightly simpler form by introducing dimensionless variables. Denoting temporarily these dimensionless variables by a prime, we define

(3.12)

Following the usual practice, we substitute the above expressions in the first-order problem (3.8)–(3.11), and then drop the primes as well as the superscripts. This yields (Fig. 1):

(3.13)

with

(3.14)

and

(3.15)

Once the solution to this problem is found, the deflection of the upper surface is found by the second dynamical condition in (3.9), viz.

. (3.16)

Fig. 1. Steady flow of a viscous ice stream down an inclined plane (physical picture) exhibits a transition from Poiseuille flow to plug flow as a result of a transition from a no-slip to a free-slip basal boundary condition at x = 0. An exact, analytic description of this idealized flow is obtained here by solving the problem shown in the mathematical picture.

4. Preliminary Results

The exact solution to the first-order problem requires a lengthy mathematical analysis, which some readers may want to skip. In order to enable all readers to get the gist of the results, we shall first deduce several fundamental properties of the exact solution, using only a modest knowledge of the mathematical properties of harmonic functions (i.e. solutions to Laplace’s equation). More particularly, we shall exhibit the nature of the singularities of the vorticity, basal stress and pressure; we shall then prove that the vorticity is non-negative and finally, rule out the possibility of recirculations, i.e. of closed stream lines which do not emanate from upstream. These a priori results complement the exact solution, once it is found, by reassuring us that the algebra and numerical calculations are correct.

Boundary conditions at x = ±∞

We start by introducing a stream function ψ such that

. (4.1)

Then, if η denotes the vorticity, i.e.

. (4.2)

we can rewrite the Equations (3.13) as follows

(4.3)

or equivalently

. Far upstream, we want the flow to tend to a Poiseuille flow. (The general Poiseuille flow has been discussed by Batchelor (1967, p. 182.)) This flow arises from the balance between the downstream component of the gravitational force and the shear stress at the basal boundary. Thus, we require that as x→–∞

. (4.4)

With the above trends, we have arbitrarily chosen the upper boundary to be the stream line ψ = 0 and the lower boundary to be the stream line ψ = −⅓. With these choices, the boundary conditions (3.14)–(3.15) become

(4.5)

and

. (4.6)

Far downstream, we anticipate a plug flow. (Here, we define plug flow as being a laminar flow which lacks shear across planes parallel to the bed.) As a result, the viscous forces become vanishingly small, and the pressure must rise to balance the component of the gravitational force in the downstream direction. Thus,

(4.7)

where C is a constant. Note that since the x-coordinate is parallel to the basal plane, which is inclined with respect to the horizontal, the total height h ⁰ + αh ¹ tends to a constant C as x→∞.

Consequences of finite viscous energy dissipation at the origin

We next turn our attention to the behaviour of the fluid in the immediate vicinity of the change from the no-slip to the stress-free boundary. Since derivatives of ψ have different values at x = z = 0 depending on how this point is approached, we anticipate that ψ will be singular at this point. Therefore, we represent the stream function there by a series which allows for the possibility of an algebraic singularity, namely

(4.8)

where r and θ are the polar coordinates, v is an unknown exponent and ψ is an analytic function at r = 0. This approach, which follows that of Richardson (1970), has been widely used in applied mathematics (see e.g. Batchelor, 1967, p. 226; Carrier and Pearson, 1976, p. 66). As a result, we can write this function as follows

. (4.9)

We should emphasize that since v is yet to be determined, there is no loss of generality in assuming that . If we substitute this representation of ψ in the biharmonic equation for the stream function, we see that

.

Consequently

. (4.10)

On account of the boundary condition (4.6), , and consequently from Equation (4.8)

(4.11)

and as a result

(4.12)

and

. (4.13)

Furthermore, the stress-free condition on θ = 0 implies that and hence

. (4.14)

Finally, the no-slip condition for θ = π yields

. (4.15)

We can check that v = 1 is not acceptable since it implies that and , which after substitution in Equation (4.10), leads to . As a result, from the two independent, linear, homogeneous Equations (4.13) and (4.14), we deduce that

(4.16)

Therefore, we are left with either

(4.17)

or

(4.18)

where β has yet to be determined. In either case, the viscous energy dissipation per unit mass Φ, near the origin is of . Therefore, if the total dissipation is to remain finite in this region, i.e. , then r ^2v–3 should be integrable over a region surrounding the origin. In other words, from energy considerations, we want to focus on solutions of Equations (4.17) or Equations (4.18) such that

i.e.

(4.19)

The smallest ν which satisfies condition (4.19) and either Equations (4.17) or Equations (4.18) is

(4.20)

As a result, we conclude that near the origin, we have

(4.21)

Note that the vorticity η is of one sign near the origin and that lines of constant η fan out of the origin. A look ahead at Figure 9 shows that the exact solution does indeed display this behavior.

We consider next the singularity of the pressure. As is the case for all Stokes flows (Batchelor, 1967, p. 229), the pressure is a harmonic function, namely

. (4.22)

This equation is a direct consequence of taking the divergence of the momentum equations in Equations (3.13) and using the continuity equation. The vorticity η, which is also a harmonic function, is closely related to the pressure. To see this, let us denote by ζ the harmonic conjugate to η. In other words, η and ζ are the real and imaginary parts of an analytic function f of the complex variable x + iz. As such, η and ζ satisfy the Cauchy-Riemann conditions, viz.

(4.23)

which allow us to write the momentum equations in Equations (3.13) as

(4.24)

As a result, we conclude that

(4.25)

where C is a constant. By using the above equation and the expression for η given in (4.21), we conclude that near the origin

. (4.26)

The existence of a singularity in the pressure is signified by the factor .

The analysis used to deduce the form of the various fields near the origin can also be used to establish that this point is the only one where these fields are singular. This result casts doubt on the accuracy of the streamwise velocity field displayed by Hutter and Olunloyo (1980) in figure 4 of their paper.

This figure shows the streamwise velocity as having strong peaks near the origin at positions above the basal boundary. While these peaks are not singularities per se, our results suggest that the transition from Poiseuille to plug flow should be gradual.

Consequences of Laplace’s equation

By combining the information we now have about the vorticity, we show next that the vorticity is positive throughout the fluid region. We shall then use this result to establish that there can be no recirculation in the flow. As already mentioned, we can prove these assertions without resorting to the explicit form of the vorticity or the stream function.

We recall that the vorticity 77 is a harmonic function in the strip 0 < z < 1 which, according to conditions (4.4)–(4.7), takes non-negative values on all the boundaries, with the possible exception of x < 0, z = 0 where its value is not known. From the maximum principle (Protter and Weinberger, 1967, p. 64), it follows that if η were positive on this no-slip half of the lower boundary, then it would be positive throughout the entire region. To prove the positivity of η(x, z), let us assume the converse, namely that η(x,0) is not positive for all x < 0, and show that this assumption leads to a contradiction. In other words, let us assume that there exists a point Ρ on the lower boundary located at xp < 0 where the vorticity vanishes:

. (4.27)

We shall restrict our attention to the case where Ρ is the only point on the no-slip boundary where η vanishes. In this case, η(x,0), which tends to 1 as x→–∞, decreases to 0 as x reaches P, and finally tends to –∞, as x approaches the origin as can be seen from approximation (4.21).

A curve, on which η(x, z) = 0, emanates out of the point P. This curve can end (i) upstream asymptotically close to the free upper surface or (ii) at a point Q, on the upper boundary or (iii) downstream. The argument which we are about to give is independent of the actual shape of the η(x, z) = 0 line. Therefore, we pick one such shape, say (ii), and proceed. Fig. 2 illustrates this curve as well as generic lines of constant vorticity. Consider two such nearby lines of constant vorticity in the upstream half of the flow region, say AB and A’B’. Clearly, since ζ is harmonic, , and consequently we can write

.

or by Green’s theorem

In view of the orthogonality of the constant η and ζ lines, the last integral is identically zero and this equation reduces to

.

Finally, since by the Cauchy-Riemann conditions the normal derivative of ζ vanishes on AB and Β’Ά’, we conclude that

. (4.28)

Fig. 2. Proof that vorticity η is everywhere positive. If we assume the converse, then η = 0 at a pointpon the lower boundary, and this point must lie upstream of the origin. (Downstream, η is zero on the lower boundary by definition.) A contour on which η = 0 emanates from Ρ and intersects the upper surface of the ice at point Q. Two points, Β and B’ shown above, can be chosen so that . Let the contours of constant η which emanate from Β and B’extend far upstream towards the points A and A’, respectively. A contradiction arises when we consider the integral of over the region bounded by these two contours (indicated by shading above). This contradiction negates the existence of point Ρ and thus assures us that η≥0 everywhere.

Since the lines of constant vorticity AB and A’B’ were chosen arbitrarily, we are forced to conclude that

(4.29)

or by the Cauchy-Riemann conditions (4.23)

(4.30)

which is clearly a contradiction. Thus

(4.31)

and the vorticity is positive throughout the region.

Returning to the equation for ψ in Equations (4.3), we can immediately infer from the fact that that ψ cannot have a local maximum in the region. Therefore

. (4.32)

In fact, we can show that ψ lies between-⅓ and 0. To that end we assume the converse, namely that there exists a region R where

. (4.33).

In this region, there exists a stream line Ψ=Ψ where Ψ is smaller than-⅓. Since the value of the stream function on the boundaries of the domain lies above-⅓, the stream line Ψ=Ψ must be closed. One such hypothetical region R enclosed by this stream line is illustrated in Figure 3. Clearly

. (4.34)

Finally, let G be the Green function for Laplace’s equation in this region, namely

. (4.35)

It is easy to show that G is positive throughout R . Furthermore, through the usual integral formulas, Ψ in R is given by

. (4.36)

This implies that ψ is larger than Ψ, which contradicts condition (4.34). Thus, not only have we shown that

(4.37)

but also that no recirculation is possible.

Fig. 3. Stream-line geometry if recirculation were possible. The positivity of η prevents ψ from having a local extremum, thus ruling out recirculation.

In summary, many of the essential properties of the flow across the transition in basal boundary condition have been derived without resort to a solution, either exact or approximate, of the governing equations. Our expectation of finite viscous energy dissipation has guided us to the realization that stress concentrations can occur only at one point, the point of transition at the bed. This result is reassuring for the fact that slow, creeping motions of viscous fluids are not expected to give rise to strong flow gradients within the interior of the fluid. The analysis of the vorticity field, which obeys Laplace’s equation, yields the simple result that the vorticity field is everywhere positive. With this result, we have ruled out the possibility of small eddies of recirculating flow. This conclusion suggests that the no-slip to free-slip transition at the bed cannot, by its own volition, give rise to an overturned stratigraphic sequence in ice cores collected downstream (e.g. Weertman, 1976).

5. The Underlying Wiener-Hopf Equation

Except for considerations pertaining to the free upper surface, our first-order problem is identical to that solved by Richardson (1970) and by Hutter and Olunloyo (1980). We shall therefore abbreviate the steps required to arrive at the solution, and give only what is necessary to present a self-contained analysis.

We focus on the stream function exclusively, and write the boundary-value problem as

(5.1)

with

(5.2)

The fact that Ψ(x,0) does not tend to zero as precludes us from immediately using Fourier transforms in the x-direction. However, this obstacle can be circumvented by replacing the constant ⅓ in Equations (5.2) by , solving the problem thus modified by means of Fourier transform, and then take the limit ∊→0 (see e.g. Richardson, 1970). We shall use this approach and replace the boundary conditions (5.2) by

(5.3)

Denoting Fourier transforms by a caret, we define

. (5.4)

The transformed stream function is a solution of the following problem:

(5.5)

together with

(5.6)

where primes denote derivatives with respect to z.

The subscripts + and — not only remind us that the range of integration is the positive or negative interval, but much more importantly, that the functions and

are analytic functions of ξ in the upper and lower halves of the complex ξ-plane (Noble, 1959, p. 12), respectively.

The most general solution of the transformed equation which satisfies the boundary conditions at z = 1 is

(5.7)

The boundary conditions at z = 0 imply that

(5.8)

Solving the last two equations for A and B, we see that

(5.9)

Substituting these expressions in the first equation in Equation (5.8), we deduce the fundamental Wiener-Hopf equation:

(5.10)

As is usual with Wiener-Hopf equations (Noble, 1959), Equation (5.10) is a single equation for two “half” unknowns, viz. and However, because these two unknowns have different domains of analyticity in the complex plane, we can split this single equation into two separate equations for o and alone. The Wiener-Hopf technique can be recognized as being conceptually similar to the separation of variables method for solving separable partial differential equations in the following sense. A single separable partial differential equation can be converted to several ordinary differential equations by writing the solution of the former equation as the product of functions of a single variable. Likewise, the Wiener-Hopf equation posed above will be converted to two equations by writing its solution as the product of two functions which have different properties in the complex ξ-plane.

Before tackling the solution of the Wiener-Hopf equation we ought to note a feature of A and B, namely that they have no poles at the zeros {ξ _n } of

(5.11)

Indeed, Equations (5.10) implies that

(5.12)

which after replacing ξ _n by , and dividing by sinh ξ _n implies that

. (5.13)

The above expression is proportional to the numerator of Β as given in Equations (5.9); thus, Β has no poles at ξ _n.

Similarly, because of Equation (5.13)

(5.14)

This shows that A has no poles at ξ _n . Our purpose in making the above observation is two-fold. First, in so doing we avoid the future computation of unnecessary residues. Secondly, we suspect that this is the source of the error in Hutter and Olunloyo’s (1980) calculation. Indeed, although they started out from the same Wiener-Hopf Equation (5.10) as we do, their approximate treatment relies on the replacement of the coefficients 2ξsinh ²ξ and sinhξcosξ–ξ by simpler functions. More specifically, they were able to achieve an approximate solution of Equation (5.10)) with much less effort by replacing the function to be split, namely

(5.15)

by the much simpler one

(5.16)

This step changes the numerators and denominators of A and Β in such a way that they no longer vanish simultaneously at the ξ _n s. Thus, in Hutter and Olunloyo’s treatment of the problem, the ξ _n are poles which must be taken into account, even though their paper made no mention of them.

6. Solution of the Wiener-Hopf Equation

In preparation for the solution of the Wiener-Hopf Equation (5.10), we examine each of the functions entering in this equation. We start with the coefficient of , which is obviously an entire function of ξ. It is well known that the hyperbolic sine of ξ can be expressed as an infinite product. In fact

(6.1)

This is just a special case of the general Hadamard’s factorization theorem (Titchmarsh, 1950, p. 250). Since we want to separate the zeros which lie in the upper half-plane from those lying in the lower half-plane, the convergence factors are needed in the representation. Note also that behavior at the origin determines the factor in front of the products. Consider next the coefficient of , viz.

. (6.2)

By writing f as follows

we see that f is an entire function of order 5 (see Titchmarsh (1950, p. 248) for a definition of “order”) of the variable . Therefore, it has an infinite number of zeros. Except for a triple zero at the origin, these zeros are complex. We shall denote those zeros which lie in the first quadrant by {ξ _n } The other zeros are then where * denotes complex conjugation.

For large n, the nth root has the simple asymptotic form

. (6.3)

This asymptotic form will prove useful in evaluating the roots of f(ξ) using iterative techniques which require an initial guess.

Making use of the previous remark about f being an entire function of order ½ of ξ², as well as the location and form of the zeros found above, we deduce from Hadamard’s factorization theorem that

(6.4)

Substituting the product representations in Equations (6.1) and (6.4) of the coefficients in the Wiener-Hopf Equation (5.10), we deduce that

(6.5)

In deriving Equation (6.5), we have succeeded in splitting the lefthand side of Equation (5.10) into two terms which are analytic and zero-free in the upper and lower halves of the ξ-plane, respectively. Incidentally, we now see that the upper half of the ξ-plane consists of the region , where stands for “imaginary part”; similarly, the lower half-plane corresponds to Thus, the two half-planes overlap over the strip

The next step in obtaining two equations out of the single Wiener-Hopf Equation (5.10) is to split the function on the righthand side in Equation (6.5). Let us denote this function by G(ξ;∊)

. (6.6)

G is a meromorphic function of ξ, i.e. a function whose only singularities are poles. We can reach our goal by writing it as a sum of partial fractions (see e.g. Titchmarsh, 1950, p. 110; Noble, 1958, p. 39). The formulas given in these references assume that the meromorphic function is bounded at infinity and has only simple poles. G is bounded at infinity but it has both single and double poles. The partial fraction formula given in the above references must be generalized to account for the double poles in the following manner

(6.7)

where a ^± are the residues at and β the residues at and –ξ _n respectively, and G_n ^–2),G_n ^–1) are the first two coefficients of the Laurent expansion of G around the double poles ξ=inπ At this stage, a very welcome miracle happens: with the exception of a ^± all the other coefficients entering in the partial fraction expansion are of order e. Therefore, they vanish in the limit ∊→0. As a result, taking that limit now, we get

(6.8)

where the subscripts ± remind us that the poles at the origin are the limits of the poles at ±i∊ As a result, Equation (6.5) becomes

(6.9)

where e(ξ) is an unknown function which is analytic in the whole ξ-plane, hence an entire function. As is usually the case with the solution of Wiener-Hopf equations, this entire function is determined by making some assumptions about the singularity in the physical domain at x = z = 0, i.e. at the location where the boundary conditions change character. From our discussion of the singularity at x = z = 0 in section 4, we know that

(6.10)

This implies (Noble, 1959, p. 36) that

(6.11)

This result, together with the assertion that the products in Equation (6.9) tend to one as ξ becomes large, imply that the entire function e(ξ) is identically zero. Consequently,

(6.12)

Equations (6.12) represent the solution to the Wiener-Hopf equation. These expressions may now be substituted into Equations (5.9) and (5.7) to obtain the Fourier transform of ψ(x,z) We remark that the above, rather lengthy analysis has provided an exact solution to the biharmonic equation for ψ.

7. Deformation of the Free Surface and the Pressure Field

Since the shear stress and x-component of the velocity on the bottom have been given by both Richardson (1970) and Hutter and Olunloyo (1980), we proceed directly to the evaluation of the free-surface deformation, which is given by Equation (3.16). To accomplish this evaluation, we require an explicit knowledge of the pressure. We recall our discussion of the pressure in section 4 and begin by evaluating ζ or rather its Fourier transform, from Equation (4.22). Clearly

. (7.1)

We next express in terms of as follows

(7.2)

Making use of the representation (5.7) for ψ, we deduce that

. (7.3)

As a result

. (7.4)

By inverting the above expression and substituting it into Equation (4.24) we get

(7.5)

In particular

(7.6)

where C is a constant picked to satisfy the boundary condition (4.4) upstream. It is easier to determine this constant from considerations regarding the upper surface. We therefore pursue the calculation of the upper surface first and defer that of the pressure field until later.

The next step is the evaluation of w_z at z = 1. This calculation is straightforward

(7.7)

Substituting Equations (7.6) and (7.7) in Equation (3.16), we see that

. (7.8)

Now, according to Equations (5.7) and (6.12)

(7.9)

As a result

(7.10)

where

. (7.11)

In evaluating the above integrals, we have, of course, to find the singularities of the integrands. In addition to the poles at nπi for the integrand of h₁ and , for that of h₂ , these two integrands do have poles at {ξ _n }, which are the zeros of

(7.12)

However, as our discussion at the end of section 5 shows, has no poles at {ξ _n } other words, the residues of and at {ξ _n } cancel each other. Consequently, we need not evaluate them.

Finally, we should note that the integrand of h₂ (1) has a pole at ξ = 0. Making use of the calculus of residues and closing the contours appropriately, we find that

. (7.13)

In the above formulas, H₁[x;] and H₂[x;] stand for the residues of h₁ and h₂ , respectively, at the poles indicated in the square brackets. The residues needed to evaluate h for x > 0 are

. (7.14)

Similarly, the residues for x < 0 are

(7.15)

To satisfy h = 0 as x→–∞,, we require

. (7.16)

To evaluate Equation (7.13) numerically, we had to determine the location of the poles ξ_k on the complex ξplane. This was accomplished by using a software package, the FORTRAN sub-routine ZANLY from the IMSL math library (IMSL, Inc., 1990), that is designed to find zeros of user-supplied functions. The initial guess of the ξ_k supplied to this sub-routine was the asymptotic expression given by Equation (6.3). One hundred poles were found using 10⁻¹⁸ as both the absolute stopping criterion and relative-error stopping criterion (these criteria are defined in the software documentation). The maximum number of iterations allowed in finding any one pole to the accuracy specified in the absolute and relative-error stopping criteria was 10 000. The first 31 of these poles are shown in Figure 4. As described below, the poles ξ_k, k > 100, were not required in the computation of the various fields due to quick convergence of the series and products we evaluated.

Fig. 4. First 31 zeros of Equation (6.2) which lie in the upper right quadrant of the complex ξ-plane.

The series and products expressed by Equations (7.13)–(7.16) were evaluated using a FORTRAN program running on a CRAY-YMP computer which performs floating-point arithmetic with 64-bit accuracy. The stopping criterion for evaluating the infinite series and products was 10⁻⁵. (If the next term in a given series, or the difference between the next term in the product and 1, is less than this criterion, the series or product is truncated at the current level.) Figure 5 displays h(x), which, as a reminder, is the first-order deviation of the free surface from its undisturbed state. As is clearly seen from the figure and by inspection of Equation (7.13), the free surface dips by about 20% over the transition in basal boundary condition, then rises linearly with x at points downstream. This linear rise is necessary to produce linear growth in the hydrostatic pressure which compensates for the gravitational force. The end result is that the free surface becomes horizontal downstream where there is no basal traction.

Fig. 5. The first-order terms of (a) the basal shear stress τ(x), (b) the basal pressure p(x,0) and (c) the free-surface elevation h(x). (a). Open circles on the graph of τ indicate the solution obtained by Hutter and Olunloyo (1980). The concentration of basal shear stress at the point of transition in basal boundary condition was cited by Hutter and Olunloyo (1981) as a possible cause of high till concentrations at the beds of some glaciers, (b). The angularity of the basal pressure represents a fundamental difficulty with the solution of the problem. We suggest that thermodynamic processes or cavitation must be considered to resolve this difficulty, (c). The linear increase of h(x) as x→∞, cancels the zero-order surface incline and renders the total free surface horizontal. This is consistent with the lack of basal traction at points x > 0.

For reference, we also show in Figure 5 the pressure p(x,0), to be determined below, and shear stress τ = u_z(x,0) at the base. Of particular interest is the correspondence between the concentration of basal shear stress and an increase in surface slope associated with the upstream limb of the free-surface dip over the origin. While the basal shear stress is singular at the origin, the free-surface elevation and slope are both finite and continuous.

We return to the pressure field as given in Equation (7.5) with C given by Equation (7.16). Following the same steps as for the surface evaluation, we write

(7.17)

where

(7.18)

Once again, the calculus of residues yields

(7.19)

where

(7.20)

and

(7.21)

As shown in Figure 5, the basal pressure (calculated following the same arithmetic procedures as for h(x)) exhibits a modest rise as the origin is approached from upstream. Directly downstream of the origin, p(x, 0) is singular. This indicates that the Taylor-series expression for p(x,z;α) given in Equation (3.1) is disordered at x = z = 0. Such disordering presents a fundamental inconsistency with the assumptions required to solve the problem.* (Ideally, we would like a solution in which the deviations of the pressure are uniformly bounded.) We shall discuss the consequences of this disorder below. Also of note in Figure 5 is the linear increase of p(x, 0) as x→+∞,. This represents the increasing hydrostatic pressure at the bed associated with the horizontal free surface.

A contour map of the first-order correction to the hydrostatic pressure p(x, z) is shown in Figure 6. To produce this contour map, we evaluated p(x,z) on a 31 by 31 grid covering the region — ½ < x < ½ and 0 < z < 1. The series representation of p for x > 0, Equation (7.19), involves summing which contains a cos(nπz) factor. Gibbs phenomena associated with this factor made the resulting contours jittery near x = 0. To eliminate this unrealistic noise, we introduced the Lanczos convergence factor, sin(nπ/m)/(πn/m), as an extra factor in the definition of . Here, n = m defines where the infinite series for p given by Equation (7.19) is truncated. The Lanczos convergence factor is a well-known remedy for Gibbs phenomena associated with truncated sine and cosine series (e.g. Arfken, 1970, p.659).

Fig. 6. Pressure near the transition from no-slip to free-slip. The strong drop in pressure downstream of the transition is of particular interest and suggests that thermodynamic or mechanical effects,which are not considered in our study, may come to bear in more realistic situations.

As with the vorticity field to be discussed below, the contours of constant pressure fan out of the origin and intersect the bed downstream of x = 0. The region where p < 0 is confined to a narrow band located just downstream of the basal boundary transition. At points far upstream, p(x, z) tends to zero as it should. Downstream, contours of constant pressure become vertical and are spaced at constant intervals in x. This reflects the fact that h(x) increases linearly with x.

As noted above, a region exists near the origin where p becomes strongly negative. Within this region, the assumptions used to justify our approximations break down. Foremost of these assumptions is the fact that we require the Taylor-series expansion (3.1) for p in terms of the small slope α be ordered as x,z→0. This series is clearly disordered when p→–∞, because p, the first-order pressure, will eventually exceed

, the zeroth-order pressure.

Another concern is that thermodynamic effects, which so far have been assumed to be unimportant, may come into play where there is variable pressure along a presumably thawed glacial bed. A strong pressure drop would cause the freezing point to elevate close to the transition in a basal boundary condition. With this elevation, a heat flux would be required to avoid basal freezing and consequent downstream migration of the basal transition (which would violate our assumption of steady state). We do not resolve this inconsistency here. However, we suspect that thermodynamic or mechanical effects associated with strong pressure drops may play a role in this region. Such effects could include: 1. cavitation (e.g. Lliboutry, 1968; Nye, 1970), 2. wake separation (e.g. Clarke and Blake, 1991), and 3. flow associated with the non-Newtonian properties of ice.

Clarke and Blake (1991) offered intriguing suggestions as to the possible consequences of an abrupt transition in the basal boundary condition which our finding concerning the pressure drop supports. Their analysis of Trapridge Glacier suggests that, in the case of a free-slip to no-slip transition (the reverse of our problem), an edge dislocation may migrate up from the bed into the interior of the ice. The flow geometry and basal conditions on Trapridge Glacier are far more complicated than what we allow for here; thus, further comment on Clarke and Blake’s (1991) suggestions is clearly beyond the scope of this study.

8. The Velocity Fields

We turn now to the tedious calculation of the velocity fields. We return to the formula (5.7) for and replace the constants A and Β by their expressions given in Equation (7.13). The result is

. (8.1)

Therefore

(8.2)

where

(8.3a)

and

(8.3b)

The integrands in both of these integrals have simple poles at ξ = 0. The subscripts + and — are a reminder that the contour must be indented above the origin in the first integral and below in the second. The contribution from these poles will yield the Poiseuille flow observed as x→–∞, and the plug flow observed as x→+∞,

The analysis proceeds along the same line as before, namely the evaluation via calculus of residues of these integrals with different expressions for x > 0 and x < 0.

(8.4)

As previously explained, U₁[x,z;] and U₂[x,z;] stand for the residues of u₁ , u₂ at the poles indicated in the square brackets. Note the Poiseuille and plug flows arising from the contributions of the poles at the origin. The other residues needed to evaluate u for x > 0 are

(8.5)

The residues needed to evaluate u(x, z) for x < 0 are

(8.6)

We record the expression for the vertical velocity even though we shall not compute it by the calculus of residues

(8.7)

where

(8.8)

The above integrands have no singularities at ξ = 0 and hence the contours need not be indented. Because of the lack of singularities at the origin, the vertical velocity tends to zero far away, as it should.

Vertical and longitudinal profiles of u(x, z) are shown in Figures 7 and 8, respectively. The vertical profiles display the transition between the Poiseuille flow and the plug flow. The longitudinal profiles display the smooth transition between no-slip and free-slip at the bed, as well as the smooth change on planes parallel to the bed. The reduction of surface velocity across the origin is smooth, unlike the result displayed by Hutter and Olunloyo (1980). This suggests that transitions in the basal boundary condition may not so easily be identified from surface-velocity measurements as was once suggested.

Fig. 7. Streamwise velocity (solid lines) as a function of z. Velocity u at various positions x_i,(x₁ =-1.0, x₂ = −0.1, x₃ =-0.1, x₄ = 0.5 and x₅ = 1.0 along the stream. At x = 1.0, the velocity profile resembles the Poiseuille flow imposed as x→–∞,. At x = 1.0, the profile is close to the plug flow expected as x→+∞,. The transition from Poiseuille to plug flow occurs over a length comparable to the depth of the stream.

Fig. 8. Streamwise velocity as a function of x. Velocity u at various for depths z_i (z₁ = 1.0, z₂= 0.5 and z₃ = 0). Open circles, crosses and open triangles denote the solutions obtained by Hutter and Olunloyo (1980). (These points were digitized from their figure 4.) The narrow peak in Hutter and Olunloyo’s solutions for z = 1.0 and z = 0.5 near x = 0 stems from an arithmetic error in their calculations. (Again, we thank Professor Hutter for bringing this error to our attention.) According to our results (solid lines), the transition from Poiseuille flow to plug flow is monotonie at all depths.

9. The Vorticity Field

Starting from Equation (7.8) and following the same steps as followed previously, we write

(9.1)

where

(9.2)

Note that the integrand of η₂ has a simple pole at ξ = 0 but not that of η₁.

The calculus of residues yields

(9.3)

The residues required to evaluate η for x> 0 are

. (9.5)

Residues needed for x < 0 are

(9.6)

A contour map of η(x, z) is displayed in Figure 9. As with the pressure map (Fig. 6), Lanczos convergence factors were used to eliminate Gibbs phenomena in the series evaluations for x > 0. Notable features seen in Figure 9 include the fanning of the contours out of the origin and the fact that η is positive throughout the entire domain. These features are confirmed by the a priori results derived in section 4. The singularity of η at the origin corresponds to the singularities of the basal shear stress, τ, and the basal pressure, p(x, 0).

Fig. 9. Vorticity near the transition from no-slip to free-slip. Note the fanning of contours which emanate from the origin (x = z = 0). At the origin, η is undefined. The zero-level contour runs along the free surface (z = 1) and along the x > 0 part of the basal surface. The η = 1 level contour emanates from the origin and asymptotically approaches the basal boundary as x→–∞,. In the vicinity of the origin, contours with value greater than 1 emanate from the origin and fold over to intersect the bed upstream of x = 0. The geometry of these contours was anticipated by the preliminary analysis of section 4.

10. The Far Fields

We have seen that as x→∞, the upper surface has the following behavior

. (10.1)

Thus for large distances, and more specifically for , the series expansion in powers of α becomes disordered. Contrary to the situation for the pressure near the origin, this second disorder manifestation can be dealt with. This requires us to investigate the flow fields at such large distances, and to that effect we introduce a downstream far-field coordinate (Hinch, 1991, p. 52)

(10.2)

as well as far-field variables

(10.3a)

Substituting these expressions in the dimensionless version of the governing Equations (2.1), we see that

(10.3b)

The boundary conditions on the bottom are those for the stress-free part, viz.

(10.4)

whereas the dynamical boundary conditions at the top are the dimensionless version of Equations (2.4)

(10.5)

Finally, the kinematic boundary condition reads

. (10.6)

We once again look for series solutions in α of the following form

(10.7)

It is easy to deduce that the zero-order pressure is hydrostatic

. (10.8)

The first-order equations yield

(10.9)

with

(10.10)

and

. (10.11)

Integrating the x-component of the momentum equation in Equations (10.9) over z, mindful of Equation (10.11), we deduce that

(10.12)

or in view of Equation

. (10.13)

In order to satisfy the dynamical condition at the upper surface, must vanish, i.e.

(10.14)

where c is an unknown constant. However, it is clear that if this far field is to match the near field as given by Equation (10.1), then c = 1 and so

. (10.15)

In summary, after considering the first-order correction to the pressure, we have found that

(10.16)

The second-order continuity equation implies that

and as a result

. (10.17)

Substituting this expression in the kinematic boundary condition, we deduce that

which is easily recognized as the integrated version of the continuity equation

. (10.18)

If these newly derived far fields are to match with the previously derived results, now looked upon as near fields, then certainly the constant in Equation (10.18) must be taken equal to the total (dimensionless) mass flux, namely ⅓. Therefore, to leading order, the far fields are

(10.19)

The far-field representations just obtained match the near-field solutions as X→0. Furthermore, they remain ordered for large X _s Therefore, they show that the previously found near-field solutions are valid as long as the downstream distance x times the slope α is not large.

11. Conclusion

We have undertaken a somewhat tedious re-evaluation of a classic theoretical problem in glaciology: the flow of a viscous ice stream across a no-slip to free-slip transition in basal boundary condition. The results are: a confirmation or correction of previous results, new results pertaining to the vorticity and the free-surface expression of the basal transition, and, most importantly, the discovery of a fundamental inconsistency with the solution which renders the problem as originally posed still open to further analysis. We have learned that the surface expression of the transition in basal boundary condition is not as marked as implied by Hutter and Olunloyo (1980). The surface velocity field changes gradually over distances of order of several ice thicknesses, from the Poiseuille flow upstream to the plug flow downstream. Thus, it may not be possible to diagnose sharp changes in basal conditions by measurement of the streamwise surface velocity. This suspicion is consistent with previous studies of the transfer of basal-sliding variations to the surface. Balise and Raymond (1985), for example, argued that sliding variations are expressed by a diffuse variation in surface velocity.

The free-surface elevation, which we derive here for the first time, provides the most clear evidence of the basal transition visible in surface measurements. The free surface dips by about 20% over the basal transition, then becomes horizontal downstream. This result is consistent with the observed topography of various ice-sheet to ice-shelf transitions. (We remark, however, that many processes not accounted for here, such as tidal flexure and transient flow effects, also affect these observations.)

The one result which we find troubling is the fact that the pressure becomes negatively infinite immediately downstream of the basal transition. This feature represents an inconsistency with the assumptions that justify our solution. Having sought a solution to this inconsistency (e.g. by attempting to find a near-field solution in which inertia or some other viscous-stress terms are able to compensate for the infinite pressure) without success, we suggest that the problem, as formulated here, remains open and unsolved. On physical considerations, the drop in basal pressure can signal the need to consider special processes associated with ice, specifically its thermodynamics and rheological behavior. Our suspicion that these special processes come into play is supported by the observations of Trapridge Glacier (Clarke and Blake, 1991) which suggest the migration of an edge dislocation into the ice from the point of basal transition.

As a final remark, we offer the results of this analysis as a means to confirm the performance of various numerical ice-flow models and inverse methods. If there is a value in the solution of highly idealized flow problems, it is in the opportunities these solutions offer to test the numerical models created to embrace more realistic flow problems.