1. Introduction

Ice-core samples taken from depth in an ice sheet reveal strong fabrics, shown by significant alignment of initially randomly distributed c axes of individual crystals, and consequent substantial differences in shear viscosities in different planes. A macroscopic viscous law for the shear stress proposed by Reference Morland and StaroszczykMorland and Staroszczyk (1998) was motivated by a simple picture of the rotational recrystallization in which individual crystal glide planes, material planes, are rotated towards planes normal to principal axes of compression, and away from planes normal to principal axes of extension, described as the principal stretch planes. The instantaneous viscous response at each stage of the deformation has reflexional symmetry in these planes;that is, the instantaneous viscous response is orthotropic with respect to the current principal stretch planes. The directional strengths of the response depend on the current deformation. The orthotropic viscous law is then a frame- indifferent relation between deviatoric stress, strain rate, deformation and the three structure tensors defined by the outer products of the three orthogonal unit vectors along the principal stretch axes. Subsequently, a set of equalities and inequalities on the instantaneous directional viscosities was derived from the rotation concepts by Reference Staroszczyk and MorlandStaroszczyk and Morland (2000), designated by one author (L.W.M.) as the ‘Staroszczyk inequalities’, which showed that at least two tensor generators in the orthotropic law were required. Further refinements (Reference Staroszczyk and MorlandStaroszczyk and Morland, 2001; Reference Morland and StaroszczykMorland and Staroszczyk, 2003) led to a model law with two tensor generators, each with a fabric response coefficient function, but with an equality from the Staroszczyk inequalities providing an explicit relation between the two functions, so only one fabric response function need be prescribed. Two such model functions were constructed by correlation with complete (idealized) monotonic uniaxial compression and shearing responses which satisfy observed limit viscosity enhancement factors associated with cold and warm ice. The initial isotropic response enters as a multiplying factor depending on a strain-rate invariant and incorporating a temperature-dependent rate factor.

While this fabric-evolution model ignores such effects as grain growth and temperature variation during the evolution, it is adopted here to illustrate, for the first time, the possible significance of fabric effects included in the determination of an ice-sheet flow, with the cold- and warm-ice models reflecting contrasting influences. We have investigated steady radially symmetric flow to restrict the system of partial differential equations to two independent variables, horizontal and vertical spatial coordinates, which allow accurate numerical solution. This incorporates a realistic lateral spreading which is lost in the simpler plane flow equations. The adopted temperature field is prescribed in terms of surface elevation and depth below surface, such that there is uniform heat flux into the base and a constant rate of temperature decrease with elevation at the surface. The applied elevation-dependent surface accumulation/ ablation increases exponentially with elevation from large ablation at the margin to smaller accumulation at the divide. The prescribed basal sliding law is a linear relation between the basal tangential traction and velocity, with coefficient proportional to the normal traction which eliminates a surface slope singularity at the margin. The theory allows a prescribed basal melting/refreezing distribution, but this is assumed zero in the illustrations.

We also consider bed topography of modest slope relative to the horizontal mean bed plane. Reference MorlandMorland (2000) derived the enhanced reduced model for steady plane isothermal linearly viscous flow over beds of moderate slope, S, which are the leading order balances in asymptotic expansions of the full system of equations and boundary conditions, depending on the very small surface slope, є, and on the allowed larger topography slope, δ (є ≪ δ ≪ 1), and on s, the topography span relative to the sheet thickness, where 1 ≲ s ≲ 1/δ. The error magnitude is of order δ/s. The expansions were generalized for prescribed temperature and non-linear viscous flow by Reference Morland, Straughan, Greve, Ehrentraut and WangMorland (2001), obtaining the same error magnitude. We now derive the analogous leading order balances, an enhanced reduced model, with the above anisotropic viscous law incorporating fabric evolution. The error magnitude is now of order δ, which is greater than before when s ≫ 1; hence the description ‘modest slope’. Illustrations, and comparisons with the isotropic theory, are presented for both the cold- and warm-ice models, for a flat bed, a bed with a single hump of maximum slope δ = 0.01 and a bed with a single hollow of maximum slope δ = 0.01.

