Proposal for a constitutive equation of temperate firn

CoM Regions Science and Technology, 13 (1986) 1-9 Elsevier Science Publishers B.V.~ Amsterdam - Printed in The Netherlands

PROPOSAL FOR A CONSTITUTIVE

EQUATION

OF TEMPERATE

1

FIRN

W. Ambach Institut fb'r Medizinische Physik Mb'llerstrasse 44, A-6020 Innsbruck (Austria)

and H. Eisner Institut fffr Experimentalphysik, Scho'pfstrasse 4 I, A-6020 Innsbruck (Austria)

(Received May 8, 1985; accepted in revisedform January6, 1986)

ABSTRACT

2. BASIC MEASUREMENT DATA

From deformation measurements at a 20 m deep firn pit, a constitutive equation for temperate firn under natural stress has been derived. The analytical and numerical formulation passes over into Glen's flow law for ice at the transition o f firn to ice. Basis o f the analysis is the determination o f the shear strain rate from the tilt rate o f the pit axis, under consideration o f the longitudinal and vertical strain rates, the compressibility o f the tim, as well as the triaxial pit deformation.

The coordinate system applied is represented in Fig. 1, the measuring quantities being as follows: strain ea of pit's cross-section in flow direction x, measured in various depths; strain eb of pit's cross-section in transverse direction z 1 measured in various depths; strain eh of the compression of layers in various depths; tilt 7x, of the pit axis to lead line (compare Fig. I); density p, dependent from the depth; surface slope a along the flow line; depth y of the markers, in relation to the surface. At 14 different depths along the wall of the pit, 6 to 7 markers each were placed, the relative distance of which was measured once a year. The pit crosssections were approximated by ellipses, with the major axis pointing in flow direction, and the minor axis pointing transverse to the flow direction. The centrepoints of the ellipses at different depths form

1. INTRODUCTION

On a 20 m deep firn pit, measurements of the cross-section deformation, of the compression of layers and of the tilt rate of the pit axis have been carried out over a period of 11 years. In the present paper, the analysis is realized with the intention of introducing a constitutive equation for temperate tim, approximating asymptotically into Glen's flow law at the transition between firn and ice. The proposed constitutive equation of temperate firn consists of 2 parts, one related to deviatoric terms, the other one related to isotropic terms due to the compressibility of the material. Essentially, the evaluation was achieved according to a procedure described by Paterson (1981, p. 108), taking into consideration the compressibility of the firn.

0165-232X/86/$03.50

© 1986 Elsevier Science Publishers B.V.

~ y

x~

PROFI,E

Y

.

yr

Fig. 1. Coordinate systems applied. Because of small angles, it holds true that arctan ~x ~ ~/x and aretan 3,x',~ 7x:

movement of the tim, the bottom of the pit is lowered from a depth of 20 m to a depth of 34 m, whereas the distance between the marker closest to the top and the marker closest to the b o t t o m decreases from

the pit axis. The strains of the ellipse axes and of the compression of layers are represented as a function of time in Fig. 2, and the tilt of the pit axis is given as a function of time in Fig. 3. Due to the immersing

~a 13

0.4

®

4/1

,.~~ 11/14

0.3

.~-~,1~''"

617/12 2/5

0.2

0.1 ---it-

67

I

I

I

I

I

I

I

I

I

I

I

68

69

70

71

72

73

74

75

76

77

78

Fig. 2a. Strains of the major ellipse axis of the pit cross-sections for individual levels over the years 1967-78 (No. 1 bottommost level).

h

Laye~

12 10 11/13

?'/

