List of Symbols
Subscripts x , y , X and Y denote partial derivatives, while overbars denote complex conjugation.
1. Introduction
The objective of this paper is to investigate some fundamental aspects of the behaviour of ice streams underlain by a perfectly plastic bed. While tLankenshiphe central parts of ice streams are often dominated by lateral shearing and can be described by relatively simple models (Reference RaymondRaymond, 1996, 2000; Reference Whillans and van derVeenWhillans and Van der Veen, 1997; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000b), the behaviour of the side margins, which can provide a significant proportion of the total drag acting on the ice stream, remains unclear (Reference Jackson and KambJackson and Kamb, 1997; Reference Harrison, Echelmeyer and LarsenHarrison and others, 1998; Reference Jacobson and RaymondJacobson and Raymond, 1998; Reference Raymond, Echelmeyer, Whillans, Doake, Alley and Bindscha-dlerRaymond and others, 2001). Understanding the interactions between bed properties, stress field and temperature in the margins of ice streams is the focus of this paper.
Reference RaymondRaymond (1996) and Reference Jacobson and RaymondJacobson and Raymond (1998) modelled the transition between the rapid sliding of an ice stream and the slow shearing flow typical of inter-stream ridges as a jump in the friction coefficient τ o in the sliding law
where u b is sliding velocity and τ b is basal shear stress. Sliding laws of this type are widely used to describe hard-bed sliding. However, recent studies of till mechanics have indicated that till can be idealized as a Coulomb-plastic material with a yield stress which is determined by the water content of the till, and is independent of strain rate (Reference Iverson, Hooyer and BakerIverson and others, 1998; Reference Tulaczyk, Mickelson and AttigTulaczyk, 1999). If sliding therefore occurs at a yield stress which is predetermined rather than dependent on sliding velocity, it is appropriate to equate basal shear stress to the yield stress in those parts of the bed where failure occurs, rather than to use a relation of the form of Equation (1).
The model used here to study the effect of plastic boundary conditions on ice-stream margins is extremely simple. We consider only variations in the stress field along a crosssection at right angles to the direction of flow of a model ice stream. Ice is assumed to have a constant viscosity, and longitudinal stresses as well as variations in ice thickness and bed topography are ignored. The ice stream is assumed to be underlain everywhere by a plastic bed. Sliding occurs where the basal shear stress equals the yield stress of the bed (which is allowed to depend on position), while no sliding is assumed to occur where the yield stress of the bed is not attained. The same sliding behaviour would also result from taking m → ∞ in the hard-bed sliding law (Equation (1)). Processes such as regelation, which are likely to give rise to finite sliding velocities even when the yield stress of basal till is not attained, are ignored by our modelling assumptions. Along with the assumption of a constant ice viscosity, this is likely to be one of the main limitations of our model.
The simplifications outlined above allow us to make considerable progress analytically, and to understand aspects of the stress field in the margins which are difficult to resolve otherwise. In particular, we are able to demonstrate that certain modelling assumptions naturally predict singularities in the stress field, and to show how such singularities can be avoided. Similar approaches were taken in the papers of Reference Hutter and OlunloyoHutter and Olunloyo (1980) and Reference Barcilon and MacAyealBarcilon and MacAyeal (1993), which concerned the problem of transitions between no-slip and free-slip regions along a flowline. Although a physically more realistic model is desirable, the added complications introduced by, for instance, a temperature-dependent Glen’s law rheology (Reference PatersonPaterson, 1994) render the kind of analysis presented here much more difficult, and therefore actually obscure our understanding of some of the fundamental aspects of the behaviour of shear margins.
One of the main reasons for studying shear margins is to improve our understanding of the processes which cause them to migrate, as changes in the width of an ice stream can have a significant effect on its discharge (Reference Van der Veen and WhillansVan der Veen and Whillans, 1996). InnReference Jacobson and RaymondJacobson and Raymond’s (1998) approach to this problem, it is unclear why there should be a discontinuity in the sliding coefficient To, and how the position of the discontinuity migrates if it exists. By contrast, the evolution of the yield stress of a plastic bed can be understood straightforwardly in terms of changes in the water content of the bed. As the bed absorbs or releases water, its porosity changes, which implies a change in effective pressure and therefore in yield stress (Reference Tulaczyk, Mickelson and AttigTulaczyk, 1999; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000a). In turn, changes in the water content of the bed are controlled by drainage and by basal melting or freezing. Although the poorly understood problem of drain- age at the base of soft-bedded ice streams is beyond the scope of the present work (cf. Reference NgNg, 1998), we will investigate how strain heating in the ice might result in melting and freezing patterns at the bed.
2. Model
2.1. The mechanical problem
The geometry of the problem is illustrated in Figure 1. We consider an infinite parallel-sided slab of ice with constant viscosity η and density ρ, inclined at some angle α > 0 to the horizontal. The ice thickness will be denoted by H. At the bed, a yield stress τc is prescribed which depends only on the coordinate x measuring distance in the cross-slope direction. Sliding occurs where basal shear stress attains this yield stress, while no slip is assumed where basal shear stress is less than the yield stress. The ice surface is assumed to be stress-free. It follows straightforwardly from these assumptions that the velocity components in the cross-slope direction and in the direction normal to the slope vanish. Furthermore, the pressure field is hydrostatic, and the only non-zero components of the deviatoric stress tensor are those describing lateral shear and shear in the downstream direction (ηux and ηuy in the notation below).
Geometry of the ice-flow problem: three-dimensional view (top) and downstream view of a cross-section (bottom).
Denoting the downslope component of velocity by u(x, y), where the coordinate y is measured perpendicularly to the bed, the only non-trivial component of the Stokes flow equations is
where g is acceleration due to gravity, and subscripts x and y represent partial derivatives; thus uxx = ∂ 2 u/∂x 2 and similarly for u yy.
Boundary conditions on u(x, y) can be stated as follows: Vanishing shear stress at the ice surface implies
At the bed, regions where sliding occurs are denoted by Γ, while regions where the yield stress is not attained are denoted by 
. To clarify our terminology, a point (x, 0) on the bed will be referred to as inside the ice stream if x ∊ Γ and as outside the ice stream if 
. The appropriate boundary conditions at the bed are then
In addition, we require u(x, y) to remain bounded as |x|→∞.
The mixed boundary conditions at the bed (Equation (4)) need to be supplemented with the inequalities
which ensure that the yield stress is not exceeded, and that sliding can only occur in the downslope direction. The relevance of these physical conditions with respect to finding a solution is that they determine where the areas of till failure Γ lie for a given yield stress distribution τc(x). The margin positions are therefore determined at all times by the spatial distribution of yield stresses, which is borne out by the analysis which follows in section 3.1. Importantly, this implies that the positions of ice-stream margins (which we identify with the boundary points of Γ) depend only on basal mechanical conditions and need not coincide with any thermal boundaries at the bed, such as the transition between an actively freezing and an actively melting bed. Of course, thermal conditions will affect the yield-stress distribution, and hence the the position of the margins, over time.
In a physically more realistic model, one would replace the constant viscosity η by an appropriate non-linear rheology (Reference PatersonPaterson, 1994, ch. 5). If processes such as regelation are to be accounted for, which allow finite sliding velocities even when the local yield stress τ c (x) is not attained, then the no-slip boundary condition (4)i needs to be replaced by an appropriate sliding law relating sliding velocity to basal shear stress, while the conditions (4)2 and (5) remain unchanged.
2.2. The thermal problem
Before proceeding to solve the mechanical problem, we outline the thermal problem which will be used to study the effect of different yield-stress distributions on strain heating in the ice, and hence on melting and freezing at the bed.
Following Reference Jacobson and RaymondJacobson and Raymond (1998), we suppose that the temperature field T(x,y) is determined by rapid diffusion of heat, and that there is no temperate ice, but un- like Jacobson and Raymond, we do not attempt to model the effect of advection due to lateral inflow of ice here. Hence the temperature field in the ice satisfies
where k is the thermal conductivity of ice, assumed constant here. The right hand side of this equation represents strain heating in the ice.
As Reference Jacobson and RaymondJacobson and Raymond (1998) point out, the residence time of ice in an ice stream is typically much shorter than the diffusive time-scale for heat transport across the thickness of the ice. The temperature field T should therefore be determined by advection along flowlines rather than by rapid diffusion. However, this observation is based on velocities near the centre of an ice stream, whereas our main interest is in the marginal areas, where ice velocities are much slower and advection is less important. Furthermore, we will find that strain heating in the margins is concentrated near the bed. The relevant length scale for diffusion is therefore less than the ice thickness, implying in turn a shorter time-scale for diffusion than might be expected. Our approach of assuming rapid diffusion may therefore yield reasonable results for heat fluxes near the bed in the marginal areas of the ice stream.
As boundary conditions on the quasi-steady heat diffusion problem (Equation (6)), we prescribe a constant temperature T s at the surface y = H, although this neglects the effect of cool air pooling in surface crevasses (Reference Harrison, Echelmeyer and LarsenHarrison and others, 1998). At the bed, it is reasonable to expect the areas in which till failure occurs to be at the melting point T = T m. As stated above, it is unlikely that the margin positions coincide exactly with thermal boundaries which exist at the bed, such as the transition between a frozen and an unfrozen bed. The bed may be temperate not only under the ice stream, but also in regions outside the margins, although basal water pressures are low there. Hence, even if there is a temperate-bed/frozen-bed transition somewhere outside the margins, it is not clear where its position should be. To avoid difficulties associated with choosing this position arbitrarily, we assume that the bed is temperate everywhere, so T = T m when y = 0. Note that these boundary conditions differ somewhat from those in Reference Jacobson and RaymondJacobson and Raymond (1998); in particular, the choice of constant temperature at the bed decouples the thermal problem in the ice from that in the bedrock, which may be taken to supply a constant geothermal heat flux q geo.
The rate of basal melting m(x) (or freezing if negative) determines, along with drainage, the rate at which the basal yield stress changes locally. m(x) is given by
where L is the latent heat of melting per unit volume of water produced, and the second term on the righthand side represents frictional heat dissipation in the bed.
The geothermal heat flux q geo and the temperature difference T m − T s, which together cause a background rate of heat gain or loss at the bed, can be prescribed independently of the mechanical problem, and are therefore of secondary interest here. The contributions to basal melting due to strain heating in the bed and in the ice consist of the heat dissipation term u(x, 0) τ c (x) and the following anomalous heat flux:
It can be shown that q strain is the heat flux kTy (x, 0) into the bed which would result from setting T m = T s in the boundary conditions for Equation (6). q strain therefore does not depend on the applied temperature difference T m − T s and purely measures the effect of strain heating in the ice, which is one of our central interests in considering the thermal problem. In terms of q strain, the rate of basal melting can be written as
The term in curly brackets on the righthand side denotes the background thermal conditions, while the last two terms, which are always positive, describe the effects of strain heating. We emphasize that the frictional heat dissipation term u(x, 0) τ c (x) vanishes where till failure does not occur; this is entirely a result of having prescribed zero slip where basal shear stress does not attain the yield stress. In a more realistic model, where sliding outside the ice stream is allowed for, frictional dissipation could be significant and contribute to a weakening of the bed there. In that case, however, the appropriate strain-heating term would be ηu y (x,0) u(x,0) as ηu y (x, 0) = τc(x) only within the ice stream.
In theory, it would be possible to make our ice-stream cross-section evolve in time simply by defining a constitutive relationship between yield stress and water content of the bed, and by ignoring water fluxes due to subglacial drainage (cf. Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000a). This, however, complicates the model and is probably not warranted in view of the simplifying assumptions that have already gone into it. Instead, we only aim to use the thermal model above to explore qualitatively the spatial distribution of heat fluxes generated by strain heating in the ice.
3. Method of Solution
3.1. The stress field
The main complication in the mechanical model of section 2.1 is that the regions in which till failure occurs are not known from the outset, but that their location is fixed implicitly by the conditions of Equation (5) and must therefore be found as part of the solution. Mathematically similar problems occur in the determination of contact areas between elastic bodies (Reference EnglandEngland, 1971), which is why they are also known as contact problems. For a bounded cross section rather than the infinite slab assumed here, the ice- flow problem can be cast as a so-called non-coercive variational inequality, which allows the existence and uniqueness of weak solutions to be deduced under appropriate conditions on the yield-stress distribution τ c (x) (Kinderlehrer and Stampacchia, 1980, ch. 3; see also Reference GillowGillow, 1998). A variational approach has the further advantage of allowing solutions to be constructed without the need to calculate explicitly the contact points at which the switch in boundary conditions occurs. In this paper, we use a conceptually simpler, though less elegant, complex variable method (Reference EnglandEngland, 1971), as this leads to some important qualitative insights into the stress field in the ice-stream margins. We further restrict ourselves to the case where the region of till failure consists of a single interval, Γ = (a, b), where a and b will be referred to as the margin positions (or just as the margins). As mentioned before, the phrase “inside the margins” will refer to x in the interval Γ, while “outside the margins” will refer to x not in Γ. The restriction to a single ice stream is not essential as our method of solution can be extended straightforwardly to the case where Γ consists of a finite number of intervals, but it simplifies the algebra involved. Inevitably, the solution procedure outlined below still involves a fair amount of analysis, which some readers may wish to skip. In order to facilitate this, we have made the presentation of results in section 4 as self-contained as possible.
Before proceeding further, we recast the model of section 2 in terms of convenient dimensionless variables by defining
This change of variables ensures that U satisfies Laplace’s equation in X and Y,
and can consequently be expressed as the real part of an analytic function of the complex variable z = X + iY. Hence
where θ(z) is an analytic function, and an overbar denotes complex conjugation. Note that the first derivative of θ(z), denoted by θ′(z), can be thought of as representing the stress field in the ice,
where 
and 
stand for real and imaginary part, respectively. The boundary conditions (3) and (4), which determine the function θ(z), become (in terms of the scaled variables chosen above)
where we have defined the following scaled difference between yield stress and driving stress for notational simplicity:
In addition, U(X, Y) is required to be bounded as |X| → ∞.
For later convenience, note that Equation (16) implies that
Using the usual differentiation rules (Reference EnglandEngland, 1971, ch. 1), the boundary conditions (14), (15) and (18) become in terms of θ(z), respectively,
where boundary values taken as z approaches the appropriate boundary from within the strip 0 < Y < π are implied. In addition, the requirement that U(X,Y) be bounded at large |X| now becomes
Next, we define the function Ω(z) as
By the Schwarz reflection principle (e.g. Reference MuskhelishviliMuskhelishvili, 1992, p.94−95), Ω(z) is analytic in the strips −π < Y < 0 and 0 < Y < π, while Equation (21) shows that it is also analytic across the real axis with the exception of the interval [a*,b*]. The remaining boundary conditions (19), (20) and (22) can be written in terms of Ω(z) as
where superscripts + and − in Equation (24) indicate limits taken as the branch cut is approached from above and below, respectively.
Equation (24) is suggestive of a classical Hilbert problem (e.g. Reference MuskhelishviliMuskhelishvili, 1992), the main complication being that Ω(z) is defined only on the strip −π < Y < π. To avoid this difficulty, we define the following conformal transformation
which maps the lines Y = π and Y = −π onto the negative and positive halves of the real axis in the ζ plane, respectively, while the branch cut Y = 0, a* ≤ X ≤ b* is mapped onto the part of the imaginary axis in the ζ plane for which exp(a*/2) ≤ = ℘ (ζ) ≤ exp(b*/2). For simplicity of notation, we let ζ = ξ + iv, where ξ and v are real, and put v a = exp(a*/2), V b = exp(b*/2) (Fig. 2).
Geometry of the ζ plane with the contour L.
The boundary value problem for Ω(z) is transformed to the ζ plane by defining the function Φ(ζ) as
which is easily seen to be single-valued and analytic in the upper half of the ζ plane cut along the line segment ξ = 0, v a ≤ v ≤ V b . Boundary conditions for Φ(ζ) follow from Equations (24−26) as
where superscripts + and − now indicate limits taken as the branch cut is approached from the left and right, respectively (Fig. 2), and 
is defined by
As before, limits taken as ζ approaches the real axis from above are implied in Equation (30).
Using the Schwarz reflection principle and Equation (30), Φ(ζ) can be continued analytically across the real axis by defining
such that Φ(ζ) is analytic in the entire ζ plane cut along the two line segments ξ = 0, v a ≤ |v| ≤ V b , where it satisfies
where the superscripts + and — again indicate limits taken as the branch cuts are approached from the left and right, respectively (Fig. 2).
Allowable solutions Φ(ζ) are further constrained by Equations (31), (33) and (23), which yield the conditions
One other constraint on Φ(ζ) arises because the boundary condition (16) was differentiated to give Equation (18) in the derivation of the boundary value problem for Φ(ζ). It therefore remains to ensure that we have not only U x = 0 but also U = 0 on both half-lines X < a* and X > b*. In order to make certain that putting U = 0 on one of these half lines (which is always possible) implies U = 0 on the other, it suffices to require that
For practical purposes (in particular, to avoid later having to evaluate the boundary values of Φ(ζ) numerically), it is convenient to reformulate this constraint as follows: First, notice that Equation (36) implies, by Cauchy’s theorem, that
for any contour L which encloses the branch cut ξ = 0, v a ≤ v ≤ V b , but does not enclose or cross the branch cut which lies in the lower half of the ζ plane. Choosing L to be the union of a line segment – R≤ ξ≤ R on the real axis and a semicircle of radius R in the upper half-plane (Fig. 2), and letting R → ∞, we have by Jordan’s lemma the following condition, equivalent to Equation (36):
Assuming that Δ(v) (and hence τc(x)) is Holder continuous, Equation (34) is a classical Hilbert problem, to which there are a variety of solutions differing in number and position of singularities (Reference MuskhelishviliMuskhelishvili, 1992, ch. 10). The nature of the original contact problem dictates that we should only consider solutions which have no singularities at all; as we shall see below, allowing singular solutions is tantamount to assuming an infinite yield stress just outside the margins. Solutions to Equation (34) which are bounded (nonsingular) in the finite part of the c plane must take the form
Where
the branch taken having cuts along ξ = 0, v a ≤ |v| ≤ v b and behaving as χ(ζ) ~ ζ2 when ζ → ∞. Evaluating χ+ (iv) on the two branch cuts yields
where 
denotes a positive square root. From this follows the equivalent expression
In addition to solving Equation (34), a solution must also satisfy the additional conditions (35) and (36). As the only free parameters in Equation (41) are the positions of the branch points iva and ivb , the additional constraints which arise serve to determine v a and v b, and hence the location of the ice-stream margins a and b. Of course, as there are only two real quantities (va and vb ) to be determined, it is important that the conditions (35) and (36) should not lead to more than two real constraints.
Evidently, Φ(ζ) defined in Equation (41) vanishes at the origin ζ = 0, so Equation (35)
1 is satisfied. Some slightly awkward manipulations, taking care with the location of branch cuts and the branches chosen, reveal that 
and consequently that Equation (35)
3 is always satisfied by Φ(ζ) defined in Equation (41). Hence only Equations (35)
2 and (36) can impose the required constraints on va
and vb
. To deal with Equation (35)
2, note that the behaviour of Φ(ζ) near the point at infinity is
Consequently Φ(∞) = 0 is equivalent to the single real constraint
Lastly, the condition (36) can be written explicitly as
which is again a single real constraint.
Equations (43) and (44) are sufficient in principle to determine the values of v a and v b , and hence the margin positions a and b. Of course, these constraints are not by themselves sufficient for a solution to be physically acceptable, as the conditions (5) must also be satisfied. A particular example of a solution to Equations (43) and (44) which can be unphysical is the choice va = vb , with vb arbitrary; this solution corresponds to the absence of ice streams. The inequality constraints (Equation (5)) must therefore be checked once a solution has been calculated. To make matters more complicated, there are yield-stress distributions for which no solution to Equations (43) and (44) is acceptable, namely those which require more than one ice stream to be present. Suffice it to say that the method above can be extended to accommodate multiple ice streams as well; details are beyond the scope of this paper.
For practical purposes, it is convenient to transform back to the original scaled coordinates X and Y before attempting to solve for the margin positions a and b. This involves a considerable amount of algebra which we omit here. It is furthermore important to recast Equation (44) such that the integral exists even when Equation (43) is not satisfied exactly. This is done by subtracting the integrand in Equation (43), multiplied by an appropriate function of ξ, from the integrand in the repeated integral (44) (note that the order of integration in Equation (44) does not commute); this is legitimate because the integral on the left hand side of Equation (43) vanishes at a solution. Suitable forms of the constraints (43) and (44), which furthermore have the expected symmetry in a* and b*, can on this basis be derived as, respectively,
Where
These constraints are solved for a given yield-stress distribution τ c (x) using a backtracking line-search modification of Newton’s method (Reference Dennis and SchnabelDennis and Schnablel, 1996), where we approximate the Jacobian by finite differences. Evaluation of the integrals on the left hand side of Equation (46) at each iteration step is carried out by first splitting the ranges of integration into the triangles a* < X′ < a* + (b* − X) and a* + (b* − X) < X′ < b* as well as the semi-infinite strips X ′ < a*, a* < X < b* and X′ > b*, a* < X < b*, which locates the singularities at X′ = a* and X ′ = b* and the rapid change in the integrand at X = X′ at the ends of the split ranges of integration. The integrals over the strips X′ < a* and X′ > b* are then transformed to a finite range of integration by using the transformations s = 1/(1 + a* − X′) and t = 1/(1 + X′ − b*), while the singularities in Equations (45) and (46) are dealt with numerically using a change of variable suggested by Reference AtkinsonAtkinson (1989, p.306-307). Integrals are evaluated using one of a variety of Gauss−Legendre quadrature rules, typically involving 64, 96 or 128 nodes; in the case of multiple integrals over rectangles and triangles, a separation of variables is carried out first (Reference Isaacson and KellerIsaacson and Keller, 1996).
The stress field proxy U X (X, Y) − iU Y (X, Y) = 2Φ(ζ) for π/2 < arg (ζ) < π can be written in terms of the real variables X and Y as
where G(X) is defined as before, and the branch of the square root outside the integral is continuous for 0 < Y < π and behaves as −exp[(2X + 2iY − a*− b*)/4]/2 when X→+∞.This expression for Φ(ζ) again preserves the expected symmetry in a* and b*. Furthermore, we have ensured that the stress field remains bounded at large |X| even when Equation (43) is not satisfied exactly by subtracting an appropriate multiple of (43) from (41). This, in turn, means that the calculated stress field converges uniformly to the actual stress field as the approximated values of a * and b* approach the values which solve Equations (43) and (44). Numerically, Equation (48) is evaluated in a similar way to the constraints (45) and (46). Given U x and U y , the velocity field can be evaluated by simple quadrature, using the boundary condition U(X, 0) = 0 on X < a* and on X > b*.
3.2. Free slip: a crack problem
One interesting result which arises from the constraint (43) is that it can only be satisfied if τ c (x) − ρgH sin a changes sign somewhere between a and b, as follows straightforwardly from the fact that the integrand in Equation (43) is of the same sign as τ c (x) − ρgH sin a. In other words, a solution to the contact problem of section 2.1 requires the yield stress τ c to exceed the driving stress somewhere inside the ice stream, where till failure occurs. τ c > ρgH sin α cannot be true if there is free slip (t c ≡ 0 every where at the bed of an ice stream). A stress field without singularities is therefore not compatible with free slip in the ice stream.
Extending the analogy with elasticity theory mentioned at the beginning of section 3.1, free slip inside the ice stream, or simply an arbitrary choice of margin positions which does not satisfy the constraints (43) and (44), corresponds to a crack problem. The Hilbert problem (34) combined with the additional constraints (35) and (36) can still be solved for free slip inside the ice stream (indeed, for any arbitrary choice of , va
and vb
). However, these solutions have singularities of the one-over-square-root type at both branch points V
a and V
b, corresponding to integrable singularities in the stress field at the margins. Specifically, these solutions have unbounded basal shear stresses just outside the margin positions. In terms of the condition (5)i, this requires an infinite yield stress at the margins.
3.3. The thermal problem
The thermal problem (Equation (6)) takes the form of the Poisson equation posed on an infinite strip with Dirichlet boundary data. A straightforward method of solution is to use a Green’s function, derived using Fourier transforms. It is important to ensure that the source term on the righthand side is integrable with respect to x. Owing to residual simple shear u
y ~ ρg sin α(H- y)/η far from the ice stream, the rate of strain heating 
which itself can be calculated using the method of section 3.1, is not integrable. An equivalent problem with an integrable source term is obtained straightforwardly if a reduced temperature Θ is defined through
So that Θ satisfies
where the source term a(x, y) on the righthand side is integrable with respect to x. Boundary conditions are Θ = 0 on y = 0, H. Defining the Fourier transform for a generic function f (x, y) which is integrable with respect to x by
Equation (50) can be turned into the ordinary differential boundary value problem
which can be solved through the use of a one-dimensional Green’s function (e.g. Reference WaltmanWaltman, 1986). Our main interest is the heat flux out of the bed, and hence the quantity kΘ y (x, 0). In terms of Fourier transforms, we have
and hence
The anomalous heat flux q strain follows as
Numerical evaluation of kΘy(x, 0) in Equation (54) is carried out by integrating with respect to s using 64-point Gauss−Legendre quadrature. Fast Fourier transforms are used to calculate the Fourier transforms, thereby approximating the infinite integrals by finite ones and computing them using the composite trapezoidal rule. This yields reasonable results provided the range of integration is chosen large enough and we are only concerned with values of q strain in the central parts of that range.
4. Results
Our analysis has shown that arbitrary assumptions about basal yield stresses under an ice stream, and about the location of the ice-stream margins, generally yield solutions with stress singularities. In particular, transitions from no slip to free slip give rise to shear stress singularities at the bed. Similar results were obtained by Reference Hutter and OlunloyoHutter and Olunloyo (1980) and Reference Barcilon and MacAyealBarcilon and MacAyeal (1993) for transitions between no slip and free slip along a flowline. Singularities are typical of abrupt changes in boundary conditions (Reference EnglandEngland, 1971), and their presence in a solution usually indicates that some physical process has been ignored. Inclusion of such a regularizing process then usually leads to a non-singular solution.
In some cases, a model with solutions containing singularities is acceptable, and no attention needs to be paid to regularizing it. Whether this is true depends on what the model is supposed to achieve. In our case, the integrable shear stress singularities predicted for free-slip/no-slip transitions correspond to integrable singularities in the strain-heating rate, so the total amount of strain heating near the margins is finite. However, the strain-heating singularities do cause singularities in the basal melt rate m(x). This is not acceptable for our purposes of understanding how strain heating in the margins and basal melt are affected by the distribution of yield stresses at the bed.
In section 3.1, we showed that solutions with nonsingular stress fields can be calculated for our model provided τc(x) varies (Hölder) continuously and the margin positions a and b are determined such that two non-trivial constraints(Equations (45) and (46)) are satisfied τ c (x) enters into these constraints through Δ(X)). A possible regularizing process which is ignored by assuming a transition from free slip to no slip is therefore the continuous change in τ c from a low value inside the margins to a high value outside.
Some concrete examples are given in Figure 3. Column 1 shows results calculated for a transition between free slip and no slip, with the corresponding singularity in basal shear stress clearly visible in panel a1. The remaining columns show results for continuous yield-stress distributions; in columns 2 and 3 these are of the form
with n = 5 and n = 50, respectively. For n ≫ 1, the yield- stress distribution (Equation (56)) approximates to free slip inside the ice streams, with a rapid increase in yield stress near x = ±5H.This corresponds to an ice stream of a width of about ten times its thickness, which is narrow compared with most real ice streams; the reason for choosing such a limited width is that numerical results are easier to display graphically. In column 4, we have used
which corresponds to a wider, asymmetrical ice stream with a central sticky spot.
Numerical calculations of stress, velocity and basal heat flux. Each column (labelled 1−4) shows results for a particular basal yield-stress distribution specified in the text. Row a shows τ c (x) as a dashed line, and basal shear stress ηuy(x,0)as a solid line. In the axis label τ0 = ρgH sin α. Row b shows the corresponding velocity profiles. u0 = ρgH2 sin α/η in the axis labels. Row c shows velocity contour maps near the left hand margins of the ice streams (a = −5H; −6.1H; −5.2H and−21.2H in c1−c4, respectively). Low-velocity contours are on the left, and high-velocity contours are on the right. Contour intervals between solid lines are ρgH2 sin α/2η in c1−c3, and twice that in c4. Dotted contours are at intervals of ρgH2 sin α/4η. Row d shows contours of the strain-heating rate η(u x 2 + u y 2). Contour intervals are 5ρ2g2H2 sin2 α/η in d1−d3, and twice that in d4. In all cases, the greatest rate of strain heating is experienced close to the bed inside the margins. Row e shows the anomalous heat flux qstrain as a solid line, and heat dissipation qdissipate(x) = τc(x) u(x,0)as a dashed line. q0 = ρ2g2H3 sin2α/η in the axis label. The margin position a is also shown as a vertical dotted line.
As panels a2−a4 show, basal shear stress inside the margins (i.e. as the margins are approached from inside the ice stream) increases in tandem with the increase in τ c (x) (as required by the boundary condition η u y (x, 0) = τ c (x)), and reaches values considerably in excess of the driving stress. This is in agreement with the analysis of section 3.1, which requires the yield stress to exceed the driving stress somewhere inside the margins, but stands in contrast to the results of Reference RaymondRaymond (1996) and Reference Jacobson and RaymondJacobson and Raymond (1998). These authors, whose assumptions regarding basal sliding are different from ours, find a stress relaxation inside the margins. This difference in results is explored further in section 5.
Corresponding to the high basal shear stresses in the margins, panels d2−d4 show that the rate of strain heating in the ice is greatest on the ice-stream side of the margins. Unsurprisingly, the anomalous heat flux q strain also peaks inside the margins (panels e2~e4). Similarly, the rate of heat dissipation in the bed u(x, 0)τ c (x) reaches a maximum close to the margins, probably due to the high yield stresses τ c (x) experienced there. Comparison of panels e2 and e3 further shows that the anomalous heat flux q strain and heat dissipation u(x, 0)τ c (x) become more sharply peaked as the yield-stress profile becomes steeper (that is, for larger n in Equation (56)); this can probably be attributed to stresses becoming nearly singular as a discontinuity in yield stress, corresponding to n → ∞, is approached. Note that there is no panel e1 because the Fourier transform method used to calculate basal heat flux cannot resolve the singularity in q strain caused by the strain-heating singularity which occurs for a free-slip/no-slip transition.
Lastly, it should be pointed out that the velocity profiles shown in Figure 3 row b agree well with simplified models of ice streams supported by lateral shear (e.g. Reference RaymondRaymond, 2000; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000b): the bulk flow of the ice streams shown in Figure 3 is dominated by lateral shear u x , while vertical shear u y is negligible by comparison, and only becomes significant near the margins. In fact, it is possible to show analytically (though we will not attempt this here for reasons of space) that the lateral shear stress ηu x calculated from Equation (48) agrees asymptotically with the stress field predicted by the simplified, depth-integrated model (e.g. Reference RaymondRaymond, 1996)
provided the location at which the stress field is evaluated lies inside the ice stream and is not close to the margins (x − a » H and b − x » H), and that τ c does not change rapidly near x. A further interesting feature of panel b4 is the concave velocity profile near the sticky spot. If the strength of the sticky spot were increased further, the ice stream would eventually split in two.
5. Discussion
Our model of an ice stream with a plastic bed shows that stress singularities can be avoided by assuming a continuous yield-stress distribution that determines where the transition from no slip to slip with a given basal shear stress should take place. The location of the margins is thus dictated entirely by mechanical conditions at the bed and need not coincide with any given thermal boundaries that may also exist. In particular, the bed just outside the margins, as well as under the ice stream, could be temperate. In that case, the margins do not coincide with a molten-bed/frozen- bed transition.
How the location of the transition between slip and no slip evolves dictates how the ice-stream margins migrate. Margin migration is therefore governed by changes in the yield-stress distribution τ c (x) over time. In theory, one could derive the rate of margin migration by differentiating the constraints (45) and (46) formally with respect to time and by relating ∂τc/∂t to the rate of basal melting or freezing m(x) through an appropriate constitutive relation while neglecting basal drainage (e.g. Reference Tulaczyk, Mickelson and AttigTulaczyk, 1999). However, the resulting expressions are lengthy and are not reported here.
Let us assume that the background rate of heat loss at the bed due to geothermal flux and the temperature difference between bed and ice-stream surface is not large enough to cause freezing everywhere. Our numerical results indicate that water is supplied to the bed and the yield stress lowered fastest on the ice-stream side of the margins, where strain heating is strongest. This suggests a tendency for the yield-stress profile τ c (x) to become steeper over time, as the yield stress is being lowered fastest inside the margins, where it is already lower than outside the margins. In turn, such a steepening in τ c (x) can be shown to slow down margin migration, and is furthermore likely to lead to strain heating and heat flux becoming increasingly concentrated (cf. panels e2 and e3 of Fig. 3).
The anomalous heat flux q strain cannot become ever more peaked as the temperature in the ice cannot exceed the melting point, and we must therefore have T y (x, 0) ≤ 0, or q strain ≤ k(T m − T s)/H. This suggests that a temperate zone may eventually form in the margins, which must then be modelled separately (e.g. Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler, 2001). The effect of such a temperate zone on water supply to the bed and hence on margin migration is unclear.
Other processes could also play an important role in ice-stream margins. For instance, the large gradients in basal yield stress which may exist near the margins could correspond to large differences in pore-water pressure, which in turn may drive an appreciable amount of water drainage towards the interstream ridges. Drainage could therefore play a role in lowering yield stresses outside the ice stream, which would accelerate margin migration in the outward direction. Likewise, allowing for finite sliding velocities outside the ice stream would lead to heat being dissipated at the bed there. As basal shear stresses are likely to be high just outside the margins (Fig. 3, row a), this heat dissipation could be significant even if sliding is slow, and would again contribute to lowering the yield stress outside the margins. Another process which may occur in practice but would be difficult to incorporate in a model is fracturing of basal ice in the margins due to the high shear stresses there. Fracturing would presumably prevent stresses from exceeding some yield stress for the ice, but how this would affect margin migration is unclear. Lastly, the effect of heat advection due to lateral inflow of ice is likely to affect thermal conditions in the margins, notably by shifting the effect of high strain heating towards the inside of the ice stream and suppressing the outward migration of the margins (see Reference Jacobson and RaymondJacobson and Raymond, 1998).
Our investigation has been based on the notion of a plastic bed, and our results differ considerably from those obtained by Reference RaymondRaymond (1996) and Reference Jacobson and RaymondJacobson and Raymond (1998), who used a hard-bed sliding law of the form of Equation (1) with a discontinuity in the sliding coefficient To. The main difference in results consists of the stress relaxation on the ice-stream side of the margins observed by these authors, compared with the stress intensification there predicted by our model. Although direct comparison of the models is made more complicated by the simplifying assumptions used in the present paper — Reference RaymondRaymond (1996) and Reference Jacobson and RaymondJacobson and Raymond (1998) use a more realistic shear-thinning rheology, as opposed to our constant viscosity — it can be shown that some of the qualitative differences between our results are a direct consequence of the difference in basal boundary conditions, and do not depend on the particular choice of rheology. For instance, the basal sliding velocity U b must be continuous in order to avoid non- integrable stress singularities (e.g. Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler, 2001). Hence, if basal sliding is described by a sliding law (Equation (1)), then τb/τ0 must be continuous. If To jumps from a low value inside the ice stream to a high value outside, as was assumed by Reference RaymondRaymond (1996) and Reference Jacobson and RaymondJacobson and Raymond (1998), then basal shear stress τb must do the same. The basal boundary conditions prescribed by these authors therefore require low basal shear stresses on the ice-stream side of the margins (defined as the location of the discontinuity in τ0). This may be contrasted with the high basal shear stresses predicted on the ice-stream side of the margins (defined as the location at which there is a change from yielding to non-yielding bed) by our model.
The definition of margin position is another important issue in comparing our model with those of Reference RaymondRaymond (1996) and Reference Jacobson and RaymondJacobson and Raymond (1998). Further qualitative insight into this problem can be gained by considering the plastic limit m →∞ in the sliding law (Equation (1)). The finite jump in the coefficient τ0 assumed by Reference RaymondRaymond (1996) and Reference Jacobson and RaymondJacobson and Raymond (1998) then corresponds to a finite jump in basal yield stress τ c in our model. Importantly, this type of yield-stress distribution differs from the two cases which we have considered above, namely the infinite jump from free slip to no slip which gives rise to stress field singularities of the one-over-square-root type, and the continuous yield-stress distributions which avoid such singularities. Although the complex variable treatment of section 3.1 implicitly assumes a (Hölder) continuous yield- stress distribution τc(x), it can be amended straightforwardly to allow for a finite number of simple discontinuities in τc(x). Two results emerge from this analysis (which we do not report here for reasons of space): Firstly, the margins, defined as the boundary points of the region of till failure Γ, do not coincide with the discontinuities in τ c (x). This is relevant because the stress intensification observed by Reference RaymondRaymond (1996) and Reference Jacobson and RaymondJacobson and Raymond (1998) occurs at these discontinuities, and would be located inside the margins according to our definition of margin position. The second interesting feature is that a finite jump in yield stress again gives rise to a stress field singularity, specifically a logarithmic singularity in the lateral shear stress ηu x , which would be difficult to detect using a purely numerical approach.
It is conceivable that allowing the sliding coefficient To in the sliding law (Equation (1)) to vary continuously with transverse position in Reference RaymondRaymond’s (1996) and Reference Jacobson and RaymondJacobson and Raymond’s (1998) model, rather than prescribing a discontinuity, would yield results more closely in line with those obtained here. Nevertheless, the disparity between their results and ours suggests that the stress field in the margins may be highly sensitive to basal conditions in a way in which the lateral-shear dominated stress field in the bulk of the ice stream is not. This clearly makes the modelling of shear margins a much more difficult task, and underlines the need for further field-based research into the basal conditions which prevail in the shear margins of ice streams.
The marginal areas of ice streams are notoriously inaccessible (Reference Harrison, Echelmeyer and LarsenHarrison and others, 1998), making drilling work hazardous and difficult. An interesting question is therefore whether measurements of surface velocity, u(x, H) in our notation, might allow bed conditions in the margins to be deduced (see also Reference Whillans and van der VeenWhillans and Van der Veen, 2001). Unfortunately, the answer is almost certainly no. Even if the margin positions a and b were known exactly, inverting velocity measurements u(x, H) to obtain basal shear stress τc(x) from, say, Equation (48) (which retains the scaled variables used in section 3.1) would still require two ill-posed mathematical procedures. Firstly, the velocity measurements u(x, H) would have to be differentiated numerically to yield the surface velocity gradient u x (x, H), and secondly, a Fredholm integral equation of the first kind with continuous kernel would have to be solved to find τ c (x) in terms of u x (x, H). In view of this ill-posedness, and because of the numerous assumptions behind our model, which render it quantitatively inaccurate, such an inversion is unlikely to produce useful results, especially if there are strong spatial variations in τc(x).
6. Conclusions
Our model shows that the location of areas in which basal sliding occurs at a basal yield stress τ c is determined in a non-trivial way by the spatial distribution of the difference between basal yield stress τc(x) and the driving stress ρgH sin α. Simplistic assumptions about basal mechanical conditions and the location of the margins, such as free slip inside the ice stream, are likely to lead to model results containing stress singularities, and consequently to singularities in the basal melt rate.
When basal yield stresses are allowed to vary continuously with position and the ice stream is allowed to choose its own margins, high stress concentrations will occur inside the margins, accompanied by correspondingly high rates of strain heating in the ice. These results contrast significantly with those of previous studies by Reference RaymondRaymond (1996) and Reference Jacobson and RaymondJacobson and Raymond (1998).
The thermal model studied in this paper further suggests that shear stresses become more concentrated over time, and the concomitant increase in strain heating may eventually lead to the formation of a temperate zone of ice near the bed in the margins.
We have also seen that a different choice of sliding law can lead to an entirely different location for the highest rate of strain heating in the margins (bearing in mind the caveat regarding the definition of margin location mentioned earlier), which in turn affects basal thermal conditions and the evolution of the margins. Moreover, a number of processes which may yet prove to be important in the problem of margin migration — basal drainage, formation of crevasses near the bed, the effect of fabric development or of a temperate zone — have yet to be incorporated into theoretical models, and their likely effects are unknown.
From a modelling perspective, one would ultimately like to have a simple, parameterized expression relating the rate of margin migration to simple bulk properties of the ice stream and the surrounding ice sheet, as this would avoid having to solve the ice-flow and thermal problems in the narrow marginal regions explicitly when trying to model the evolution of an ice stream as a whole. However, the physical complexity of the margins makes finding such an expression a difficult, if not impossible, endeavour.
Acknowledgements
I would like to thank R. Greve (Scientific Editor), T. Jóhan- nesson and an anonymous referee in particular for making this paper more readable. G. K. C. Clarke commented on an earlier draft. Financial support in the form of a Killam postdoctoral research fellowship at University of British Columbia is gratefully acknowledged.