The numerical procedure is a non-trivial iteration between the hyperbolic partial differential equations governing the fabric evolution for a given velocity field, starting with the isotropic flow solution, and the solution of the two-point boundary-value problem given by the enhanced reduced model for the flow with given fabric. The final solution criterion is that the relative changes of margin radius, divide height and mean thickness are all less than 10^-5. Full details of the numerical methods are presented in section 7.

2. Fabric Evolution Anisotropic Viscous Law

3. Axially Symmetric Flow

Steady axisymmetric flow of a large ice sheet is described by the radial and vertical velocity fields, u(r,z) and w(r,z), respectively, where r is a horizontal radial coordinate and z is a vertically upward distance coordinate of cylindrical polar coordinates (r, θ, z), and all the flow variables are independent of the polar angle, θ. The corresponding nonzero physical components of the strain rate are

subject to the assumption of incompressibility, equivalent to the mass-balance condition,

which is satisfied identically by expressing the velocity components in terms of a stream function, ω(r,z),

The horizontal radial and vertical momentum balances for the very slow flow in which acceleration terms are negligible are, respectively,

where ρ = 918kgm^–3 is the constant ice density, g =9.81 ms^–2 is the constant gravity acceleration, and the circumferential momentum balance is automatically satisfied.

It is supposed that the atmospheric pressure at the surface is uniform, and stress is measured relative to this pressure. Then the ice-sheet surface z = h(r) is traction-free, and is subjected to a normal net ice flux q_n , which is positive over accumulation zones and negative over ablation zones. It is expected that accumulation will occur over a central zone at higher elevations, and ablation below a snowline. In general, q_n will depend on elevation h and location r. Let the surface have unit tangent and outward normal vectors s and n in a righthand sense, then

The zero traction is conveniently expressed in terms of vanishing normal and tangential tractions t _n = n ×σ× n and t _s = s ×σ× n in the Orz plane,

The kinematic condition prescribing the surface accumulation flux becomes

where v _n is the normal ice velocity and is the equivalent vertical flux which can be identified with equivalent vertical snowfall. In general, will depend on surface elevation and location, but the later illustrations assume dependence on elevation only - dependence on location only is a simpler problem, but artificial. Equation (32) is the additional boundary condition required to determine the unknown surface h(r). A prescribed temperature field incorporates a prescribed surface temperature.

The bed z = f(r) is prescribed. The vectors s and n are defined by

Normal and tangential tractions t _n and t _s are given by

and normal and tangential velocities v _n and v _s by

The kinematic condition prescribing the normal basal melt flux, bn, negative if freezing occurs, is

A linear sliding relation between the tangential traction and velocity is assumed, given by

where γ is a friction coefficient depending in general on the overburden pressure –t _n. The proportionality of t _s to t _n ensures bounded surface slope at the margin where the ice thickness, and overburden pressure, approach zero (Reference Morland and JohnsonMorland and Johnson, 1980). The case of no-slip gives rise to unbounded surface slope at a margin at which there is ablation in the reduced model (Reference Morland and JohnsonMorland and Johnson, 1980; Reference MorlandMorland, 1997), which violates the expansion scheme. Figure 1 shows an ice-sheet profile with the above conditions displayed, and gives notation for physical coordinates, velocity components and the later dimensionless variables.

Fig. 1. Ice-sheet geometry.

4. Enhanced Reduced Model

We now introduce typical magnitudes for the physical variables so that normalized dimensionless variables can be defined. Let d ₀ be a thickness magnitude and q ₀ a vertical velocity magnitude, defined by a typical surface accumulation, then pgd ₀ is an overburden pressure magnitude, d₀/q₀ is a flow time magnitude and q ₀/d ₀ is a vertical strain-rate magnitude. Typical values are adopted here,

then

Define

then Z, H, P, Q, B are order unity variables, F is of order of its amplitude, ĸ, which may be order unity or less, and W is yet to be estimated. The viscous law (14) now becomes

where is defined by

and with the chosen parameters,

We formally treat the product as an order unity quantity, since it was noted by Reference Morland and SmithMorland and Smith (1984) that as a(T) becomes small in cold upper regions of an ice sheet, so does the strain rate. є² is the magnitude of the dimensionless viscosity in law (42).

