Introduction

Reference GlenGlen’s law (1955), which represents ice as a power-law fluid, has proven to be a useful and effective tool for modeling the shearing behavior of ice in simple flow configurations. However, Reference GlenGlen (1958) himself explicitly recognized the limited nature of the power law, and discussed some of the implications of more general constitutive equations. In particular, he considered an expression for the rate of deformation as a function of stress for an isotropic and incompressible fluid, specialized to exclude dependence on the spherical part of the stress. Glen’s law is a further specialization of the resulting equation that neglects terms of second order in the stress. Reference GlenGlen (1958) showed that the power law cannot consistently represent results from both uniaxial compression and direct shear tests, nor can it account for measurements made under simultaneous shear and compression achieved in torsion of an annular sample. The higher-order theory, however, is able to reconcile these observations.

This paper considers the application of the “second-order fluid” model to polycrystalline ice. Reference MorlandMorland (1979, Reference Morland and Spring1981) has discussed this possibility for a more general constitutive equation, of which the second-order fluid is a special case. The second-order fluid has close affinities to the general flow law treated by Reference GlenGlen (1958), but there are several fundamental differences. First, the second-order fluid arises in an expansion of stress as a function of the deformation-rate, while glaciologists often work with the inverse of such a relation. Non-linear constitutive equations often cannot be inverted, so that direct comparisons between these different approaches may be difficult. For example, the procedure discussed by Reference GlenGlen (1958), when applied to a relationship for stress as a function of deformation rate, leads to the definition of the “Reiner-Rivlin fluid” (see, e.g. Reference Bird, Bird, Armstrong and HassagerBird and others, 1977), which exhibits behavior significantly different from that of the material defined by Glen’s general expression. A second major difference between the second-order fluid and earlier formulations is that it involves an expansion not only in higher powers of the deformation rate, but also in higher order derivatives of the deformation rate. This extension was critical in developing a constitutive model capable of representing the unequal normal stress differences commonly observed in non-linear fluids in shearing flows, a phenomenon central to the discussion to follow. Finally, the coefficients in the general flow law posed by Reference GlenGlen (1958) are allowed to be functions of the invariants of the stress. Similarly, the coefficients in the constitutive equations considered by Reference MorlandMorland (1979, Reference Morland and Spring1981) are allowed to depend upon the invariants of the rates of deformation. The coefficients in the expansion leading to the second-order fluid are taken to be material constants.

The last feature of the second-order fluid model bears elaboration in the present context. An important consequence of this assumption is that the relationship between shear stress and rate of shearing deformation for the second-order fluid is identical to that for a linearly viscous fluid. Thus, the second-order fluid is in no way presented here as a replacement for, or generalization of, Glen’s law. Clearly, it fails to represent the non-linear rate dependence of ice in shear. However, just as the power law is a useful artifice for capturing non-linear shearing behavior, the second-order fluid is likewise a convenient idealization for representing certain other classes of non-linear phenomena. Prominent among these are normal stress effects.

The term “normal stress effects” refers here to various measurable kinematic and dynamic phenomena that are a consequence of rate-dependent normal stresses associated with shearing deformations. In particular, materials such as the second-order fluid, when subjected to simple shear, exhibit non-hydrostatic normal stresses perpendicular to the plane of shear and in the direction of the flow that are quadratic in the shear rate. Similarly, the material treated by Reference GlenGlen (1958), when subjected to a simple shear stress, exhibits a non-zero dilation rate normal to the plane of shear. It is important to emphasize that the term “normal stress effect,” in this context, is completely unrelated to the question of whether the rate of shear deformation of ice is affected by the spherical part of the stress. Indeed, as noted previously, the shear rate of the second-order fluid depends only on the shear stress.

Both theoretical considerations and experimental evidence indicate that normal stress effects are among the first manifestations of non-linearity in many fluids. Power-law constitutive models cannot properly represent normal stresses, and thus cannot be used to predict an entire class of important non-linear phenomenology. The purpose of this paper is to study the possibility of normal stress effects in ice. The second-order fluid, which is an idealization commonly used by polymer rheologists, is adopted as a tentative model. The three coefficients appearing in this model are determined from data for polycrystalline ice in biaxial creep. An excellent representation of primary and secondary creep is obtained. Normal stress effects are illustrated by examining the free-surface depression predicted for open-channel flow. Using the material properties determined for ice, it appears that such a depression could be measurable on glaciers on sufficiently steep slopes. Furthermore, the principal stresses at the surface of the channel are changed in both magnitude and direction, which could influence the development of crevasses.

