1. Introduction

The solution of the complete ice-sheet flow equations is a complex problem due to the thermomechanical coupling of the equations, anisotropy of ice associated with the formation and evolution of fabric, and the presence of moving boundaries changing the geometry of an ice sheet. For this reason, a number of simplifications are adopted in ice-sheet models in order to reduce the complexity of calculations. The most common simplifications are: (1) the reduction of geometrical dimensionality by taking advantage of symmetries in a flow field and solving either a plane or an axi-symmetric problem, and (2) the assumption of isotropic behaviour for the ice. Such an approach, in which Stokes’equations for an incompressible viscous fluid are solved with the constitutive relation in the form of Glen’s flow law, was employed, for instance, by Reference Hooke, Raymond, Hotchkiss and GustafsonHooke and others (1979), Reference RaymondRaymond (1983), Reference HodgeHodge (1985), Reference HansonHanson (1995) and Reference HvidbergHvidberg (1996).

A fundamentally different method is to construct an approximate solution to the problem equations by exploiting the long aspect ratio of natural ice masses, reflecting varying conditions in lateral directions compared to the normal direction through the ice thickness. Such an approach, which constitutes the shallow-ice approximation (SIA) theory (Reference Fowler and LarsonFowler and Larson, 1978; Reference Morland and JohnsonMorland and Johnson, 1980; Reference HutterHotter, 1981,Reference Hutter1983; Reference MorlandMorland, 1984), uses scaled variables based on typical physical magnitudes and provides explicit relations for the velocity and stress fields in ice sheets with slowly varying bed topography and free-surface elevation. The first application of the SIA to numerical modelling of ice sheets is due to Reference Hutter, Yakowitz and SzidarovszkyHutter and others (1986); since then the SIA has been widely implemented, with various modifications, in many large-scale ice-sheets models (Reference Hindmarsh, Morland, Boulton and HutterHindmarsh and others, 1987; Reference HerterichHerterich, 1988; Reference Dahl-JensenDahl-Jensen, 1989; Reference HuybrechtsHuybrechts, 1990; Reference Fabre, Letréguilly, Ritz and MangeneyFabre and others, 1995).

Although the SIA has proved to be a very effective tool in ice-sheet modelling, its application is restricted due to the assumption of isotropic ice response. For this reason, an attempt has been made recently by Reference Mangeney and CalifanoMangeney and Califano (1998) to extend the SIA applicability to anisotropic ice. These authors have analysed the problem of ice flow in the ice-divide region by adopting a transversely isotropic fabric corresponding to that in the Greenland Icecore Project (GRIP) ice core. The same fabric has also been used with the complete set of mechanical equations by Reference Mangeney, Califano and CastelnauMangeney and others (1996, Reference Mangeney, Califano and Hutter1997) to solve numerically the problem of a steady-state flow of ice under isothermal conditions. In the latter papers, however, the empirically derived anisotropic fabric is, in fact, a function of the ice depth only, since no constitutive law that relates the fabric evolution to the flow field is included in the model.

A proper approach requires that the ice response is described by a frame-indifferent constitutive law that relates the stress to a limited number of variables representing the deformation and the evolving structure (fabric) of ice. Recently several constitutive models describing strain-induced anisotropy of polar ice have been developed, but here we focus on the theory developed by Reference Morland and StaroszczykMorland and Staroszczyk (1998) and Reference Staroszczyk and MorlandStaroszczyk and Morland (2000), in which it is assumed that the macroscopic response of ice can be described in terms of fabric induced entirely by macroscopic deformation, and all microscopic processes occurring at the grain level are ignored. It is also supposed that the induced anisotropy depends only on the current strains and does not depend on the deformation history, which ignores the effects of crystal interactions which may depend on the nature of the deformation process, and therefore induce different fabric for different deformation paths. It is believed, however, that this approximation is the most simple approach to an evolving anisotropic viscous law which can be implemented in large-scale ice-sheet modelling, since it requires that only current deformation gradients are calculated in addition to the velocity and pressure fields. The adopted form of an orthotropic constitutive law expresses the deviatoric stress in terms of the strain rate, strain and three structure tensors defined by the current principal stretch axes. The relation is separable in the isotropic dependence on strain rate and fabric dependence on deformation, and in its simplified form involves two fabric-response functions with dependence on the principal stretches and an invariant measure of total deformation. The constitutive law is used to analyse a plane, gravity-driven, steady flow of an isothermal ice sheet without making the shallow-ice approximation. Assuming a horizontal bedrock with no-slip basal conditions and a fixed, free-surface elevation, ice-accumulation rates needed to maintain the prescribed geometry, together with the velocities and stresses, are calculated by using a finite-element method (FEM). Numerical simulations have been conducted for different combinations of model parameters in order to investigate their effect on the flow field. The results obtained for anisotropic ice are compared with those for isotropic ice (modelled as a Newtonian fluid) to show how the anisotropy influences the velocities and stresses in the ice sheet.