In the reduced model, leading order equations are developed by introducing the coordinate and variable stretchings

where R, U and Ω are order unity, and

We are allowing a bed topography with modest slope δ, and with amplitude k and span ĸ relative to the thickness d ₀, where

as treated by Reference Morland, Straughan, Greve, Ehrentraut and WangMorland (2000, 2001) and Reference Cliffe and MorlandCliffe and Morland (2002). Due to the bed topography, variables change with r over a dimensionless length scale s near the bed, so that derivatives leave the order unchanged. Thus F’(R) is order ĸ/(єs) = S/є ≫ 1, and it follows from the basal melt condition (36) that W is also that order. All R derivatives of order unity variables will be of order (єs)^–1, with higher derivatives multiplied successively by this factor. We expect the surface to vary on the R scale, i.e. the scaled surface slope H′(R) = Γ(R) to be order unity, but will allow it to be as great as the scaled bed slope F′(R) = ß = O(δ/є) in the later asymptotic expansions. For valid expansions at the margin, it must be assumed that the bed slope in a margin region must not be greater than є in magnitude, and hence F’(R) is order unity there. The analysis for the isotropic viscous law shows that the error of the enhanced reduced model is of order δ/s. We will now show that with the anisotropic viscous law (42) the error becomes of order S, which is greater when s ≫ 1.

First we must estimate the magnitudes of the dimensionless scaled stress components and their derivatives which appear in the momentum balances and boundary conditions. Note that all components of never exceed unity in magnitude because of property (17) and definition (13). Immediately from (46), but and are greater than order unity, so by the incompressibility condition (27) their leading order terms must be equal and opposite, namely of order (єs)^–1. This implies that where R derivatives need the factor (єs)^–1 near the bed, Z derivatives are increased by the factor ĸ¯¹ to reflect the topography effects, and are therefore of order є/S smaller than the corresponding R derivative. Thus, from (46),

The relation applies throughout the depth, even when Z derivatives away from the bed leave magnitudes unchanged. Further, the order of the R derivative of W in is a factor δє or δ², near or away from the bed, respectively, of the Z derivative of U.

All components of each M ^(l) are order unity. Define

Note that e ⁽²⁾ = (0,1,0), so that for l = 1,3, = 0, and

The non-zero component of A ⁽²⁾ is

and for l = 1,3 the non-zero components of A ^(l) are

all dominated by with relative error of order δ by (49).

The non-zero components of C are

where the product does not exceed unity in magnitude, so again all components are dominated by with relative error of order δ.

The usual scaled dimensionless deviatoric stress, is defined by

then from (42), (49), (50) and (52–54),

all with relative error δ at most, where

The expressions for c_rz and c_zz are order unity with relative error δ. Then the shear stress invariant, J, defined in (2), becomes

with the chosen parameters (4), (38) and (39). With the isotropic viscous law, c_rz ≡ 1, c_zz ≡ 0. By (49), is of order (δєs)^–1, and hence the deviatoric stress components given by (42) are of order є/(δs), which is much smaller than unity in view of the strong inequalities (48);in fact, not greater than δ if we assert s ≳ єδ^–2, which is a reasonable restriction, consistent with (48), and is assumed here. Now all the components of are order δ/є.

From (30) and (33) it follows that ∆_h and ∆_f are both unity, with maximum error δ², so from (31) and (32) the surface conditions become, to lead order, neglecting terms of order δ,

Similarly, neglecting order δ terms, on the bed, from (34) and (35),

so that the basal boundary conditions (from (36) and (37)) become

and we suppose Λ is order unity or greater. The momentum

balances (29) are exactly

Since the components of are at most of order δ/є, and recalling the deduction above Equation (49) that R derivatives are not greater than a factor δ/є of Z derivatives, the leading order balances of (63) and (64), neglecting terms of order δ in the former, and order δ and δ² in the latter, become

Equations (60), (62) and (65) are those of the enhanced reduced model, but the fabric changes the deviatoric stress expressions (56) and the invariant J defined by (59) through the factors c _rz, c _zz and c defined by (57) and (58).

5. Flow Solution

