A viscoelastic Rivlin-Ericksen material model applicable to glacier ice

We present a viscoelastic constitutive relation which describes transient creep of a modified second grade fluid enhanced with elastic properties of a solid. The material law describes a Rivlin-Ericksen material and is a generalization of existing material laws applied to study the viscoelastic properties of ice. The intention is to provide a formulation tailored to reproduce the viscoelastic behaviour of ice ranging from the instantaneous elastic response, to recoverable deformation, to viscous, stationary flow at the characteristic minimum creep rate associated with the deformation of polycrystalline ice. We numerically solve the problem of a slab of material shearing down a uniformly inclined plate. The equations are made dimensionless in a form in which elastic effects and/or the influence of higher order terms (i.e., strain accelerations) can be compared with viscous creep at the minimum creep rate by means of two dimensionless parameters. We discuss the resulting material behaviour and the features exhibited at different parameter combinations. Also, a viable range of the non-dimensional parameters is estimated in the scale analysis.


Introduction
Creep of ice is an inevitable issue when dealing with glaciers and terrestrial ice masses.From the pioneering work in the 1950s (Glen, 1952;Nye, 1953;Steinemann, 1958), the viscous constitutive relation for stationary creep of ice has emerged, which, in the glaciology community, is referred to as Glen's flow law.The Glen flow law is a generalized Newtonian material model (e.g.Crochet et al., 1984) with powerlaw viscosity and is widely used when modeling the flow of glaciers, ice streams and ice sheets at any spatial and tempo-Correspondence to: P. Riesen (riesen@vaw.baug.ethz.ch)ral scale, in spite of its constraint to only predict a stationary stress/strain rate relation at the minimum creep rate observed in laboratory creep experiments.However, the laboratory experiments on the creep of ice conducted in the past consistently revealed non-stationary creep of ice (e.g.Jellinek and Brill, 1956;Mellor and Cole, 1982;Jacka, 1984).
Several attempts were made to describe the transient creep of polycrystalline ice observed in laboratory creep experiments by viscoelastic constitutive equations, relating strain rates, strain, stress and time in different ways.Nonlinear, time-dependent constitutive relations describing creep behaviour observed in uni-axial compression tests at constant load were given by Sinha (1978) and Le Gac and Duval (1980), and reviewed in Ashby and Duval (1985).The work of Szyszkowski and Glockner (1985) considered a nonlinear constitutive equation based on spring and dash-pot elements.A description of the transient strain rate as a nonlinear power-law function of stress and strain was proposed by Azizi (1989).Shyam-Sunder and Wu (1989a,b) published a differential flow model, and in Shyam-Sunder and Wu (1990) it was compared with the models of Sinha, and Le Gac and Duval.Later on, Meyssonnier and Goubert (1994) took up the models of Le Gac and Duval, and Shyam-Sunder and Wu and proposed some modifications.While the early attempts of the aforementioned models describe only uni-axial creep responses, and lack an obvious generalization to multiaxial creep states (e.g.Sinha), the more recent models of Le Gac and Duval; Shyam-Sunder and Wu, and Meyssonnier and Goubert, include one or more state variables which must be modeled by additional differential equations.Morland (1979) and Morland and Spring (1981) proposed constitutive equations of rate type (Lockett, 1972;Hutter, 1983) to describe the viscoelastic responses of isotropic polycrystalline ice.They considered constitutive relations which relate stress and stress rates to either strain rates and strain accelerations, which is a fluid type model (Morland and Spring, 1981), or else to strain and strain rates, i.e. a solid type model P. Riesen et al.: A viscoelastic Rivlin-Ericksen material model applicable to glacier ice (Spring and Morland, 1982).These models are capable of reproducing idealized decelerating (primary), stationary (secondary) and accelerating (tertiary) creep responses, but were not applied to real creep data.
If stress rates are excluded from the constitutive relation and the stress is considered only a function of the strain rates, strain accelerations, and possibly higher order derivatives, the relation then describes a fluid of differential type (Rivlin and Ericksen, 1955).These fluids are also known as fluids of grade N (e.g.Joseph, 1990), or, order fluids (Owens and Phillips, 2002).The incompressible second grade fluid model (fluid of grade 2) was first applied to creep experiments and ice mechanics by McTigue et al. (1985), who investigated the relevance of normal stress differences in shearing flow.However, the second grade fluid model has no strain rate-dependent viscosity.Man and Sun (1987) thereupon postulated a modified second grade fluid with power-law viscosity.The relation of Man and Sun describes primary and secondary creep, and the material behaviour asymptotes to a Glen-like power-law fluid model for vanishing second order terms.For ice, Sun estimated the phenomenological coefficients for the modified second grade fluid model with powerlaw viscosity from triaxial laboratory creep experiments.
In recent field observations on an Alpine glacier (Gornergletscher, Switzerland), repeated near-reversal of flow, accompanied by reversed displacement direction was observed (Sugiyama et al., 2007(Sugiyama et al., , 2008)).The change of motion takes place within a few days during the periodical drainage of a supra-glacial lake, and is possibly related to the unloading and stress redistribution during the rapid drainage.The question was raised whether the retrograde movement may be attributed to viscoelastic recovery properties of the ice.To elaborate on this hypothesis, an appropriate viscoelastic constitutive model for glacier ice is needed, preferably applicable to multi-axial deformations.No attempt has been made so far to take into account transient, recoverable deformation effects when an external load (e.g. a forming and draining lake) is applied to and removed from a glacier.This was our motivation to construct a simple constitutive formulation with a material response varying smoothly in between the limits of ice behaving as a viscous fluid on long time scales, and ice behaving as an elastic solid on short time scales.We here propose a constitutive model able to predict instantaneous elastic strain followed by recoverable, transient strain fading into a steady creep response associated with the stationary minimum creep rate of ice, based on a further generalization of the modified second grade fluid of Man and Sun (1987).A tertiary response with increasing strain rate after the minimum strain rate is not considered.
Ultimately, we target to corroborate whether glacier ice allows for deformations such as those observed.By now, we develop the constitutive relation aimed for, and investigate it in a numerical example.The possible range of parameters is collected and the characteristics which can be exhibited by the rheological model are illustrated.

