1. Introduction
The constant stress and constant strain-rate responses of ice at constant temperature in uniaxial compressive stress exhibit non-viscous behaviour: the strain rate is not a unique function of stress, and vice versa. At constant stress the initial primary strain rate decreases with time to a minimum, described as secondary creep. It then increases and approaches an asymptotic limit, called tertiary creep. Analogously, at constant strain rate the initial stress increases to a maximum then decreases to an asymptotic limit.
Reference Mellor and TrydeMellor (1980) noted that the initial elastic strain when a stress is applied is much smaller than the creep strains, so this solid response will be neglected. Such variations with time at constant temperature can then be described by a viscoelastic fluid model, and are used here to construct a simple constitutive law of differential type.
The timescales of the primary, secondary and tertiary response increase significantly with a decrease in temperature. Following Reference MorlandMorland (1979) it is proposed that the ice is thermo-rheologically simple (Reference Morland and LeeMorland and Lee, 1960); that is, there is a common response on a pseudo-timescale depending on temperature, analogous to the temperature-dependent rate factor applied in the conventional viscous law. There is a considerable variation of temperature through an ice sheet. The time that different ice elements have been stressed for ranges from very short near the surface in an accumulation zone, to very long as the ice passes through the sheet. There are concerns about whether a flow solution based on the conventional viscous response can realistically reflect the time-dependent viscoelastic behaviour, which itself has such significant variation through the sheet. Only comparison flow solutions can determine the significance of such viscoelastic effects. The present constitutive model, which is the simplest consistent with the uniaxial stress response, could be a basis for preliminary comparisons.
The idealized responses adopted above have been widely discussed by Reference HookeHooke (1980), Reference JackaJacka (1984a,b) and Reference Jacka and MaccagnanJacka and Maccagnan (1984), where experimental results show times to minimum strain rate at different temperatures, and by Reference Budd and JackaBudd and Jacka (1989) and Reference Jacka, Li and HondohJacka and Li (2000). However, the responses are also shown as functions of the evolving strain, and interpreted as variations with evolving fabric. Fabric evolution models predict distinct behaviour in the uniaxial response, with viscosity enhancement in tertiary creep greater than unity. This is in contrast to the decreased viscosity described above, suggesting that explicit time-dependent viscoelastic behaviour is required to model the above responses instead of, or in addition to, fabric evolution. There has been no clear separation of the two processes to date, and this is discussed in section 8. The influence of such time-independent fabric evolution accompanied by induced anisotropy is ignored in the present, preliminary, simple construction.
Following Reference Mellor and TrydeMellor (1980) and the other references above, Figure 1 shows qualitatively the axial compressive strain-rate r(t) response in time t at constant uniaxial compressive stress σ and constant temperature. Figure 2 shows qualitatively the uniaxial compressive stress σ(t) response at constant axial compressive strain rate r and constant temperature. As the constant temperature increases, the strain rate at each constant stress increases and the primary and tertiary creep timescales decrease. In Figure 1, the curve rises and shifts left with increasing temperature. In addition, the stress at each constant strain rate decreases, so in Figure 2 the curve lowers and shifts left with increasing temperature.
Fig. 1. Axial strain-rate r response in time t at constant axial stress σ.
Fig. 2. Axial stress σ response in time t at constant axial strain rate r.
In Figure 1, the initial strain rate is r 0(σ) which then decreases during primary creep to a minimum strain rate (secondary creep) rm (σ) at time tm (σ). It then increases during tertiary creep to an asymptotic limit rate re (σ). The time τ(σ) denotes the time at which the strain rate has increased a factor 1 − δ of re − rm from rm , where 0 < δ ≪ 1 is a chosen parameter. The parameters r 0(σ), tm (σ), re (σ) and τ(σ) describe significant features of the responses for varying σ, and may be the main inferences obtained from experimental data.
In Figure 2, the initial stress is σ 0(r) which then increases to a maximum stress σM (r) at time tM (r), then decreases to an asymptotic limit σE (r). The parameters σ 0(r), tM (r) and σE (r) describe significant features of the responses for varying r, and may be the main inferences obtained from experimental data. The present analysis supposes that
but analogous analysis applies if the inequalities are reversed. The reciprocal relations
are reasonable assertions (Reference Morland and SpringMorland and Spring, 1981), and Reference Mellor and TrydeMellor (1980) also suggests
These define unique boundaries in the r − σ domain, shown in Figure 3 as the three curves r expressed as a function of σ and conversely, separating regions in which the coefficient functions of the viscoelastic differential relation have specific values.
Fig. 3. Domain of axial strain rate r vs stress σ, with primary, secondary and tertiary creep boundaries.
It is conjectured that as the constant stress σ → 0 then r 0, rm , re → 0 and as the constant strain rate r → 0 then σ 0, σM , σE → 0. However, their general dependence on σ and r, respectively, must be determined from test data. Reference Mellor and TrydeMellor (1980) proposed that the responses from both tests must reflect the same material properties, but that no existing constitutive model simulates both responses. Reference MorlandMorland (1979) and Reference Morland and SpringMorland and Spring (1981) have proposed viscoelastic fluid models of differential type which exhibit the main features of both responses. Correlation with one type of response still required information from the second type, however, and would not allow exact correlation with an independent second type response.
In this paper, after presenting the axial differential relations which describe the constant stress and constant strainrate responses shown in Figures 1 and 2, a corresponding frame-indifferent viscoelastic fluid law is constructed. This law is simpler than the earlier laws, and is correlated with independent constant stress and constant strain-rate responses. That is, data from constant stress and constant strain-rate tests determine the viscoelastic fluid law of this form. An idealized family of responses is adopted to demonstrate the correlation process and the determination of the constitutive (material) functions which define the law. A simple illustration is also presented.
2. Constant Stress Response
It is first necessary to show how the primary, secondary and tertiary creep responses with time at constant stress, then at constant strain rate, can be described by differential equations which can be correlated with a viscoelastic fluid law of differential type.
Following Reference Morland and SpringMorland and Spring (1981), note that the strain rate r(t) at constant stress shown in Figure 1 repeats the range rm ≤ r ≤ re in the primary and tertiary creep zones, with strain acceleration r negative and positive, respectively. The latter has the simplest differential equation description
where the strain acceleration r is a material time derivative of r; that is, the time derivative of r following the ice particle. The quadratic relation, Equation (5), has negative and positive roots:
which describe monotonic primary and tertiary creep, respectively, provided that
There are then unique inverse relations
For 
and 
as r → rm
(σ), and for 
to switch from 
to 
as t passes through tm
(σ), it is necessary that
The asymptotic limit r → re
(σ) as t → ∞ requires 
as r → re
(σ), which implies
Since there is only primary creep 
in r ≥ re
(σ), that is, only a single root of Equation (5),
Further, neither root given by Equations (6) and (7) applies in r < rm (σ) nor in r > r 0(σ), so it is convenient to define the extensions, with f(r, σ) discontinuous at r = r 0(σ),
These explicit properties are noted in Figure 3. The time tm (σ) to minimum strain rate rm (σ) is bounded and the time to the tertiary limit re (σ) is unbounded, satisfied if
A complete family of constant stress responses at constant temperature over an appropriate stress range is required to determine the functions f(r, σ) and F(r, σ). Equations (13) and (14) cover the ranges r < rm
(σ) and r > re
(σ). Given
in Equations (6) and (7) and using the inverse relations of Equations (9) and (10),
3. Constant Strain-Rate Response
There is an analogous analysis for the constant strain-rate response shown in Figure 2. The governing differential equation for σ(t) is
where 
is the material time derivative of σ. The primary and tertiary roots are
subject to G[r, σE (r)] = 0 and
There are then unique inverse relations:
Since there is only primary creep 
in σ ≤ σE
(r),
Again, we define the continuous extensions
which are noted in Figure 3. The time tM (σ) to maximum stress σM (r) is bounded and time to the tertiary limit σE (r) is unbounded, satisfied if
A complete family of constant stress responses at constant temperature over an appropriate strain-rate range is required to determine the functions g(r, σ) and G(r, σ). Equations (24) and (25) cover the ranges σ < σE
(r) and σ > σM
(r), and then, given 
in Equations (18) and (19) and using the inverse relations of Equations (22) and (23),
4. Viscoelastic Fluid Law
With the usual incompressibility assumption, a constitutive viscous law for ice determines the deviatoric stress
where σ is the Cauchy stress and p is the mean pressure. In terms of the instantaneous strain rate,
where v ( x , t) is the velocity field in spatial coordinates x , grad denotes the spatial gradient and T denotes tensor transpose. To incorporate the strain acceleration required in the constant uniaxial stress response, Reference MorlandMorland (1979) introduced the frame-indifferent second Rivlin–Ericksen tensor. To describe the constant strain-rate response, Reference Morland and SpringMorland and Spring (1981) introduced a frame-indifferent deviatoric stress-rate tensor. However, these responses can both be simulated without such rate tensors, by simply allowing dependence on material time derivatives of strain rate and deviatoric stress invariants. This yields a simpler form of tensor relation than that proposed by Reference Morland and SpringMorland and Spring (1981). Such a form is now analysed to show how it can be determined by independent data from constant stress and constant strain-rate tests.
The simple form proposed, to describe the response at constant temperature T, is given by the equivalent coaxial relations
where I and J are the second invariants of
D
and
S
, respectively, and 
and 
are their respective material time derivatives:
Φ(I, J, , ) and Ψ(I, J, , ), related, are the response functions which describe the particular viscoelastic properties of the ice. The first-order differential equations relating
D
and
S
, Equations (30), must be completed by initial conditions on I and J. Dependence on the third invariants of
D
and
S
, and/or their time derivatives, may be required for correlations with more general stress responses.
For uniaxial compressive stress σ, with axial compressive strain rate r accompanied by lateral (tensile) strain rates r /2 in all lateral directions to satisfy incompressibility,
At constant σ, each of the diagonal components of Equation (30) becomes
and at constant r,
For comparison, Equations (5) and (17) can be written with expressed in terms of and expressed in terms of :
Define
where by Equations (32) and (33), σ and r are replaced by σ = (3J)1/2, r = 2(I/3)1/2.
That is, f(r, σ), F(r, σ), g(r, σ) and G(r, σ) given by the primary and tertiary response data relations (13), (16), (24) and (27) determine the functions 
, 
, 
and 
.
Now Φ(I, J, , ) must be chosen so that Equation (34) at constant σ, , yields Equation (36) and Equation (35) at constant r, , yields Equation (37). A simple construction is the additive form
where A(I, J) and B(I, J) are functions of I and J subject to the restriction
consistent with the vanishing of F and G in the same zones.
This is not a unique prescription of a consistent Φ, but appears to be the simplest form. Here, given independent 
and 
, there is still a single arbitrary function A(I, J) or B(I, J) which could be used to ‘improve’ correlation with a different response, such as simple shear. If the family of constant stress responses contain the same information as the family of constant strain-rate responses, then 
and are not independent of 
and but the above relations still apply. There is an analogous construction for the response function Ψ.
5. Temperature Dependence
The significant temperature dependence is described by assuming that the ice is thermo-rheologically simple, satisfying the above time-dependent relations with real time t replaced by a pseudo-time ζ( X , t) (Reference Morland and LeeMorland and Lee, 1960; Reference MorlandMorland, 1979) for each material element X , depending on the temperature history of that element, defined by
Here t 0 is an initial time before which the ice is not stressed or when a steady state has been maintained for sufficient time. It states that all rate processes are increased at each time for each element by a factor a[T ( X , t)] which is the analogue of the temperature-dependent rate factor applied in the usual viscous law. Reference Smith and MorlandSmith and Morland (1981) derived an accurate representation for a(T) from data presented by Reference Mellor and TestaMellor and Testa (1969) for temperatures between melting and 60 K below melting. A good simplifying approximation over the range of practical significance from melting to 40 K below melting is (Reference MorlandMorland, 1997)
where [20 K] is a typical temperature change magnitude over an ice-sheet depth and the dimensionless temperature 
is zero at melting and −2 at 40 K below melting where a = 7.9×10−
4. The present uniaxial stress analysis supposes the ice is near the melting point, when a = 1 and ζ = t with respect to t
0 = 0. At different constant temperatures, t is replaced by ζ = at.
For a general temperature field with rate factor a[T ( X , t)], with time t replaced by ζ in the viscoelastic fluid relations,
and Equations (30), (41) and (42) become
At low temperatures, the rate factor a ≪ 1 significantly affects the relative weights of the terms in the square brackets and of the outside terms of Φ. There still remains an arbitrary function 
or B to improve other correlations, but recall that the uniaxial responses do not determine , , or over a complete I − J domain.
Note that the material time derivative is required in Equation (43), as in Equation (31) for the invariants. The former involves the gradient in spatial coordinates necessary to describe an ice-sheet domain subject to surface accumulation/ ablation and basal melting, steady or evolving. This adds a numerical complexity to the flow equations, adding hyperbolic equations to the elliptic equations of the slow flow. However, this situation has been addressed successfully for a steady radial flow incorporating fabric evolution (Reference Morland and StaroszczykMorland and Staroszczyk, 2006).
6. Idealized Responses
The form of strain-rate response r(t) at each constant stress σ and constant temperature T shown in Figure 1 can be simulated by
where b and c are functions of σ, satisfying r = rm
and 
at 
and r → re
as 
. The strain acceleration is then
which must be positive during tertiary creep 
and negative during primary creep 
. This requires that
The condition that re − r = δ Δ at t = τ, that is, at 
where 
(a measure of how fast the tertiary creep approaches its limit) is
Setting r 0 − re = k Δ, which relates the strain-rate decrease between its initial value and final tertiary limit value to the tertiary increase from the minimum value, provides another measure of the overall response. This is expressed
Finally, following Reference Mellor and TrydeMellor (1980) and Reference JackaJacka (1984a), the strain ∊m at tm is 0.01 for all σ, so
where
The Equation (54) restriction is adopted in an illustration, but it is noted that this magnitude of ∊m applied to the responses at or near the melting temperature yields timescales that are too short.
Setting
the relations (52–54) with Equation (55) become
subject to the Equation (51) restrictions:
Since β(γ) is monotonically increasing for each k, and β(1) < 0 and β → ∞ as γ → ∞, there is a unique γm (k) such that the left inequality is satisfied for γ > γm (k). Further, since β decreases with increasing k at fixed γ, γm (k) increases with k. The right inequality is expressed by
Now h(γ) is monotonically decreasing for each k, and h(1) > 0 and h → −∞ as γ → ∞, so there is a unique γx (k) such that the right inequality is satisfied for γ < γx (k). Further, since h increases with increasing k at fixed γ, γx (k) increases with k. Also, h(γm ) is positive by the second equality of Equation (61), so γx > γm ; that is, there is a finite range (γm , γx ) of γ for each k which satisfies the Equation (60) restrictions. A valid set of parameters requires that the root γ of Equation (57) for given k(σ) and τ(σ) falls in this range. With a valid set, Equation (59) then determines tm (σ) given re (σ).
As a simple example, it is supposed that there is a response 
independent of σ:
The secondary response rm (σ) is the usual viscous relation, and Equation (62) assumes that the strain rate relative to rm (σ) varies on a timescale stretched by tm (σ) independent of σ. While the latter is consistent with all timescales increasing as σ decreases, it introduces an indeterminancy as σ → 0. Then Equation (49) gives
where Re and the coefficients γ and β are constant, independent of σ. Re is the strain-rate enhancement factor, or reciprocal of the viscosity enhancement factor, relating tertiary and secondary responses. Equations (57) and (58) for the constant γ and β, subject to the Equations (60) and (61) validity restrictions, are unchanged, but now
Equation (59) becomes
where the constant ∊⋆ is given by
The tm variation with σ is the reciprocal of that of rm , certainly consistent with behaviour as σ increases and decreases. Given rm (σ), the constants Re , R 0, δ and τ can be specified and ∊⋆ , hence tm (σ), determined. Alternatively, the constant ∊⋆ can be specified and R 0 or τ determined.
The differential relation Equation (5) now becomes
where the symbol ° denotes differentiation with respect to 
or 
, and Equations (13) and (16) give
with 
and 
zero elsewhere. By Equation (63),
and can be expressed in terms of y by the unique =
and 
for the monotonic primary and tertiary creep responses, respectively, defined by Equation (63). The constitutive functions and in the viscoelastic response function Φ of the model relation Equation (41), using Equations (38) and (67), are then
The domains of 
and 
defined by uniaxial stress response are given by
However, other stress configurations may require 
and 
over different (I, J) domains, and may imply dependence on third invariants as the shear experiments of Reference SteinemannSteinemann (1954) did for the viscous law.
These functions are completed by prescribing the secondary-creep minimum strain rate rm (σ) defined by the standard viscous law. This three-term polynomial representation was constructed by Reference Smith and MorlandSmith and Morland (1981) with finite viscosity at zero stress, closely correlated with Reference GlenGlen’s (1955) uniaxial compression data at near melting temperature, a = 1, over a stress range 0–8.66×105 N m− 2. With a stress unit of 105 N m− 2 (0.1MPa) and a strain-rate unit of a− 1, their relation is
With Equations (32) and (33), this gives
The choice of idealized constant stress responses (Equation (62)) also determines the constant strain-rate responses, since it applies for all σ and is an implicit equation for σ(t) at each constant r. That is,
using Equation (4) to invert rm
(σ). 
can be expressed in terms of r and y, hence σ by Equation (62), on the separate primary and tertiary monotonic creep responses 
and t > 1. Now differentiating Equation (62) with respect to t at constant r gives
which vanishes when 
, that is, at 
, t = tm
(σ). Further, 
has the opposite sign to 
and is therefore positive in 
and negative in 
. This verifies the primary monotonically increasing σ in t < tm
(σ) and the tertiary monotonically decreasing σ in t > tm
(σ). The functions g (r, σ) and G(r, σ) are then determined by Equation (27) and 
and 
by Equation (39).
7. Illustration
Sets of parameters satisfying Equations (57–59) and consistent with the restriction Equation (60) have been calculated, and the following are chosen to illustrate the response and determine the corresponding constitutive model response functions:
Figure 4 shows three associated constant stress responses at
Fig. 4. Example strain rates r(t) at three constant stresses.
For σ = 1 the initial strain rate is r = 1.18, but then decreases too rapidly in time to be shown. The short timescale is due both to the chosen idealized response given by Equations (62) and (63) and the Mellor suggestion and ∊m = 0.01 applied near melting, with a = 1. An ∊m increasing with temperature as does the strain rate, with larger rate of increase near melting, increases ∊⋆ given by Equation (65) and hence tm (σ) given by Equation (66).
Figure 5 shows four constant strain-rate responses, again with very short timescales, at
Fig. 5. Example stresses σ(t) at four constant strain rates.
The constitutive functions 
and 
are formally determined by Equations (68) and (70–72), and 
and 
by Equations (24), (27) and (39) once the time in the separate primary and tertiary creep responses is expressed in terms of r and σ. For the example Equation (62), this is simply determining the primary and tertiary roots 
and 
of
for each r and σ, hence I and J, over the chosen ranges. Here, these are 0 ≤ σ ≤ 3 and 0 ≤ r ≤ 3 corresponding to 0 ≤ J ≤ 3 and 0 ≤ I ≤ 6.75. There is an indeterminancy as r, σ → 0, no root , for this example. The maximum function values are very large due to the short timescale of the example and the very small ∊⋆
. They are shown as surfaces of 
and 
in Figures 6 and 7 respectively, and surfaces of ∊⋆g(I,J) and (∊⋆
)2
G(I,J) in Figures 8 and 9 respectively. Figures 6–9 also show their common non-zero (I, J) domains which are covered by uniaxial stress. Analyzing actual data from families of constant stress and constant strain-rate responses to determine these constitutive functions will be much more complicated, but will necessarily yield the correct timescales associated with the stress, strain rate and temperature ranges covered.
Fig. 6. The response function .
Fig. 7. The response function 
.
Fig. 8. The response function ∊⋆g (I, J).
Fig. 9. The response function (∊⋆ )2 G(I, J).
8. Viscoelasticity, Fabric and Induced Anisotropy
As noted in section 1, constant uniaxial stress response due to fabric evolution (based on crystal rotation arguments and a consequence of fabric models) shows an increased tertiary viscosity over the minimum (secondary) viscosity. This is in contrast to the decreased tertiary viscosity shown by Reference Mellor and TrydeMellor (1980) and many subsequent experimental results. Reference Jacka, Li and HondohJacka and Li (2000) carried out experiments which, at least to the attained strains of approximately 10%, indicate that at low temperatures and stresses the fabric remains random and there is no change of viscosity in the tertiary response.
The fabric seems to be dominated by rotation at low temperatures and higher stresses, while at high temperatures it is dominated by recrystallization. Reference Jacka, Budd, Jones, McKenna, Tillotson and JordaanJacka and Budd (1991) note that theoretical considerations of crystal rotation at high stress by Reference Azuma and HigashiAzuma and Higashi (1985) and Reference AlleyAlley (1988) ignore recrystallization. At the relatively high temperatures and stresses of laboratory tests, it is possible that recrystallization is the cause of the decreased tertiary viscosity rather than the time-dependent viscoelastic response modelled here, though this is not implied directly.
While Reference Budd and JackaBudd and Jacka (1989) indicate that primary creep exhibits time-dependent recoverable strain (a viscoelastic effect), Reference Gao and JackaGao and Jacka (1987) emphasize that further compression of samples previously compressed to the tertiary limit do not exhibit a change of strain rate from the previous tertiary limit; this is a rejection of a viscoelastic effect once the fabric is fully developed. Further, observed fabrics in the tertiary limit of uniaxial compression are not consistent with those predicted by crystal rotation alone, but are consistent with a combination of rotation and recrystallization which can explain the decrease in viscosity.
However, the above inferences do not imply that viscoelastic processes, with explicit time dependence, are not occurring. Existing evidence does not clearly distinguish fabric and viscoelastic effects. Both processes may be taking place, with one or the other dominating in different regimes. Reference Gao and JackaGao and Jacka (1987) also note that applying a different stress configuration in the further compression of a sample previously compressed to its tertiary limit results in a different strain rate. This confirms that fabric evolution does indeed induce anisotropy. If viscoelastic processes are also occurring (given that induced anisotropy is significant), the present isotropic viscoelastic fluid model must be generalized, much as the conventional viscous law has been generalized, to incorporate fabric effects.
It is a major task to perform the necessary experimental tests required to infer an appropriate constitutive model, but a preliminary approach would be to correlate classes of simple idealized models to existing data. There still remains the need to extrapolate experimentally inferred responses to the wide range of temperatures and stresses arising through an ice sheet. However, it is better to extrapolate a model consistent with experimental behaviour than to continue applying the extrapolated viscous fluid model.
9. Conclusions
The present analysis has focused on the constant uniaxial stress and constant strain-rate responses of ice in time, and shown how a simple viscoelastic fluid model can simulate, and be determined by, the two families of responses. The specific simplification adopted is probably too simple to correlate with real families of responses, and the families chosen as examples too artificial. However, the analysis demonstrates how an actual correlation can be constructed given the necessary real data. Different stress configurations may show that the viscoelastic model needs further generalization. More importantly, there is a need to separate explicit time-dependent viscoelastic effects and fabric evolution effects, and it may prove necessary to incorporate both into a model describing real response.
Acknowledgements
I thank T.H. Jacka for advising me, with supporting references, that the assumed time dependence of the experimental responses may be interpreted entirely as strain dependence, and that the decreased tertiary viscosity can be explained by fabric evolution. I also appreciate the constructive comments made by the two referees (K. Hutter, R. LeB. Hooke) and the scientific editor (R. Greve).
177