- 0.4 (~

~'::"

0.3

.

0.2

..;;~f'"

O.1

?;;4

67

I"~.

9/8

.a....~ elf

I..~.///

2

~...~-~

, ~...'~'-'~

I

I

I

I

I

I

i

!

I

i

i

68

69

70

71

72

73

74

75

76

77

78

Fig. 2b. Strains of the compression of individual layers over the years 1967-78. Layer 1 is situated between levels 1 and 2 (No. 1 bottom-most layer).

Level

2/4

No. ~

~b 7

-0.5

9 8/3/5/12 6 10

® -0.4

/ \d'

14/1/11

-0.3

13

/

-0.2 It'.'" s."

-0.1

i

I

I

I

I

I

I

I

I

I

I

I

67

68

69

70

71

72

73

74

75

76

77

78

Fig. 2c. Strains of the minor ellipse axis of the pit cross-sections for individual levels over the years 1967-78 (No. 1 bottommost level).

"~x = _15.5x10-3

0.05

I

67

I

68

69

~' ~

70

71

I

I

72~...

a-1

II

I

1

I

I

74

75

76

77

78

M

measured - 0.10

Formation of c r e v a s s e

.j.~%

~

"

- 0.20

- 0.30

Fig. 3. Tilt 3'x of the pit axis as a function of time: data of measurements (crosses) and calculated curve (full line). 17.6 m to 11.9 m. The range of density of the entire pit amounted to between 631 kg/m a at the topmost level in 1967 and 864 kg/m a at the b o t t o m m o s t level in 1978. The slope of the surface was measured along

the movement distance; an averaged value of 5.14" was applied in the analysis. Temperature measurements of the firn down to a depth of 15 m show that during winter negative

temperatures occur down to a depth of 7 m from the surface, whereas between June and November, the temperature amounts to 0°C throughout. The negative temperatures within the topmost layers do not cause any significant disturbances with respect to the strain rates.

X(p*) = bxy/rxy 3.2.2

Calculation of deviatoric stresses by means of eqn. (4) with the deviatoric strain rates (compare paragraph 3.3.3) yields • 1

, 3. D E V I A T O R I C TERMS OF THE C O N S T I T U T I V E EQUATION

(6)

•

Oxx

=

exx ~k(p*)

,t

,

, eyy ° Y Y = x (p*)

'

, azz

ot

-

ezz k(p*)

(7)

3.2.3

Calculation of D(p*) by eqn. (5) gives

3.1. Conception

The following deviatoric constitutive equation proposed by Mellor (1975, p. 275) is introduced: eeff

=

O(p*)reaff with p* = P / (Pi - P)

(1)

and 2r2eff = ox2x+ Oy~ + Oz2z+ 2r2xy,

(2)

2~2ff

(3)

.,2 + t:lTy 0,2 + ~:zz "2 + = ~2XX

2~x2y

Because of the orientation of the coordinate system in the direction of the flow, it holds true that rxz = ryz = 0 and exz = ~:yz = 0. Contrary to Glen's flow law, at constant temperatures, D is a function of p*. The setting of 4eft proportional to reaffis supported by the indication, that in a graph with logarithmic scales on both coordinate axes, a system of parallel straight lines results for various snow densities (Mellor, 1975, Fig. 6, p. 274). The same slope of the straight lines is an indication that for the term reff the same exponent can be set for ice and firn. From the parallel displacement of the straight line it can be deduced further that the factor of proportionality depends on the density. When applied to deviatoric terms, it holds true that • t

el/

,

t

= X(p )off

with

X(0 *) = D(p* )r2ff.

(4) (5)

3.2. Analysis of the deviatoric terms 3.2.1

Calculation of k(p*) by means of eqn. (4) with shear strain rate and shear stress (compare paragraphs 3.3.1 and 3.3.2) yields

D(p*)

=

X(p*)/re2ff

.

Furthermore, the function D(p*) is approximated by O(p*) = A + O l e x p ( - d , p * ) + O2exp(-d2p*).

(8)

with A, D1, D2, d], d2 being constants of the constitutive equation. 3.3. Evaluation of the deviatoric terms 3.3.1 Shear strain rate

By definition, it holds true that •

1

exy=~"

(au -~y + ~--xV) .

(9)

According to measurements, made so far, it is true that ~v/ax .~ 0.04 ~u/ay, so that av/ax can be neglected, compared with ~u/~y (Eisner et al, 1984a). Calculation of au/~y is achieved according to Paterson (1981, p. 108) and Raymond (1971) by au _ a3'x + v

0y

at

a~/x

ay

- ~/x

(au

ax

a~)

- 7z

au

g -z

(10)

In this equation, 7x, ")'z are the tilts of the pit axis in the xy-plane and yz-plane respectively, and u, o are the velocity components. The profiles of the pit's tilt have been approximated by means of the best fitting straight lines, therefore it is true that a'yx/ay = 0. Moreover, it is true that 7z = 0. The slopes of the best fitting straight lines 7x are represented as a function of time in Fig. 3. According to eqn. (10) for the case ~'x = 0, ~'z = 0 and aTx/ay = 0, one obtains

au

aTx

- --

ay

(11)

at