Constitutive relations
In our description we use the Eulerian notation.The kinematic measures used are summarized in Table 1, following Lockett (1972) and Hutter and Jöhnk (2004).Before we discuss the proposed generalized constitutive relation, we now quote the relations for the already mentioned material laws.In general, the ice is assumed an incompressible material, which introduces the pressure as an independent field, absorbing any isotropic stress contribution.

Glen's flow law
In Glen's flow law, the Cauchy stress tensor t is given by the relation Here, p is the pressure due to the incompressibility constraint, I is the identity tensor, η(D) the strain rate-dependent power-law viscosity, B a constant, and n the power-law exponent.For n = 1, the material model becomes non-Newtonian.This relation is valid for all times; obviously it can merely describe stationary creep.For the monotonic secondary creep regime of ice, a value of n = 3 is commonly used.The parameter B is in the range of 1.2 to 2.9 MPa d 1 3 for ice at temperatures between 0 • C and −10 • C, as recommended by Paterson (1994).

The modified second order fluid model
The modified second order fluid model with power-law viscosity (MSOFM) was introduced by Man and Sun (1987) as where A (1,2) are the first and second Rivlin-Ericksen tensors, describing the current strain rate and strain accelerations (Rivlin and Ericksen, 1955).The first Rivlin-Ericksen tensor in terms of the spatial velocity gradient is A (1) = L+L T = 2D.The strain acceleration tensor A (2) , and higher order tensors are constructed via the recurrence relation (viii) in Table 1.The coefficients α 1,2 are termed normal stress coefficients.The viscosity η(A (1) ) follows a power-law relation of the form (3) In Sun (1987), the MSOFM was proposed as an improvement on Glen's flow law to include non-stationary creep.Glen's flow law is contained in Eq. (2) as the asymptotic limit after transient creep has died out.Or else, with use of (vii) and (viii) (Table 1), Eq. ( 1) is recovered when α 1 = α 2 = 0 in Eq. ( 2), m = (1 − n)/n, and µ = 2 −1/n B, is substituted in Eq. ( 3).If m = 0, the MSOFM becomes the second order The trace of a measure, i.e. tr(D) = D ii fluid model.Clearly, Eq. ( 2) is a generalization of the second order fluid model.A good summary on related generalizations of the second order model and the MSOFM is given in Massoudi and Vaidya (2008).Sun (1987) performed the exploitation of the second law of thermodynamics (Clausius-Duhem inequality) for the MSOFM and deduced the restrictions on the stress function and its parameters µ, α 1 and α 2 .There is much controversy on the second order fluid models.Experimental observation and mathematical analyzes on the fluid model and its stability properties do not share the same consequences on the sign of the normal stress coefficient α 1 (see e.g.Dunn and Fosdick, 1974;Joseph, 1976;Müller and Wilmanski, 1986;Joseph, 1990;Rajagopal and Srinivasa, 2008).A review on this topic is given by Dunn and Rajagopal (1995).Here, we follow the work of Sun (1987), who used the restrictions (i) α 1 + α 2 = 0, and (ii) m = −2/3, corresponding to n = 3 in Glen's flow law, Eq. (1).Sun determined mean values of µ = 2.41 MPa d 1 3 and α 1 = 161 MPa d 2 from fitting the model to the data of triaxial (McTigue et al., 1985), and pressure-meter (Kjartason, 1986) creep experiments.