2. Problem Formulation

We consider a plane-strain, gravity-driven, steady flow of an ice sheet over a rigid horizontal bedrock (Fig. 1). The problem is solved in rectangular Cartesian coordinates Oxyz, with the horizontal plane Oxy coinciding with the flat bed, the z axis directed upward and passing through the ice divide, and the x axis in the direction of flow to the right of the divide. For simplicity, the ice-sheet profile is assumed symmetric about the plane x = 0.

Fig. 1. Cross-section of an idealized symmetric ice sheet.

Let σ be the Cauchy stress tensor, with components σ_xx , σ_zz and σ_xz in the plane of flow, and v be the ice velocity vector, with components u, v and w in the x, y and z directions respectively (here, for the flow in the O_xz plane, v = 0). The deviatoric stress σ ′ is defined in terms of σ and the mean pressure p by

where I is the identity tensor, and tr σ denotes the trace of σ . The deviatoric stress σ ′ is determined by the viscous constitutive law, described in section 3, while the pressure p, being here a workless constraint stress compatible with the incompressible response, is not determined by the flow law, but by the momentum balance and boundary conditions.

In the absence of inertia forces, due to a negligibly small Reynolds number for the gravity-driven ice-sheet flows, the momentum-balance equations become the equilibrium relations

where ρ is the ice density, assumed constant, and g is the gravitational acceleration. The mass balance is here the incompressibility condition div v = 0, which in components reads

The equilibrium and incompressibility relations (2.2)-(2.4) have to be supplemented by the boundary conditions. At the undeformable bedrock it is assumed that the ice temperature is significantly below the melting point so that there is no melt and no basal sliding at the interface. Thus the conditions at z = 0 are

In view of the symmetry of the flow, at the ice divide x = 0 the relations

hold. At the free surface z = h(x,t) the kinematic condition is

where t denotes time, U_s and w_s are, respectively, the free-surface horizontal and vertical velocities, and q = q(x,t) is the accumulation rate (ablation if q ≤ 0). The free surface is assumed to be traction free, that is the atmospheric pressure is considered negligibly small compared to typical stress levels in the ice sheet. Hence,

where n_s is a unit outward vector normal to the free surface. In the flow problem considered here we assume that the free-surface elevation is prescribed and stationary; that is h(x,t) = h(x) and ∂h/∂t = 0 in (2.7), and we calculate the ice velocities and the associated accumulation rate q(x) required to maintain the fixed ice-sheet geometry. These quantities will then be compared for the isotropic and anisotropic solutions.

3. Constitutive Law

The anisotropic viscous behaviour of ice is described by an orthotropic model formulated by Reference Morland and StaroszczykMorland and Staroszczyk (1998) and further extended by Reference Staroszczyk and MorlandStaroszczyk and Morland (2000). The adopted flow law expresses the deviatoric stress σ ′ in terms of the current strain rate D, the left Cauchy-Green strain tensor B and three structure tensors M r ( r = 1,2,3) defined by the current (evolving) principal stretch axes. Here we denote the spatial rectangular Cartesian coordinates by Ox_i, (i = 1, 2, 3), with the equivalence x ₁ = x, x ₂ = y, x ₃ = z, and, for the velocity components, v ₁ = u, v ₂ = v, v ₃ = w. Let the particle reference coordinates be OX_i (i = 1, 2, 3 ), with the equivalence X ₁ = X, X ₂ = Y and X ₃ = Z. Then the deformation gradient, F, spatial velocity gradient, L, and the strain rate, D, have components

