1 Introduction
Several studies have approached questions on flow patterns in ice sheets near the summit (Reference ReehReeh, 1988; Reference Dahl-JensenDahl-Jensen, 1989) and the transition from ice-sheet flow to ice-shelf flow (Reference HerterichHerterich, 1987) by including the role of deviatoric stress gradients. Van der Reference Van der VeenVeen and Whillans (1989) and Reference Van der VeenVan der Veen (1989) have presented a numerical scheme that solves the equations of force balance without any mathematical approximations. They start their iteration from a “first approximation” that formally solves the equations of force balance in a form today referred to the shallow-ice approximation and improve on this approximate solution. Their fixed-point iteration scheme seems to work efficiently for ice sheets, though no account is given about its performance in the case of smaller and steeper glaciers. Reference MullerMuller (1991) proposed a similar approach for solving the first-order equations for stress and velocity components in the two-dimensional plane-flow approximation. This algorithm, based on finite-difference discretization along the x axis (method of lines) and a fixed-point iteration scheme to meet the surface-boundary conditions, proved to be very simple and flexible, and could be readily extended to the three-dimensional case (Blatter, 1995).
The fixed-point iteration applied by Reference MullerMuller (1991) and Blatter (1995) converges rapidly for large ice-sheet configurations with a small aspect ratio. However, for small glaciers, where corrective terms become important, the number of necessary iteration steps grows rapidly with the increasing aspect ratio. Furthermore, the number of necessary iteration steps also increases with increasing horizontal grid resolution. These factors limit the applicability of fixed-point iteration in the achievable horizontal grid resolution and exclude small and steep glaciers, for which convergence cannot be achieved at all. The fixed-point iteration scheme can handle prescribed basal velocity and it can handle mixed basal boundary conditions, prescribed basal velocity in some parts and basal shear traction in other parts of the bed. However, with mixed boundary conditions, convergence of the iteration is so poor that fixed-point iteration in practice is restricted to prescribing the basal velocity along the entire bed. However, this procedure then suffers from obvious deficiencies of adequate imposition of boundary conditions.
In this paper, an improved algorithm is presented. The replacement of fixed-point iteration by a Newton iteration scheme substantially reduces the number of necessary iteration steps. Although a single Newton iteration step requires more computations, the total computation time is considerably reduced. Furthermore, the new technique easily allows the imposition of mixed basal boundary conditions. The limitation in the horizontal grid resolution can be reduced by replacing the single- by a multiple-shooting scheme, With this multiple-shooting scheme, the convergence criterion is no longer the limiting factor in many applications. However, the large memory required for handling the linear algebra limits its applicability.
The problem treated in this paper does not handle the complete glacier-modelling problem. The acceleration term can be omitted in the momentum-balance equations, though not, as sometimes staled, because acceleration is negligibly small (Whillans, 1987). On the contrary, acceleration is large enough that the lime needed to adjust the velocity field to a new situation is negligibly short (seconds) compared with the time-scales that aie relevant for glaciological studies (hours to years). Therefore, stress and velocity fields are treated as quasi-stationary, and that has an important consequence. The stress and velocity fields only depend on the instantaneous stale of the geometry and other conditions that are relevant for the rheology but not their rates of changes, In this sense, the mechanical problem treated in section 2 constitutes a complete mathematical problem that can be solved independently of mass balance and consequent changes in geometry.
This paper presents the methods of multiple-shooting and Newton iteration schemes, and the handling of different types of basal boundary conditions. Here, we discuss the two-dimensional equations, even though the three-dimensional problem closely follows the same line as the simpler case. The mathematical details, important information on the implementation and performance of the routines arc presented by Colinge (paper in preparation), together with practical examples. A practical application of the three-dimensional, first-order scheme on Haut Glacier d’Arolla has been presented by Reference Hubbard, Blatter, Nienow, Mai and HubbardHubbard and others (1998). The performance of multiple-shooting the Newton iteration scheme has been demonstrated by Reference Blatter, Clarke and ColingeBlatter and others (1998). In this paper, referred to as Paper II in the remainder of this paper, the relation between basal stress distribution and basal movement is investigated in detail.
2 Stress and velocity fields
Governing equations
Here, a summary of the model equations in the plane-flow approximation is presented. A scale analysis is conducted that suggests a perturbation-type solution procedure of the governing equations in terms of a small parameter. Correspondingly, equations are called zeroth, first- and higher-order equations depending upon whether they contain the small parameter in the zeroth, first or higher powers. The second-order equations of the conservation of mass and momentum and the How law (constitutive relationship) are listed and a detailed description of the integration scheme as well as the Newton iteration scheme for its solution in the two-dimensional plane-flow approximation are given.
The geometry of the ice mass is defined by the upper free surface z = S(x) and the basal surface z = B(x), in which (x,z) are orthogonal Cartesian coordinates with the z axis pointing opposite to the direction of gravity. As is common in glaciology, the strain-rate tensor D and the Stress deviator Σ are related by
where the rate factor A is, in general, a function of the temperature but is held constant here where considerations are limited to isothermal conditions. The fluidity F depends on the second invariant of the stress deviator, the traditional choice being F(III) = (III)(n-1)/2 with a proven law exponent n≈3. Regularity of the integration procedure requires a finite viscosity law, requesting F(0)≠0, different from Glens flow law.
The total stress T will be decomposed into pressure P and the stress deviator Σ,T = -P1 + Σ, normal components will be denoted by σH (ηα summation, k = i\ z) and shear-stress components by ru :l.j = χ,ζ), while τν, arc ι be total normal stress components.
To obtain an objective hierarchy of approximations, we scale the spatial variables x, z and S, the horizontal and vertical velocity components u and W, the stresses, the rate factor and the fluidity by
where the tilded quantities are dimensionless and {Φ} is the dimensional order of magnitude for a dimensional quantity Φ. Typical vertical and horizontal extents of the ice mass are {H} and {L}, respectively, and the aspect ratio is ε = {H}/{L}; ρ and g are the constant density of ice and acceleration of gravity. The horizontal velocity component U is scaled with the shear velocity in a homogeneous parallel-sided slab of thickness {H} when How law in Equation (1) with its power-law fluidity, F, is imposed: the vertical velocity component w is taken to be a factor f smaller. Pressure and total normal stress, P and Tkk are scaled with the hydrostatic pressure al the bed of a parallel-sided slab, whilst shear stress τij (i ≠ j) and normal deviatoric stresses σkk are sealed with the basal driving stress.
For glaciers, the aspect ratio? is small and can be used as an order-of-magnitude estimate of the various terms in the mass- and force-balance equations and the constitutive equations. In plane flow, these laws take the form (see Blatter, 1995)
in which for n = 3
and t0 is proportional to the inverse of the viscosity ai vanishing effective stress. Note that Equation (8) is obtained from the horizontal force balance by an integration over depth and incorporating the vertical force balance and free-surface stress-free boundary conditions (sec e.g. Hut ter 1983, chapter 5).
We now apply a terrain following coordinate transformation by mapping the local ice thickness h on to the unit interval
We also employ the notation
Journal of Glaciology where f is the shear traction parallel to the surface defined by ζ = constant. In terms of the new independent variables, the transformed T term takes the form
where
With the above transformation rules, the field Equations (7)—(10) may be written in the form
This is a set of four equations for the variables iï, w. f and σ and, if it were not for the T term, it would be a system of first-order partial differential equations. This T term makes it a system of integro-dififerential equations, which is of second order in the derivative. Solutions to it are sought in the domain bounded by the free and basal surfaces, for which boundary conditions need to be imposed. In the transformed terrain following coordinates, the locations of then imposition arc ζ = 0 and ζ = 1.
Physically, the boundary condition at the free surface consists of a vanishing shear traction (vanishing horizontal air-pressure gradient is already incorporated in Equation (8))
At the basal surface, the form of the boundary conditions depends on the local conditions that arc postulated to apply. Either one assumes no-slip conditions and then has
or a certain sliding law. In this latter case, the basal shear traction η, and the sliding velocity ûb are functionally related via
But, in addition, the tangency condition of the flow must hold
The simplest case is perfect slip for which tb=0. It is obvious that a total of three boundary conditions must be imposed, either Equations (23) and (24) or else (23), (25) and (26).
These facts suggest the field Equations (19)—(22) should be regarded not as a system of partial differential equations but rather as though they were ordinary differential equations to be integrated in the ζ direction. (This will be formally achieved by replacing all derivatives by finite differences in the following section.) We have made this interpretation implicitly apparent in Equations (19) -(22) by writing the terms involving the ζ derivatives on the lefthand side and all other terms on the righthand side. Viewed this Way, Equations (19)-(22) are three first-order ordinary differential equations (ODEs) in ζ and one algebraic equation for v, W, Tand σ, the Ç integration of which indeed requires three boundary or initial conditions as surmised above. Because one condition is imposed at ζ — 1 and two are imposed at ζ = 0. i.e. at the two end points of the interval ζ € [0.1], the problem is called a two-point boundary-value problem.
How this integration is performed will be explained below. Here, we simply mention that, as e -> 0, the above equations reduce to the zeroth-order equations, generally known as the shallow-ice approximation, which is commonly used in models of large ice sheets. Eor smaller mountain glaciers, the shallow-ice approximation misses significant physics. Deleting terms that contain ε2 yields the first-order approximation which already captures the essential patterns of the stress and velocity fields, and gives realistic results for small ice masses.
Line integration
The idea of discretizing the partial differential equations in every dimension except one to obtain a system ol ordinary differentia] equations is called the method of lines (Reference Verwer and Sanz-SmiaVerwer and Sanz-Serna, 1984). By introducing a discrete grid on the | axis and approximating the derivatives by finite differences, Equations (19)—(21) can be rewritten as ODEs, and Equation (22) becomes algebraic. For each vertical gridline i, this establishes a sot of three ODEs for the three unknowns f i, ti; and i/'i and one algebraic equation for if·,. This large set of ODEs can be integrated by using a standard numerical integrator I e.g. a Runge-Kutta scheme). The integration begins at the base (Ç = 0), with starting values for f|, and «b at each gridpoint. At each step of the numerical integration, the algebraic equation is solved explicitly or with a numerical root-finder. It is important that the algebraic équation always has a unique real solution. This is the case for small aspect ratios generally relevant for glaciers and for the How law in Equation (11) with a flow-law exponent n = 3 (paper in preparation by J. Colinge).
To arrive at a proper finite-difference scheme, several points must be taken into account. The shear stress t is computed from Equation (20) and depends on ?σ/? and σ. Furthermore, σ is computed from Equation (22) and depends on dù/?. This cascade of dependence on derivatives reduces the order of the whole difference scheme to p- 1, even if all difference schemes for single derivatives are of the order p. Thus, to obtain a consistent difference scheme of order p, the Θ?/? in Equation (22) must be dis-cretized to order/; + 1 (paper in preparation by. J. Colinge). The integration from Ç = 0 to ζ = I with the starting values for basal shear stress and velocity components, i.e. three conditions, does not automatically satisfy the surface boundary conditions in Equation (23). In order to solve the boundary-value problem, the proper basal values for ti,b can be found iteratively. A good initial choice is the shallow-ice approximation of the basal shear stress, which in many eases is already close to the solution:
If basal conditions have a large variability, then an initial f? = () is often the belter choice to start with. With the basal shear stress and the values for the basal velocity components, Equation (22) is first used to calculate the basal value for σ·φ. Integrating upwards from the base yields surface values ffs =£ 0. By using this result, a correction to fUt must be found to obtain t\fi closer to the required boundary condition. This method is called shooting, or single shouting (Reference Ascher, Mathrij and RussellAscher and others, 1988). The type of basal boundary conditions that must be met is defined by the sliding parameterization. Prescribing basal velocity is the simplest and numerically most stable case, except for the no-slip conditions of the limited value for studies of basal sliding. Mixed boundary conditions, vanishing basal velocity at non-sliding locations and basal friction represented by basal shear traction at sliding local ions may be useful for specific studies on basal conditions. The application of a sliding law which relates basal shear traction to basal velocity requires a nested iteration of this sliding law within the shooting procedure.
The T term (see Equation (17)) on the righthand side of Equation (20) requires different treatment. This term contains the integral off and, because it is of the second order in the derivatives, it can only be computed after the hue integration has been completed. In practice, the T term is computed by using the f from the previous iteration step and then including it in Equation (20). Its inclusion makes the computation lime per iteration step slightly longer but the number of necessary iteration steps is hardly affected. Numerical experiments show that the convergence becomes more difficult if nested iterations such as for sliding laws or for the T term must be considered.
Single-shooting fixed-point iteration
A simple way to define a single-shooting procedure is to write the problem, Equations (20)-(22), as a fixed-point problem (Reference MullerMuller, 1991; Blatter, 1995). We introduce the following notation, where M is the number of gridpoints for the method of lines:
A converging fixed-point iteration is equivalent to a constructive proof of Banach's fixed-point theorem (Ascher and others. 1988). Thus, for (a > 0. a sufficiently small; the fixed point fb of G, G(f *,) = f*b, is the correct solution. In the simplified ease of an infinite slab with a parallel Hat surface and base, a simple criterion
can be rigorously proven (paper in preparation by J, Colinge) to ensure the existence and uniqueness of the desired fixed point (Banach's theorem). This exact criterion lor the slab is approximately valid for the case of realistic glacier geometries and one can expect to be forced to decrease a either when Δ decreases or the aspect ratio f increases (see also Blatter, 1995).
This criterion is very restrictive and explains the impossibility of applying the method to small and steep glaciers. Moreover, it uses global quantities of ice geometry and local conditions may require a smaller ? achieve convergence. This is especially true if the surface slope locally displays large longitudinal variations, such as in icefalls or at steep glacier snouts. This criterion provides information regarding the stability of the numerical process. An unstable process is one which assigns largely different results to pairs of very close starting conditions. Since the idea of single shooting is to use the final results to learn something about the. starting conditions, it is essential to avoid instability.
The restrictive criterion in Equation (311 also indicates the impossibilily of applying (he line integration in a single-shot inverse mode. One may be tempted to start the integration at the surface with zero shear traction and measured velocity components to compute basal shear traction and basal velocity (Reference Van der VeenVan der Veen and Whillans, 1989). Although the mathematical problem is well posed, it is ill-conditioned and thus basal conditions obtained by this integration are extremely sensitive to the input at the surface and may be far away from the correct solution {see e.g. Reference LliboutryLliboutry, 1987; Reference Bahr, Pfeffer and MeierBahr and others, 1994).
Single-shooting Newton iteration
From the above consideration, we look for T\, such that Ts(fb) = 0. Colinge (paper in preparation] has shown that fs is infinitely differentiable with respect to ti This is only true if the flow law F is differentiable and invertible in the entire domain; otherwise, the uniqueness of the solution is no longer guaranteed, although the numerical integration may produce a result. Thus, Glen's How law needs lo be modified to avoid infinite viscosity for stress-free conditions (Reference HutterHutter, 1981, Reference Hutter1983; Reference Smith and MorlandSmith and Morland. 1981).
Another way of implementing a single-shooting procedure is by solving f s(f b) = 0 with a New ton iteration:
Then
is a series expansion of Ts, where
is the jacobian matrix. To define the Newton iteration, we need to truncate Equation (33) and solve the remaining linear equation
and thus find the basal f,'' for the next from the old iteration step
The Newton iteration uses all gridpoints at the surface to correct the stress values at each gridpoint at the base and converges quadratically close to the solution (Ascher and others. 1988). On the other hand, the fixed-point iteration only uses the gridpoints al the surface to correct the basal stress in the saine vertical profile and converges linearly. T> compute the components of the jacobian
where
we have to integrate the set of equations M - 1 times, once for each choice of the indices j = 1,…, M - 1. The small increment δT in the corresponding basal shear stress should be chosen as small as possible to obtain accurate values of the jacobian components. The smallest and best choice for δT is the square root of the machine precision.
For large M, the above computation of the jacobian may become very expensive in computing time, since we only get one column of the matrix per integration. However, there is a great potential for reducing the number of necessary integrations. For physical reasons, we can anticipate that the components of the jacobian become smaller if the difference | j — i | grows, thus thejacobian is approximately a banded matrix. The Newton Kantorovitch theorem (Reference Ascher, Mathrij and RussellAscher and others, 1988) ensures the convergence of the Newton iteration even with an approximate Jacobian, provided the approximation is close enough. This suggests that we assume the jacobian is a k diagonal matrix, where k is odd. a nd then calculate the approximate components of the jacobian by
for I = 1,…, k and j =l,l + k,l + 2k,…. This enables us to compute M/k columns of the jacobian matrix at each integration and thus reduces the number of necessary integrations to k. independently of the generally much larger number of gridlines M, Moreover, the linear system, Equation 135), can be solved faster by using a banded matrix. The bandwidth k must be chosen such that k ? > 5/ι, i.e. beyond the range of influence of basal perturbations.
When the convergence threshold is not reached with the first iteration step, we repeat the above procedure of computing the jacobian and integrate the equations with an improved choice of the basal shear traction. Sometimes the Newton iteration can be continued by keeping the old jacobian, although convergence of the iteration may become slower than with an updated Jacobian or may not be achieved at all.
It is also possible to compute a number C single, called the condition, which describes the stability of the process (Reference Ascher, Mathrij and RussellAscher and others, 1988). For a slab, this number is
A large condition number indicates an unstable algorithm which is very sensitive to initial (basal) values. In many cases, we observed numerically that, for a chosen geometry (synthetic or realistic),Δε could be chosen about half the size than for the fixed-point iteration. This is a clear improvement, interesting for many studies, but the spatial resolution is still limited.
Multiple-shooting Newton iteration
The computation becomes unstable as ? decreases. One possibility of stabilizing the iteration scheme is to divide the glacier into N horizontal layers, applying the single shooting with Newton iterations to each layer and require that the final solution has no discontinuity. More precisely, with the starting shear stresses and velocities f{),p (/.,, 3 = 1,···,N (f7n.| is given) and the linal values T\.j{T{ij. UQ,J), UIJÇTQJ, U()j) we have to solve the following problem
in terms of Toj, UQJ. This approach is called multiple shoaling (Reference Ascher, Mathrij and RussellAscher and others, 1988) and it is possible to prove for the slab that
This is a clear improvement, since Cmultiple < Csingle and the stability improves for increasing N. It is also possible to prove fora slab with a linear flow law that such a process converges to the unique solution of the partial differential equation (paper in preparation by J. Colinge).
Again, the analysis of the simplified case of the slab has been relevant in the general ease. So, this algorithm theoretically allows for any spatial resolution but requires much more computing time and memory (if Δ is halved, then N must be doubled). The fact (hat the stability of the algorithm suffers, if Δ decreases is a property of the method of" lines discretization. The condition number of the ODE system is
which obviously grows as Δ becomes smaller.
Mixed basal boundary conditions
The above algorithms and analyses do not only apply for prescribed zero basal velocity; any reasonable basal velocity can be imposed. However, in some cases it may be preferable to prescribe the basal shear traction in a given area of the glacier instead of the basal velocity itself. The algorithm still works in this more general application with little modification.
Again, it is possible to give a fixed-point formulation of the problem
where Yih = τ,.ι, if i G E where i K denotes the domain in which basal shear traction is prescribed) and YJj, = ?,Λ, if i £ E, Ta = (f i.s. ? ? ·. fM-i,s) · The iteration parameter a must be chosen differently for the two domains with either prescribed basal velocity or prescribed basal shear traction. However, the fixed-point iteration was found to converge very slowly and, thus, its application in the case of mixed basal boundary conditions is not practical.
The resolution with a Newton iteration is preferable and, because the problem for mixed basal boundary conditions is less stable, multiple shooting is necessary. The algorithm must be prepared to include ?/,.ι, as an unknown if i g E and f/.i, il i £ E. This implies computing a mixed Jacobian matrix by considering part ial derivatives with respect to «,.ι, il"/' € E and with respect to τ,.ι, if/ £ E.
In the basal section, where shear traction is prescribed. the basal velocity becomes a solution of the computation, This fact poses an additional problem with the discretization scheme. Using, for example, a symmetric second-order difference for ??/? makes the σ at odd gridpoints dependent on σ at even gridpoints and vice versa. Additionally, a symmetric second-order scheme for ?σ/? transfers this uncoupling to the solution off.
This uncoupling is less effective as long as the velocity component u is prescribed along the entire bed and solutions are acceptable as long as variations in the basal velocity are small. However, if basal shear traction f\, is prescribed instead al parts of the bed, the decoupling produces solutions for even and odd gridpoints which oscillate with a wavelength of two gridcells, There are several possibilities for removing part of the oscillation, for example, by smoothing the solution with a filter function. However, the meaning of such a smoothed solution is not clear. A variation of this method is to smooth the solution for the basal velocity u b, within the region of prescribed basal shear traction and I hen to recompute a solution by using the algorithm prescribing this u b, instead. This furnishes a well-defined numerical solution to an approximate problem and the well-posedness of the problem suggests that this solution should be close to the true one. In practice, however, this method was not satisfactory. After smoothing u b the corresponding tb was not as close to the originally prescribed one as desirable.
A preferable approach is to avoid the symmetries that disconnect the odd and even points, for example, by using non-symmetric finite differences to approximate ?ιι/θ. There are several ways to do this but some ways remove only part of the oscillation (the odd and even points become only-weakly coupled). The following scheme removes the oscillation to a large degree and also maintains the consistency in the order (second order) of the whole difference scheme. A fourth-order symmetric scheme for
is first used to compute cre(Xi,Ç). The corresponding ^(^ivC)jfi(C)) 's 'hen used to compute the value of ù,(C + Δ<″). The value of f;(C + ?ζ) is dependent on à, as well as on dà/?. For à, the above σν{χ,. ζ) is used but
where crteft,i+i (C) is recomputed from Equation (22) with the asymmetric order-three scheme
and σ,?,ι,,./ ι (0 with
This makes three separate computations ????/? and solutions of Equation (22) necessary. On the other hand, this scheme has a well-defined order and meets the boundary conditions correctly. If Equation (22) is solved iteratively with a numerical root finder, the first solution, e.g. for σc is a very close approximation for the other two solutions σleft and σright, and thus allows a substantial reduction in the number of iteration steps. Using the average σ,. = (iTriw|M + à\,.(t)/2 produces the same solution (within the precision of the numerical scheme) as using the above suggested a,.
3 Flow in a parallel-sided plane slab
Second- and first-order slab equations
The geometry of the ice slab is defined by the upper free surface S = H and the basal surface B = 0 in Cartesian coordinates (x.z), with the x axis aligned parallel to the direction of the surface slope. The z axis and the direction of gravity include an angle ?. The equations are written in dimensionless form using the scaling
where t,? and w are the shear stress, deviatoric normal stress, longimdinal and transverse (vertical) velocity component, respectively. The symbols with the "tilde" are the corresponding dimensionless quantities; the rate factor, A is scaled with A 0 and Glen's exponent n= 3.
The above sealing serves the purpose of writing the equations in dimensionless form and, thus, does not contain the scaling factor ε. However, in the following equations we formally re-introduce factors ε and ε2 with ε = 1 to mark the first- and second-order terms corresponding to the scaling in Equations (2) (6).
The dimensionless, two-dimensional second-order slab equations for the four unknowns f. <r, u and w are
with the flow law in Equation (11) and the surface-boundary condition fs = 0.
In the first-order approximation, Equations (55) and (56) reduce to
Where as Equations (57) and (54) remain as in the second-order form. Note that in the first-order approximation, the solution for to from Equation (54) does not feed back to the solutions of the other three equations; thus, w can be determined by quadrature of Equation (54) after the other three variables are computed.
The foregoing scaling reduces the number office parameters for a homogeneous slab to two, namely i0 and AQ. The rate factor A0 is simply a multiplier for the velocity held and, as long as A 0 is homogeneous over the whole domain, the stress field is independent of Ay,. This permits normalization of the velocity field for different choices of tu and makes A0 a function of iy, reducing the number of free parameters to the single parameter t(V In all examples, we choose tQ> and AQ such that the asymptotic surface velocity
The asymptotic stress and velocity profiles are defined as the solution for the entirely non-sliding slab, i.e. the solution in the absence of any sliding perturbations. The limit of A0 —> 0 constitutes a Newton fluid, whereas the limit-2 t0 —? 0 constitutes the classical Glen How law with infinite viscosity for vanishing effective stress. In all presented numerical examples, we choose tn = 0.1 and Arj = 5/.'î. With this assumption, the only freedom rests in the choice of the basal stress pattern over a chosen area and/or the velocity over the remaining areas.
To simulate an infinitely long plane-strain slab, the model domain of finite length is numerically connected end to end. With homogeneous basal boundary conditions in the finite domain, this produces a homogeneous solution for the infinite slab; with inhomogeneous boundary conditions, this yields a periodic pattern with a wavelength of the length of the finite domain.
Comparison between first- and second-order solutions
The differences between first- and second-order equations arc terms that contain derivatives of f and w with respect to x. These terms may become significantly large if sharp longitudinal transitions occur, either abrupt changes in the surface of bottom topography or sliding-non-sliding transitions. Only the latter are possible for parallel-sided ice slabs.
We have performed a number of computations with the first- and second-order models for a long non-sliding slab with a shorter zone of perfect slip. One representative example is illustrated in Figures 1-6. The chosen grid size is Δx = 1.0 and the length of the zone with frictionless sliding is defined by seven gridpoints with prescribed vanishing basal shear traction, tb= 0, and the no-slip condition is imposed in the rest of the slab by ub = Q The transition occurs over one gridcell, i.e. over the distance of one slab thickness. The non-dimensional basal shear traction scales to τ>, = 1 for non-sliding conditions far away from the influence of the sliding zone (asymptotic behaviour) and the corresponding shear velocity is normalized to Ûs = 1.
In all cases, the basal shear traction reaches a peak value at the non-sliding gridpoint adjacent to the sliding zone and then drops sharply to the prescribed zero shear traction at the first sliding gridpoint (Fig. 1). Correspondingly, the longitudinal velocity component increases sharply, also changing from a typical shear-flow profile to nearly plough-flow in the interior of the sliding area (figs 2 and 3), The transitions are also accompanied by a peak in downward vertical velocity where the ice enters sliding, and upwards, where ice leaves the sliding area (Fig. 4), Of course, (he surface of the slab would change corresponding to the non-zero vertical velocity component and given. surface mass flux. In this study, we do not investigate surface response to variations in basal conditions but wc look at snapshots of the quasi-stationary stress and velocity fields with sharp basal transitions for a given geometry. Such transitions may be unrealisti-cally abrupt (Reference HutterHutter and Olunloyo, 1980, 1981; Reference Banunn and MacAyealBarcilon and MacAyeal, 1993) and, in reality, ice may break at the peak stress locations. On the other hand, this illustrates the extent to which the first-order solution captures the essential patterns of stress-strain relations and produces reliable results relatively close to the second-order solutions.
Fig. 1. Longitudinal profiles of basal shear traction in the first- and second-order approximations (terms linear and quadratically in e are accounted for) across the transitions between non-sliding areas and a limited zone of perfect sliding for a parallel-sided ice slab. The sliding is specified for seven gridpoinls (2-8) by prescribed vanishing basal shear traction. The non-sliding parts are defined by vanishing basal velocity. The grid size is Ax = I which corresponds to the thickness of the slab.
Fig. 2. Longitudinal profiles of basal and surface longitudinal-velocity component in the first-and second-order approximations across the transitions between non-sliding areas and a limited zone of perfect sliding as in Figure 1. The subscripts S and b refer to surface and base, respectively, ike solid and dashed lines represent first-order solutions and the dotted and clash-dotted lines the second-order solutions.
The second-order solutions show profiles with negative shear stress in the interior of the sliding area [Fig. 5}. It is not exactly clear whether this is physically realistic or whether it is a numerical artefact, or both. In the solution for f, this looks like a small oscillation around zero-stress values with a wavelength of two gridcells (Fig. 6). This fact strongly indicates that the oscillation is numerically induced. although the exact reason is not clear. The problem stems from the T term (sec Equation (20) and thus far no solution is available.
Fig. 3. Vertical profiles of longitudinal velocity inside and adjacent to an area of perfect sliding as described in Figure 1. The labels refer to the corresponding gridpoints as illustrated in Figure 1. The decrease in velocity with distance to the bed is assumed to be a numerical artefact in the second-order solutions.
Fig. 4. Longitudinal profiles of vertical velocity component cil the surface in the second-order approximation across the transitions between non-sliding areas and a limited zone of perfect sliding as in Figure 1.
Fig. 5. Vertical profiles of shear stress inside and adjacent to an area of perfect sliding as described in Figure I. The labels refer-to the corresponding gridpoints as illustrated in Figure i. 7 he negative shear-stress profiles are assumed to be a numerical artefact in the second-order solutions.
Fig. 6. Longitudinal profiles of shear stress at various levels of z in the second-order approximation across the transitions between non-sliding areas and a limited zone of perfect sliding as in Figure L The negative shear stresses are assumed to be a numerical artefact in the second-order solutions.
4 Conclusions and prospects
The governing equations for the two-dimensional stress and velocity Gelds in glaciers are presented in Cartesian and terrain-following coordinates. A sealing analysis allows the definition of a hierarchy of approximations based on powers of the small-aspect ratio of the iee geometry. Deleting terms that contain c to powers 2 and higher leads to a set of three coupled non-linear partial differential equations, referred to as first-order approximation in this paper. We present discretization schemes based on the method of lines and discuss various integration and iteration schemes.
The single-shooting fixed-point iteration is easy to program and, for large ice masses with small aspect ratios, it converges reasonably fast (Blatter, 1995). At present, this is the method of choice for three-dimensional finite-difference modelling of higher-order glacier mechanics and dynamics. In this paper, we introduce two-dimensional (plane-strain) models with more elaborate numerical methods: Newton iteration insiearl of fixed-point iteration and multiple shooting instead of single shooting.
The Newton iteration scheme is faster than the fixed-point scheme Cor the same glacier and the same horizontal résolution Ax. For large iee sheets with smooth surface topography, the Newton iteration allows for substantially smaller Ax than the fixed-point iteration. For small glaciers with smooth topography, this improvement is only small but, for glaciers with sleep sections and variable surface profile, the fixed-point iteration performed better. Fixed-point iteration only works reasonably well, if basal velocity is prescribed along the entire bed. On the other hand, Newton iteration can adequately handle mixed basal conditions, where the velocity is prescribed only at parts of the bed and basal shear traction at the rest of the bed. This allows the modelling of glaciers with floating parts such as on sub-glacial lakes or ice shelves.
Journal of Glaciology The multiple-shooting scheme allows much better spatial resolution than the single-shouting scheme, with Δx as small as 5% of the average ice thickness for the first-order approximation. However, this is at the price of long computation times and large work-memory requirements. An application, which corroborates these statements is presented in Paper 11. This scheme may still be practical, if single situations tor the stress and velocity fields arc to be modelled. If the mechanical model must be used repeatedly for coupling a surface-evolution model, the single shooting with a well-implemented Newton iteration is certainly the method of choice. In the second-order approximation, the T term must be iterated parallel to the Newton iteration. This reduces the convergence radius and, thus far, the best horizontal resolution that could be achieved with the second-order code was about 75% of the average ice thickness. More elaborate discretization schemes, which promise faster computation with higher spatial resolution suitable for three-dimensional ice masses, are presently being investigated.
The physics that is included in the glacier model presented here poses a novel problem concerning the required information on the glacier to exploit fully the capabilities of the model. The bed topography becomes available for a growing number of glaciers. The general lack of information on sliding patterns poses the most serious limitation, since the longitudinal coupling of stresses strongly depends on the spatial variations of the basal velocity; see Paper II. On the other hand, these requirements point to the desirable field measurements thai provide the necessary input for their subsequent interpretation based on model computation (Reference Harbor, Copland, Hubbard and NiatrHarbor and others, 1997; Reference Hubbard, Blatter, Nienow, Mai and HubbardHubbard and others, 1998).
5 Acknowledgements
We thank A. Ohmura and G. Wanner for their support and stimulating discussions. K. Hotter commented on an earlier version of this paper and helped to improve the manuscript substantially.