The Second-Order Fluid

A quite general model for material behavior is the simple material. This material is defined by the constitutive assumption that the stress is a functional of the history of the strain (Reference NollNoll, 1958)Footnote ^* , and it includes as special or asymptotic cases many of the classical theories of material behavior, such as linear and non-linear viscosity, linear and non-linear visco-elasticity, and linear and non-linear elasticity. This model is not only highly useful for conceptual purposes (e.g. Reference Truesdell, Noll and FlüggeTruesdell and Noll, 1965) but also may be used to solve certain boundary-value problems exactly (Reference Coleman, Coleman, Markovitz and NollColeman and others, 1966). However, usually the model is too general to allow analytical solutions. A useful expediency is, then, to choose special simple materials, explicit enough to allow solution of boundary-value problems, but not so special as to exclude the physical phenomena one wishes to model. For such choices, there are two obvious options (1) arbitrary choice of constitutive equations for special simple materials, and (2) development of a hierarchy of more and more complex materials by rational approximation schemes. option 1 is not discussed here. Within the realm of option 2, the choice that gives answers that are most acceptable to those who do experimental work on fluid-like materials is the method of retarded motions (Reference Coleman and NollColeman and Noll, 1960). There, the functional is approximated by a sequence of functionals of motions more and more retarded, i.e. stretched out, in time. For incompressible materials, the first-order approximation is the linearly viscous fluid. The second-order approximation isFootnote ^†

(1)

where T_ij is the stress, P is an isotropic pressure, μ₁, μ₂, and μ₃ are material constants,

(2)

is twice the deformation rate (“strain-rate”) D_ij used in fluid mechanics, v_i is the velocity, and

(3)

The quantities Aij ⁽¹⁾ and Aij ⁽²⁾ were introduced by Reference Rivlin and EricksenRivlin and Ericksen (1955), and are now commonly called Rivlin–Ericksen tensors. The material model (1), properly referred to as the “incompressible fluid of grade 2”, is commonly called the second-order fluid. The special case μ₂ = μ ₃ = 0 recovers the linearly viscous fluid, while setting μ₂ = 0 recovers the Reiner–Rivlin fluid (Reference RivlinRivlin, 1948). As with any good putative material model, the basis of the second-order fluid is a strong theoretical argument, but its ultimate acceptance as a constitutive equation for any real material depends upon comparisons between exact solutions to boundary-value problems and experimental data. Steps in that process for ice are presented here.

One distinguishing feature of the second-order fluid model is that it allows for non-zero and unequal normal stress differences in shearing flows. The normal stress differences are defined by

(4a)

(4b)

For a shear rate au ∂u ₁/∂x ₂ = K, N ₁ = −2μ₂k² and N ₂ = (2μ₂ + μ₃)k². It is easily seen that a power-law constitutive equation of the form

(5)

where cannot yield non-zero values for N ₁ and N ₂ in rectilinear shearing flows. Reference GlenGlen’s law (1955), commonly used to represent the behavior of ice, is of the form of Equation (5), with n generally found to be about 3. Normal stress effects, if present, thus cannot be predicted by Glen’s law.

Direct measurements on a great variety of non-linear fluids, most of them polymer solutions or melts, have shown that most exhibit non-zero and unequal normal stress differences in shearing (e.g. Reference Bird, Bird, Armstrong and HassagerBird and others, 1977). A power law can accurately represent the secondary creep behavior of these materials, but cannot predict any of the non-linear phenomena associated with normal stresses. Ice is not a polymer in the classical sense, and for that reason it might be inferred that this experimental evidence carries little weight. However, it is significant that few materials have ever been found that exhibit non-linear shear viscosity that do not also exhibit normal stress effects. It is interesting to note that non-linear kinetic theories for fluids (e.g. Reference Bird, Bird, Armstrong and HassagerBird and others, 1977, part 2; Reference Truesdell and MuncasterTruesdell and Muncaster, 1980) arrive at a similar conclusion based on detailed considerations of momentum transfer on a molecular scale.