and the deformation gradient is determined by the kinematic relation

where the superposed dot (’) denotes the material time derivative and the summation convention for repeated subscripts is applied. In practice, (3.2) must be solved simultaneously with the momentum balance, incompressibility condition and constitutive law, and is subject to an initial condition that F = I when the ice is first deposited at the surface. The strain B, its principal stretch axes given by unit vectors e ^(r) ( r = 1,2,3), and the squares of the principal stretches br(r=1,2,3) are defined by

where λr are principal stretches along the vectors ε r, and, due to the ice incompressibility,

The orthotropic structure tensors M (r), describing reflexional symmetries in the principal stretch planes normal to the axes ε (r), are defined by

Following Reference Staroszczyk and MorlandStaroszczyk and Morland (2000) we express the orthotropic viscous law in a form for a’ in terms of D and B and involving two fabric-response functions, each depending on a single deformation invariant argument:

where μo = μo(trD²,T) defines the isotropic fluid viscosity at temperature T when B = I, and R is an invariant measure of total deformation

The normalization

has been introduced to yield the isotropic fluid viscous law a’ = 2μ0D when B = I, that is when b₁ = b₂ = b₃ = 1 and R = 3. Although the adopted form (3.6) is a considerable simplification of a general orthotropic law (Reference BoehlerBoehler, 1987), this reduced form still retains ample flexibility to correlate with observed viscous responses. This has been demonstrated in Reference Staroszczyk and MorlandStaroszczyk and Morland (2000), where the model was used to simulate the ice behaviour in simple stress and strain configurations, corresponding to those occurring in um-axial compression and simple shear tests conducted in a laboratory. The simulated responses given by the model have also been compared with the predictions of a micro-macroscopic model of Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier (1999) and a good agreement between the results of the two theories, based on very different approaches, has been found (Reference Staroszczyk and GagliardiniStaroszczyk and Gagliardini, 1999).

The viscous response defined by the law (Equation (3.6)) is determined by the two fabric-response functions f(b_r) and G(K). It is expected that these should be simple, smooth functions, and it has been deduced (Reference Staroszczyk and MorlandStaroszczyk and Morland, 2000) that f(b_r) should be positive and monotonic increasing, that G(K) should be negative, and that f(b_r) and G(K) are explicitly related so the model is defined by a single-fabric function. This last result follows from an equality between instantaneous viscosities in plane flow at different axial stretches, being one of a set of equalities and inequalities deduced for valid responses at all possible differential stretches. These are based on the assumption that the alignment of ice-crystal c axes towards the direction of compression (and away from the direction of extension) depends on the relative magnitudes of the three principal stretches λi. The smaller a given principal stretch is compared to the other two stretches, the stronger is the alignment of c axes towards the direction of this stretch, and hence the easier is the crystal basal gliding on the plane normal to this principal stretch axis. Thus,

where

which, due to the normalization (3.8), imposes a restriction on f(br) at b r = = 1:

Further, the limit of (3.9) as b1 → ∞ yields the relation

Two additional relations between limit values of f(b_r) and G(K) are obtained by considering um-axial compression and simple shear tests, in which the limit ratios of a fabric-induced viscosity to the isotropic viscosity are measured at large strains. The reciprocal values of these limit ratios are conventionally described as enhancement factors for compression and shear. By applying the viscous law (3.6) to the um-axial, unconfined compression (when both lateral principal stretches are equal), we obtain the relation

where A is the reciprocal of the axial-enhancement factor. For simple shear (when the principal stretch normal to the shear plane is unity), the law (3.6) gives

where S is the reciprocal of the shear-enhancement factor. These two relations (3.13) and (3.14), together with (3.12), determine /(0), f(∞) and G(∞):