The elastic modified second order isotropic material model (EMSOIM)
In the past, the role of primary creep and elastic effects on glacier flow and observed flow anomalies on a scale of hours to a few days has scarcely been investigated.To strike a new path in this direction, we adopt the MSOFM, in which the ice is able to reproduce both primary and secondary creep effects with good agreement on the creep experiments analyzed by Sun (1987).However, we further require the material to exhibit two contrasting properties: (1) viscous stationary creep of a fluid on long time scales, and (2) elasticity of a solid when the time scale under consideration becomes short, allowing elastic strain jumps and reversible creep.The material law needs to include properties of an isotropic elas-tic solid.We propose to extend the constitutive form of the MSOFM by an explicit dependence on the Finger strain tensor (Table 1), postulating as a frame-invariant functional relation for the stress.We call the corresponding material an elastic modified second order isotropic material (EMSOIM).The functional form ( 4) is a further generalization of the MSOFM, and belongs to the class of Rivlin-Ericksen materials (Rivlin and Ericksen, 1955).
We now restrict the constitutive model ( 4) by ad-hoc assumptions to make the mathematical proof of the thermodynamic behaviour performed by Sun (1987) applicable to the EMSOIM.In this way, we preserve the essential properties, namely (i) inclusion of elasticity effects, which, paired with the viscous effects, allow for relaxation phenomena, and (ii) use of the MSOFM concept to account for the primary and secondary creep regimes, in the context of an extended Glen flow law.
For an isothermal process with a body at uniform temperature the thermodynamic analysis of Sun (1987) is recovered for the EMSOIM (Eq.4) if the following postulates hold: 1.The Cauchy stress tensor t can be additively composed as where tD is a dissipative stress component that does not depend on e, and tE is an elastic stress contribution.
2. The functional dependence of the Helmholtz free energy may be additively decomposed as 3. The component ψ2 of the Helmholtz free energy is a convex function of its argument A (1) .
We refrain from carrying out the complete thermodynamic analysis here as the problem which we consider in the following is an isothermal process and does not require an energy balance to be solved.An admissible form of the elastic stress contribution t E for an incompressible material with elastic deformation limited to shearing is the relation where is the deviatoric strain tensor, and β is a variable shear modulus of the form with initial rigidity β 0 = 2G, where G = 3500 MPa corresponds to the shear modulus of ice (Schulson and Duval, 2009).The constant c ≥ 0 is referred to as "fading elasticity factor".The purpose of the exponential dependence of the shear modulus on strain evolution is the ability to destroy elasticity on long time scales.It introduces an exponentially fading strength of elasticity with increasing strain accumulation.This is analogous to an exponentially fading memory of the material's elastic properties with increasing deformation.This should not be confused with the concept of materials with fading memory (Coleman and Noll, 1960).Thus, in the EMSOIM, the Cauchy stress tensor t takes the form where α = α 1 = −α 2 was used.This representation is a fairly general form of a material law containing the modified second order fluid model by Man and Sun (1987), the Glen flow law and a Hooke-type elasticity relation with vanishing bulk modulus for an incompressible isotropic material of the class of materials of differential type (Rivlin and Ericksen, 1955).A similar, even more general constitutive relation is discussed by Zhou (1991).
In Table 2 we summarize values of the parameters m, µ, α and β 0 which are available from the literature and experimental data.The values of n and B in Glen's flow law, given in Sect.2.1, are expressed in terms of m and µ.

