Introduction

On the large timescales of ice-sheet flow, the ice is assumed to be incompressible and to obey a non-linearly viscous constitutive law for the shear response, neglecting the shorter-timescale viscoelastic effects. In all large-scale ice-sheet modelling to date, the behaviour is supposed to be that of a simple viscous incompressible fluid for which, at constant temperature, the deviatoric stress depends only on the strain rate. The pressure is a workless constraint, not given by any constitutive law, but determined by the momentum balances and boundary conditions. Such a viscous law, necessarily isotropic by material frame indifference, has a general quadratic representation, with alternative, but equivalent, stress and strain-rate formulations, as discussed by Reference MorlandMorland (1979). However, it is still common practice to ignore the quadratic term and adopt a simple relation proposed by Reference NyeNye (1953) in which the deviatoric stress is coaxial with the strain rate, and which depends on only one of the two stress (or strain-rate) invariants.

Reference GlenGlen (1958) (acknowledging F. Ursell) presented the quadratic viscous relation for the strain rate, and noted that Reference SteinemannSteinemann’s (1954) combined compression and shear data were inconsistent with the simple form, so that dependence on a second invariant, or the quadratic term, or both, is necessary. Standard single stress component tests, uniaxial compression or simple shear, are not sufficient to determine the general form, nor can they verify the validity of the simpler parallel form. The determination of two response functions with two invariant arguments requires biaxial or combined compression and shear tests. While combined compression and shear tests have been conducted (see, e.g., Reference Li, Jacka and BuddLi and others, 1996; Reference Warner, Jacka, Li, Budd, Hutter, Wang and BeerWarner and others, 1999) there has been no attempt to correlate data with other than the simple parallel form. Reference Morland and EarleMorland and Earle (1983) demonstrate that conventional triaxial stress tests (identical lateral axial stresses) provide no more information than uniaxial stress tests, and explore the domains of the two-invariants plane covered by uniaxial, biaxial and triaxial stress tests, but not that for combined compression and shear tests. The analyses presented here of confined and unconfined compression combined with shear show that both can check the consistency of the general quadratic isotropic form, independent of the actual response functions, assess the significance of the contribution made by a quadratic term and then determine the two response functions of two invariant arguments. Laterally confined or unconfined biaxial stress tests are also analyzed;the former is inadequate to determine a general relation, and the latter is more cumbersome than the combined compression and shear tests, and does not yield a check on the consistency of a viscous fluid relation assumption.

A general viscous relation that includes the quadratic term has yet to be incorporated into the flow equations for an ice sheet, but would certainly yield changes from the simple coaxial form. The commonly adopted reduced model equations are the leading-order balances of an asymptotic expansion in a small parameter arising from a dimensionless viscosity magnitude in the parallel form, which in turn defines the small surface slope magnitude or sheet aspect ratio (shallow-ice approximation). Allowing that the quadratic term is significant, two small factors arise in the corresponding dimensionless viscous relation, but can be related. It is now shown that the relative magnitudes of the stress components are changed from those with the parallel form, but that the same leading-order momentum balances, and their explicit depth integrations to yield the same expressions for the pressure and horizontal shear stresses, are obtained, though with different error magnitudes. However, the resulting expressions for the velocity gradients with depth cannot be explicitly integrated, so no expressions for the velocity fields are available to incorporate into the surface and bed kinematic conditions to obtain a differential equation for the surface profile. This same problem arises in the simple coaxial form if the response function depends on two invariants, which is the minimum generalization implied by Reference SteinemannSteinemann (1954). In any case, the reduced model applies only when the bed topography also has small slopes, which is not reality, and a quadratic term in the viscous relation will have a significant effect on the full flow equations. It seems crucial that the nature of the viscous relation, whether it is linear or quadratic and whether it depends on one or two invariants, is explored further to justify or improve current ice-sheet modelling.

The Viscous Relation

Let σ denote the Cauchy stress tensor, and the deviatoric stress tensor, then

(1)

where p is the mean pressure and l is the unit tensor, and let D denote the strain-rate tensor, the symmetric part of the velocity-gradient tensor. Then the most general frame-indifferent relation between and D can be expressed in two alternative, but equivalent, forms of the Rivlin–Ericksen representation between tensors with zero trace:

(2)

(3)

where I ₁, I ₂, I ₃ are the principal invariants of D and J ₁, J ₂, J ₃ are the principal invariants of , defined by (omitting the minus sign in the second invariant for convenience),

(4)

(5)

The vanishing of I ₁ is the incompressibility condition, and the vanishing of J ₁ is by the definition of the deviator, so the response functions, ϕ ₁, ϕ ₂, ψ ₁, ψ ₂, each depend on only two non-trivial invariants. While expansions (2) and (3) are equivalent, there is no explicit algebraic inversion. Note that these expansions for a simple viscous fluid are necessarily isotropic in all reference configurations, and cannot describe induced anisotropy associated with the fabric developed as the ice elements deform and crystal glide planes are reoriented. It is the stress formulation, (2), which is required for substitution in the momentum balances of a general ice-sheet flow. Note that ϕ ₂ = 0 implies, and is implied by, ψ = 0.

The dependence on temperature, T, is assumed to be given by a rate factor, a(T), by replacing D by an effective strain rate, , defined by

(6)

in all the above relations, where a(T) is an increasing function of T; that is, the actual strain rate at a given stress increases with temperature. As noted by Reference MorlandMorland (1979), this is the assumption of a thermorheologically simple response (Reference Morland and LeeMorland and Lee, 1960) in which the same processes occur, but on a timescale factored by a(T).

The pioneering experimental work of Reference GlenGlen (1952, Reference Glen1953, Reference Glen1955, Reference Glen1958) and Reference SteinemannSteinemann (1954) was used to construct power laws for a simpified form of (3), proposed by Reference NyeNye (1953), with an equivalent form of (2), namely

(7)

in which D is coaxial with , and there is only one response function, ψ₁, depending on only one invariant, J ₂, which is a measure of the shear stress squared, or one response function, ϕ ₁, depending on one invariant, I ₂, which is a measure of the strain rate squared. The tests were mainly on polycrystalline ice at constant temperature with randomly oriented crystals, in either unconfined compression or simple shear, so that only one function of one argument could be inferred. Comparisons of known datasets at the time were made by Reference Smith and MorlandSmith and Morland (1981), which demonstrated wide differences, and the form (7) was constructed from the Reference GlenGlen (1955) data, showing that a three-term fifth-order polynomial representation, with finite viscosity at zero stress, was a much closer correlation than a power law with infinite viscosity at zero stress. Reference LliboutryLliboutry (1969) and Reference Colbeck and EvansColbeck and Evans (1973) had constructed similar polynomial representations from different data, but without error estimates. Reference Smith and MorlandSmith and Morland (1981) also constructed exponential representations for the rate factor, a(T), over different temperature ranges from the data of Reference Mellor and TestaMellor and Testa (1969).

Combined Shear and Compression Analysis

I now analyze the stress and strain-rate patterns for shear combined with both confined and unconfined compression, and derive the various relations following from the general isotropic viscous fluid relations (2) and (3). In particular, universal relations independent of the two response functions are derived, from which data can confirm the consistency of an isotropic fluid assumption, or reject the assumption, and from which data can estimate the significance of the quadratic tensor terms in the general isotropic relations if isotropy is valid. The relations allow the two functions of two invariant arguments to be determined by these tests, or, if the quadratic term is small, one function of two invariants, or of one invariant if the simple form (7) proves to be adequate.

Longitudinally confined compression

For combined shear and compression with no longitudinal extension and free lateral extension, the strain-rate tensor and quadratic tensor in Equation (2) are

(8)

(9)

where – is the rate of axial compression, is the shear stain rate and the invariants are

(10)

The associated stress tensor, mean pressure and deviatoric stress tensor are

(11)

(12)

and the invariants are

(13)

(14)

While σ_z and τ are known applied stresses, the constraint stress, σ_x , necessary in general to impose the zero strain rate, D _xx , is given by the constitutive law, and is not an applied stress.

It is clear, then, that the stress formulation, (2), is the natural form to correlate with the test data, since the strainrate formulation, (3), implicitly involves the unknown σ_x . Substituting the strain-rate components and invariants given by (8–10) into (2) yields the following relations:

(15)

(16)