")'x = 0 is complied with, if the pit axis is perpendicular to the surface (3'x' = tan a). In this particular case, the strain rates au/ax and av/ay do not have any influence on the tilt rate of the pit axis 7x. The tilt rate of the pit axis is then brought about by the term au/ay only. The temporary course of "Ix(t) was disturbed from 1975 onwards by the formation of an up to 4 m wide transverse crevasse, 33.5 m to the upper side of the pit (Fig. 3). For the further evaluation therefore, the values of the pit's tilt up to 1974 are of interest only. For the. calculation of 7x, the function 7x(t) for the period from 1967 to 1974 was approximated by means of the best fitting polynomial (Fig. 3) 7x = a + bt + ct 2 .

(12)

The position of the pit axis, vertical to the surface, is obtained by to--- 4.74 a, corresponding to a period of time between 1971 and 1972 (Fig. 3). Before the time when the pit axis reaches a perpendicular position, this fitting curve does show a flatter slope than afterwards. This is in agreement with the influence of the longitudinal and the vertical strain rates. The tilt rate 7x of profiles with an upward slope is less negative because of the longitudinal strain rate, whereas the tilt rate of profiles with a downward slope is a more negative one. For 7x = 0, the value of 3'x was calculated from eqn. (12) by differentiation. As a further step, it results from eqn. (9) with au/ay >>av/ax, that 1

(13)

~xy = -i ~/x .

The value bxy is assumed to be constant with time for the period from 1967 till 1974. The same assumption is applied in Fig. 3, which shows a good agreement between the measured and calculated data. 3.3.2 Shear stress

According to Paterson (1981, p. 109) it is true that 7"xy =

-Ff

Y

pg sin ot d y , o

(14)

with F allowing for the influence of the transverse profile of the glacier bed. Subsequently, it is true that F = 1/2 (Nye, 1965). This value holds true for all shapes of cross-sections in layers close to the surface and along the centre line of the glacier, with the exception of an infinitely extended plate. Further correction terms for Zxy have been neglected (terms T, G according to Paterson (1981), p. 100). 3.3.3 Deviatoric

s t r a i n rates

The deviatoric strain rates can be calculated from the deformation of the pit by the following steps: At first, the strains for each level are determined according to (Fig. 2) Aa i ea=~ " , ai

Ah i eh=Z

Ab i ,

eb=Z

hi

(15)

bi

The values Aai/a i are added up step by step, from year to year, always using the value of the respective previous year for ai; Similarly for eh, eb. Contrary to the temporary course of 7x (Fig. 3), no disturbing influence due to formation of crevasses is to be recognized with Ca, eh, eb. The time functions ea, eh, eb are fitted by means of best fitting straight lines; consequently the strain rates ea, eh, eb are temporarily constant. Because of the tri-axial state of stress, the following equations hold true for the strain rates (Bader et al., 1951): =

+

Ch = - - P~x x + eYy'~y -- P~zzzz , :

-

+

(16)

+

with ~xx being the uni-axial strain rate in the x-direction (subscript), exclusively originating from Oxx (superscript); similarly for b:~, ~z. v is the viscous Poisson's ratio, to be calculated as a function of density by (Bader et al., 1951) ['1.48 (1 - 0.80 × 10 -3 p)(1 + 3.87 × 10 -3 p)] -1 v

L

p • 10-3

J (17)

After determination of &x.~, by,y, ~zz from eqn. (16), the strain rates exx' eyy' ezz for the further evaluation result from

=

~yy

+

=

_

=

-

-

-

l~xxx + ~

(18)

p e"zz zz ,

_

-

3.3.5 Approximation From eqn. (5), D(p*) is calculated and approximated by eqn. (8), namely

gz,

+

,

D ~ * ) = A +D l e x p ( - d l p * ) + D~exp(-d2p*)

and the corresponding deviatoric strain rates from

By the best fitting curve one obtains (Fig. 4): A = 6.04 X 10 -15 kPa-3s -1, D1 = 3.94 X 10-1°kPa-3s -1, D2 = 7.07 × 10-1akPa-as -1, d l = 2.071, and d2 = 0.419, where dl and d2 are independent of the system of units. As an attempt, the evaluation was repeated for the limiting values v = 0 and v = 0.5 in order to judge the influence of v on the final result. Generally, the value of v depends both from the density and from the state of stress (Salm, 1977). However, from the re-evaluation it was found that the value v only slightly influences the result (Fig. 5a). For ice at 0°C it holds that A = 5.3 X 10-1SkPa-3s -1 (Paterson, 1981, p. 39, Table 3.3). In order to obtain the numerical value of A with a higher degree of accuracy, the approximation was additionally made for p* > 6, by using the best fitting curve through the points of measurements (Fig. 5b). The deviatoric constitutive equation for temperate firn therefore reads numerically for natural stress conditions with eeff (s -1) and reff(kPa) for p* > 2.5