Unidirectional flow of the EMSOIM
We solve the balance equations for incompressible, isothermal Stokes (creeping) flow with the EMSOIM law.The governing equations are Equation ( 12) is the momentum balance, where ρ is the (constant) material density and f an external force.Equation ( 13) is the mass balance, which requires the velocity v to be solenoidal.The remaining equations describe the constitutive relation of the EMSOIM.To solve the system forward in time, the set of equations must be complemented by an evolution equation for the Finger strain tensor e from the present time t.Differentiation of e, using (vi) from Table 1 which follows as an identity.
As a benchmark example, we consider one-directional shearing flow of a slab of material down a uniformly inclined plate (UIP), as illustrated in Fig. 1.Congruent to the two-dimensional spatial Cartesian coordinate system (x, z) depicted in Fig. 1, we define the particular reference configuration X = (X, Z).The plate is considered infinite in the downstream direction (X, x), and the slab thickness is monotonic.The flow is driven by a gravitational force f = (gsinφ,−gcosφ) where φ is the inclination angle of the plate.We seek the velocity profile v = (v(z),0) which varies with the thickness of the slab.In this rectilinear flow problem, the following relations exist between present and reference configuration: where u is a spatial displacement in the flow direction.The deformation gradient is then We note that ∂u ∂Z = ∂u ∂z ∂z ∂Z = ∂u ∂z , so that the Finger strain tensor is given by It may appear appropriate to use a geometrically linearized strain tensor, i.e., neglecting terms less than O((∂u/∂z) 2 ).However, this would render the use of a variable modulus inconsistent, as tr(e 2 ) ∝ (∂u/∂z) 2 .Thus, the geometrical linearization implies explicitly that one only moves little away from the reference configuration; the material's elastic properties would remain unchanged in that range of deformation.
Since we require degradation of material elasticity, terms less than O((∂u/∂z) 2 ) should be retained in the formulation for consistency.
The first Rivlin-Ericksen tensor is determined from the spatial velocity gradient as The evaluation of Eq. ( 17) brings up the relation between velocity v and displacement u as which, in this unidirectional case, is simply v = ∂u/∂t.We keep v and u separate and insert Eq. ( 14), together with Eqs. ( 20) and ( 21) into the momentum equation ( 12).The stress divergence yields two equations as with the velocity v, displacement u, and pressure p as unknowns.
Note, the pressure equation ( 24) is decoupled from the flow equation ( 23).Equation ( 24) is interesting as it obviously describes deviations from a hydrostatic stage (linear pressure variation with slab thickness) which arise from second order effects and elasticity.These normal stress contributions due to second order effects (term associated with α) and elasticity (term associated with β 0 ) in Eq. ( 24) carry opposite signs.If we would have used a geometrically linearized strain tensor, elastic normal stresses, i.e., the second term in Eq. ( 24) would be absent.
In the following, we will only be concerned with the solution of the flow problem (Eq.23).
We now non-dimensionalize problem (Eq.23) by replacing the relevant fields on the basis of a characteristic time scale as where each bracketed term represents a characteristic scale of the respective field.For Eq. ( 23), upon inserting Eq. ( 25) and dividing by the load ρ g and rearranging, we obtain with the -coefficients The coefficients 1−3 contain ρ g to the first power.This means that only two -products are independent, namely H = 2 / 1 and K = 3 / 1 , evaluated as The non-dimensional number H measures the significance of strain accelerations, while K is the initial rigidity modulus.
v(z = 0, t) = 0, ( 31) where Eq. ( 31) is the Dirichlet boundary condition, assuming that the slab adheres to the plate, and Eq. ( 32) describes a set of initial conditions.
We have now reduced the set of coefficients to the parameters H and K as measures on viscoelastic, transient effects (strain acceleration and/or elasticity) relative to the purely viscous power-law material, i.e., the first term of Eq. ( 29).