The above values ensure that the viscosity ratios for um-axial compression and simple shear at large strains reach the limits which are equal to the prescribed model parameters A and S, but the rate at which the limit viscosities are approached during ice deformation depends on specific cd forms of the response functions. We have explored two sets of simple functions:

where α is determined by the restriction (3.11), m is a free parameter, and r = 1,2,3. Below we illustrate the results predicted by the viscous law (Equation (3.6)) for um-axial, unconfined compression (Fig. 2), and simple shear (Fig. 3). In Figure 2, A_× ≥ 1 is a lateral stretch measured in the direction normal to the axis of compression, and in Figure 3, R ≥ 0 denotes a shear strain measured in the plane of shearing. The curves labelled (1) and (2) correspond to the function (3.16) with m = 1.5 and 2 respectively, and the labels (3) and (4) refer to the function (3.17) with m = 1 and 1.5 respectively These particular values of m have been chosen as giving good correlation with the micro-macroscopic model of Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier (1999); see Reference Staroszczyk and GagliardiniStaroszczyk and Gaghardini (1999). The results presented in the figures have been obtained for the limit viscosity ratios A = 3 and S = 0.2, corresponding to the enhancement factors 1/3 for compression and 5 for shear.

Fig. 2. Evolution of the normalized axial viscosity with increasing stretch λ1 in uniaxial compression for different response functions f(br).

Fig. 3. Evolution of the normalized shear viscosity with increasing strain R in simple shear started from an isotropic state for different response functions f(b_r).

For other combinations of A > 1 and S ≤ 1 qualitatively similar results are obtained; that is, a continuous strengthening of ice (increase in the axial viscosity) during axial compression, and an initial strengthening (for R ≤ 2), followed then by a significant softening of ice, during simple shearing. Such a qualitative behaviour agrees with the observed responses of cold ice subjected to stress levels typically occurring in large ice sheets, for which, according to Reference Pimienta, Duval and LipenkovPimienta and others (1987), the value of A can be as high as 10 for a single-maximum fabric. We have chosen here A = 3, since this smaller value seems to be more relevant to polar ice sheets and corresponds to the value obtained by Reference Mangeney, Califano and CastelnauMangeney and others (1996) for the ice near the bottom of the GRIP ice core. It should be pointed out here that values of A ≤ 1, corresponding to the enhancement factors for compression greater than unity, have been measured for warm ice (near melting) subjected to large strain rates (Reference Budd and JackaBudd and Jacka, 1989), and as such are appropriate for modelling the flows of temperate glaciers, rather than large cold ice sheets. However, as will be demonstrated later, the influence of the axial viscosity ratio A (as opposed to the shear factor S) on the global flow of large ice sheets is very small.

4. Dimensionless Variables

Before solving numerically the equilibrium and mass-conservation equations (2.2)-(2.4) together with the constitutive law (3.6), we eliminate physical dimensions from the equations in a way suggested by Reference MorlandMorland (1984). We adopt a typical ice depth, h*, as a length scale, and a typical accumulation rate, if, as a velocity scale. These two characteristic magnitudes determine other scaling parameters which (assuming p = 918 kg m⁻³ and g = 9.81 m s⁻²) are collected in Table 1.

Using the adopted scales, we introduce dimensionless variables, indicated by a superposed bar, defined as follows:

An essential feature of the ice-sheet analysis is that lateral derivatives with respect to are much smaller than normal derivatives with respect to . In order that maximum and derivatives have equal status, we apply coordinate stretching in the direction, an approach on which the shallow-ice approximation theory is based (Reference HutterHutter, 1983; Reference MorlandMorland, 1984). Hence, we adopt a small parameter ε ≪ 1 and introduce differentially scaled coordinates

By choosing ε to be the aspect ratio defined by the ratio of the characteristic ice depth h* to the characteristic lateral dimension of the glacier L*

both and coordinates become order unity. Since w is order unity by construction, we re-normalize the velocities by

so that is also order unity (which follows from the mass-conservation balance). Typically, the aspect ratio varies from ε = 10⁻³ for Antarctic to ε = 5 × 10⁻³ for Greenland ice sheets.