Here we decouple the radial flow boundary-value problem defined by (55–59), (60), (62) and (65) from the fabric evolution, by assuming that the fabric through the ice-sheet domain is known;that is, c_rz(r,z) defined by (57), c_zz(r,z) defined by (58) and c(r,z) defined by (59) are given functions. Subsequently we solve the fabric-evolution equations assuming that the velocity field is known, and then construct the coupled flow by an iteration procedure. is defined by (43) with a(T) given by (7) and ψ(J) by (6). Solutions of the momentum balances (65) subject to the surface conditions (60), and the invariant J defined by (59), have the simple expressions

Denote basal values on Z = F(R) by a subscript ‘b’, then by (66) and (62),

where ∆(R) is the ice thickness at R. From (57), (60), (67) and (47),

Integrating twice with respect to Z,

Now the surface accumulation condition (60) and basal melt condition (62) with the relations (47) become

where subscript ‘s’ denotes evaluation on the surface Z = H(R). Differencing the relations in (72), and evaluating Ω_s by (70), gives

where U _b is given by (67) in terms of Γ(R) and Λ(D), and

and g and g_2s depend on Γ, hence dg_2s/dR includes Γ’(R). It follows that (73) is a second-order differential equation for H(R).

This differential equation applies over an unknown span, the margin radius R _M, which must be determined as part of the flow solution. The subscript ‘M’ will denote evaluation at the margin. In section 7 the explicit second-order equation for H(R) will be transformed into an equivalent system of first-order equations, and an asymptotic limit for Γ’(R) as R → 0, since Γ’(R) is indeterminate at the divide, will be given (see (90)). As R → R_M , so Δ → 0, and recalling the assumption that ß is order unity near the margin, an asymptotic analysis (Reference MorlandMorland, 1997) determines both Γ_M and ΓM in terms of R_M and prescribed functions evaluated at R _M; thus

Note that Γ_M < ß _M for positive thickness inside the margin, so the Γ′_м factor is non-zero, and and denote derivatives with respect to R, and not their given arguments. While the second-order equation could be integrated from the margin for a sequence of trial R _м, with the margin slope given by (76), until the R _м is found such that the surface slope becomes zero at R = 0, the indeterminancy there in Γ′ is not given by the finite asymptotic limit for each incorrect trial R_м, and the numerical procedure is unstable (Reference MorlandMorland, 1997). The stable approach (Reference Cliffe and MorlandCliffe and Morland, 2002) treats the problem as a two-point boundary-value problem with both R_м and the divide elevation H(0) = H _D as unknown parameters, where the subscript ‘D’ denotes evaluation at R = 0. The differential equation is integrated from both the margin where H = F and H′ = Γ_м, and from the divide, R = 0, where H = H _D and H′ = 0, for a sequence of trial R _м and H _D until both surface elevation H and slope Γ match at a chosen interior radius. The method is described further in section 7. Once the profile H(R) with R _м and H _D are determined, the velocity fields and derivatives can be determined to apply to the fabric-evolution equations.

6. Fabric Evolution

Now we assume that the domain and velocity fields and derivatives required in the evolution equations for B are known. Let be a dimensionless time with unit d ₀/q ₀, then the kinematic relation (9) gives the non-zero relations

The relations are retained in full since some components of B increase indefinitely as the stretches normal to the dominant compression increase indefinitely; we cannot set bounds. The incompressibility condition (11) is satisfied identically by using

and ignoring the evolution equation for B_zz . It is expected that B _rr ≥ B _θ ≥ 1, and B_zz ≤ 1, or close to this situation, and both B_rr and B_θθ will approach infinity, and B_zz will approach zero, as the ice is compressed vertically as it approaches the bed, while B_rz , zero in the undeformed state, may be positive or negative as the ice flows through the sheet, but not greater in magnitude than B _rr. Numerical solution is most sensibly conducted in order unity variables, and after some trials we found that introducing the following normalized variables yielded satisfactory evolution equations,

The three evolution equations for ζ₁, ζ₂ and ζ₃, by (78), are

where, in steady flow, the material time derivative reduces to the convective terms

for all ζ, and by the incompressibility relation (79),

7. Numerical Method