1

1

~)y =~yy - ~- J, ,

(19)

1 :

-

7

'

with J l = exx + eyy + ezz •

(20)

Equations (16) and (18) differ in signs. Thereby allowance is made for the fact that the deformation measurements have been carried out on a cylindrical cavity and not on a cylindrical solid body (Eisner et al. 1984b). The dependence of the viscous Poisson's ratio from the state of stress remains unconsidered (Salm, 1977). 3.3.4 Deviatoric stresses For the calculation of Oxx, ' a'yy, Ozz, eqn. 7 is applied and reff results from eqn. (2). 5 0 -

m I 4 0 .bg t~ I O A 3 0 -

C¢" v

2 0 -

.tO-

0

0

i

i

i

]

2

4

6

0

Fig. 4. Curve of D(o*) for t~* > 2.5 by approximation from eqn. (8).

.......

1

10

.........

i,

12

i ::I. 4

i ~.B

50.)4 u) equ.

#.

40

30

®

17

ED

"Y = 0 . 0

+

'V=O. 5

O,-

c,

20-

10-

0 2

0

4

6

10

8

12

14

1~

7

®

I 0.

Q ,

5

A(Ice,OOC)

,~*

o

0

2'

.4 '

'

6'

8'

1' 0

I' 2

1 4'

"t6

Fig. 5. Curve of D(p*) by appro×imation from eqn. (8). a. Evaluation with u = O, v = 0.5 and v depending on the density [eqn. (17)]. b. Evaluation/or p* > 6 in order to determine the asymptotic value as p* -* oo.

~eff = [6.04 X 10-1s+ 3.94 × 10-1°exp(-2.071p *) + 7.07 × 10-13exp(_0.419p,)] reif3 .

(21)

For p* ~ *% the deviatoric constitutive equation for firn (eqn. (21)) passes over into Glen's flow law for ice in its analytical and numerical formulation. This is considered to be the confirmation for the expediency of eqn. (1).

4. ISOTROPIC TERMS OF THE CONSTITUTIVE EQUATION 4.1. Conception The following i~tropic constitutive equation was introduced by applying the f'trst invariants of the strain rate and the stress tensor

with p* = p/(,o i - P),

J; = H(p*)Ia ,

~ o. is fulfilled. The numerical formulation of the isotropic term of the constitutive equation therefore reads with Ja (s -a) and 11 (kPa) for p* > 2.5

(22)

p*

and with the condition H(,o*) ~ 0 at the transition of compressible firn to incompressible ice.

Jl = [4.74 X 10-nexp(-1.081p *) 4.2. Evaluation of the isotropic terms of the constitutive equation

+ 9.64 × 10-Bexp(-0.131p*)] 11 •

Figure 6 shows the ratio J1/11 as a function of p*, where Jx and 11 are calculated by J; =

+

+

zz,

5. F I N A L CONCLUSIONS

(23)

t

11 -- 3(Oxx + Oyy + Ozz) ,

The constitutive equation of temperate firn under natural state of stress is composed of deviatoric and isotropic terms as follows:

(24)

where b~ are known from eqn. (18), Oxx, O'zz from eqn. (7), and where Oyy is calculated from the density profile by

1

+ 7 J' where with

y

Oyy

=

-f

p g cos a d y .

(25)

0 The data are plotted in Fig. 6, where the bars are the standard deviations, obtained from individual layers. The plot H(p*) was approximated by H(p*) = Haexp(-h aP*) + H2exp(-h2p*)

8ti

(28)

means the Kronecker symbol. According

b'q = D(p*)r2effa b ,

[see eqns. (4, 5)] and

J1

[see eqn. (22)]

= H(,o*)/l

eqn. (28) reads finally 1 ei/ : D(p * )'reffOtj ' + 7 H(p*)I18i/

(26)