Applying now the normalizations (4.1), (4.3) and (4.4), the momentum-balance and mass-conservation equations (2.2)-(2.4) become

In order to solve the above system of equations, the deviatoric stresses need to be related to the velocity components and , which reduces the problem to solving the system of three equations with three unknown functions , and . This requires that the components of the Cauchy-Green strain tensor B and the strain-rate tensor D are expressed in the stretched coordinates (4.2). Accordingly, the deformation gradient components (3.1), become

and the related Cauchy-Green tensor components (3.3), are given by

Further, the components of the velocity gradient, defined by (3.1)₂, are

and the strain-rate tensor components are given by

From the above relations it follows that and , and hence, for the isotropic ice flow, governed by the law , the relations and hold. On the other hand, since the normalized pressure gradients are order unity, Equation (4.5) implies that , and hence it follows that and . The latter relations imply, in view of and being order unity, that . We anticipate that as the ice fabric evolves from its initial isotropic state, none of the different axial and shear viscosities change by more than one order of magnitude, so all the above relations still hold true for anisotropic ice. This prompts the re-scaling of viscosities by introducing a new isotropic viscosity by

This viscosity re-normalization corresponds equivalently to defining ε ² as a dimensionless viscosity (Reference Morland and JohnsonMorland and Johnson, 1980; Reference MorlandMorland, 1984), and is that adopted by Reference HutterHutter (1983) and then by Reference Mangeney and CalifanoMangeney and Califano (1998).

Now, with the relations (4.8)-(4.12), the constitutive law (3.6), after its normalization by using (4.1), yields

where

In the isotropic state, when B = I, that is b1, = b₂ = b₃ = 1 and R = 3, in view of the identities (3.5)₂, (3.5)₄ and the normalization (3.8), we have

and the Relations (4.13)-(4.15) for the deviatoric stresses reduce to the isotropic viscous fluid relations

The shallow-ice approximation, not adopted here, is given by the leading order balances when all terms of order ε and smaller are neglected.

5. Numerical Model and Results