We apply the following iteration procedure. First solve the flow problem for the profile and velocity field, and the strain-rate components and velocity derivative ∂W/∂R which are required in the fabric-evolution equations (81), for a given fabric distribution defined by functions c_rz (r,z) and c _zz r,z). Then solve the fabric-evolution problem with this flow field to determine new fabric functions c_rz (r,z) and czz(r,z), which are then used to solve a new flow problem, and continue this sequence. We start by solving the flow problem with isotropic ice for which c_rz (r,z) ≡ 1 and

c_zz(r,z) ≡ 0, and continue the iteration until the final flow solution profile is sufficiently close to the previous one. To move between the two problems – a second-order ordinary differential equation with two end-point boundary conditions, on an unknown span R _M changing at each iteration, and the hyperbolic partial differential equations governing the fabric evolution – we map each ice-sheet domain in the sequence on to the same fixed domain, a unit square. Introducing new variables (t,y) by the transformations

the surface Z = H(R) becomes y = 0, the bed Z = F(R) becomes y = 1, and the divide R = 0 becomes t = 0. The margin R = R _M, where the thickness Δ(R_M) = 0, is expanded to 0 < y < 1 at t = 1, but no depth integrals in the flow equations are required when Δ(R) = 0. It is supposed that there is a single snowline in the later illustrations, at t_E, so that there is surface accumulation for 0 ≤ t < t _E and surface melt for t _E < t < 1; that is

The mapped domain and a typical ice particle path are shown in Figure 2.

Fig. 2. Typical ice particle path in the (t,y) domain.

With the change of variables

the second-order differential equation defined by (73–75), with definition (71) for g₂(R,Z), and unknown parameters RM and H _D, can be expressed as four first-order differential equations for R _M and H _D:

subject to the boundary conditions

where Γ_M is given by (76). The complexity lies in the surface slope derivative given by

where

and

The integrand functions of (R,Z) must be transformed to (t,y) dependence before integration, and in particular the partial derivative ∂c/∂R at constant Z must be expressed in terms of partial derivatives with respect to t and γ, then the latter is integrated out by parts.

The expression (89) for G is indeterminate as t → 0 at the divide, and as t → 1 when at the margin. The latter asymptotic limit is simply , where Γ_M is given by (76). The former asymptotic limit is

The velocity fields U and W are then given by (47), and the strain-rate components, required by the fabric-evolution equations (81), are given by (46) and (56) with Ω given by (70). These involve integrals from y- to 1 at t_i for each mesh point (t _i,y_j ). In addition, (81) requires the velocity derivative ∂W/∂R at constant Z, which is determined numerically by cubic expansions of W in the t and y directions using values at (t_i ,y_j ) and three adjacent points.

A uniform mesh is placed on the square with n intervals of length h _t in the t direction, and m intervals of length h _y in the y direction, defining mesh points (t_i,y_j ), (i = 0, ...,n; j= 0, ...,m). Given that the fabric-evolution solution has determined c_rz, c_zz and hence c, at each mesh point (t_i ,y_j ), the required t derivatives are obtained by cubic interpolation based on four adjacent mesh points in the t direction. For the flow problem, the depth integrations are performed using a trapezoidal approximation, with error and the differential equations (88) are solved using a second-order Runge–Kutta scheme with error Newton’s root-finding method is used to determine the parameters R _M and H _D which match the surface elevation and slope at the interior point t = 0.5. Using double precision, the matching error for all the later examples was initially of order 10^–8, retaining eight decimal places in R _M and H _D, but the final solutions were obtained retaining only six decimal places. The iteration between flow and fabric evolution was continued until a measure of the change from the previous flow solution defined by

was acceptably small, where [ ] denotes the jump from the previous solution in the iteration sequence. In all the later examples er did not exceed order 10^–5, and the rate of decrease in the later stages suggested that a further iteration would have yielded a much smaller value, so this final solution can be expected to change less than 10^–5 with a further iteration.

Figure 2 shows a typical particle path in the (t, y) domain, starting from an accumulation point on the surface (t _s,0), where 0 < t _s < t _E, and returning to the surface at an ablation point (t _r,0) where t _E < t _r < 1. Let for all ζ, and let s denote arc length on a particle path in the (t, y) domain defined by , for which

The fabric-evolution equations (81), with the relation (82), now become