Numerical implementation
The field equations are numerically implemented using the DOLFIN/FFC finite element software (Kirby and Logg, 2006;Logg and Wells, 2010).We solve the flow problem of Eqs. ( 29), (30) together with the homogeneous boundary condition (31) and the initial conditions (32) in a mixed problem.The weak form of the discrete Galerkin formulation is to find for all admissible (b h , w h ) ∈ U h × V h .Here, F is a bilinear form, where (•,•) is an L 2 ( ) inner product for scalars defined with respect to the partition of the bounded domain of R 1 into finite elements, and the subscript h is a discretization parameter.The finite element spaces U h and V h are appropriate spaces of square integrable basis functions with derivatives also being square integrable.We use continuous Galerkin elements with Lagrange polynomials of second order.The time derivatives are discretized using the implicit backward Euler scheme.We apply Newton's method to solve the resulting system of nonlinear algebraic equations.
Note that Eq. ( 33) results in a system of nonlinear differential-algebraic equations (DAEs) of the form G(t, u, v, ∂v/∂t) = 0.It contains the initial value problem of satisfying a consistent set of initial conditions for the N-dimensional vectors of unknowns u(t 0 ) = u 0 , ∂u/∂t| t 0 = v(t 0 ) = v 0 , and possibly ∂v/∂t| t 0 (e.g. Brown et al., 1998).Initially, the creep rate of ice is larger than in the stationary creep regime and it is decelerating (primary creep).We thus require v 0 > v steady .To compute a consistent v 0 , we solve one time step with H = 0 using a very small step size.In this way, we obtain a valid steady solution (u,v).The velocity was then multiplied by a constant a, i.e. (u 0 ,v 0 ) = (u,a v).We use a = 2.5 and then start the actual computation with H specified, and (u 0 ,v 0 ) as initial guesses.If K = 0, we reset the displacement u 0 = 0 as no initial instantaneous elastic displacement occurs.If K > 0, u 0 also served as initial guess.

