UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Accurate modelling of glacier flow 1981

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
UBC_1982_A1 W34.pdf [ 20.48MB ]
Metadata
JSON: 1.0052458.json
JSON-LD: 1.0052458+ld.json
RDF/XML (Pretty): 1.0052458.xml
RDF/JSON: 1.0052458+rdf.json
Turtle: 1.0052458+rdf-turtle.txt
N-Triples: 1.0052458+rdf-ntriples.txt
Citation
1.0052458.ris

Full Text

ACCURATE MODELLING OF GLACIER FLOW by EDWIN DONALD WADDINGTON B . S c , U n i v e r s i t y of Toronto , 1971 M . S c , U n i v e r s i t y of A l b e r t a , 1973 THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n THE FACULTY OF GRADUATE STUDIES (Department of Geophysics and Astronomy) We accept t h i s t h e s i s as conforming to the r e q u i r e d s tandard THE UNIVERSITY OF BRITISH COLUMBIA November, 1981 © Edwin Donald Waddington, 1981 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e f o r reference and study. I further agree that permission for extensive copying of t h i s thesis for s c h o l a r l y purposes may be granted by the head of my department or by h i s or her representatives. I t i s understood that copying or p u b l i c a t i o n of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of Geophysics and Astronomy The University of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Da t e February 20, 1982, /7cn i i ABSTRACT Recent i n t e r e s t i n c l i m a t i c change and i c e .sheet v a r i a t i o n s p o i n t s out the need for accura te and n u m e r i c a l l y s t a b l e models of time-dependent i c e masses. L i t t l e a t t e n t i o n has been pa id to t h i s t o p i c by the g l a c i o l o g i c a l community, and there i s good reason to b e l i e v e tha t much of the p u b l i s h e d l i t e r a t u r e on numer ica l m o d e l l i n g of the f low of g l a c i e r s and i c e sheets i s q u a n t i t a t i v e l y i n c o r r e c t . In p a r t i c u l a r , the importance of the n o n l i n e a r i n s t a b i l i t y has not been w i d e l y r e c o g n i z e d . The purposes of t h i s t h e s i s are to develop and to v e r i f y a new numer ica l model for g l a c i e r f l o w , compare the model t o another w i d e l y accepted model, and to demonstrate the model i n s e v e r a l g l a c i o l o g i c a l l y i n t e r e s t i n g a p p l i c a t i o n s . As i n e a r l i e r work, the computer model so lve s the c o n t i n u i t y equat ion together w i t h a f low law for i c e . Thickness p r o f i l e s a long flow l i n e s are ob ta ined as a f u n c t i o n of time for a temperate i c e mass w i t h a r b i t r a r y bed topography and mass ba lance . A set of necessary t e s t s to be s a t i s f i e d by any numer ica l model of g l a c i e r f low i s p r e s e n t e d . The numer ica l s o l u t i o n s are compared w i t h a n a l y t i c a l s o l u t i o n s ; these i n c l u d e a s imple t h i c k n e s s - v e l o c i t y r e l a t i o n to check terminus m o b i l i t y , and Burgers* equat ion to check c o n t i n u i t y and dynamic behaviour w i t h f u l l n o n l i n e a r i t y . An attempt has been made to v e r i f y the accuracy of the computer model of Budd and Mclnnes (1974), Rudd (1975) and Mclnnes ( u n p u b l i s h e d ) . These authors have repor ted problems w i t h numer ica l i n s t a b i l i t y . I f the e x i s t i n g documentation i s i i i a c c u r a t e , the Budd-Mclnnes model appears to s u f f e r from mass c o n s e r v a t i o n v i o l a t i o n s both l o c a l l y and g l o b a l l y . The new numer ica l model developed i n t h i s t h e s i s can be used to r e c o n s t r u c t the v e l o c i t y f i e l d w i t h i n the g l a c i e r at each time s t ep ; t h i s v e l o c i t y f i e l d s a t i s f i e s c o n t i n u i t y and G l e n ' s f low law for i c e . I n t e g r a t i o n of t h i s v e l o c i t y f i e l d y i e l d s the t r a j e c t o r i e s of i n d i v i d u a l i ce elements f l o w i n g through the t i m e - v a r y i n g i ce mass. The t r a j e c t o r i e s and v e l o c i t y f i e l d are checked by comparison w i t h an a n a l y t i c a l s o l u t i o n for a steady s t a t e i c e sheet (Nagata, 1977). The model i n t h i s t h e s i s i s not r e s t r i c t e d to steady s t a t e , and i t avo ids the v i o l a t i o n s of mass c o n s e r v a t i o n , and the approx imat ions about the v e l o c i t y f i e l d found i n some p u b l i s h e d t r a j e c t o r y models . T h e ' f e a s i b i l i t y of us ing s t a b l e i sotopes to i n v e s t i g a t e p r e h i s t o r i c surg ing of v a l l e y g l a c i e r s has been s t u d i e d w i t h a model s i m u l a t i n g the S tee le G l a c i e r , Yukon T e r r i t o r y . A s l i d i n g < v e l o c i t y and surge d u r a t i o n were s p e c i f i e d , based on the o b s e r v a t i o n s of the 1966-67 surge . A surge p e r i o d of roughly 100 years gave the most r e a l i s t i c i c e t h i c k n e s s throughout the surge c y c l e . By c a l c u l a t i n g i c e t r a j e c t o r i e s and us ing two p l a u s i b l e r e l a t i o n s h i p s between 6 ( 0 1 B / 0 1 6 ) and p o s i t i o n or h e i g h t , l o n g i t u d i n a l s e c t i o n s and sur face p r o f i l e s of 6 were c o n s t r u c t e d for t imes be fo re , d u r i n g , and a f t e r a surge . D i s c o n t i n u i t i e s of up to 0 . 8 ° / O o were found acros s s e v e r a l sur faces d i p p i n g upstream i n t o the g l a c i e r . Each of these sur faces i s the present l o c a t i o n of the i ce which formed the i c e - a i r i n t e r f a c e at the time a p r e v i o u s surge began. I t may be d i f f i c u l t to observe these sur faces on the S t e e l e G l a c i e r due to i v the l a rge and poor ly -under s tood background v a r i a b i l i t y of 6. The genera t ion of wave og ives has been examined f o l l o w i n g the theory of Nye ( I 9 5 8 [ b ] ) r wherein waves are caused by a combinat ion of seasonal v a r i a t i o n i n mass balance and p l a s t i c deformat ion i n an i c e f a l l . The wave t r a i n generated on a g l a c i e r i s shown i n t h i s t h e s i s to be a c o n v o l u t i o n of the v e l o c i t y g r a d i e n t w i t h an i n t e g r a l of the mass balance f u n c t i o n . T h i s i n t e g r a l i s the impulse response of the g l a c i e r sur face to a s t ep i n the v e l o c i t y f u n c t i o n . S p a t i a l v a r i a t i o n s i n the g l a c i e r w i d t h and mass balance a l s o c o n t r i b u t e to the wave t r a i n . T h i s f o r m u l a t i o n i s used to e x p l a i n why many i c e f a l l s do not generate wave og ives i n s p i t e of l a rge seasonal balance v a r i a t i o n s and l a r g e p l a s t i c de format ions . V TABLE OF CONTENTS A b s t r a c t i i L i s t Of Ta b l e s x i v L i s t Of F i g u r e s xv Acknowledgements x i x Chapter 1: Beginnings 3 1.1 I n t r o d u c t i o n 3 1.1.1 Ends And Means 3 1.1.2 Conventions Used 9 1.2 Previo u s Work 11 1.2.1 Four C e n t u r i e s Of G l a c i e r Flow Theory 11 1.2.2 Previous Ice P r o f i l e Models 13 1.2.3 Previous Ice T r a j e c t o r y Models 15 1.3 The C o n t i n u i t y Equation For An Ice Mass 17 1.3.1 The General G l a c i e r Flow Problem 17 1.3.2 Rectangular C r o s s - s e c t i o n Flow Model 18 1.3.3 The C o n t i n u i t y Equation 20 1.3.4 P h y s i c a l I n t e r p r e t a t i o n 22 1.4 P h y s i c s Of Ice Deformation 23 1.4.1 I n t r o d u c t i o n 23 1.4.2 S t r e s s E q u i l i b r i u m Equations 24 1.4.3 C o n s t i t u t i v e R e l a t i o n For Ice Deformation 26 1.4.4 Shear S t r e s s And Ice Flux 30 1.5 B a s a l S l i d i n g 35 1.5.1 I n t r o d u c t i o n 35 1.5.2 Basal Ice Temperature And S l i d i n g 35 1.5.3 Philosophy Of S l i d i n g In T h i s Study 36 v i Chapter 2: Models And Tests ; 38 2.1 I n t r o d u c t i o n 38 2.1.1 O u t l i n e 38 2 .1 .2 Importance Of Model T e s t i n g 38 2.2 C o n t i n u i t y Equat ion P r o f i l e Model 42 2.2.1 I n t r o d u c t i o n 42 2 .2 .2 The Numerica l Scheme 43 2 . 2 . 3 Boundary C o n d i t i o n s 46 2 .2 .4 Numerica l S t a b i l i t y 48 2 .2 .5 Accuracy 55 2.3 T e s t i n g The C o n t i n u i t y Model 57 2.3.1 I n t r o d u c t i o n 57 2 .3 .2 C o n t i n u i t y Test With Terminus Motion 57 2 . 3 . 3 C o n t i n u i t y Test With Burger s ' Equat ion 62 2.4 The Ice T r a j e c t o r y Computer Model . .' 67 2.4.1 I n t r o d u c t i o n 67 2 .4 .2 The V e l o c i t y And Displacement F i e l d s 67 2 . 4 . 3 The Ice P a r t i c l e T r a j e c t o r i e s 70 2 .4 .4 Accuracy Of The T r a j e c t o r y Model 70 2.5 T e s t i n g The T r a j e c t o r y Model 74 2.5 .1 I n t r o d u c t i o n 74 2 .5 .2 Nagata Ice Sheet Test 74 2 . 5 . 3 Surface Mass Conserva t ion Test 79 Chapter 3: Can S tab le Isotopes Reveal A H i s t o r y Of Surging? 83 3.1 I n t r o d u c t i o n 83 3.2 S t e e l e G l a c i e r 85 3.2.1 General D e s c r i p t i o n 85 v i i 3.2.2 G l a c i e r Surges 86 3.2.3 Observations Of The S t e e l e G l a c i e r 88 3.2.4 P e r i o d Of S t e e l e Surges 91 3.3 Numerical Model 1 93 3.3.1 Flow Law Constants And Shape F a c t o r 93 3.3.2 Bed Topography 94 3.3.3 Channel Width 96 3.3.4 Mass Balance 97 3.3.5 C y c l i c Surge P a t t e r n For The Model 100 3.4 S t e e l e G l a c i e r Model 2 103 3.4.1 Problems With S t e e l e Model 1 103 3.4.2 S i m p l i f i c a t i o n s 108 3.5 S t a b l e Isotopes In G l a c i o l o g y 108 3.5.1 D e f i n i t i o n Of The Del Scale 108 3.5.2 F a c t o r s A f f e c t i n g Del 110 3.5.3 Previous I s o t o p i c S t u d i e s 112 3.5.4 Del R e l a t i o n s For The Model 116 3.6 Model R e s u l t s : Surge P e r i o d And T r a j e c t o r i e s 119 3.6.1 I n t r o d u c t i o n 119 3.6.2 P e r i o d i c a l l y Repeating S t a t e 119 3.6.3 Ice T r a j e c t o r i e s 126 3.7 Model R e s u l t s : D i s t r i b u t i o n Of Isotopes — 127 3.7.1 I n t r o d u c t i o n 127 3.7.2 Model 61: L o n g i t u d i n a l 6 S e c t i o n s 129 3.7.3 Model 61: Surface 6 P r o f i l e s 133 3.7.4 Model 62: S e c t i o n s And Surface P r o f i l e s 135 3.7.5 Are The P r e d i c t e d E f f e c t s Observable? 138 3.7.6 C o n c l u s i o n s 140 v i i i Chapter 4: Wave Ogives 142 4.1. I n t r o d u c t i o n 142 4.1.1 D e s c r i p t i o n Of Ogive Systems 142 4.1.2 Theories Of Wave Ogive Formation 145 4.1.3 Disappearance Of The Waves 146 4.2 Nye's Theory Of Wave Ogives 148 4.2.1 O u t l i n e Of The Nye Theory 148 4.2.2 Nye's Annually Repeating State S o l u t i o n 148 4.2.3 Ogives Above The F i r n L i n e 151 4.2.4 An Unanswered Question 152 4.3 A New S o l u t i o n For Ogives 153 4.3.1 Using Method Of C h a r a c t e r i s t i c s 153 4.3.2 Separable Mass Balance 155 4.3.3 The G e n e r a l i z e d V e l o c i t y 157 4.3.4 The Upstream Boundary C o n d i t i o n 157 4.3.5 The Terms Of The Flux S o l u t i o n 158 4.3.6 P h y s i c a l I n t e r p r e t a t i o n 160 4.3.7 The Green's F u n c t i o n For Ogives 167 4.3.8 A Co n v o l u t i o n Formulation For Ogives 167 4.4 Some Simple Examples 168 4.4.1 I n t r o d u c t i o n 168 4.4.2 Example: L i n e a r V e l o c i t y Gradient 169 4.4.3 Example: Double Step I c e f a l l Model 174 4.5 A u s t e r d a l s b r e e n 176 4.5.1 I n t r o d u c t i o n 176 4.5.2 Estimated Wave Generation 177 4.5.3 Ogive S o l u t i o n For A u s t e r d a l s b r e e n 178 4.5.4 F i n d i n g The Wave Generating Region 181 ix 4.6 C o n c l u s i o n s 184 L i s t Of Symbols 186 L i t e r a t u r e C i t e d 201 Appendix 1: C o n t i n u i t y Model 235 A1.1 The Numerical Scheme 235 A1.1.1 The C o n t i n u i t y Equation 235 A1 .1.2 A Ma t r i x Formulation 238 A1.2 N o n l i n e a r i t y 239 A1.3 Boundary C o n d i t i o n s 243 A 1.3.1 The Upper Boundary 243 A1.3.2 The Downstream Boundary 245 A1.3.3 Nonzero Flux Leaves Downstream Boundary 245 A1.3.4 Moving Wedge Terminus 247 A1.4 Numerical S t a b i l i t y 254 A1.4.1 I n t r o d u c t i o n 254 A1.4.2 The L i n e a r Computational I n s t a b i l i t y 255 A1.4.3 The-Nonlinear I n s t a b i l i t y 257 A1.4.4 V e l o c i t y Smoothing 260 A1.4.5 Numerical D i s s i p a t i o n 262 A1.4.6 D i s s i p a t i o n From The V e l o c i t y Equation 265 A1.4.7 Wavenumber S p e c t r a l T r u n c a t i o n 267 A1.5 Accuracy 269 A1.5.1 6 Parameter: Accuracy 269 A1.5.2 Phase E r r o r s 271 A1.5.3 T r u n c a t i o n E r r o r 274 Al.5.4 I n t e r p o l a t i o n E r r o r 278 Appendix 2: Ice T r a j e c t o r y Model 282 A2.1 I n t r o d u c t i o n 282 X A2.2 The V e l o c i t y F i e l d 284 A2.2.1 The Rectangular Flow Model 284 A2.2.2 The Downslope V e l o c i t y 285 A2.2.3 The L o n g i t u d i n a l S t r a i n Rate 286 A2.2.4 V e l o c i t y Normal To The Bed 288 A2.3 Ice Displacement F i e l d 290 A2.3.1 Four Point I n t e r p o l a t i o n 290 A2.3.2 Displacements At Meshpoints 291 A2.4 Ice P a r t i c l e T r a j e c t o r i e s 295 A2.4.1 T r a c k i n g Procedure 295 A2.4.2 P a r t i c l e s Which Reach Ice Surface 296 A2.4.3 T r a c k i n g Backwards In Time 296 A2.4.4 Boundary C o n d i t i o n At Upstream End 297 Appendix 3: Aspects Of D i s c r e t e Data S e r i e s 298 A3.1 The Z Transform 298 A3.2 A l i a s i n g 299 Appendix 4: Densit y Of G l a c i e r Ice 302 A4.1 F i r n As E q u i v a l e n t Ice Thickness 302 A4.2 Constant Density Assumption 303 Appendix 5: C o n t i n u i t y Equation For An Ice Mass 306 A5.1 Mass Conservation In A Moving Continuum 306 A5.2 In A S t a t i o n a r y G l a c i e r C r o s s - s e c t i o n 307 A5.3 In An A r b i t r a r y Channel 310 A5.4 In A Rectangular Channel 311 A5.5 In Bed-normal Coordinates 313 Appendix 6: Equations For P e r t u r b a t i o n s 316 Appendix 7: V e l o c i t y Equation For An Ice Mass 318 A7.1 I n t r o d u c t i o n 318 A7.2 The Shear S t r e s s Equation 319 A7.3 Approximations 322 A7.4 Shape F a c t o r s 328 A7.5 Shear S t r a i n Rate 330 A7.6 Ice Flux And Average V e l o c i t y 332 Appendix 8: G l a c i e r S l i d i n g 334 A8.1 Measurements . 334 A8.2 P h y s i c a l Processes In S l i d i n g 335 A8.3 Computer Models Of S l i d i n g 338 A8.3.1 Using Weertman S l i d i n g 338 A8.3.2 Budd-Mclnnes Model 338 A8.3.3 S l i d i n g In T h i s Study 342 Appendix 9: Burgers' Equation 343 Appendix 10: Matrix C o e f f i c i e n t s 346 Appendix 11: Convergence C r i t e r i a 350 Appendix 12: Machine Roundoff E r r o r s 354 Appendix 13: D i f f e r e n c i n g Scheme For The Flux Gradient .... 356 A13.1 I n t r o d u c t i o n 356 A13.2 P e r t u r b a t i o n Equations 357 A13.3 Space D i f f e r e n c i n g Scheme 358 A13.4 T r a n s f e r Function 360 A13.5 C o n d i t i o n s On The Mesh I n t e r v a l 361 Appendix 14: Ice Surface E l e v a t i o n s 364 Appendix 15: A n a l y t i c Models Of Ice Sheets 366 A15.1 I n t r o d u c t i o n 366 A15.2 Nye Ice Sheet Model 366 A15.3 Nagata Ice Sheet Model 367 A15.3.1 Basic Equations 367 x i i A l 5 . 3 . 2 Ice Depth, Mass Ba lance , And V e l o c i t y 368 A15 .3 .3 S t reaml ines 370 A15.4 H a e f e l i - P a t e r s o n Ice Sheet Model 374 Appendix 16: Tests Of The Budd-Mclnnes Model 375 A16.1 I n t r o d u c t i o n 375 A16.2 Ice Flow In The Budd-Mclnnes Model 376 A16.3 V a t n a j o k u l l Model : Non l inear I n s t a b i l i t y 377 A16.4 V a t n a j o k u l l (Model 1) : Mass Conserva t ion 381 A16.5 B r u a r j S k u l l (Model 2 ) : Steady State F lux 386 A16.6 Fedchenko G l a c i e r : Steady Sta te F lux 392 A16.7 Fedchenko G l a c i e r : N o n s l i d i n g Model 395 A16.8 Fedchenko G l a c i e r : Dynamic Behaviour 398 A16.9 Conc lus ions 411 Appendix 17.: Four C e n t u r i e s Of G l a c i e r Flow Theory . 414 A17.1 I n t r o d u c t i o n 414 A17.2 The Years 1570 To 1840 415 A17.2 .1 E a r l i e s t P ioneers 415 A17 .2 .2 H . B. De Saussure 416 A17 .2 .3 Rendu 418 A17.3 1840 To 1915 419 A17.3 .1 L o u i s Agas s i z 419 A17 .3 .2 J . D. Forbes 421 A 1 7 . 3 . 3 John T y n d a l l 424 A17 .3 .4 Many Wondrous Theor ie s 427 A 1 7 . 3 . 5 Ice Deformation Experiments 429 A17 .3 .6 Mathemat ica l G l a c i o l o g y . . . - 431 A17 .3 .7 The B r i n y Depths Of G l a c i e r s 434 A17.4 1915 To 1953 436 x i i i A17.4.1 I n t r o d u c t i o n 436 A17.4.2 Shear Plane S l i p Or V i s c o u s Flow? 437 A17.4.3 Continuum Mechanics For G l a c i e r s 439 A17.4.4 The G l a c i e r A n t i c y c l o n e 441 A17.4.5 F i e l d S t u d i e s : 1934 Spit z b e r g e n E x p e d i t i o n 444 A17.4.6 F i e l d S t u d i e s : J u n g f r a u j o c h Research Party .... 446 A17.4.7 E x t r u s i o n Flow 449 Appendix 18: S t a b i l i t y C o n d i t i o n For A Surge Bulge 454 Appendix 19: S t e e l e G l a c i e r T r i b u t a r i e s 457 xiv LIST OF TABLES 2.1. Parameters For Vatnajokull (Figure 2.3) 53 2.2. Parameters For Test With Moving Terminus 58 2.3. Model Parameters For Burgers' Equation Test 66 2.4. Parameters For Nagata Ice Sheet 78 2.5. Residence Times In Nagata Ice Sheet 80 2.6. Nagata Ice Sheet Surface Boundary Test 82 3.1. Parameters For Steele Steady State 104 3.2. Velocity Pattern For Steele Surge 120 3.3. Numerical Time Steps For Steele Surge 124 4.1. Odinsbreen: Linear Velocity Approximations 177 A1 5.1 . Parameters For Nagata Ice Sheet 374 A16.1. Parameters For Vatnajokull Model 379 A16.2. Parameters For Fedchenko Model 398 A16.3. Parameters For S l i d i n g Models 405 A19.1. Ice Flux From Steele Glacier Tributaries (a) 457 A19.2. Ice Flux From Steele Glacier Tributaries (b) 459 A19.3. T r i b u t a r i e s : Effect On Mass Balance 460 X V LIST OF FIGURES F r o n t i s p i e c e .. 2 1.1. Rectangular C r o s s - s e c t i o n Flow 20 1.2. V e r t i c a l Prism For C o n t i n u i t y I n t e r p r e t a t i o n 23 2.1. Numerical Scheme At Ice D i v i d e 47 2.2. The Wedge Terminus 48 2.3. Example Of The Nonlinear I n s t a b i l i t y 52 2.4. F i l t e r To Suppress The N o n l i n e a r I n s t a b i l i t y 55 2.5. C o n t i n u i t y Test with Moving Terminus 59 2.6. C o n t i n u i t y Test With Moving Terminus 61 2.7. N o n l i n e a r Test With Burgers' Equation 65 2.8. Meshpoints For Ice V e l o c i t y C a l c u l a t i o n s 68 2.9. Nagata Steady Ice Sheet 75 2.10. Growth Of Nagata Ice Sheet 76 2.11. V e l o c i t y F i e l d For Nagata Model 77 2.12. T r a j e c t o r i e s In Nagata Model 79 2.13. Sur f a c e Mass Cons e r v a t i o n Test 81 3.1. I c e f i e l d Ranges L o c a t i o n Map 86 3.2. S t e e l e G l a c i e r And T r i b u t a r i e s 89 3.3. Model 1 For S t e e l e G l a c i e r 98 3.4. S l i d i n g Model For S t e e l e G l a c i e r 101 3.5. S t e e l e Model 1 Growth To Steady State 104 3.6. Steady S t a t e S t r e a m l i n e s For Model 1 105 3.7. Model 2 For S t e e l e G l a c i e r 106 3.8. Reference Surface For 6-x F u n c t i o n 118 3.9. S l i d i n g V e l o c i t y : S t e e l e G l a c i e r Model 121 3.10. Pre- And Post-surge P r o f i l e s : 47 Year P e r i o d 122 x v i 3 .11 . P r e - And Post-surge P r o f i l e s : 97 Year P e r i o d 123 3 .12 . Pre - And Post-surge P r o f i l e s : 147 Year P e r i o d 124 3 .13 . S tee le G l a c i e r T h i c k n e s s : One Surge Cyc le 125 3 .14 . Ice T r a j e c t o r i e s For 97 Year Surge P e r i o d 126 3 .15 . L o n g i t u d i n a l 6 S e c t i o n s : Model 61 130 3 .16 . Model 61 : Surface 6 P r o f i l e s 134 3 .17 . Model 62: L o n g i t u d i n a l 6 Sec t ions 136 3 .18 . Model 62: Surface 6 P r o f i l e s 137 4 . 1 . Aus terda l sbreen V e l o c i t y And Mass Balance 151 4 . 2 . The C h a r a c t e r i s t i c s In T- t Space 154 4 . 3 . Double Step I c e f a l l Model 162 4 . 4 . S i n g l e V e l o c i t y Step Model 164 4 . 5 . Three Fac tor s Generat ing Waves 166 4 . 6 . Ogives From A V e l o c i t y Grad ient 170 4 . 7 . Flow Past A V e l o c i t y G r a d i e n t : Numerica l S o l u t i o n . . . . 173 4 . 8 . Double Step I c e f a l l Model 174 4 . 9 . Odinsbreen: G e n e r a l i z e d V e l o c i t y Per U n i t Width 176 4 .10 . Aus te rda l sbreen Wave Ogives 179 4 . 1 1 . Aus te rda l sbreen Ice Thickness 180 4 .12 . V a r i a t i o n s On Odinsbreen I c e f a l l 183 A 1 . 1 . Mesh Increment On Bed 236 A 1 . 2 . Model Terminus 248 A 1 . 3 . A l i a s i n g And The Nonl inear I n s t a b i l i t y 260 A 1 . 4 . Trans fer Funct ions Of Smoothing Schemes 261 A 1 . 5 . Trans fer F u n c t i o n s : Slope-dependent Damping 266 A 1 . 6 . F i l t e r To Suppress N o n l i n e a r I n s t a b i l i t y 267 A 1 . 7 . Trans fer F u n c t i o n Modulus For V a r i o u s 6 270 A 1 . 8 . Trans fer F u n c t i o n Phase Comparison 272 x v i i A 2 . 1 . Mesh For Ice Displacement C a l c u l a t i o n s 283 A 2 . 2 . The Rectangular Flow Model 284 A 2 . 3 . Four P o i n t I n t e r p o l a t i o n Scheme 290 A 2 . 4 . I n t e r p o l a t i o n Surface f (P ) 292 A 2 . 5 . C e l l Ver tex N o t a t i o n For Downward V e l o c i t y 294 A 2 . 6 . Displacement F i e l d In A Steady State 295 A 2 . 7 . C e l l Ver tex N o t a t i o n For Negat ive Time 297 A3.1 . The Z Plane 299 A 3 . 2 . S i g n a l s With Wavelengths 1.5Ax And 3 . 0AX 300 A 4 . 1 . Force Balance On An Ice Element 303 A 5 . 1 . Surfaces for D e r i v a t i o n of C o n t i n u i t y Equat ion 307 A 5 . 2 . The Thin C r o s s - s e c t i o n L i m i t 309 A 5 . 3 . Coord inates And V a r i a b l e s In Rectangular Channel . . . . 313 A 1 0 . 1 . Q u a n t i t i e s In Slope C a l c u l a t i o n 348 A 1 3 . 1 . Damping Using The Ice Surface Slope 359 A14.1 . Ice Surface E l e v a t i o n 364 A 1 5 . 1 . Nagata Steady Ice Sheet 373 A 1 6 . 1 . V a t n a j o k u l l : I n s t a b i l i t y And Growth Rate 378 A 1 6 . 2 . G l a c i e r Mass As A F u n c t i o n Of Time 383 A 1 6 . 3 . Rate Of Growth As A F u n c t i o n Of Length L 384 A 1 6 . 4 . B r u a r j d k u l l Ice P r o f i l e s (Model 2) 387 A 1 6 . 5 . B r u a r j d k u l l F lux Test 389 A 1 6 . 6 . Fedchenko G l a c i e r F l u x Test 393 A 1 6 . 7 . Fedchenko Steady Sta te Ice P r o f i l e s 396 A 1 6 . 8 . Fedchenko N o n s l i d i n g Model 397 A 1 6 . 9 . Fedchenko n o n s l i d i n g model w i t h n=3 399 A16 .10 . Growth Of Fedchenko G l a c i e r 400 A 1 6 . 1 1 . Fedchenko Growth: Other N o n s l i d i n g Models 402 x v i i i A16.12 . Fedchenko Growth: Moderate S l i d i n g 404 A16 .13 . Fedchenko G l a c i e r : 0-V Plane 406 A16 .14 . Fedchenko Growth: S l i d i n g Model (2) 407 A16 . 1 5 . Fedchenko Growth: S l i d i n g Model (3) 409 A16.16. Fedchenko Growth: S l i d i n g Model (4) 410 A 1 8 . 1 . Advancing Surge Bulge 454 x i x ACKNOWLEDGEMENTS In the e a r l y stages of t h i s work, encouragement and support from my a d v i s o r G. K. C. C l a r k e helped me keep my research g o i n g , e s p e c i a l l y dur ing the per iods when the computer programs were uncoopera t ive . H i s enthusiasm i s contag ious , and i t i s d e l i g h t f u l to work w i t h h im. My a p p r e c i a t i o n of him grows a long w i t h my knowledge of g l a c i o l o g y . D i s c u s s i o n s w i t h W. S. B. P a t e r s o n , B. B. Narod, W. H . Mathews, D. W. Oldenburg, C. F . Raymond, R. D. R u s s e l l , M. C. Quick, R. A . B i n d s c h a d l e r , and W. D. H i b l e r I I I have shed l i g h t on aspects of g l a c i e r dynamics and numer ica l methods. Dur ing the l a t e stages of t h i s p r o j e c t , P . K. F u l l a g a r , complet ing h i s P h . D . t h e s i s i n the adjacent o f f i c e , has been a frequent companion on " the n ight s h i f t " exchanging ideas and mutual encouragement. B. B. Narod put i n long hours p roo f read ing the manuscript and o f f e r i n g h e l p f u l sugges t ions , and he, together w i t h a group of co l l eagues and f r i e n d s , helped me through that l a s t h e c t i c n ight when the manuscript was produced. J u l i a Forbes , through her encouragement, i n t e r e s t , and suppor t , has helped me to formulate and to work toward my p e r s o n a l g o a l s . I have made good f r i e n d s at U . B . C . In p a r t i c u l a r , Barry Narod, J . G . N a p o l e o n i , Bo Chandra, and Peter F u l l a g a r have shared many happy non-academic exper iences w i t h me. I was supported i n part by a N a t i o n a l Research C o u n c i l of Canada postgraduate s c h o l a r s h i p , and by an H. R. M a c M i l l a n Fami ly F e l l o w s h i p from the U n i v e r s i t y of B r i t i s h Columbia . The computations have been c a r r i e d out at the Computing Centre at U . B . C . A u s t i n P o s t , of the U . S. G e o l o g i c a l Survey, p rov ided the f r o n t i s p i e c e photograph of the Tr imble G l a c i e r . FRONTISPIECE: Wave og ives on the North Branch, Tr imble G l a c i e r , Ala ska Range, 6 1 ° 4 0 ' N , 152°18'W. Photo by A u s t i n P o s t , U . S. G e o l o g i c a l Survey, 1965. 2 3 CHAPTER j_ : BEGINNINGS ' " I t s the job t h a t ' s never s t a r t e d as takes longest to f i n i s h " as my o l d ga f fer used to s a y . ' 1 1.1 INTRODUCTION 1.1.1 ENDS AND MEANS Al though g l a c i e r s and i c e sheets o f ten appear to be far-removed from most day-to-day mat te r s , there are a number of c o m p e l l i n g reasons to study t h e i r behav iour . Advances of some g l a c i e r s would threa ten roads , dams and mines. Berendon G l a c i e r , B r i t i s h Columbia was s t u d i e d by U n t e r s t e i n e r and Nye (1968) and by F i s h e r and Jones (1971) fo r t h i s reason. Ice avalanches from g l a c i e r s have caused a long h i s t o r y of d e s t r u c t i o n . For example, a s e r i e s of four i c e avalanches from the Randa G l a c i e r , S w i t z e r l a n d ( A g a s s i z , 1840, p . 158) between 1636 and 1819 des t royed many b u i l d i n g s and f i e l d s and k i l l e d dozens of c i t i z e n s . Ice avalanches i n t o moraine-dammed lakes i n Peru caused damaging f loods ( L l i b o u t r y and o t h e r s , 1977). Advancing g l a c i e r s can dam streams or r i v e r s ; the r e s u l t i n g l a k e s o f ten d r a i n c a t a s t r o p h i c a l l y (Clarke and Mathews 1981; C l a r k e , i n press ) when the i c e dam f a i l s . Cunningham (1854 ( r e p r i n t e d 1970), p . 100) repor ted damaging f l o o d s on the Indus R i v e r i n the n ineteennth c e n t u r y , and Forbes (1845, p . 262) d e s c r i b e d the 1818 d i s a s t e r when the Getroz g l a c i e r dammed the 1 Sam Gamgee, i n The F e l l o w s h i p of the R i n g . J. R. R. T o l k i e n . 4 Dranse i n the V a l , de Bagnes i n S w i t z e r l a n d . Severa l glacier-dammed lakes threa ten a proposed p i p e l i n e . route i n the Yukon T e r r i t o r y (Canada, u n p u b l i s h e d ) . There i s s t i l l no u n i v e r s a l l y accepted theory on the cause of i ce ages and c o n t i n e n t a l g l a c i a t i o n ; an unders tanding of i ce sheet dynamics helps to s e l e c t and t e s t hypotheses . To c o r r e c t l y i n t e r p r e t the geomorphological record of P. le is tocene i c e sheets , we must understand the processes of g l a c i a l e r o s i o n and d e p o s i t i o n . Thi s r e q u i r e s a knowledge of g l a c i e r mechanics (e.c>. B o u l t o n , 1979; H a l l e t , 1979). The volume, d i s t r i b u t i o n , and rate of growth and decay of the P l e i s t o c e n e i c e sheets (§_.£. P a t e r s o n , 1972) are important data for the d e t e r m i n a t i o n of the v i s c o s i t y of the upper mant le , v e r t i c a l c r u s t a l movements, and sea l e v e l changes ( e . £ . Andrews, 1974). The i s o t o p i c compos i t ion of p o l a r i c e sheets has been used to r e c o n s t r u c t temperature changes and c l i m a t e over the past 10 s years (Dansgaard and o t h e r s , 1969). To c o r r e c t l y date deep c o r e s , i t i s necessary to determine the flow p a t t e r n w i t h i n the i c e sheet (Dansgaard and Johnsen, I969[a ] ; P h i l b e r t h and Federer , 1971; Hammer and o t h e r s , 1978). A c u r r e n t ques t ion of some concern i s the p o s s i b i l i t y of g l o b a l atmospheric warming due to combustion of f o s s i l f u e l s and c l e a r i n g of temperate f o r e s t s (e«g_« SMIC, 1971). A m u l t i d i s c i p l i n a r y study (NOAA, unpubl i shed) i s underway i n B o u l d e r , Colorado to i n v e s t i g a t e the e f f e c t increased atmospheric C0 2 would have on the A n t a r c t i c i c e sheets . D i s i n t e g r a t i o n of the East A n t a r c t i c Ice Sheet c o u l d r a i s e sea l e v e l by 75 m and s u b s t a n t i a l l y reduce the albedo of the ea r th 5 ( W i l s o n , 1969). Of more immediate concern i s the p o s s i b i l i t y of a surge and d i s i n t e g r a t i o n of the West A n t a r c t i c i c e sheet ; t h i s c o u l d r a i s e sea l e v e l by seven metres i n l e s s than 100 years (Thomas and o t h e r s , 1979). A group at NASA (NASA, unpubl ished) i s u s ing s a t e l l i t e r ad iometry , a l t i m e t r y and radar imaging to monitor and to h e l p model v a r i a t i o n s of the Greenland i c e sheet . Nye (1951, I952[a ] , 1953, 1957) made the f i r s t q u a n t i t a t i v e s tud ie s of the steady flow of g l a c i e r s and i ce sheets us ing a n a l y t i c a l models, and Weertman (1958), L l i b o u t r y ( I958 [b ] ) , and Nye (1960, 1961, I963[a ] , I963[b] , 1963[C ] , I965[a ] , I965[b]) developed the theory of g l a c i e r v a r i a t i o n s , k inemat ic waves, and response to c l i m a t e , by us ing p e r t u r b a t i o n methods. Many i n t e r e s t i n g g l a c i o l o g i c a l problems have l a rge temporal v a r i a t i o n s or compl ica ted boundary c o n d i t i o n s ; the a n a l y t i c a l s o l u t i o n s cannot be used. Answers to some of these more compl i ca ted problems can be found by numer ica l methods us ing f i n i t e d i f f e r e n c e s on d i g i t a l computers. However, numer ica l s o l u t i o n s have t h e i r own s p e c i a l p i t f a l l s . A numer ica l s o l u t i o n of a d i f f e r e n t i a l equat ion may d i f f e r from the c o r r e c t s o l u t i o n for many reasons (e.<j. Richtmyer and Morton, 1967; Gary, 1975). I t i s extremely d i f f i c u l t to prove that a numer ica l model has c o r r e c t l y so lved a p a r t i c u l a r d i f f e r e n t i a l equat ion w i t h c o m p l i c a t e d boundaries i f no a n a l y t i c a l check i s a v a i l a b l e ; yet t h i s i s p r e c i s e l y the type of problem for which numer ica l models are neces sary . I t i s e s s e n t i a l to f i r s t check numer ica l models aga ins t a n a l y t i c a l s o l u t i o n s fo r a v a r i e t y of s imple r problems. The g l a c i o l o g i c a l l i t e r a t u r e c o n t a i n s very l i t t l e d i s c u s s i o n of 6 model v e r i f i c a t i o n i n s p i t e of i t s obvious importance. There i s some i n d i c a t i o n t h a t much of the p u b l i s h e d l i t e r a t u r e on numerical m o d e l l i n g may be q u a n t i t a t i v e l y i n c o r r e c t . The major t h r u s t of my work has been aimed at understanding the problems of numerical models, f i n d i n g ways to a v o i d the problems, and d e v i s i n g t e s t s to v e r i f y the accuracy of the model r e s u l t s . With t h i s i n mind, I have developed a new computer model of g l a c i e r flow (Appendix 1, and Chapter 2, S e c t i o n 2.2). L i k e s e v e r a l p r e v i o u s models (Budd and Jenssen, 1975; Bi n d s c h a d l e r , u n p u b l i s h e d ) , t h i s model uses f i n i t e d i f f e r e n c e s to s o l v e the mass c o n s e r v a t i o n equation together with a flow law fo r i c e , to g i v e the time-dependent g l a c i e r s u r f a c e f o r a temperate i c e mass i n a channel of a r b i t r a r y width and bed topography, with an a r b i t r a r y mass balance, assuming the flow i s d r i v e n by g r a v i t a t i o n a l s t r e s s e s . I have analyzed the numerical s t a b i l i t y and accuracy of t h i s model as thoroughly as i s p o s s i b l e f o r n o n l i n e a r equations (Appendix 1). In Chapter 2 I present a set of t e s t s comparing the numerical s o l u t i o n s to a n a l y t i c a l s o l u t i o n s t o check terminus m o b i l i t y and both l o c a l and g l o b a l mass c o n s e r v a t i o n , i n c l u d i n g a case with a n o n l i n e a r flow law. I f , i n a computer model, the g l a c i e r terminus moves i n c o r r e c t l y , i t can s e r i o u s l y a f f e c t the i c e t h i c k n e s s and the v e l o c i t y throughout the g l a c i e r ( S e c t i o n A1.3.4). The p h y s i c s of the deformation of a g l a c i e r snout i s complicated (Nye, 1967) f o r any r e a l i s t i c i c e rheology. In the standard numerical approximation ( e . S . Budd and Jenssen, 1975; B i n d s c h a d l e r , unpublished, p. 105), the terminus i s simply a wedge-shaped 7 volume w i t h s lope and apex chosen so as to conserve mass (see S e c t i o n A 1 . 3 . 4 ) . The e r r o r i n u s ing t h i s k inemat ic approximat ion cannot be determined? the c o r r e c t genera l s o l u t i o n for the motion of a g l a c i e r terminus on an a r b i t r a r y s lope w i t h G l e n ' s f low i n tensor form i s s t i l l an unsolved problem. However, I have t e s t ed the numerica l implementat ion of the wedge terminus by comparing the computed s o l u t i o n to a time-dependent a n a l y t i c a l s o l u t i o n w i t h a s impler " r h e o l o g y " (Sec t ion 2 . 3 . 2 ) . Many standard numer ica l schemes for l i n e a r equat ions break down when a p p l i e d to n o n l i n e a r e q u a t i o n s . I t i s important to t e s t a numer ica l model w i t h a n o n l i n e a r problem. Burgers ' equat ion (Sec t ion 2 .3 .3) i s a n o n l i n e a r h y p e r b o l i c equat ion w i t h an a n a l y t i c a l s o l u t i o n ; i t i s a l s o r e l a t e d to the mass c o n s e r v a t i o n e q u a t i o n . I have compared numer ica l r e s u l t s w i t h the a n a l y t i c a l s o l u t i o n to Burgers ' equat ion to show that my model c o r r e c t l y so lves n o n l i n e a r problems. For some g l a c i e r f low problems, such as d a t i n g i c e cores (Dansgaard and Johnsen, 1969[a]) and f i n d i n g the temperature d i s t r i b u t i o n of c o l d i ce masses ( e . £ . Jenssen, 1977), i t i s necessary to know the t r a j e c t o r i e s of i c e p a r t i c l e s . My computer model can c a l c u l a t e the v e l o c i t y f i e l d on a v e r t i c a l l o n g i t u d i n a l mesh for a time-dependent g l a c i e r , by us ing G l e n ' s f low law to f i n d the h o r i z o n t a l v e l o c i t y , and us ing the c o n t i n u i t y equat ion to then d e r i v e the v e r t i c a l v e l o c i t y . P a r t i c l e t r a j e c t o r i e s are found by a numer ica l i n t e g r a t i o n of the time-dependent v e l o c i t y . I t e s t e d t h i s par t of the computer model aga in s t an a n a l y t i c a l s o l u t i o n by Nagata (1977) fo r p a r t i c l e paths i n a steady i c e sheet . 8 T h i s numer ica l model i s probably the most thoroughly and a c c u r a t e l y t e s t e d of i t s type . The set of t e s t s which I have assembled, or o thers s i m i l a r to them, should be used to v e r i f y any numer ica l model of g l a c i e r f l o w . Only then can the models be used w i t h conf idence to so lve more compl i ca ted problems. I have used t h i s new computer model i n two s t u d i e s . Prev ious e f f o r t s have concentra ted on us ing v a r i a t i o n s i n i s o t o p i c r a t i o s i n i c e cores to i n v e s t i g a t e c l i m a t e change, assuming steady s ta te f l o w . In Chapter 3 I have eva lua ted the f e a s i b i l i t y of us ing s t a b l e i so tope measurements to study the surge h i s t o r y of v a l l e y g l a c i e r s , assuming a constant c l i m a t e and unsteady f l o w . The example I cons idered was the S tee le G l a c i e r , Yukon T e r r i t o r y . I found that surg ing leaves a c h a r a c t e r i s t i c p a t t e r n i n the i so tope d i s t r i b u t i o n , but p r e l i m i n a r y measurements of 6 ( 0 1 8 / 0 1 6 ) suggest that t h i s p a t t e r n may be masked by other e f f e c t s . F i n a l l y , i n Chapter 4 I have d e r i v e d a l i n e a r c o n v o l u t i o n r e l a t i n g the ampli tude of wave og ives to the v e l o c i t y , channel wid th and mass balance i n i c e f a l l s . T h i s work i s an ex tens ion of s t u d i e s by Nye ( I 9 5 8 [ b ] ) . I used the computer model to v e r i f y the c o n v o l u t i o n f o r m u l a t i o n and to determine which fea tures of the i c e f a l l on A u s t e r d a l s b r e e n , Norway are most c r i t i c a l to the f o r m u l a t i o n of i t s l a rge wave o g i v e s . 9 1.1.2 CONVENTIONS USED There are many d i v e r s e views on the most a p p r o p r i a t e s t y l e fo r a P h . D . t h e s i s . My aim has been to produce a document which f u l l y d e s c r i b e s my work, and which can be understood on i t s own by those w i t h a bas i c knowledge of p h y s i c s or p h y s i c a l g l a c i o l o g y . I have documented a l l my numer ica l methods i n d e t a i l , and summarized the r e l e v a n t work of o t h e r s ; re ferences s u b s t a n t i a t e the t ex t r a ther than s u b s t i t u t e for i t . T h i s r e s u l t s i n a lengthy manuscr ip t . To keep the main t e x t as short as p o s s i b l e , I have p laced the numer ica l methods and the background m a t e r i a l i n appendices . The work of o ther s should be c l e a r l y i d e n t i f i a b l e . I hope that t h i s l e v e l of d e t a i l w i l l be a p p r e c i a t e d by some readers , s ince b r e v i t y w i l l be r e q u i r e d i n the v e r s i o n of t h i s work i n p r e p a r a t i o n for p u b l i c a t i o n . Unless s t a t ed o t h e r w i s e , I have used a r ighthanded l o c a l l y orthonormal coord ina te system such that the x a x i s l i e s a long the g l a c i e r bed down the c e n t r e l i n e of the c h a n n e l . The y a x i s i s t r ansver se and h o r i z o n t a l , and the z a x i s i s normal to the bed and p o s i t i v e upward i n the v e r t i c a l p lane c o n t a i n i n g the c e n t r e l i n e . The v e l o c i t y components are (u ,w,v) a long the ( x , y , z ) axes . T h i s n o t a t i o n d i f f e r s from the s tandard convent ion ( i . . e . (u ,v ,w) ) due to h i s t o r i c a l developments i n the t h e s i s . Underscores are used t o i n d i c a t e t e n s o r s . The number of underscores i n d i c a t e s the rank of the t e n s o r , e . £ . v i s the v e l o c i t y v e c t o r , and A i s a c o e f f i c i e n t m a t r i x . The dot symbol when l o c a t e d above a v a r i a b l e , i n d i c a t e s 10 a time d e r i v a t i v e . When l o c a t e d between two v e c t o r s , i t represents the standard s c a l a r inner product or dot p roduct , e.cj . Malvern ( 1969, p . 17) . A bar above a v a r i a b l e i n d i c a t e s an average v a l u e , u s u a l l y over a depth range (z d i r e c t i o n ) , or over time (§_.£. annual averages ) . A l i s t of symbols, together w i t h t h e i r meanings and the s e c t i o n i n which each i s i n t r o d u c e d , can be found f o l l o w i n g Chapter 4. Equat ion numbers, and t e x t u a l re ferences to equat ion numbers, are enclosed in round parentheses , e . £ . ( 2 . 2 . 5 ) , or ( A 5 . 6 ) , or ( A 1 . 1 . 3 ) . The c h a r a c t e r s preceding the f i r s t decimal p o i n t are the chapter or appendix number. The middle number ( i f present) i d e n t i f i e s the chapter subsec t ion where the equat ion i s g i v e n , and the f i n a l number i s the consecut ive equat ion number i n that s u b s e c t i o n . References to chapter s e c t i o n s are always i d e n t i f i e d as such, and the numbers are not enc losed i n parentheses . The LITERATURE CITED i s i n the s t y l e of The J o u r n a l of G l a c i o l o g y . I was not ab le to o b t a i n the use of some of the very e a r l y l i t e r a t u r e , and some of the l i t e r a t u r e i n languages other than E n g l i s h . In those cases , where I c o u l d not v e r i f y the c i t a t i o n s of o t h e r s , I have i n c l u d e d the c i t i n g author i n square b r a c k e t s . 11 1.2 PREVIOUS WORK 1.2.1 FOUR CENTURIES OF GLACIER FLOW THEORY The framework w i t h i n which we c u r r e n t l y understand and i n v e s t i g a t e g l a c i e r flow has been assembled i n the past three decades. Deep c o r i n g techniques (Hansen and Langway, 1966) and rad io-echo sounders (Evans, 1963), coming s h o r t l y a f t e r important experiments on the deformat ion of i c e (Glen , 1952, 1955) and mathematical treatment of the flow of i c e sheets and g l a c i e r s (Nye, 1951, I952[a] , 1957) s t a r t e d the r a p i d growth of g l a c i o l o g i c a l r e s e a r c h . However, i n v e s t i g a t i o n of the f low of g l a c i e r s goes back hundreds of y e a r s . In Appendix 17, I review the development of ideas on g l a c i e r f low i n that e a r l y p e r i o d . Some e a r l y works were based on f e r t i l e imag inat ion and l i m i t e d o b s e r v a t i o n s ; o ther s were conc i s e and l u c i d statements on t o p i c s which are sub jec t s of research today . Some misconcept ions about g l a c i e r f low were r a i s e d , debated, and re so lved more than once d u r i n g the p e r i o d . The main t h r u s t of research d u r i n g the second h a l f of the n ine teenth century was d i r e c t e d by p h y s i c i s t s ; i n the e a r l y par t of t h i s c e n t u r y , g e o l o g i s t s dominated the f i e l d ( w i t h a few notab le e x c e p t i o n s ) , and research p r i o r i t i e s and t h e o r i e s r e f l e c t e d t h i s d i f f e r e n c e . Contemporary reviews of g l a c i e r f low theory were g iven by C r o l l (1875, Chapter XXX, p . 495) , by G e i k i e (1894, Chapter 3, p . 25) , by R u s s e l l (1897, Chapter 9, p . 160), by Hawkes (1930), by Matthes (1942), by P eru tz (1947) , and by Orowan (1949) . S ince 1950, research on g l a c i e r flow has progressed 12 r a p i d l y . The flow law for i c e was e s t a b l i s h e d for many p r a c t i c a l purposes by Glen (1952, 1955) and Nye (1953). Weertman (1957), L l i b o u t r y d 968[a] , I968[b] ) , Kamb (1970), Nye d 9 6 9 [ b ] , 1970) and Morland d 976 [ a ] , 1976[b3) c o n t r i b u t e d to the theory of g l a c i e r s l i d i n g ; some aspects of t h i s que s t ion are s t i l l u n r e s o l v e d . Papers by Nye (1951, I952[a ] , I952[b] , 1 9 5 2 [ C ] , 1957, 1959[C ] ) e s t a b l i s h e d r e a l i s t i c a n a l y t i c a l s o l u t i o n s for steady g l a c i e r and i c e sheet p r o f i l e s , v e l o c i t i e s , and s t r e s s f i e l d s , w h i l e i d e n t i f y i n g many u s e f u l approx imat ions . V i a l o v (1957) used G l e n ' s f low law to d e r i v e a steady i c e sheet p r o f i l e which matched the flow l i n e through M i r n y , A n t a r c t i c a . Weertman (1961[b]) examined the e f f e c t s of l o n g i t u d i n a l s t r a i n ra tes on steady i c e sheet p r o f i l e s , and i n c l u d e d i s o s t a t i c depress ion of the bed. Weertman (1963) cons idered the e f f e c t s of f r i n g i n g mountain ranges on steady i c e sheet s , and (1966) the e f f e c t of a basa l water l a y e r . Al though temporal v a r i a t i o n s of i c e masses are d i f f i c u l t to study f u l l y , some u s e f u l a n a l y t i c a l r e s u l t s have been d e r i v e d . Kinemat ic waves on g l a c i e r s were observed by V a l l o t (1900) and were s t u d i e d by L l i b o u t r y d 9 5 8 [ b ] ) , Nye (1958[a] ) f and by Weertman (1958) us ing p e r t u r b a t i o n methods. Nye, i n a s e r i e s of papers , (1960, 1961, I963[a ] , I963[b] , 1963[C ] , I965[a ] , I965[b]) extended the method to analyze the response of g l a c i e r s to c l i m a t i c change, and to es t imate past c l i m a t e from the record of advance and r e t r e a t of g l a c i e r s . Bodvarsson (1955) d e r i v e d equat ions for a t h i n i c e sheet and analyzed i t s s t a b i l i t y to c l i m a t i c change. T h i s model i s not w i d e l y used due to i t s assumed r e l a t i o n between basa l s t r e s s and i c e f l u x . Weertman 13 (1961[a]) performed a s i m i l a r s t a b i l i t y a n a l y s i s assuming Weertman (1957) s l i d i n g . He a l s o d e r i v e d (I964[a]) the time s ca le s for the growth or decay of a p e r f e c t l y p l a s t i c i c e sheet . Jenssen and Radok (1963) obta ined a numer ica l s o l u t i o n for the temperature f i e l d i n the c e n t r a l reg ion of an i ce sheet undergoing t h i n n i n g . The study of f u l l y time dependent i c e masses w i t h a r b i t r a r y boundaries and source terms i s a recent development made p o s s i b l e by h igh speed computers. Numerical s o l u t i o n s of the equat ions governing i c e masses have been obta ined for a range of problems by, e.g_. Shumskiy (1963), Budd and Jenssen (1975), Mahaffy (1976), Jenssen (1977), B indschad ler (unpublished) and C l a r k e (1976). The complete s o l u t i o n of the equat ions of mot ion , the c o n s t i t u t i v e equat ions and the equat ions of s t a t e for a time v a r y i n g i c e mass w i t h a r b i t r a r y sources and boundaries i s an out s t and ing problem. 1.2.2 PREVIOUS ICE PROFILE MODELS Computer models which f i n d the sur face he ight of t i m e - v a r y i n g g l a c i e r s and i c e sheets are a r e l a t i v e l y new research t o o l . Campbell and Rasmussen (1969, 1970) and Rasmussen and Campbell (1973) developed a model which found i c e depth at p o i n t s on a h o r i z o n t a l mesh. They assumed tha t the i c e was a v i s c o u s mater ia l- w i t h a ba sa l f r i c t i o n c o e f f i c i e n t determined by mass f l u x . By a r b i t r a r i l y l ower ing the ba sa l f r i c t i o n c o e f f i c i e n t they s imula ted g l a c i e r surges . 14 Budd and others (1971) and Budd and Jenssen (1975) developed a f i n i t e d i f f e r e n c e model to so lve the c o n t i n u i t y equat ion for the g l a c i e r t h i c k n e s s p r o f i l e a long a f l o w l i n e . Whi le the equat ions are s i m i l a r to those I have used, our numer ica l methods d i f f e r i n some respects (see Appendix 16). These authors have i n c l u d e d a d d i t i o n a l important p h y s i c a l p r o p e r t i e s of i c e f low (e.g_. the e f f e c t of l o n g i t u d i n a l s t r e s s d e v i a t o r s ) and I have devoted more e f f o r t toward o b t a i n i n g an accurate and s t a b l e numer ica l scheme; the a d d i t i o n a l p h y s i c a l p r o p e r t i e s w i l l be i n c l u d e d l a t e r . Thi s model was l a t e r developed by Budd (1975) and Budd and Mclnnes (1974, 1978, 1979) to generate p e r i o d i c surges . Working from the assumption that ba sa l meltwater can cause s l i d i n g , they used the s t r a i n energy d i s s i p a t i o n to r e d i s t r i b u t e the basa l shear s t r e s s , caus ing l a rge l o n g i t u d i n a l s t r a i n r a te s and r a p i d f l o w . The s l i d i n g behaviour of the model, viewed as a q u a l i t a t i v e phenomenon, may be i t s most important c o n t r i b u t i o n to our ideas on surges . B indschad le r (unpublished) developed another f i n i t e d i f f e r e n c e p r o f i l e model s i m i l a r to the one I d e s c r i b e i n t h i s t h e s i s . B indschad le r a l s o d i d a c a r e f u l a n a l y s i s of numer ica l s t a b i l i t y and used a numer ica l scheme s i m i l a r to the one I d i s c u s s i n Appendix 1. He used t h i s model to i n v e s t i g a t e the changes i n the surge-type V a r i e g a t e d G l a c i e r , A l a s k a , d u r i n g i t s qu ie scent phase. Mahaffy (Mahaffy, u n p u b l i s h e d ; Mahaffy 1976; Mahaffy and Andrews 1976; Andrews and Mahaffy, 1976) used a two-dimensional f i n i t e d i f f e r e n c e model to study the i c e t h i c k n e s s and l a t e r a l extent of the Laurent ide i c e sheet and the Barnes i c e cap . 15 Jenssen (1977) p u b l i s h e d the only f u l l y three d imens iona l model of i c e sheets . T h i s model c a l c u l a t e d f l o w l i n e s and temperature as w e l l as i c e surface h e i g h t . I t s accuracy was l i m i t e d by severe computer s i z e l i m i t s (Jenssen used a mesh of 12 x 12 x 10 p o i n t s to represent the whole Greenland Ice Sheet ) . However, t h i s model i s an ambit ious development, and probably w i l l be fo l l owed by other models of t h i s k i n d . 1.2.3 PREVIOUS ICE TRAJECTORY MODELS Q u a n t i t a t i v e attempts to c a l c u l a t e s t r eaml ine s date back to the l a t e n i n e t e e n t h , c e n t u r y . Nansen ( c i t e d by Shumskiy, 1978, p . 133) c a l c u l a t e d flow l i n e s near the i ce d i v i d e i n c e n t r a l Greenland by assuming (1) steady s t a t e , (2) constant i c e t h i c k n e s s , d e n s i t y and mass ba lance , and (3) h o r i z o n t a l v e l o c i t y independent of depth . H a e f e l i (I963[b]) independent ly d e r i v e d the same s o l u t i o n . Re id (1894) and F i n s t e r w a l d e r (1897) independently developed a method of c a l c u l a t i n g s t r e a m l i n e s i n steady g l a c i e r s us ing the concept of f low tubes and p r o p e r t i e s of smooth vec tor f i e l d s . A l l these methods were on ly q u a l i t a t i v e ; none made any use of the c o n s t i t u t i v e p r o p e r t i e s of i c e . H a e f e l i (1961) d e r i v e d the v e l o c i t y f i e l d i n the c e n t r a l p o r t i o n of a steady i so therma l i c e sheet , assuming (1) no s l i d i n g , and (2) deformation by shear p a r a l l e l to the bed, us ing G l e n ' s f low law for i c e ( G l e n , 1955). Nagata (1977) developed an a n a l y t i c steady s t a t e i c e sheet model assuming no h o r i z o n t a l shear de format ion ; the f low was a l l 16 basa l s l i d i n g f o l l o w i n g a Weertman-type r e l a t i o n (see A 8 . 3 . 1 ) . By assuming that the v e r t i c a l v e l o c i t y was constant a long the i ce sheet s u r f a c e , Nagata a l s o d e r i v e d the s t r e a m l i n e s , and used the r e s u l t s to model the c o n c e n t r a t i o n of meteor i te s by g l a c i e r flow i n A n t a r c t i c a (Nagata, 1978). I use t h i s model as a t e s t for my numer ica l s t reaml ine c a l c u l a t i o n s i n Sec t ion 2 . 5 . 2 . N i e l s o n and Stockton ( 1956) derived^ the flow f i e l d i n v a l l e y g l a c i e r s of constant v a l l e y c r o s s - s e c t i o n assuming steady p l a s t i c f l o w , and Shumskiy (1967) found a s o l u t i o n for s t r e s s and v e l o c i t y i n a steady g l a c i e r w i t h n o n l i n e a r v i s c o s i t y . Severa l t r a j e c t o r y models have been d e r i v e d for regions near i c e d i v i d e s on steady i c e sheets i n order to date i c e cores (Dansgaard and Johnsen, I969[a ] ; P h i l b e r t h and Federer , 1971; Hammer and o t h e r s , 1978), and to model temperature w i t h depth (Weertman, 1968). A l l these models make some assumptions about v e r t i c a l s t r a i n ra te s or temperature g r a d i e n t s . Budd and others (1971) c a l c u l a t e d t r a j e c t o r i e s a long steady s t a t e A n t a r c t i c flow l i n e s assuming that the v e r t i c a l s t r a i n ra te was constant i n any v e r t i c a l column, (p. 51) or weighted by the h o r i z o n t a l v e l o c i t y v a r i a t i o n w i t h depth (p. 55) . T h i s process does not appear to s a t i s f y c o n t i n u i t y l o c a l l y . The i ce sur face e l e v a t i o n was c a l c u l a t e d by the numer ica l model d e s c r i b e d i n Sec t ion 1 .2 .2 . Jenssen (1977) c a l c u l a t e d t r a j e c t o r i e s i n a f i n i t e d i f f e r e n c e three d imens iona l time-dependent i c e sheet model, i n order to so lve for the temperature f i e l d . Jenssen (1978) a l s o found the t r a j e c t o r i e s of i c e p a r t i c l e s for a surg ing model of a f l o w l i n e through M i r n y , A n t a r c t i c a by 1 7 Budd and Mclnnes (1978). He d i d not d e s c r i b e the method used to o b t a i n the f l o w l i n e s . Th i s i s the only prev ious model, of which I am aware, to c a l c u l a t e the t r a j e c t o r i e s i n a t i m e - v a r y i n g i c e mass. 1.3 THE CONTINUITY EQUATION FOR AN ICE MASS 1.3.1 THE GENERAL GLACIER FLOW PROBLEM The equat ion of c o n t i n u i t y expresses the manner i n which the i c e mass changes i t s shape over t ime , i n response to mass input (accumulat ion or a b l a t i o n ) , sub ject to the phys i c s of deformat ion and s l i d i n g of i ce (Sec t ions 1.4 and 1.5) , and w i t h the assumptions about the flow f i e l d imposed below. I have d e r i v e d the c o n t i n u i t y equat ion for an i c e mass from f i r s t p r i n c i p l e s i n Appendix 5, u s ing s tandard methods of continuum mechanics (e.c[. T r u e s d e l l and Toupin , 1960; M a l v e r n , 1969; Prager , 1973). In t h i s S e c t i o n , I w i l l g ive the r e s u l t i n g e q u a t i o n , the assumptions i n v o l v e d i n i t s d e r i v a t i o n , and the p h y s i c a l i n t e r p r e t a t i o n of i t s terms. The c o o r d i n a t e system I have used i n t h i s study i s d e s c r i b e d i n S e c t i o n 1 .1 .2 . The p o s i t i o n vec tor x i s the t r i p l e t [ x , y , z ] , and the v e l o c i t y v e c t o r v(x) i s the t r i p l e t [ u ( x ) , w ( x ) , v ( x ) ] . Three ba s i c assumptions of the g l a c i e r model a r e : 1. matter i s conserved . 2. momentum i s conserved, i . e . a c c e l e r a t i o n n e g l i g i b l e . 3. i c e i s incompres s ib le (see Appendix 4 ) . 18 To so lve the genera l set of c o n s e r v a t i o n equat ions , c o n s t i t u t i v e equa t ions , and equat ions of s t a t e for the temperature d i s t r i b u t i o n , i n t e r n a l energy c o n t e n t , and a l l the components of the s t r e s s tensor and the v e l o c i t y f i e l d , w i t h boundary c o n d i t i o n s on a p o s s i b l y a r b i t r a r y boundary, i s a problem to i n s t i l l a sense of h u m i l i t y i n even the most e n t h u s i a s t i c and o p t i m i s t i c g l a c i o l o g i s t or- numer ica l a n a l y s t . A l l a t tempts , of which I am aware, to f i n d s o l u t i o n s to i ce flow problems s t a r t by making some a d d i t i o n a l assumptions about the channel geometry (boundaries and symmetry), or about the temperature d i s t r i b u t i o n (e.g_. i s o t h e r m a l ) , and/or about the flow f i e l d i t s e l f ( e . £ . plane s t r a i n , s imple shear , uniform s t r a i n r a t e s , e t c . ) . The model I de sc r ibe i n t h i s study i s no e x c e p t i o n . 1.3.2 RECTANGULAR CROSS-SECTION FLOW MODEL The g l a c i e r flow volume being modelled (see F i g u r e 1.1) i s assumed to have a r e c t a n g u l a r c r o s s - s e c t i o n , and a width W(x) . The two d imens iona l model i n c l u d e s v a r i a t i o n s i n the t h i r d dimension i n an approximate way by the assumptions that the v e l o c i t y components u and v are independent of y , and that the l a t e r a l component w v a r i e s l i n e a r l y w i t h y ( l a t e r a l s t r a i n ra te independent of y ) , such that the net v e l o c i t y y_(x) at the margins i s p a r a l l e l to the margins , i . . e . 19 dv = 0 dy (1 .3 .2 ) dv = u ( x , z , t ) dW dy W(x) dx (1 .3 .3 ) The g l a c i e r t h i c k n e s s h ( x , t ) i s a l s o independent of y . d h ( x , t ) = 0 dy (1 .3 .4 ) Nye ( 1 9 5 9 [ C ] , equat ion (33)) suggested t h i s method of i n c l u d i n g l a t e r a l v a r i a t i o n s , and Budd and Jenssen (1975, equat ion 3.34) i n c o r p o r a t e d i t i n t o t h e i r model . F igure 1.1 i l l u s t r a t e s the form of t h i s model. For an i ce sheet , W(x) can be the d i s t a n c e between two p o s s i b l y n o n p a r a l l e l f l o w l i n e s ; the ' w a l l s ' of the channel are a mathematical f i c t i o n , and the assumptions (1 .3 .1 ) through (1 .3 .4 ) are reasonable . For a v a l l e y g l a c i e r , drag from the v a l l e y w a l l s i s impor tant , and the i c e t h i c k n e s s and g l a c i e r bed vary w i t h y . I f I attempt to i d e n t i f y W(x) w i t h the v a l l e y w i d t h , the assumptions (1 .3 .1 ) to (1 .3 .4 ) may be g r o s s l y v i o l a t e d . I f , however, I l e t W(x) be the d i s t a n c e between two f l o w l i n e s near the g l a c i e r c e n t r e l i n e , §_.£. W(x) may be a few percent of the v a l l e y w i d t h at the l e v e l of the i c e s u r f a c e , then a l l the assumptions are reasonable , and I o b t a i n a c e n t r a l f l o w l i n e s o l u t i o n , but w i t h the l a t e r a l v a r i a t i o n i n v a l l e y width i n c l u d e d to a good a p p r o x i m a t i o n . The e f f e c t of the v a l l e y s i d e w a l l drag can be i n c l u d e d i n an approximate way by us ing shape f a c t o r s (Nye, 1 9 6 5 [ C ] ) to modify the shear s t r e s s (see equat ion (1 .4 .25) be low) . 20 FIGURE 1 . 1 . Rectangular C r o s s - s e c t i o n F low. The t r i a d x -y-z shows the coord ina te axes , and the bo ld arrows u , v , and w are the vec tor components of the v e l o c i t y v . The example shows the v e l o c i t y i n the accumulat ion zone (v i s n e g a t i v e ) . 1 . 3 . 3 THE CONTINUITY EQUATION With these assumptions, the well-known c o n t i n u i t y equat ion (the d e r i v a t i o n of which I show i n Appendix 5 ) i s d h ( x , t ) + 1 g_Q(x,t) = A ( x , t ) at wTx) dx. ( 1 . 3 . 5 ) where h ( x , t ) i s the i c e t h i c k n e s s normal to the bed, Q(x , t ) i s the i c e f l u x through a c r o s s - s e c t i o n from bed to s u r f a c e , and the source term A ( x , t ) i s the net mass balance normal to the bed, the net accumulat ion or a b l a t i o n ra te i n i ce e q u i v a l e n t t h i c k n e s s u n i t s per u n i t t i m e , i n c l u d i n g s n o w f a l l , sur face m e l t i n g , and basa l m e l t i n g and r e f r e e z i n g . The l a s t two 21 ' c o n t r i b u t i o n s are u s u a l l y n e g l i g i b l y sma l l (e.g.. R o t h l i s b e r g e r , 1972). The i c e f l u x Q(x , t ) i s d e f i n e d by /%h(x,t) Q ( x , t ) = W(x) I u ( x , z , t ) dz (1 .3 .6) J 0 where u ( x , z , t ) i s the v e l o c i t y component p a r a l l e l to the bed. I t w i l l be d e r i v e d i n the next s e c t i o n . Q ( x , t ) can a l s o be w r i t t e n Q(x , t ) = V ( x , t ) h ( x , t ) W(x) (1 .3 .7) where V ( x , t ) i s u ( x , t ) , i . . e . the downslope v e l o c i t y u ( x , z , t ) averaged between the bed and the s u r f a c e . J h ( x , t ) u ( x , z , t ) dz (1 .3 .8) . . . . . . 0 Assuming that the upstream end of the g l a c i e r s e c t i o n under c o n s i d e r a t i o n i s at x=0, the boundary c o n d i t i o n i s Q(0 , t ) = Q ( t ) 0 (1 .3 .9) I f x=0 a c t u a l l y represents the p h y s i c a l upper extent of the i c e mass, Q 0 ( t ) i s i d e n t i c a l l y z e r o . For the case of an i c e d i v i d e , t h i s i s achieved by s e t t i n g V ( 0 , t ) = 0 (1 .3 .10) by a v a n i s h i n g i c e surface s lope angle (see ( 1 .4 . 3 8 ) ) , and l e t t i n g i c e t h i c k n e s s h ( 0 , t ) v a r y . For the case of a v a l l e y g l a c i e r o r i g i n a t i n g on a s l o p e , the boundary c o n d i t i o n (1 .3 .9) i s ach ieved by s e t t i n g h ( 0 , t ) = 0 (1 .3 .11) I f the lower end x=L(t) i s the g l a c i e r t e rminus , then L ( t ) i s d e f i n e d i m p l i c i t l y by 22 h ( L ( t ) , t ) = 0 (1 .3 .12) (This i s not a mathematical boundary c o n d i t i o n , but a p h y s i c a l c o n d i t i o n d e f i n i n g the l i m i t s of the i c e mass.) An i n i t i a l c o n d i t i o n of the form h (x ,0 ) = H (x) 0 (1 .3 .13) i s a l s o r e q u i r e d . For example, one s imple i n i t i a l c o n d i t i o n i s H o (x )=0 , i_.e. u n g l a c i e r i z e d ground. 1.3.4 PHYSICAL INTERPRETATION Equat ion (1 .3 .5 ) may be i n t e r p r e t e d i n the f o l l o w i n g manner. Consider a v e r t i c a l pr i sm of i c e as shown i n F igure 1.2, extending from the bed to the surface h ( x , t ) , w i t h wid th W(x) i n the y d i r e c t i o n , and t h i c k n e s s 6x i n the x d i r e c t i o n . Let p be the constant d e n s i t y of g l a c i e r i c e . When (1 .3 .5 ) i s m u l t i p l i e d by the constant />W(x)6x6t, the f i r s t term i s the net change i n mass i n the pr i sm i n a time 6t (the pr i sm i s then t a l l e r or s h o r t e r ) . The second term on the l e f t i s the d i f f e r e n c e i n mass between that which flowed out of the pr i sm through the downslope f ace , and that which flowed i n t o the pr i sm through the upstream face , i n the time 6 t . T h i s i s the net l o s s of mass from the pr i sm i n t o the downstream f l o w . The term on the r i g h t s ide i s j u s t the t o t a l mass added to the pr i sm i n time 6t by s n o w f a l l or m e l t i n g . Thus, (1 .3 .5 ) s t a t e s that the t o t a l l a y e r of mass added t o the top of the pr i sm at any p o s i t i o n x i s the sum of the s n o w f a l l onto the pr i sm and the net mass l e f t i n s i d e the pr i sm by s p a t i a l f low v a r i a t i o n s . 23 1.4 PHYSICS OF ICE DEFORMATION 1.4.1 INTRODUCTION In t h i s S e c t i o n , I w i l l o u t l i n e the d e r i v a t i o n of a second equat ion r e l a t i n g g l a c i e r t h i c k n e s s and i ce f l u x so that the c o n t i n u i t y equat ion (1 .3 .5 ) can be s o l v e d . There are three steps i n t h i s d e r i v a t i o n . Th i s i s a s tandard procedure i n m o d e l l i n g g l a c i e r flow (e.g.. Pa te r son , 1980; Raymond, 1980). Newton's second law e s t a b l i s h e s r e l a t i o n s h i p s between the sur face and body forces and the a c c e l e r a t i o n s i n any continuum. The s t r e s s e q u i l i b r i u m equat ions are o u t l i n e d i n S e c t i o n 1 .4 .2 . Second, observa t ions and theory of the deformat ion of g l a c i e r i c e e s t a b l i s h c o n s t i t u t i v e r e l a t i o n s h i p s between the s t re s se s a p p l i e d to i c e , and the r e s u l t i n g de format ion . G l e n ' s flow law 24 for i c e ( G l e n , 1955, 1958) i s presented i n S e c t i o n 1 .4 .3 . F i n a l l y , i n Sec t ion 1.4.4, a f t e r making s i m p l i f y i n g assumptions, s u b s t i t u t i n g the s t r e s s equat ions of Sec t ion 1.4.2 i n t o G l e n ' s flow law to get the s t r a i n r a t e , and i n t e g r a t i n g over a p p r o p r i a t e c o o r d i n a t e s , I express the downslope component of the i ce v e l o c i t y i n terms of the i c e t h i c k n e s s . S ince i ce f l u x , t h i c k n e s s , and v e l o c i t y are r e l a t e d through ( 1 . 3 . 7 ) , t h i s w i l l complete the d e r i v a t i o n of a second equat ion needed to so lve the c o n t i n u i t y equat ion (1..3.5) for i c e t h i c k n e s s and f l u x . 1.4.2 STRESS EQUILIBRIUM EQUATIONS Since g l a c i e r i c e deforms s l o w l y , the a c c e l e r a t i o n term i n Newton's second law can be n e g l e c t e d , l e a v i n g the r e s u l t t h a t , for any volume V of a s l o w l y deforming cont inuous medium, B + T = 0 (1 .4 .1) where*B i s the t o t a l body fo rce found by i n t e g r a t i n g the s p e c i f i c body forces ( force per u n i t mass) f.(r) at p o s i t i o n £ over a l l r_ throughout the volume V . I t s components are B V (1 .4 .2 ) and T i s the t o t a l sur face t r a c t i o n on the sur face S e n c l o s i n g the volume V w i t h surface normal v e c t o r n . The components of T are T = | | n (r ) c ( r ) dS jk (1 .4 .3) S tf.k i s the s t r e s s t e n s o r , i_.e. the force per u n i t area on a surface normal to the X j a x i s a c t i n g i n the x k d i r e c t i o n . The 25 E i n s t e i n c o n v e n t i o n , whereby repeated i n d i c e s are summed, i s used i n t h i s s e c t i o n . Orthogonal axes are assumed. A p p l y i n g Gauss' Theorem (e.g.. Prager , 1973, p . 29) to (1 .4 .3 ) and s u b s t i t u t i n g (1 .4 .2 ) and (1 .4 .3 ) i n t o (1 .4 .1 ) by components g ives d 3 r = 0 dx V 1 (1 .4 .4 ) S ince the volume V i s a r b i t r a r y , the g l o b a l equat ion (1 .4 .4 ) has a l o c a l counterpar t Ic ( r ) (r) f ( r ) + k i = 0 dx (1 .4 .5) k Assuming that the net angular a c c e l e r a t i o n i s z e r o , and s e t t i n g to zero the sum of moments a c t i n g on the volume V , g ives the r e s u l t c = <t (1 .4 .6 ) i j J i . i . e . the s t r e s s tensor i s symmetric . The development i s very s i m i l a r to (1 .4 .1 ) through ( 1 . 4 . 5 ) , and i s g iven i n Prager (1973, p . 47 ) . For a g l a c i e r , the only body force i s g r a v i t y , so f ( r ) = g (1 .4 .7 ) When the x , a x i s ( to be c a l l e d x) i s taken a long the g l a c i e r bed which i s at an angle fix) to the h o r i z o n t a l , the x 3 a x i s (z) i s normal to i t and upward, and the x 2 a x i s (y) i s h o r i z o n t a l , then the s t r e s s equat ions (1 .4 .5 ) are 26 off <̂ff he x x + x y + x z + pgsin(p)= 0 (1 .4 .8 ) dx dy Tz ^ff <5ff dff x z + y z + z z - />gcos(*) = 0 (1 .4 .9) <̂x by £z do 5ff f̂f xy_ + Yl + v z = 0 (1 .4 .10) dx dy T z 1.4.3 CONSTITUTIVE RELATION FOR ICE DEFORMATION The c o n s t i t u t i v e equat ions for deformation r e l a t e the s t r e s se s a p p l i e d to i c e to the r e s u l t i n g deformation r a t e . The components of the s t r a i n ra te tensor are 1 D — i + —j dx dx • j i (1 .4 .11) where u^ i s the i t h component of the i c e v e l o c i t y . Glen (1958) showed that the most genera l r e l a t i o n between the s t r e s s tensor and the s t r a i n ra te tensor for a n o n l i n e a r , o r i g i n a l l y i s o t r o p i c m a t e r i a l had the form € = A(T ,T ,T )6 + B(T ,T ,T )c i j 1 2 3 i j 1 2 3 i j + C(T ,T ,T )<s a 1 2 3 ik k j (1 .4 .12) where 6-tj i s the Kroenecker d e l t a 27 6 = 0 i ^ j i j = 1 i=j (1 .4 .13) and T 1 f T 2 , and T 3 are the f i r s t , second, and t h i r d s c a l a r i n v a r i a n t s of the s t r e s s tensor (e.g.. Prager , 1973, p . 22) T = c 1 i i (1 .4 .14) T = J_( <r e ~ c a ) 2 2 i j j i i i j j (1 .4 .15) T - M2<y a c - 3* c e + t f * t f ) 3 6 i j jk k i i j j i kk i i j j kk (1 .4 .16) Terms i n h igher powers of • can be e l i m i n a t e d by the Hami l ton-Cay ley equat ion (§_.g_. Prager , 1973, p . 25 ) . The stress-dependent c o e f f i c i e n t s must be f u n c t i o n s o f , at most, the s c a l a r i n v a r i a n t s , because the r e l a t i o n (1 .4 .12) i s independent of the cho ice of axes . The c o e f f i c i e n t s A , B, and C may be dependent on temperature . Rigsby (1958) showed t h a t , to a good approx imat ion , the deformation ra te of i c e c r y s t a l s was independent of the h y d r o s t a t i c pressure p , where P = 1 « 3 i i (1 .4 .17) when the i c e temperature was measured r e l a t i v e to the p r e s s u r e - m e l t i n g p o i n t . T h i s r e s u l t means that the c o n s t i t u t i v e r e l a t i o n can be w r i t t e n more s imply i n terms of the d e v i a t o r i c s t r e s s * | j 28 cr ' = e - _1_ 6 cr i j i j 3 i j kk (1 .4 .18) By d e f i n i t i o n , the f i r s t s c a l a r i n v a r i a n t T', of cr'j i s z e r o . By f u r t h e r assumptions, f i r s t tha t the d e n s i t y of g l a c i e r i ce i s constant (see Appendix 4 ) , second, t h a t , for a g iven s t r e s s , the components of s t r a i n ra te are p r o p o r t i o n a l t o the components of the s t r e s s d e v i a t o r t e n s p r , and t h i r d , that the second i n v a r i a n t of the s t r a i n r a te tensor i s a f u n c t i o n of T 2 only (see G l e n , 1958) the f low law presented by Nye (1953) reduced the genera l r e l a t i o n (1 .4 .12) to € = B(T' ) <r' i j 2 i j (1 .4 .19) When € and T are the square root s of the second s c a l a r i n v a r i a n t s of and (T = JT2 ) , Nye (1953) showed that a p l a u s i b l e r e l a t i o n s h i p was n t = A T (1 .4 .20) Using the case of a s imple s l a b deforming by shear p a r a l l e l to the x a x i s , combining (1 .4 .20) w i t h (1 .4 .19) i m p l i e s that n-1 B (T ' ) = A T (1 .4 .21) 2 and n-1 e = A T * ' (1 .4 .22) i j i j The exponent n i n (1 .4 .20) i s independent of temperature. Values i n the l i t e r a t u r e vary from 1.5 (Gerrard and o t h e r s , 1952) to 4.2 (G len , 1955), and the va lue u s u a l l y used for g l a c i e r m o d e l l i n g i s n=3 (e.cj . P a t e r s o n , 1980). The f a c t o r A 29 f o l l o w s the e x p o n e n t i a l A r r h e n i u s temperature dependence A = A exp(-Q/RT) 0 (1 .4 .23) where A 0 i s a c o n s t a n t , R i s the gas constant (8.314 J K ^ m o l - 1 ) , Q i s the a c t i v a t i o n energy fo r c reep , and T i s the temperature ( ° K ) . Values of Q for secondary creep of p o l y c r y s t a l l i n e i c e are 60 k j m o l " 1 f o r T < - 1 0 ° C , and approximate ly 139 kJ m o l " 1 for T > -10°C (Pa ter son , 1981, p . 34 ) . The presence of s m a l l amounts of water causes g r a i n boundary s l i d i n g (Barnes and o t h e r s , 1971? Jones and Brunet , 1978) above - 1 0 ° C . The deformat ion below -8 °C i s dominated by basa l g l i d e (Barnes and o t h e r s , 1971). In t h i s s tudy, the i c e i s assumed to be i s o t h e r m a l at 0 ° C , and the flow law parameters used are (Pa ter son , 1981, Table 3 . 3 , P- 39) n=3 A = 5.3 1 0 " 1 5 s" 1 k P a " 3 (1 .4 .24) These va lues apply only for secondary c reep , a f t e r the i n i t i a l t r a n s i e n t response to l o a d i n g has d i e d away. Other c o n s t i t u t i v e r e l a t i o n s have been proposed for g l a c i e r i c e , such as a h y p e r b o l i c s ine r e l a t i o n (Barnes and o t h e r s , 1971), or a po lynomia l w i t h odd order s t r e s s terms (Meier , 1958, 1960; L l i b o u t r y , I969[a ] ; Colbeck and Evans, 1973). However, l a b o r a t o r y experiments (e.g.. G l e n , 1952, 1955; Steinemann, 1958) and f i e l d measurements of c l o s u r e of boreholes and tunnel s and deformat ion of boreholes (e.cj . Gerra rd and o t h e r s , 1952; Nye, 1953; Mathews, 1959; M e i e r , 1960; Paterson and Savage, I963[a] , 1963[b]; H a e f e l i , 1963[a] ; Shreve and Sharp, 1970; Raymond, 1971; P a t e r s o n , 1977) and the f low of i c e she lves (Thomas, 1973) 30 i n d i c a t e that ( 1 . 4 . 2 2 ) , known as G l e n ' s flow law, i s a u s e f u l and s a t i s f a c t o r y c o n s t i t u t i v e r e l a t i o n for g l a c i e r i c e . 1.4.4 SHEAR STRESS AND ICE FLUX In t h i s s e c t i o n , I w i l l o u t l i n e the d e r i v a t i o n of the shear s t r e s s p a r a l l e l to the bed, and how i t can be i n t e g r a t e d to g ive the downslope v e l o c i t y component and the i ce f l u x . The e r r o r s and assumptions are e x p l i c i t l y shown. The d e t a i l s are g iven i n Appendix 7. The s t r e s s e q u i l i b r i u m equat ion (1 .4 .8 ) can be i n t e g r a t e d from the sur face to a he ight z above the g l a c i e r bed to g ive the shear s t r e s s c ( x , z ) = s / > g ( h - z ) s i n « xz 1 + 0 9*' 3*' 2h xx + h yy dx dx pqha (1 .4 .25) which i s d e r i v e d from (A7.3 .21) i n Appendix 7. The l e a d i n g f a c t o r i s the standard formula for shear s t r e s s i n a p a r a l l e l - s i d e d s l a b deforming by s imple shear p a r a l l e l to the bed (e.c[. P a t e r s o n , 1969, p . 91) when the shape f a c t o r s i s u n i t y (see Appendix 7, S e c t i o n A 7 . 4 ) . The surface s lope a i s an e f f e c t i v e s lope averaged over a d i s t a n c e of at l e a s t the order of 4h (Budd , l968 ; I970 [a ] ) . The average s t r e s s d e v i a t o r s i n the second term are de f ined by 31 h(x) be' = _ 1 _ r d<r' z<h xx ( h - z ) l xx dz' dx J dx (1 .4 .26) z = 0 z=h with a corresponding d e f i n i t i o n f o r the yy component. If I' co n s i d e r the x- and y - d i r e c t e d f o r c e s on an i c e column from the i c e s u r f a c e to the bed, the c o r r e c t i o n terms i n (1 .4 .25) are, very roughly speaking, r a t i o s of the normal f o r c e s to the b a s a l shear f o r c e . These r a t i o s are u s u a l l y very small f o r g l a c i e r s and i c e sheets. L l i b o u t r y ( I 9 5 8[b]), Shumskiy (1961), and Robin (1967) used a c o r r e c t i o n term s i m i l a r to t h i s to account f o r l o n g i t u d i n a l s t r a i n , and C o l l i n s (1968) p u b l i s h e d a mathematical j u s t i f i c a t i o n of i t . Nye ( I 9 6 9[a]) s i m p l i f i e d the a n a l y t i c a l f o r m u l a t i o n by an a p p r o p r i a t e c h o i c e of axes. My f o r m u l a t i o n d i f f e r s i n some d e t a i l s , p a r t l y because I use the axes of the numerical model (see S e c t i o n 1 .3) . Budd (1968; 1970[a]; I970[b]; 1971) gave a d e t a i l e d d i s c u s s i o n of s t r e s s v a r i a t i o n s i n g l a c i e r s , i n c l u d i n g c o r r e c t i o n terms and the wavelength ranges f o r which they may be important. Hutter ( i n press) g i v e s the most recent and r i g o r o u s treatment of s t r e s s i n g l a c i e r s . The shape f a c t o r i s an approximate c o r r e c t i o n between zero and u n i t y f o r the r e d u c t i o n i n the shear s t r e s s « r x z along the channel c e n t r e l i n e when some of the weight of the g l a c i e r i s supported by the v a l l e y s i d e w a l l s . I t was f i r s t used q u a n t i t a t i v e l y i n work on r e c t i l i n e a r flow i n r e c t a n g u l a r , p a r a b o l i c , and e l l i p t i c a l channels by Nye ( I 9 6 5 [ c ] ) . I d e s c r i b e shape f a c t o r s i n more d e t a i l i n Appendix 7, S e c t i o n A7.4. 32 In a r r i v i n g at (1 .4 .24) i n Appendix 7, I assumed that the s lope angles o of the i c e s u r f a c e , and p of the g l a c i e r bed were s m a l l , i . e . | o ( x ) | « 1 ( 1 .4 .27) |f(x) | « 1 (1 .4 .28) I a l s o assumed that a was never n e g l i g i b l y smal l compared to p. T h i s may not be true near an i c e d i v i d e . Al though (1 .4 .8) c o n t a i n s s i n $ , the f i n a l r e s u l t (1 .4 .25) for the shear s t r e s s depends on ly on s i n e . Nye ( l952[b] ) f i r s t p o i n t e d out t h i s r e s u l t . When the smal l angle assumptions (1 .4 .27) and (1 .4 .28) h o l d , the l o n g i t u d i n a l s t r e s s grad ient term d * ^ /dx i n (1 .4 .8 ) i n t r o d u c e s a term i n io-p) which cance l s the bed s lope dependence, l e a v i n g only the sur face slope dependence of ( 1 . 4 . 2 5 ) . With the assumption that the major shear deformat ion occurs p a r a l l e l to the bed, . i . e . iv /bu "57/ bz « 1 (1 .4 .29) the shear s t r a i n ra te * = 1 xz 2 du + ^v dz dx (1 .4 .30) i s approx imate ly € = 1 l u xz 2 bz (1 .4 .31) which can be i n t e g r a t e d d i r e c t l y , from the bed to he ight z , to g ive 33 u ( x , z ) = u (x) + z e dz ' xz 0 (1 .4 .32) where u_(x) i s the basa l s l i d i n g v e l o c i t y d i s cus sed i n Sec t ion 1.5. I f I f u r t h e r assume that T, the square root of the second i n v a r i a n t of the s t r e s s d e v i a t o r tensor (1 .4 .19) i s approx imate ly equal to the shear s t r e s s exz , i_.e. shear s t r e s s p a r a l l e l to the bed i s by far the l a r g e s t s t r e s s d e v i a t o r component, or 6T T - e e = xz xz e xz << 1 (1 .4 .33) then I can s u b s t i t u t e G l e n ' s f low law (1 .4 .22) for the s t r a i n ra te « x Z i n ( 1 . 4 . 3 2 ) , u s ing <rxz from (1 .4 .25) for both T and * , to get the v e l o c i t y component u ( x , z ) p a r a l l e l to the bed. u ( x , z ) - u (x) = n n+1 n+1 2 A [ s ( x ) P q s i n ( g ( x ) ) ] [h - (h-z) ] [1 + e ( x ) ] (n+T5 (1 .4 .34) where the e r r o r e ( x ) , i . e . the terms not i n c l u d e d i n the computer model , has the form e(x) = 0 he' de' 2h xx + h yy dx dx + (n-1) 6T e - b_v /fru X Z bx/ bz pgho (1 .4 .35) where the symbol 0 [x ] means " i s of the order of x", i . e . the f u n c t i o n goes to zero at the same ra te as x . I have assumed that 34 the g l a c i e r i s i s o t h e r m a l , i . e . temperate, so that the c o e f f i c i e n t A of G l e n ' s flow l aw(1 .4 .22 ) can be t r e a t e d as a c o n s t a n t . I f the i c e temperature v a r i e s w i t h z , the i n t e g r a l can be eva lua ted n u m e r i c a l l y . The downslope i c e f l u x for use i n the c o n t i n u i t y equat ion (1 .3 .5 ) i s h ( x , t ) Q(x , t ) =1 u ( x , z , t ) dz ^0 n n+2 = u ( x , t ) h ( x , t ) + 2A[s/pgsinq] [ h ( x , t ) ] [1+ e (x ) ] s (n+2) (1 .4 .36) With the assumptions d i scus sed above, the e r r o r term i n v o l v i n g e(x) i s s m a l l ; i t i s neg lec ted i n the computer model i n i t s present form. The average v e l o c i t y V ( x , t ) used i n S e c t i o n 1.3 and Appendix 1 i s de f ined as V ( x , t ) = Q ( x , t ) / h ( x , t ) (1 .4 .37) which i s n n+1 V ( x , t ) = u ( x , t ) + 2A[spgsino] [ h ( x , t ) ] [1+ e (x ) ] s (n+2) (1 .4 .38) The term on the r i g h t due to the i n t e r n a l deformation i s j u s t (n+1)/(n+2) times the downslope v e l o c i t y component at the i c e sur face u ( x , h ( x ) , t ) from ( 1 . 4 . 3 4 ) . 35 1.5 BASAL SLIDING 1.5.1 INTRODUCTION The one q u a n t i t y s t i l l r e q u i r e d to complete the d e r i v a t i o n of the downslope v e l o c i t y u ( x , z ) and the f l u x Q(x) i s the basa l s l i d i n g v e l o c i t y u 5 ( x ) which appeared as an i n t e g r a t i o n constant in ( 1 . 4 . 3 2 ) . Raymond (1980), i n a recent rev iew, gave a summary of s l i d i n g behaviour , measurements, and the p h y s i c a l processes p o s s i b l y i n v o l v e d , and p o i n t e d out some d i f f i c u l t i e s of q u a n t i t a t i v e m o d e l l i n g of s l i d i n g . In Appendix 8, I summarize c u r r e n t ideas on the phys i c s of s l i d i n g , and review the use of s l i d i n g i n computer models. In t h i s s e c t i o n , I d i s c u s s the importance of s l i d i n g , the way I t r e a t s l i d i n g i n Chapter 3, and the reason for my c h o i c e . 1.5.2 BASAL ICE TEMPERATURE AND SLIDING Ice masses which are c o l d at the base, have temperatures below the pressure m e l t i n g p o i n t , do not appear to s l i d e . The basa l i c e i s e f f e c t i v e l y f rozen to the g l a c i e r bed, and u (x) = 0 s (1 .5 .1 ) Ice masses which have temperatures at the pressure m e l t i n g p o i n t at the i c e - r o c k i n t e r f a c e e x h i b i t s l i d i n g v e l o c i t i e s which range from zero to va lues much g rea te r than the v e l o c i t i e s due to i n t e r n a l de format ion . The model I d e s c r i b e i n t h i s study assumes an i so therma l i c e mass. Due to the e x i s t e n c e of the 36 geothermal heat f l u x , and the r e s u l t i n g geothermal temperature g r a d i e n t , the only p o s s i b l e e s s e n t i a l l y i so therma l i c e mass i s one at the pressure m e l t i n g p o i n t throughout i t s volume ( n e g l e c t i n g a p o s s i b l y c o l d sur face l a y e r caused by d i f f u s i v e p e n e t r a t i o n of the winter c o l d wave), because only then can the geothermal f l u x be absorbed at the base by being transformed i n t o energy of f u s i o n . T h i s means that s l i d i n g v e l o c i t i e s can be an important component of motion for my mode l l ing s i t u a t i o n s . 1.5.3 PHILOSOPHY OF SLIDING IN THIS STUDY C o r r e c t l y m o d e l l i n g the p h y s i c a l processes of g l a c i e r s l i d i n g i s , at the p re sen t , very d i f f i c u l t , due to inadequate o b s e r v a t i o n s , and the l a rge number of u n c o n t r o l l e d p h y s i c a l v a r i a b l e s p o s s i b l y i n v o l v e d i n s l i d i n g proces ses . In Appendix 8, I have o u t l i n e d the problems of measurements, the p h y s i c a l c o m p l i c a t i o n s of s l i d i n g proces ses , the present s t a t e of s l i d i n g theory and i t s q u a n t i t a t i v e a p p l i c a t i o n i n computer models. In the models presented i n t h i s s tudy , I do not attempt to i n v e s t i g a t e or to s imula te the p h y s i c s of g l a c i e r s l i d i n g . My a im, i n Chapter 3, i s to i n v e s t i g a t e the consequences of surg ing (de f ined by a p e r i o d i c s l i d i n g h i s t o r y ) on s t r u c t u r e s w i t h i n a g l a c i e r , g iven that p e r i o d i c surg ing occurs i n the d e f i n e d manner. I do not attempt to induce surges i n the model by any p a r t i c u l a r p h y s i c a l mechanism. T h i s approach to i n v e s t i g a t i n g e f f e c t s of surg ing was a l s o used by Campbell and Rasmussen (1969) and by C l a r k e (1976). I c o u l d e a s i l y i n c o r p o r a t e the t h e o r i e s of Weertman (1957), 37 Nye (I969[b]y 1970), Kamb (1970), Morland (1976[a l ; I976[b] ) , or Budd (1975) to c a l c u l a t e the s l i d i n g v e l o c i t i e s , but the r e s u l t s would be n u m e r i c a l l y suspect , due to the problems d i scus sed i n Appendix 8, and would add noth ing to my i n v e s t i g a t i o n of the consequences of s u r g i n g . My approach i s , i n s t e a d , to use a predetermined s l i d i n g f u n c t i o n u 5 ( x , t ) (based as c l o s e l y as p o s s i b l e on the reasonably w e l l i n f e r r e d s l i d i n g h i s t o r y of a surg ing g l a c i e r such as the S tee l e ) as a d r i v i n g f u n c t i o n f o r p e r i o d i c surges i n the computer model . I c a l c u l a t e the response of the g l a c i e r model to t h i s d r i v i n g f u n c t i o n by us ing c o n t i n u i t y and G l e n ' s flow law to f i n d the i n t e r n a l de format ion . For my purpose of f i n d i n g the e f f e c t s of surg ing on the i n t e r n a l s t r u c t u r e , t h i s approach i s no worse than us ing a n u m e r i c a l l y inadequate s l i d i n g theory , and i t has the d i s t i n c t advantage that I can c o n t r o l the s l i d i n g at w i l l w h i l e I r e l a t e s l i d i n g pa t t e rns to r e s u l t i n g changes w i t h i n the i c e mass. 38 CHAPTER 2: MODELS AND TESTS "Things are seldom what they seem; Skim m i l k masquerades as c r e a m . " 1 2.1 INTRODUCTION 2.1.1 OUTLINE In t h i s chapter , I o u t l i n e the opera t ion of the computer models and I de sc r ibe t e s t s used to v e r i f y t h e i r c o r r e c t o p e r a t i o n . In t h i s i n t r o d u c t o r y s e c t i o n , I e x p l a i n why I t h i n k t e s t s are impor tant . In Sec t ion 2 . 2 , I d e s c r i b e the c o n t i n u i t y equat ion g l a c i e r p r o f i l e model , and i n Sec t ion 2 . 3 , how I have t e s t e d i t . In s e c t i o n 2.4 I d e s c r i b e the p a r t i c l e t r a j e c t o r y c a l c u l a t i o n s , and i n S e c t i o n 2 . 5 , how they were t e s t e d . 2 .1 .2 IMPORTANCE OF MODEL TESTING A n a l y t i c a l s o l u t i o n s of i n i t i a l va lue problems are most d e s i r a b l e , because the c o r r e c t n e s s of the s o l u t i o n can be v e r i f i e d fo r a l l space and time s imply by s u b s t i t u t i n g the s o l u t i o n i n t o the d i f f e r e n t i a l e q u a t i o n . U n f o r t u n a t e l y , a n a l y t i c a l s o l u t i o n s to i c e flow problems are r e s t r i c t e d to a few cases w i t h s imple boundary c o n d i t i o n s , uncomplicated r h e o l o g i e s , and, o f t e n , steady s t a t e s . A f i n i t e d i f f e r e n c e numer ica l model uses a set of a l g e b r a i c 1 H . M . S . P i n a f o r e . G i l b e r t and S u l l i v a n . 39 equat ions whose s o l u t i o n c l o s e l y approximates a d i g i t i z e d v e r s i o n of the t rue s o l u t i o n of the d i f f e r e n t i a l e q u a t i o n , to w i t h i n a t r u n c a t i o n e r r o r (see Appendix 1, S e c t i o n A 1 . 5 . 3 ) . Numerical s o l u t i o n s can extend the domain of s o l v a b l e problems to i n c l u d e those w i t h q u i t e genera l boundary c o n d i t i o n s , i c e rheo logy , and temporal v a r i a t i o n s . The p r i c e which i s pa id for t h i s increa sed g e n e r a l i t y , however, i s a new inherent u n c e r t a i n t y i n the v a l i d i t y of the s o l u t i o n o b t a i n e d . S u b s t i t u t i o n of a numer ica l s o l u t i o n i n t o the f i n i t e d i f f e r e n c e equa t ions , or i n t o the d i f f e r e n t i a l e q u a t i o n , can g i v e , at most, the r e s i d u a l e r r o r s i n the most recent time s t ep . The long term i n t e g r a t e d e r r o r i s unknown. The danger i s that a numer ica l s o l u t i o n w i l l behave i n a q u a l i t a t i v e l y reasonable manner, yet q u a n t i t a t i v e l y may be, over some time s c a l e s , g r o s s l y wrong. For i n s t a n c e , i c e v e l o c i t i e s and t h i c k n e s s e s may be i n e r r o r by tens of pe rcent , and the phase of c y c l i c phenomena such as surg ing may become t o t a l l y u n r e l a t e d to the phase of the t rue s o l u t i o n . These p o s s i b i l i t i e s should make us q u i t e c a u t i o u s about any p r e d i c t i v e c l a i m s made for numer ica l models. We c o u l d ob ta in a r e s u l t no more a c c u r a t e , at be s t , than an educated guess, yet be l u l l e d i n t o b e l i e v i n g that i t was an q u a n t i t a t i v e p r e d i c t i o n of g l a c i e r behav iour . There are two p o s s i b l e sources of e r r o r i n numer ica l s o l u t i o n s . F i r s t , the computer program may not c o r r e c t l y so lve the set of a l g e b r a i c e q u a t i o n s . Programming e r r o r s , such as i n c o r r e c t cons tants or m i s s i n g minus s i g n s , sometimes go undetected . Spurious numer ica l " s o l u t i o n s " of t h i s k i n d have on occas ion found t h e i r way i n t o the s c i e n t i f i c l i t e r a t u r e . By 40 c a r e f u l program design and the use of c o n s i s t e n c y checks , these e r r o r s can be e l i m i n a t e d , a l though not a l l programmers have taken the time to do so . Assuming then , tha t the computer program works c o r r e c t l y , there i s s t i l l another source of e r r o r . The numer ica l scheme used i n the computer model may not adequately represent the d i f f e r e n t i a l equat ion at a l l time s c a l e s . -One example of t h i s i s the l i n e a r computat iona l i n s t a b i l i t y (Appendix 1, S e c t i o n A 1 . 4 . 2 ) . The s o l u t i o n of the f i n i t e d i f f e r e n c e equat ions may i n c l u d e an e x p o n e n t i a l l y growing h igh wavenumber o s c i l l a t i o n complete ly u n r e l a t e d to the d i f f e r e n t i a l e q u a t i o n . F o r t u n a t e l y , t h i s e r r o r i s u s u a l l y easy to r e c o g n i z e ! The other cause fo r concern i s the p o s s i b i l i t y that the numer ica l s o l u t i o n may d r i f t away from the t rue s o l u t i o n , yet s t i l l look " p h y s i c a l l y r e a s o n a b l e " . T h i s c o u l d be caused by any of s e v e r a l f a c t o r s . For example, i n a p p r o p r i a t e mesh i n t e r v a l s (Appendix 1, Sec t ion A 1 . 5 . 2 ) may cause i n c o r r e c t d i s p e r s i o n and s p e c t r a l a t t e n u a t i o n , d i s t o r t i n g the s o l u t i o n to an unacceptable degree. I n t r o d u c i n g smoothing schemes in attempts to suppress numer ica l i n s t a b i l i t i e s can cause s i m i l a r d i s t o r t i o n s (and may s t i l l f a i l to t o t a l l y remove the i n s t a b i l i t i e s ) . For example, the w i d e l y repor ted i c e sheet model of the Melbourne group (Budd and Jenssen , 1975; Budd and Radok, 1971) was used by Mclnnes ( u n p u b l i s h e d ) , who appears to have encountered a l l of the above d i f f i c u l t i e s . In an attempt to model surg ing g l a c i e r s i n the presence of growing numer ica l i n s t a b i l i t y , Mclnnes (p. 64) r e p o r t e d : 41 " D i f f e r e n t time steps g ive s l i g h t l y d i f f e r e n t surge times and per iods and there fo re at the same growth t ime , p r o f i l e s are not d i r e c t l y comparable. A l s o due to the d i f f e r e n t number of s teps , and the smoothing scheme not being p e r f e c t , the lower the time s tep the more i t e r a t i o n s , which leads to more smoothing which tends to lower e i t h e r the depth or the base v e l o c i t y , and there fore a f f e c t s the surg ing t i m e s . " and f u r t h e r , d i s c u s s i n g a c r i t e r i o n to in t roduce automatic smoothing at the appearance of i n s t a b i l i t y problems (p. 64) : "Us ing t h i s t e s t before automatic smoothing, l e s sens the number of t imes smoothing i s used, and there fo re decreases the e f f e c t smoothing has on the exact p r o f i l e s . " There i s o b v i o u s l y l i t t l e cause for optimism i n the b e l i e f that t h i s numer ica l model , for i n s t a n c e , was p r o v i d i n g a s o l u t i o n tha t c l o s e l y matched the t rue s o l u t i o n at a l l t imes . The r e s u l t s of the Mclnnes study were p u b l i s h e d by Budd (1975) and by Budd and Mclnnes (1974; 1978; 1979). I mention t h i s example, not to c r i t i c i z e any p a r t i c u l a r i n d i v i d u a l s , but to i l l u s t r a t e the l ack of a t t e n t i o n p a i d to t h i s s e r ious ques t ion by most members of the g l a c i o l o g i c a l m o d e l l i n g community. Even a major g l a c i e r m o d e l l i n g program has apparent ly had se r ious d i f f i c u l t i e s w i t h accuracy , c o n s i s t e n c y , and mass conserva t ion (see Appendix 16), yet model v e r i f i c a t i o n has not been given p r i o r i t y d i s c u s s i o n i n the p u b l i s h e d l i t e r a t u r e . Because numer ica l model r e s u l t s cannot be v e r i f i e d for the c o m p l i c a t e d problems the numer ica l models are c rea ted to s o l v e , i t i s impera t ive that numer ica l schemes be v e r i f i e d by comparing t h e i r numer ica l s o l u t i o n s w i t h a n a l y t i c a l s o l u t i o n s fo r s impler problems before the numer ica l models are used for new problems. There i s always a temptat ion w i t h a new model to rush i n t o the s o l u t i o n of compl ica ted g l a c i o l o g i c a l problems. T h i s urge to 42 break new ground i n a hurry must be c o n t r o l l e d u n t i l the numer ica l model has been demonstrated to work a c c u r a t e l y on known ground. Only i n t h i s way can there be any conf idence i n the r e s u l t s . 2.2 CONTINUITY EQUATION PROFILE MODEL 2.2.1 INTRODUCTION I have w r i t t e n a FORTRAN IV computer program which n u m e r i c a l l y so lves the c o n t i n u i t y equat ion ( 1 . 3 . 5 ) , when prov ided w i t h a subrout ine to c a l c u l a t e the v e l o c i t y V ( x , t ) averaged through the i c e t h i c k n e s s . The flow equat ion (1 .4 .38) based on G l e n ' s flow law (1 .4 .22) i s used for g l a c i e r s i m u l a t i o n s (Chapter 3 ) . I d e s c r i b e the computat iona l scheme i n d e t a i l i n Appendix 1. I summarize the computat iona l aspects i n t h i s S e c t i o n , i n c l u d i n g the numer ica l scheme (Sec t ion 2 . 2 . 2 ) , the boundary c o n d i t i o n s ( Sec t ion 2 . 2 . 3 ) , numer ica l s t a b i l i t y (Sec t ion 2 . 2 . 4 ) , and accuracy (Sec t ion 2 . 2 . 5 ) . . • Inputs to the model are bedrock e l e v a t i o n b ( x ) , mass balance A ( x , t ) , the i n i t i a l i c e t h i c k n e s s H 0 ( x ) , the i ce f l u x Q 0 ( t ) through the upslope boundary, and parameters to s p e c i f y the s l i d i n g and flow p r o p e r t i e s of the deforming medium. Output from the model i s the t i m e - v a r y i n g i c e t h i c k n e s s p r o f i l e . 43 2 .2 .2 THE NUMERICAL SCHEME Complete d e t a i l s of the numer ica l scheme are g iven i n Appendix 1, Sec t ions A1.1 and A 1 . 2 . I so lve the c o n t i n u i t y equat ion (1 .3 .5 ) by a f i n i t e d i f f e r e n c e method (e.g_. Richtmyer and Morton, 1967). The p a r t i a l d i f f e r e n t i a l equat ion (1 .3 .5 ) i s approximated by a set of a l g e b r a i c equat ions for the i c e t h i c k n e s s (h- | j=1, J } at a set of mesh p o i n t s at equal h o r i z o n t a l i n t e r v a l s of Ax . S t a r t i n g from an i n i t i a l c o n d i t i o n { h ^ l j ^ j j } , the s o l u t i o n i s obta ined by time marching, us ing a p o s s i b l y v a r i a b l e time s tep A t . The f i n i t e d i f f e r e n c e equat ions are n+1 n n+1 n+1 n n h - h + 6 (Q - Q ) + (1-6) (Q - Q ) _ j i " j + 1/2 j - 1 / 2 j + 1/2 j - 1 / 2 At W Ax W Ax j j j j n+1 n = 6A + (1-6)A (2 .2 .1) j j 1 £ j < J 1 < n < N where s u p e r s c r i p t s i n d i c a t e the time s t ep , and s u b s c r i p t s i n d i c a t e the s p a t i a l mesh p o i n t . Mass balance A- and i c e t h i c k n e s s h- are measured normal to the bed, and the mesh w increments AXj are measured a long the bed. The weight f ac tor 6 i s a constant between zero and u n i t y , used to s t a b i l i z e the scheme. I d i s c u s s 6 f u r t h e r i n S e c t i o n 2 . 2 . 4 . The i c e f l u x Q j +1/2 between the meshpoints i s r e l a t e d to the i c e t h i c k n e s s hj at the meshpoints , the channel wid th Wj + 1 /z and the v e r t i c a l l y averaged downslope v e l o c i t y v j • 1/2 between the meshpoints by 44 Q j ± l / 2 = V W (h J ± 1 / 2 J ± 1 / 2 j±1 + h ) i (2 .2 .2) 2 Most of the v a r i a b l e s i n (2 .2 .2 ) are shown i n F igure 2 . 1 . The f l u x i s c a l c u l a t e d midway between mesh p o i n t s because of numer ica l s t a b i l i t y c o n s i d e r a t i o n s . The set of equat ions (2 .2 .1 ) can be w r i t t e n i n matr ix form as where the components of the v e c t o r H are the unknown i ce t h i c k n e s s e s at the mesh p o i n t s at the fu ture time step {h-* 1 | j = 1,J} , the r i g h t s ide vec tor D conta in s q u a n t i t i e s from J the prev ious time s tep , and the matr ix M c o n t a i n s c o e f f i c i e n t s i n v o l v i n g v e l o c i t y at the future time s t ep . S ince the v e l o c i t y i s r e l a t e d to the i c e t h i c k n e s s (Sect ions 1.3 and 1.4) , t h i s makes (2 .2 .3 ) n o n l i n e a r , and the system must be so lved i t e r a t i v e l y . For the f i r s t i t e r a t i o n , the v e l o c i t y p r o f i l e at the prev ious time s tep {V? | j=3/2, 5/2,...J+1/2} i s used as an es t imate of { 0 VP* 1 | j = 3/2 , 5/2,...J+1/2} , to c a l c u l a t e the f i r s t e s t imate of the i c e t h i c k n e s s { 0 h?* 1 | j = 1 ,J } . P r e s c r i p t s i n d i c a t e i t e r a t i o n number. These i c e t h i c k n e s s es t imates are then used to c a l c u l a t e an es t imate { , Vp + 1 | j = 3/2 , 5 /2 , . . . J+1/2} of the l o n g i t u d i n a l v e l o c i t y p r o f i l e at the fu ture s t ep . The r e s i d u a l s (2 .2 .4 ) M H = D (2 .2 .3 ) 45 n+1 n+1 n+1 n n r = 2p 6[Q - Q ] + 2p (1-6)[Q - Q ] j " j j+1/2 j - 1 / 2 j 3+1/2 j - 1 / 2 n+1 n n+1 n + h - h - 6A At - (1-6)A At j j j j 1 ^ j < J 1 < n < N (2 .2 .4 ) where p = At/(2AxW ) (2 .2 .5 ) j j then measure the degree to which the c u r r e n t es t imates of i c e t h i c k n e s s and v e l o c i t y f a i l to s a t i s f y the c o n t i n u i t y equat ion ( 1 . 3 . 5 ) . By us ing a l i n e a r i z e d equat ion (2 .2 .6 ) to r e l a t e r e s i d u a l s to i c e t h i c k n e s s c o r r e c t i o n s 6h- r e q u i r e d to make the r e s i d u a l s go to z e r o , J n+1 \-~\ br** 1 r = > —j 6h (2 .2 .6 ) j *—* dh* + 1 k k=1 k or A 6h = r (2 .2 .7 ) I o b t a i n e s s e n t i a l l y a m u l t i - d i m e n s i o n a l f o r m u l a t i o n of the Newton-Raphson method (e.cj . Carnahan and o t h e r s , 1969, p . 319) to so lve ( 2 . 2 . 4 ) . The i t e r a t i o n s terminate when the l a r g e s t r e s i d u a l i n abso lute va lue i s sma l le r than a preset t e s t c r i t e r i o n . I d i s c u s s the cho ice of a c r i t e r i o n i n Appendix 11. 4 6 2 . 2 . 3 BOUNDARY CONDITIONS The set of equat ions (2 .2 .1 ) w i t h (2 .2 .2 ) c o n s i s t s of J equat ions i n J+2 unknowns { h?* 1 | j = 0,J+1}. S ince equat ion (1 .3 .5 ) i s a f i r s t order d i f f e r e n t i a l e q u a t i o n , i t r e q u i r e s one boundary c o n d i t i o n s e l e c t e d from (1 .3 .10) through ( 1 . 3 . 1 2 ) . Th i s g ives one of the r e q u i r e d two e x t r a equa t ions . At an i c e d i v i d e , zero input f l u x i s modelled by i n c l u d i n g an image p o i n t h 0 at a d i s t a n c e Ax out s ide the boundary j=1, w i t h h = h O 2 V = - v 1 ' 2 3 / 2 (2 .2 .8 ) T h i s forces the surface s lope to be zero at the boundary, but the i c e t h i c k n e s s can vary w i t h t i m e . T h i s s i t u a t i o n i s i l l u s t r a t e d i n F igure 2 . 1 . For a v a l l e y g l a c i e r o r i g i n a t i n g on a bedrock s l o p e , h = 0 1 (2 .2 .9) and the i c e surface s lope can ad jus t to any a p p r o p r i a t e v a l u e . To model only a p o r t i o n of an i c e mass such that the upstream end of the model i s some d i s t a n c e downslope from the bergschrund or d i v i d e ( e . £ . Chapter 4 , where an i c e f a l l i s m o d e l l e d ) , the i c e f l u x must be s p e c i f i e d at 1/2 mesh increment above the boundary, i . . e . n Q = Q (nAt) 1 / 2 0 (2 .2 .10) The second e x t r a equat ion r e q u i r e d to balance the number of 47 FIGURE 2 . 1 . Numerica l Scheme At Ice D i v i d e . equat ions and unknowns comes from the treatment of the downslope boundary. I f I t r e a t the boundary as f i x e d i n space, and s imply a l l o w i c e to flow through i t and out of the model, I can get the equat ion by w r i t i n g h j + 1 as an a p p r o p r i a t e e x t r a p o l a t i o n of the i ce t h i c k n e s s at the meshpoints . I use the second order Newton's d i v i d e d d i f f e r e n c e s po lynomia l (e.g.. Carnahan and o t h e r s , 1969, p . 12) . I f I choose to f o l l o w the a c t u a l motion of the g l a c i e r snout , d e f i n e d by x = L ( t ) , I must keep t r a c k of L ( t ) when i t f a l l s between the meshpoints , and be ab le to add or sub t rac t p o i n t s to or from the mesh as the terminus moves. I assume tha t the terminus i s wedge shaped from the l a s t mesh p o i n t J to the snout at L ( t ) , as shown i n F i g u r e 2 . 2 . Then, I apply the p r i n c i p l e of c o n s e r v a t i o n of mass to the shaded s e c t i o n of the wedge, o b t a i n i n g the r e q u i r e d f i n a l equat ion by b a l a n c i n g the 48 4 . — A X / 2 — H FIGURE 2 . 2 . The Wedge Terminus. volume change of the wedge aga ins t i c e f l u x i n t o the wedge from upslope p l u s net mass balance on the upper s u r f a c e , dur ing each time s t e p . The d e t a i l s may be found i n Appendix 1, Sec t ion A 1 . 3 . 4 . As the terminus moves, the i c e t h i c k n e s s h^. at the upstream end, the snout p o s i t i o n L ( t ) , and the s lope of the i c e surface are a l l free to ad jus t to changes i n the f l o w . T h i s terminus model i s s i m i l a r i n many respects to that used by B i n d s c h a d l e r (unpub l i shed ) . 2 .2 .4 NUMERICAL STABILITY I f the f i n i t e d i f f e r e n c e numer ica l scheme and the mesh increments are chosen u n w i s e l y , the f i n i t e d i f f e r e n c e equat ions may admit a s o l u t i o n q u i t e d i f f e r e n t from the s o l u t i o n of the d i f f e r e n t i a l e q u a t i o n . Th i s u s u a l l y i n v o l v e s a spur ious e x p o n e n t i a l growth of some wavenumber component which q u i c k l y 49 dominates the d e s i r e d bounded, p h y s i c a l l y reasonable s o l u t i o n . Rigorous s t a b i l i t y a n a l y s i s of n o n l i n e a r equat ions i s u s u a l l y not f e a s i b l e , but s t a b i l i t y c r i t e r i a for l i n e a r i z e d analogues g e n e r a l l y g ive u s e f u l g u i d e l i n e s and i n s i g h t s i n t o s t a b i l i t y of the n o n l i n e a r forms. The f i r s t type of s t a b i l i t y problem, the l i n e a r computa t iona l i n s t a b i l i t y , i n v o l v e s the cho ice of the mesh increments At and Ax. For many systems of f i n i t e d i f f e r e n c e equa t ions , i n s t a b i l i t y can occur when the time s tep At i s too l a rge r e l a t i v e to the space s tep Ax . I f Ax/At i s much greater than the m a t e r i a l v e l o c i t y V , m a t e r i a l t r a v e l s many mesh i n t e r v a l s per time s tep , and the system tends to ' f o r g e t ' the p h y s i c a l s o l u t i o n . The von Neumann method (e.g^. Richtmyer and Morton , 1967, p . 70) , i s one s tandard s t a b i l i t y a n a l y s i s for l i n e a r or l i n e a r i z e d e q u a t i o n s . The method i n v o l v e s f i n d i n g the t r a n s f e r f u n c t i o n T(m) i n the wavenumber domain which m u l t i p l i e s the F o u r i e r t rans form of the s o l u t i o n at the prev ious time step n , to g ive the F o u r i e r t rans form of the s o l u t i o n at the future t ime s tep n+1. I f T(m) < 1 (2 .2 .11) at a l l wavenumbers m, no wavenumber can grow, so no i n s t a b i l i t y can e x i s t . A f t e r l i n e a r i z i n g (1 .3 .5 ) by s e t t i n g the v e l o c i t y V to a cons tant i n ( 1 . 3 . 7 ) , I f i n d tha t by choosing 0 ^ 1 / 2 (2 .2 .12) I o b t a i n u n c o n d i t i o n a l l i n e a r computat iona l s t a b i l i t y for any cho ice of Ax and A t . The second type of numer ica l i n s t a b i l i t y i s c a l l e d ' the n o n l i n e a r i n s t a b i l i t y ' ( e . £ . P h i l l i p s , 1959; Mesinger and 50 Arakawa, 1976, p . 35) . T h i s i s a problem which a r i s e s i n any numer ica l s o l u t i o n of a d i f f e r e n t i a l equat ion us ing a d i s c r e t e mesh, and having terms which are n o n l i n e a r i n some combinat ion of the dependent and the independent v a r i a b l e s . In ( 1 . 3 . 5 ) , dQ/dx has t h i s form. The n o n l i n e a r i t y pumps energy (squared ampl i tude of the wavenumber spectrum) from the low wavenumber end to the h igh wavenumber par t of the wavenumber spectrum of the dependent v a r i a b l e , and the a l i a s i n g (Appendix 3) due to d i s c r e t e sampling f o l d s t h i s energy back to the lower wavenumbers, where i t d i s t o r t s the s o l u t i o n . Since the n o n l i n e a r i n s t a b i l i t y has an importance which i s not wide ly recognized i n the g l a c i e r m o d e l l i n g community, I d i s c u s s i t i n d e t a i l i n Appendix 1, Sect ion A 1 . 4 . 3 . I f a f u n c t i o n i s known only at d i s c r e t e i n t e r v a l s Ax, a wellknown r e s u l t from sampling theory (see Appendix 3) i s the f ac t that i t s F o u r i e r spectrum can be found on ly up to a wavenumber m^, c a l l e d the Nyquis t wavenumber, m = rr N Ax (2 .2 .13) T h i s i s a sampling ra te of two samples per c y c l e . Wavenumbers above t h i s l i m i t are m i s i n t e r p r e t e d as lower wavenumbers (see F i g u r e A 3 . 2 ) , by being ' f o l d e d ' back i n t o the spectrum s y m m e t r i c a l l y about mN (see F i g u r e A 1 . 3 ) . I t i s easy to show (see (A1 .4 .7 ) through (A1 .4 .9 ) ) that m u l t i p l y i n g two b a n d - l i m i t e d F o u r i e r s e r i e s together g ives a product b a n d l i m i t e d to the sum of the bandwidths of the two s i g n a l s . T h i s happens w i t h the i c e f l u x Q=hV. I f both h and V 51 are always b a n d l i m i t e d to 2/3m^, t h e i r product i s b a n d l i m i t e d to 4/3m^, which i s a l i a s e d back onto the i n t e r v a l from 2/3m^ to m^, but the lower 2/3 of the spectrum remains c o r r e c t (see F i g u r e A 1 . 3 ) . I t i s e v i d e n t , t h e n , that to avo id the non l inear i n s t a b i l i t y , the a l i a s e d s i g n a l at h igh wavenumbers, i . e . above 2/3m^, must be h e a v i l y a t t e n u a t e d . At the same t ime , the a t t e n u a t i o n must not d i s t o r t the low wavenumbers which c o n t a i n i n f o r m a t i o n about the g l a c i e r . Budd and Jenssen (1975) encountered an i n s t a b i l i t y which they a t t r i b u t e d to machine roundoff e r r o r . I t h i n k i t was a c t u a l l y the n o n l i n e a r i n s t a b i l i t y . Budd and Jenssen attempted to cope w i t h the i n s t a b i l i t y by smoothing the v e l o c i t y p r o f i l e whenever i t began to o s c i l l a t e s p a t i a l l y . Mclnnes (unpubl i shed , p . 58; p . 102) used the same computer programs to s imulate surg ing at B r u a r j o k u l l , I c e l a n d . The broken curves i n F igure 2.3 (redrawn from Mclnnes , (unpub l i shed ) , p . 58) show the i n s t a b i l i t y which arose as he attempted to b u i l d up the g l a c i e r to a steady s t a t e w i t h no s l i d i n g , and no smoothing of the p r o f i l e s . In Appendix 1 S e c t i o n A 1 . 4 . 4 , I d i s c u s s the v e l o c i t y smoothing method used by Budd and Jenssen (1975), and a l s o , the a d d i t i o n of a p u r e l y numer ica l d i s s i p a t i o n term to (1 .3 .5 ) to p r e f e r e n t i a l l y damp h igh wavenumbers. I conclude that the use of e i t h e r of these methods i s hard to j u s t i f y . The two methods I use i n t h i s numer ica l model are s u p e r i o r on p h y s i c a l grounds. F i r s t , i f the flow equat ion for the continuum (e.g.. ( 1 .4 .34 ) ) depends on the l o c a l i c e sur face s lope o ( x ) , then l a r g e ampl i tude bumps i n the s o l u t i o n p r o f i l e should tend to 52 1500 o UJ X 1000 500 Vatnajokull model Plots at 50 years 200 years 250 years 10 20 DISTANCE (km) 30 AO FIGURE 2 . 3 . Example Of The N o n l i n e a r I n s t a b i l i t y . The dashed curves redrawn from Mclnnes (unpubl i shed , p . 58) show the onset of the n o n l i n e a r i n s t a b i l i t y on i c e surface p r o f i l e s at 50 year i n t e r v a l s for a f l o w l i n e on V a t n a j o k u l l ( I c e l a n d ) . The Budd-Mclnnes model d i d not c a l c u l a t e f l u x at the midpoint s of the mesh i n t e r v a l s for h ( x ) , and t h e r e f o r e was v u l n e r a b l e to the n o n l i n e a r i n s t a b i l i t y . The s o l i d curves are the 50 year p r o f i l e s from my computer model us ing parameters i n Table 2 . 1 . The Budd-Mclnnes model a l s o appears to c rea te mass (see Appendix 16). d i f f u s e out r a p i d l y , due to the p h y s i c a l p r o p e r t i e s of the medium (See Appendix 6, which r e l a t e s p e r t u r b a t i o n s of (1 .3 .5 ) to a d i f f u s i o n equat ion f o l l o w i n g Nye (1960)) . T h i s same process should a l s o e f f i c i e n t l y smooth out h i g h wavenumber i n s t a b i l i t i e s . In Appendix 13, I have shown t h a t , when the i c e 53 f l u x i s c a l c u l a t e d at the midpoint s of the mesh i n t e r v a l s , . i . e . at J ± 1 / 2 , as shown i n F igure 2 . 2 , th i s , p h y s i c a l d i f f u s i o n process i s incorpora ted i n t o the numer ica l model and c o n t r o l s the n o n l i n e a r i n s t a b i l i t y at a l l wavenumbers. For example, the broken curves i n F igure 2.3 show the onset of the n o n l i n e a r i n s t a b i l i t y i n the Budd-Mclnnes (1974) computer model (Budd, 1975; Mclnnes , u n p u b l i s h e d ) . The o r i g i n a l c a p t i o n on F igure 4.3 of Mclnnes , from which these p r o f i l e s are redrawn, was " P r o f i l e s from the V a t n a j o k u l l model at f i f t y year i n t e r v a l s , showing the increase i n the magnitude of the o s c i l l a t i o n s w i t h t ime , due to the two p o i n t f i n i t e d i f f e r e n c e approx imat ion . In t h i s case , no smoothing was u s e d . " These authors d i d not c a l c u l a t e the i c e f l u x at the midpoint s of t h e i r mesh i n t e r v a l s . As a r e s u l t , the d i f f u s i v e mechanism of G l e n ' s f low law was unable to prevent s e r ious a l i a s i n g at the Nyquis t wavelength ( i n t h i s case , 2 km). The t r i g g e r for the i n s t a b i l i t y c o u l d have been a l a rge t r u n c a t i o n e r r o r (Appendix 1, S e c t i o n A 1 . 5 . 3 ) r e s u l t i n g from the use of a forward d i f f e r e n c e at the boundary (Budd and Jenssen, 1975). The s o l i d n A s g p Ax At -n -1 bar a ms" 2 kg m" 3 m a 2 .225 1.0 9.8 910. 1000. 1.0 TABLE 2 . 1 . Parameters for V a t n a j o k u l l (F igure 2.3) curves are the 50 year p r o f i l e s u s ing my computer model w i t h the parameters i n Table 2 . 1 . As far as I can t e l l , Mclnnes a l s o used 54 these v a l u e s . The n o n l i n e a r i n s t a b i l i t y i s removed at a l l wavelengths due to the c h o i c e of numer ica l scheme. As w e l l as the h igh wavenumber o s c i l l a t i o n , the Budd-Mclnnes model appears to s u f f e r from a mass c o n s e r v a t i o n e r r o r . I d i s c u s s t h i s f u r t h e r i n Appendix 16. A second method must be used when the i c e f l u x i s not a f u n c t i o n of the l o c a l i c e s l o p e . For example, the g r a v i t a t i o n a l s t r e s s may be c a l c u l a t e d us ing an in termedia te or l a rge sca le s lope (e.a.. B i n d s c h a d l e r , u n p u b l i s h e d , p . 92 ) . In t h i s case , I remove the n o n l i n e a r i n s t a b i l i t y wi thout d i s t o r t i n g the low wavenumber s i g n a l at a l l , by, at the complet ion of each time s t ep , t a k i n g the F o u r i e r t rans form of both the v e l o c i t y p r o f i l e and the i c e t h i c k n e s s p r o f i l e , and m u l t i p l y i n g by the lowpass f i l t e r i n F i g u r e 2 . 4 ; t h i s f i l t e r has a c u t o f f at 2/3m^. I then perform the inver se F o u r i e r t rans form to o b t a i n the p r o f i l e s which are b a n d l i m i t e d at 2/3 m^, and unaf fec ted by a l i a s i n g . The n o n l i n e a r i n s t a b i l i t y cannot grow, and cannot a f f e c t the s o l u t i o n below 2/3m v . By a s u i t a b l y smal l cho ice of the mesh increment Ax , m.. can be made l a r g e enough so that a l l p h y s i c a l l y i n t e r e s t i n g wavenumbers i n the g l a c i e r p r o f i l e are w e l l below the c u t o f f wavenumber. P h i l l i p s (1959), who o r i g i n a l l y i d e n t i f i e d the n o n l i n e a r i n s t a b i l i t y , suppressed i t i n t h i s manner. However, the procedure was q u i t e c o s t l y to implement, because h i s work predated the Fast F o u r i e r Transform a l g o r i t h m (Cooley and Tukey, 1965). In my work, u s ing the f i l t e r of F i g u r e 2.4 w i t h the Fast F o u r i e r Transform method d i d not d r a m a t i c a l l y increa se the computer execut ion time for the t e s t s I performed. 55 IT(m)l m AX TT/2 2T T /3 FIGURE 2.4. F i l t e r To Suppress The Nonlinear I n s t a b i l i t y , The Nyquist wavenumber i s at mAx=ir. 2.2.5 ACCURACY Having achieved a stable scheme, I must now ask how accurately i t solves the p a r t i a l d i f f e r e n t i a l equation. In Section A1.5, I examine the accuracy of the numerical scheme by two somewhat complementary methods. The f i r s t method follows from the von Neumann s t a b i l i t y analysis in Section 2.2.4. At a l l wavenumbers m, I compare both amplitude and phase of the transfer function T(m) of the the numerical scheme with those of the p a r t i a l d i f f e r e n t i a l eqation. Errors in amplitude represent incorrect attenuation, and errors in phase represent incorrect propagation speeds and dispersion. For the li n e a r equation, the conclusion i s to use 6 = 1 / 2 (2.2.14) for the most accurate amplitudes (see Figure A1.7), and to take Ax as small as feasibly possible to get accurate phase (see 56 F igure A1.8). I use these c o n c l u s i o n s as a guide when s e l e c t i n g parameters for the n o n l i n e a r model . The second method i s an e s t i m a t i o n of the ' t r u n c a t i o n e r r o r ' . Th i s i s the d i f f e r e n c e between an exact s o l u t i o n h ( x , t ) of the d i f f e r e n t i a l e q u a t i o n , and an exact s o l u t i o n {h M | j = l , J } of the f i n i t e d i f f e r e n c e equa t ions . I t i s a s p a t i a l domain e r r o r e s t i m a t e , i . e . i t es t imates the t o t a l e r r o r at each x p o s i t i o n , r a ther than e s t i m a t i n g the e r r o r i n each s ine wave component. A f t e r assuming that i c e t h i c k n e s s h ( x , t ) and i c e f l u x Q(x , t ) are i n f i n i t e l y d i f f e r e n t i a b l e , the f i n i t e d i f f e r e n c e s o l u t i o n can be expressed as t runca ted Tay lor expansions about the s o l u t i o n h ( x , t ) and Q(x , t ) of the d i f f e r e n t i a l e q u a t i o n . The e r r o r i s expressed i n terms of the f i r s t neg lec ted d e r i v a t i v e s . The r e s u l t i s that the t r u n c a t i o n e r r o r has the form J n € = At j (1-26) * 2h 2 at 2 1 a 3 h n 1 b 3Q n + A t 2 — + A x 2 — 6 a t 3 j 24 dx 3 j ( 2 . 2 . 1 5 ) The l e a d i n g term vanishes w i t h the cho ice 6 = 1 / 2 ( 2 . 2 . 1 6 ) T h i s i s the same r e s u l t obta ined from the wavenumber a n a l y s i s ( 2 . 2 . 1 4 ) . The t r u n c a t i o n e r r o r i s a l s o minimized by keeping the mesh increments as sma l l as p o s s i b l e . The c o e f f i c i e n t s can be es t imated from the t h i r d d e r i v a t i v e s of h and Q i n the numer ica l model to get a q u a n t i t a t i v e es t imate of the e r r o r . 57 2.3 TESTING THE CONTINUITY MODEL 2.3 .1 INTRODUCTION Two t e s t s are d e s c r i b e d i n t h i s subchapter . The f i r s t v e r i f i e s tha t the model s a t i s f i e s c o n t i n u i t y w i t h nonconstant w i d t h , ba lance , bed, v e l o c i t y , and i ce t h i c k n e s s , and that the terminus moves c o r r e c t l y . The second t e s t v e r i f i e s that the programs a c c u r a t e l y so lve a r e a l i s t i c n o n l i n e a r problem i n which the f low v e l o c i t y depends on both h ( x , t ) and i t s g rad ient i n x . Al though other t e s t s of numer ica l models are a v a i l a b l e ( e . £ . Waddington, 1979), the t e s t s presented here are a s i m p l e , reasonably comprehensive, and s t r i n g e n t t r i a l of any numer ica l f l o w l i n e model based on the c o n t i n u i t y equat ion ( 1 . 3 . 5 ) . I t h i n k that i t i s a reasonable proposa l tha t these t e s t s be used as a minimum standard for v e r i f i c a t i o n and comparison of a l l numer ica l models of g l a c i e r f low p r i o r to attempts to use them to model compl ica ted i c e masses. 2 .3 .2 CONTINUITY TEST WITH TERMINUS MOTION The a b i l i t y of the g l a c i e r model terminus to move c o r r e c t l y , i . . e . at a ra te c o n s i s t e n t w i t h the f low law being used, i s c r i t i c a l to accura te s i m u l a t i o n of g l a c i e r f low ( e . £ . N y e ( l 9 6 3 [ a ] , I963[b] ) ; Nye i n d i s c u s s i o n of a paper by Mahaffy and Andrews (1976)) . I f the model terminus advances too s l o w l y , i t w i l l ac t as a dam, r e s u l t i n g i n a g l a c i e r s o l u t i o n which i s too t h i c k , s low, and s h o r t , even though c o n t i n u i t y i s s a t i s f i e d everywhere. I f the terminus of the model advances more r a p i d l y 58 h s • s w w c AX At 0 0 0 0.1 -1 .0 0.01 1 .0 1 .0 -0 .02 0.01 0.5 0.1 -0.1 -0.01 1.0 1 .0 -0.02 0.01 0.5 TABLE 2 . 2 . Parameters for c o n t i n u i t y t e s t w i t h moving terminus . The f i r s t l i n e g ives the advancing model, and the second l i n e g ives a r e t r e a t i n g model . See F i g u r e 2.5 and F i g u r e 2 . 6 . than the c o r r e c t r a t e , the g l a c i e r s o l u t i o n w i l l be too t h i n , f a s t f l o w i n g , and extended. No s imple approximat ions are a v a i l a b l e to s i m p l i f y the s t r e s s equat ions (1 .4 .8 ) through (1 .4 .10) near the g l a c i e r t e rminus . Nye (1967) p u b l i s h e d a s o l u t i o n for the shape of a g l a c i e r t e rminus . S ince t h i s s o l u t i o n assumed steady s t a t e , a h o r i z o n t a l bed, and p e r f e c t l y p l a s t i c i c e , i t i s not s u i t a b l e for i n c l u s i o n i n my numerica l model . The terminus model being t e s t e d i s d e s c r i b e d i n S e c t i o n 2 . 2 . 3 and A 1 . 3 . To check that the numer ica l model s a t i s f i e s c o n t i n u i t y everywhere and advances or r e t r e a t s c o r r e c t l y when i t i s a l lowed to choose i t s own terminus p o s i t i o n , I use a c h a n n e l , mass ba l ance , and a f low law g i v i n g an a n a l y t i c a l s o l u t i o n to the c o n t i n u i t y equat ion (A1.1 .1) and show that the model reproduces t h i s s o l u t i o n a c c u r a t e l y . My c h o i c e i s to some extent a r b i t r a r y . Other s o l u t i o n s c o u l d e a s i l y be found. However, the one used below i s a good one because i t i s ev ident at a g lance whether the model r e s u l t s are c o r r e c t , and i t t e s t s the model w i t h a number of nonconstant and n o n t r i v i a l input f u n c t i o n s . Let h 0 , So* s, W 0 , W , and c be c o n s t a n t s , l e t the channel width W(x) be 59 Distance FIGURE 2 . 5 . C o n t i n u i t y Test w i t h Moving Terminus. Flow i s to the r i g h t . The s o l u t i o n p r o f i l e s are shown at i n t e r v a l s of 10 time u n i t s up to 40 u n i t s , then at 5 u n i t s . The model i s d e s c r i b e d by equat ions (2 .3 .1 ) through (2 .3 .5 ) w i t h the constants i n Table 2 . 2 . W(x) = W +" W x o (2 .3 .1 ) and the mass balance A(x) be A(x) = sx + c/W(x) (2 .3 .2 ) I use a v e l o c i t y (averaged over depth) g iven by ( 2 . 3 . 3 ) . Th i s has no p h y s i c a l s i g n i f i c a n c e ; i t i s merely a numer ica l t e s t . 60 V ( x , t ) = c x / [ W ( x ) ( h ( x , t ) - h )] 0 x*0 = c / [ s ( t )W.(x) ] x = 0 (2 .3 .3 ) where the surface s lope s i s g iven by s ( t ) = s + s t 0 (2 .3 .4 ) (During the numerica l s o l u t i o n procedure , care should be taken to a v o i d the s i n g u l a r i t y i n (2 .3 .3 ) when h ( x , t ) approaches h 0 ; t h i s can be done by approximat ing (2 .3 .3 ) by a s u i t a b l y smooth a n a l y t i c a l f u n c t i o n i n the reg ion of h ( x , t ) near h 0 and h ( x , t ) > h 0 . ) S u b s t i t u t i o n of (2 .3 .1 ) through (2 .3 .5 ) i n t o (A1 .1 .1 ) v e r i f i e s that the t h i c k n e s s s o l u t i o n h ( x , t ) i s h ( x , t ) = h + s ( t ) x 0 (2 .3 .5 ) D i s t ance x i s measured a long the g l a c i e r bed. Any bed p r o f i l e may be used i n the numer ica l model . The t h i c k n e s s h 0 at x=0 i s constant for a l l t ime , and h ( x , t ) v a r i e s l i n e a r l y w i t h x at a l l t imes , r e s u l t i n g i n a wedge-shaped " g l a c i e r " w i t h s lope s ( t ) . T h i s s lope s ( t ) changes l i n e a r l y w i t h time at the constant ra te s. By s o l v i n g (2 .3 .5 ) for the va lue of L such that h ( L ( t ) , t ) = 0 , the c o r r e c t g l a c i e r terminus p o s i t i o n i s found to be at L ( t ) = - h / s ( t ) (2 .3 .6 ) o F i g u r e 2.5 shows the numer ica l r e s u l t s for the advancing model us ing the cons tants i n the f i r s t l i n e of Table 2 . 2 , and the c u r v i l i n e a r bed shown i n F i g u r e 2 . 5 . The wedge s o l u t i o n advances to the r i g h t and i s shown at equal time i n t e r v a l s of 10 nondimens iona l i zed u n i t s for the f i r s t 40 u n i t s , and at 5 u n i t s 61 FIGURE 2 . 6 . C o n t i n u i t y Test With Moving Terminus. R e s u l t s i n F igure 2.5 w i t h bed e l e v a t i o n removed and d i s t a n c e measured a long the bed. t h e r e a f t e r . F i g u r e 2.6 shows the same r e s u l t s w i t h the bed e l e v a t i o n s u b t r a c t e d , and d i s t a n c e x measured a long the model bed r a t h e r than h o r i z o n t a l l y , i . e . i t should be the s o l u t i o n h ( x , t ) i n ( 2 . 3 . 5 ) . In f a c t , i t reproduces (2 .3 .5 ) to w i t h i n one par t i n 1 0 3 . The same degree of accuracy i s ob ta ined w i t h the r e t r e a t i n g model (Table 2.2) which d u p l i c a t e s the curves of F i g u r e 2 . 6 , but i n the reverse o r d e r . T h i s v e r i f i e s that the numer ica l model s a t i s f i e s c o n t i n u i t y everywhere and moves the terminus c o r r e c t l y under q u i t e genera l c o n d i t i o n s ; the width 62 v a r i e s w i t h x , the mass balance v a r i e s w i t h x , and the v e l o c i t y v a r i e s w i t h h , x , and t . 2 . 3 . 3 CONTINUITY TEST WITH BURGERS' EQUATION The t e s t i n S e c t i o n 2 .3 .2 v e r i f i e s that the model works c o r r e c t l y w i t h genera l geometr i ca l input and a moving terminus . In t h i s s e c t i o n I show that the i t e r a t i v e procedure i n the numer ica l model works a c c u r a t e l y by c o r r e c t l y s o l v i n g a f u l l y n o n l i n e a r equat ion w i t h a r e a l i s t i c form of v e l o c i t y depending on both i c e t h i c k n e s s and s l o p e . The problem so lved a l s o i n c l u d e s k inemat ic waves. The theory of propagat ion of shock waves i n a gas w i t h d i f f u s i o n has been i n v e s t i g a t e d by many a u t h o r s . C o n t r i b u t o r s to the l i t e r a t u r e of gas dynamics have i n c l u d e d Stokes (1848), Rankine (1870) and Tay lo r (1910). Methods d i scus sed i n the t ex t on n o n l i n e a r waves by Whitham (1974) are c l o s e l y fo l l owed here . One s tandard approach to f i n d the gas d e n s i t y as a f u n c t i o n space and time i n a shock f r o n t i s to so lve a c o n t i n u i t y equat ion analogous to (1 .3 .5 ) (but w i t h no source term on the r i g h t hand s ide ) together w i t h an equat ion analogous to a flow law r e l a t i n g gas f l u x Q to gas d e n s i t y H. One of the s imples t f l u x r e l a t i o n s i s Q = Q (H) - v oH (2 .3 .7 ) The f i r s t term g ives the tendency of f l u x to increase w i t h t h i c k n e s s , and the second term i s d i f f u s i v e damping w i t h d i f f u s i o n c o e f f i c i e n t i/>0. When the v a r i a b l e change 63 dQ (H) ( 2 .3 .8 ) c(H) = — dH i s in t roduced to the c o n t i n u i t y equat ion (1 .3 .5 ) w i t h mass balance A ( x , t ) set to zero and width W(x) set to u n i t y , the r e s u l t i s a n o n l i n e a r d i f f u s i o n equat ion known as Burgers ' equat ion (Burgers , 1948). When the f l u x law (2 .3 .7 ) has the form the n o n l i n e a r Cole-Hopf t r a n s f o r m a t i o n (Co le , 1951; Hopf, 1950), reduces (2 .3 .9 ) to a l i n e a r d i f f u s i o n e q u a t i o n , to which the a n a l y t i c a l s o l u t i o n i s we l l -known. A p p l y i n g the inverse t r a n s f o r m a t i o n to the s o l u t i o n g ive s the a n a l y t i c a l s o l u t i o n to Burger s ' e q u a t i o n , and thus to (1 .3 .5 ) w i t h the f l u x r e l a t i o n ( 2 . 3 . 1 0 ) . The d e t a i l s of the s o l u t i o n are g iven i n Appendix 9. S ince i t i s rare to f i n d n o n t r i v i a l a n a l y t i c a l s o l u t i o n s to n o n l i n e a r p a r t i a l d i f f e r e n t i a l equat ions , t h i s r e s u l t i s remarkable . I t p rov ides an e x c e l l e n t o p p o r t u n i t y for an exact t e s t of the model w i t h a n o n l i n e a r f low law. Burger s ' equat ion has appeared p r e v i o u s l y i n the g l a c i o l o g i c a l l i t e r a t u r e i n papers on kinematic waves of f i n i t e ampl i tude (Johnson, 1968; L i c k , 1970; H u t t e r , 1980). When the i n i t i a l c o n d i t i o n i s = v a 2 c(x,t) bx 7 (2 .3 .9) Q(H) = oH 2 + pH + r - vdH dx (2 .3 .10) c ( x , 0 ) = A 6(x) (2 .3 .11) and the boundary c o n d i t i o n s are 64 c( co , t ) = c ( - oo,t) = 0 (2.3.12) the s o l u t i o n i s (see Appendix .9) V X 2 / ( 4 l / t ) R t e (e - 1 ) c ( x , t ) = (2.3.13) J ? + ( e R - 1) J OO _ z 2 e dz xA/TTTF The parameter R i s equal to h/2v. Th i s - so lu t ion (2.3.13) i s a s i n g l e asymmetric decaying hump which propagates i n the p o s i t i v e x d i r e c t i o n . A shock ( v e r t i c a l f r o n t ) tends to form on the l e a d i n g edge, but i s prevented from doing so by d i f f u s i o n . When o=1/2 and 0=7=0 , we ob ta in from (2.3.8) c ( x , t ) = H ( x , t ) (2.3.14) so that the s o l u t i o n (2.3.13) i s a l s o the s o l u t i o n for H ( x , t ) . I t i s worth not ing that equat ion (2.3.9), as w e l l as being a n o n l i n e a r t e s t , i s a l s o a t e s t of • k inemat ic wave behaviour ( L i g h t h i l l and Whitham, 1955). For a p p l i c a t i o n to waves on g l a c i e r s , see Nye (1960). I t i s r e a d i l y apparent from (2.3.9) that the d i f f u s i v e proper ty of c ( x , t ) d c ( x , t ) = i / d 2 c ( x , t ) (2.3.15) dt dx"2" i s c a r r i e d as a k inemat ic wave at v e l o c i t y dx = c ( x , t ) (2.3.16) dt which i s thus the k inemat ic wave v e l o c i t y . (The l e f t s ide of (2.3.15) i s a t o t a l d e r i v a t i v e . ) t h i s k inemat ic wave behaviour i s i n c l u d e d i n the complete s o l u t i o n (2.3.13). To t e s t the numer ica l s o l u t i o n of (1.4.1) aga ins t the a n a l y t i c a l s o l u t i o n (2.3.13) and (2.3.14), I s e t the m a t e r i a l 65 1 .0 FIGURE 2 . 7 . Nonl inear Test With Burger s ' E q u a t i o n . Dashed l i n e s i n d i c a t i n g the a n a l y t i c a l s o l u t i o n (2 .3 .13) have been superimposed on the numer ica l model s o l u t i o n for a s i n g l e hump. Because of the c l o s e agreement (one part i n 1 0 3 ) , the curves are i n d i s t i n g u i s h a b l e i n t h i s F i g u r e . P l o t s are at i n t e r v a l s of 2 time u n i t s . The parameters of the model are g iven i n Table 2 . 3 . v e l o c i t y V ( x , t ) to be, u s ing ( 2 . 3 . 1 0 ) , V ( x , t ) = Q ( x , t ) H ( x , t ) H*0 = 0 H=0 (2 .3 .17) I so lve the f i n i t e d i f f e r e n c e equat ions on the i n t e r v a l [ - L , L ] , choos ing L s u f f i c i e n t l y l a r g e that H(-L , t )=0 adequately approximates ( 2 . 3 . 1 2 ) . The c o n d i t i o n at +co i s not needed s ince 66 A V L to o e y Ax At 1.0 0.1 7.5 2.0 1/2 0.0 0.0 0. 125 0.05 TABLE 2 . 3 . Parameters for Burger s ' equat ion t e s t . (1 .4 .1 ) i s f i r s t o r d e r . S ince f i n i t e d i f f e r e n c e equat ions cannot represent the impulse i n i t i a l c o n d i t i o n ( 2 . 3 . 1 1 ) , I s t a r t i n s t e a d w i t h an i n i t i a l c o n d i t i o n (2 .3 .13) c ( x , t 0 ) at some sma l l t o >0 . In F i g u r e 2 . 7 , the a n a l y t i c a l s o l u t i o n (2 .3 .13) i s superimposed as dashed l i n e s on the numer ica l model r e s u l t s . The parameters are g iven i n Table 2 . 3 . The agreement i s so c l o s e ( g e n e r a l l y to b e t t e r than three f i g u r e s ) that the dashed l i n e s cannot be d i s t i n g u i s h e d . T h i s t e s t shows that the numer ica l model conserves mass w i t h f u l l y n o n l i n e a r e q u a t i o n s . The i t e r a t i v e scheme for n o n l i n e a r i t i e s works c o r r e c t l y . S ince the numer ica l s o l u t i o n f i n d s the c o r r e c t t ime response for the hump, the behaviour of k inemat ic waves i n the numer ica l s o l u t i o n i s a l s o c o r r e c t . 67 2.4 THE ICE TRAJECTORY COMPUTER MODEL 2.4.1 INTRODUCTION I have w r i t t e n a FORTRAN IV computer program which l o c a t e s the t r a j e c t o r i e s of s p e c i f i e d p a r t i c l e s of i c e , as they flow through a t i m e - v a r y i n g g l a c i e r . I have d e s c r i b e d the model i n d e t a i l i n Appendix 2. The input s to the model are the bedrock topography, the cons tant s for G l e n ' s f low law ( 1 . 4 . 2 2 ) , and the number N, and the i n i t i a l p o s i t i o n s , of the i c e p a r t i c l e s to be t r a c k e d . The t r a j e c t o r y model uses the same mesh increments Ax and At as the c o n t i n u i t y model ( Sec t ion 2 . 2 ) , and the same assumptions about the channel geometry. At each time s tep , the i c e t h i c k n e s s p r o f i l e {h-| j=1,J} and the basa l s l i d i n g p r o f i l e { us. | j=1,J} J J used by the" c o n t i n u i t y model are a l s o used as input to t h i s model . 2 .4 .2 THE VELOCITY AND DISPLACEMENT FIELDS The t r a j e c t o r y of an i c e p a r t i c l e i s g iven by P ( t ) , i t s displacement vec tor as a f u n c t i o n of t i m e . For a p a r t i c l e at p o s i t i o n P 0 at time t 0 , where v ( x , t ) i s the i c e v e l o c i t y f i e l d ( u , w , v ) . To so lve ( 2 . 4 . 1 ) , I s t a r t at each time s tep by f i n d i n g the v e l o c i t y f i e l d t P ( t ) P (2 .4 .1 ) o 68 FIGURE 2 . 8 . Meshpoints For Ice V e l o c i t y C a l c u l a t i o n s . v ( x , t ) throughout the v e r t i c a l plane i n the g l a c i e r c e n t r e l i n e , on the x-z mesh shown i n F i g u r e 2 . 8 , and d e s c r i b e d i n Appendix 2 , Sec t ion A 2 . 1 . F i r s t , the downslope v e l o c i t y component u ( x , z ) i s found at each meshpoint us ing the i n t e g r a t e d form (1 .4 .34) of G l e n ' s f low law, w i t h h ( x ) , o ( x ) , and u g ( x ) g iven by the c o n t i n u i t y model s o l u t i o n . The l o n g i t u d i n a l s t r a i n ra te du/dx i s es t imated by a f i n i t e d i f f e r e n c e u ( i + 1 , j ) - u ( i - 1 , j ) o u ( i , j ) = dx DX + DX i - 1 , j i , j (2 .4 .2 ) of the downslope v e l o c i t y from ( 1 . 4 . 3 4 ) , and the l a t e r a l s t r a i n ra te dw/^y i s g iven i n terms of u ( x , z ) and the channel width W(x) by ( 1 . 3 . 3 ) . The i n c o m p r e s s i b i l i t y c o n d i t i o n w i t h mass c o n s e r v a t i o n 69 g ives dv = - 51 du - dw d"x 5y (2 .4 .3 ) Using the approximat ion ( u s u a l l y very good, e.g_. R o t h l i s b e r g e r , 1972) tha t basa l m e l t i n g i s n e g l i g i b l e as far as mass balance i s concerned, i . e . (2 .4 .3 ) i s i n t e g r a t e d n u m e r i c a l l y by Simpson's r u l e (e.c[. Carnahan and o t h e r s , 1969, p . 73) from the bed to l e v e l z to g ive v ( x , z , t ) . The l a t e r a l v e l o c i t y component w ( x , y , z ) i s e x a c t l y zero on the f l o w l i n e down the centre of the channe l , and i s very smal l i n a narrow flow volume (Figure 1.1) cent red on t h i s f l o w l i n e . The average value of w i s zero i n t h i s volume. The computer model uses w=0. When the v e l o c i t y f i e l d has been complete ly determined at the meshpoints i n F i g u r e 2 . 8 , the displacement f i e l d of the i c e l e a v i n g each meshpoint P 0 and going to a p o i n t P i n a time i n t e r v a l At i s found by e s t i m a t i n g the i n t e g r a l (2 .4 .1 ) by i . e . , the average of the v e l o c i t i e s at the beg inning and end of the time i n t e r v a l A t . The v e l o c i t y v (P , t+At ) i s e s t imated from the va lues of v at the four surrounding meshpoints at time t+At us ing an i n t e r p o l a t i o n scheme d e s c r i b e d i n Appendix 2, S e c t i o n A 2 . 3 . 1 . v(x , 0,t) = 0 (2 .4 .4 ) P(t+At) - P 0 ( t ) = 1 v ( P 0 , t ) 2 L At (2 .4 .5 ) 70 2 .4 .3 THE ICE PARTICLE TRAJECTORIES At each time s tep , and for each of the N p a r t i c l e s being t r a c k e d , the coord ina te s (6x ,6z) r e l a t i v e to a meshpoint ( i , j ) are recorded . The di sp lacements of the p a r t i c l e s i n the time At are i n t e r p o l a t e d from the d i sp lacements of i c e at the four surrounding meshpoints , and the new coord ina te s of the p a r t i c l e s are then saved. The program checks whether the new p o s i t i o n s are s t i l l w i t h i n the g l a c i e r mass, and saves the i n t e r p o l a t e d times and p o s i t i o n s at which i c e p a r t i c l e s reach the i c e s u r f a c e . By us ing At<0, the program i s e a s i l y adapted to t r ack p a r t i c l e s backward i n time and upslope (e.g.. from a b o r e h o l e ) , to f i n d where and when they entered the i c e mass as p r e c i p i t a t i o n . 2 .4 .4 ACCURACY OF THE TRAJECTORY MODEL The v e l o c i t y f i e l d i n t h i s model i s obta ined by assuming that shear ing p a r a l l e l to the g l a c i e r bed i s the dominant component of de format ion . The downslope v e l o c i t y u ( x , z ) i s then found us ing the s t r e s s equat ions and the mechanical p r o p e r t i e s of i c e . T h i s i s the only p lace where the force equat ions are used. The other v e l o c i t y components are determined p u r e l y k i n e m a t i c a l l y , u s ing c o n t i n u i t y and geometric assumptions about the l a t e r a l v a r i a t i o n of the f low f i e l d . When the assumption of predominant ly shear deformation h o l d s , t h i s approach works very w e l l . The f r a c t i o n a l e r r o r i n l o n g i t u d i n a l s t r a i n r a t e , and i n the v e r t i c a l v e l o c i t y , i s q u i t e 71 s m a l l . The l e a d i n g term i s approximate ly the r a t i o of the unbalanced l o n g i t u d i n a l fo rces on a v e r t i c a l column, to the basa l shear f o r c e . The d e t a i l s are g iven i n Appendix 7, equat ion ( A 7 . 5 . 9 ) , which i s repeated here as ( 2 . 4 . 6 ) . e(x) = 0 dtf' dtf' 2h xx + h yy 6x 6x />gho + (n-1) 6 T a xz d v /ou o~x/ dz max (2 .4 .6) When the assumption of predominantly shear deformation i s not a c c u r a t e , the downslope v e l o c i t y component u ( x , z ) i s r e l a t i v e l y i n a c c u r a t e . The f i n a l term i n the e r r o r est imate (2 .4 .6 ) may be of order u n i t y , or l a r g e r , when du/oz i s s m a l l . The second term i n the same equat ion may a l s o be of order u n i t y , or l a r g e r , when s t r e s s d e v i a t o r s other than the shear <JXZ c o n t r i b u t e to the e f f e c t i v e s t r e s s . Furthermore, the presence of any a d d i t i o n a l s t r e s s components other than * x Z always softens the i c e , so there i s no p o s s i b i l i t y of compensating approximat ions i n t h i s term. T h i s s i t u a t i o n may a r i s e at an i c e d i v i d e , where the predominant motions may be v e r t i c a l s i n k i n g , and l o n g i t u d i n a l e x t e n s i o n as the i c e f lows downslope i n both d i r e c t i o n s . F o r t u n a t e l y , the i c e v e l o c i t y i s very smal l at an i c e d i v i d e , so that the t o t a l e r r o r i n the t r a j e c t o r i e s i s not l a r g e . To be s a f e , the t r a j e c t o r i e s near an i c e d i v i d e should be i n t e r p r e t e d only q u a l i t a t i v e l y . L o n g i t u d i n a l s t r e s se s are a l s o important i n i c e f a l l s . The" use of G l e n ' s f low law, and t h e . i n c l u s i o n of 72 l o n g i t u d i n a l s t r a i n ra te s ( to the approximat ion of ( 2 . 4 . 6 ) ) i s an improvement over the t r a j e c t o r y models of Nagata (1977), who assumed that du/dz was z e r o , and Dansgaard and Johnsen (I969[a]) who set du/dz to a constant at each x i n the lowest 400 metres , and zero above that l e v e l , and assumed that du/dx was constant for a g iven z , i n a model for the Camp Century , Greenland b o r e h o l e . The fac t that my ' model i s time dependent i s a l s o an improvement over these models . One advantage of t h i s model i s tha t the mass c o n s e r v a t i o n law i s s t i l l obeyed g l o b a l l y and l o c a l l y at a l l t i m e s . T h i s has not been the case w i t h some other t r a j e c t o r y models. For i n s t a n c e , the Weertman (1968) a n a l y t i c a l model for i c e v e l o c i t y and temperature at Camp Century , Greenland , assumed a constant v e r t i c a l s t r a i n r a t e , but used a h o r i z o n t a l v e l o c i t y g iven independently by i n t e g r a t i n g G l e n ' s f low law (1 .4 .22) to get a r e s u l t l i k e ( 1 . 4 . 3 4 ) . Dansgaard and Johnsen (1969[b]) showed that t h i s v i o l a t i o n of c o n t i n u i t y was the l i k e l y cause of a d i sc repancy of 2°C between p r e d i c t e d and measured temperature at the bed at Camp Century , Green land . The f low model used by Dansgaard and Johnsen (1969[a]) s a t i s f i e d c o n t i n u i t y everywhere, but assumed the form of the h o r i z o n t a l v e l o c i t y component, in s t ead of us ing G l e n ' s f low law. Whi le the model I d e s c r i b e d i n Sec t ion 2.2 does not yet i n c l u d e temperature v a r i a t i o n s , that i s a s imple development that w i l l not v i o l a t e mass c o n s e r v a t i o n . Budd and o t h e r s , (1971, p . 42) assumed a constant v e r t i c a l s t r a i n ra te dv/dz independent of depth , at each p o s i t i o n x , 73 r e g a r d l e s s of the va lues of the other terms i n ( 2 . 4 . 3 ) . They a l s o d e s c r i b e d another model (p. 43) i n which the v e r t i c a l s t r a i n ra te was weighted by the downslope v e l o c i t y u ( x , z ) at each p o i n t . While t h i s makes the v e r t i c a l s t r a i n ra te curves more c l o s e l y resemble the expected shape i n a r e a l i ce mass, i t s t i l l does not s a t i s f y c o n t i n u i t y l o c a l l y through ( 2 . 4 . 3 ) . An incompres s ib le continuum f lowing w i t h such a v e l o c i t y f i e l d would have to l o c a l l y c r e a t e or des t roy mass. Constant v e r t i c a l s t r a i n ra te (or a s t r a i n ra te s p e c i f i e d a p r i o r i ) may be a reasonable approximat ion when d e r i v i n g a n a l y t i c a l s o l u t i o n s to some flow problems to i l l u s t r a t e the p h y s i c s i n v o l v e d (e.cj . R o b i n , 1955; f o r the e f f e c t of advec t ion on temperature p r o f i l e s ) , and i t a r i s e s n a t u r a l l y i n some simple k inds of f l o w , i_.e. h o r i z o n t a l shear ing r e s t r i c t e d to the basa l l a y e r (e.c=. H i l l , 1950, p . 233; Nye, 1951; or Nye, 1957). However, g iven the s o p h i s t i c a t i o n of c u r r e n t computer models a c c e p t i n g otherwise q u i t e genera l i n p u t s , t h i s assumption now appears to be a needless l i m i t a t i o n . Whi le the e r r o r s i n v o l v e d may t u r n out , i n some s i t u a t i o n s , to be s m a l l , i t i s p r e f e r a b l e to a v o i d them at no a d d i t i o n a l c o m p l i c a t i o n . The only model of which I am aware which attempts to use the c o n t i n u i t y equat ion i n a manner s i m i l a r to the way I d e s c r i b e d i n S e c t i o n 2 .4 .2 to f i n d the v e r t i c a l v e l o c i t y component, and a l s o makes no a p r i o r i assumptions about the flow f i e l d other than those d e s c r i b e d i n S e c t i o n 2 . 4 . 2 , i s the p r e l i m i n a r y t h r e e - d i m e n s i o n a l model of Jenssen (1977) for the Greenland i c e cap. T h i s model a l s o so lve s for the temperature d i s t r i b u t i o n , for use i n the temperature dependent flow law. The 74 Jenssen model i s s t i l l under development to improve the boundary treatment and and to reduce inaccuracy r e s u l t i n g from the coarse g r i d imposed by computer l i m i t a t i o n s (Greenland was modelled by a 12x12 map g r i d , 10 p o i n t s deep) ; however, the Jenssen approach i s the l o g i c a l next s tep that must be taken towards the numer ica l s o l u t i o n of the complete set of equat ions for i c e sheet s . 2.5 TESTING THE TRAJECTORY MODEL 2.5 .1 INTRODUCTION In S e c t i o n 2 .5 , I d e s c r i b e two t e s t s of the t r a j e c t o r y model . In S e c t i o n 2 . 5 . 2 , I compare the steady s t a t e t r a j e c t o r i e s and v e l o c i t y f i e l d c a l c u l a t e d by my computer model to an a n a l y t i c a l s o l u t i o n for a steady s t a t e i c e sheet (Nagata,1977). In S e c t i o n 2 . 5 . 3 , I use the balance c o n d i t i o n (A5.8) at the i ce sur face to show that c o n t i n u i t y i s s a t i s f i e d by the v e l o c i t y f i e l d . 2 .5 .2 NAGATA ICE SHEET TEST Nagata (1977) d e r i v e d an a n a l y t i c a l s o l u t i o n for the sur face p r o f i l e , v e l o c i t y f i e l d , and s t r e a m l i n e s of a steady two-dimens iona l i c e sheet r e s t i n g on a f l a t bed. For the steady s t a te ca se , s t r e a m l i n e s and i c e t r a j e c t o r i e s c o i n c i d e . I d e s c r i b e the Nagata model i n d e t a i l i n Appendix 15, S e c t i o n A 1 5 . 3 . Nagata (1978) used t h i s i c e sheet model to 75 E CD O c o o CD 1.0 0.5 Nagata (1977) Model Steady state ice sheet Basal sliding m=3 Velocity f i e l S and streamlines •6723 0.5 x/L 0 -2 1.0 FIGURE 2.9. Nagata Steady Ice Sheet. The a n a l y t i c a l s o l u t i o n f o r the i c e t h i c k n e s s , mass balance, p a r t i c l e paths, and v e l o c i t y f i e l d u s i ng the parameters i n equation (2.5.2). The v e l o c i t i e s have been m u l t i p l i e d by 250 y e a r s . The numbers are the time i n years f o r i c e to flow the l e n g t h of each p a r t i c l e path. e x p l a i n the c o n c e n t r a t i o n of m e t e o r i t e s at the M e t e o r i t e Ice F i e l d i n A n t a r c t i c a . The Nagata model assumes that the forward i c e v e l o c i t y u(x) i s constant throughout a v e r t i c a l column and depends on the b a s a l shear s t r e s s (1.4.25) through a Weertman-type r e l a t i o n (see Appendix 8, S e c t i o n A8.3.1) of the 76 T X (km) FIGURE 2 .10 . Growth Of Nagata Ice Sheet . The i c e surface e l e v a t i o n i s shown at i n t e r v a l s of 500 y e a r s , s t a r t i n g from i c e - f r e e c o n d i t i o n s . The flow parameters, are g iven i n Table 2.4 and the mass balance i n F i g u r e 2 . 9 . The steady s t a te i c e t h i c k n e s s agrees w i t h the a n a l y t i c a l s o l u t i o n by Nagata (1977) to one par t i n 1 0 3 . form m u(x) = A T (2 .5 .1 ) Al though i t i s not c l e a r i n the o r i g i n a l paper , the Nagata model a l s o assumes that the v e r t i c a l v e l o c i t y component v ( x , z ) i s equal to a constant b a long the upper sur face of the i c e sheet . The s ca l e of the i c e sheet model i s set by H , the t h i c k n e s s at the i c e d i v i d e . Whi le the form of the mass balance and the v e l o c i t y f i e l d are not a c l o s e r e p r e s e n t a t i o n of r e a l i c e sheets , the Nagata model i s w e l l s u i t e d to t e s t i n g the accuracy of my numerica l 77 Distance (km) FIGURE 2 . 1 1 . V e l o c i t y F i e l d For Nagata Model . The v e l o c i t y v e c t o r s c a l c u l a t e d by the Waddington t r a j e c t o r y model for the steady s t a te p r o f i l e i n F i g u r e 2.10 and us ing flow parameters i n Table 2 .4 . The v e l o c i t i e s have been m u l t i p l i e d by 250 years (The u n i t s of the v e c t o r s are d i s p l a c e m e n t ) . t r a j e c t o r y program. I compare my numer ica l r e s u l t s to the Nagata model shown i n F igure 2.9 for the cons tants i n Table 2 .4 . F i g u r e 2.10 shows the growth of the Nagata i c e sheet to steady s t a t e us ing the p r o f i l e model d e s c r i b e d i n S e c t i o n 2 . 2 ) . The steady s t a t e p r o f i l e agrees w i t h the a n a l y t i c a l model i n F i g u r e 2.9 to w i t h i n one par t i n 1 0 3 . The s o l i d v e c t o r s i n F i g u r e 2.11 show the steady s t a t e v e l o c i t y f i e l d c a l c u l a t e d by my numer ica l model . The v e l o c i t i e s have been m u l t i p l i e d by a f a c t o r of 250 y e a r s . The v e l o c i t y 78 m A b H L p g b a r " 2 a'1 m a" 1 m km kg m" 3 m s" 2 2 100. 1.0 3000 454.6 910. 9.8 AX km A t a 7.215 10.0 TABLE 2 . 4 . Parameters for Nagata i c e sheet . f i e l d agrees w i t h the a n a l y t i c a l s o l u t i o n (F igure 2.9) to w i t h i n a few p a r t s i n 10 3 (except near the terminus where the agreement i s on ly to two par t s i n 1 0 2 , because the mass balance i n the numer ica l model remains f i n i t e w h i l e the a n a l y t i c a l mass balance, goes to -oo ) . F i g u r e 2.12 shows the t r a j e c t o r i e s for i c e e n t e r i n g the i c e sheet at time t=0. ; these are the same f i v e p o i n t s shown for the s t r e a m l i n e s i n F i g u r e 2 . 9 . The arrowheads i n d i c a t e 250 year i n t e r v a l s . The t o t a l res idence t imes a long each of the f i v e s t r e a m l i n e s are compared to the a n a l y t i c a l va lues i n Table 2 . 5 . The good agreement between the numer ica l r e s u l t s and the a n a l y t i c a l s o l u t i o n i n d i c a t e s tha t the v e l o c i t y f i e l d i s r e c o n s t r u c t e d a c c u r a t e l y and the i n t e g r a t i o n of the v e l o c i t y f i e l d to get f l o w l i n e s i s c o r r e c t and a c c u r a t e . 79 0 100 200 300 X (km) 400 500 FIGURE 2 .12 . T r a j e c t o r i e s In Nagata Mode l . The s o l i d curves are the p a r t i c l e paths c a l c u l a t e d by the Waddington t r a j e c t o r y model by i n t e g r a t i n g the v e l o c i t y f i e l d i n F igure 2 . 1 1 . The arrowheads i n d i c a t e 250 year i n t e r v a l s . In Table 2.5 the t o t a l res idence times i n the i c e sheet are compared to va lues for the a n a l y t i c a l s o l u t i o n . 2 . 5 . 3 SURFACE MASS CONSERVATION TEST Conserva t ion of mass at the g l a c i e r surface w i t h normal v e c t o r n i m p l i e s tha t the normal v e l o c i t y of the i c e - a i r i n t e r f a c e i s equal to the sum of the normal i c e v e l o c i t y p lu s the sur face accumulat ion r a t e , i . . e . — • n = v • n + a • n (2 .5 .2 ) dt T h i s equat ion i s d e r i v e d i n Appendix 5, S e c t i o n A 5 . 2 . When i c e 80 S t reaml ine Number Numerical model years A n a l y t i c a l model years 1 6746. 6723. 2 4625. 4606. 3 3335. 3322. 4 2350. 2346. -> 5 1477. 1466. TABLE. 2 . 5 . Residence times i n Nagata i c e sheet . t h i c k n e s s change rate d h / d t and mass balance vec tor a are measured normal to the bed, and the v e l o c i t y components (u ,v) are p a r a l l e l to and normal to the bed, (2 .5 .2) becomes oh = v ( x , h ) - u (x ,h ) djh + A(x) S t 6x (2 .5 .3 ) where A(x) i s now a s c a l a r g i v i n g the mass ba lance . My d e r i v a t i o n (Sec t ion 2 .4 .2 and Appendix 2, Sec t ion A2.2) of the v e l o c i t y f i e l d uses on ly i n c o m p r e s s i b i l i t y (2 .4 .3 ) and the basa l boundary c o n d i t i o n ( 2 . 4 . 4 ) ; (2 .5 .3 ) can- be used as an independent check of the degree to which the v e l o c i t y f i e l d conserves mass. Models which use (2 .5 .3 ) to d e r i v e the v e r t i c a l v e l o c i t y ( e . £ . Weertman, 1968; Budd and o t h e r s , 1971, p . 42) do not have t h i s c o n s i s t e n c y check. F i g u r e 2.13 (a) shows the mass balance A(x) and the surface r i s e dh/dt for a t y p i c a l time s t ep , w i t h parameters d e s c r i b e d by Table 2 . 6 , d u r i n g the growth of the Nagata i c e sheet model . I ' a p p l y the t e s t at a time (2500 years ) when the i c e sheet i s growing v i g o r o u s l y because c o n d i t i o n s at tha t t ime impose the most s t r i n g e n t c o n d i t i o n s on accuracy . F igure 2.13 (b) shows the r e s i d u a l of ( 2 . 5 . 3 ) , i . . e . the d i f f e r e n c e between the l e f t and r i g h t s ide s of the equat ion on s u b s t i t u t i n g the va lues of h , A , 81 FIGURE 2 .13 . Surface Mass Conserva t ion T e s t . Time=2500 years d u r i n g the growth of the Nagata i c e sheet model to a steady s t a t e , (a) shows the mass balance ( s o l i d curve) and the ra te of sur face r i s e dh/dt (broken c u r v e ) , (b) shows the r e s i d u a l of ( 2 . 5 . 3 ) . I t i s very smal l (note s ca l e change of 10~ 2 , i n d i c a t i n g that the v e l o c i t y f i e l d s a t i s f i e s i n c o m p r e s s i b i l i t y . The terminus i s at x=37.8 km. The r e s i d u a l increase s near the terminus because there are few meshpoints i n each v e r t i c a l column; v e r t i c a l i n t e g r a t i o n i s i n a c c u r a t e . u , and v from the numer ica l s o l u t i o n . The r e s i d u a l e r r o r i s three orders of magnitude l e s s than the average magnitude of the mass ba lance , except near the terminus where there are i n s u f f i c i e n t mesh p o i n t s i n any v e r t i c a l column to guarantee accura te v e r t i c a l i n t e g r a t i o n of ( 2 . 4 * 3 ) . T h i s reg ion has no 82 Time Ax DZ At a km m a 2500. 7.215 150. 10. TABLE 2 . 6 . Nagata i ce sheet surface boundary t e s t . e f f e c t on the t r a j e c t o r i e s i n v e s t i g a t e d i n t h i s work. S i m i l a r t e s t s (not shown) on the S tee le G l a c i e r Model 1 (Figure 3.3) r o u t i n e l y g ive r e s i d u a l s of the order of one par t i n 10 2 of the average mass ba lance , even though the mass balance, i s d i s c o n t i n u o u s ( t r i b u t a r i e s ) , the width i s v a r i a b l e , and the flow law i n c l u d e s a height-dependent l o n g i t u d i n a l s t r a i n r a t e . R e s u l t s of these t e s t s i n d i c a t e that the v e l o c i t y f i e l d s a t i s f i e s the c o n t i n u i t y equat ion (2 .3 .3) to a very good a p p r o x i m a t i o n . Th i s t e s t v e r i f i e s the accuracy of the t r a j e c t o r y model under t i m e - v a r y i n g c o n d i t i o n s . 83 CHAPTER 3: CAN STABLE ISOTOPES REVEAL A HISTORY OF SURGING? 3.1 INTRODUCTION The s t a b l e heavy i sotopes O 1 8 of oxygen and D (deuterium) of hydrogen i n g l a c i e r i c e have been wide ly used as i n d i c a t o r s of c l i m a t i c change (e . c i . Dansgaard and o t h e r s , 1969; 1971). Thi s procedure r e q u i r e s assumptions about the p a t t e r n of g l a c i e r f l o w . In t h i s chapter , I i n v e s t i g a t e a r e l a t e d problem; assuming that the past c l i m a t e i s known, can the s t a b l e i sotope d i s t r i b u t i o n be used to r e v e a l the flow h i s t o r y of t i m e - v a r y i n g g l a c i e r s ? The example I cons ider i n t h i s chapter i s the S tee le G l a c i e r , Yukon T e r r i t o r y . Th i s g l a c i e r was observed to surge i n 1966-1967 (S t an ley , 1969). I used the computer models d e s c r i b e d in Chapter 2 to f i n d the i c e surface and the v e l o c i t y f i e l d throughout the surge c y c l e , and to c a l c u l a t e the i c e t r a j e c t o r i e s . Knowing the t r a j e c t o r i e s and the i s o t o p i c compos i t ion (from c l i m a t e ) at the time and p lace the m a t e r i a l was p r e c i p i t a t e d as snow on the g l a c i e r s u r f a c e , has a l lowed me to c o n s t r u c t l o n g i t u d i n a l c r o s s - s e c t i o n s and sur face p r o f i l e s showing the i s o t o p i c d i s t r i b u t i o n at a s e r i e s of t imes dur ing the surge c y c l e ; t h i s p a t t e r n would f a c i l i t a t e the s e l e c t i o n of optimum borehole and surface sampling l o c a t i o n s for an i s o t o p i c study of past f low p a t t e r n s . In S e c t i o n 3 .2 , I d e s c r i b e the S tee le G l a c i e r and i t s surge h i s t o r y . In Sec t ions 3.3 and 3.4 I de sc r ibe the computer models I used to s imula te the S tee le G l a c i e r . Model 1 i n S e c t i o n 3.3 i s 84 a d e t a i l e d model based on a l l the a v a i l a b l e d a t a . Because of l i m i t e d s l i d i n g data and an approximat ion i n the mass ba lance , i t i s not at present the best model for p a r t i c l e t r a j e c t o r y c a l c u l a t i o n s . With more complete data and some computer model re f inements , i t may become so . Model 2 i n Sec t ion 3.4 i s a s i m p l i f i e d v e r s i o n b e t t e r matched to the r e s o l u t i o n of the s l i d i n g o b s e r v a t i o n s . Model 2 i s used for the t r a j e c t o r y and i s o t o p i c c a l c u l a t i o n s i n Sec t ions 3.6 and 3 .7 . In Sec t ion 3 .5 , I d e s c r i b e the use of s t a b l e i so topes i n g l a c i o l o g y , and present two p o s s i b l e i s o t o p i c r e l a t i o n s for p r e c i p i t a t i o n at S tee le G l a c i e r . The i c e surface p r o f i l e c a l c u l a t i o n s in Sec t ion 3.6 i n d i c a t e that a surge p e r i o d of approximate ly 100 years i s a p p r o p r i a t e for the S tee le G l a c i e r i f the present mass balance and the v e l o c i t y of the 1966-67 surge are r e p r e s e n t a t i v e of the average long- term c l i m a t e and of the g l a c i e r f low p a t t e r n . In S e c t i o n 3 .7 , I present computed l o n g i t u d i n a l c r o s s - s e c t i o n s and surface p r o f i l e s of the i s o t o p i c d i s t r i b u t i o n 6 ( 0 1 8 / 0 1 6 ) for a model w i t h a surge p e r i o d i c i t y of 97 y e a r s . I s o t o p i c d i s c o n t i n u i t i e s occur i n the i c e a long sur faces which were at the i c e - a i r i n t e r f a c e i n the accumulat ion reg ion when a surge began. Even for the i s o t o p i c - p r e c i p i t a t i o n model most favourable to the format ion of d i s c o n t i n u i t i e s , the d i s c o n t i n u i t i e s do not exceed one DEL u n i t on the S t e e l e G l a c i e r ; t h i s amount may be hidden by noi se i n the S tee le G l a c i e r environment. 85 3.2 STEELE GLACIER 3.2.1 GENERAL DESCRIPTION S tee le G l a c i e r ( 6 1 ° 1 0 ' N , 140°15'W) i s a surg ing v a l l e y g l a c i e r on the northeas t s lopes of the I c e f i e l d Ranges (see F igure 3.1) of the S t . E l i a s Mountains , Yukon T e r r i t o r y , Canada. P r i o r to 1963 i t was c a l l e d the Wolf G l a c i e r . (Sharp (1951) and Bostock (1948, p . 99) s t a t e that t h i s i s more proper than the of ten-used n o t a t i o n Wolf Creek G l a c i e r . ) . The g l a c i e r l eng th v a r i e s from 34 to 44 km, and the width of the main channel i s one to two km. The main channel f lows down from 3000 m e l e v a t i o n on the nor th s ide of M t . S tee le (F igure 3.2) to 1200 m e l e v a t i o n i n S tee le Creek, where the i c e i s stagnant and moraine-covered between surges . The c o n t i n e n t a l s lope of the I c e f i e l d Ranges i s s e m i - a r i d . Annual p r e c i p i t a t i o n drops from 300 cm a " 1 at Yakutat on the Gul f of A l a s k a (see F i g u r e 3 . 1 ) , to about 35 cm at Kluane Lake (Wood, 1972). The f i r n l i n e on the S t e e l e G l a c i e r i s p r e s e n t l y very h igh (2400 m or more) . A major source of accumulat ion i s avalanches from the n o r t h face ( l a b e l l e d (0) i n F i g u r e 3.2) of Mt . S tee le (5070 m). In i t s upper reaches , the S tee le i s an ex tens ive system of accumulat ion bas ins and converg ing t r i b u t a r i e s . These t r i b u t a r y i c e streams (see (1) through ( 5 ) , F i g u r e 3.2) a l s o c o n t r i b u t e a s u b s t a n t i a l f r a c t i o n of the S tee le G l a c i e r mass ba lance . In i t s lower reaches , the S tee le G l a c i e r makes a 90° bend to the e a s t , and enter s a s t r a i g h t and narrow v a l l e y formed by the Wolf Creek monocline (Sharp, 1943). 8 6 FIGURE 3 . 1 . I c e f i e l d Ranges L o c a t i o n Map. N o n s t i p p l e d areas are major g l a c i e r s and i c e f i e l d s . The t r i a n g l e s i n d i c a t e major summits of the I c e f i e l d Ranges (S t . E l i a s Mounta ins ) . The S tee le G l a c i e r i s on the nor th-ea s t ( c o n t i n e n t a l ) s lope (top c e n t r e ) . 3 .2 .2 GLACIER SURGES A g l a c i e r surge (e.cj . Meier and P o s t , 1969) i s a short p e r i o d of very r a p i d f l o w , d u r i n g which i c e i s t r a n s f e r r e d from an i c e r e s e r v o i r area to an i c e r e c e i v i n g area downstream. A surge i s f o l l o w e d by a longer p e r i o d of s t agna t ion and a b l a t i o n 87 i n the r e c e i v i n g a rea , w i t h renewed i c e b u i l d u p in the i c e r e s e r v o i r ; these areas do not n e c e s s a r i l y correspond to the accumulat ion and a b l a t i o n zones d e f i n e d by mass ba lance . For a l a rge v a l l e y g l a c i e r l i k e the S t e e l e , the maximum v e l o c i t y d u r i n g a surge may be 500 m a " 1 to 10 km a " 1 , w i t h downstream i c e displacement of one to 10 km. Surging appears to be a p e r i o d i c phenomenon. For g l a c i e r s l i k e , the S t e e l e , the surge d u r a t i o n i s t y p i c a l l y one to two y e a r s , w i t h a qu iescent phase l a s t i n g from 20 to 150 y e a r s . G l a c i e r s of a l l s i z e s have been observed to surge. Examples are the Trapr idge G l a c i e r ( i n the S t e e l e Creek watershed) which i s on ly three km l o n g , and the Muldrow G l a c i e r on Mt . M c K i n l e y , A la ska which i s over 50 km l o n g . There has been s p e c u l a t i o n that the A n t a r c t i c i c e sheet may surge ( H o l l i n , 1969; W i l s o n , 1969). Surging g l a c i e r s are found i n many p a r t s of the w o r l d , e«2. Alaska (Tarr and M a r t i n , 1914, p . 168), B r i t i s h Columbia and the Yukon T e r r i t o r y . (Pos t , 1969), the A r c t i c I s l ands ( H a t t e r s l e y - S m i t h , 1964; L0ken, 1969), I ce l and (Thorar in s son , 1969), the Karakoram (Hewi t t , 1969), the Pamirs , T ien Shan, Caucasus, and Kamchatka (Dolgoushin and Os ipova , 1975). Surges occur i n both temperate and c o l d or subpolar g l a c i e r s . Surging does not appear to be t r i g g e r e d by c l i m a t e v a r i a t i o n s . The h i g h v e l o c i t y d u r i n g surg ing i s g e n e r a l l y a t t r i b u t e d to r a p i d ba sa l s l i d i n g . V a r i o u s hypotheses on the mechanism of surg ing have been put forward . The more p l a u s i b l e ones i n c l u d e thermal r e g u l a t i o n (Robin , 1955; C l a r k e , 1976; L l i b o u t r y , I969[b ] ) , s t r e s s i n s t a b i l i t i e s ( e . g . P o s t , 1960; R o b i n , 1969) and ba sa l water f i l m i n s t a b i l i t i e s (Weertman, 1962, 1969; Robin 88 and Weertman, 1973; Budd, 1975). A concensus on the cause of g l a c i e r surg ing has yet to be reached. In t h i s chap te r , I i n v e s t i g a t e the e f f e c t of p e r i o d i c surg ing on the s t a b l e i sotope d i s t r i b u t i o n i n the S tee le G l a c i e r . The amount of ba sa l s l i d i n g i s important to t h i s q u e s t i o n ; the mechanism which causes i t i s n o t . For t h i s reason, I s p e c i f y a basa l s l i d i n g v e l o c i t y u s ( x , t ) c o n s i s t e n t w i t h the observa t ions of the S tee le G l a c i e r surge of 1966-1967, and I use t h i s as a boundary c o n d i t i o n for the computer model . 3 .2 .3 OBSERVATIONS OF THE STEELE GLACIER The flow p a t t e r n before and d u r i n g a surge i s b e t t e r known for the S tee le G l a c i e r than for most other surg ing g l a c i e r s . S c i e n t i f i c e x p e d i t i o n s sponsored by the American Geographica l S o c i e t y and l e d by W. A . Wood exp lored the S tee le Creek watershed i n 1935, 1936, 1939, and 1941. Work i n c l u d e d the exper imenta l use of ob l ique a e r i a l photography for mapping areas of h igh r e l i e f . These a i r photos and the panoramas from ground c o n t r o l p o i n t s g ive i n f o r m a t i o n on the i ce surface e l e v a t i o n and s t a t e of f low dur ing the quiescent phase. Sharp (1943) d e s c r i b e d the geology of S tee le Creek, observed the g l a c i e r s of the S tee le Creek bas in (Sharp, 1947), and i n t e r p r e t e d the. g l a c i a l h i s t o r y (Sharp, 1951). The Surveys and Mapping Branch of the Department of Energy, Mines and Resources , Government of Canada, obta ined v e r t i c a l a e r i a l photography of the I c e f i e l d Ranges from 10 700 m i n the summer of 1951. F l i g h t l i n e s A13232/33 covered the S tee le 89 61°20 61°05 140° 30 U0° FIGURE 3 . 2 . S tee le G l a c i e r And T r i b u t a r i e s . The c e n t r a l f l o w l i n e used i n the computer model i s marked at 2.5 km i n t e r v a l s . The Hodgson G l a c i e r i s i n c l u d e d i n the channel width f u n c t i o n (F igure 3.3 (b)) for the main i c e s tream. For Model 1, the mass c o n t r i b u t i o n s from the minor t r i b u t a r i e s (0) through (5) are inc luded i n the mass balance f u n c t i o n (F igure 3.3 ( a ) ) . G l a c i e r . T h i s photography was used to prepare the government map 115F at the s ca le 1:250,000; i t was a l s o used to prepare a map of the S tee le G l a c i e r at the s ca l e 1:25,000 (Topographica l Survey, 1967). Obl ique a e r i a l photography of the S tee le G l a c i e r was ob ta ined p e r i o d i c a l l y , between 1960 and 1965 by W. A . Wood (American Geographica l S o c i e t y ) and by A . Post (U .S . G e o l o g i c a l S u r v e y ) . In 1960, Post p r e d i c t e d an imminent surge for the 90 S tee le G l a c i e r (Wood, 1972), and i n 1965, noted the f i r s t s igns of increa sed flow on the upper S tee le (Wood, 1972). Dur ing the 1966-1967 surge , h igh ob l ique a e r i a l photos were obta ined by Wood, pre-1941 survey s t a t i o n s were r e - o c c u p i e d , and the Surveys and Mapping Branch obta ined v e r t i c a l a i r photo coverage on August 13 and September 15, 1966. Wood (1972) and S tan ley (1969) measured displacements of i d e n t i f i a b l e surface markings . Features o r i g i n a l l y between the Hodgson conf luence and the 90° bend (F igure 3.2) a l l advanced at roughly the same v e l o c i t y of 5 km a " 1 d u r i n g 1966, w i t h a t o t a l displacement of roughly 8 km when the surge ended l a t e i n 1967. There was no s i g n i f i c a n t l a t e r a l v a r i a t i o n of v e l o c i t y beyond 200 m from the g l a c i e r marg ins . Bayrock (1967) observed the d e t a i l s of the terminus advance, and the r e a c t i v a t i o n of stagnant i c e . A bulge of a c t i v e i c e 30 m h igh moved forward at about 10 m d " 1 (3650 m a " 1 ) . T h i s bulge sometimes has been c a l l e d the terminus i n the l i t e r a t u r e on the 1966-1967 surge. A l f o r d from the Whitehorse o f f i c e of the Water Survey of Canada obta ined monthly a i r photographs of the advancing bulge d u r i n g the winter of 1966-1967. He observed an advance of 6000 feet (1830 m) between September 10, 1966 and January 15, 1967 (Roots , 1967). By August 1967, the a c t i v e bulge had slowed to 2 m d " 1 (730 m a - 1 ) (Thomson, 1972). S t a n l e y (1969) i d e n t i f i e d three zones based on sur face e l e v a t i o n changes d u r i n g the surge . An upper zone, e n t i r e l y i n the accumulat ion a rea , apparent ly was not i n v o l v e d i n the 1966-1967 surge . In a middle zone (the i c e r e s e r v o i r zone) from a p o i n t above the f i r n l i n e (based on S t a n l e y ' s d e s c r i p t i o n , I 91 l o c a t e t h i s po in t at about x=7 km) to a po in t about 3 km above the 90° bend (at x=27 km i n my model ) , there was a net lower ing of the i c e sur f ace , w i t h a maximum value of 130 m above the Hodgson-Steele conf luence . In the lower remaining zone ( i c e r e c e i v i n g area) the surface rose by up to 100 m. Dur ing the winter of 1966-1967, the Hodgson G l a c i e r began a y e a r - l o n g surge dur ing which i t pushed the s t i l l - s u r g i n g S tee le i c e stream to one t h i r d of i t s normal width at the i c e s u r f a c e . The e f f e c t of the Hodgson surge on deep i c e i s unknown. The Hodgson i c e formed a l a rge lobe extending three km down the main S tee le c h a n n e l . Th i s t r i b u t a r y surge may have been t r i g g e r e d by a reduced c o n f i n i n g s t r e s s r e s u l t i n g from the r a p i d lower ing of the S tee le G l a c i e r at t h e i r c o n f l u e n c e . 3 .2 .4 PERIOD OF STEELE SURGES Sharp (1951) showed t h a t , based on 1941 o b s e r v a t i o n s , the S t e e l e G l a c i e r below the bend had been stagnant s i n c e 1916 or e a r l i e r . Based on b i o l o g i c a l r e c o l o n i z a t i o n r a te s i n s i d e the most recent t r i m l i n e , he e s t imated that the l a s t advance ended i n the p e r i o d 1840-1890, j . . e . 115 to 75 years before the 1966-1967 surge . Wood (1972) p o i n t e d out evidence for i c e d i sp lacements of over two km between 1935 and 1941 near the Steele-Hodgson c o n f l u e n c e . The p r e v i o u s l y smooth i c e i n that r eg ion was h e a v i l y crevassed i n 1941, and the crevasses looked s e v e r a l years o l d . The evidence suggests tha t the S tee le G l a c i e r had a minor surge 92 or a f a i l e d attempt at a surge i n 1938 or 1939. Sharp ' s (1951) es t imate of the time of the prev ious terminus advance would then imply a b u i l d u p time of 50 to 90 y e a r s . Sharp (1951) p o i n t e d out moraines 100 m to 150 m above the 1941 i c e surface below the bend. An a b l a t i o n est imate of 2 m a " 1 i m p l i e s the maximum i c e e l e v a t i o n occurred 50 to 75 years p r i o r to 1941 . There are no d i r e c t observa t ions on the S tee le G l a c i e r of r e g u l a r looped moraines , or of d i s r u p t e d ogive p a t t e r n s which c o u l d r e v e a l p r e - h i s t o r i c a l surge ep i sodes . The assumption of p e r i o d i c i t y i s based on o b s e r v a t i o n of surg ing g l a c i e r s elsewhere ( e . £ . V a r i e g a t e d G l a c i e r , A l a s k a , (B indschadler and o t h e r s , 1977)) . The assumption of exact p e r i o d i c i t y for the computer model i s , at be s t , an approx imat ion ; s i g n i f i c a n t c l i m a t e v a r i a b i l i t y on the time s ca l e of g l a c i e r surges has been w i d e l y documented ( e . £ . Bryson and Goodman, 1980; G r i b b o n , 1979). I f the S tee le G l a c i e r surges p e r i o d i c a l l y , the p e r i o d i s i n the range 50 to 150 y e a r s . 93 3.3 NUMERICAL MODEL 1 3.3.1 FLOW LAW CONSTANTS AND SHAPE FACTOR The computer model i n i t s present form assumes that the i ce i s i s o t h e r m a l . Temperatures i n the upper 100 metres of the S tee le G l a c i e r ( J a r v i s and C l a r k e , 1974; C l a r k e and J a r v i s , 1976) f o l l o w i n g the 1966-1967 surge were i n the range -1 °C to - 7 ° C . I t i s l i k e l y that the basa l i c e , w i t h i n which most of the shear deformation takes p l a c e , i s at or near the pressure m e l t i n g p o i n t . The assumption of temperate i c e i s not unreasonable . I use the constants (Pa ter son , 1981, p . 39) A = 5.3 1 0 - 1 5 s " 1 k P a ' 3 n=3 (3 .3 .1) i n G l e n ' s f low law (1 .4 .22) for the S tee le G l a c i e r s i m u l a t i o n . S ince most of the motion i n a surge i s due to s l i d i n g , changes i n the f low law cons tant s have l i t t l e e f f e c t on the i c e movement. The shape of the v a l l e y c r o s s - s e c t i o n i s unknown. I use a shape f a c t o r (see S e c t i o n A7.4) of s = 0.8 independent of p o s i t i o n x . S ince the shape f a c t o r i s r a i s e d to the nth power (see equat ion ( 1 . 4 . 3 4 ) ) , u n c e r t a i n t y i n s i n t r o d u c e s s u b s t a n t i a l u n c e r t a i n t y i n t o the deformation v e l o c i t y . However, as p o i n t e d out above, the i n t e r n a l de format ion of the S tee le G l a c i e r i s a smal l component of the t o t a l m o t i o n . 94 3 .3 .2 BED TOPOGRAPHY There are no p u b l i s h e d i ce depth data for the S tee le G l a c i e r . J a r v i s and C la rke (1974) and C la rke and J a r v i s (1976) reached depths exceeding 100 metres w i t h a h o t - p o i n t d r i l l wi thout any i n d i c a t i o n of bottom; the g l a c i e r i s l i k e l y much t h i c k e r than t h i s . T y p i c a l depths of l a rge v a l l e y g l a c i e r s can exceed 600 metres (e.g.. L o w e l l G l a c i e r , S t . E l i a s Mountains , based on monopulse radio-echo sounding, 1977; G. K. C. C l a r k e , per sona l communicat ion) . A rough est imate of the depth e.cj . above the confluence of the S tee le and the Hodgson can be obta ined by assuming a ba sa l shear s t r e s s of one bar and measuring the i c e surface s lope to be c=0.03 from the 1951 map (Topographica l Survey, 1969). Assuming a shape f a c t o r of s=0.8, the s t r e s s r e l a t i o n (1 .4 .25) g ives the depth as (approximately) h = xz = 472 metres s/>go (3 .3 .2 ) A l t e r n a t i v e l y , assuming flow by s imple shear w i t h no s l i d i n g , and us ing the 1951 v e l o c i t y of 25 m a " 1 measured at t h i s p o i n t by Wood (1972) (1 .4 .38) w i t h the flow law cons tant s (3 .3 .1 ) g ives the depth est imate l / ( n + D = 445 metres h = (n+2)V n 2A( Syoga) (3 .3 .3 ) S ince the motion of the S tee le G l a c i e r i s o b v i o u s l y not an example of steady n o n s l i p f low i n a c y l i n d r i c a l c h a n n e l , these are merely rough e s t imate s . The presence of ba sa l s t r e s s of l e s s than one bar , or of nonzero ba sa l s l i d i n g would r e s u l t i n an 95 overes t imate of the t rue t h i c k n e s s . N e g l e c t i n g the s t r e s s p e r t u r b a t i o n s and l a t e r a l asymmetry caused by bends i n the c h a n n e l , I assume the c e n t r a l f l o w l i n e (x a x i s ) f o l l o w s the broken curve shown i n F igure 3 . 2 . For the S tee le G l a c i e r bed topography i n F i g u r e 3.3 ( c ) , I use an e x p o n e n t i a l f u n c t i o n having the form -bx h(x) = ae + c (3 .3 .4 ) where a , b , and c are constants determined by f i t t i n g the three p o i n t s : ( 1 ) 2900 metres e l e v a t i o n at the bergschrund of the main i c e stream (x=0). (2) 1200 metres at the 1978 terminus p o s i t i o n (x=42 km). (3) 1650 metres at the conf luence of the Hodgson and S tee le G l a c i e r s (X=18 km); t h i s va lue would g ive an i c e depth of about 400 metres at t h i s p o i n t i n 1951 (Topographica l Survey, 1969)., The 1951 l o n g i t u d i n a l sur face p r o f i l e i s shown i n F i g u r e 3.3 ( c ) . The a c t u a l topography beneath the S tee le G l a c i e r i s undoubtedly more c o m p l i c a t e d . Th i s approximat ion means that the computer model cannot be expected to q u a n t i t a t i v e l y reproduce the observed i c e surface e l e v a t i o n s . 96 3 .3 .3 CHANNEL WIDTH The width of the main i c e stream of S tee le G l a c i e r can be es t imated by measuring the s epara t ion of the l a t e r a l moraine r i d g e s on the topographic map at the s c a l e 1:25,000 (Topographica l Survey, 1969). The p u b l i s h e d government mapsheet 115G and 1 1 5 F ( E l / 2 ) , at the s ca le 1 :250 ,000 , (from which F i g u r e 3.2 i s drawn), i s not a r e l i a b l e i n d i c a t i o n of the channel w i d t h . I d e a l l y , t r i b u t a r i e s should be i n c l u d e d i n the numer ica l model as separate i c e streams w i t h t h e i r own bed and mass balance f u n c t i o n s , and coupled to the main i c e stream by t h i c k n e s s and f l u x c o n d i t i o n s at the j u n c t i o n s ; t h i s op t ion i s not a v a i l a b l e i n the computer model i n i t s present form. I n s t e a d , I can i n c l u d e the e f f e c t s of t r i b u t a r i e s i n approximate ways through the width or mass balance f u n c t i o n s . S ince the d i scharge and depth of the Hodgson G l a c i e r are probably comparable t o those of the S tee le at t h e i r c o n f l u e n c e , I use the sum of t h e i r widths above t h i s p o i n t . T h i s i s a s imple a p p r o x i m a t i o n . In f a c t , t h e i r surface g r a d i e n t s are not equal everywhere i n t h i s r e g i o n , t h e i r mass balance f u n c t i o n s probably d i f f e r , and the two g l a c i e r s do not surge together (the Hodgson l a s t surged i n 1967-68, one year a f t e r the S t e e l e ) . Because of t h e i r much sha l lower depths , to i n c l u d e the minor t r i b u t a r i e s (1) through (5) (F igure 3.2) i n t h i s way would tend to g r o s s l y reduce the t h i c k n e s s and v e l o c i t y of the main c h a n n e l . The mass c o n t r i b u t i o n of these t r i b u t a r i e s i s i n c l u d e d i n Model 1 by an a d d i t i o n to the mass balance f u n c t i o n (F igure 3.3 (a)) at the reg ions of c o n f l u e n c e . 97 3 .3 .4 MASS BALANCE No comprehensive mass balance measurements have been made on the S t e e l e G l a c i e r . Between an average date of J u l y 19, 1974, and an average date of J u l y 2, 1975, C o l l i n s (unpubl ished) measured a b l a t i o n at e ight survey t a r g e t s near the r i g h t angle bend (x=30 km). The measurements ranged from 1.57 m to 2.77 m of i c e , w i t h a mean of 2.15 m. A l l o ther evidence i s i n d i r e c t . On the lower S tee le G l a c i e r , S tan ley (1969) es t imated the a b l a t i o n to be 1.5 m a " 1 based on downwasting of i ce shown to be stagnant by Sharp (1951). S ince s n o w f a l l can depend s t r o n g l y on l o c a l c o n d i t i o n s , e x t r a p o l a t i n g mass balance i n f o r m a t i o n even from adjacent v a l l e y s i s at best a r i s k y procedure . However, the q u a l i t a t i v e pa t t e rns of mass ba lance , and i t s order of magnitude throughout the I c e f i e l d Ranges g ive some c o n t r o l on reasonable es t imates for the S t e e l e G l a c i e r . Marcus and Ragle (1970) repor ted winter accumulat ion measurements on a t r a v e r s e across the I c e f i e l d Ranges from the lower Seward G l a c i e r to the Kaskawulsh G l a c i e r . The va lues for the Kaskawulsh are i n s t r u c t i v e because both Kaskawulsh and S tee le G l a c i e r s occupy s i m i l a r p o s i t i o n s on the east s ide of the I c e f i e l d Ranges (F igure 3 . 1 ) . The p r e c i p i t a t i o n on the Kaskawulsh increases w i t h e l e v a t i o n . The 1964-65 winter accumulat ion was 1.7 metres of i c e e q u i v a l e n t at 2640 metres e l e v a t i o n , decreas ing to 0.35 metres of i c e e q u i v a l e n t at 1615 metres e l e v a t i o n . These authors a l s o reported n e g l i g i b l e summer m e l t i n g on the i c e p l a t e a u at D i v i d e S t a t i o n 30 km west of the Kaskawulsh G l a c i e r . D i v i d e S t a t i o n at 2620 metres 98 1 1 r ol 1 i L 10001 1 1 0 10 20 30 40 50 x (km) FIGURE 3.3..Model 1 For Steele G l a c i e r . (a) Mass balance. The dashed curve i s balance for the main ice stream. The s o l i d curve includes mass contributions from the t r i b u t a r i e s (0) through (5) i d e n t i f i e d in Figure 3.2. (b) Glacier width. The broken curve is the width of the main ice stream only. The s o l i d curve includes the width of the Hodgson Glacier . (c) Bed topography. The ice surface p r o f i l e in 1951 i s also shown. 99 e l e v a t i o n rece ive s about 2.2 metres ( i c e e q u i v a l e n t ) net accumulat ion each year ; the record there i n d i c a t e s tha t 1965-66 was a low s n o w f a l l year by about 40% i n the Kaskawulsh reg ion (Marcus and Ragle , (1970, F i g u r e 7 ) ) . K e e l e r (1969) repor ted t h a t , on Mount Logan, 60 km south of the S tee le G l a c i e r , e l e v a t i o n has l i t t l e e f f e c t on p r e c i p i t a t i o n above 2500 metres e l e v a t i o n . The net .accumulat ion i s about 0.8 metres of i ce a n n u a l l y . S tan ley (1969, F igure 1) l o c a t e d the f i r n l i n e at about 2400 metres on the S tee le G l a c i e r based on the 1951 a e r i a l photography. Wood (1972) put the f i r n l i n e i n the higher e l e v a t i o n range of 2750 to 2900 metres . T h i s would leave almost no accumulat ion a r e a . Sharp (1947) gave the es t imate of 8000 to 9100 feet (2440 to 2775 metre s ) . The broken l i n e i n F i g u r e 3.3 (a) shows an e s t imated mass balance f u n c t i o n fo r the main i c e stream of the S t e e l e G l a c i e r c o n s i s t e n t w i t h C o l l i n s (unpubl ished) and w i t h the i n d i r e c t o b s e r v a t i o n s . The i c e f l u x from the i t h sma l l t r i b u t a r y g l a c i e r i s i n c l u d e d through a p e r t u r b a t i o n 6b^ to the mass balance f u n c t i o n . In Appendix 19 I d e s c r i b e two methods of e s t i m a t i n g the i c e f l u x from each t r i b u t a r y and how I use t h i s to es t imate the 6 b £ . The s o l i d curve i n F i g u r e 3.3 (a) shows the mass balance w i t h the t r i b u t a r y terms 6 b ^ i n c l u d e d . Whi le F igure 3.3 (a) represent s the best a v a i l a b l e es t imate of the S tee le G l a c i e r mass ba lance , i t should be kept i n mind tha t mass balance may change s i g n i f i c a n t l y w i t h longterm c l i m a t i c change. The decades 1930-1950, upon which much of the data fo r the S tee le G l a c i e r i s based, appear t o have been an 100 e x c e p t i o n a l l y warm p e r i o d (§_.£. Hansen and o t h e r s , 1981; Schneider and Mesirow, 1976, Chapter 3 ) . 3 .3 .5 CYCLIC SURGE PATTERN FOR THE MODEL Because I am i n v e s t i g a t i n g consequences of s u r g i n g , r a ther than surge mechanisms, I s p e c i f y a p r i o r i the s l i d i n g v e l o c i t y u 5 ( x , t ) . T h i s aspect of surge m o d e l l i n g i s d i s cus sed i n Sec t ion 1 .5 .3 . I have chosen to use a s l i d i n g v e l o c i t y having the form u ( x , t ) = X ( x , t ) T ( t ) (3 .3 .5 ) s because i t can represent the observed s l i d i n g of the S tee le G l a c i e r reasonably w e l l . Other f u n c t i o n a l forms of u s ( x , t ) c o u l d e q u a l l y w e l l f i t the observa t ions of S tan ley (1969) and Wood (1972). The time dependent term T ( t ) i s a nondimensional we ight ing f a c t o r between zero and u n i t y . F i g u r e 3.4 (b) shows the form of T ( t ) for a surge c y c l e of l e n g t h t « . Dur ing the quiescent s tage , T has some smal l constant va lue f ( £ . £ . f=0 g ive s no s l i d i n g ) . The surge s t a r t s at time t 0 and the v e l o c i t y r i s e s to the peak value by time t 1 f remains at the maximum u n t i l t 2 , then f a l l s back to the nonsurge l e v e l by time t 3 . The term X ( x , t ) shown i n F igure 3.4 (a) g ive s the normal ized s p a t i a l d i s t r i b u t i o n of the s l i d i n g v e l o c i t y at a time ( t - t 0 ) . Observat ions on the S tee le G l a c i e r ( S t a n l e y , 1969) and on other surg ing g l a c i e r s , and some t h e o r e t i c a l work on surg ing ( e . g . R o b i n , 1969; Robin and Weertman, 1973) suggest that r a p i d s l i d i n g s t a r t s i n a sma l l r e g i o n , and the boundaries of t h i s 101 FIGURE 3 .4 . S l i d i n g Model For S tee le G l a c i e r . (a) X ( x , t ) g ives the normal ized s p a t i a l dependence of the s l i d i n g v e l o c i t y at time t . (b) T ( t ) i s the temporal we ight ing f u n c t i o n for the s l i d i n g v e l o c i t y i n (a) . zone of r a p i d s l i d i n g then propagate down (and p o s s i b l y up) the g l a c i e r . U 0 i s the maximum s l i d i n g v e l o c i t y d u r i n g the surge . In t h i s model , each t r a n s i t i o n from a zone of r a p i d s l i d i n g to a zone of no s l i d i n g i s g iven by one h a l f c y c l e of a c o s i n e . The d i s t r i b u t i o n of s l i d i n g can change d u r i n g the surge as the four p o i n t s x . 2 ( t ) , x . , ( t ) , x , ( t ) , and x 2 ( t ) move at v e l o c i t i e s c . 2 ( t ) , c . , ( t ) , c , ( t ) , and c 2 ( t ) r e s p e c t i v e l y . The data from the surge of the S tee le G l a c i e r are not s u f f i c i e n t l y complete to a l l o w d e t a i l e d e s t i m a t i o n of the c - ( t ) . I use constant va lues fo r these four v e l o c i t i e s , so that the v e l o c i t y t r a n s i t i o n p o i n t s move a c c o r d i n g to 1 02 x ( t ) = x (t ) + ( t - t ) c 2 2 0 0 2 (3 .3 .6 ) w i t h s i m i l a r equat ions for the other three p o i n t s . The constants x . 2 ( t 0 ) , x . , ( t 0 ) , x , ( t 0 ) , and x 2 ( t 0 ) de f ine the extent of the t r i g g e r zone. Observat ions by Raymond and o ther s (unpubl i shed , F i g u r e 9) on the Var i ega ted G l a c i e r , A l a s k a , i n d i c a t e a regu lar increase of s l i d i n g v e l o c i t y , from 0.05 m d " 1 (18 m a " 1 ) i n 1973, to 0.3 m d " 1 (110 m a " 1 ) i n 1979 i n the upper reaches of the g l a c i e r . The V a r i e g a t e d G l a c i e r appears to have a surge p e r i o d of about 20 year s , and i s expected to surge sometime i n the mid-1980 ' s . The s i m u l a t i o n s which I have c a r r i e d out for t h i s chapter do not i n c l u d e a pre-surge increase i n basa l s l i d i n g , a l though i t probably c o u l d be model led s a t i s f a c t o r i l y by the Weertman (1957) s l i d i n g mechanism; t h i s o p t i o n i s a v a i l a b l e i n the computer program. However, f o r the S t e e l e G l a c i e r , the observed pre-surge v e l o c i t i e s (Wood, 1972) are low, and are not s u f f i c i e n t l y d e t a i l e d to warrant t h i s a d d i t i o n a l c o m p l i c a t i o n . The l e a d i n g edge of the zone of r a p i d s l i d i n g ( i . e . x , ( t ) to x 2 ( t ) ) i s a reg ion of compressive f l o w . The i n s t a b i l i t y of reg ions of compressive f low i s w ide ly recognized ( e . g . P a t e r s o n , 1969, p . 207) . For nonsurging g l a c i e r s , the surface r i s e r e s u l t i n g from compressive flow i n the a b l a t i o n area i s l a r g e l y balanced by sur face m e l t i n g . A surge l a s t i n g only one or two years occurs too q u i c k l y for m e l t i n g to have any a p p r e c i a b l e c o n t r o l on surface e l e v a t i o n ; to a v o i d a growing bulge or shock wave moving w i t h the l e a d i n g edge of the zone of r a p i d s l i d i n g , the v e l o c i t i e s c , and c 2 of the v e l o c i t y d i s tu rbance must be 103 s u f f i c i e n t l y greater than the m a t e r i a l v e l o c i t y U 0 that the v e l o c i t y change can move ahead of any bulge forming i n the i c e . For the elementary case of i c e w i t h s l i d i n g v e l o c i t y U 0 and t h i c k n e s s h 0 advancing i n t o t h i n n e r stagnant i c e of t h i c k n e s s h , i n a channel of constant w i d t h , a s imple c o n t i n u i t y argument i n Appendix 18 shows that when c , = c 2 , the v e l o c i t y d i s tu rbance must move at l e a s t at the speed c £ U i o to prevent the development of a shock f ron t i n the advancing v e l o c i t y t r a n s i t i o n . Equat ion (3 .3 .7 ) can be used to es t imate the i c e t h i c k n e s s h 0 and h , ; for the S tee le G l a c i e r , I have used i t o n l y to choose reasonable va lues of c , and c 2 . 3.4 STEELE GLACIER MODEL 2 3.4.1 PROBLEMS WITH STEELE MODEL 1 I ran the S tee le G l a c i e r model 1 (F igure 3.3) us ing the flow law parameters (3 .3 .1 ) w i t h no s l i d i n g . The numer ica l parameters are given i n Table 3 . 1 . The sur face p r o f i l e s at 50 year i n t e r v a l s s t a r t i n g from i c e - f r e e c o n d i t i o n s are shown i n F i g u r e 3 . 5 . The i n d i v i d u a l t r i b u t a r i e s coa le sced by 250 y e a r s . The maximum v e l o c i t y (averaged over depth) i n steady s t a te i s 33 m a " 1 at x=12.5 km. The corresponding v e l o c i t y at the i ce sur face i s (see Sec t ion 1.4.4) (n + 2) / (n+l ) t imes t h i s , j . . e . 41 m a " 1 . The f i n a l steady s t a t e l e n g t h i s 35.5 km, and the h - h o 1 (7. 7. 7 ) 104 n A 5 g fi Ax At -n -1 kPa s ms" 2 kg irr 3 m a 3 5.3 1 0 " 1 5 0.8 9.8 900. 500. 1 .0 TABLE 3 . 1 . Parameters for S tee le steady s t a t e . AOOOl T Steele Glacier Growth to steady state No sliding Profiles every 50 years 3000 2 0 0 0 h 1000 0 10 20 30 x (km) FIGURE 3 . 5 . S tee le Model 1 Growth To Steady S t a t e . There was no s l i d i n g . The model parameters are g iven i n Table 3 . 1 . The i c e surface p r o f i l e s are shown at 50 year i n t e r v a l s s t a r t i n g from i c e - f r e e c o n d i t i o n s . The i c e lobes at 11 km and 13 km are caused by i c e from t r i b u t a r i e s (1) and ( 2 ) . maximum depth i s 444 m at x=25 km. The sur face p r o f i l e s from Model 1 appear to be reasonable , 105 and the model coped w i t h the r a p i d l y - v a r y i n g mass balance and width f u n c t i o n s i n F i g u r e 3 .3 . However, I c a l c u l a t e d some of the steady s t a t e s t r eaml ine s for Model 1; these are shown i n F i g u r e 3.6 (note that the ab sc i s s a i s i ce t h i c k n e s s r a ther than e l e v a t i o n ) . The v e r t i c a l mesh increment for the v e l o c i t y 1 1 r Steele Glacier Model 1 Steady streamlines No sliding Arrows every 100 years FIGURE 3 .6 . Steady Sta te S t reaml ines For Model 1. T h i s c r o s s - s e c t i o n shows i c e dep th . Arrows on the s t r eaml ine s i n d i c a t e each 100 years of f l o w . There was no s l i d i n g . The p e r t u r b a t i o n s to the smooth s t r eaml ine s are due to the a d d i t i o n of i c e f l u x from t r i b u t a r i e s through mass balance terms. c a l c u l a t i o n s (see F i g u r e 2.8) was DZ = 15m The p e r t u r b a t i o n s i n the s t r e a m l i n e s are caused by the a d d i t i o n s of t r i b u t a r y i c e f l u x through the mass balance f u n c t i o n . I am l o o k i n g for i r r e g u l a r i t i e s i n s t r e a m l i n e shape and i n i sotope d i s t r i b u t i o n r e s u l t i n g from s u r g i n g ; to have fea tures of t h i s 106 £ co o c a a CO. u) a 3 2 1 0 -1 •2 3 0 0 0 h £ 2 0 0 0 h •g £ 1000 3000 2000 CXI X 1000 Steele Glacier Model 2 Growth to steady state No sliding Profiles every 50 years 0 10 20 30 AO 50 (km) FIGURE 3 . 7 . Model 2 For S tee le G l a c i e r . (a) Mass ba lance . (b) G l a c i e r w i d t h . (c) Bed topography and i c e sur face p r o f i l e s at 50 year i n t e r v a l s d u r i n g growth to steady s t a te w i t h no s l i d i n g . form in t roduced through the mass balance f u n c t i o n i s 107 u n d e s i r a b l e . Adding the t r i b u t a r y f l u x at the g l a c i e r surface ra ther than at the channel margins i s an adequate approximat ion i f on ly the i ce surface shape i s d e s i r e d . I f p a r t i c l e t r a j e c t o r i e s are a l s o d e s i r e d , t h i s approximat ion i s unacceptab le . An a d d i t i o n a l problem a r i s e s w i t h Model 1 when s l i d i n g i s added. The width f u n c t i o n (F igure 3.3 (b)) i s a d e t a i l e d r e p r e s e n t a t i o n of the v a l l e y of S tee le Creek. The s p a t i a l d i s t r i b u t i o n of the surg ing v e l o c i t y of the S tee le G l a c i e r i s not known w i t h the same r e s o l u t i o n ; Wood (1972) and S tan ley (1969) obta ined only s p a t i a l l y and tempora l ly averaged v e l o c i t i e s . F o r c i n g the i c e to s l i d e at n e a r l y constant v e l o c i t y through a channel of h i g h l y v a r i a b l e width can r e s u l t i n some u n r e a l i s t i c surface c o n f i g u r a t i o n s . For i n s t a n c e , i n reg ions where the channel width g rad ient i s l a rge and n e g a t i v e , the i c e can t h i c k e n r a p i d l y and o b t a i n a reversed surface s l o p e . More d e t a i l e d data on how the g l a c i e r a c t u a l l y changes speed to prevent t h i s s i t u a t i o n are not a v a i l a b l e for the 1966-67 surge . I t i s necessary to use a width f u n c t i o n which has the same degree of s p a t i a l d e t a i l as the s l i d i n g d a t a . When a network computer model for t r i b u t a r i e s i s developed, and when more d e t a i l e d s l i d i n g o b s e r v a t i o n s are a v a i l a b l e , the amount of d e t a i l i n Model 1 can be j u s t i f i e d . 108 3 .4 .2 SIMPLIFICATIONS F i g u r e 3.7 (a) and (b) show a s i m p l i f i e d model for the S tee le G l a c i e r . Thi s model resembles Model 1 ' i n i t s gross f e a t u r e s , yet avoids the d i f f i c u l t i e s I d e s c r i b e d i n the prev ious s e c t i o n . F i g u r e 3 .7(c) shows the growth of Model 2 to steady s ta te w i t h no s l i d i n g , us ing the parameters i n Table 3 . 1 , and s t a r t i n g from i c e - f r e e c o n d i t i o n s . The steady s t a t e l e n g t h i s 36.5 km, the maximum i c e t h i c k n e s s i s 452 m at x=21 km, and the maximum v e l o c i t y (averaged over depth) i s 30.4 m a " 1 at x=11 km. These va lues are c l o s e to the va lues for Model 1 (see S e c t i o n 3 . 4 . 1 ) . 3.5 STABLE ISOTOPES IN GLACIOLOGY 3.5.1 DEFINITION OF THE DEL SCALE The s tandard method of d e s c r i b i n g the i s o t o p i c compos i t ion of oxygen and hydrogen i n water i s the DEL s ca l e ( 6 ) . The r a t i o R of the c o n c e n t r a t i o n s of the heavy and l i g h t i so topes b 1 8 / 0 1 6 and D/H can be measured w i t h a mass spectrometer ; a p r a c t i c a l c o n c e n t r a t i o n sca le should be based on R. The 6 va lue (3 .5 .1 ) of an i c e sample i s the r e l a t i v e d i f f e r e n c e between R^ of the sample and R of a re ference sample known as Standard Mean Ocean Water (SMOW) ( C r a i g , 1961). 1 0 9 R - R 6 = S SMOW 1 0 3 R SMOW -I ( 3 . 5 . 1 ) A drawback of t rue SMOW i s the f ac t that no samples are a v a i l a b l e . Samples of other 0 1 8 / 0 1 6 i s o t o p i c s tandards from the U . S . N a t i o n a l Bureau of Standards have been d i s t r i b u t e d i n the past by the I n t e r n a t i o n a l Atomic Energy Agency, V i e n n a . In September 1 9 7 6 , at V i e n n a , the C o n s u l t a n t s ' Meet ing on S tab le Isotope Standards and I n t e r c a l i b r a t i o n in Hydrology and Geochemistry set up a s tandard sample c a l l e d Vienna SMOW. The d i f f e r e n c e between true SMOW ( C r a i g , 1 9 6 1 ) and Vienna SMOW i s - 0 . 0 5 ° / o o - Th i s d i f f e r e n c e , i s not s i g n i f i c a n t for my g l a c i o l o g i c a l a p p l i c a t i o n s . The r e p o r t s of i s o t o p i c data for the I c e f i e l d Ranges (West and Krouse , 1 9 7 2 ; West, unpub l i shed ; Ahern , unpubl i shed [b]) predate t h i s change. Dansgaard ( 1 9 6 9 ) repor ted a r e p r o d u c i b i l i t y of ± 0 . 1 2 ° / o o (per m i l ) for r o u t i n e mass spectrometer measurements of 6 . Ahern (unpubl i shed [ b ] , p . 1 5 8 ) r epor ted ± 0 . 0 3 ° / O o for samples from the S t e e l e G l a c i e r . Dansgaard's l a b o r a t o r y i n Copenhagen a l s o now achieves t h i s l e v e l . T h i s i s adequate for work on g l a c i e r s , s i n c e 6 may vary by s e v e r a l p a r t s ° / 0 0 to tens of p a r t s ° / 0 o for samples from any one g l a c i e r . 1 10 3 .5 .2 FACTORS AFFECTING DEL The nonzero va lue of DEL (6) for an i c e sample from a g l a c i e r i s the r e s u l t of a long s e r i e s of processes i n the h y d r o l o g i c a l c y c l e s ince the water l e f t the w e l l - m i x e d ocean (where 6 i s c l o s e to z e r o ) . At 0 ° C , the vapour pressures of the three major i s o t o p i c forms of water have the approximate r a t i o s ( e . g . Dansgaard, 1964) H O 1 6 : H O 1 8 : HDO = 1.000 : 0.989 : 0.904 2 2 (3 .5 .2 ) and the d i f f e r e n c e s increase w i t h decreas ing temperature . The r e s u l t i n g d i f f e r e n c e s i n v o l a t i l i t y lead to temperature dependent i s o t o p i c f r a c t i o n a t i o n i n evapora t ion and condensat ion proces se s . Under fa s t evapora t ion or condensat ion c o n d i t i o n s (i_.e. e q u i l i b r i u m c o n d i t i o n s do not e x i s t between the vapour and the l i q u i d phases) the f r a c t i o n a t i o n f a c t o r (the r a t i o of the c o n c e n t r a t i o n s of a p a r t i c u l a r i s o t o p i c spec ies i n the two phases) i s c o m p l i c a t e d . The process which tends to c o n t r o l the 6 va lues i n g l a c i e r p r e c i p i t a t i o n i s the condensat ion of d r o p l e t s from c l o u d vapour; f o r t u n a t e l y , t h i s c an , i n most cases , be adequately model led by a R a y l e i g h condensat ion proces s , . i . e . a slow condensat ion ( q u a s i - e q u i l i b r i u m of the vapour and l i q u i d phases) w i t h immediate removal of the condensate (Dansgaard, 1964). For slow condensat ion or e v a p o r a t i o n , the f r a c t i o n a t i o n f a c t o r i s j u s t a r a t i o of the vapour pressures of the d i f f e r e n t i s o t o p i c spec ies of water at the ambient temperature . These r a t i o s are w e l l known above 0°C from l a b o r a t o r y measurements, and have been e x t r a p o l a t e d to -20 °C (Dansgaard, 1964, Table 1) us ing a formula of Zhavoronkov and o t h e r s , (1955). When the 111 R a y l e i g h condensat ion model i s a p p l i c a b l e , the 6 value of the p r e c i p i t a t i o n i s p r i m a r i l y an i n d i c a t i o n of the condensat ion temperature . In genera l terms, 6 va lues tend to decrease w i t h a l t i t u d e and w i t h l a t i t u d e , and, at any one s i t e , to be more negat ive i n winter than i n summer. A c o n t i n e n t a l e f f e c t i s a l s o sometimes observed (Dansgaard, 1964); the 6 va lues of p r e c i p i t a t i o n at constant condensat ion tejnperature may decrease w i t h d i s t a n c e from the ocean due to d e p l e t i o n of heavy i sotopes from the storm systems through p r e c i p i t a t i o n , and due to d i l u t i o n w i t h i s o t o p i c a l l y l i g h t vapour from freshwater sources . Fac to r s other than temperature can i n f l u e n c e the 6 va lue of p r e c i p i t a t i o n . Dansgaard (1964) d i scus sed the e f f e c t s of evapora t ion from f a l l i n g d r o p l e t s , i s o t o p i c exchange between drops and a i r through which they f a l l , n o n - e q u i l i b r i u m phase changes, and v a r i a t i o n s i n the frequency and i s o t o p i c compos i t ion of the source storms. While these processes can cause v a r i a t i o n s i n 6 from storm to s torm, t h e i r e f f e c t on the average summer or winter 6 va lue tends to be constant from year to y e a r . 6 ( 0 1 8 / 0 1 6 ) and 6(D/H) are l i n e a r l y r e l a t e d under R a y l e i g h c o n d i t i o n s (Dansgaard, 1964); s imultaneous measurement of 6 ( 0 1 8 / 0 1 6 ) and 6(D/H) can be used to r e v e a l the presence of n o n - e q u i l i b r i u m condensa t ion . • S e v e r a l processes i n snow and f i r n tend to homogenize the i s o t o p i c d i s t r i b u t i o n , o b l i t e r a t i n g d i f f e r e n c e s between i n d i v i d u a l s torms, an"d sometimes the summer to winter d i f f e r e n c e . In reg ions w i t h summer m e l t i n g , r e c r y s t a l l i z a t i o n i n the presence of p e r c o l a t i n g meltwater can b r i n g the whole v e r t i c a l column of snowpack to the average 6-value (Dansgaard 1 12 and o t h e r s , 1973). However, the e f f e c t s of meltwater are not always s i m p l e ; Ahern (unpubl i shed [ a ] , p . 109) found that p e r c o l a t i n g meltwater c o u l d enhance r a ther than smooth the i s o t o p i c v a r i a t i o n s i n a c o l d snowpack w i t h v a r i a b l e d e n s i t y . In reg ions w i t h no summer m e l t i n g , some smoothing of the i s o t o p i c p r o f i l e occurs due to s u b l i m a t i o n and r e c r y s t a l l i z a t i o n i n the f i r n . V e r t i c a l vapour motion i s most pronounced i n reg ions of h igh v e r t i c a l g rad ien t s of temperature i n the f i r n ( e . £ . due to l a rge seasonal temperature v a r i a t i o n s ) , or in stormy regions w i t h frequent barometr ic pressure changes (Dansgaard, and o t h e r s , 1973). Below the depth at which f i r n reaches a d e n s i t y of 550 kg n r 3 , the " vapour spaces i n the f i r n are i s o l a t e d . The i s o t o p i c p r o f i l e i s smoothed only by s o l i d d i f f u s i o n . Th i s process i s too slow to have an a p p r e c i a b l e e f f e c t on i c e i n the S tee le G l a c i e r . 3 .5 .3 PREVIOUS ISOTOPIC STUDIES Assuming (1) that the p r e c i p i t a t i n g a i r masses f o l l o w s i m i l a r t r a c k s w i t h s i m i l a r frequency the year around and from year to year , (2) that R a y l e i g h condensat ion o c c u r s , (3) that sur face and condensat ion temperatures can be s imply r e l a t e d , (4) tha t the 6-temperature r e l a t i o n s h i p i s constant w i t h t ime , and (5) tha t the flow p a t t e r n of the i c e mass can be c a l c u l a t e d , then 6 va lues i n i c e cores can be r e l a t e d to c l i m a t e at the time of p r e c i p i t a t i o n . T h i s was f i r s t p o i n t e d out by Dansgaard (1954). A thorough d i s c u s s i o n of i c e core s t u d i e s can be found 1 13 i n Chapter 15 of Paterson (1981). The Greenland i c e sheet prov ides the most s u i t a b l e i ce flow and m e t e o r o l o g i c a l c o n d i t i o n s for a s imple c l i m a t i c i n t e r p r e t a t i o n of a deep i c e core (Dansgaard and o t h e r s , 1973). The f i r s t major d r i l l i n g program was undertaken i n 1956 by S . I . P . R . E . (U.S . Army Snow, I c e , and Permafrost Research E s t a b l i s h m e n t , now c a l l e d CRREL, Co ld Regions Research and Eng ineer ing L a b o r a t o r y ) ; a 411 m core was recovered at S i t e 2, i n northwest Greenland . T h i s was fo l l owed by deep cores at Camp Century , Greenland i n 1966 (1387 m), and at Byrd S t a t i o n , A n t a r c t i c a , i n 1968 (2164 m). The Camp Century core has been used to d e r i v e c l i m a t e v a r i a t i o n s over the past 100,000 years (Dansgaard and Johnsen, I969[a ] , I969[b] ; Dansgaard and o t h e r s , 1969, 1971). The Byrd S t a t i o n core a l s o shows long term c l i m a t e v a r i a t i o n s (Eps te in and o t h e r s , 1970; Johnsen and o t h e r s , 1972), but i s more d i f f i c u l t to date a b s o l u t e l y , because the annual v a r i a t i o n s of 6 were not preserved d u r i n g the i c e format ion p roce s s . The l e n g t h of the f l o w l i n e , and the time s c a l e for t h i s hole are i n d i s p u t e (Robin , 1977). Cor ing programs and c l i m a t i c i n t e r p r e t a t i o n s have a l s o been undertaken at other s i t e s i n Greenland (Dansgaard and o t h e r s , 1973) , i n A n t a r c t i c a , i n c l u d i n g Vostok (Barkov and o t h e r s , 1974, 1975, 1977), Dome C ( L o r i u s and o t h e r s , 1979), L i t t l e America V (Dansgaard and o t h e r s , 1977), and Terre A d e l i e ( L o r i u s and M e r l i v a t , 1977), and at s i t e s i n the Canadian A r c t i c , i n c l u d i n g Meighen Ice Cap (Koerner and o t h e r s , 1973; Koerner and P a t e r s o n , 1974) , Devon Ice Cap (Paterson and o t h e r s , 1977; Paterson and C l a r k e , 1978; F i s h e r , 1979), and Agas s i z Ice Cap, E l le smere 1 1 4 I s l a n d (D. F i s h e r , per sona l communicat ion) . Ice cores for i s o t o p i c a n a l y s i s have been obta ined from the p l a t e a u at 5400 m on Mount Logan, Yukon T e r r i t o r y by G. Holdswor th . Dur ing per iods of ex tens ive g l a c i a t i o n , deep sea sediments are enr i ched i n O 1 8 ; the i s o t o p i c compos i t ion of sea water i s a l t e r e d because of the l a rge volume of 0 1 8 - d e p l e t e d i c e on l a n d . Measurements of the i s o t o p i c compos i t ion of deep sea sediments ( e . g . Hayes and o t h e r s , 1976) have complemented the c l i m a t i c s t u d i e s of deep i c e c o r e s . The v a l i d i t y of the c l i m a t e i n t e r p r e t a t i o n of i c e cores has a l s o been supported by other s t u d i e s . Robin (1976) and Johnsen (1977) found that the temperature h i s t o r y d e r i v e d from the i s o t o p i c records was compat ib le w i t h the present v e r t i c a l d i s t r i b u t i o n of temperature i n boreho le s . Paterson and C l a r k e (1978) used the i s o t o p i c a l l y - d e r i v e d temperature h i s t o r y as a boundary c o n d i t i o n for a time-dependent heat flow model for the Devon I s l a n d boreho le s . P i c c i o t t o and o t h e r s , ( i960) demonstrated the ex i s t ence of a l i n e a r r e l a t i o n s h i p between c l o u d temperature and 6 on the coast of East A n t a r c t i c a , and L o r i u s and Delmas (1975) found a l i n e a r r e l a t i o n s h i p between 6 (D/H) and ten metre f i r n temperatures . P i c c i o t t o and o t h e r s , (1968), and L o r i u s and o t h e r s , (1970) used the annual v a r i a t i o n of 6 i n near- sur face samples to measure the recent accumulat ion r a te at A n t a r c t i c s i t e s . L o r i u s and o t h e r s , (1969) measured r e g i o n a l v a r i a t i o n s of 6 (D/H) i n A n t a r c t i c p r e c i p i t a t i o n . West and Krouse (1972) measured i s o t o p i c r a t i o s at s e v e r a l s i t e s i n the S t . E l i a s Mounta ins , Yukon T e r r i t o r y , i n c l u d i n g 115 Mount l o g a n , D i v i d e S t a t i o n , and Rusty G l a c i e r , o b t a i n i n g mass balance es t imates and r e l a t i n g i s o t o p i c compos i t ion to weather p a t t e r n s . A l o n g i t u d i n a l surface sampling of Rusty G l a c i e r (a smal l surge-type g l a c i e r i n the S tee le Creek bas in) showed a sys temat ic increase of 6 w i t h he ight i n the a b l a t i o n zone, demonstrat ing that the genera l i s o t o p i c p a t t e r n of the g l a c i e r was not des t royed by s u r g i n g . Ahern (unpubl i shed [ b ] , p . 162) obta ined an i s o t o p i c p r o f i l e to a depth of 36 m in a borehole at x=15 km (F igure 3.2) on the S tee le G l a c i e r . Th i s p r o f i l e appeared to show p e r i o d i c o s c i l l a t i o n s of wavelength 7 metres and amplitude ± 1 . 5 ° / 0 o i n 6 ( 0 1 8 / 0 1 6 ) . I s o t o p i c s tud ie s on g l a c i e r s and i c e sheets are c l o s e l y r e l a t e d to s t u d i e s of i c e f l o w . The age of the i c e must be determined, and annual l a y e r s cannot always be detected i s o t o p i c a l l y . Time s c a l e s fo r i c e cores have been d e r i v e d by Dansgaard and Johnsen ( I969 [a ] ) , by P h i l b e r t h and Federer (1971), and by Hammer and others (1978), us ing assumptions of steady s t a t e , p r o x i m i t y to an i c e d i v i d e , and s p e c i f i c forms of v e r t i c a l s t r a i n ra te or temperature g r a d i e n t . Paterson and o thers (1977) measured the v e r t i c a l deformation i n the Devon boreho le s ; t h i s a l lowed them to e l i m i n a t e the s t r a i n ra te assumptions when d e r i v i n g a time sca le from the f l o w . For boreholes a t l a rge d i s t a n c e s from i c e d i v i d e s , the h o r i z o n t a l v e l o c i t y i n f l u e n c e s the t ime s ca l e c a l c u l a t i o n , and two-dimens iona l f low s i m u l a t i o n s , such as g iven by Budd and o thers (1971, p . 148) for the Byrd s t a t i o n f l o w l i n e , must be used to date the i c e . 1 16 A conceptua l i n c o n s i s t e n c y can a r i s e w i t h the use of steady s t a t e f low models; the e x i s t e n c e of the c l i m a t i c v a r i a t i o n s , which the i s o t o p i c p r o f i l e s attempt to determine , c o u l d prec lude the e x i s t e n c e of steady s ta te c o n d i t i o n s . In a d d i t i o n , the i c e sur face e l e v a t i o n may vary i n the absence of c l i m a t i c v a r i a t i o n s . For t h i s reason, time-dependent flow models are e s s e n t i a l to the i n t e r p r e t a t i o n of some i s o t o p i c da t a . Jenssen (1978) i n v e s t i g a t e d the e f f e c t of i ce sheet e l e v a t i o n changes on the i s o t o p i c p r o f i l e i n the Vostok c o r e . He c a l c u l a t e d i ce t r a j e c t o r i e s and s imula ted i s o t o p i c p r o f i l e s on the A n t a r c t i c i c e sheet f l o w l i n e from Vostok to M i r n y , assuming p e r i o d i c surges ; the time-dependent i c e sur face p r o f i l e s were computed by Budd and Mclnnes (1978, 1979). T h i s i s the on ly prev ious work of which I am .aware i n which t r a j e c t o r i e s and i s o t o p i c d i s t r i b u t i o n are c a l c u l a t e d i n a t i m e - v a r y i n g i c e mass. Jenssen (1978) was i n t e r e s t e d i n the e f f e c t of surg ing on the c l i m a t i c i n t e r p r e t a t i o n ; i n t h i s c h a p t e r , I i n v e s t i g a t e the f e a s i b i l i t y of u s ing i s o t o p i c i n f o r m a t i o n to r e v e a l the surge h i s t o r y . Jenssen (1978) used a l i n e a r i s o t o p i c - e l e v a t i o n r e l a t i o n s i m i l a r to equa t ion (3 .5 .3 ) below. The d e t a i l s of the f l o w l i n e c a l c u l a t i o n s were not d e s c r i b e d . 3 .5 .4 DEL RELATIONS FOR THE MODEL In the computer model, I assume t h a t : (1) a l l p r e c i p i t a t i o n has the annual average va lue for the g iven e l e v a t i o n or l o c a t i o n , (!•£_• r a p i d i s o t o p i c homogenization i n the f i r n ) , and 1 17 (2) i s o t o p i c d i f f u s i o n i s n e g l i g i b l e i n the s o l i d i c e . I use two 6 - p r e c i p i t a t i o n models, r e p r e s e n t i n g opposi te extremes of topographic c o n t r o l of p r e c i p i t a t i o n . The f i r s t , which I c a l l Model 61 , i s based on the assumption that the i s o t o p i c compos i t ion of s n o w f a l l i s complete ly c o n t r o l l e d by the surface e l e v a t i o n of the g l a c i e r h ( x , t ) , and can be approximated by a l i n e a r r e l a t i o n s h i p over the range of e l e v a t i o n s c o n s i d e r e d , i . e . 6 ( x , t ) = 6 0 + c h ( x , t ) (3 .5 .3) Values of the constants 6 0 and c c o n s i s t e n t w i t h the few measurements i n the S tee le Creek ba s in (West and Krouse , 1972) and on the S tee le G l a c i e r (Ahern, unpubl i shed [b]) are 6 0 = -15.0 °/oo c = -0.004 m- 1 % o (3 .5 .4) The g rad ien t c i s c l o s e to the va lue of -0.005 m" 1 °/oo found by Sharp and o t h e r s , (1960) for c at Blue G l a c i e r , Washington, U . S . A . I t i s w i t h i n a f a c t o r of two of the va lue of -0.002 nr 1 ° / 0 0 r epor ted by Dansgaard (1961). Jenssen (1978) used c = -0.006 irr 1 °/oo i n a m o d e l l i n g study of an A n t a r c t i c f l o w l i n e . I s o t o p i c va lues on i c e sheets and i c e caps can be d e s c r i b e d by Model 61, because the i c e sheet sur face height s t r o n g l y i n f l u e n c e s the v e r t i c a l motion of incoming a i r masses. The second model, c a l l e d Model 62, assumes that the i s o t o p i c compos i t ion of the s n o w f a l l i s complete ly c o n t r o l l e d by the r e g i o n a l mountain topography and by d i s t a n c e from the storm source i n the Gul f of A l a s k a ; 6 i s assumed to be a f u n c t i o n only of p o s i t i o n x . Approximat ing the 1951 S tee le G l a c i e r surface (Topographica l Survey, 1969) by a s t r a i g h t l i n e i n the 1 18 FIGURE 3 .8 . Reference Surface For 6-x F u n c t i o n . The accumulat ion area of the S tee le G l a c i e r i s assumed to extend (see F igure 3.7 (a)) to x=13 km ( v e r t i c a l broken l i n e ) . The curved l i n e i s the observed i c e surface e l e v a t i o n i n 1951. A p p l y i n g the 6 - e l e v a t i o n r e l a t i o n (3 .5 .3 ) to the l i n e a r approximat ion (dashed l i n e ) to t h i s s u r f a c e , and us ing the cons tant s (3 .5 .4 ) g ives a l i n e a r 6-x r e l a t i o n s h i p w i t h the cons tant s i n ( 3 . 5 . 6 ) . accumulat ion a rea , as shown i n F i g u r e 3 .8 , and a p p l y i n g the 6-h r e l a t i o n s h i p to t h i s l i n e , w i t h the constants ( 3 . 5 . 4 ) , g ives the x-6 r e l a t i o n s h i p 6 ( x , t ) = 6 0 + k x (3 .5 .5 ) w i t h the cons tants 6o = — 26.6 ° / o o k = 0.00022 m- 1 ° / o o (3 .5 .6 ) Model 62 i s more a p p r o p r i a t e than Model 61 for sma l l v a l l e y g l a c i e r s . The S tee le G l a c i e r , being a l a rge v a l l e y g l a c i e r , probably f a l l s between the two extremes. Because Model 61 a l l o w s the i s o t o p i c compos i t ion of s n o w f a l l to vary w i t h both p o s i t i o n and t ime , w h i l e Model 62 a l l o w s o n l y v a r i a t i o n w i t h p o s i t i o n , Model 61 can produce l a r g e r v a r i a t i o n s or d i s c o n t i n u i t i e s i n 6 ( 0 1 8 / 0 1 6 ) . Models 61 and 62 119 can be cons idered to g ive the maximum and minimum s t r u c t u r e r e s p e c t i v e l y to the i s o t o p i c d i s t r i b u t i o n w i t h i n the S tee le G l a c i e r . 3.6 MODEL RESULTS: SURGE PERIOD AND TRAJECTORIES 3.6.1 INTRODUCTION The surge p e r i o d i c i t y of the S tee le G l a c i e r i s unknown. In S e c t i o n 3 . 6 . 2 , I show three computer s i m u l a t i o n s of the S tee le G l a c i e r us ing surge per iods spanning the range of 50 to 150 years suggested by f i e l d o b s e r v a t i o n s . I used the same v e l o c i t y p a t t e r n d u r i n g the surge i n a l l three ca ses . The p e r i o d (97 years ) which gave the most reasonable i c e t h i c k n e s s p r o f i l e s at a l l t imes was s e l e c t e d for t r a j e c t o r y and i s o t o p i c c a l c u l a t i o n s . In S e c t i o n 3 . 6 . 3 , I present t y p i c a l t r a j e c t o r i e s for i c e p a r t i c l e s i n t h i s model . These t r a j e c t o r i e s show p e r i o d i c abrupt changes of d i r e c t i o n and speed. 3 .6 .2 PERIODICALLY REPEATING STATE To ob ta in a p e r i o d i c a l l y r epea t ing surge c y c l e , I used the cons tant s i n Table 3.2 to generate the s l i d i n g v e l o c i t y . Th i s s l i d i n g p a t t e r n i s shown at i n t e r v a l s of 0.25 years i n F i g u r e 3 .9 . The surge d u r a t i o n was two y e a r s . I used three surge p e r i o d s ; 47, 97, and 147 years span the range of p o s s i b l e p e r i o d s d e s c r i b e d i n S e c t i o n 3 . 2 . 4 . S t a r t i n g from the n o n s l i d i n g steady s t a te i n F i g u r e 3 .6 , I ran the S tee le G l a c i e r Model 2 through f i f t e e n surge c y c l e s of l e n g t h t 4 . By that t ime , the 120 to (a) 0.0 x . 2 < t 0 ) (km) 8.0 t , (a) 0.75 X . i ( t 0 ) (km) 18.0 t 2 (a) 1.5 X i ( t 0 ) (km) 19.0 t 3 (a) 2.0 X 2 ( t 0 ) (km) 26.0 t « (a) 1 ) 47. 0 2) 97. 0 c . : ( t 0 ) (m a - 1 ) -1000.0 3) 147. 0 c . i ( t 0 ) (m a' 1 ) 7500.0 f 0.0 C i ( t 0 ) (m a' 1 ) 15000.0 U 0 (m a- 1 ) 5000. c 2 ( t 0 ) (m a" 1 ) • 15000.0 TABLE 3 .2 . V e l o c i t y p a t t e r n for S tee le surge. The three values of surge p e r i o d t « were used for F igures 3 .10, 3 .11 , and 3.12 r e s p e c t i v e l y . model no longer "remembered" the i n i t i a l steady s t a t e c o n d i t i o n ; i t repeated the same sur face p r o f i l e sequence to one par t i n 10' w i t h each new c y c l e for 97 and 147 year p e r i o d s , and to a few par t s i n 10 3 "with a 47 year p e r i o d . The time steps i n the computer model must be very smal l when the g l a c i e r i s surg ing i n order to mainta in accuracy (see Sec t ion A 1 . 5 ) . Table 3.3 summarizes the numerica l time step sequence used. The other numer ica l and p h y s i c a l cons tant s had the va lues shown i n Table 3 . 1 . F i g u r e s 3.10 through 3.12 show the pre-surge g l a c i e r sur face ( s o l i d l i n e ) at t 0 , and the post-surge sur face (broken l i n e ) at t 3 for the three surge p e r i o d s of 47, 97, and 147 years r e s p e c t i v e l y . The model w i t h a 47 year surge p e r i o d ( f i g u r e 3.10) has a pre-surge p r o f i l e ( s o l i d l i n e ) which i s l e s s than 100 m t h i c k beyond x=20 km. The post-surge p r o f i l e i s l e s s than 100 m t h i c k beyond x=25 km. I t seems u n l i k e l y that a surge of i c e always l e s s than 100 m t h i c k c o u l d advance 12 km i n two y e a r s . The p r o f i l e s shown i n t h i s diagram are u n r e a l i s t i c . Th i s c o n c l u s i o n i n d i c a t e s that the S tee le G l a c i e r accumulat ion area cannot p r o v i d e s u f f i c i e n t mass i n j u s t 47 years to generate 121 8000 6000 o "54000 > cn c CO 2000 T Steele Glacier Model 2 Surge period 97 year* Sliding velocity at intervals of 0.25 years t=1.0 0 10 20 30 40 x (km) FIGURE 3 .9 . S l i d i n g V e l o c i t y : S tee le G l a c i e r Mode l . The s l i d i n g v e l o c i t y p r o f i l e i s shown at i n t e r v a l s of 0.25 years d u r i n g the surge of two years d u r a t i o n . The end of each curve i n d i c a t e s the p o s i t i o n of the advancing terminus ( for the model w i t h 97 year surge p e r i o d . ) surges which move as q u i c k l y or as f a r as the 1966-67 example. E i t h e r some surges must be l e s s v i g o r o u s , or the surge p e r i o d must be s u b s t a n t i a l l y longer than 47 y e a r s . The model w i t h a p e r i o d of 147 years (F igure 3.12) i s reasonably t h i c k at a l l t i m e s . I t c o u l d be c r i t i c i z e d on the ba s i s of the e x c e p t i o n a l l y t h i c k i c e lobe below x=25 km f o l l o w i n g the surge . Radio echo sounding ( e . £ . Narod and C l a r k e , 1980) of the S tee le G l a c i e r would r e s o l v e the v a l i d i t y of t h i s c r i t i c i s m . T h i s model c o u l d be acceptab le for the t r a j e c t o r y a n a l y s i s , but I d i d not use i t because the obse rva t ions of S t an ley (1969), Wood (1972) and 1 22 40001 | 1 | 1 1 Steele Glacier Model 8 Surge period 47 years Pre—surge profile Post-surge profile 3000 J - 4 - » c n 2000 1000 I I I I I J 0 10 20 30 40 50 x (km) FIGURE 3 .10 . P r e - And Post-surge P r o f i l e s : 47 Year P e r i o d . The s o l i d curve shows the g l a c i e r as the surge beg ins . The broken curve shows the sur face e l e v a t i o n as the surge ends two years l a t e r . The s l i d i n g f u n c t i o n i s g iven i n Table 3.2 and F igure 3 .9 . The i c e depths of the lower g l a c i e r are unreasonably t h i n ; the S tee le has i n s u f f i c i e n t accumulat ion to surge as v i g o r o u s l y as the 1966-67 event as f r e q u e n t l y as every 47 y e a r s . Sharp (1947), d e s c r i b e d i n S e c t i o n 3 . 2 . 4 , i n d i c a t e that 150 years i s an upper l i m i t for the time between surges . The surge p e r i o d of 97 years i s the most acceptab le of the t h r e e . There i s l i t t l e change i n the i c e surface e l e v a t i o n above x=8 km. Between 8 km and 20 km, the surface drops by up to 110 m, w h i l e below 20 km, the surface r i s e s d u r i n g the surge by up to 120 m. T h i s p a t t e r n agrees w e l l w i t h the zones d e s c r i b e d by S t an ley (1969). I used t h i s model for the t r a j e c t o r y and 123 4000 3000 CD X 2000 1000 T T T Steele Glacier Model 2 Surge period 97 years Pre-surge profile Post—surge profile - - 0 10 20 (km) 30 40 50 FIGURE 3 .11 . P r e - And Post-surge P r o f i l e s : 97 Year P e r i o d . The s o l i d curve shows the g l a c i e r as the surge beg ins . The broken curve shows the surface e l e v a t i o n as the surge ends two years l a t e r . The s l i d i n g f u n c t i o n i s g iven i n Table 3.2 and F i g u r e 3 .9 . T h i s model i s used i n S e c t i o n 3 .7 . i s o t o p i c c a l c u l a t i o n s . F i g u r e 3.13 shows an or thographic view of the i c e t h i c k n e s s h ( x , t ) d u r i n g one complete surge c y c l e for the case w i t h a 97 year p e r i o d . The r a p i d decrease i n surface e l e v a t i o n between x=8 km and x=20 km and between time zero and two years i s hidden by the pre-surge p r o f i l e , but the r a p i d surface r i s e of the lower g l a c i e r can be seen. The subsequent slow a b l a t i o n of the terminus r e g i o n , and the b u i l d u p and advance of the i c e i n the source reg ion are e v i d e n t . Three steps or waves can be seen i n the a b l a t i o n zone. One wave forms d u r i n g each surge . The i c e 124 4000 3000 J C O l X 2000 1000 (km) FIGURE 3 .12 . P r e - And Post-surge P r o f i l e s : 147 Year P e r i o d . The s o l i d curve shows the g l a c i e r as the surge beg ins . The broken curve shows the surface e l e v a t i o n as the surge ends two years l a t e r . The s l i d i n g f u n c t i o n i s g iven i n Table 3.2 and F i g u r e 3 . 9 . T h i s model i s p h y s i c a l l y reasonable , but 147 years i s an upper l i m i t on the surge p e r i o d based on the obse rva t ions of S tan ley (1969) , Wood (1972) , and Sharp (1947) . I n t e r v a l Time (a) At (a) t 0 - t , 0.0 - 0.75 0.01 1 t i - t 2 0.75 - 1.5 0.01 t 2 - t 3 1.5 - 2.0 0.01 t 3 ~ t a 2 .0 - 7.0 0.1 7.0 - 97.0 1.0 TABLE 3 . 3 . Numerica l time steps for S t e e l e surge . 125 FIGURE 3 .13 . S tee le G l a c i e r T h i c k n e s s : One Surge C y c l e . Ice t h i c k n e s s h ( x , t ) i n o r thographic view from a p o i n t 45° above the x a x i s , and 45° around the t h i c k n e s s a x i s from the time a x i s . The model used the constants i n Tables 3 . 1 , 3 . 2 , and 3 .3 . The surge s t a r t e d at time=0.0 and ended at 2 y e a r s . The surge p e r i o d was 97 y e a r s . The t ransver se l i n e s i n d i c a t e 1 km i n t e r v a l s , and the l o n g i t u d i n a l l i n e s are spaced at i n t e r v a l s of 5 y e a r s . t h i c k e n s as a r e s u l t of being forced i n t o the converg ing channel (see F i g u r e 3.7 ( b ) ) . The observed p r o f i l e of the S tee le G l a c i e r i n 1951 (F igure 3.3 (c)) has long surface u n d u l a t i o n s , but these c o u l d r e s u l t from bed topography. Repeated surface mapping and a r ad io-echo depth survey would re so lve t h i s q u e s t i o n . The u n d u l a t i o n s i n F igure 3.13 may be an a r t i f i c e of the model due to the poorly-known s p a t i a l v a r i a t i o n of the surge v e l o c i t y , and the i n c l u s i o n of the Hodgson G l a c i e r through the wid th f u n c t i o n 126 1 1 1 1 Steele Glacier Model 2 x (km) FIGURE 3.14. Ice T r a j e c t o r i e s .For 97 Year Surge P e r i o d . T r a j e c t o r i e s for 5 p a r t i c l e s s t a r t i n g on the i c e surface as a surge beg ins . The dashed p r o f i l e i s the i ce t h i c k n e s s at t 0 . The dashed-dotted p r o f i l e i s the post-surge i c e t h i c k n e s s at t 3=2 y e a r s . The arrows on the t r a j e c t o r i e s i n d i c a t e time i n t e r v a l s of 48.5 year s , ! • £ . the midpoint and end of each 97 year surge c y c l e . Beyond x=22 km, the i c e i s v i r t u a l l y stagnant between surges , and the arrows for mid-per iod and fo r the end of the p e r i o d are superimposed. shown i n F i g u r e 3.3 ( b ) . 3 .6 .3 ICE TRAJECTORIES F i g u r e 3.14 shows the t r a j e c t o r i e s of p a r t i c l e s s t a r t i n g on the i c e sur face at t 0 as a surge b e g i n s . The arrows on the t r a j e c t o r i e s show the p a r t i c l e l o c a t i o n s each 48.5 year s , i.e. at the midpoint and at the t e r m i n a t i o n of each 97 year surge c y c l e . The pre-surge (dot ted curve) and the post-surge 1 27 (dotted-dashed curve) i c e t h i c k n e s s are a l s o shown. The i c e beyond x=22 km i s v i r t u a l l y stagnant between surges , and the arrows i n d i c a t i n g the end of each surge c y c l e are superimposed on the m i d - c y c l e arrow. Between x=8 km and x=20 km, the t r a j e c t o r i e s descend a b r u p t l y d u r i n g the surge; the h igh l o n g i t u d i n a l s t r a i n ra te r e s u l t i n g from the grad ient of the s l i d i n g v e l o c i t y i n F igure 3.9 causes a genera l t h i n n i n g i n t h i s zone. The t r a j e c t o r i e s r i s e d u r i n g the surge i n the reg ion below x=20 km because the v e l o c i t y g rad ien t i s sma l l and the i c e t h i c k e n s as i t f lows i n t o a channel of decreas ing w i d t h . The endpoint of each t r a j e c t o r y i n d i c a t e s the p o s i t i o n of i c e sur face at the time the i c e p a r t i c l e came to the s u r f a c e . 3.7 MODEL RESULTS: DISTRIBUTION OF ISOTOPES 3.7.1 INTRODUCTION In t h i s s e c t i o n , I show the i s o t o p i c d i s t r i b u t i o n s found by c a l c u l a t i n g i ce t r a j e c t o r i e s for the surge model w i t h 97 year p e r i o d (F igures 3.11 and 3.13) and the p r e c i p i t a t i o n - 6 models 61 and 62. I have chosen to look at the i so tope d i s t r i b u t i o n at four times dur ing the c y c l e : (1) when the surge begins at t=0.0 , (2) a t mid-surge , t=1.0 y e a r s , (3) at t=1u.O y e a r s , when i t i s aga in p o s s i b l e to walk over most of the g l a c i e r s u r f a c e , and (4) at t=50.0 y e a r s , near the midpoint of the qu iescent p e r i o d . S ince the v e l o c i t y f i e l d i s smooth ly-vary ing even d u r i n g a surge ( S t a n l e y , 1969), d i s c o n t i n u i t i e s cannot be in t roduced i n t o the i so tope d i s t r i b u t i o n w i t h i n the i c e by the r a p i d flow a l o n e ; 1 28 t h i s would r e q u i r e shear f r a c t u r i n g and d i s l o c a t i o n on a l a rge s c a l e . D i s c o n t i n u i t i e s or sharp changes i n the g rad ien t of 6 can be in t roduced only by abrupt changes i n the 6-value of p r e c i p i t a t i o n f a l l i n g onto a g iven element of i c e at the g l a c i e r s u r f a c e . When 6 i s cons idered to be a f u n c t i o n of p o s i t i o n and sur face e l e v a t i o n (Sec t ion 3 . 5 . 2 ) , the 6-value of new p r e c i p i t a t i o n f a l l i n g on an element of i ce at the surface can change a b r u p t l y only i f (1) the surface e l e v a t i o n changes a b r u p t l y , or (2) the sur face element i s r a p i d l y moved to a new p o s i t i o n at which s n o w f a l l has a d i f f e r e n t 6 - v a l u e . Thus, i s o t o p i c d i s c o n t i n u i t i e s can be formed only i n i c e which (1) i s at the i c e - a i r i n t e r f a c e at the moment a surge beg ins , (2) i s i n the accumulat ion r e g i o n , and (3) p a r t i c i p a t e s i n the surge through forward motion or a decrease i n i c e t h i c k n e s s . The reg ion above x=8 km on the S tee le G l a c i e r d i d not take part i n the 1966-67 surge ( S t a n l e y , 1969); t h i s fea ture i s i n c o r p o r a t e d i n t o the computer model. S ince a l l the i c e downstream from x=25, and the deep i c e throughout the g l a c i e r comes from the reg ion upstream from x=8 km, there cannot be any d i s c o n t i n u i t i e s or abrupt v a r i a t i o n s i n 6 i n t h i s i c e . (The surges do, of course , a l t e r the p o s i t i o n s and d i s t o r t the g rad ien t s of the i s o d e l l i n e s , but a l l q u a n t i t i e s remain s l o w l y - v a r y i n g ; t h i s would be a d i f f i c u l t and s u b t l e matter to i n t e r p r e t i n a f low regime as c o m p l i c a t e d as that of a surg ing v a l l e y g l a c i e r ) . In S e c t i o n 3 . 7 . 2 , I show l o n g i t u d i n a l v e r t i c a l s e c t i o n s of the S t e e l e G l a c i e r between x=5 km and x=25 km at the four times metioned above us ing the height-dependent p r e c i p i t a t i o n model 61 . I present d e t a i l e d sur face p r o f i l e s of 6 at the same 129 four t imes i n Sec t ion 3 . 7 . 3 . In S e c t i o n 3 . 7 . 4 , I show s e c t i o n s and p r o f i l e s at the same times us ing the the p r e c i p i t a t i o n model 62. 3 .7 .2 MODEL 61: LONGITUDINAL 6 SECTIONS Would a w e l l - p l a c e d b o r e h o l e , or s e r i e s of boreholes i n the S tee le G l a c i e r recover c o n v i n c i n g i sotop ' ic evidence of past surges? In F igure 3 .15 , I show the expected i s o t o p i c d i s t r i b u t i o n s i n a l o n g i t u d i n a l s e c t i o n at the four t imes 0, 1, 10, and 50 years measured from the i n i t i a t i o n of a surge . The reg ion from x=5 km to x=25 km c o n t a i n s a l l the i ce which i s capable of c o n t a i n i n g d i s c o n t i n u i t i e s i n 6 and i t s g rad ient (see S e c t i o n 3 . 7 . 1 ) . A 6-value was as s igned to each p o i n t on a mesh ( i n d i c a t e d by dots) w i t h a h o r i z o n t a l spacing of 500 m, and a v e r t i c a l spacing of 5 m. The 6-value was c a l c u l a t e d by a p p l y i n g the i s o t o p i c model 61 to the s t a r t i n g coord ina te s of the i c e O v e r l e a f : FIGURE 3 .15 . L o n g i t u d i n a l 6 s e c t i o n s : Model 61. The 6-values w i t h i n the S tee le G l a c i e r are contoured at i n t e r v a l s of 0 . 1 ° / o o « T n e dots show the 5x500 m mesh at which 6 was eva lua ted us ing Model 61 . (Sec t ion 3 . 5 . 4 ) . The surge model (F igures 3.11 and 3.13) had a p e r i o d of 97 y e a r s . Each l i n e of s teep 6 grad ient which i n t e r s e c t s the i c e surface and i s a t tenuated w i t h d i s t a n c e upstream and i n t o the i c e o u t l i n e s the r e l i c t i c e - a i r i n t e r f a c e at the s t a r t of a p rev ious surge . a) t=0 y e a r ; pre-surge p r o f i l e . b) t=1 y e a r ; mid-surge . c) t= l0 y e a r s ; post-surge p r o f i l e . Foot t r a v e l would again be p o s s i b l e . d) t=50 y e a r s ; t h i s i s approx imate ly the midpoint of the qu ie scent p e r i o d . In Model 61, p r e c i p i t a t i o n i s a l i n e a r f u n c t i o n of i c e sur face e l e v a t i o n ; t h i s tends to over-es t imate the magnitude of i s o t o p i c d i s c o n t i n u i t i e s due to s u r g i n g . 130 -C cn "5 100 10 15 x (km) 131 t r a j e c t o r i e s which passed through the meshpoints at the time of the c r o s s - s e c t i o n . The Model 61 ( Sec t ion 3.5.4) assumes that the i s o t o p i c r a t i o of p r e c i p i t a t i o n i s a l i n e a r f u n c t i o n of the i c e surface e l e v a t i o n . I c a l c u l a t e d the t r a j e c t o r i e s us ing the v e l o c i t y eva luated on a mesh (F igure 2.8) w i t h DZ = 15 m Ax = 500 m (3.7.1) (The v e l o c i t y mesh can be coarser than the i s o t o p i c mesh because there are no d i s c o n t i n u i t i e s i n the v e l o c i t y f i e l d ) . The 6-va lues for the i s o t o p i c mesh were contoured a u t o m a t i c a l l y w i t h an i n t e r v a l of 0 . 1 ° / o o « The only fea tures i n the otherwise smooth i s o t o p i c d i s t r i b u t i o n s are curved l i n e s which i n t e r s e c t the g l a c i e r sur face and d i p upstream i n t o the i c e . C r o s s i n g one of these l i n e s from above to below, 6 decreases by an amount which v a r i e s from 0 . 8 ° / o o ° r more at the surface near the f i r n l i n e (X=13 km) to near zero at depth . Each of these l i n e s r e v e a l s the present l o c a t i o n of i ce which was at the surface i n the accumulat ion reg ion at the time a p rev ious surge began. The l i n e generated by the surge c y c l e shown i n F i g u r e 3.15 i s not v i s i b l e unti& (d) at 50 y e a r s . From zero to ten year s , the d i s c o n t i n u i t y i s too c l o s e to the i c e surface and to the edge of the 5 m mesh to be d i s p l a y e d i n the contour p l o t s . I d i s c u s s the i c e surface f u r t h e r i n the next s e c t i o n . To understand the c r e a t i o n and e v o l u t i o n of one of these r e l i c t pre-surge sur face l i n e s , cons ider the i c e at x=12 km where the - 2 4 ° / 0 o contour i n t e r s e c t s the i c e sur face i n F i g u r e 3.15 ( a ) . Dur ing the surge, t h i s element of i c e moves r a p i d l y downstream and to a lower e l e v a t i o n . In ( b ) , i t i s at 1 32 x=8 km. The combined e f f e c t of bed s lope and i c e t h i n n i n g has lowered t h i s i c e element by approximate ly 100 m, so , u s ing the value of c from ( 3 . 5 . 4 ) , the snow f a l l i n g on t h i s element of i ce has a 6-value of -23 .6°/oo» A t the end of the surge , t h i s i ce element i s f u r t h e r downstream and s t i l l l ower ; the 6-value of new snow on t h i s element i s near -23°/oo. A 6 d i s c o n t i n u i t y of the order of 1°/oo has been c r e a t e d . As snow accumulates dur ing the qu iescent phase, the d i s c o n t i n u i t y i s b u r i e d ; i t i s at a depth of 20 m i n ( d ) . The magnitude of the 6 t r a n s i t i o n across the r e l i c t surface decreases w i t h d i s t a n c e up the g l a c i e r , because the s l i d i n g v e l o c i t y d u r i n g the surge and the amount of s u r f a c e - l o w e r i n g a l s o decrease i n t h i s d i r e c t i o n (see F i g u r e s 3.9 and 3 .11 ) . A f t e r one complete surge , t h i s l i n e i s c l e a r l y v i s i b l e i n F i g u r e 3.15 ( a ) ; i t i n t e r s e c t s the i c e surface at X=15 km. During the subsequent surge (b) and ( c ) , t h i s l i n e i s seen to move down the g l a c i e r . S ince i t then i n t e r s e c t s the g l a c i e r surface i n the a b l a t i o n zone (at x=20 km i n ( c ) ) , the downstream p a r t , w i t h the l a r g e s t 6 c o n t r a s t , i s r a p i d l y des troyed by m e l t i n g (d ) , u n t i l , at the s t a r t of the t h i r d surge c y c l e ( a ) , a l l tha t remains i s a s l i g h t p e r t u r b a t i o n of the i s o d e l l i n e s near x=17.5 km. An i s o t o p i c p r o f i l e i n a s i n g l e borehole would show, at be s t , a s i n g l e abrupt decrease of about 0 . 8°/ Oo to 0 . 5 ° / o o at the most recent or second most recent r e l i c t pre-surge s u r f a c e . Because t h i s i c e i s d i s p l a c e d by three to f i v e km d u r i n g each surge , i t i s not p o s s i b l e to see more than one d i s c o n t i n u i t y i n a s i n g l e b o r e h o l e . A l a rge number of boreholes (5-30) would be r e q u i r e d to i n t e r p r e t the 6 p a t t e r n w i t h r e l i a b i l i t y and 1 33 p r e c i s i o n . T h i s would be an expensive f i e l d p r o j e c t . 3 .7 .3 MODEL 61: SURFACE 6 PROFILES A cheap a l t e r n a t i v e to borehole d r i l l i n g i s d e t a i l e d sur face sampling i n the a b l a t i o n zone. In F i g u r e 3.16, I show the i s o t o p i c 6-values a long the g l a c i e r sur face for the same models and the same times as for the c r o s s - s e c t i o n s i n S e c t i o n 3 . 7 . 2 . The 6-values were p o o r l y determined by the computer model where the p r o f i l e s are broken . When an i c e p a r t i c l e i s very near the s u r f a c e , and i t s v e l o c i t y i s n e a r l y p a r a l l e l to the s u r f a c e , smal l t r u n c a t i o n e r r o r s i n the normal v e l o c i t y component (see S e c t i o n 2 .5 .3 ) can cause the t r a j e c t o r y to i n t e r s e c t the i c e surface at a l a rge d i s t a n c e from the c o r r e c t p o s i t i o n ; t h i s g ives an i n c o r r e c t va lue for 6. I e s t imate the e r r o r i n 6 to be ± 0 . 1 ° / O o i n t h i s r e g i o n . E l sewhere , the t r a j e c t o r i e s i n t e r s e c t the i c e surface at l a r g e r ang le s , g i v i n g r e l i a b l e i n t e r s e c t i o n p o s i t i o n s and 6-va lues . At any g iven t ime , i s o t o p i c d i s c o n t i n u i t i e s due to three p rev ious surges can be seen, d i m i n i s h i n g i n ampli tude from 0 . 8 ° / o o to 0 . 2 ° / o o w i t h age. As i n the prev ious s e c t i o n , i t i s p o s s i b l e to f o l l o w the motion of any r e l i c t pre-surge s u r f a c e , t h i s t ime see ing on ly where i t i n t e r s e c t s ' the present i c e s u r f a c e . The four i s o t o p i c peaks and d i s c o n t i n u i t i e s l a b e l l e d t=0, t=1, t= lO, and t=50 i n F i g u r e 3.16 show the p o s i t i o n and the ampl i tude change across a r e l i c t sur face d u r i n g the c y c l e i n which i t forms. The d i s c o n t i n u i t i e s l a b e l l e d s imply 1, 10, and 50 show how the p a t t e r n moves downstream w i t h the i c e i n the 134 -23 to O \ G O -27 t=50 Steele Glacier Model 2 Surface isotopic profiles Precipitation Model <J 1 Surge ' Surge eriod 97 years uration 2 years 10 15 20 25 (km) 30 FIGURE 3.16. Model 61: Surface 6 P r o f i l e s . 6 was c a l c u l a t e d at 50 m i n t e r v a l s a long the surface of the S tee le G l a c i e r model w i t h 97 year surge p e r i o d (F igures 3.11 and 3 .13 ) . The p r o f i l e s correspond to the l o n g i t u d i n a l s e c t i o n s i n F igure 3 .15 . The d i s c o n t i n u i t i e s i n 6 occur where the r e l i c t pre-surge surfaces i n F i g u r e 3.15 i n t e r s e c t the i c e s u r f a c e . The 6-values may be i n e r r o r by ± 0 . 1 ° / O o where the l i n e s are broken, due to sma l l v e r t i c a l v e l o c i t y e r r o r s i n t r a j e c t o r i e s n e a r l y p a r a l l e l to the s u r f a c e . The f i r n l i n e i s at x=13 km. second surge (1 and 10), then how the i n t e r s e c t i o n p o i n t r e t r e a t s back upstream (50) as the sha l low downstream p o r t i o n of the r e l i c t sur face i s melted o f f . In format ion on e i t h e r surge p e r i o d i c i t y , or on i c e displacement d u r i n g prev ious surges c o u l d be d e r i v e d from a sur face i s o t o p i c sampling program i f Model 61 i s a p p r o p r i a t e for 1 35 the S tee le G l a c i e r , and i f the background