Introduction
Civil engineering problems involving ice as a damaging material have been recently addressed. At very high strain rates brittle behaviour is predominant and existing models for continuous damage of rocks or ceramics (and possibly concrete) can be adapted without a big risk. At moderate strain rates, of the order of 10−1s−1, the ductile behaviour of ice cannot be neglected. If the steady creep law of undamaged ice is now well established, at least in the range of stresses relevant to engineering problems, including it in a damage model is not straightforward (Reference SinhaSinha, 1988: Reference Meyssonnier, Duval, Axelsson and FranssonMeyssonnier and Duval, 1989). Concerning the ductile brittle transition zone, which is probably involved in most of the ice-structure interaction events, the transient behaviour of undamaged ice has been the object of extensive experimental and theoretical studies (Reference Traetteberg, Gold, Fiederking and FrankensteinTraettebergand others, 1975; Reference DuvalDuval. 1978; Reference SinhaSinha, 1978; Reference Cole and SodhiCole, 1991) and still needs to be investigated (Reference GoldGold, 1994). So far, two kinds of models have been published which describe the transient creep of ice by means of creep functions (Reference SinhaSinha. 1978) or internal state variables (Reference Le Gac, Duval and TrydeLe Gac and Duval, 1980; Reference Sunder and WuSunder and Wu, 1989a). Although these models are quite different from each other, their respective authors have shown them to work well in simple loading cases. In our opinion, their ability to reproduce ice behaviour correctly in complex loading situations remains to be proven.
The present paper focuses only on models involving internal variables, since they lead to constitutive equations which can be easily implemented into numerical codes. They can summarize and bring out, at the macroscopic level, the microscopic processes that control the ice deformation, a brief review of which is given below.
Deformation mechanisms of polycrys-talline ice
Single crystal deformation
The ice Ih single crystal deforms essentially by slip on the basal plane (0001) normal to the hexagonal symmetry-axis <0001>. Crystals well oriented for basal glide exhibit an accelerating transient creep until a steady creep rate, proportional to the square of the stress, is reached. The accelerating transient creep is explained by the increase in the dislocation density up to a steady value (Reference WeertmanWeertman, 1983). According to Reference Duval, Ashby and AndermanDuval and others (1983) the non-basal deformation under a prescribed strain rate at –10°C, requires stresses 60 times larger than that for an easy-glide oriented crystal. Data reviewed by these authors show that the steady creep rate of easy-glide oriented crystals, subject to a constant stress, can be more than 1000 times that of polycrystalline ice.
Polycrystal deformation
Under stresses larger than 0.05 MPa at –10°C (Reference Duval, Ashby and AndermanDuval and others, 1983), the deformation of polycrystalline ice is due to the movement of dislocations in each grain. Because of the strong anisotropy of ice single-crystals, the grains of polycrystalline ice have only two slip systems available for deformation, and these are in the basal plane. Reference Duval, Ashby and AndermanDuval and others (1983) suggest that dislocation climb on planes normal to the basal plane, and possibly on the prismatic planes, is the most likely process to accommodate an imposed deformation. Since the dislocations experience a small resistance to glide motion, the steady creep rate is then controlled by dislocation climb.
According to Reference Duval, Ashby and AndermanDuval and others (1983), the initial strain rate of primary creep corresponds to that of the easy-glide oriented grains. Then, because of the grain anisotropy, the macroscopically isotropic polycrystal has to cope with the strain incompatibilities of easy-glide versus hard oriented crystals, which cause an increasingly non-uniform stress field to develop at the grain scale. These internal stresses oppose forward deformation, giving kinematic hardening. The recoverable strain observed when unloading a polycrystal corresponds to the relaxation of the internal stresses. It is many times the pure elastic strain corresponding to the same loading conditions (Reference DuvalDuval, 1978; Reference SinhaSinha, 1978; Reference Ashby and DuvalAshby and Duval, 1985; Reference Cole and SodhiCole, 1991). The energy storage inside the grains, the release of which permits strain recovery, is generally explained by the bowing of dislocation segments in a substructure network or by the piling of dislocations against obstacles such as grain boundaries and sub-boundaries (Reference PoirierPoirier, 1977). The recovery processes limiting the stored energy level are, beside recrystalliz-ation and microcracking which result in irreversible changes of the structure, dislocation climb and possibly grain-boundary sliding.
According to Reference SinhaSinha (1978), the delayed elastic strain has its origin in the latter mechanism. Sliding occurs at grain boundaries until local stresses, due to irregularities and triple junctions, increase to equilibrate the applied stress. Crystal accommodation being purely elastic, the resulting recoverable strain should be of the order of the pure elastic strain. This reasoning is consistent with Reference SinhaSinha's (1978) observations on columnar ice for short loading durations (10 min).
Along with the long-range processes which yield kinematic hardening, short-range interactions between the moving dislocations and the crystal lattice generate friction forces which act independent of the direction of motion. According to Reference Duval, Ashby and AndermanDuval and others (1983), this isotropic hardening is responsible for the zero creep rate periods they observed during stress drop experiments. The relaxation of the internal stresses associated with isotropic hardening is linked to the recovery of the dislocation substructures.
Internal-state variable models for ice
An interesting feature of such models is that they can be handled using the general framework of thermodynamics. They consist of constitutive and evolution equations which involve stress, strain and state variables which characterize the current state of work-hardening. The present study focuses on the models of Reference Le Gac, Duval and TrydeLe Gac and Duval (1980) and Reference Sunder and WuSunder and Wu (1989a). For con-venience they are referred to as LGD and SSW models in the following, and the same notation is used, when possible, for the two models.
Common bases for the models
In the absence of microcracking, the minimum creep rate is given by Norton-Hoff's law, expressed in a multiaxial form as
where 
is the deviatone part and 
the second invariant of the applied stress tensor 
. J2 is defined by
The parameter 
is temperature dependent and the value of exponent n is close to 3.
In the following we consider only deformation at a constant temperature. Also, we consider neither micro-cracking nor the dynamic recrystallization processes which occur after the minimum creep rate has been reached. In this respect the minimum creep rate is then seen as a steady creep rate.
The total strain 
is split into a pure elastic 
and a visco-plastic 
component as
By convention the transient creep strain 
is defined as
where 
is the creep strain resulting from Equation (1). The internal state variables of the models are the visco-olastic strain 
and the variables 
and r which describe, respectively, the kinematic and isotropic hardening.
The thermodynamic forces are the stress 
, the kinematic internal stress 
and the isotropic internal stress R associated with 
and r respectively. They are defined by the set of equations
where 
, is the Helmotz free energy potential.
If an expression for the dual dissipation potential, 
can be found, the evolution equations of the model are derived by applying the normality rule. Hence,
Clausius-Duhem's inequality is then automatically verified, provided that 
is positive, convex and null at the origin 
.
LGD model
The dissipation potential is written as the sum of a flow potential and a recovery term as
where 
if 
.
The free energy, Ψ, taken as
where Ψe denotes the component relative to the purely elastic (Hookean) deformation. The model parameters A, B, C, Η, Κ are positive constants (which may depend on the temperature T).
Since Norton's law (1) has to be satisfied at steady state, the following relations hold:
in which 
and 
are the steady values of 
and R.
Adopting the value n = 3, the multiaxial model derived from Equations (7) and (8) reduces, in the case of uniaxial compression, to
where 
are the tensor components relative to the direction of loading, and E is the Young's modulus.
SSW model
The dissipation potential is decomposed into a component 
which corresponds to the steady-state flow, and a transient flow potential, 
, as
The last expression implies that the isotropic stress R is strictly positive. Using the normality rule, it is shown that the internal variable 
is identical to the transient strain 
as defined by Equation (4).
The free energy function is
where Η is a positive constant and g is a function of r only. As a consequence of 
and Equation (12), the transient strain is totally recoverable.
The rate equation for R is taken independently as
Since J2() is strictly positive, the condition R > 0 is fulfilled as long as the initial value, R0 of R is positive.
Equation (13) is essentially the same as Equation (13) in Reference Sunder and WuSunder and Wu (1989b), on the obvious condition that the equivalent quantities σeq and εeq, defined in the original paper, are effectively tensor invariants (that reduce to |σ| and |ε| in the uniaxial case, instead of σ and ε, as stated by these authors). The uniaxial equations derived from Equations (11), (12) and (13) are, with n = 3,
The last equation does not match identically Equation (11) in Reference Sunder and WuSunder and Wu (1989a), which is given for arbitrary loading histories and in which the absolute values are dropped. Nevertheless the two formulations are equivalent if n = 3 and as long as R0 > 0 since 
.
Evaluation of the models
The LGD and SSW models have been tested by Reference Sunder and WuSunder and Wu (1990) against the master curves drawn by Reference Ashby and DuvalAshby and Duval (1985) from the dala of Reference JackaJacka 1984) obtained from constant stress experiments. The two models have a sufficient number of parameters to provide a good fit to these relatively simple and smooth master curves. The present analysis examines their ability to reproduce relatively complex loading cases.
Experimental data
Since published data concerning varying load experiments are scarce, a series of uniaxial compression tests was performed. Samples of granular isotropic ice were prepared by freezing a mixture of sieved crystals and deaerated-dcionized water. The ice density was 0.914 ± 0.005. Specimens of two average grain-sizes, 1 and 3.5 mm, were prepared. The mean grain-size was determined by counting the number of grains in a given area on photographs of thin sections cut from each specimen. In order to obtain cylindrical specimens approximately 120 mm long and 65 mm in diameter, the moulded samples were frozen on a circular steel platen, designed to fit the upper platen of the testing machine, then machined on a lathe. The specimens were stored at –10°C during 2 days before testing.
The tests were performed in a cold room at the same temperature with a lever operated machine. The axial force was measured by a load transducer mounted under the lower fixed steel platen. The axial strain was measured by two or three LVDTs mounted on the middle part of the specimens. The error in the strain measurements was within ±2 × 10 −5. The stress and strain measurements were stored in a data acquisition system.
Although the specimens were centered carefully, most tests exhibited a large difference between the strains measured for each specimen (i.e. during the same test, at the same time). This reveals a lack of uniformity in the strain measurements, which may have several causes, e.g. departure of the load from uniaxial, the testing machine being too compliant at high stresses; misalignment of the strain transducers; or small inhomogeneities in the specimens, increasing during the test. This discrepancy can be quantified by the ratio 
, where εmax and εmin are the maximum and minimum strains, measured for a specimen, at the end of the first loading phase. 
was less than 10% for only one half of the tests. Compared with this, the error associated with the measuring device (±2 × 10−5) was negligible, since the total strains were always higher than 3 × 10−3.
For comparison with model predictions, two tests were used consisting of constant stress loading followed by total unloading; during Test 2 two stress jumps and drops were applied during unloading. The corresponding measured
Fig. 1. Test 1 LVDΤ creep and recovery measurements; the loading conditions are shown by shaded area (1.65–0.07 M Pa)
Fig. 2. Test 2 LVDΤ measurements with stress jumps during unloading and mean curve used for model comparisons; the loading conditions are shown by shaded area (1.54, 0.05, 0.98, 0.05, 0.93, 0.05 M Pa)
strain curves are shown in Figures 1 and 2, along with the loading histories. Test 1 (grain-size 1 mm) was selected because of its very low ratio 
. Test 2 was taken as representative of “acceptable” tests, with 
, for specimens of average grain diameter 3.5 mm.
The model evaluations were made by using, for each test, the mean arithmetic strain curve calculated as the point-by-point average (at the same time) of the strains measured with the extensometers. No smoothing procedure was used in the calculation of this curve. The mean curve relative to Test 2 is shown in Figure 2; that of Test 1 was not drawn on Figure 1. for clarity.
Numerical solutions
Both models were applied using an implicit scheme for the integration of the sets of differential equations. An adaptive time step was used so as to ensure numerical Stability. The adjustment of model parameters was done using an iterative procedure which examined each parameter in turn and computed the best estimate so as to minimize the mean square deviation between the computed strain and the mean strain curve. This procedure does not require that the mean strain curve be smoothed before processing, since it constitutes a regression procedure by itself.
The values of parameters H and Κ (Equations (10)–(14)) displayed on the figures are scaled with the conventional Young's modulus E = 9500 M Fa (Reference SinhaSinha, 1979). The fluidity parameters 
A, B, C (Equations (1),(10),(14)) are scaled with 
where th stands for theoretical. This value was derived from our experimental data at –10°C and is a little less than that given by Reference JackaJacka (1984).
Evaluation of LGD model
When optimized to fit only the loading part of Test 1, the model exhibits a relaxation time that is too long for the
Fig. 3. Comparison of Le Gac and Duval's model with Test 1 data, a, computed strain curve fitted on the entire observed curve; b, computed recovered strain when the parameters are optimized to fit only the loading branch
recoverable part of the transient strain (see Fig. 3, curve (b)). This trend is removed when fitting the model over the whole curve, but then the steady state is reached too early at a strain near 4 × 10−3 (instead of 1% as is generally agreed) with a strain rate a little too high (see Fig. 3, curve (a)).
Fig. 4. Comparison of Le Gac and Duval's model with Test 2 data. Optimized strain curve (thick curve)
Figure 4 shows that the model fails to simulate the ice response to the stress jumps applied during the unloading phase of Test 2. The corresponding strain increments are lower than expected. Moreover since the model has to respond quickly to applied stress increments, its response on first loading is also very fast and does not give a good description of primary creep (as in Test 1).
This behaviour is explained by the fact that all the rapid strain variations following stress changes are controlled by the fluidity parameter, A, in Equation (10), which characterizes the mobility of the dislocations in the basal planes of the crystals.
Evaluation of SSW model
The ability of the model to describe the loading unloading cycle of Test 1 was evaluated by computing the optimized sets of parameters which fit either the whole cycle, or only the strain recovery branch. The results
Fig. 5. Comparison of Sunder and Wu's model with Test 1 data, a, computed strain nine, fitted on entire observed curve, along with its transient and steady-creep components; b, computed primary creep branch with parameters optimized to fit only the observed strain recovery
shown in Figure 5 indicate that this model cannot reproduce adequately both the loading and unloading parts of the strain curve. It can be argued that this result has no general value since it concerns a single experiment. Nevertheless, the shape shown in Figure 1 is always observed when granular ice is unloaded after the steady creep rate has been reached. The two simulations shown in Reference Sunder and WuSunder and Wu (1989a, Fig.8), concern tests in which the steady state was not reached when unloading occured
Fig. 6. Comparison of Sunder and Wu's model with Test 2 data. Optimized strain curve and transient and steady creep components
at a strain of 10−4, and seem reasonably good because the creep parameter 
was overestimated by a factor 3.
When tested against Test 2 (see Fig. 6), the model shows a good ability to reproduce the fast response to the increments of stress which were made during unloading. In contrast, the steady creep rate is too high and attained too quickly, at a strain of about 1.5 × 10−3, which renders the computed primary creep curve quite unrealistic.
For the SSW model the viscoplastic strain is decomposed as 
, and all the rapid variations in the total strain rate are governed primarily by 
. Because 
reaches its limiting value very soon on first loading, the steady state is achieved too quickly. In addition, since the observed recoverable strain is small, and the model assumes that 
is totally recoverable, the steady state is obtained at loo low a value of the total strain.
Alternative model
From this brief evaluation it appears that LGD and SSW models are unable to provide a correct description of both fast responses to stress jumps and primary creep.
In order to uncouple the first loading and long-term behaviour from the short-term response to fast changes in stress, we propose to split the viscoplastic strain, 
, into two components as
where 
is related to the kinematic hardening processes and provides the fast variations of strain, and 
is related io the isotropic hardening. The dissipation potential is then decomposed as
with 
if 
, and the free energy potential is written as
where 
and 
are conjugate internal variables relative to kinematic and isotropic hardening.
Using Equations (16) and (17) and applying the normality rule, leads to the multiaxial equations of the model. In the uniaxial case, with n = 3, they reduce to
Fig. 7. Comparison of the proposed model with Test 1 data. Optimized strain curve and strain components
where 
are the tensor components along the direction of loading.
Since the model has to be equivalent to Norton's law (1) at steady-state 
the following constraint on the fluidities can be derived
and the steady values of the internal stresses are
The physical interpretation which can be given is that the 
component, which includes all the kinematic effects, is produced only by a fraction of the total mobile dislocations. These dislocations are trapped in the networks and pile-ups capable of storing energy. The
Fig. 8. Comparison of the proposed model with Test 2 data. Optimized strain curve and strain components
rate at which this energy is stored and released is controlled by glide in the basal planes. A part of 
is unrecoverable as dislocations can escape (by climb) from the substructure traps. In that sense, the model for 
is the same as in LGD model, except that the isotropic hardening acting on the basal movement of the dislocations is neglected. The second component 
can he seen as the deformation resulting from the fraction of dislocations whose basal motion is impeded by local kinematic stresses and which escape directly by climb, taking no part in the anelastic processes. It is subject to isotropic hardening, which increases with the dislocation density.
Comparison of this model with experimental data from Test 1 and Test 2 is shown in Figures 7 and 8. Test 1 is very well reproduced by the model, with parameters optimized to fit only the loading part of the strain curve. For Test 2, the strain response to stress jumps during total unloading is quite satisfactory as well as the simulation of the primary creep curve. Clearly, better results should be achieved if 
and 
were described by a full LGD model, but with a major price to pay as regards to the great number of parameters to deal with.
Discussion
Contrary to the SSW model, the LGD and the proposed models do not involve explicitly the transient creep strain 
in their formulations. For these models, 
must be calculated as 
. In the LGD model, ε is the solution of Equations (10) and 
is given by Equation (9). In the present model the corresponding equations arc (18) and (19).
Reference SinhaSinha's (1979) equation gives the total strain as 
where 
is the (steady) creep strain resulting from Equation (1). The delayed elastic strain 
is thus formally identical to 
as defined by Equation (4). This component, as well as the anelastic moduli studied by Reference Gold, Traetteberg and WeaverGold and Traetteberg (1975), exhibits an exponential time dependence in the form Exp (–at1/3) (where a is a temperature-dependent parameter). In the absence of dislocation creep and recovery processes, the present model, with Ai = Β = 0 in Equations (18), does not involve the same time dependence, since it then reduces to 
and 
. Nevertheless this model was shown to be able to give acceptable responses to varying loads, at the scale of the total strains achieved during the tests 
and for loading durations higher than three hours. To conform to Reference SinhaSinha's (1979) or Reference Gold, Traetteberg and WeaverGold and Traetteberg's (1975) results, while keeping the framework of internal state variable models, more complexity must be added in the formulation of the kinematic strain component. The SSW model and Reference SinhaSinha's(1979) equation, describe the transient creep strain, or the delayed elastic strain, as totally recoverable. This is not the case for the LGD and the proposed models, in which recovery processes for kinematic hardening, and the existence of isotropic hardening, render part of the transient strain unrecoverable. This latter point seems to be in accordance with our observations (i.e. for total strains higher than 5 × 10−3).
Conclusion
An evaluation of the models of Reference Le Gac, Duval and TrydeLe Gac and Duval (1980) and Reference Sunder and WuSunder and Wu (1989a) has been presented. It was based on the simulation by the models of two uniaxial compression tests performed on isotropic granular ice under varying load at –10°C. Examination of the computed optimized strain curves showed that the two models fail to reproduce correctly both the fast variations of strain following stress jumps, and primary creep.
A new model, involving six parameters (instead offne for LGD model and four for SSW model), based on a decomposition of the viscoplastic strain into two components which account separately for the kinematic and isotropic hardenings, was shown to improve the quality of the simulations.
Much work remains to be done to obtain the reliable data required to confirm the validity of this model for uniaxial situations involving relatively high total strains, and to establish its correspondence with rheological models using delayed clastic and plastic terms.
Acknowledgements
This work was supported by the Centre National de la Recherche Scientifique, CNRS/SPI, France. The authors are particularly grateful to Dr P. Duval. Dr L.W. Gold and the anonymous reviewers for very helpful comments and discussion. Special thanks are due to Dr L. W. Gold and an anonymous reviewer for English corrections.