The problem defined by Equations (4.5)-(4.7) together with (4.13)-(4.15) is solved for the velocities , and the pressure by using the finite-element method. A weighted residual, or Galerkin, approach is adopted (Reference Zienkiewicz and TaylorZienkiewicz and Taylor, 1989), in which the problem equations are satisfied in an integral mean sense. The plane domain of interest is discretized by using triangular finite elements with six nodes (three vertices and three mid-side points). In each element the velocity field is approximated by using six nodal values of and and bi-quadratic shape functions, while the pressure field is approximated by three nodal values of (at vertices) and bi-linear shape functions. All surface integrals are evaluated numerically by applying Gauss-Legendre quadrature with seven sampling points within the element. Simultaneously with solving the FEM equations for , and Equation (3.2) describing the evolution of the deformation gradient F has to be solved to determine the fabric associated with the deformation field. For the plane problem considered here it is equivalent to solving a system of four first-order differential equations for the deformation gradient components F₁₁, F₁₃₎ F₃₁ and F33 at each node of the discrete system. We have assumed bi-linear variation of the Fij components over the element and hence the above system of differential equations is solved at each vertex nodal point. In order to find a steady solution of the problem discussed, an iterative procedure is followed. At the beginning of calculations it is assumed that ice is isotropic, that is F = I. For this initial isotropic state the velocities are determined by solving the FEM equations, and they are used to calculate from (3.2) the rate at which the deformation gradient changes. By imposing a restriction on maximum variation of the deformation gradient over one iteration, a maximum time-step At (which here has little physical meaning since the steady problem is analysed) is evaluated for solving (3.2) and updating the components F_ij. In computations, the above restriction has been defined in terms of the maximum relative change in the Euclidean norm of a vector containing the components of F at all discrete nodes, and it has been required that this maximum change over one iteration does not exceed 10–3. The current deformation-gradient components are used to find the values and directions of the principal stretches λr (r = 1, 3 ; λ2 = 1), which then define, through (3.5), the structure tensors M (r) and the values of the response functions f(b_r) and G(R that is a new fabric. For this new fabric, the FEM equations are solved to provide new values of the velocities and pressures, which are then used again to update the fabric, and the whole iteration process described above is continued until steady flow is reached. In our simulations, the calculations have been terminated when the relative change between two consecutive steps of the Euclidean norm of a vector containing all nodal velocities is smaller than 10–5. Usually, depending on the model parameters defining the strength of anisotropy, the stationary solution has been reached for .

A crucial part of a numerical model for the ice-sheet-flow analysis is the solution of Stokes’ Equations (4.5)- (4.7), since the incompressibility constraint (4.7) may give rise to numerical difficulties which are due to the absence of a pressure term in (4.7). Hence, a common method of treating this problem consists in adding an artificial pressure term to the continuity equation. There are a number of possible ways to do this. In the pseudo-compressibility method, formulated by Reference ChorinChorin (1967), a compressibility term in the form of a pressure time derivative is added to Relation (4.7), and when a steady state is reached, this artificial term vanishes. In an alternative approach, so-called pressure-correction method (Reference HirschHirsch, 1992), which can be applied to both stationary and time-dependent flows, an iterative procedure between the velocity and pressure fields is applied. This results in solving the momentum equations together with a Poisson equation for the pressure correction, obtained by taking the divergence of the momentum equations. The performance of this scheme depends very much on the accuracy of solving the Poisson equation. The pressure-correction method is closely related to a fractional-step method, also called projection method, which, with some modification, has recently been applied by Reference Mangeney, Califano and CastelnauMangeney and others (1996) to the steady ice-sheet-flow problem. Yet another approach consists in the adaptation of a method used for solving the problems encountered in incompressible elasticity and plasticity (Reference Zienkiewicz and TaylorZienkiewicz and Taylor, 1989). In this method, a term involving the pressure, times some small parameter, is subtracted from both sides of the incompressibility relation (4.7) and the resulting system of the momentum and modified continuity equations is solved by an iteration procedure, in which the initially adopted pressure is gradually updated until the stationary solution is found. In some ways, the latter method can be regarded as an improvement on the usual penalty method, in which the pressure times some penalty parameter is simply inserted in the continuity relation. Finally, it is also possible to solve the FEM equations resulting from the incompressible Stokes’ equations by using a special solver for semi-definite systems of equations (personal communication from R. C. A. Hindmarsh, 1998). In our model we have tested two methods in order to compare their numerical performance. The first method has been the pseudo-compressibility method, in which the system of first-order differential equations is solved. As an initial solution the velocities and pressures given by the SIA for isotropic ice have been used. The other method has been a combination of the penalty method, with a very small penalty parameter, and the application of a special solver for the semi-definite matrices. The small penalty parameter has been introduced since the available solver has still caused some numerical problems (oscillating solutions appeared) in the vicinity of ice-sheet margins with a steep, free-surface elevation, for instance for Vialov’s profiles, and at places with a rapid change in the finite-element mesh size. The comparison of the results given by the two methods has shown that the times, needed by both schemes to reach a stationary flow do not differ by more than 10 per cent, therefore the second method has been preferred in computations as being much cheaper (algebraic instead of differential equations are solved).

The numerical results presented below have been obtained for a simple geometry defined by a parabolic, free-surface elevation

where H is the ice thickness at the divide, and L is the lateral extent of the glacier (see Fig. 1). By adopting the latter geometrical dimensions as the typical scales introduced in section 4, that is by assuming H = h* and L = L*, Equation (5.1) can be written as

where is the normalized free-surface elevation, and . The simulations have been performed for the aspect ratio ε = H/L = l0–2. The ice anisotropy has been modelled by the response function (3.16) with m = 2, and we have adopted the constant isotropic dimensionless viscosity ; that is, for the needs of this study, we assume that the isotropic viscosity is independent of the strain rate and temperature.

In Figure 4 we show the variation of the normalized horizontal and vertical velocities at the free surface, and as well as the accumulation rate required to maintain the fixed geometry of the ice sheet. The results obtained for different combinations of the parameters A and S, describing the limit viscous response of ice in compression and shear, have been compared in order to illustrate the effect of these parameters on ice-sheet flow. We note immediately that the influence on the global flow of the parameter A, defining (for A > 1) the hardening of ice in uni-axial compression (recall Fig. 2), is very limited. Comparison of the results obtained for A = 3 and S = 0.2 (dashed lines) with the results for A = 10 and S = 0.2 (dotted lines) shows that the differences between the free-surface velocities are small and confined to the region of positive accumulation. Another feature worth noting is that for a large value of A meaning a high resistance of the anisotropic ice to axial stresses, the accumulation rate, q, necessary to maintain the fixed free-surface profile may be smaller than that for the isotropic case (though this applies only to a small region near the divide). A general conclusion for practical applications to be drawn at this point is that the significance of the parameter A (the reciprocal of the enhancement factor for compression), in terms of the global flow of the ice sheet, is very small. Hence, since there is still some disagreement among glaciologists on the proper magnitudes of A for polar ice, mainly due to the lack of realistic experimental data, our numerical results show that even the use of precisely determined enhancement factors for compression will have a negligibly small effect on the global flow-rate estimation for large ice sheets. On the other hand, as follows from Figure 4, the significance of the limit-viscosity factor S (being the reciprocal of the shear-enhancement factor) is crucial. It is clearly seen that the decrease of S (the increase of the enhancement factor) considerably accelerates the ice-sheet flow. Comparison of the results obtained for A = 3 and S = 0.2 (dashed lines) with those for A = 3 and S = 0.4 (dashed-dotted lines) indicates that both horizontal and vertical velocities, practically over the whole extent of the ice sheet, are about twice as high for S = 0.2 as for S = 0.4. This means that the global flow of ice is approximately proportional to the shear-enhancement factor 1/S. This feature is confirmed by Figure 5, in which we illustrate how the free-surface horizontal velocity, obtained for anisotropic ice (S ≤ 1) increases with respect to the corresponding velocity for the isotropic ice (S = 1); the results have been obtained for A = 3. The horizontal-velocity ratio is plotted against the enhancement factor 1/S for different locations We observe that at locations distant from the ice divide (curves for x/L = 0.4 and 0.6) the increase in the flow velocity is almost exactly proportional to the shear-enhancement factor. Only in the vicinity of the ice divide is this increase slightly smaller (curves for x/L = 0.1 and 0.2), which is obviously due to relatively high levels of the axial deviatoric stresses here, as opposed to the region far from the divide, where the flow is dominated by the shear stresses. To some extent, the conclusions drawn here can be related to the particular geometry of the adopted ice-sheet profile, therefore they should be further confirmed by the results of simulations carried out for other geometries, as well as by results obtained for the problem in which the free-surface elevation is calculated for a prescribed accumulation rate.

Fig. 4. (a) Free-surface horizontal velocities, (b) vertical velocities, and (c) the accumulation rates, for different combinations of the limit-viscosity factors A and S. The isotropic-^ results correspond to A = S=1 (solid lines).

Fig. 5. Variation of the free-surface horizontal velocity with increasing shear-enhancement factor 1/S at different locations x/L.

In Figure 6 we present depth profiles of the horizontal and vertical velocities at different locations , calculated for anisotropic ice (A = 3, S = 0.2) and isotropic ice (A = S = 1). The velocities are normalized by the respective free-surface values, that is the ratios and are plotted against the normalized elevation . We note that the horizontal-velocity profiles (Fig. 6a) calculated for the anisotropic ice (lines) vary with , contrary to the isotropic ice flow for which the normalized profile (indicated by circles) is common for all locations. This velocity pattern for isotropic ice differs slightly from those predicted by the models in which Glen’s flow law with the power factor n = 3 is used (Reference Dahl-JensenDahl-Jensen, 1989; Reference HvidbergHvidberg, 1996) due to the assumption of the constant isotropic viscosity made here. The most significant differences between the horizontal velocities obtained for anisotropic and isotropic ice occur in the near-divide region, while for x/L > 0.6 the profiles for both cases coincide. As regards the vertical velocity profiles (Fig. 6b) a more complicated pattern is observed, and the profiles for both anisotropic and isotropic ice vary with x/L; for the sake of clarity only the profiles for x/L = 0 (circles) and x/L = 0.6 (triangles) are plotted for the isotropic ice. The largest discrepancies between the velocities for anisotropic and isotropic ice take place, as for the horizontal velocities, near the ice divide, and again the corresponding profiles for anisotropic and isotropic ice coincide at locations distant from the divide.

Fig. 6. (a) Depth profiles of the normalized horizontal and(b)vertical velocities at different locations x/L for the anisotropic ice (lines) and the isotropic ice (symbols).

Finally, in Figure 7 we illustrate the depth variation of the shear and axial deviatoric stresses at different locations x/L, and again the results obtained for anisotropic (A = 3, S = 0.2) and isotropic ice are compared. Since the stresses and are of different orders, O(e) and O(e2) respectively, they have been re-normalized to become order unity, and hence the values and are plotted against the normalized elevation z/h. We observe in Figure 7a that the shear stresses are not affected by the anisotropy of ice, since the results for both cases are practically the same (lines for the anisotropic and symbols for the isotropic ice), which agrees with the shallow-ice approximation result that the direct shear stress is independent of the material properties. On the other hand, as seen in Figure 7b, the ice anisotropy dramatically changes the axial deviatoric stresses in a region close to the ice divide. An interesting result, which can be, however, a consequence of the particular ice-sheet profile adopted here, is that the maximum axial deviatoric stresses occur not at the ice divide, which is the case for isotropic ice, but at some distance (x/L ∼ 0.2) from the divide.

Fig. 7. (a) Depth profiles of the normalized shear stresses and (b) the longitudinal deviatoric stresses axz/ε at different locations x/L for the anisotropic ice (lines) and the isotropic ice (symbols).

Conclusions

In the paper we have solved the problem of a plane, steady-state, isothermal flow of a large polar ice sheet with a prescribed free-surface elevation, for which the accumulation rates required to maintain the fixed geometry are calculated. The strain-induced anisotropy of ice has been modelled by applying an orthotropic constitutive law, in which the viscous response of ice depends on current strain rates and deformations. The results of numerical simulations carried out for different sets of rheological parameters A and S, describing the limit viscous behaviour of ice in compression and shear, have shown the following:

1. (1) The global ice-sheet flow is little sensitive to the magnitude of the enhancement factor for compression and the influence of this parameter on the flow is confined to the region in the vicinity of the ice divide. For large values of the parameter A (reflecting significant hardening of ice in axial compression) the accumulation rates predicted for anisotropic ice flow may be smaller in the near-divide region than the rates obtained for isotropic ice;

2. (2) The ice-sheet flow is very sensitive to the magnitude of the enhancement factor for shear. It has been found that in most of the ice sheet (apart from a small domain adjacent to the divide) the flow rate is practically proportional to the value of this factor;

3. (3) In the near-divide region the normalized depth profiles of the horizontal and vertical velocities evaluated for anisotropic ice differ from those for isotropic ice, whereas in a region far from the ice divide the velocity profiles obtained for both anisotropic and isotropic cases coincide;

4. (4) The anisotropy of ice practically does not the affect the shear stresses, but significantly influences the distribution of longitudinal deviatoric stresses in the vicinity of the ice divide.

Since the above results have been obtained for a particular, idealized ice-sheet profile, further simulations involving a more realistic situation, in which a steady free-surface elevation is calculated for given accumulation rates, should be conducted before the general conclusions on the effect of the ice anisotropy on ice-sheet flow have been drawn. Nevertheless, the results obtained in this work have shown that anisotropy significantly affects the flow of ice in the near-divide region, where the deviatoric axial and shear stresses are of comparable magnitudes, while in the region where axial stresses are small and flow is dominated by shear stresses, anisotropy plays practically no role and isotropic ice solutions can be applied.