The second-order fluid model successfully represents data from viscometric measurements for numerous materials. In addition, its application in various boundary-value problems leads to predictions of normal stress effects that are qualitatively consistent with laboratory observations. Included among these are the phenomena of rod climbing (the Weissenberg effect) (Reference Joseph and FosdickJoseph and Fosdick, 1973), die swell (Reference Bagley, Schreiber and EirichBagley and Schreiber, 1969), secondary flows in non-circular conduits (Reference Langlois and RivlinLanglois and Rivlin, 1963), cross-stream migration of particles (Reference Ho and LealHo and Leal, 1976), and free-surface bulging in open-channel flows (Reference Wineman and PipkinWineman and Pipkin, 1966). All of these effects are due to non-zero normal stress differences. The second-order fluid is the simplest model that can represent this phenomenology.

Determination of Material Parameters

Although the second-order fluid has been suggested previously (Hobbs, 1973; Reference Rundle and PassmanRundle and Passman, 1982) as a model for geological materials in creep, there has been no simple and direct method available for determining the three coefficients in Equation (1) from standard laboratory tests. However, Reference PassmanPassman (1982) has recently presented an exact solution to the boundary-value problem corresponding to a cylinder of second-order fluid in biaxial creep. By fitting the analytical expression for axial strain versus time to laboratory data, the material properties can be determined directly. Indeed, this has been done successfully for polycrystalline halite (Reference PassmanPassman, 1982) and polycrystalline olivine at high temperature (a paper in preparation by D.F. McTigue and S.L. Passman).

Consider a right circular cylinder of ice, as commonly tested in a biaxial device (Fig. 1), and choose a rectangular Cartesian coordinate system with x ₁ along the axis of the cylinder. Let the cylinder have a length L in a reference configuration. Assume the ice is incompressibleFootnote ^* and that it undergoes a homogeneous motion, with its length as a function of time being l(t). Let

(6)

Clearly, α is related to the infinitesimal (“engineering”) axial strain ∊ by α = 1 + ∊ and to the natural strain e by e = Inα. Then the deformation gradient tensor (see e.g. Reference SchowalterSchowalter, 1978) is diagonal, and its components F_iI have the form

(7)

the second two terms being equal because of the choice of coordinate system, and having the value 1/α because incompressibility requires that the determinant of F_iI be 1. The Rivlin–Ericksen tensors (1) (2) A_ij and A_ij (Equations (2) and (3)) are easily computed from F_iI . They take the forms

(8)

(9)

where a = α/α, and a superposed dot indicates differentiation with respect to time.

Fig. 1. Definition sketch for the biaxial creep test.

Now by Equation (1), if⁽¹⁾ A_ij and⁽²⁾ A_ij have the forms of Equations (8) and (9), then T_ij must have the form

(10)

where v is the confining pressure, and σ is the axial stress in excess of the confining pressure. It then follows from Equation (1) that

(11)

and

(12)

If one eliminates v from these equations, one obtains

(13)

(14)

with the initial conditions

(15)

(16)

Reference PassmanPassman (1982) details a number of solutions to this system, corresponding to different types of fluids in polymer rheology. The only situation of interest here is

(17)

The solution is facilitated by the introduction of the following parameters:

(18)

(19)

(20)

The parameters λ and μ may be either positive or negative, but may not vanish or be infinite. Likewise, ξ² may vanish or be positive or negative. Discussion is limited here to the case ξ² > 0, which yields physically reasonable results and has been found to provide the best fits to laboratory creep data for geological materials (Reference PassmanPassman, 1982). For this case, the solution is

(21)

(22)

In the following section. data for ice are fitted by Equation (22) in order determine the material constants μ₁, μ₂, and μ₃.

Experimental Procedure and Data