Only two of the axial relations are independent since the traces of both sides of (2) vanish. From (15)₂ and (16)₁,

(17)

and τ is given directly by (16)₂. It now follows that

(18)

which determine both ϕ₁ and ϕ₂ as functions of and from measure , , σ_z and τ in the tests, and in turn as functions of I ₂ and I ₃ using (10). Substituting these expressions for φ ₁ and ϕ ₂ in (17)₁ shows that

(19)

Note that the ratio, τ/ arising in (19), by (16)2, is bounded as → 0. Now (19), a universal relation independent of the response functions ϕ ₁ and ϕ ₂, is a necessary condition that a general isotropic viscous relation (2), or its equivalent form (3), holds. That is, measured σ_x would provide a check on the consistency of the simple viscous fluid assumption, and, while not sufficient to confirm its validity, failure of (19) would show this assumption is not valid. In particular, it could demonstrate that subsequent tertiary creep associated with recrystallization is not an isotropic viscous response, and would necessarily involve dependence on more variables describing the fabric evolution.

Given that law (2) is valid, then a measure of the significance of the quadratic term relative to the commonly assumed parallel relation is the ratio

(20)

where measures the magnitude of D . The common parallel relation is only satisfactory if |Q _cc| << 1 at all of the test points covering the range of strain rates arising in ice-sheet flows. In this case, though, (18)₁ still determines ϕ ₁ as a function of the two invariants, I ₂ and I ₃ which is more general than the inverse of the commonly assumed form (7), already noted to be inconsistent with Reference SteinemannSteinemann’s (1954) combined shear and compression data.

Unconfined compression

With no restraint stresses on flow in any direction normal to the compression axis, the stress, deviatoric stress and the quadratic tensor in (3) have the forms

(21)

(22)

(23)

where –σ_z is the axial compressive stress, τ is the shear stress and the pressure and invariants are

(24)

The distinct xx and yy components in (23) imply that with the quadratic term present in (3), there must be corresponding distinct strain-rate components, so the strain-rate tensor has the form

(25)

where the latter relation is the incompressibility condition, and the invariants are

(26)

In this case, with _z and being the measured strain rates, it is clear that the strain-rate formulation (3) is the most natural for correlation, allowing the unknown _x and _y to be given by the constitutive law in terms of known stresses σ_z and τ. Thus

(27)

(28)

where only two of the axial relations are independent since the traces of both sides of (3) vanish. From (28),

(29)

which, using (24), determines ψ ₁ and ψ ₂ as functions of J ₂ and J ₃. Then from (27),

(30)

which are universal relations independent of the response functions, so measured ɛ_x or ɛ_y would provide a check on the consistency of the simple viscous fluid assumption. The significance of the quadratic term is measured by the ratio

(31)

where measures the magnitude of the deviatoric stress , and the parallel relation is only satisfactory if |Q _cu| << 1 at all of the test points covering the range of expected deviatoric stresses arising in ice-sheet flows.

Biaxial Stress Analysis

In the absence of shear, the stress and strain-rate tensors are simply diagonal with diagonal components σ_x , σ_y , σ_z and ε_x , ε_y , ε_z , respectively. It will be seen that both laterally confined and unconfined longitudinal compression configurations do not yield as much information as the combined compression and shear configurations. It is supposed that σ_y and σ_z are the known applied stresses.

Laterally confined compression

With no compression in the x direction, the strain-rate tensor and quadratic tensor in (2) have the diagonal forms

(32)

where – is the rate of axial compression and the invariants are

(33)

Immediately, there is only a single strain-rate component and one non-trivial invariant, so this configuration covers only one axis of the two-invariant plane. The associated diagonal stress and deviatoric stress tensors, with a necessary constraint stress, σ_x , have the forms

(34)