with Ha = 4.74 X 10-1akPa-ls-a,H2 =9.64 X 10 -13 kPa-ls -1, hi = 1.081, h2 = 0.131, where the values ha, h2 are independent from the system of units. By this approximation, the condition H(p*) ~ 0 for 6

(27)

or

(29)

(30)

~i/= [A + D l e x p ( - d , p * ) + D2exp(~/2p*)] r2effob 1

+-- [Hlexp(-hlp*) + H2exp(-h2p*)] I18t/. 3

-

(31) T ta

I

4

~

2

o o

2;

4;

6;

8;

.'t ; 0

,l ; 2

"t '.,a

I 6

Fig. 6 Curve of H(p*) for p* > 2.5 by approximation from eqn. (26). Bars indicate the standard deviation, obtained from individual layers.

with A = 6.04 × 10-1SkPa-as -1, D1 = 3.94 X 10 -1° kPa-as -1, dl = 2.071, D2 = 7.07 X 10-1akPa-as -1, d2 = 0.419, H1 = 4.74 × 10-11kPa-I s-I, hi = 1.081, /-/2 = 9.64 × 10-13kPa-I s -1, h2 = 0.131. This numerical formulation is valid for p* > 2.5. At the transition of firn to ice it holds p* ~ 0% thus Glen's Law is fulfilled by •

2

t

ei/= AreffOi/

and

"

_

3

eeff - A'reff .

hi (Ahi) J1 t u,u x,y,z x', y', z' ol

(32)

ACKNOWLEDGEMENTS

The authors would like to express their sincere gratitude to the Austrian Academy of Sciences (Akademie der Wissenschaften), Vienna, for the financial support of the field works, to the Federal Ministry for the Interior (Bundesministerium fiJr Inneres), Vienna, for the transportation of the material; many friends and colleagues for their assistance with the field works, and Mr. K. Leitl for his co-operation in evaluating the data.

thickness of a firn layer (differences). first invariant of strain rate tensor. time. velocity components. coordinates. coordinates. surface slope. tilt of pit axis in system x, y and x', , dy dy' y respectively (Tx = ~--, 7x' = ~ ) .

~z ea eb eh

xf:*) P

P, Pi, P*

oi/(ai/) I

LIST OF SYMBOLS

A

a,b,c

ai, ( /'i)

b~, (AbD D1, D2, dl, d2 D(p*)

F g

H~*) H1,/'/2 hi, h2

constant of constitutive equation. coefficients of an approximating polynomial. major axis of ellipse of pit cross-section (differences). minor axis of ellipse of pit cross-section (differences). constants of constitutive equation. parameter of constitutive equation, depending on p*. components of strain rate tensor (deviator). uni-axial strain rate in x-direction (subscript), caused by Oxx (superscript), analogously in y, z directions. proportionality factor. acceleration of gravity. parameter of constitutive equation, depending on p*. constants of constitutive equation.

1"eft

transverse tilt of pit axis. strain of major ellipse axis of pit crosssection. strain of minor ellipse axis of pit crosssection. strain of compression of layers. effective strain rate. parameter of constitutive equation, depending on 0". viscous Poisson's ratio. # density (Pi = density of ice; p* = ) Pt-P

components of stress tensor (deviator). effective stress. Kronecker's symbol.

REFERENCES Eisner, H., Ambach, W. and Schneider, H. (1984a). Time dependent tilt of a 20 m deep f'lrn pit. Polafforschung,

54(2): 85-93. Eisner, H., Ambach, W. and Schneider, H. (1984b). Evaluation of strainrate measurements of a 20 m deep firnpit, applying a Newtonian model (Kesselwandferner, Oetztal Alps 1967-1978). Zeitschrift fiir Gletscherkunde und G1azialgeologie,20(2), in press. Mellor, M. (1975). A review of basic snow mechanics. International Association of Hydrological Sciences Publ. No. 114, Proceedings of the Grindelwald Symposium April 1974, pp. 251-291. Nye, J.F. (1965). The flow of a glacier in a channel of rectangular, elliptic or parabolic cross-section. J. Giaciol. 5(41): 661-690. Paterson, W.S.B. (1981). The Physics of Glaciers. 2nd Edition, Pergamon Press Oxford. Raymond, C.F. (1971). Determination of the three-dimensional velocity field in a glacier.J. Glaciol. 10(58): 3943 Salm, B. (1977). Snow forces. J. Glaciol. 19(81): 67-100.