The creep curves given here were obtained from ice grown in a cold room at about −8° C. A mold, similar to that described by Reference ColeCole (1979), was about 80% filled with snow or crushed ice which had been passed through a #30 sieve but had been retained on a #40 sieve. It was evacuated with a vacuum pump for about 15 min and then flooded with de-aerated, de-ionized, and doubly distilled water at 0° C. While still under vacuum, the mold was tapped lightly to release any trapped air bubbles. It was allowed to freeze for about 20 min, and then the vacuum was released and the top cap of the mold unscrewed to allow for expansion of the freezing sample. The whole was then left in the cold room, usually overnight, to freeze solid.

The samples were removed by gently warming the mold and allowing the ice to slide out. At this point the samples were approximately 140 mm long and 85 mm in diameter. There was usually a thin core of air bubbles running down the center of the sample. However, because the sample was cut longitudinally into four pieces and then each piece turned on a lathe, the bubbly center part was removed. The final samples tested were 12 mm in diameter and 40 mm long, with stainless steel caps frozen on to the ends. Thin sections of each sample were examined to check that the crystals were randomly oriented. The average grain-size was 1 mm.

The samples were deformed under a constant uniaxial dead-weight compressive load in a stainless steel jig. Friction was minimized by using linear ball bearings, and the ensuing deformation was measured with an LVDT. The whole jig was placed in a high-pressure cell with two electrical feedthroughs, one for the LVDT connections and one for a thermocouple. The pressurizing fluid was Dow Corning silicone fluid 200, which served also to prevent evaporation of the sample. The whole was placed in a cold room, and the compressive creep deformation was recorded as a function of time on a strip-chart recorder. The data were converted to give a normal creep curve of strain against time for a constant hydrostatic pressure. Some results of these tests have already been discussed (Reference Jones and ChewJones and Chew, 1983).

Four tests were done at the same temperature (−9.5° C to −9.8° C) and axial stress in excess of the confining pressure (−0.47 MPa), but at different confining pressures (0, 37, and 50 MPa). The tests were carried out for 23 to 35 d. The data reveal the onset of tertiary creep at late time; i.e. there is evidence of strain-softening behavior. It is known for other materials that the second-order fluid can represent primary and secondary, but not tertiary, creep (Reference PassmanPassman, 1982)Footnote ^* . Thus, for the purpose of fitting Equation (22) to the measured strains, the data have been truncated at the onset of tertiary creep.

The three coefficients characterizing the second-order fluid were determined by multiple regression, fitting Equation (22) to the data as shown in Figure 2. The results are given in Table I. The properties inferred from tests 2, 3, and 4 are quite consistent. Because these tests were run at essentially the same temperature, possible temperature dependence of the coefficients is not revealed. A wide range of confining pressures was used, and there is no discernible effect on the material properties. This observation lends credence to the assumption of incompressibility. The fits shown in Figure 2 and the corresponding coefficients recorded in Table I were determined using all of the available data. The values obtained for μ₂ and μ₃ are highly sensitive to the data at early time. This was substantiated by repeating the regression analyses with the first point in each data set deleted. The results of this procedure, which are not shown here, were values for all three coefficients in all four tests within a few percent of their respective means. Thus, the scatter in the values recorded in Table I can be ascribed primarily to uncertainty in the creep response at early time. Test 1 in particular is strongly affected by the first data point, and the resulting parameters are not included in the calculated averages.

Fig. 2. Fits of Equation (22) to data for primary and secondary creep. Parameters determined by regression are shown in Table I.

Table I. Constitutive Coefficients, Temperature, and Confining Pressure Forfour Biaxial Creep Tests

It should be emphasized that all of these tests were carried out at the same deviatoric stress; only the confining pressure was varied. In general, of course, the coefficients are expected to be different for different shear rates (e.g. μ₁ may vary in the fashion of the term multiplying A_ij ⁽¹⁾ in Equation (5)). this type of generalization has been discusesd by Reference MorlandMorland (1979,Reference Morland and Spring1981). The results presented here, then, are not intended to suggest that the second-order fluid model gives a complete representation for the behavior of ice. Rather, attention is confined to a limited test condition in order to obtain a first estimate of possible normal stress effects.

For shearing flows, the average values obtained for ice yield normal stress differences of

(23a)

(23b)