(92–94) are five ordinary differential equations for subject to the initial surface conditions at s = 0 on y = 0, t = t _s where 0 ≤ t _s ≤ t _E,

Values of to and in turn the values of c_rz , c_zz , c and their derivatives, must be determined at each mesh point (t _i,y _j), where they are required for the subsequent flow solution.

The divide t = 0 is a particle path with s = y, but as t → 0, so does U, and and hence y _d → ∞ when W ≠ 0, but

and ν is bounded. On t = 0, (positive), and we suppose the bed slope does not exceed order є near the divide, so the derivative ∂W/∂R is order unity there. Then, neglecting the term of order є, the three non-trivial equations for to on the divide become

which can be integrated numerically. However, if W→ 0 as well, then y _d is indeterminate, as is the case in our later examples which assume zero basal melting, and hence W = 0 at the base of the divide. However, the ν limit given in (95) is still the correct asymptotic behaviour, the other terms of y_d being smaller, and since is positive, then is negative by incompressibility, and hence W behaves like –(1 – y) and ν like (1 – y)^–1, positive, as y → 1. Immediately so behaves like a positive power of (1 – y), and like a larger positive power of (1 – y), so B _rr → ∞, B _rz → 0 and B _zz → 0, by incompressibility. On the divide, the unit vectors e ^(l) are along the coordinate axes so b ₁ = B _rr, b ₂ = B _θθ, b ₃ = B _zz, and hence, by (13), ξ₁ → ∞, ξ₃ → –∞, η ~ K ~ 2b₁ → ∞, and by (12) each is zero. Thus by (21), (57) and (58),

This asymptotic limit is approached smoothly by the numerical integration to the last point before the bed y = 1. The bed is also a particle path when B = 0, but direct numerical integration in s = t from the margin is not possible because of the singular, but not explicit, behaviour of ν there.

The first approach was to integrate the five differential equations for to from the surface mesh points on (0,t_E) using a second-order Runge–Kutta method with integration interval h_t , until the closest point to (t _r,0) was reached. Values at each mesh point were then interpolated by using nine or ten closest path points in both rectangular and polar expansions. It was found that all such interpolations became erratic in some neighbourhoods of the bed, y = 1, which we believe is because the paths become parallel and bunch close together there.

An accurate method is to integrate the path equations (92) for and backwards from each mesh point (t _i,y_j ), to determine the corresponding start value t _s on y = 0. The final value is determined by using repeated tangent extrapolations over shorter intervals from the last point before the line y = 0 is crossed, stopping when the last interval is less than 10^–6. Then all five evolution equations (92) and (93) are integrated forward from (t _s,0) to (t _i,y _j). The high accuracy is confirmed by the accurate return to (t,•,yj). With the necessarily fine mesh using second-order integration routines, computing time for the reverse and return integrations for the large number of mesh points is extremely large (many hours on an Apple G4), and must be repeated many times for the iteration between evolution and flow. An alternative strategy was therefore devised. The backward integration from each (t_i,y_i ) was performed to determine the point (t₀,y₀ ) where the path first crossed a previous mesh line: t = t_i–1 or y = y_j–1 or y = y_j+1; the four possibilities are shown in Figure 2 with typical mesh points denoted by a cross and the previous mesh line intersections by circled points. The above tangent extrapolation method was again used to determine accurate (t₀,y₀). The values of at (t ₀,y₀) are determined by quadratic interpolation using values at three appropriate adjacent mesh points, then the five equations are integrated forward to (t_i,y_j). The different possible configurations for the location of the point (t ₀,y ₀), which require known values at different sets of adjacent mesh points, complicates the procedure for covering the domain to provide all required values for the interpolations, but this strategy has proved very fast and accurate. The accuracy is confirmed by applying the full reverse and return integrations to selected sets of mesh points: adjacent to the divide, adjacent to the bed, and in the ablation zone of the surface;in particular the bed neighbourhood could be expected to yield the worst errors. In all the later examples, it was only at a few places in the bed, towards the margin, that significant differences arose, up to a few per cent in small zones, but these should make little difference to the depth integrals in the flow equations. The solution for mesh points on the bed was obtained by integrating the five differential equations from a surface point (10^–7,0) for which the particle hugs the divide, bed and expanded margin line, and linearly interpolating between values at integration points t on either side of t_i. The final values at (1,0) are used for the complete line t = 1, 0 ≤ y ≤ 1, since this is the expanded margin.

