Primary, secondary and tertiary creep of ice modelled as a viscoelastic fluid

Abstract

As an ice sheet evolves, there are ice elements near the surface only recently subjected to stress following deposition, and others that have been subjected to stress over many ranges of time. The constant stress and constant strain-rate responses of ice in uniaxial compressive stress exhibit non-viscous behaviour, that is, the strain rate is not fixed by the stress (and conversely) but both vary with time. 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, described as tertiary creep. Analogously, at constant strain rate the initial stress increases to a maximum then decreases to an asymptotic limit. These responses are used to construct a simple viscoelastic fluid constitutive law of differential type. Such a time-dependent law, with timescales changing widely with temperature, can be expected to yield a flow field in an ice sheet that is very different from that obtained from the viscous law. Only comparison solutions for both constitutive laws can determine the differences and significance of the non-viscous behaviour, and the simple law constructed would be a candidate for such comparisons.

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.

Mellor (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 Morland (1979) it is proposed that the ice is thermo-rheologically simple (Morland 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 Hooke (1980), Jacka (1984a,b) and Jacka and Maccagnan (1984), where experimental results show times to minimum strain rate at different temperatures, and by Budd and Jacka (1989) and Jacka 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 Mellor (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 ₀(σ) which then decreases during primary creep to a minimum strain rate (secondary creep) r_m (σ) at time t_m (σ). It then increases during tertiary creep to an asymptotic limit rate r_e (σ). The time τ(σ) denotes the time at which the strain rate has increased a factor 1 − δ of r_e − r_m from r_m , where 0 < δ ≪ 1 is a chosen parameter. The parameters r ₀(σ), t_m (σ), r_e (σ) 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 σ ₀(r) which then increases to a maximum stress σ_M (r) at time t_M (r), then decreases to an asymptotic limit σ_E (r). The parameters σ ₀(r), t_M (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

((1))

but analogous analysis applies if the inequalities are reversed. The reciprocal relations

((2))

((3))

are reasonable assertions (Morland and Spring, 1981), and Mellor (1980) also suggests

((4))

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 ₀, r_m , r_e → 0 and as the constant strain rate r → 0 then σ ₀, σ_M , σ_E → 0. However, their general dependence on σ and r, respectively, must be determined from test data. Mellor (1980) proposed that the responses from both tests must reflect the same material properties, but that no existing constitutive model simulates both responses. Morland (1979) and Morland 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 Morland and Spring (1981), note that the strain rate r(t) at constant stress shown in Figure 1 repeats the range r_m ≤ r ≤ r_e in the primary and tertiary creep zones, with strain acceleration r negative and positive, respectively. The latter has the simplest differential equation description

((5))

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:

((6))

((7))

which describe monotonic primary and tertiary creep, respectively, provided that

((8))

There are then unique inverse relations

((9))

((10))

For and as r → r_m (σ), and for to switch from to as t passes through t_m (σ), it is necessary that

((11))

The asymptotic limit r → r_e (σ) as t → ∞ requires as r → r_e (σ), which implies

((12))

Since there is only primary creep in r ≥ r_e (σ), that is, only a single root of Equation (5),

((13))

Further, neither root given by Equations (6) and (7) applies in r < r_m (σ) nor in r > r ₀(σ), so it is convenient to define the extensions, with f(r, σ) discontinuous at r = r ₀(σ),

((14))

These explicit properties are noted in Figure 3. The time t_m (σ) to minimum strain rate r_m (σ) is bounded and the time to the tertiary limit r_e (σ) is unbounded, satisfied if

((15))

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 < r_m (σ) and r > r_e (σ). Given in Equations (6) and (7) and using the inverse relations of Equations (9) and (10),

((16))

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

((17))

where is the material time derivative of σ. The primary and tertiary roots are

((18))

((19))

subject to G[r, σ_E (r)] = 0 and

((20))

((21))

There are then unique inverse relations:

((22))

((23))

Since there is only primary creep in σ ≤ σ_E (r),

((24))

Again, we define the continuous extensions

((25))

which are noted in Figure 3. The time t_M (σ) to maximum stress σ_M (r) is bounded and time to the tertiary limit σ_E (r) is unbounded, satisfied if

((26))

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),

((27))

4. Viscoelastic Fluid Law

With the usual incompressibility assumption, a constitutive viscous law for ice determines the deviatoric stress

((28))

where σ is the Cauchy stress and p is the mean pressure. In terms of the instantaneous strain rate,

((29))

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, Morland (1979) introduced the frame-indifferent second Rivlin–Ericksen tensor. To describe the constant strain-rate response, Morland 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 Morland 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

((30))

where I and J are the second invariants of D and S , respectively, and and are their respective material time derivatives:

((31))

Φ(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,

((32))

((33))

At constant σ, each of the diagonal components of Equation (30) becomes

((34))

and at constant r,

((35))

For comparison, Equations (5) and (17) can be written with expressed in terms of and expressed in terms of :

((36))

((37))

Define

((38))

((39))

where by Equations (32) and (33), σ and r are replaced by σ = (3J)^1/2, r = 2(I/3)^1/2.

((40))

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

((41))

where A(I, J) and B(I, J) are functions of I and J subject to the restriction

((42))

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) (Morland and Lee, 1960; Morland, 1979) for each material element X , depending on the temperature history of that element, defined by

((43))

Here t ₀ 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. Smith and Morland (1981) derived an accurate representation for a(T) from data presented by Mellor 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 (Morland, 1997)

((44))

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. 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,

((45))

and Equations (30), (41) and (42) become

((46))

((47))

((48))

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 (Morland 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

((49))

where b and c are functions of σ, satisfying r = r_m and at and r → r_e as . The strain acceleration is then

((50))

which must be positive during tertiary creep and negative during primary creep . This requires that

((51))

The condition that r_e − r = δ Δ at t = τ, that is, at where (a measure of how fast the tertiary creep approaches its limit) is

((52))

Setting r ₀ − r_e = 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

((53))

Finally, following Mellor (1980) and Jacka (1984a), the strain ∊_m at t_m is 0.01 for all σ, so

((54))

where

((55))

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

((56))

the relations (52–54) with Equation (55) become

((57))

((58))

((59))

subject to the Equation (51) restrictions:

((60))

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

((61))

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 t_m (σ) given r_e (σ).

As a simple example, it is supposed that there is a response independent of σ:

((62))

The secondary response r_m (σ) is the usual viscous relation, and Equation (62) assumes that the strain rate relative to r_m (σ) varies on a timescale stretched by t_m (σ) independent of σ. While the latter is consistent with all timescales increasing as σ decreases, it introduces an indeterminancy as σ → 0. Then Equation (49) gives

((63))

where R_e and the coefficients γ and β are constant, independent of σ. R_e 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

((64))

Equation (59) becomes

((65))

where the constant ∊^⋆ is given by

((66))

The t_m variation with σ is the reciprocal of that of r_m , certainly consistent with behaviour as σ increases and decreases. Given r_m (σ), the constants R_e , R ₀, δ and τ can be specified and ∊^⋆ , hence t_m (σ), determined. Alternatively, the constant ∊^⋆ can be specified and R ₀ or τ determined.

The differential relation Equation (5) now becomes

((67))

where the symbol ° denotes differentiation with respect to or , and Equations (13) and (16) give

((68))

with and zero elsewhere. By Equation (63),

((69))

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

((70))

((71))

((72))

The domains of and defined by uniaxial stress response are given by

((73))

((74))

However, other stress configurations may require and over different (I, J) domains, and may imply dependence on third invariants as the shear experiments of Steinemann (1954) did for the viscous law.

These functions are completed by prescribing the secondary-creep minimum strain rate r_m (σ) defined by the standard viscous law. This three-term polynomial representation was constructed by Smith and Morland (1981) with finite viscosity at zero stress, closely correlated with Glen’s (1955) uniaxial compression data at near melting temperature, a = 1, over a stress range 0–8.66×10⁵ N m^{− 2}. With a stress unit of 10⁵ N m^{− 2} (0.1MPa) and a strain-rate unit of a^{− 1}, their relation is

((75))

With Equations (32) and (33), this gives

((76))

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,

((77))

using Equation (4) to invert r_m (σ). 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

((78))

which vanishes when , that is, at , t = t_m (σ). Further, has the opposite sign to and is therefore positive in and negative in . This verifies the primary monotonically increasing σ in t < t_m (σ) and the tertiary monotonically decreasing σ in t > t_m (σ). 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:

((79))

Figure 4 shows three associated constant stress responses at

((80))

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 t_m (σ) given by Equation (66).

Figure 5 shows four constant strain-rate responses, again with very short timescales, at

((81))

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

((82))

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 (∊^⋆ )² 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 (∊^⋆ )² 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 Mellor (1980) and many subsequent experimental results. Jacka 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. Jacka and Budd (1991) note that theoretical considerations of crystal rotation at high stress by Azuma and Higashi (1985) and Alley (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 Budd and Jacka (1989) indicate that primary creep exhibits time-dependent recoverable strain (a viscoelastic effect), Gao 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. Gao 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