where the shear rate k is in units of s⁻¹. It is interesting to compare these results with those for other materials modeled as second-order fluids. Polymer rheologists generally find from direct viscometric measurements that N ₁ > 0 and N ₂ < 0 (e.g. Reference Bird, Bird, Armstrong and HassagerBird and others, 1977), and observe that | N ₂ | is about one order of magnitude smaller than | N ₁ |. The results for ice which are inferred here from non-viscometric measurements, contrast strongly with data for polymers; N ₂ is positive and two orders of magnitude greater than N ₁. These differences are not surprising because the mechanisms that give rise to normal stresses in polymers are interactions of long-chain molecules. That is clearly not the case for ice. For salt (Reference PassmanPassman, 1982) and olivine (a paper in preparation by D.F. McTigue and S.L. Passman), relationships similar to those observed in ice are found: N ₁ > 0, N ₂ > 0, and N ₂ » N ₁ Footnote ^*

Free-Surface Depression in Open-Channel Flow

A significant qualitative difference between the flow of a second-order fluid and that of a power-law fluid is the shape of the free surface in steady flow in an open channel. Reference Wineman and PipkinWineman and Pipkin (1966) showed that a negative second normal stress difference N ₂ will result in a free surface that is convex upward in this type of flow. For slow flow, the kinematics are dominated by the linear term in Equation (1). An approximation to the non-linear normal stresses can then be determined based on the velocity field for a linearly viscous fluid. This results in a non-equilibrated normal stress on the nominal (horizontal) free surface, which is balanced by a proportional rise or fall of the surface. Reference TannerTanner (1970) implemented this observation as an experimental method for determining the normal stress difference N ₂. Reference SchowalterSchowalter (1978) has presented a general expression relating the maximum surface deflection to N ₂ for flow in a semicircular channel:

(24)

where ρ is the density, g the acceleration due to gravity, h the central (maximum) rise or fall, β the channel slope, and R the channel radius(Fig.3). The velocity field, to this order of approximation, has the same parabolic profile as is found for a Newtonian fluid. Thus, by Equations (1) and (4b), the normal stress difference for the second-order fluid is found to be

(25)

where r is the radial coordinate. Subtituting where r is the radial coordinate. Substituting Equation (25) into Equation (24) and integrating gives

(26)

Because 2μ₂ + μ₃ > 0, h is negative, representing a concave surface, i.e. a depression.

Fig. 3. Definition sketch for channel flow.

It is recognized, of course, that most glacial surfaces exhibit a concave surface in the accumulation region and a convex bulge in the ablation region. These are generally well understood in terms of the equation of conservation of mass in an axially extending or shortening flow region (Reference RaymondRaymond, 1971), where vertical flow results in surface displacements over the thicker central part of the channel. Nonetheless, it is interesting to consider whether the normal stress effect examined here could contribute measurably to observed surface deformation. Using the average values for μ₁, μ₂, μ₃ and μ₃ given in Table I and taking ρ = 9 × 10² kg/m³, g = 9.8 m/s²,h can be estimated for various channel slopes and radii. Results for such a calculation are shown in Figure 4. Because the normal stress effects are of second-order in the rate of shear deformation, they become significant only at higher rates, i.e. for large channels and/or steep slopes.

Fig. 4. Calculated surface depression as a function of channel slope and radius for a semicircular channel, using mean values from Table I.

The principal stresses at the free surface are affected by the non-linear normal stress terms, and this may influence the development of crevasses. In the case of a power-law fluid, the normal stress differences in a steady uniform flow are zero, and the principal stresses are a tension and a compression of equal magnitude, oriented at 45° to the glacier margin (Reference NyeNye, 1952). Crevasses normal to the principal tension are favored. For the second-order fluid, the normal stress difference causes a change in magnitude and orientation of the principal stresses. In a semicircular channel, the principal stresses T ₁ and T ₂ are given by

(27)

(28)

where B = ρg sin β/2 μ1. As a sample calculation, consider a flow for which p = 9 × 10² kg/m³, β = 2°, and μ₁, μ₂, and μ₃ are given the average values shown in Table I, At the margin of a 500 m radius flow, the compressive principal stress T ₁ is decreased by 12% and the tensile principal stress T ₂ is increased by 14% relative to the purely viscous case (μ₂ = μ₃ = 0). Furthermore, the angle θ at which T ₁ meets the boundary is given by