Creep under step function load
The UIP problem may be interpreted in analogy to a shear creep experiment.At time t < 0, the load is zero and the material has been at rest for a long time.For t ≥ 0, we assume the load to be constant, i.e., equal to sinφ.However, as we are also interested in the unloading phase, we artificially remove the (gravitational) load after some time, so we define the modified load sinφ(t) with where φ c = 12 • , H (•) is the Heaviside function, and t r is a portion of the experiment run-time (60% of the total runtime).Of the solutions, we display velocity (creep rate) and displacement (creep) at the surface of the slab as functions of time.In Figs. 2 to 4 we elucidate results for different parameter combinations.
If second order (creep accelerations) and elasticity effects are absent, i.e.H = K = 0, the purely viscous flow problem, equivalent to the Glen power-law is solved.In that case, the material simply shears down the plate with steady creep rate (velocity).As the load is removed at t = 0.9, the creep rate instantaneously drops to zero.This case is depicted by the black solid curve in Fig. 2a.In Fig. 2b, the corresponding creep curve is depicted.The displacement increases linearly with time with a slope corresponding to the creep rate (Fig. 2a).At unloading (t = 0.9), the creep curve remains at a constant level of permanent creep experienced so far.
For increasing (decreasing) power-law exponent m, the constant creep rate increases (decreases) proportionally.However, the creep rate is fixed for a constant load, and no transient behaviour occurs.In all the following computations, we used m = −2/3 (n = 3).
In Fig. 2a, b, we have 0 ≤ H ≤ 10, while K = 0. Here, the material initially deforms rapidly with high creep rate.For H ≤ 1.0, the creep rate decays rapidly and asymptotically reaches the steady creep rate (the solution with H = 0,  Fig. 2a, black solid).If H = 5 or 10, the creep rate decelerates significantly, but does not reach stationary creep until t = 0.9.When the load is removed at that time, the creep rates do not drop instantaneously to zero, but further decay for all H > 0. For small H, the creep curves in Fig. 2b are composed of four sections; an initial interval of decelerating (primary) creep followed by monotonically increasing, stationary creep, and a recurring interval of decelerating creep, which then goes over into u = const.However, as H is increased, the alternation of these different creep intervals is blurred and the displacement function transforms into a single interval of almost permanently decelerating (primary) creep (see H = 10, Fig. 2b).Note that the decelerating creep at t > 0.9 is particularly interesting since the material now exhibits, after removal of the load, continued decelerating creep.This situation indicates a retardation of the material's response to the removal of the load; with increasing H the creep rate decays increasingly slower.The creep undergone by the material is not recoverable, so the material is a viscous fluid with the ability to experience transient creep, according to the MSOFM.
In Fig. 2c, we now impose elastic properties with K > 0 and set H = 0.For any K the material immediately starts to creep with initial velocity as the solution with K = 0.All creep rates then decay asymptotically to zero; the time it takes to reach zero depends on K.At the removal of the load (t = 0.9), the creep rates instantaneously jump to negative values, which indicates that creep starts to recover.In Fig. 2d, the increase of displacements decelerates with increasing time t and eventually the individual displacement solution approaches an asymptotic limit (e.g.K = 10 3 in Fig. 2d).At unloading time t = 0.9, the displacements start to decrease again quasi-exponentially and re-approach zero.Thus, for K > 0, the material becomes elastic and rigid, which prevents permanent creep but allows complete recovery of the displacement.
In Fig. 3, we display the material behaviour when elastic and second order effects are both activated with nonzero H and K.We show solutions for H = 2.0 and varying K (Fig. 3a, b), and for K = 10 2 with variable H (Fig. 3c,  d).When H > 0, the slab starts to creep with initially high creep rate.If K is very small, i.e., K = 1.0, the decay of the creep rate slows down and becomes almost constant (Fig. 3a).At removal of the load (t = 0.9), the decay of the creep rate speeds up again and approaches zero fast.If K = 10 2 , the creep rate initially decays rapidly but slows www.nonlin-processes-geophys.net/17/673/2010/ Nonlin.Processes Geophys., 17, 673-684, 2010 down considerably after t ∼ 0.3.If K = 10 3 , or larger, the creep rates decay very rapidly and even drop below zero and then increase again.So, the creep rates begin to oscillate and fade to zero afterwards.If the load is removed, the creep rates drop, oscillate and fade again to zero.Thus, since K 0 and H = 0, the material is strongly elastic and able to respond to creep accelerations; the material behaviour becomes resilient.The corresponding creep curves show decelerating creep reaching a maximum at increasingly earlier time, at which time recovery is also activated (Fig. 3b).In Fig. 3c, we observe how the creep rate decays increasingly less rapidly with increasing H.The material response gets strongly delayed.If H is small, the creep rate quickly decays asymptotically to zero.At unloading (t = 0.9), it jumps to negative values and again asymptotically returns to zero.If H = 1 or 5, the creep rate first decreases fast and then the decay slows down rather quickly.For large H ≥ 10, the creep rate decreases slowly and unloading at t = 0.9 has almost no effect.In this case, for increasing H and fixed K, the material experiences increasingly more creep until at unloading, the creep decreases (recovers) again (Fig. 3d).