where p = –(σ_x + σ_y + σ_z /3. The two independent relations of (2) give

(35)

where (35)₁ determines ϕ ₂, and

(36)

These determine both ϕ ₁ and ϕ ² as functions of from measured , σ_z , σ_y , σ_x in the tests, but cannot distinguish dependence on I ₂ and I ₃. Note that the constraint, σ_x , is not determined by the viscous relation in the form (2), but the two independent relations for given by the form (3) would implicitly determine σ_x in terms of σ_y and σ_z . Further, there is no relation independent of the response functions to check the consistency of the simple viscous fluid assumption. Given that the law (2) is valid, then a measure of the significance of the quadratic term relative to the commonly assumed parallel relation is the ratio

(37)

Laterally unconfined compression

With no constraint in the x direction, the stress and deviatoric stress tensors now have the diagonal forms

(38)

and the quadratic tensor in (3) has the form

(39)

where J ₂ and the strain-rate tensor are given by

(40)

By (3),

(41)

(42)

from which

(43)

(44)

but there is no relation independent of the response functions to check the consistency of the viscous fluid assumption. The significance of the quadratic term is measured by the ratio

(45)

Reduced Model for Ice-Sheet Flow

With two distinct strain-rate dependent terms in the deviatoric stress law (2) it is sensible to express this in dimensionless variables appropriate to the ice-sheet flow equations, following Reference Morland and JohnsonMorland and Johnson (1980) and Reference MorlandMorland (1984), to expose the two small parameters which arise, instead of introducing directly a small aspect ratio magnitude (Reference HutterHutter, 1983). First express (2) in the normalized form:

(46)

(47)

(48)

where the parameters σ₀ and D ₀ are typical deviatoric stress and strain-rate magnitudes in ice-sheet flows, so now ϕ ₁ and ϕ ₂, retaining the same response function notation, are order unity. It is supposed that the relation is not dominated by the quadratic term, so that ϕ ₁ is certainly order unity, though ϕ ₂ could be smaller in magnitude. The ice is assumed incompressible with density ρ = 918 kg m^–3. Now let the coordinates, x, y, z, and the velocity components, u, v, w, be dimensionless with units d ₀ and q ₀, respectively, which are respectively an ice-sheet thickness magnitude and surface accumulation magnitude, so the dimensionless velocity gradient and strain rate, have unit q ₀/d ₀, and define dimensionless stress, , with unit a typical overburden pressure ρgd ₀, where g is the gravity acceleration, then

(49)

(50)

The viscous relation (46) is now expressed by

(51)

(52)

and with the parameters above,

(53)

(54)

Both ɛ and the ratio v ₂/v ₁ are inversely proportional to d ₀, so the ratio v ₂/ɛ ³ is independent of d ₀, and, since d ₀ is the only parameter above which may change much between ice sheets, then in general v ₂ = O(ɛ³); that is,

(55)

In relation (51), ɛ ² still defines the very small dimensionless viscosity magnitude, not changed by the ϕ ₂ term, and again (Reference Morland and JohnsonMorland and Johnson, 1980; Reference MorlandMorland, 1984) the small parameter ɛ is the basis for the asymptotic expansions. To obtain a finite-span ice-sheet profile from leading-order balances, it is necessary to introduce the coordinate and velocity scalings

(56)

where the upper-case variables are order unity, and derivatives with respect to X, Y, Z do not change orders of magnitude from that of the dependent variable. The corresponding asymptotic analysis when bed topography of small slope, but larger than ɛ, induces greater local magnitudes, presented by Reference MorlandMorland (2000, Reference Morland, Straughan, Greve, Ehrentraut and Wang2001) for plane flow and Reference Cliffe and MorlandCliffe and Morland (2002) for radially symmetric flow, follows for three-dimensional flow, and will extend to the general viscous relation (51). Here it is supposed that the bed topography is as smooth as the surface and does not induce greater magnitudes, when the reduced model for the simple viscous relation (7) has error O(ɛ ²). With the scalings (56),

(57)

are all order unity, while

(58)

are both O(ɛ ^–1) with the leading-order expressions shown having relative errors O(ɛ ²). It now follows, that, to leading order with relative errors O(ɛ ²),

(59)

where two of the relations are exact. That is, the four leading-order components shown are all O(ɛ ^–2), while the two exact components are both O(ɛ ^–1).

Recall that with temperature variation, the strain rate, D, must be replaced by the effective strain rate, , defined by (6), incorporating the rate factor, a(T ), which can become very small in cold upper regions. However, Reference MorlandMorland (1984) argued that this was compensated by small strain rates in the upper regions so that the formal expansion in ɛ, ignoring the magnitude of a, is valid. The stress estimates presented now ignore the factor a, which would need to be incorporated into the final explicit relations through (6). Now from (51), with relative error O(ɛ ²),

(60)

are all O(ɛ), but are O(ɛ ²) when ϕ ₂ = 0. Also, with relative error O(ɛ) from the neglected ϕ ₂ contribution,

(61)

are both O(ɛ) and independent of φ ₂, and are identical to the expressions with the parallel relation. By (60) and (61), with relative error O(ɛ ²),

(62)

where P is the dimensionless pressure corresponding to Equation (1) ₂, which is again the parallel relation result.

Now examine the momentum balances, equilibrium equations due to the very small inertia terms. Let the z axis point vertically upwards, then

(63)

With the above leading-order stress expressions, the momentum balances become, with error O(ɛ),

(64)

where the latter vertical balance has error O(ɛ ²). These are the same leading-order balances as when the simple parallel relation is used, and by (61) are independent of ϕ ₂, so the balances depend only on ϕ ₁. However, if ϕ ₁ is correlated with single stress test data when ϕ ₂ is actually significant, these balances lead to invalid velocity fields, and hence to an invalid sheet profile.

To lead order, now neglecting O(ɛ), the traction-free surface, Z = H(X, Y, t), conditions are

(65)

where t denotes a dimensionless time with units d ₀/q ₀. Equation (64), then (63), can be explicitly integrated with respect to Z from the surface to yield the familiar stress expressions

(66)

then applying relations (61) and (58),

(67)

The rate factor, a(T), simply enters as a multiplying factor in both velocity gradient expressions (67). Note that can only be replaced by when the parallel relation holds, and in general depends also on φ ₁ in relation (3);that is, the quadratic term contribution enters the reduced model equations. Since depends on all the strain-rate components, and hence all velocity gradients, to leading order, the two equations (67) cannot be formally integrated with respect to Z to determine expressions for U and V. Further, if is expressed in terms of and where and by numerically inverting (2), and depend on all the deviatoric stress components to leading order, not simply on σ_xz and σ_yz given by (66), so again (67) cannot be formally integrated. That is, the reduced model proceeds as far as expressing the pressure and horizontal shear stress components in terms of the surface profile, but allows no further explicit integration to determine velocities and apply surface and bed kinematic conditions. The same situation arises with the coaxial relation if the response functions ϕ ₁ and ψ ₁ both depend on two invariants, which was indicated as a minimum generalization of the common relation (7) by Reference SteinemannSteinemann (1954).

Conclusions

It has been shown that the two types of combined shear and compression tests, longitudinally confined and unconfined, can determine the two response functions of two invariant arguments arising in a general viscous fluid relation for an incompressible material. Further, they yield universal relations independent of the response functions which check the consistency of the viscous fluid assumption, though not its verification, and which can reject its validity. They also provide a direct assessment of the significance of the quadratic tensor term compared to the linear tensor term. Biaxial stress tests are not as useful. The combined shear and compression tests can be used to determine, or reject, a general viscous relation for the deviatoric response of ice. In particular, the commonly adopted relation in which the deviatoric stress is coaxial with the strain rate, the linear tensor term, can be judged as an adequate approximation, or not. Allowing that the quadratic term is significant, it would make significant contributions in the momentum and energy balances of ice-sheet flow. Current modelling adopts the simple coaxial viscous relation. The commonly adopted reduced model (shallow-ice approximation) equations are the leading-order balances of an asymptotic approximation using the small surface slope magnitude or aspect ratio as small parameter, and are valid only when the bed topography also has small slope; the full equations are necessary otherwise. This asymptotic analysis has now been extended to the general viscous relation to show that the relative stress magnitudes are changed, but that the leading-order momentum balances are the same, allowing explicit depth integration to yield the same expressions in terms of the surface profile for the pressure and horizontal shear stresses, but with greater error magnitude. However, the resulting expressions for the velocity gradients with depth can no longer be integrated explicitly, so the velocities and surface and bed kinematic conditions cannot be applied to yield a reduced equation for the surface profile. The same problem arises if the viscous relation is a coaxial form, but with dependence on two invariants. When the bed topography has finite slope, which is the real situation, and the full flow equations are required, a quadratic term will have a significant influence, so its assessment, and determination if significant, is crucial for ice-sheet flow modelling.