8. Illustrations

A temperature distribution used by Reference MorlandMorland (1997) is adopted,

where H and Δ denote H(R) and Δ(R), respectively, which depends on surface elevation and depth below the surface. The corresponding thermal boundary conditions are

with the realistic properties that the surface temperature decreases at a rate of 0.8 K per 100 m rise, and that there is a uniform heat flux into the base equivalent to a temperature gradient of 0.5 K per 100 m, and the Laplacian of T, arising in the energy balance, is bounded.

The surface accumulation/ablation distribution adopted is an elevation-dependent example with large margin ablation, –Q _M, used by Reference MorlandMorland (1997),

with a decay height scale H = 0.25 corresponding to h = 500m. The equilibrium height (snowline) where Q = 0 is at H = 0.64 corresponding to h = 1280 m. Reference Cliffe and MorlandCliffe and Morland (2004) have investigated the dynamic stability of radial flow for linearly viscous ice with varying weighting of the accumulation dependence on elevation and location, and demonstrate in particular that dependence on elevation alone is unstable. However, direct solution of the steady flow equations with this function has always proved numerically stable, and it describes a simple, physically plausible distribution. The basal melting B ≡ 0, and a constant friction coefficient Λ ≡ 25, are adopted in all the examples presented.

We compare solutions for a flat bed, for a bed with a single symmetric hump of amplitude 0.4, corresponding to 800 m, maximum slope δ = 0.01, and slope span s = 40, corresponding to 80 km, centred at R _c = 0.4, corresponding to r _c = 479 km, and for a single symmetric basin which is the image of the above hump. All three cases are described by the smooth bed form

with

where the R span of the slope sections is e, and of the flat section is 2e, and is a monotonic polynomial

F(R) and its first three derivatives are continuous at R _c–2e, R_c–e, R_c + e and R_c + 2e, and F(R) has a maximum slope at χ = 0.5 and where

with the chosen parameters, corresponding to a maximum physical slope δ = 0.01 and a hump/basin spanning 174 km. The flat, hump and basin topographies are given, respectively, by κ = 0, 0.4 and –0.4.

The solution computations were performed on a mesh with n = m = 1000 so that h _t = h _y = 0.001. In all three cases, the fabric parameter c_zz remained very small, not exceeding a few per cent, so the parameter c in (59) remained close to unity. The fabric influence on the flow is therefore mainly through the fabric parameter c _rz, the viscosity factor arising in the relation (56), and defined by (57). We now present comparisons of the isotropic-ice, cold-ice, and warm-ice solutions for the three bed topographies. Tables 1–3 show H _D, R _M, R _E, (c_rz )_DB and (c_rz )_M, where R _E is the scaled radius where Q(H) is zero, and (c_rz )_DB denotes the value at the base of the divide, here R = 0, Z = 0, given by the asymptotic limit (97), for the flat, hump and basin beds. At the base of the divide, where the vertical compression is large, (c_rz)_DB has the same limit values for all three topographies, but this greatly exceeds unity for cold ice, and is significantly less than unity for warm ice, associated with the limit compression factor A, given in (20), exceeding unity for cold ice and less than unity for warm ice. The viscosity factor at the margin, (c_rz )_M, decreases from its unit isotropic value for the cold ice, and more dramatically for the warm ice, in the flat- and hump-bed cases, but in the basin-bed case it still decreases dramatically for warm ice, but is close to unity for cold ice. For all three bed forms, the margin span R_M increases from isotropic to cold to warm ice, while the divide height H_D increases slightly then decreases more for the same sequence, but, as expected, increases from basin to flat to hump beds.

Figures 3–5 show the bed and surface profiles for the isotropic, cold and warm ice, for the flat, hump and basin beds, respectively. These illustrate the distinct, and significant, effects of the two forms of fabric evolution treated, particularly on the span R _M. For a given R _M and fabric law, the divide height H _D is essentially governed by the steady-state assumption, that there is no net flux in or out of the sheet.