Fading elasticity
As already mentioned, when K > 0 (Figs.2c, d and 3), the material is essentially a viscoelastic solid which can only experience limited, though fully recoverable creep.It is obvious that these properties of a viscoelastic solid dominate the material behaviour for all times t.However, viscous creep of a fluid should be the dominating material behaviour for t 0, thus with increasing time t, the material needs to forget its elastic properties.This is achieved by adjusting the fading elasticity factor c. In our creep experiment, if K > 0 and the load sinφ(t) is applied, the displacement approaches some asymptotic limit for increasing time.In that case, i.e., t 0, the last term of Eq. ( 29), associated with K, becomes dominant at t 0, while due to v → 0 the other terms diminish.Setting c ≥ 0 exponentially attenuates the growth of that term with increasing displacement (i.e., strain).This is equivalent to an exponential decay of the initial dimensionless modulus K with increasing displacement/strain.Note that this results in a temporal response of the material behaviour, however it is not an explicit time-dependent response, as the change in the phenomenological parameter is linked to the deformation.This response can be physically  interpreted as a change in resistance of the material, which e.g.Ashby and Duval (1985) and Castelnau et al. (2008) associate with a change in the internal stress field of the material.
We display the associated material behaviour in Fig. 4. If second order effects are left aside (H = 0) we observe the material creep rate to decay towards zero and then increase again (Fig. 4a).The larger c is, the earlier the decay of the creep rate will be interrupted and the creep rate increases again, approaching some constant nonzero value.For such cases, the creep decelerates and then accelerates again, increasing gradually with time.When the load is removed at t = 0.9, some recovery of creep takes place where the displacements decrease and approach some constant value again, corresponding to the amount of permanent viscous creep acquired (Fig. 4b).The larger c is, the more permanent viscous creep and the less recovery of creep occurs.If second order effects are taken into account, e.g.H = 3, primary creep is activated again and oscillating creep rates (re-)appear.The increase of c in such cases dampens the oscillations and prevents the creep rate from decaying to or below zero, maintaining increasing viscous creep with increasing c. for the MSOFM, under the conditions of m = −2/3, (n = 3 in Glen's flow law) and α 1 + α 2 = 0 (thermodynamic constraint).The resulting mean estimates of Sun were listed in Table 2.We now use m and α as determined by Man and Sun (1987, see Table 2) and insert them into the definition of H (first equation of Eq. 28).In Fig. 5, we plot the variation of H as a function of time scale.The limits of the abscissa encompass the approximate range over which the experiments of Mellor and Cole (1982) and Jacka (1984) lasted.The time scales shaded in light grey ranges from 1 to 10 days, and corresponds to the durations of the creep experiments analyzed by Sun (1987), from which he estimated the parameters µ Eq. ( 28) 2 with 1.4 ¡ £ ¡ 2.4 MPad 1/3 , ¤ 0 = 7000 MPa Fig. 5. Change of the dimensionless numbers H and K of Eq. ( 28) as functions of the time scale and the parameters from Table 2. and α.The value of µ deduced by Sun is in the range of µ expected from the Glen flow law (Paterson, 1994).Considering a variation of µ according to the range given in Table 2, the change of H with increasing time scale and fixed α = 161 MPa d 2 is given by the dark grey bar in Fig. 5a.The range of expected H is the intersection of the grey bar with the time scale interval shaded in light grey, which is 1.5 ≤ H ≤ 115.For a time scale of 5 d, H is ≈5.In Fig. 5b, the change of K with increasing time scale, according to the second equation of ( 28), is shown.Here, the dark grey bar indicates the value of K for β 0 = 7000 MPa, and 1.4 ≤ µ ≤ 2.4 MPa d 1 3 .Expected values of K for a time scale between 1 to 10 d lie in the interval of 2.9 × 10 3 ≤ K ≤ 2 × 10 4 , with K ≈ 8 × 10 3 at [T] = 5 d.
We note that there is a considerable variation of H of two orders of magnitude with increasing [T], whereas K changes only one order of magnitude across the time scales of 1 to 10 d.In this range, the absolute value of K is about two orders larger than that of H at a given time scale.Thus, the solid elastic properties are the prominent feature of the ice at any time scale.Since the value of H shows strong variation across [T], the behaviour of a strongly elastic material (large K) is expected to vary quite significantly, depending on H.The occurrence of oscillating creep (Sect.5) is a striking indication to this.

Relevance of acceleration effects (H > 0)
As pointed out, the magnitude of H obviously varies strongly, depending on the scales considered.The interpretation is that the strain acceleration term in the constitutive relation can alter the material behaviour considerably, as demonstrated by the appearance of oscillating creep rates with increasing H and large K.The oscillating creep questions consideration of second order effects (i.e., strain accelerations) in a creeping (Stokes) flow problem.Nevertheless, if no creep accelerations are considered, i.e., α (H) is zero, it is not pos-sible to reproduce primary creep rates which are up to two orders larger than the steady secondary creep rate (e.g.Jacka, 1984;Castelnau et al., 2008).For K > 0, but H = 0, there is no solution with creep rates larger than the viscous, stationary creep rate (K = 0).The strain accelerations with H > 0 are necessary to capture the primary creep regime with the deceleration of the creep rate.However, in the EMSOIM constitutive equation, the significant decrease of strain rate in the primary creep regime will also be influenced by the evolving elasticity of the material and not only by strain accelerations, as in the MSOFM.Presumably, the actual value of α as a material parameter in the EMSOIM is smaller than it can be expected on the basis of the MSOFM.Thus, strain accelerations and decay of elasticity in the EMSOIM should be designed in such way that the interference does not result in oscillating creep behaviour.

Final remarks
We applied a viscoelastic material law in a simple unidirectional flow problem and studied the influences of the various material parameters.The EMSOIM relates stress to strain and its derivatives up to second order, i.e. strain rates and strain accelerations.As the EMSOIM is a generalized material law incorporating the Glen flow law, the MSOFM and a nonlinear elasticity relation, depending on the choice of parameters, it reproduces material behaviour with respect to either a single or multiple of the incorporated constitutive relations.
For unidirectional flow considered here, the set of field equations was substantially reduced, however, the constitutive model was capable of producing complex material responses.Such characteristics, including total or partial recovery of deformation and enhanced viscous deformation with non-stationary creep rates, may be possibly encountered in the observations on Gornergletscher.Such a multi-dimensional situation makes the problem much more challenging and besides the difficulties associated with the viscoelastic theory and material behaviour, numerically dealing with the problem in the context of finite elements becomes more involved.Nevertheless, we target on considering further applications of the EMSOIM constitutive model and the investigation of viscoelastic properties of glacier ice.

Fig. 1 .
Fig.1.Problem geometry of one-directional shearing flow of a slab of material down an uniformly inclined plate (UIP).The coordinate systems of the reference (X) and present (x) configurations are oriented concordantly, inclined by some angle φ relative to the gravitational force f acting in the vertical.

P.
Riesen et al.: A viscoelastic Rivlin-Ericksen material model applicable to glacier ice The non-dimensionalization emphasizes the parameters' relationship to transient effects.Dropping the superscript bars, the non-dimensional form of the initial/boundary value problem takes the form

Fig. 2 .Fig. 2 .
Fig. 2. (a, c) Velocity (creep rate) v at the slab surface as function of time.(b, d) Displacement (creep) u at the slab surface as function of time.In (a, b), solutions of v and u are displayed for different values of the dimensionless number H while K = 0.In (c, d), solutions of v and u are displayed with H = 0 and varying K.Other parameters are as displayed in the plot headers.All fields (u, v, t) are non-dimensional.

680P.Fig. 3 .Fig. 3 .
Fig. 3. (a, c) Velocity (creep rate) v at the slab surface as function of time.(b, d) Displacement (creep) u at the slab surface as function of time.In (a, b) solutions of v and u are displayed for H = 2 and variable K, while in (c, d) solutions of v and u are shown for K = 10 2 and different values of H.Other parameters are as displayed in the plot headers.All fields (u, v, t) are non-dimensional.

Fig. 4 .Fig. 4 .
Fig. 4. (a, c) Velocity (creep rate) v at the slab surface as a function of time.(b, d) Displacement (creep) u at the slab surface as a function of time.In (a, b) solutions of v and u are displayed for different values of the memory damping factor c H = 0.In (c, d) solutions of v and u are shown for various values of c and H = 3.0.Other parameters are as displayed in the plot headers.All fields (u, v, t) are non-dimensional.
ranges of H and K The EMSOIM includes three relevant material parameters, i.e., µ, α, and β 0 .The parameters µ and α have been determined based on laboratory experiments and creep function fitting.McTigue et al. performed the parameter identifications of µ, α 1 , and α 2 for a second order fluid model with m = 0. Sun (1987) re-fitted the creep data of McTigue et al.

Table 1 .
Definition of kinematic tensor quantities in Eulerian notation.In index notation, the tensor quantity is non-bold and indexed (i.e.L =L ij ).The material derivative () is given as∂()∂t + v j • () ,j , where v =v j is the velocity.Note that for (viii), we have A (0) = I, and A (1) = 2D.

Table 2 .
Values for the parameters of the Glen flow law and MSOFM, available from the referenced literature, expressed for the generalized EMSOIM constitutive relation.(*) The value for the modulus β 0 corresponds to 2G, where G = 3500 MPa is the shear modulus of ice