(29)

which, for large values of the right-hand side, leads to the approximation

(30)

For the example cited in the foregoing, ө = 41°; i.e. the normal stresses cause a small rotation of the direction most conducive to crevasse formation.

Conclusions and Discussion

Glacier ice is most often represented by a power-law constitutive model that provides good estimates of shear stresses and overall kinematics in simple flow configurations. However, such a model cannot represent the non-zero normal stress differences that have been observed in virtually all non-linear fluids for which the appropriate measurements have been made. The simplest constitutive model that arises from strong theoretical arguments and embodies these effects is the second-order fluid. The model provides an excellent fit to data for primary and secondary creep of ice in biaxial tests in a single deviatoric stress state, and the inferred material constants are reasonably consistent among the four tests considered. The results indicate positive first and second normal stress differences, which contrasts with other materials. This provides a point of departure for inquiry into possible normal stress effects in natural glacier flows. The free-surface deformation in a steady channel flow has been considered briefly as an example of such phenomena, and it appears that it could be discernible for thick glaciers on steep slopes. In addition, the normal stresses alter the principal stresses in a fashion that favors crevasse formation in uniform flow.

A limitation of the present analysis that should be noted is that the data-fitting procedure is somewhat insensitive to changes in the constants μ ₂ and μ₃; i.e. good fits of Equation (22) to the creep data can be found for broad ranges of these parameters. The coefficients μ₂ and μ₃, are determined principally by the curvature of the strain-time plot at early time. Most standard biaxial tests, including those reported here, provide only sparse data in this critical region, and regression analyses are consequently ill-constrained. This could be rectified, of course, in tests run explicitly for the purpose of examining the early time response in detail. A second reason for the low parameter sensitivity is that the response of a material in biaxial creep is dominated by its shearing propertiesFootnote ^* . While the biaxial creep test is not ideally suited to the measurement of normal stress coefficients, it has the advantage that it is a standard test for which experimental equipment is easily available and highly refined, and for which there is a large body of existing data. Despite these limitations, it is emphasized that the four data sets considered here all yield consistent values for μ₂ and μ3, implying that they are reasonably well constrained. Of course, the four tests cited were performed in the same deviatoric stress state, so that any non-linear rate dependence of the shear stress is not revealed.

The Reiner-Rivlin fluid, which can be obtained from Equation (1) by setting μ₂ = 0, also embodies certain normal stress effects. However, all fluids for which normal stresses have been measured exhibit a non-zero value of N₁ (cf. Equation (4a)), contradicting the Reiner–Rivlin model (see e.g. Reference Bird, Bird, Armstrong and HassagerBird and others, 1977). Furthermore, it can be shown that the behavior of the Reiner–Rivlin fluid in a biaxial creep test is indistinguishable from that of the linearly viscous fluid. Thus, it, too, is incapable of representing both primary and secondary creep. The second-order fluid, in contrast, does capture the transient primary creep phenomenon. However, it should be noted that the observation of primary creep in a biaxial test does not necessarily imply the existence of normal stresses. While they are closely linked in the present analysis, it is important to recognize that the representation of primary creep is not unique to models of this type. For example, the four-element visco-elastic material known as the Burgers substance, like the second-order fluid, can capture an initial elastic displacement, primary creep, and asymptotic steady secondary creep, but does not predict non-zero normal stress differences.

Another limitation, which is not unique to the second-order fluid model, lies in the extension of laboratory measurements on a small homogeneous sample to large-scale field situations in which ice may be highly inhomogeneous and extensively fractured. It is not known whether such circumstances may enhance or suppress normal stress effects in the bulk material when it is viewed as a continuum.

It remains to be seen whether normal stress effects of any sort are discernible in natural glacier flows. Experience with other non-linear fluids in other contexts suggests that such effects can be of practical significance. Calculations based on power-law models can never point toward these phenomena. The search for normal stress effects in ice must begin with a constitutive model capable of representing them in a self-consistent fashion, and the second-order fluid is the simplest such model. An important step in the assessment of this sort of phenomenology is the quantification of material properties. This has now been initiated for ice, and the results serve as a point of departure for further investigation.