Fig. 3. Free surface profiles for isotropic, cold and warm ice; for flat bed.

Fig. 4. Free surface profiles for isotropic, cold and warm ice on the illustrated basin bed. illustrated hump bed.

Fig. 5. Free surface profiles for isotropic, cold and warm ice on the illustrated basin bed.

The viscosity factor c_rz is identically unity for the isotropic ice. Figure 6 shows the distribution of c_rz , first down the divide as a function of H – Z, then along the bed as a function of R, for cold and warm ice, in the flat-bed case. This gives a measure of how the fabric is affecting the shearing rate of the ice, down the symmetry axis and along the bed. Figure 7 shows the same for the hump bed and Figure 8 for the basin bed. There is a significant effect from the bed topographies. While c_rz for the warm ice is mainly decreasing along these two boundaries, consistent with the increase of span R _M from the isotropic span, the effect of decreasing c_rz along the bed, and in the sheet interior, for the cold ice, must offset the increase down the divide, since R _M is again greater than the isotropic span for each bed form, though less than the warm span.

Fig. 6. Distributions of the viscosity factor, c_rz , down the divide (as a function of H – Z) and along the bed (as a function of R) for cold and warm ice; for flat bed.

Fig. 7. Distributions of the viscosity factor, c_rz , down the divide (as a function of H — Z) and along the bed (as a function of R) for cold and warm ice; for flat bed.

Fig. 8. Distributions of the viscosity factor, c_rz , down the divide (as a function of H — Z) and along the bed (as a function of R) for cold and warm ice; for basin bed.

Figure 9 shows, for the flat bed, the distributions of the horizontal scaled velocity, U, with height at locations R = 0.2, R = 0.4, R = 0.6 for isotropic, cold and warm ice. Recall that the hump/basin topographies are centred at R _c = 0.4. The effects of fabric are increased with distance from the divide as expected, particularly near the bed. Figure 10 shows the same for the hump bed, and Figure 11 for the basin bed, with similar effects, and with increased changes near the bed due to the topographies.

Fig. 9. Horizontal velocity depth profiles at R = 0.2, R = 0.4 and R = 0.6, for isotropic, cold and warm ice; for flat bed.

Fig. 10. Horizontal velocity depth profiles at R = 0.2, R = 0.4 and R = 0.6, for isotropic, cold and warm ice; for hump bed.

Fig. 11. Horizontal velocity depth profiles at R = 0.2, R = 0.4 and R = 0.6, for isotropic, cold and warm ice; for basin bed.

9. Conclusions

An enhanced reduced model for steady radially symmetric flow of an ice sheet with fabric evolution, over bed topographies with a modest maximum slope, S, has been derived, where є ≪ δ ≪ 1 and є is the surface slope magnitude. It is described by the leading order balances of the mass and momentum equations, boundary conditions and the constitutive law for the ice governing the fabric (induced anisotropy) evolution, with errors of order δ.

Numerical solutions for a flat bed, a bed with a single symmetric hump and a bed with a single symmetric basin have been constructed, for two forms of the fabric evolution, described as cold ice and warm ice. Tables of the key parameters – divide height, margin radius, equilibrium-line radius and the fabric-induced viscosity factor – are presented for the isotropic, cold and warm ice for each of the bed forms. The radial span, margin radius, is seen to be strongly influenced by the induced fabric. Figures illustrate the distinct ice-sheet profiles, the distribution of the fabric- induced viscosity factor down the divide and along the bed, and the distributions with height of the horizontal velocity, for the different ice laws and bed topographies. All illustrate the significant effect of fabric evolution compared to the conventional isotropic ice model.

The adopted fabric-evolution model does not account for possible fabric weakening due to migration recrystallization occurring at high (near-melting) temperature or high strain rate, which may arise in thin zones near an ice-sheet bed. Since bed conditions have a major influence on the flow and profile, such weakening could make a significant difference to the solutions and conclusions. If, however, such zones are not extensive, but confined to limited areas of the bed, then the magnitudes of the effects derived here may still apply. To reach a firmer conclusion, a constitutive model which realistically describes such weakening effects must first be established and incorporated in the ice-sheet flow equations, and then solutions for simple flow configurations obtained.