UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Accurate modelling of glacier flow Waddington, Edwin Donald 1981

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

Item Metadata

Download

Media
[if-you-see-this-DO-NOT-CLICK]
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
Original Record: 1.0052458 +original-record.json
Full Text
1.0052458.txt
Citation
1.0052458.ris

Full Text

ACCURATE MODELLING OF GLACIER FLOW by EDWIN DONALD WADDINGTON B.Sc,  U n i v e r s i t y of T o r o n t o ,  1971  M.Sc,  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 in THE FACULTY OF GRADUATE STUDIES (Department of G e o p h y s i c s and  We a c c e p t t h i s t h e s i s as to the  Astronomy)  conforming  required standard  THE UNIVERSITY OF BRITISH COLUMBIA November,  1981  © Edwin Donald Waddington,  1981  In p r e s e n t i n g  this thesis in p a r t i a l  f u l f i l m e n t of  requirements f o r an advanced degree at the  the  University  of B r i t i s h Columbia, I agree t h a t  the L i b r a r y s h a l l make  it  and  f r e e l y a v a i l a b l e f o r reference  study.  I  further  agree t h a t p e r m i s s i o n f o r e x t e n s i v e copying of t h i s t h e s i s f o r s c h o l a r l y purposes may  be  department or by h i s or her understood t h a t  granted by  representatives.  s h a l l not  be  Geophysics and  The U n i v e r s i t y of B r i t i s h 2075 Wesbrook P l a c e Vancouver, Canada V6T 1W5  /7cn  t  e  It is  allowed without my  permission.  Da  my  copying or p u b l i c a t i o n of t h i s t h e s i s  f o r f i n a n c i a l gain  Department of  the head of  February 20,  1982,  Astronomy  Columbia  written  i i ABSTRACT  Recent  interest  i n c l i m a t i c change and i c e .sheet  p o i n t s out the need f o r a c c u r a t e and n u m e r i c a l l y of  time-dependent  to  believe  that  much  of  nonlinear  incorrect.  instability  purposes  of  this  accepted  glaciologically As  in  thesis  not  been  are  and  there  to  and  ice  good  is  the importance of  the  recognized.  The  d e v e l o p and t o v e r i f y a new  f l o w , compare the model  interesting  is  sheets  widely  model, and to demonstrate  earlier  continuity  In p a r t i c u l a r ,  has  n u m e r i c a l model f o r g l a c i e r widely  models  the p u b l i s h e d l i t e r a t u r e on  n u m e r i c a l m o d e l l i n g of the flow of g l a c i e r s quantitatively  stable  i c e masses. L i t t l e a t t e n t i o n has been p a i d t o  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, reason  variations  to  another  the model i n  several  applications.  work,  the  equation together  computer  model  w i t h a flow law f o r  solves ice.  the  Thickness  p r o f i l e s a l o n g flow l i n e s a r e o b t a i n e d as a f u n c t i o n of time f o r a temperate i c e mass w i t h balance.  A  set  of  necessary  n u m e r i c a l model of g l a c i e r solutions  arbitrary  bed  tests  flow  is  to  topography be  and  satisfied  presented.  The  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 ;  mass by any  numerical  these include  a s i m p l e 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 t o check t e r m i n u s m o b i l i t y , and B u r g e r s *  e q u a t i o n t o check c o n t i n u i t y and dynamic  behaviour  with f u l l n o n l i n e a r i t y . An  attempt  has  been  computer model of Budd Mclnnes ( u n p u b l i s h e d ) . numerical  made  and  t o v e r i f y the a c c u r a c y of (1975)  and  These a u t h o r s have r e p o r t e d problems  with  instability.  Mclnnes  If  the  (1974),  existing  Rudd  the  documentation  is  iii  accurate,  the Budd-Mclnnes model appears  conservation The  violations  new  each  time  Glen's  step;  the  model  this  suffer  from  mass  globally.  developed  the v e l o c i t y f i e l d  f l o w law f o r  yields  both l o c a l l y and  numerical  used t o r e c o n s t r u c t  to  i n t h i s t h e s i s can be  within  the  glacier  at  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  ice.  Integration  trajectories  of  of  this  individual  velocity  ice  field  elements flowing  through the t i m e - v a r y i n g i c e 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 a  steady  state  ice  sheet  t h e s i s i s not r e s t r i c t e d violations  of  mass  w i t h an a n a l y t i c a l (Nagata,  t o steady  and  it  avoids  and the a p p r o x i m a t i o n s  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 The'feasibility prehistoric  of u s i n g  surging  stable  isotopes  and  surge  gave  surge c y c l e . By plausible height,  duration  were  the most r e a l i s t i c calculating  relationships  to  specified,  ice  between  surge  investigate  for  times  D i s c o n t i n u i t i e s of up surfaces  dipping  before, to  upstream  0.8°/Oo into  trajectories  during, were the  difficult  interface  at the time a p r e v i o u s  ice  roughly the  using  two  position  profiles  of  after  glacier.  the  throughout  and  and  on  of  and  found  s u r f a c e s i s the p r e s e n t l o c a t i o n of the ice-air  period  6(01B/016)  sliding  based  ice thickness  l o n g i t u d i n a l s e c t i o n s and s u r f a c e  constructed  about  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  o b s e r v a t i o n s of the 1966-67 s u r g e . A 100 y e a r s  the  models.  model s i m u l a t i n g the S t e e l e G l a c i e r , Yukon T e r r i t o r y . A < velocity  for  1977). The model i n t h i s  state,  conservation,  solution  6 a  which  were surge.  across Each  or  several of  these  formed  the  surge began. I t may be  t o observe t h e s e s u r f a c e s on the S t e e l e G l a c i e r due  to  i v  the l a r g e and p o o r l y - u n d e r s t o o d background v a r i a b i l i t y of The g e n e r a t i o n of wave o g i v e s has been the  theory  of  Nye  (I958[b])r  wherein  c o m b i n a t i o n of s e a s o n a l v a r i a t i o n i n mass d e f o r m a t i o n i n an i c e f a l l . is  shown  in  this  thesis  is  the  impulse  following  waves are caused by a balance  The wave t r a i n g e n e r a t e d  and  on a g l a c i e r  mass  balance  function.  response of the g l a c i e r  w i d t h and mass b a l a n c e a l s o c o n t r i b u t e t o the wave  wave  This  surface to a  s t e p 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  formulation  plastic  t o 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 integral  examined  6.  glacier  train.  This  i s used t o e x p l a i n why many i c e f a l l s do not g e n e r a t e  ogives  large p l a s t i c  in  s p i t e of l a r g e s e a s o n a l b a l a n c e v a r i a t i o n s and  deformations.  V  TABLE OF CONTENTS  Abstract  i i  List  Of T a b l e s  xiv  List  Of F i g u r e s  xv  Acknowledgements Chapter  xix  1: B e g i n n i n g s  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 C o n v e n t i o n s U s e d  9  1.2 P r e v i o u s  Work  11  1.2.1 F o u r C e n t u r i e s  Of G l a c i e r Flow T h e o r y  1.2.2 P r e v i o u s  Ice P r o f i l e  1.2.3 P r e v i o u s  I c e T r a j e c t o r y Models  1.3 The C o n t i n u i t y  Equation  Models  11 13 15  F o r An I c e Mass  17  1.3.1 The G e n e r a l G l a c i e r F l o w P r o b l e m  17  1.3.2 R e c t a n g u l a r  18  Cross-section  1.3.3 The C o n t i n u i t y  Flow Model  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 1.4 P h y s i c s 1.4.1  22  Of I c e D e f o r m a t i o n  23  Introduction  23  1.4.2 S t r e s s E q u i l i b r i u m E q u a t i o n s  24  1.4.3 C o n s t i t u t i v e R e l a t i o n F o r I c e D e f o r m a t i o n  26  1.4.4 S h e a r S t r e s s  30  1.5 B a s a l 1.5.1  And I c e F l u x  Sliding  35  Introduction  1.5.2 B a s a l  35  I c e T e m p e r a t u r e And S l i d i n g  1.5.3 P h i l o s o p h y  Of S l i d i n g  In T h i s  Study  35 36  vi Chapter 2.1  2: Models And T e s t s  ;  Introduction  2.1.1  Outline  2.1.2  Importance  38 38 38  Of Model T e s t i n g  38  2.2 C o n t i n u i t y E q u a t i o n P r o f i l e Model  42  2.2.1  Introduction  42  2.2.2  The N u m e r i c a l Scheme  43  2.2.3  Boundary C o n d i t i o n s  46  2.2.4  Numerical 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  Introduction  57  2.3.2  C o n t i n u i t y Test With Terminus M o t i o n  57  2.3.3  C o n t i n u i t y Test With B u r g e r s '  62  2.4 The Ice T r a j e c t o r y  Equation  Computer Model . .'  67  2.4.1  Introduction  67  2.4.2  The V e l o c i t y And Displacement  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  Accu racy Of The T r a j e c t o r y  70  Fields  Model  2.5 T e s t i n g The T r a j e c t o r y Model  67  74  2.5.1  Introduction  74  2.5.2  Nagata Ice Sheet Test  74  2.5.3  S u r f a c e Mass C o n s e r v a t i o n Test  79  Chapter  3: Can S t a b l e I s o t o p e s R e v e a l A H i s t o r y Of S u r g i n g ? 83  3.1  Introduction  3.2 S t e e l e G l a c i e r 3.2.1  General D e s c r i p t i o n  83 85 85  vii  3.2.2  G l a c i e r Surges  3.2.3  O b s e r v a t i o n s Of  3.2.4  Period  3.3  Of  86 The  Steele  N u m e r i c a l Model  Steele Glacier  88  Surges  91  1  93  3.3.1  F l o w Law  3.3.2  Bed  3.3.3  Channel Width  96  3.3.4  Mass B a l a n c e  97  3.3.5  Cyclic  3.4  C o n s t a n t s And  Shape F a c t o r  93  Topography  94  Surge P a t t e r n  For  The  Model  100  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  3.4.2  Simplifications  3.5  Stable  Isotopes  D e f i n i t i o n Of  3.5.2  Factors  3.5.3  Previous  3.5.4  Del  108  The  Del  108  Scale  108  A f f e c t i n g Del Isotopic  Relations  Model R e s u l t s :  For  110  Studies  112  The  116  Model  Surge P e r i o d  3.6.1  Introduction  3.6.2  Periodically  3.6.3  Ice  3.7  103  In G l a c i o l o g y  3.5.1  3.6  1  And  Trajectories  119 119  Repeating  State  119  Trajectories  Model R e s u l t s :  126  D i s t r i b u t i o n Of  3.7.1  Introduction  3.7.2  M o d e l 61:  Longitudinal  3.7.3  M o d e l 61:  Surface  3.7.4  M o d e l 62:  Sections  3.7.5  Are  3.7.6  Conclusions  The  Predicted  Isotopes  —  127 127  6 Sections  6 Profiles And  Surface  129 133  Profiles  E f f e c t s Observable?  135 138 140  vi i i  Chapter 4.1.  4: Wave O g i v e s  142  Introduction  142  4.1.1  D e s c r i p t i o n Of O g i v e S y s t e m s  142  4.1.2  T h e o r i e s Of Wave O g i v e F o r m a t i o n  145  4.1.3  Disappearance  146  Of The Waves  4.2 Nye's T h e o r y Of Wave O g i v e s 4.2.1  O u t l i n e Of The Nye T h e o r y  4.2.2  Nye's A n n u a l l y  4.2.3  Ogives  4.2.4  An Unanswered Q u e s t i o n  4.3  Repeating  Above The F i r n  148 148 State Solution  Line  148 151 152  A New S o l u t i o n F o r O g i v e s  153  4.3.1  U s i n g Method Of C h a r a c t e r i s t i c s  153  4.3.2  Separable  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 U p s t r e a m Boundary C o n d i t i o n  157  4.3.5  The Terms Of The F l u x  158  4.3.6  Physical Interpretation  160  4.3.7  The G r e e n ' s F u n c t i o n F o r O g i v e s  167  4.3.8  A Convolution Formulation  167  4.4 Some S i m p l e  Mass B a l a n c e  Solution  Examples  168  4.4.1  Introduction  4.4.2  Example: L i n e a r V e l o c i t y  4.4.3  Example: D o u b l e S t e p  4.5  For Ogives  168 Gradient  Icefall  Model  Austerdalsbreen  169 174 176  4.5.1  Introduction  176  4.5.2  E s t i m a t e d Wave G e n e r a t i o n  177  4.5.3  Ogive S o l u t i o n For Austerdalsbreen  178  4.5.4  F i n d i n g The Wave G e n e r a t i n g  181  Region  ix  4.6 C o n c l u s i o n s List  Of  Symbols  Literature Appendix A1.1  184 186  Cited  201  1: C o n t i n u i t y  The N u m e r i c a l  A1.1.1  Model  Scheme  The C o n t i n u i t y  A1 .1.2 A M a t r i x  235 235  Equation  235  Formulation  238  A1.2 N o n l i n e a r i t y  239  A1.3 B o u n d a r y C o n d i t i o n s  243  A 1.3.1  The Upper Boundary  243  A1.3.2 The Downstream B o u n d a r y  245  A1.3.3 N o n z e r o F l u x  245  L e a v e s Downstream Boundary  A1.3.4 M o v i n g Wedge T e r m i n u s A1.4 N u m e r i c a l A1.4.1  247  Stability  254  Introduction  A1.4.2 The L i n e a r  254  Computational  A1.4.3 T h e - N o n l i n e a r A1.4.4 V e l o c i t y  Instability  Instability  Smoothing  255 257 260  A1.4.5 N u m e r i c a l 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 E q u a t i o n  265  A1.4.7 Wavenumber S p e c t r a l  267  Truncation  A1.5 A c c u r a c y  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 Al.5.4 Appendix A2.1  Error  Interpolation  Error  2: I c e T r a j e c t o r y  Introduction  274  Model  278 282 282  X  A2.2 The V e l o c i t y F i e l d A2.2.1  284  The R e c t a n g u l a r  Flow M o d e l  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 R a t e  286  A2.2.4 V e l o c i t y Normal To The Bed  288  A2.3 I c e D i s p l a c e m e n t A2.3.1  Four  Point  Field  290  Interpolation  290  A2.3.2 D i s p l a c e m e n t s A t M e s h p o i n t s A2.4  Ice P a r t i c l e  Trajectories  295  Procedure  295  A2.4.1  Tracking  A2.4.2  P a r t i c l e s Which Reach I c e S u r f a c e  A2.4.3 T r a c k i n g  Backwards In Time  A2.4.4 Boundary C o n d i t i o n Appendix A3.1  A t U p s t r e a m End  3: A s p e c t s Of D i s c r e t e D a t a  Series  A4.1  Firn  298  299 Of G l a c i e r I c e  As E q u i v a l e n t  A4.2 C o n s t a n t D e n s i t y Appendix  297  298  A3.2 A l i a s i n g 4: D e n s i t y  296 296  The Z T r a n s f o r m  Appendix  291  Ice Thickness  Assumption  5: C o n t i n u i t y E q u a t i o n  A5.1 Mass C o n s e r v a t i o n  F o r An I c e Mass  In A M o v i n g C o n t i n u u m  302 302 303 306 306  A5.2 I n 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 C h a n n e l  310  A5.4  311  In A R e c t a n g u l a r  Channel  A5.5 I n B e d - n o r m a l C o o r d i n a t e s Appendix  6: E q u a t i o n s  Appendix  7: V e l o c i t y E q u a t i o n  A7.1  Introduction  For Perturbations F o r An I c e Mass  313 316 318 318  A7.2 The Shear S t r e s s  Equation  319  A7.3 A p p r o x i m a t i o n s  322  A7.4 Shape F a c t o r s  328  A7.5 S h e a r S t r a i n R a t e  330  A7.6 I c e F l u x  332  Appendix  And A v e r a g e V e l o c i t y  8: G l a c i e r S l i d i n g  334  A8.1 Measurements  .  A8.2 P h y s i c a l P r o c e s s e s  334  In S l i d i n g  335  A8.3 Computer M o d e l s Of S l i d i n g A8.3.1  Using  338  Weertman S l i d i n g  A8.3.2 Budd-Mclnnes A8.3.3 S l i d i n g  338  Model  In T h i s  338  Study  342  Appendix  9: B u r g e r s ' E q u a t i o n  343  Appendix  10: M a t r i x  346  Appendix  11: C o n v e r g e n c e  Appendix  12: M a c h i n e  Appendix  13: D i f f e r e n c i n g Scheme F o r The F l u x  A13.1  Coefficients Criteria  350  Roundoff E r r o r s  354  Introduction  A13.2 P e r t u r b a t i o n  Gradient  .... 356 356  Equations  357  A13.3 S p a c e D i f f e r e n c i n g Scheme  358  A13.4 T r a n s f e r  360  Function  A13.5 C o n d i t i o n s  On The Mesh I n t e r v a l  Appendix  14: I c e S u r f a c e  Appendix  15: A n a l y t i c M o d e l s Of I c e S h e e t s  A15.1  Elevations  Introduction  361 364 366 366  A15.2 Nye I c e Sheet M o d e l  366  A15.3 N a g a t a  367  A15.3.1  I c e S h e e t Model  Basic  Equations  367  xi i A l 5 . 3 . 2 Ice D e p t h , Mass B a l a n c e ,  And V e l o c i t y  A15.3.3 Streamlines A15.4 H a e f e l i - P a t e r s o n  368 370  Ice Sheet Model  Appendix 16: T e s t s Of The Budd-Mclnnes Model  374 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 M o d e l : N o n l i n e a r 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 C o n s e r v a t i o n  381  A16.5 B r u a r j S k u l l (Model 2 ) : Steady S t a t e F l u x  386  A16.6 Fedchenko G l a c i e r : Steady S t a t e F l u x  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 C o n c l u s i o n s  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  415  1570 To 1840  A17.2.1 E a r l i e s t  Pioneers  415  A 1 7 . 2 . 2 H . B. De Saussure  416  A 1 7 . 2 . 3 Rendu  418  A17.3 1840 To 1915  419  A17.3.1 Louis Agassiz  419  A17.3.2 J .  D. Forbes  421  A 1 7 . 3 . 3 John T y n d a l l  424  A 1 7 . 3 . 4 Many Wondrous T h e o r i e s  427  A17.3.5 Ice Deformation Experiments  429  A17.3.6 Mathematical G l a c i o l o g y . . . -  431  A 1 7 . 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  xi i i A17.4.1  Introduction  436  A17.4.2 Shear P l a n e S l i p A17.4.3 C o n t i n u u m  Or V i s c o u s  Flow?  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  Studies:  1934 S p i t z b e r g e n  A17.4.6 F i e l d  Studies:  Jungfraujoch  A17.4.7 E x t r u s i o n  Expedition  444  Research Party  .... 446  Flow Condition  437  449  Appendix  18: S t a b i l i t y  F o r A Surge Bulge  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  454 457  xiv LIST OF TABLES 2.1. Parameters For V a t n a j o k u l l  ( F i g u r e 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  82  Boundary Test  3.1. Parameters For S t e e l e Steady State  104  3.2. V e l o c i t y P a t t e r n For S t e e l e Surge  120  3.3. Numerical Time Steps For S t e e l e Surge  124  4.1. Odinsbreen: L i n e a r V e l o c i t y Approximations  177  A1 5.1 . Parameters For Nagata Ice Sheet  374  A16.1. Parameters For V a t n a j o k u l l 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 S t e e l e G l a c i e r T r i b u t a r i e s (a)  457  A19.2. Ice Flux From S t e e l e G l a c i e r T r i b u t a r i e s (b)  459  A19.3. T r i b u t a r i e s : E f f e c t On Mass Balance  460  XV  L I S T OF  FIGURES  Frontispiece  .. Cross-section  1.1.  Rectangular  1.2.  Vertical  2.1.  N u m e r i c a l Scheme At I c e D i v i d e  47  2.2.  The Wedge T e r m i n u s  48  2.3.  Example  2.4.  Filter  2.5.  C o n t i n u i t y Test  with Moving  Terminus  59  2.6.  Contin uity Test  With Moving  Terminus  61  2.7.  Nonlinear  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  75  Prism  Flow  2  For C o n t i n u i t y  Of The N o n l i n e a r  20 Interpretation  Instability  To S u p p r e s s The N o n l i n e a r  Test  With Burgers'  52  Instability  Equation  Steady I c e Sheet Ice Sheet  55  65  2.10.  Growth  2.11.  Velocity Field  2.12.  T r a j e c t o r i e s In N a g a t a M o d e l  79  2.13.  Surface  81  3.1.  Icefield  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  98  3.4.  Sliding  3.5.  S t e e l e Model  1 Growth  3.6.  Steady State  Streamlines  3.7.  Model  3.8.  Reference Surface  3.9.  Sliding  3.10.  Of N a g a t a  23  76  F o r Nagata Model  Mass C o n s e r v a t i o n Ranges L o c a t i o n  77  Test  Map  86  1 For Steele G l a c i e r Model  For S t e e l e G l a c i e r  101  To S t e a d y S t a t e F o r Model  1  2 For S t e e l e G l a c i e r F o r 6-x  Post-surge  105 106  Function  V e l o c i t y : S t e e l e G l a c i e r Model  P r e - And  104  Profiles:  47 Y e a r P e r i o d  118 121 122  xvi 3.11.  P r e - And P o s t - s u r g e P r o f i l e s : 97 Year P e r i o d  123  3.12.  P r e - And P o s t - s u r g e P r o f i l e s :  124  3.13.  S t e e l e G l a c i e r T h i c k n e s s : One Surge C y c l e  125  3.14.  Ice T r a j e c t o r i e s  126  3.15.  Longitudinal 6 Sections:  3.16.  Model 6 1 : S u r f a c e 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 S e c t i o n s  136  3.18.  Model 62: S u r f a c e 6 P r o f i l e s  137  4 . 1 . Austerdalsbreen  147 Year P e r i o d  For 97 Year Surge P e r i o d Model 61  130  V e l o c i t y And Mass Balance  4.2.  The C h a r a c t e r i s t i c s  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 F a c t o r s  166  4.6.  Ogives From A V e l o c i t y G r a d i e n t  170  4.7.  Flow Past A V e l o c i t y G r a d i e n t : N u m e r i c a 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.  O d i n s b r e e n : 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.  154  G e n e r a t i n g Waves  Austerdalsbreen  4.11. Austerdalsbreen 4.12.  In T - t Space  151  Wave Ogives  179  Ice T h i c k n e s s  180  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 N o n l i n e a r I n s t a b i l i t y  260  A1.4. Transfer  F u n c t i o n s Of Smoothing Schemes  261  A 1 . 5 . Transfer  F u n c t i o n s : Slope-dependent  266  Damping  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  A1.7. Transfer  F u n c t i o n Modulus For V a r i o u s 6  270  A1.8. Transfer  F u n c t i o n Phase Comparison  272  xvi 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 R e c t a n g u l a r 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  290  A2.4. Interpolation  Scheme  Surface f(P)  292  A 2 . 5 . C e l l V e r t e x N o t a t i o n For Downward V e l o c i t y  294  A 2 . 6 . D i s p l a c e m e n t F i e l d In A Steady S t a t e  295  A 2 . 7 . C e l l V e r t e x N o t a t i o n For N e g a t i v e Time  297  A3.1 . The Z Plane  299  A 3 . 2 . S i g n a l s W i t h Wavelengths  1.5Ax And  300  3.0AX  A 4 . 1 . F o r c e Balance On An Ice Element  303  A 5 . 1 . S u r f a c e s f o r D e r i v a t i o n of C o n t i n u i t y E q u a t i o n  307  A 5 . 2 . The T h i n C r o s s - s e c t i o n L i m i t  309  A 5 . 3 . C o o r d i n a t e s And V a r i a b l e s A10.1. Quantities  In R e c t a n g u l a r Channel  In Slope C a l c u l a t i o n  ....  313 348  A 1 3 . 1 . Damping U s i n g The Ice S u r f a c e Slope  359  A14.1 . Ice S u r f a c e E l e v a t i o n  364  A 1 5 . 1 . Nagata Steady Ice Sheet  373  A16.1. Vatnajokull: 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  387  (Model 2)  A16.5. B r u a r j d k u l l Flux Test  389  A 1 6 . 6 . Fedchenko G l a c i e r F l u x T e s t  393  A 1 6 . 7 . Fedchenko  396  Steady S t a t e Ice P r o f i l e s  A 1 6 . 8 . Fedchenko N o n s l i d i n g Model  397  A 1 6 . 9 . Fedchenko  399  n o n s l i d i n g model w i t h n=3  A16.10.  Growth Of Fedchenko G l a c i e r  400  A16.11.  Fedchenko Growth: Other N o n s l i d i n g Models  402  xvi i i A16.12 . Fedchenko G r o w t h : Moderate S l i d i n g  404  A16.13.  406  Fedchenko G l a c i e r :  0-V P l a n e  A 1 6 . 1 4 . Fedchenko G r o w t h : S l i d i n g Model  (2)  407  A 1 6 . 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  A18.1.  454  Advancing  Surge Bulge  xix ACKNOWLEDGEMENTS In the e a r l y s t a g e s of t h i s work, encouragement from my a d v i s o r going,  G. K. C. C l a r k e  helped  me  and support  keep  my  research  e s p e c i a l l y d u r i n g the p e r i o d s when the computer  were u n c o o p e r a t i v e . H i s enthusiasm  is  contagious,  programs  and  it  is  delightful  t o work w i t h h i m . My a p p r e c i a t i o n of him grows a l o n g  with  knowledge  my  W. S. B. P a t e r s o n ,  of  glaciology.  W. D. H i b l e r I I I  have  M. C. Quick, shed  light  R. A . B i n d s c h a d l e r ,  on  a s p e c t s of  dynamics and n u m e r i c a l methods. D u r i n g the l a t e project,  P . K. F u l l a g a r ,  adjacent o f f i c e , shift"  helpful and  completing  has been a f r e q u e n t  exchanging  put i n l o n g  ideas  hours  his  stages  companion  on  the  of  manuscript  this i n the  "the  and mutual encouragement.  proofreading  glacier  Ph.D. thesis  night  B. B. Narod  and  s u g g e s t i o n s , and h e , t o g e t h e r w i t h a group of  friends,  with  B. B. N a r o d , W. H . Mathews, D. W. O l d e n b u r g ,  C. F . Raymond, R. D. R u s s e l l , and  Discussions  offering colleagues  h e l p e d me t h r o u g h t h a t l a s t h e c t i c n i g h t when the  m a n u s c r i p t was p r o d u c e d . J u l i a Forbes, support,  has  personal  goals.  particular,  through  helped I  me have  Barry Narod,  her to made  encouragement,  formulate good  and  interest,  and  t o work toward my  friends  at  U.B.C.  In  J . G . N a p o l e o n i , Bo Chandra, and P e t e r  F u l l a g a r have shared many happy  non-academic  experiences  with  me. I Canada  was  supported  postgraduate  i n p a r t by a N a t i o n a l Research C o u n c i l of scholarship,  and  by  an  H. R. M a c M i l l a n  F a m i l y 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 C o l u m b i a .  The C e n t r e at  computations  have  been  c a r r i e d out a t the Computing  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 S u r v e y , f r o n t i s p i e c e photograph of the T r i m b l e G l a c i e r .  provided  the  FRONTISPIECE: Wave o g i v e s on the North Branch, Trimble G l a c i e r , Alaska Range, 61°40'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 S u r v e y , 1965.  2  3 CHAPTER j _ : BEGINNINGS ' " I t s the j o b 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 g a f f e r used t o say.'1  1.1  INTRODUCTION  1.1.1  ENDS AND MEANS  Although far-removed  glaciers  and  ice  sheets  often  from most d a y - t o - d a y m a t t e r s ,  British  destruction.  and  For  example, a s e r i e s of f o u r i c e a v a l a n c h e s  1819  destroyed  caused  damaging  often  drain  floods  (reprinted  1970),  p.  the  Sam Gamgee,  (Lliboutry  streams  the  158)  between  or  ice  dam  and  others,  rivers;  the  in  1977).  resulting  ( C l a r k e and Mathews 1981; fails.  Cunningham  (1854  100) r e p o r t e d damaging f l o o d s on the Indus  R i v e r i n the n i n e t e e n n t h described  p.  from  i n t o moraine-dammed l a k e s  catastrophically  i n p r e s s ) when  1840,  of  many b u i l d i n g s and f i e l d s and k i l l e d  A d v a n c i n g g l a c i e r s can dam  1  some  from g l a c i e r s have caused a l o n g h i s t o r y  dozens of c i t i z e n s . Ice a v a l a n c h e s  Clarke,  of  reason.  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 ,  lakes  Advances  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  Ice a v a l a n c h e s  Peru  be  dams and m i n e s . Berendon G l a c i e r ,  by F i s h e r and Jones (1971) f o r t h i s  1636  to  t h e r e are a number of  c o m p e l l i n g reasons t o study t h e i r b e h a v i o u r . g l a c i e r s would t h r e a t e n r o a d s ,  appear  century,  1818 d i s a s t e r  and  Forbes  (1845,  p.  262)  when the G e t r o z g l a c i e r dammed the  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  in  the  glacier-dammed  Val,  de  lakes  of  in  Switzerland.  Several  t h r e a t e n a proposed p i p e l i n e . r o u t e i n the  Yukon T e r r i t o r y (Canada, There i s s t i l l  Bagnes  unpublished).  no u n i v e r s a l l y a c c e p t e d t h e o r y on the  cause  i c e ages and c o n t i n e n t a l g l a c i a t i o n ; an u n d e r s t a n d i n g of  sheet dynamics h e l p s t o s e l e c t  and t e s t h y p o t h e s e s .  ice  To c o r r e c t l y  i n t e r p r e t the g e o m o r p h o l o g i c a l r e c o r d of P . l e i s t o c e n e i c e s h e e t s , we  must  erosion  and  d e p o s i t i o n . T h i s r e q u i r e s a knowledge of g l a c i e r mechanics  (e.c>.  Boulton, of  understand  1979; H a l l e t ,  growth  Paterson,  and  processes  of  glacial  1979). The volume, d i s t r i b u t i o n , and r a t e  decay  of  the  Pleistocene  ice  sheets  (§_.£.  1972) are i m p o r t a n t d a t a f o r the d e t e r m i n a t i o n of  viscosity sea  the  of  the  the upper m a n t l e , v e r t i c a l c r u s t a l movements, and  l e v e l changes  ( e . £ . Andrews,  1974).  The i s o t o p i c c o m p o s i t i o n of p o l a r i c e s h e e t s has been  used  t o r e c o n s t r u c t temperature changes and c l i m a t e over the past years  (Dansgaard  cores,  i t i s necessary  ice  sheet  Federer,  and  others,  1969).  To  c o r r e c t l y date deep  t o determine the flow p a t t e r n w i t h i n  (Dansgaard  and  Johnsen,  1971; Hammer and o t h e r s ,  10 s  I969[a];  the  Philberth  and  possibility  of  1978).  A c u r r e n t q u e s t i o n of some c o n c e r n i s  the  g l o b a l atmospheric warming due t o combustion of f o s s i l f u e l s and clearing  of  temperate  multidisciplinary Boulder, atmospheric  study  Colorado C02  Disintegration  to  would of  forests (NOAA,  (e«g_«  unpublished)  investigate have  on  the  the  SMIC, is  1971).  A  underway  in  effect  Antarctic  increased  ice  sheets.  the E a s t A n t a r c t i c Ice Sheet c o u l d r a i s e  l e v e l by 75 m and s u b s t a n t i a l l y reduce the a l b e d o of  the  sea  earth  5 (Wilson,  1969). Of more immediate c o n c e r n 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 s h e e t ; c o u l d r a i s e sea (Thomas  and  l e v e l by seven metres  others,  i s using s a t e l l i t e  in  less  than  1979). A group a t NASA (NASA,  radiometry,  a l t i m e t r y and  100 y e a r s  unpublished)  radar  imaging  m o n i t o r and t o h e l p model v a r i a t i o n s of the G r e e n l a n d i c e Nye (1951, studies  of  I952[a],  1953,  steady  flow  the  a n a l y t i c a l models, Nye (1960, developed  1961,  1957)  to  I963[a],  climate,  interesting variations  I963[b],  solutions complicated  be  can  on  variations,  be  complicated boundaries this  Many  temporal  the a n a l y t i c a l  found by n u m e r i c a l methods  using  computers. pitfalls.  particular  of  However,  numerical  A numerical s o l u t i o n  from the c o r r e c t  i s e x t r e m e l y d i f f i c u l t t o prove t h a t a a  large  some  many reasons (e.<j. Richtmyer and M o r t o n ,  solved  methods.  more  digital  to  I965[b])  these  have t h e i r own s p e c i a l  correctly  are  have  Answers  and  k i n e m a t i c waves, and  boundary c o n d i t i o n s ;  of a d i f f e r e n t i a l e q u a t i o n may d i f f e r  It  quantitative  I965[a],  perturbation  problems  used.  problems  f i n i t e differences  for  sheet.  and i c e s h e e t s u s i n g  1963[C],  using  complicated  cannot  solutions  by  glaciological or  to  and Weertman ( 1 9 5 8 ) , L l i b o u t r y ( I 9 5 8 [ b ] ) ,  the t h e o r y of g l a c i e r  response  made the f i r s t  of g l a c i e r s  this  solution  1967; G a r y ,  numerical  differential  1975).  model  equation  has with  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 ;  i s p r e c i s e l y the type of problem f o r which n u m e r i c a l  yet models  necessary. It  is  essential  analytical solutions glaciological  to for a  literature  f i r s t check n u m e r i c a l models variety contains  of very  simpler little  against  problems. discussion  The of  6  model v e r i f i c a t i o n some  indication  numerical The  in spite that  of i t s obvious  much  of  thrust  models,  p r o b l e m s , and d e v i s i n g t e s t s With  this  in  flow  Like  previous  Bindschadler,  for  temperate  driven  and C h a p t e r (Budd  have  channel  analyzed  1).  In  numerical  of  mass b a l a n c e ,  the  numerical  Chapter  2 I present to  analytical  computer  any  approximation unpublished,  local  model,  can s e r i o u s l y the g l a c i e r  of a g l a c i e r  realistic (e.S. p.  ice  differences  surface for a width  assuming  affect  and  bed  the flow i s  and a c c u r a c y o f  the  equations comparing to  check  mass c o n s e r v a t i o n ,  glacier  terminus  moves  t h e i c e t h i c k n e s s and t h e  ( S e c t i o n A 1 . 3 . 4 ) . The p h y s i c s o f  snout  and  and g l o b a l  the  rheology.  Budd  105),  finite  solutions  in  the deformation  computer  1975;  a s e t of t e s t s  w i t h a n o n l i n e a r flow law.  for  and b o t h  throughout  a new  as i s p o s s i b l e f o r n o n l i n e a r  a case  velocity  o f t h e model  Jenssen,  stability  including  i t  the  2.2).  arbitrary  mobility  incorrectly,  avoid  stresses.  solutions  a  to  glacier  terminus  If,  understanding  t o g e t h e r w i t h a f l o w law  time-dependent  t h i s model a s t h o r o u g h l y (Appendix  on  2, S e c t i o n  and  t h i s model u s e s  w i t h an a r b i t r a r y  by g r a v i t a t i o n a l  I  the  1,  models  the  i c e mass i n a  topography,  literature  the accuracy  s o l v e t h e mass c o n s e r v a t i o n e q u a t i o n give  ways  I have d e v e l o p e d  (Appendix  unpublished),  i c e , to  is  incorrect.  finding  to v e r i f y  mind,  model o f g l a c i e r several  There  o f my work h a s been aimed a t  the problems of n u m e r i c a l  to  published  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 major  results.  the  importance.  is  complicated  In  the  Jenssen,  terminus  (Nye,  standard  1975;  i s simply  1967)  numerical  Bindschadler, a wedge-shaped  7 volume w i t h s l o p e and apex chosen so as t o  conserve  mass  (see  S e c t i o n A 1 . 3 . 4 ) . The e r r o r i n u s i n g t h i s k i n e m a t i c a p p r o x i m a t i o n cannot  be  determined?  motion of a g l a c i e r flow  in  tensor  the  correct  general  s o l u t i o n for  t e r m i n u s on an a r b i t r a r y s l o p e  form  is s t i l l  comparing  the  Glen's  an u n s o l v e d p r o b l e m . However, I  have t e s t e d the n u m e r i c a l i m p l e m e n t a t i o n of the by  with  computed  solution  to  a n a l y t i c a l solution with a simpler "rheology"  wedge a  terminus  time-dependent  (Section  2.3.2).  Many s t a n d a r d n u m e r i c a l schemes f o r l i n e a r e q u a t i o n s down  when  test a  numerical  equation an  applied  analytical  with  a  solution;  it  is  solution  For some g l a c i e r and  to  flow p r o b l e m s ,  Johnsen,  calculate  1969[a])  the  law  continuity Particle  to  find  equation  to  trajectories  the time-dependent model  the  against  Burgers'  an  the  mass  results  with  dating  and f i n d i n g the  (e.£. of  horizontal  are  to  with  e q u a t i o n t o show t h a t my  Jenssen,  field  glacier, velocity,  derive  the  ice  on  cores  temperature  1977),  ice p a r t i c l e s .  velocity  then  related  such as  l o n g i t u d i n a l mesh f o r a time-dependent flow  problem.  to  problems.  n e c e s s a r y t o know the t r a j e c t o r i e s can  i s important  numerical  Burgers'  d i s t r i b u t i o n of c o l d i c e masses  model  also  I have compared  model c o r r e c t l y s o l v e s n o n l i n e a r  (Dansgaard  nonlinear  It  break  2.3.3) i s a nonlinear hyperbolic equation  conservation equation. the  nonlinear equations.  model  (Section  analytical  to  the  it  is  My computer a  vertical  by u s i n g and  vertical  Glen's  using  the  velocity.  found by a n u m e r i c a l i n t e g r a t i o n of  v e l o c i t y . I t e s t e d t h i s p a r t of the analytical  p a r t i c l e p a t h s i n a steady i c e  solution sheet.  by  Nagata  computer (1977)  for  8  T h i s n u m e r i c a l model i s p r o b a b l y the accurately  tested  assembled,  of  its  most  thoroughly  and  t y p e . The set of t e s t s which I have  or o t h e r s s i m i l a r t o them, s h o u l d be used  to  verify  any n u m e r i c a 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 c o n f i d e n c e t o s o l v e more c o m p l i c a t e d p r o b l e m s . I  have  used  Previous  efforts  isotopic  ratios  this  new  have in  concentrated  ice  cores  assuming steady s t a t e f l o w . feasibility  of  computer  using  model on  using  unsteady  Glacier,  flow.  Yukon  The  example  Territory.  characteristic  pattern  I in  p r e l i m i n a r y measurements of 6 ( 0 may be masked by o t h e r Finally, relating  the  variations  in  have  evaluated  the  i s o t o p e measurements t o 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 and  two s t u d i e s .  t o i n v e s t i g a t e c l i m a t e change,  In Chapter 3 I  stable  in  I  1 8  /0  that isotope  1 6  constant  c o n s i d e r e d was the  found the  a  )  surging  climate Steele  leaves  distribution,  a but  suggest t h a t t h i s p a t t e r n  effects.  i n Chapter 4 I have d e r i v e d a  linear  convolution  a m p l i t u d e of wave o g i v e s t o the v e l o c i t y , c h a n n e l  w i d t h and mass b a l a n c e i n i c e f a l l s . T h i s work i s an e x t e n s i o n of s t u d i e s by Nye ( I 9 5 8 [ b ] ) . the  convolution  I used the computer  model  to  verify  f o r m u l a t i o n and t o d e t e r m i n e which f e a t u r e s  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 f o r m u l a t i o n of i t s l a r g e wave o g i v e s .  to  of the  9  1.1.2 There for  CONVENTIONS USED  are many d i v e r s e views on the most a p p r o p r i a t e  a Ph.D. thesis.  fully by  describes  those  detail,  a  I  basic  have  knowledge  documented  of  all  physics  my  results  the  text  rather  than  or  which  physical  numerical  and summarized the r e l e v a n t work of  substantiate  as  My aim has been t o produce a document  my work, and which can be understood on i t s own  with  glaciology.  style  methods i n  others;  substitute  references  for  i t . This  i n a l e n g t h y m a n u s c r i p t . To keep the main t e x t as  possible,  I  have  placed  the  background m a t e r i a l i n a p p e n d i c e s . clearly  identifiable.  appreciated  numerical  methods  The work of o t h e r s  short and the  should  I hope t h a t t h i s l e v e l of d e t a i l w i l l  by some r e a d e r s , s i n c e b r e v i t y w i l l  be  required  be be in  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 f o r p u b l i c a t i o n . Unless  s t a t e d o t h e r w i s e , I have used a r i g h t h a n d e d l o c a l l y  o r t h o n o r m a l c o o r d i n a t e system such t h a t the x the  glacier  and  lies  along  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 transverse bed  axis  and h o r i z o n t a l , and the z a x i s  positive  upward  is  normal  to  the  i n the v e r t i c a l p l a n e c o n t a i n i n g the  centreline. The v e l o c i t y components are This (i..e.  notation (u,v,w)  from  the  are  used  to  i n d i c a t e s the rank of  the  dot symbol  tensor,  axes.  convention  i n the  indicate tensors.  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 The  standard  ) due t o h i s t o r i c a l developments  Underscores underscores  differs  ( u , w , v ) a l o n g the ( x , y , z )  thesis.  The number of  e.£. v  is  the  matrix.  when l o c a t e d above a v a r i a b l e ,  indicates  10 a  time  derivative.  represents  the  When  located  standard  e.cj. M a l v e r n ( 1969, p .  scalar  a  depth  two  vectors,  it  i n n e r product or dot  product,  i n d i c a t e s an average v a l u e ,  usually  17).  A bar above a v a r i a b l e over  between  range  (z  d i r e c t i o n ) , or over time (§_.£. annual  averages). A l i s t of symbols, section  in  Chapter  which  together  each  is  with  meanings  numbers,  the  i n t r o d u c e d , can be found f o l l o w i n g  are  numbers,  and  enclosed  in  textual  references  to  round p a r e n t h e s e s , e . £ .  ( A 5 . 6 ) , or ( A 1 . 1 . 3 ) . The c h a r a c t e r s  equation (2.2.5),  p r e c e d i n g the f i r s t  subsection  where the e q u a t i o n  g i v e n , and the f i n a l number i s the c o n s e c u t i v e that  subsection.  identified  as  References  such,  and  the  to chapter numbers  are  equation  sections not  are  or  decimal  a r e the c h a p t e r or appendix number. The m i d d l e number  p r e s e n t ) i d e n t i f i e s the c h a p t e r  in  and  4.  Equation  point  their  (if is  number always  enclosed  in  parentheses. The  LITERATURE  CITED  is  in  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 a b l e t o o b t a i n the use of some of the v e r y 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 than  English.  In  c i t a t i o n s of o t h e r s , brackets.  those  cases,  where  other  I c o u l d not v e r i f y the  I have i n c l u d e d the c i t i n g a u t h o r i n square  11 1.2 PREVIOUS WORK 1.2.1  FOUR CENTURIES OF GLACIER FLOW THEORY  The framework w i t h i n investigate  glacier  which  flow  we  currently  has been assembled  understand i n the p a s t  d e c a d e s . Deep c o r i n g t e c h n i q u e s (Hansen and Langway, radio-echo  sounders  (Evans,  1963),  coming  i m p o r t a n t e x p e r i m e n t s on the d e f o r m a t i o n 1955)  and  glaciers  (Nye, 1951,  I952[a],  research.  1957)  1966)  and  shortly  after  (Glen,  1952,  ice  s t a r t e d the r a p i d growth  However,  development  of  i n v e s t i g a t i o n of the flow of  g l a c i e r s goes back hundreds of y e a r s .  In Appendix 17,  I  review  of i d e a s on g l a c i e r flow i n t h a t 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 i m a g i n a t i o n  and  o b s e r v a t i o n s ; o t h e r s were c o n c i s e and l u c i d statements which  three  m a t h e m a t i c a l t r e a t m e n t of the flow of i c e s h e e t s and  glaciological  the  of  and  limited on t o p i c s  are s u b j e c t s of r e s e a r c h t o d a y . Some m i s c o n c e p t i o n s about  g l a c i e r flow were r a i s e d , d e b a t e d , and r e s o l v e d 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 r e s e a r c h d u r i n g the second half  of  the  n i n e t e e n t h c e n t u r y 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 p a r t of t h i s c e n t u r y , g e o l o g i s t s (with  a  few  notable  theories reflected glacier  flow  difference.  were  given  Contemporary  by C r o l l  p.  Perutz (1947), Since  160),  by  Hawkes  and by Orowan  1950,  field  research  (1930),  reviews  of  (1875, Chapter XXX,  p . 4 9 5 ) , by G e i k i e (1894, Chapter 3, p . 2 5 ) , by Chapter 9,  the  e x c e p t i o n s ) , and r e s e a r c h p r i o r i t i e s and  this  theory  dominated  Russell  (1897,  by Matthes (1942),  by  (1949). on  glacier  flow  has  progressed  12 r a p i d l y . The flow law f o r i c e was e s t a b l i s h e d purposes  by  Lliboutry  d 968[a],  and  Glen (1952,  Morland  glacier  unresolved. 1957,  and Nye ( 1 9 5 3 ) . Weertman (1957),  I 9 6 8 [ b ] ) , Kamb ( 1 9 7 0 ) , Nye  d 976[a],  sliding;  1955)  f o r many p r a c t i c a l  1976[b3)  some  Papers  Nye  1959[C]) e s t a b l i s h e d and  contributed  aspects  by  steady  glacier  fields,  while i d e n t i f y i n g  of  (1951,  realistic  this  question  are  I952[b],  analytical  still 1952[C],  solutions  velocities,  for  and s t r e s s  approximations.  Vialov  (1957) used G l e n ' s flow law t o d e r i v e a steady i c e sheet  profile  which  useful  1970)  t o the t h e o r y of  I952[a],  i c e sheet p r o f i l e s , many  d969[b],  matched the flow l i n e t h r o u g h 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 steady  i c e sheet p r o f i l e s ,  the bed. Weertman (1963)  and i n c l u d e d i s o s t a t i c  considered  the  rates  depression  effects  of  of  fringing  mountain ranges on steady i c e s h e e t s , and (1966) the e f f e c t b a s a l water  on  of a  layer.  A l t h o u g h 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 t o study  fully,  some u s e f u l a n a l y t i c a l r e s u l t s  K i n e m a t i c waves on g l a c i e r s were  studied  by  Lliboutry  were observed by V a l l o t d958[b]),  Nye  Weertman (1958) u s i n g p e r t u r b a t i o n methods. papers,  (1960,  1961,  have been d e r i v e d .  I963[a],  I963[b],  (1900)  (1958[a])f  I 9 6 5 [ b ] ) extended the method t o a n a l y z e the response of t o c l i m a t i c change,  for  a  thin  of g l a c i e r s .  relation  Bodvarsson  (1955)  i c e sheet and a n a l y z e d i t s  c l i m a t i c change. T h i s model assumed  of  I965[a], glaciers  and t o e s t i m a t e p a s t c l i m a t e from the r e c o r d  of advance and r e t r e a t equations  and by  Nye, i n a s e r i e s 1963[C],  and  between  is basal  not  widely  used  derived  stability due  to  to its  s t r e s s and i c e f l u x . Weertman  13 (1961[a])  performed  Weertman  (1957)  a  similar  sliding.  He  stability also  analysis  d e r i v e d ( I 9 6 4 [ a ] ) the time  s c a l e s f o r the growth or decay of a p e r f e c t l y p l a s t i c Jenssen and Radok temperature  assuming  ice  (1963) o b t a i n e d a n u m e r i c a l s o l u t i o n  field  in  the  central  region  of  an  sheet.  for ice  the sheet  undergoing t h i n n i n g . The study of f u l l y time dependent boundaries possible  and by  source  high  terms  speed  is  i c e masses w i t h  a  computers.  recent  development  Numerical s o l u t i o n s  e q u a t i o n s g o v e r n i n g i c e masses have been o b t a i n e d problems by, e.g_.  Shumskiy  Mahaffy  Jenssen ( 1 9 7 7 ) , B i n d s c h a d l e r  (1976),  (1963),  C l a r k e ( 1 9 7 6 ) . The complete  arbitrary  Budd  and  made of  the  f o r a range of  Jenssen  (1975),  (unpublished)  and  s o l u t i o n of the e q u a t i o n s of m o t i o n ,  the c o n s t i t u t i v e e q u a t i o n s and the e q u a t i o n s of s t a t e f o r 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 s o u r c e s outstanding  models  time-varying glaciers  an  which and  ice  find  the  sheets  surface  height  of  are  a  relatively  new  tool.  Campbell  and  Rasmussen  Campbell (1973) developed a on  a  horizontal  (1969, model  1970)  which  and  found  flux.  coefficient  By  Rasmussen and ice  depth  at  mesh. They assumed t h a t the i c e was a  viscous material- with a basal f r i c t i o n c o e f f i c i e n t mass  is  PREVIOUS ICE PROFILE MODELS  Computer  points  boundaries  problem.  1.2.2  research  and  arbitrarily  lowering  they s i m u l a t e d g l a c i e r  surges.  the  determined by  basal  friction  14 Budd  and  others  (1971)  developed  a f i n i t e difference  equation  for  the  numerical These  methods  authors  properties  have  of  to  thickness  similar  differ  Budd  model  glacier  W h i l e the e q u a t i o n s are  and  in  included  and solve  to  those  properties developed to  stable  will  I  additional  Appendix  16).  important  be  to  physical  of l o n g i t u d i n a l toward  included  later.  This  stress  obtaining  model  an  physical  was  later  1978,  1979)  s u r g e s . Working from the assumption  b a s a l m e l t w a t e r can cause s l i d i n g , they used the redistribute  large longitudinal strain behaviour  our  n u m e r i c a l scheme; the a d d i t i o n a l  periodic  dissipation  flowline. used,  by Budd (1975) and Budd and Mclnnes (1974,  generate  (1975)  continuity  have  some r e s p e c t s (see  i c e flow (e.g_. the e f f e c t  and  the  p r o f i l e along a  d e v i a t o r s ) and I have devoted more e f f o r t accurate  Jenssen  the  rates  basal  and  that  strain  energy  shear s t r e s s ,  causing  rapid  flow.  of the model, viewed as a q u a l i t a t i v e  The  sliding  phenomenon,  may  be i t s most important c o n t r i b u t i o n t o our i d e a s on s u r g e s . Bindschadler difference  profile  t h e s i s . Bindschadler stability  and  used  (unpublished)  developed  another  finite  model s i m i l a r t o the one I d e s c r i b e i n t h i s also did a careful a  numerical  analysis  scheme  d i s c u s s i n Appendix 1. He used t h i s  model  of  numerical  s i m i l a r t o the one I to  investigate  changes i n the s u r g e - t y p e 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  the its  quiescent phase. Mahaffy Andrews 1976; finite  (Mahaffy,  unpublished;  Andrews and M a h a f f y ,  difference  Mahaffy  1976)  used a  1976; Mahaffy and two-dimensional  model t o study the i c e t h i c k n e s s  e x t e n t of the L a u r e n t i d e  i c e sheet and the Barnes i c e  and  lateral  cap.  15  Jenssen (1977) p u b l i s h e d the o n l y f u l l y model  of  ice  temperature limited  sheets.  This  as w e l l as  by  ice  model  calculated  surface  severe computer  three  height.  dimensional  flowlines  and  accuracy  was  Its  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 t o r e p r e s e n t the whole Greenland Ice However, t h i s model i s an a m b i t i o u s will  development,  be f o l l o w e d by o t h e r models of t h i s  1.2.3  p.  late  kind.  streamlines  flow l i n e s near the  Greenland  by  thickness,  d e n s i t y and mass b a l a n c e ,  independent  assuming  (1)  steady  of d e p t h . H a e f e l i  same  solution.  ice state,  (2)  of  glaciers  smooth  fields.  (1894)  All  1978,  central  constant  independently  and  ice  these  derived  Finsterwalder  (1897)  streamlines  u s i n g the concept of flow tubes and  vector  qualitative;  in  to  and (3) h o r i z o n t a l v e l o c i t y  (I963[b])  Reid  divide  i n d e p e n d e n t l y developed a method of c a l c u l a t i n g steady  d a t e back  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,  133) c a l c u l a t e d  the  probably  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 the  and  Sheet).  methods  in  properties were  only  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  ice. Haefeli p o r t i o n of sliding,  (1961) a  steady  derived  the v e l o c i t y f i e l d i n the c e n t r a l  isothermal  ice  sheet,  assuming  and (2) d e f o r m a t i o n by shear p a r a l l e l t o the bed,  G l e n ' s flow law f o r i c e ( G l e n , Nagata  (1)  no  using  1955).  (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 d e f o r m a t i o n ;  the flow was a l l  16 b a s a 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 By  assuming  (see  A8.3.1).  that  the v e r t i c a l v e l o c i t y was c o n s t a n t a l o n g the  i c e 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 t o model 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 by flow  in  Antarctica  (Nagata,  1978).  I use t h i s model as a t e s t  f o r my n u m e r i c a l s t r e a m l i n e c a l c u l a t i o n s i n S e c t i o n N i e l s o n and S t o c k t o n  ( 1956)  glacier  derived^  the  2.5.2.  flow  field  in  v a l l e y g l a c i e r s of c o n s t a n t v a l l e y c r o s s - s e c t i o n assuming  steady  plastic  stress  flow,  and  Shumskiy (1967) found a s o l u t i o n f o r  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 . S e v e r a l t r a j e c t o r y models have  been  derived  for  regions  near i c e d i v i d e s on steady  i c e s h e e t s i n order t o date i c e c o r e s  (Dansgaard  I969[a];  and  Johnsen,  Hammer and o t h e r s , (Weertman,  P h i l b e r t h and F e d e r e r ,  1978), and t o model  1968).  temperature  with  s t a t e A n t a r c t i c flow l i n e s assuming  that  r a t e was c o n s t a n t i n any v e r t i c a l column, velocity  variation  p r o c e s s does not appear t o s a t i s f y surface  elevation  described in Section Jenssen difference  (1977)  was  about  gradients.  Budd and o t h e r s (1971) c a l c u l a t e d t r a j e c t o r i e s a l o n g  horizontal  depth  A l l these models make some assumptions  v e r t i c a l s t r a i n r a t e s or temperature  the  1971;  with  the  vertical  by  strain  (p. 51) or weighted by depth  (p. 5 5 ) . T h i s  continuity locally.  calculated  steady  the  The  numerical  ice  model  1.2.2. calculated  trajectories  in  a  finite  t h r e e d i m e n s i o n a l time-dependent i c e sheet model, i n  o r d e r t o s o l v e f o r the temperature  field.  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 f o r a s u r g i n g model of a f l o w l i n e t h r o u g h M i r n y ,  Antarctica  by  17 Budd  and Mclnnes ( 1 9 7 8 ) . He d i d not d e s c r i b e  the method used  o b t a i n the f l o w l i n e s . T h i s i s the o n l y p r e v i o u s model, of I  am aware,  t o c a l c u l a t e the t r a j e c t o r i e s  to  which  in a time-varying ice  mass.  1.3 THE CONTINUITY EQUATION FOR AN ICE MASS  1.3.1  THE GENERAL GLACIER FLOW PROBLEM  The e q u a t i o n of c o n t i n u i t y e x p r e s s e s the the input  ice  mass  subject  and s l i d i n g of i c e ( S e c t i o n s  the assumptions derived  the  about the  continuity  flow  field  equation  Prager,  (e.c[.  1973).  equation,  the  Truesdell  In  this  and  for  described [x,y,z],  which  coordinate  system  physics  an  below.  and w i t h I  have  i c e mass from f i r s t  methods  of  1960;  will  of  continuum  Malvern,  give  the  1969;  resulting  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 The  I  the  imposed  Toupin,  Section,  assumptions  to  1.4 and 1 . 5 ) ,  p r i n c i p l e s i n Appendix 5, u s i n g s t a n d a r d mechanics  in  changes i t s shape over t i m e , i n response to mass  ( a c c u m u l a t i o n or a b l a t i o n ) ,  deformation  manner  terms. I  have  used  in  this  study  is  i n S e c t i o n 1 . 1 . 2 . The p o s i t i o n v e c t o r x i s the t r i p l e t and  the  velocity  vector  v(x)  is  the  triplet  [u(x),w(x),v(x)]. Three b a s i c assumptions 1. m a t t e r i s  of the g l a c i e r model a r e :  conserved.  2. momentum i s c o n s e r v e d , i . e . a c c e l e r a t i o n 3. i c e i s i n c o m p r e s s i b l e  (see  Appendix 4 ) .  negligible.  18 To  solve  the  constitutive  general  equations,  set  of  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  all  the  components  of  t e n s o r and the v e l o c i t y f i e l d ,  with  stress  boundary c o n d i t i o n s on problem  to  instill  enthusiastic  start  channel  geometry  strain  sense  humility  boundary, in  (boundaries  itself  etc.).  and  or- n u m e r i c a l  analyst.  symmetry),  or  flow  about  the  about  the  i s o t h e r m a l ) , a n d / o r about strain,  a  the most  to f i n d s o l u t i o n s to ice  (e.g_.  ( e . £ . plane  is  even  by making some a d d i t i o n a l assumptions  distribution  rates,  of  arbitrary  of  of which I am aware,  problems  flow f i e l d  a  possibly  and o p t i m i s t i c g l a c i o l o g i s t  A l l attempts,  temperature  a  equations  equations,  state  the  and  conservation  simple  The model I d e s c r i b e  shear,  the  uniform  i n t h i s study i s no  exception.  1.3.2  RECTANGULAR CROSS-SECTION FLOW MODEL  The g l a c i e r flow volume b e i n g m o d e l l e d (see 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 ,  The two d i m e n s i o n a l dimension  in  an  model  includes  approximate  way  F i g u r e 1.1)  and a w i d t h W ( x ) .  variations  in  the  by the assumptions  v e l o c i t y components u and v are independent of y , and lateral  is  third  that  the  that  the  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  rate  independent of y ) , such  that  the  margins i s p a r a l l e l t o the m a r g i n s ,  net i..e.  velocity  y_(x)  at  the  19 dv = 0 (1.3.2)  dy dv = u ( x , z , t ) dW dy The g l a c i e r  thickness  W(x)  h(x,t)  dx  (1.3.3)  i s a l s o independent  dh(x,t)  of  y.  = 0  dy Nye  (1.3.4)  ( 1 9 5 9 [ C ] , equation  lateral  variations,  incorporated  it  form of t h i s For possibly  (33))  and Budd and Jenssen (1975, into  their  an  ice  sheet,  ( 1 . 3 . 4 ) are  reasonable.  For  a  important,  valley  attempt  to  and  the  the  glacier,  identify  (1.3.1)  1.1  to  'walls'  the  of the c h a n n e l are a  assumptions  drag  W(x)  from  with  the  may  (1.3.1)  the  and g l a c i e r  (1.3.4)  centreline,  valley  are but  reasonable, with  through  walls  is  bed v a r y w i t h y .  If  valley  width,  the  the and  lateral  ice I  the  be g r o s s l y v i o l a t e d . flowlines  §_.£. W(x) may be a few p e r c e n t  v a l l e y w i d t h at the l e v e l of assumptions  3.34)  illustrates  I l e t W(x) be the d i s t a n c e between two  glacier  solution,  equation  can be the d i s t a n c e between two  and the i c e t h i c k n e s s  assumptions  the  W(x)  nonparallel flowlines; fiction,  however,  model. Figure  including  model.  mathematical  I  suggested t h i s method of  surface,  then  obtain a central  variation  in  If, near  of  all  the the  flowline  valley  width  i n c l u d e d t o a good a p p r o x i m a t i o n . The  effect  an approximate  of the v a l l e y s i d e w a l l d r a g can be i n c l u d e d i n  way by u s i n g  modify the shear s t r e s s  (see  shape  factors  equation  (Nye,  (1.4.25)  1965[C])  below).  to  20  FIGURE 1 . 1 . R e c t a n g u l a r C r o s s - s e c t i o n F l o w . The t r i a d x - y - z shows the c o o r d i n a t e a x e s , and the bold arrows u, v, and w are the v e c t o r components of the v e l o c i t y v . The example shows the v e l o c i t y in the a c c u m u l a t i o n zone (v i s n e g a t i v e ) .  1 . 3 . 3 THE CONTINUITY EQUATION With  these a s s u m p t i o n s , the well-known c o n t i n u i t y e q u a t i o n  (the d e r i v a t i o n of which I show i n Appendix 5 ) i s dh(x,t)  +  at  where h ( x , t ) the  ice  1 wTx)  g_Q(x,t) = A ( x , t )  i s the i c e t h i c k n e s s  flux  (1.3.5)  dx.  normal t o the bed,  through a c r o s s - s e c t i o n  the source term A ( x , t )  i s the net mass  bed,  accumulation  the  equivalent surface  net  thickness  u n i t s per u n i t  or  Q(x,t)  is  from bed t o s u r f a c e ,  and  balance  normal  to  the  rate  in  ice  ablation  time,  including  m e l t i n g , and b a s a l m e l t i n g and r e f r e e z i n g .  snowfall,  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 s m a l l (e.g..  Rothlisberger,  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)  = W(x) I  Q(x,t) where u ( x , z , t )  u(x,z,t)  J 0  dz  (1.3.6)  i s the v e l o c i t y component p a r a l l e l t o the bed.  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 ) Q(x,t) where  V(x,t)  is  can a l s o be w r i t t e n  = V ( x , t ) h ( x , t ) W(x)  u(x,t),  (1.3.7)  i . . e . the downslope v e l o c i t y  averaged between the bed and the  . . .  h(x,t)  u(x,z,t)  Q(0,t)  dz  (1.3.8)  0  Assuming t h a t the upstream end of the g l a c i e r c o n s i d e r a t i o n i s at x=0,  u(x,z,t)  surface.  J . . .  s e c t i o n under  the boundary c o n d i t i o n i s = Q (t) (1.3.9)  0  I f x=0 a c t u a l l y r e p r e s e n t s the p h y s i c a l upper e x t e n t of the mass,  Q0(t)  It  ice  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 a c h i e v e d by s e t t i n g V(0,t) =0 by a v a n i s h i n g i c e letting  ice  surface  thickness  (1.3.10) slope  h(0,t)  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 ,  angle  vary.  (see  (1.4.38)),  and  For the case of a v a l l e y  the boundary  condition  (1.3.9)  i s a c h i e v e d by s e t t i n g h(0,t) If  the  lower  =0  (1.3.11)  end x=L(t) i s the g l a c i e r t e r m i n u s , then L ( t )  d e f i n e d i m p l i c i t l y by  is  22 h(L(t),t) =0 (This  (1.3.12)  i s not a m a t h e m a t i c a l  boundary c o n d i t i o n , but  a  physical  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) (1.3.13)  0  is  also  required.  Ho(x)=0,  i_.e. u n g l a c i e r i z e d 1.3.4  y  (1.3.5)  may  be  interpreted  d i r e c t i o n , and t h i c k n e s s  the  in  shorter).  is  following 1.2,  w i t h w i d t h W(x) i n  6x i n the x d i r e c t i o n . L e t p be i c e . When ( 1 . 3 . 5 ) i s  c o n s t a n t />W(x)6x6t, the f i r s t  mass i n the p r i s m i n a time 6t  the  i c e as shown i n F i g u r e  from the bed t o the s u r f a c e h ( x , t ) ,  the c o n s t a n t d e n s i t y of g l a c i e r by  i n i t i a l condition  ground.  C o n s i d e r a v e r t i c a l p r i s m of  extending the  one s i m p l e  PHYSICAL INTERPRETATION  Equation manner.  For example,  (the  term i s the net change i n  prism  The second term on the l e f t  multiplied  is  then  taller  i s the d i f f e r e n c e  or  i n mass  between t h a t which flowed out of the p r i s m through the  downslope  f a c e , and t h a t which flowed i n t o the p r i s m through the  upstream  face,  in  the  time  6t.  p r i s m i n t o the downstream  T h i s i s the net l o s s of mass from the f l o w . The term on the  right  side  j u s t the t o t a l mass added t o the p r i s m i n time 6t by s n o w f a l l m e l t i n g . Thus, to  the  top  ( 1 . 3 . 5 ) s t a t e s t h a t the t o t a l of  the  by s p a t i a l  flow  variations.  or  l a y e r of mass added  p r i s m at any p o s i t i o n x i s the sum of  s n o w f a l l onto the p r i s m and the net mass l e f t  is  inside  the  the  prism  23  1.4 PHYSICS OF ICE DEFORMATION 1.4.1 In  INTRODUCTION  this Section,  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  equation r e l a t i n g g l a c i e r c o n t i n u i t y equation in  this  glacier  t h i c k n e s s and i c e  flux  derivation.  T h i s i s a s t a n d a r d procedure  flow (e.g.. P a t e r s o n ,  and body f o r c e s and the a c c e l e r a t i o n s  observations  ice establish applied  to  the  and  steps  in modelling  1980; Raymond, 1980).  The s t r e s s e q u i l i b r i u m e q u a t i o n s Second,  that  ( 1 . 3 . 5 ) can be s o l v e d . There a r e t h r e e  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 surface  so  between  the  i n any c o n t i n u u m .  are o u t l i n e d i n S e c t i o n  1.4.2.  t h e o r y of the d e f o r m a t i o n of  glacier  constitutive relationships  between  the  stresses  i c e , and the r e s u l t i n g d e f o r m a t i o n . G l e n ' s flow law  24 for ice (Glen, Finally,  in Section  substituting flow  1955,  law  the  to  appropriate  1958)  is  1.4.4, after  presented  the  strain  coordinates,  I  of S e c t i o n  rate,  and  1.4.2  and  1.4.3.  assumptions, into Glen's  integrating  over  e x p r e s s the downslope component of  the i c e 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 . thickness,  Section  making s i m p l i f y i n g  s t r e s s equations  get  in  Since ice  v e l o c i t y are r e l a t e d through ( 1 . 3 . 7 ) ,  flux,  this  will  complete the d e r i v a t i o n of a second e q u a t i o n needed t o s o l v e c o n t i n u i t y equation  1.4.2  (1..3.5) f o r i c e t h i c k n e s s and f l u x .  STRESS EQUILIBRIUM EQUATIONS  Since g l a c i e r Newton's  the  second  i c e deforms  s l o w l y , the a c c e l e r a t i o n  law can be n e g l e c t e d ,  term  l e a v i n g the r e s u l t  in  that,  f o r any volume V of a s l o w l y deforming c o n t i n u o u s medium, B + T = 0 where*B i s specific  the body  total forces  body  force  (1.4.1) found  by  integrating  ( f o r c e per u n i t mass) f.(r)  over a l l r_ throughout the volume V . I t s components  the  at p o s i t i o n £ are  B V  (1.4.2)  and T i s the t o t a l s u r f a c e t r a c t i o n on the s u r f a c e the  volume  S  enclosing  V w i t h s u r f a c e normal v e c t o r n . The components of T  are T  =||  n (r) S  tf.k  i s the s t r e s s t e n s o r ,  surface  normal  to  i_.e. the f o r c e  c  jk  (r)  per  dS  unit  (1.4.3) area  on  a  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 r e p e a t e d  indices  used i n t h i s s e c t i o n . O r t h o g o n a l axes are Applying  Gauss'  Theorem  ( 1 . 4 . 3 ) and s u b s t i t u t i n g components  (e.g..  ( 1 . 4 . 2 ) and  are  summed,  assumed.  Prager, (1.4.3)  1973, into  V  (1.4.1)  1  by  (1.4.4) (1.4.4)  has  counterpart (r)  Assuming  the  to  d3r = 0  S i n c e the volume V i s a r b i t r a r y , the g l o b a l e q u a t i o n  setting  p . 29)  gives  dx  a local  is  f (r)  that  the  +  Ic  dx  net  k i(r) = 0  (1.4.5)  k angular a c c e l e r a t i o n  i s z e r o , and  t o z e r o the sum of moments a c t i n g on the volume V , g i v e s  result c ij  . i . e . the s t r e s s t e n s o r similar (1973, p .  to  (1.4.1)  (1.4.6)  = <t J i  i s symmetric. through  The  (1.4.5),  development and  is  very  i s g i v e n i n Prager  47).  For a g l a c i e r ,  the o n l y body f o r c e  When the x , a x i s  f(r)  = g  (to  be  is gravity,  (1.4.7) called  x)  glacier  bed which i s at an a n g l e fix)  a x i s (z)  i s normal t o i t and upward, and  horizontal,  so  then the s t r e s s e q u a t i o n s  is  taken  along  the  t o the h o r i z o n t a l , the x 3 the  x2  (1.4.5)  are  axis  (y)  is  26 <^ff  he  dx  dy  Tz  ^ff  <5ff  dff  off  x z + pgsin(p)= 0  (1.4.8)  xz + yz + z z - />gcos(*) = 0 <^x by £z  (1.4.9)  do dx  1.4.3 The  xx +  xy_ +  xy +  5ff dy  Yl  +  ^ff vz  =  Tz  0  (1.4.10)  CONSTITUTIVE RELATION FOR ICE DEFORMATION  constitutive  equations  for  deformation  s t r e s s e s a p p l i e d t o i c e t o the r e s u l t i n g d e f o r m a t i o n components of the s t r a i n r a t e t e n s o r  1  relate  the  rate.  The  are  — i + —j dx dx • j i  D  (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 t h a t the most g e n e r a l  the s t r e s s t e n s o r and the s t r a i n r a t e t e n s o r  r e l a t i o n between  for  a  nonlinear,  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 €  ij  = A(T ,T ,T )6 + B(T ,T ,T )c 1 2 3 ij 1 2 3 ij + C(T ,T ,T )<s a 1 2 3 ik kj  where 6- j i s the Kroenecker t  delta  (1.4.12)  27 6  i j  and  T1f  T2,  and  =0  i^j  = 1  i=j  T3  are  (1.4.13)  the f i r s t ,  i n v a r i a n t s of the s t r e s s t e n s o r T  - M2<y  6  3  higher  scalar  a  c  i j jk k i  powers  Hamilton-Cayley stress-dependent  (1.4.14)  = J_( <r e ~ c a ) 2 ij ji i i jj  2  in  1973, p . 2 2 )  ii  T  Terms  (e.g.. P r a g e r ,  scalar  = c 1  T  second, and t h i r d  of  equation  - 3*  •  (1.4.15) c  + t f * t f )  e  i j j i kk can  i i j j kk  be  eliminated  (§_.g_. P r a g e r ,  1973,  (1.4.16) by  the  p. 25).  The  at most,  the  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 ,  i n v a r i a n t s , because the r e l a t i o n ( 1 . 4 . 1 2 ) i s  independent  of the c h o i c e of a x e s . The  coefficients  A,  B,  and  C  may  be  dependent  on  temperature. R i g s b y (1958) showed t h a t , deformation  rate  of  h y d r o s t a t i c pressure  ice  to  a  crystals  good was  approximation, independent  of  the the  p , where P = 1 «  3 when  the  ice  pressure-melting  temperature  ii was  (1.4.17) measured  relative  the  p o i n t . T h i s r e s u l t means t h a t 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 i m p l y i n terms of stress *|j  to  the  deviatoric  28 cr '  ij  By  e  =  definition,  - _1_ 6  ij  3  the  ice  is  constant  (see  ij  first  z e r o . By f u r t h e r a s s u m p t i o n s ,  cr  kk  (1.4.18)  scalar  i n v a r i a n t T', of cr'j  f i r s t t h a t the d e n s i t y of  Appendix  4),  second, t h a t ,  s t r e s s , the components of s t r a i n r a t e are components  of  the  glacier  for a given  proportional  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 ,  (see  Glen,  1958)  the  flow  When €  and  invariants  T  are  of  the  = B(T' )  ij  and  of  T2  to <r'  ij  2  square  the  law p r e s e n t e d by Nye (1953)  reduced the g e n e r a l r e l a t i o n ( 1 . 4 . 1 2 ) €  to  t h a t the  second i n v a r i a n t of the s t r a i n r a t e t e n s o r i s a f u n c t i o n only  is  roots  (1.4.19) of  the  second  scalar  (T = JT 2 ) , Nye (1953) showed t h a t 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)  U s i n g the case of a s i m p l e s l a b d e f o r m i n g by shear the x a x i s ,  parallel  to  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  *' i j  i j The  exponent  n i n ( 1 . 4 . 2 0 ) i s independent of  V a l u e s i n the l i t e r a t u r e v a r y 1952)  to  4.2  (1.4.22)  (Glen,  1955),  temperature.  from  1.5  (Gerrard  and  the  v a l u e u s u a l l y used f o r  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).  and  The  others,  factor  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) (1.4.23)  0  where  A0  is  a  constant,  R  is  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 f o r c r e e p , and is  the  temperature  polycrystalline approximately  ice  (°K). are  Values  for  p . 3 4 ) . The presence of s m a l l boundary  sliding  of Q f o r secondary c r e e p of  60 k j m o l " 1  139 kJ m o l " 1  (Barnes  for  T > -10°C  amounts  and  of  others,  In  and  (Paterson,  1981,  causes  grain  1971? Jones and B r u n e t , is  dominated  by  1971).  t h i s s t u d y , the i c e i s assumed t o be i s o t h e r m a l at 0 ° C ,  and the flow law parameters P-  T < -10°C,  water  1978) above - 1 0 ° C . The d e f o r m a t i o n below - 8 ° C b a s a l g l i d e (Barnes and o t h e r s ,  T  used are  (Paterson,  1981, Table  3.3,  39) n=3  A = 5.3 1 0 " 1 5  s" 1  kPa"3  (1.4.24)  These v a l u e s a p p l y o n l y f o r secondary c r e e p , a f t e r t r a n s i e n t response  the  initial  t o 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 f o r g l a c i e r ice,  such  as  a  hyperbolic  s i n e r e l a t i o n (Barnes and o t h e r s ,  1971), or a p o l y n o m i a l w i t h odd o r d e r s t r e s s terms ( M e i e r , 1960; L l i b o u t r y ,  I969[a];  Colbeck  l a b o r a t o r y e x p e r i m e n t s (e.g.. G l e n , and  field  1952,  Evans,  1973).  However,  1955; Steinemann,  1958)  measurements of c l o s u r e of b o r e h o l e s and t u n n e l s and  d e f o r m a t i o n of b o r e h o l e s 1953;  and  1958,  (e.cj. G e r r a r d and  Mathews, 1959; M e i e r ,  1963[b]; H a e f e l i , 1971; P a t e r s o n ,  1963[a];  others,  1952;  1960; P a t e r s o n and Savage, Shreve  and  Sharp,  1977) and the flow of i c e s h e l v e s  1970;  Nye,  I963[a], Raymond,  (Thomas,  1973)  30 indicate  that  (1.4.22),  and s a t i s f a c t o r y  1.4.4  known as G l e n ' s flow l a w , i s a u s e f u l  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  ice.  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 t o the bed, and how i t can be i n t e g r a t e d to g i v e the downslope v e l o c i t y component and the i c e and  assumptions  Appendix  flux.  The  are e x p l i c i t l y shown. The d e t a i l s  errors  are g i v e n i n  7.  The s t r e s s e q u i l i b r i u m e q u a t i o n  ( 1 . 4 . 8 ) can  be  integrated  from the s u r f a c e to a h e i g h t z above the g l a c i e r bed t o g i v e shear  stress  c (x,z) xz  which factor  is  = s/>g(h-z)sin«  derived  is  the  parallel-sided  (see  effective of  from  standard slab  bed (e.c[. P a t e r s o n , unity  the  1 +  0  (A7.3.21) formula  deforming  2h  in for  9*'  3*'  dx  dx  xx + h pqha  yy (1.4.25)  Appendix 7. The l e a d i n g shear  stress  in  a  by s i m p l e shear p a r a l l e l t o the  1969, p . 91) when  the  shape  factor  s  is  Appendix 7, S e c t i o n A 7 . 4 ) . The s u r f a c e s l o p e a i s an  s l o p e averaged over a d i s t a n c e of at l e a s t  4h ( B u d d , l 9 6 8 ;  the  order  I 9 7 0 [ a ] ) . The average s t r e s s d e v i a t o r s  i n the  second term a r e d e f i n e d by  31  be'  dx  h(x) d<r' (h-z)l xx d z ' _  =  xx  ice  r  z<h  (1.4.26)  dx  z  0  a corresponding  consider  _  J  = with  1  z=h definition  for  the  yy  t h e x- and y - d i r e c t e d f o r c e s on an i c e column  surface  t o t h e bed, t h e c o r r e c t i o n t e r m s  very  roughly  shear  f o r c e . These r a t i o s a r e u s u a l l y v e r y  and  speaking,  ratios  of the normal  (1967) u s e d a c o r r e c t i o n t e r m s i m i l a r t o longitudinal  strain,  justification formulation  some  model  details,  discussion  of  1970[a]; stress  from t h e  small  for glaciers  ( 1 9 6 1 ) , and R o b i n  this  to  account  choice  of  axes.  the a n a l y t i c a l My  formulation  1.3).  I970[b];  1971)  gave  a  detailed  in  glaciers,  including  variations  t e r m s and t h e w a v e l e n g t h r a n g e s f o r w h i c h t h e y  important.  Hutter  unity  channel  (in press)  gives  t h e most r e c e n t  and  may be  rigorous  of s t r e s s i n g l a c i e r s .  The and  for  (1968) p u b l i s h e d a m a t h e m a t i c a l  correction  treatment  are,  p a r t l y b e c a u s e I use t h e a x e s o f t h e  (see S e c t i o n  (1968;  Budd  I'  f o r c e s t o the basal  (I969[a]) s i m p l i f i e d  i t . Nye  by an a p p r o p r i a t e  in  numerical  of  and C o l l i n s  If  (1.4.25)  in  L l i b o u t r y ( I 9 5 8 [ b ] ) , Shumskiy  i c e sheets.  differs  component.  shape f a c t o r i s an a p p r o x i m a t e c o r r e c t i o n between  zero  f o r the reduction  the  centreline  supported  by  quantitatively parabolic,  the in  when  i n the shear some  valley work  and e l l i p t i c a l  « r  x  z  along  of the weight of the g l a c i e r i s  sidewalls. on  stress  rectilinear  It  was  flow  first  i n rectangular,  c h a n n e l s by Nye ( I 9 6 5 [ c ] ) . I  shape f a c t o r s i n more d e t a i l  i n Appendix  used  describe  7, S e c t i o n A7.4.  32  In  arriving  at  (1.4.24)  i n Appendix 7, I assumed  s l o p e a n g l e s o of the i c e s u r f a c e , small,  that  the  and p of the g l a c i e r bed were  i.e. |o(x)| «  1  ( 1 .4.27)  |f(x) | «  1  (1.4.28)  I a l s o assumed t h a t a was never n e g l i g i b l y s m a l l compared t o T h i s may not be t r u e near an i c e Although  (1.4.8)  divide.  c o n t a i n s s i n $ , the f i n a l  f o r the shear  s t r e s s depends o n l y on s i n e . Nye  pointed  this  out  result.  When  ( 1 . 4 . 2 7 ) and ( 1 . 4 . 2 8 ) h o l d ,  the  the  small  result  bed  slope  dependence of  dependence,  angle  longitudinal  leaving  only  (1.4.25)  (l952[b])  the  first  assumptions  stress  term d * ^ / d x i n ( 1 . 4 . 8 ) i n t r o d u c e s a term i n io-p) the  p.  gradient  which c a n c e l s surface  slope  (1.4.25).  W i t h the assumption t h a t the major shear d e f o r m a t i o n o c c u r s parallel  t o the bed, . i . e . iv  "57/ the shear  strain  /bu bz  «  1  (1.4.29)  rate * = 1 du + ^v dx xz 2 dz  (1.4.30)  i s approximately € xz  =  1 lu 2 bz  which can be i n t e g r a t e d d i r e c t l y , give  (1.4.31) from the bed t o h e i g h t  z,  to  33 z u(x,z)  e  = u (x) + 0  where  u_(x)  dz'  xz  (1.4.32)  i s the b a s a l s l i d i n g v e l o c i t y d i s c u s s e d i n S e c t i o n  1.5. If  I f u r t h e r assume t h a t  invariant  of  the  the square r o o t of  T,  stress  deviator  tensor  a p p r o x i m a t e l y e q u a l t o the shear s t r e s s exz parallel  to  component,  the  bed  is  e  *  , t o get  in  =  e e  (1.4.19)  i_.e. shear  the l a r g e s t s t r e s s  is  stress deviator  << 1  xz  (1.4.33)  xz  G l e n ' s flow law ( 1 . 4 . 2 2 )  (1.4.32),  u s i n g <rxz  the  strain  f o r both T and  p a r a l l e l t o the b e d .  - u (x) =  2A[s(x)Pqsin(g(x))]  n  (n+T5  error  for  from ( 1 . 4 . 2 5 )  the v e l o c i t y component u ( x , z ) u(x,z)  where the  T-  xz  then I can s u b s t i t u t e «xZ  far  second  or  6T  rate  by  ,  the  [h  n+1  - (h-z)  n+1  ] [1 +  e(x)] (1.4.34)  e(x),  i.e.  the  terms  not  included  in  the  computer m o d e l , has the form e(x)  = 0  2h  he'  dx  xx + h  de'  dx  yy  +  6T  (n-1) e  XZ  -  b_v bx/  /fru bz  pgho  (1.4.35) where  the  symbol  0[x]  means " i s of the o r d e r of x " , i . e .  the  f u n c t i o n goes t o z e r o a t the same r a t e as x . I have assumed t h a t  34 the  glacier  coefficient constant.  is A  isothermal, of  Glen's  i . e . temperate,  flow l a w ( 1 . 4 . 2 2 )  so  that  the  can be t r e a t e d as a  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 e v a l u a t e d n u m e r i c a l l y . The downslope i c e f l u x (1.3.5)  f o r use i n the  continuity  equation  is h(x,t)  Q(x,t)  u(x,z,t)  =1  dz  ^0 = u (x,t) s  h(x,t)  + 2A[s/pgsinq] (n+2)  n  [h(x,t)]  n+2  [1+  e(x)] (1.4.36)  With e(x)  the  assumptions d i s c u s s e d above, the e r r o r term i n v o l v i n g  is small; i t is neglected in  present  the  computer  model  in  its  S e c t i o n 1.3  and  form.  The  average  velocity  V(x,t)  used  in  Appendix 1 i s d e f i n e d as V(x,t) = Q(x,t)/h(x,t) which i s V(x,t) = u (x,t) s  + 2A[spgsino] (n+2)  (1.4.37) n  [h(x,t)]  n+1  [1+  e(x)] (1.4.38)  The term on the r i g h t due t o the i n t e r n a l (n+1)/(n+2) surface  times  u(x,h(x),t)  the from  deformation  is  just  downslope v e l o c i t y component a t the i c e (1.4.34).  35 1.5 BASAL SLIDING  1.5.1  INTRODUCTION  The one q u a n t i t y s t i l l of  the  the downslope v e l o c i t y u ( x , z ) and the f l u x Q(x)  s l i d i n g velocity u5(x) in  r e q u i r e d t o complete  sliding  possibly  quantitative In  behaviour,  involved,  pointed  m o d e l l i n g of  out  constant  r e v i e w , gave a  measurements, and the p h y s i c a l  and  Appendix 8,  i s the b a s a l  which appeared as an i n t e g r a t i o n  ( 1 . 4 . 3 2 ) . Raymond ( 1 9 8 0 ) , i n a r e c e n t  of  derivation  some  summary processes  difficulties  of  sliding.  I summarize c u r r e n t  i d e a s on the p h y s 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  this  of s l i d i n g , the way I  section,  I  discuss  t r e a t s l i d i n g i n Chapter  1.5.2 Ice  slide.  In  3, and the reason f o r my c h o i c e .  BASAL ICE TEMPERATURE AND SLIDING  masses  temperatures  the importance  models.  which  are  cold  at  the  base,  have  below the p r e s s u r e m e l t i n g p o i n t , do not appear  The b a s a l i c e i s e f f e c t i v e l y  f r o z e n t o the  glacier  to  bed,  and u (x) s Ice p o i n t at  = 0 (1.5.1)  masses which have t e m p e r a t u r e s at the p r e s s u r e m e l t i n g the i c e - r o c k  interface  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 z e r o t o v a l u e s much g r e a t e r than the to  internal  deformation.  assumes an i s o t h e r m a l  The  model  i c e mass. Due  to  velocities  I describe in this the  existence  of  due study the  36 geothermal  heat  f l u x , and the r e s u l t i n g geothermal  gradient,  the o n l y p o s s i b l e  one  the  at  pressure  essentially  melting  isothermal ice  point  throughout  (neglecting a possibly cold surface layer penetration  temperature  caused  mass  its by  is  volume  diffusive  of the w i n t e r c o l d wave), because o n l y 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 t h a t s l i d i n g v e l o c i t i e s can be an i m p o r t a n t component of motion f o r my m o d e l l i n g s i t u a t i o n s .  1.5.3  PHILOSOPHY OF SLIDING IN THIS STUDY  Correctly sliding  modelling  the  i s , at the p r e s e n t ,  observations,  and  variables possibly I have o u t l i n e d  the  physical  processes  very d i f f i c u l t ,  large  number  due  to  aim,  the  problems  the models p r e s e n t e d  of  measurements,  i n t h i s study,  or t o s i m u l a t e the p h y s i c s of  i n Chapter 3, i s t o i n v e s t i g a t e  (defined  by  the  given that  manner.  I  do  periodic  not attempt  I do not attempt  glacier  sliding.  surging  was  occurs  in  the  to My  surging within a defined  t o 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 of  models.  the consequences of  surging  physical  s t a t e of s l i d i n g  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  glacier,  effects  inadequate  i n v o l v e d i n s l i d i n g p r o c e s s e s . In Appendix 8,  t h e o r y 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  investigate  glacier  of u n c o n t r o l l e d 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 p r o c e s s e s , the p r e s e n t  In  of  also  used  approach  to  investigating  by Campbell and  Rasmussen  (1969) and by C l a r k e ( 1 9 7 6 ) . I could easily  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 ( 1 9 7 6 [ a l ;  Budd (1975) t o 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  would be n u m e r i c a l l y s u s p e c t , due to the problems Appendix 8,  and  consequences of  function u 5 ( x , t ) inferred  S t e e l e ) as  a  add  results  discussed  in  n o t h i n g t o my i n v e s t i g a t i o n of  the  surging.  My approach  well  would  I 9 7 6 [ b ] ) , or  is,  instead,  to  use  a  predetermined  (based as c l o s e l y as p o s s i b l e  on the  s l i d i n g h i s t o r y of a s u r g i n g g l a c i e r driving  function  computer m o d e l . I c a l c u l a t e  for  periodic  sliding reasonably  such as  surges  in  the the  the response of the g l a c i e r model t o  t h i s d r i v i n g f u n c t i o n by u s i n g c o n t i n u i t y and G l e n ' s flow law t o find  the  internal  deformation.  For my purpose of f i n d i n g the  e f f e c t s of s u r g i n g on the i n t e r n a l s t r u c t u r e , no worse than u s i n g a n u m e r i c a l l y inadequate it  this  approach  is  s l i d i n g t h e o r y , and  has the d i s t i n c t advantage t h a t I can c o n t r o l the s l i d i n g at  w i l l while I relate the i c e mass.  s l i d i n g patterns  t o r e s u l t i n g changes w i t h i n  38 CHAPTER 2: MODELS AND TESTS " T h i n g s 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 c h a p t e r , models  and  operation. tests  are  I  I o u t l i n e the o p e r a t i o n  describe  tests  used  to  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 m p o r t a n t . In S e c t i o n 2 . 2 ,  of  the  verify  their  I e x p l a i n why  I describe  it.  In  calculations,  s e c t i o n 2.4  I describe  and i n S e c t i o n 2 . 5 ,  2.1.2  because  the  correctness and  solution  differential  rheologies,  into  the  with  simple  and, o f t e n ,  A f i n i t e difference 1  how  the p a r t i c l e  time  value of  simply  problems the by  I  have  trajectory  boundary  steady  are  solution  most can be  substituting  equation.  a n a l y t i c a l s o l u t i o n s t o i c e flow problems are cases  think  IMPORTANCE OF MODEL TESTING  v e r i f i e d f o r a l l space  few  I  how they were t e s t e d .  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 desirable,  correct  the c o n t i n u i t y  e q u a t i o n g l a c i e r p r o f i l e model, and i n S e c t i o n 2 . 3 , tested  computer  the  Unfortunately,  restricted  conditions,  to  a  uncomplicated  states.  n u m e r i c a l model uses a set  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 .  of  algebraic  39  equations version  whose of  solution  the  closely  truncation  Numerical  solutions  error  (see  this  uncertainty  Appendix 1,  general  temporal v a r i a t i o n s .  increased  digitized  generality,  in  the  boundary  validity  conditions,  is  of  a  the  the  solution  obtained. difference  will  danger  is  that  a  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  of p e r c e n t , may  numerical  behave i n a q u a l i t a t i v e l y r e a s o n a b l e manner,  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 instance,  most,  i n the most r e c e n t time s t e p . The l o n g term  i n t e g r a t e d e r r o r i s unknown. The solution  for  inherent  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  residual errors  ice  new  S u b s t i t u t i o n of a n u m e r i c a l s o l u t i o n i n t o the f i n i t e equations,  problems  The p r i c e which i s p a i d  however,  to  Section A1.5.3).  can extend the domain of s o l v a b l e  t o i n c l u d e those w i t h q u i t e and  a  t r u e 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 ,  within a  rheology,  approximates  yet  wrong. For  may be i n e r r o r by t e n s  and the phase of c y c l i c phenomena  such  as  surging  become t o t a l l y u n r e l a t e d t o the phase of the t r u e s o l u t i o n .  These p o s s i b i l i t i e s  s h o u l d make  predictive  made f o r n u m e r i c a l models. We c o u l d o b t a i n a  claims  r e s u l t no more a c c u r a t e ,  quite  cautious  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  are  two  possible  sources  of  error  s o l u t i o n s . F i r s t , the computer program may not set  of  incorrect  algebraic  constants  undetected. occasion  any  of  behaviour.  There  the  about  at b e s t , than an educated g u e s s , yet be  l u l l e d into b e l i e v i n g that glacier  us  or  equations. missing  in numerical  correctly  Programming e r r o r s , minus  signs,  solve such  sometimes  as go  S p u r i o u s n u m e r i c a l " s o l u t i o n s " of t h i s k i n d have on  found t h e i r way  into  the  scientific  literature.  By  40 careful  program d e s i g n and the use of c o n s i s t e n c y  e r r o r s can be e l i m i n a t e d , taken the time t o do Assuming  that  there  is s t i l l  another  used  in  computer  differential the  not  all  programmers  the computer  source of model  program works  error.  The  correctly,  numerical  computational  may i n c l u d e an e x p o n e n t i a l l y unrelated  instability  (Appendix  growing h i g h wavenumber  t o the d i f f e r e n t i a l  i s u s u a l l y easy t o  still of  solution  may  equation.  (Appendix  factors.  1,  oscillation Fortunately,  the  possibility  that  d r i f t away from the t r u e s o l u t i o n ,  For example,  1, S e c t i o n A 1 . 5 . 2 )  is  equations  l o o k " 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  several  the  recognize!  The o t h e r cause f o r concern i s numerical  scheme  may not a d e q u a t e l y r e p r e s e n t  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  this error  have  e q u a t i o n at a l l time s c a l e s . -One example of t h i s  linear  completely  these  so.  then,  the  although  checks,  inappropriate  mesh  the yet any  intervals  may cause i n c o r r e c t d i s p e r s i o n  and  spectral  attenuation,  d i s t o r t i n g the s o l u t i o n t o an u n a c c e p t a b l e  degree.  Introducing  smoothing  numerical i n s t a b i l i t i e s still  fail  to  (unpublished),  can cause s i m i l a r d i s t o r t i o n s  t o t a l l y remove the i n s t a b i l i t i e s ) .  the w i d e l y r e p o r t e d and J e n s s e n ,  schemes i n a t t e m p t s t o suppress  1975; Budd and Radok,  1971)  was  who appears t o have encountered In an attempt t o model  presence  growing  reported:  For  may  example,  i c e sheet model of the Melbourne group (Budd  difficulties. of  (and  numerical  surging  instability,  used  by  Mclnnes  a l l of the above glaciers Mclnnes  in  the  (p.  64)  41 " D i f f e r e n t time s t e p s g i v e s l i g h t l y d i f f e r e n t surge times and p e r i o d s and t h e r e f o r e at the same growth t i m e , p r o f i l e s are not d i r e c t l y comparable. A l s o due t o the d i f f e r e n t number of steps, and the smoothing scheme not being p e r f e c t , the lower the time s t e p the more i t e r a t i o n s , which l e a d s t o more smoothing which tends t o lower e i t h e r the depth or the base v e l o c i t y , and t h e r e f o r e affects the surging times." and  further,  discussing  a  criterion  to  introduce  smoothing at the appearance of i n s t a b i l i t y problems  automatic  (p.  64):  " U s i n g t h i s t e s t before a u t o m a t i c smoothing, lessens the number of times smoothing i s u s e d , and t h e r e f o r e d e c r e a s e s the e f f e c t smoothing has on the e x a c t 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 f o r optimism i n the that  this  numerical  model,  for  instance,  was  s o l u t i o n t h a t c l o s e l y matched the t r u e s o l u t i o n The  results  providing  at  all  a  times.  of the Mclnnes study were p u b l i s h e d by Budd (1975)  and by Budd and Mclnnes example,  belief  not  to  (1974;  criticize  1978;  1979).  I  mention  this  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 a c k of a t t e n t i o n p a i d t o t h i s  serious  question  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  glacier  difficulties (see  modelling  program  with accuracy,  Appendix 16),  p r i o r i t y discussion  yet  consistency,  model  it  is  apparently and  mass  had s e r i o u s conservation  v e r i f i c a t i o n has not been g i v e n  i n the p u b l i s h e d  Because n u m e r i c a l model r e s u l t s complicated  has  literature. cannot be v e r i f i e d f o r  problems the n u m e r i c a l models are c r e a t e d t o  the  solve,  i m p e r a t i v e t h a t n u m e r i c a l schemes be v e r i f i e d by comparing  t h e i r numerical solutions with a n a l y t i c a l solutions problems  before  for  the n u m e r i c a l models are used f o r new p r o b l e m s .  There i s always a t e m p t a t i o n w i t h a new model t o r u s h solution  of  simpler  complicated  glaciological  into  the  p r o b l e m s . T h i s urge to  42 break new ground  in  numerical  has  model  a  hurry been  must  be  controlled  demonstrated  the  t o work a c c u r a t e l y  known ground. Only i n t h i s way can t h e r e be the  until  any  confidence  on in  results.  2.2 CONTINUITY EQUATION PROFILE MODEL  I  2.2.1  INTRODUCTION  have  written  numerically provided  solves  with  a  a  FORTRAN  the  continuity  subroutine  to  on  Glen's  flow  equation  which  (1.3.5),  when  the v e l o c i t y V ( x , t )  The flow  (1.4.22)  program  is  equation used  (1.4.38)  for  (Chapter  detail  Appendix 1. I summarize the c o m p u t a t i o n a l a s p e c t s i n  this Section, the  boundary  (Section  i n c l u d i n g the conditions  Q0(t) the  numerical (Section  2 . 2 . 4 ) , and a c c u r a c y  Inputs to balance  I d e s c r i b e the c o m p u t a t i o n a l  glacier  simulations in  3).  law  computer  calculate  averaged through the i c e t h i c k n e s s . based  IV  A(x,t),  the the  model  initial  (Section  2 . 2 . 3 ) , numerical  (Section  are  scheme  scheme  2.2.5).  bedrock  2.2.2), stability  . •  elevation  b(x),  mass  i c e t h i c k n e s s H 0 ( x ) , the i c e  through the upslope boundary, and s l i d i n g and flow p r o p e r t i e s  in  parameters  to  flux  specify  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  profile.  43 2.2.2  THE NUMERICAL SCHEME  Complete d e t a i l s  of  the  numerical  scheme  are  given  in  Appendix 1, S e c t i o n s A1.1 and A 1 . 2 . I  solve  the  continuity  equation  (1.3.5)  d i f f e r e n c e method (e.g_. R i c h t m y e r and M o r t o n , differential algebraic mesh  equations  points  an i n i t i a l time  equation  (1.3.5)  is  by  1967).  approximated  f o r the i c e t h i c k n e s s  a  The p a r t i a l  by  a  condition { h ^ l j ^ j j } , using  a  the  possibly  f i n i t e d i f f e r e n c e equations  of  ( h - | j = 1 , J } at a set  solution  is  variable  of  obtained  are n n (Q - Q ) j + 1/2 j-1/2 Ax j  n+1 n = 6A + (1-6)A j j  (2.2.1)  1 £ j < J superscripts the  indicate  spatial  t h i c k n e s s h- are  mesh  measured  by  time s t e p A t . The  n+1 n n+1 n+1 h - h + 6 (Q - Q ) + (1-6) _j i " j + 1/2 j-1/2 At W Ax W j j j  indicate  set  at equal h o r i z o n t a l i n t e r v a l s of A x . S t a r t i n g from  marching,  where  finite  the  time  point.  normal  1 < n < N step,  Mass  to  the  and  subscripts  balance bed,  A- and  and  the  ice mesh  w  increments  AXj  i s a constant scheme.  are measured a l o n g the bed. The weight f a c t o r 6  between z e r o and  unity,  I discuss 6 further in Section  used  to  t h i c k n e s s hj at the m e s h p o i n t s ,  the v e r t i c a l l y averaged downslope meshpoints by  the  2.2.4.  The i c e f l u x Q j +1/2 between the meshpoints ice  stabilize  i s r e l a t e d t o the  the c h a n n e l w i d t h  velocity  v  j • 1/2  Wj /z +1  between  and the  44 Q  j±l/2  = V  Most  of  the  flux  is  calculated  J±1/2  variables  W  J±1/2  (h  j±1 2  + h ) i  (2.2.2)  i n ( 2 . 2 . 2 ) are shown i n F i g u r e 2 . 1 . The  midway  between  mesh  points  because  of  numerical 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 e q u a t i o n s  ( 2 . 2 . 1 ) can be w r i t t e n i n m a t r i x form  as M H = D where the components  of  the  at  are  {h-* 1 | j = 1,J}, J  the r i g h t s i d e v e c t o r D c o n t a i n s time s t e p ,  points  H  at  previous  mesh  vector  thicknesses  the  the  (2.2.3)  the  related  makes  to  (2.2.3)  iteratively.  the  ice thickness  nonlinear, For  the  future  step  quantities  from  and  first  Since  (Sections  the  system  estimate  {0VP* 1 | j = 3/2, of  the  ice  used  5/2,...J+1/2} thickness  ice  t o c a l c u l a t e an e s t i m a t e  (2.2.4)  the  velocity  1.3 and 1 . 4 ) , must  be  this solved  i s used  as  , t o c a l c u l a t e the  {0h?* 1 | j = 1 , J } thickness  .  an first  Prescripts  estimates  are  { , V p + 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 residuals  coefficients  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  i n d i c a t e i t e r a t i o n number. These then  ice  time  the p r e v i o u s time s t e p {V? | j = 3 / 2 , 5 / 2 , . . . J + 1 / 2 } e s t i m a t e of  unknown  and the m a t r i x M c o n t a i n s  i n v o l v i n g v e l o c i t y at the f u t u r e time s t e p . is  the  the  future  step.  The  45  r  n+1  j  =  n+1 n+1 n n 2p 6[Q - Q ] + 2p ( 1 - 6 ) [ Q - Q ] " 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 then  measure  the  j  =  At/(2AxW ) j  (2.2.5)  degree t o which the c u r r e n t e s t i m a t e s of  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  continuity  ice  equation  (1.3.5). By  using a l i n e a r i z e d equation  (2.2.6) to relate  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 t o make the  residuals residuals  go t o z e r o , r  J \-~\ br** = > —j 6h *—* d h *+ 1 k k=1 k  n+1  1  j  A 6h = r  or I  obtain  (2.2.6)  essentially  a  (2.2.7)  multi-dimensional  f o r m u l a t i o n of  Newton-Raphson method (e.cj. Carnahan and o t h e r s , to  solve  residual in criterion.  (2.2.4). absolute  The  iterations  value  is  1969,  p.  the 319)  t e r m i n a t e when the l a r g e s t  smaller  than  a  preset  test  I d i s c u s s the c h o i c e of a c r i t e r i o n i n Appendix 11.  46  2.2.3 The  set  equations (1.3.5)  in  BOUNDARY CONDITIONS of  equations  J+2  (2.2.1) with (2.2.2) consists  unknowns  { h?* 1 | j = 0,J+1}.  i s a f i r s t order d i f f e r e n t i a l equation,  boundary c o n d i t i o n s e l e c t e d  point  h0  equation  i t requires  one  equations.  At an i c e d i v i d e , z e r o i n p u t f l u x image  J  from ( 1 . 3 . 1 0 ) through ( 1 . 3 . 1 2 ) . T h i s  g i v e s one of the r e q u i r e d two e x t r a  an  Since  of  at a d i s t a n c e  i s m o d e l l e d by  including  Ax o u t s i d e the boundary  j=1,  with h  = h O  2  V 1  This forces the  ice  -v  =  '2  3  /2  (2.2.8)  the s u r f a c e s l o p e t o be z e r o at  thickness  can  i l l u s t r a t e d in Figure  vary  with  the  time.  boundary,  This  situation  h  slope,  = 0 (2.2.9)  1  and the i c e s u r f a c e s l o p e can a d j u s t t o any a p p r o p r i a t e To model o n l y a p o r t i o n end  of  of  an  ice  mass  (e.£.  Chapter 4 ,  where  the i c e f l u x must be s p e c i f i e d a t  above the boundary,  such  value. that  the  the model i s some d i s t a n c e downslope from the  b e r g s c h r u n d or d i v i d e modelled),  is  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  upstream  but  an  icefall  is  1/2 mesh increment  i..e. Q  n  = Q (nAt)  (2.2.10) The second e x t r a e q u a t i o n r e q u i r e d t o b a l a n c e the number of 1 / 2  0  47  FIGURE 2 . 1 . N u m e r i c a l Scheme At Ice  Divide.  e q u a t i o n s 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 s p a c e , and  simply  a l l o w i c e t o flow through i t and out of the model, I can get  the  equation  the  by w r i t i n g h j + 1 as an a p p r o p r i a t e  i c e t h i c k n e s s at the m e s h p o i n t s . divided differences p.  e x t r a p o l a t i o n of  I use the second o r d e r Newton's  p o l y n o m i a l (e.g.. Carnahan and o t h e r s ,  1969,  12). If  I  choose  to  follow  the a c t u a l motion of the  s n o u t , d e f i n e d by x = L ( t ) , I must keep falls  between  the  meshpoints,  track  of  L(t)  glacier when  and be a b l e t o add or  subtract  p o i n t s t o or from the mesh as the t e r m i n u s moves. I assume the  terminus  snout at L ( t ) , principle  of  it  that  i s wedge shaped from the l a s t mesh p o i n t J t o the as  shown  in  Figure 2.2.  Then,  I  apply  c o n s e r v a t i o n of mass t o the shaded s e c t i o n of  wedge, o b t a i n i n g the r e q u i r e d f i n a l e q u a t i o n  by  balancing  the the the  48  4.—  FIGURE 2 . 2 .  volume  AX/2—H  The Wedge T e r m i n u s .  change of the wedge a g a i n s t i c e f l u x i n t o the wedge from  upslope p l u s net mass balance time  step.  The  details  on the upper s u r f a c e ,  may  be  A 1 . 3 . 4 . As the terminus moves, upstream  found i n Appendix 1,  the  ice  thickness  terminus  model  Bindschadler  is  similar  in  h^. at  many  the  flow.  ice This  r e s p e c t s t o t h a t used by  NUMERICAL STABILITY  the f i n i t e d i f f e r e n c e  numerical  scheme  and  the  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  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 differential  in  the  (unpublished).  2.2.4  increments  changes  each  Section  end, the snout p o s i t i o n L ( t ) , and the s l o p e of the  s u r f a c e are a l l f r e e t o a d j u s t t o  If  during  equation.  This  from the  usually  e x p o n e n t i a l growth of some wavenumber  solution  involves  component  a  which  mesh  equations of  the  spurious quickly  49 dominates the d e s i r e d bounded, Rigorous usually  stability  not f e a s i b l e ,  analogues stability The  analysis  but  generally  useful  of the n o n l i n e a r  forms.  first  type  of  At and A x . For  equations,  instability  large r e l a t i v e the  nonlinear criteria  guidelines  stability  instability,  increments  of  stability  give  computational  than  p h y s i c a l l y reasonable  involves many  can  systems  of  V,  material  the  finite  per time s t e p , and the system tends  physical  solution.  to  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  n,  to  i s one  standard  stability  i n the wavenumber  the  for  finding  the  domain which m u l t i p l i e s  of the s o l u t i o n at  time the  step future  If (2.2.11)  at a l l wavenumbers m, no wavenumber exist.  mesh  analysis  T(m) < 1  can  greater  many  of the s o l u t i o n at the p r e v i o u s  g i v e the F o u r i e r t r a n s f o r m  t ime s t e p n+1.  too  von Neumann method (e.g^. Richtmyer and  linear  the F o u r i e r t r a n s f o r m  difference  'forget'  1967,  f u n c t i o n T(m)  the mesh  much  Morton,  transfer  into  linear  of  travels  intervals  p. 70),  and i n s i g h t s  is  is  linearized  occur when the time s t e p At i s  velocity  The  for  choice  to the space s t e p A x . I f A x / A t  material  equations  problem, the  solution.  can grow, so no  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  to a c o n s t a n t i n ( 1 . 3 . 7 ) ,  I f i n d t h a t by  instability  the v e l o c i t y V  choosing  0^1/2  (2.2.12)  I obtain unconditional l i n e a r computational  stability  for  any  is called  'the  c h o i c e of Ax and A t . The nonlinear  second  type  instability'  of  numerical  (e.£.  instability  Phillips,  1959;  Mesinger  and  50 Arakawa,  1976,  p. 35).  This  i s a problem which a r i s e s i n any  n u m e r i c a l s o l u t i o n of a d i f f e r e n t i a l e q u a t i o n u s i n g mesh,  and  a  discrete  h a v i n g terms which are n o n l i n e a r i n some c o m b i n a t i o n  of the dependent and  the  independent  variables.  In  (1.3.5),  dQ/dx has t h i s form. The  nonlinearity  wavenumber spectrum) wavenumber  part  pumps  the  wavenumber  v a r i a b l e , and the a l i a s i n g  distorts  the  importance modelling  back  is  to  spectrum of the  the  the  lower  wavenumbers,  the high  dependent sampling where  it  S i n c e the n o n l i n e a r i n s t a b i l i t y has an  not  community,  end  (Appendix 3) due t o d i s c r e t e  to  solution.  which  (squared a m p l i t u d e of  from the low wavenumber  of  f o l d s t h i s energy  energy  I  widely  recognized  discuss  it  in  in  the  glacier  d e t a i l i n Appendix 1,  Sect i o n A 1 . 4 . 3 . I f a f u n c t i o n i s known o n l y at wellknown fact  that  result  from  i t s Fourier  discrete  sampling t h e o r y (see spectrum  can  be  intervals  Ax,  Appendix 3) i s  found  only  up  a the  to  a  wavenumber m^, c a l l e d the N y q u i s t wavenumber, m  =  N This  is  a  rr Ax  sampling r a t e 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 Figure A3.2),  by  being  s y m m e t r i c a l l y about m N (see It  'folded'  two  band-limited  (A1.4.7) Fourier  product b a n d l i m i t e d t o the sum of signals.  This  happens  lower back  wavenumbers into  the  (see  spectrum  Figure A1.3).  i s easy t o show (see  multiplying  (2.2.13)  the  through  (A1.4.9))  s e r i e s together bandwidths  of  that  gives a the  two  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  is bandlimited  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^ t o m^, but  the  Figure  lower  2/3  of  the  spectrum  remains  correct  (see  A1.3).  It  is  evident,  instability, 2/3m^,  be  heavily  i n f o r m a t i o n about the and  they a t t r i b u t e d actually  to  attenuated.  must not d i s t o r t  Budd  that  avoid  the  the  i.e.  same  the low wavenumbers  (1975) encountered  machine  roundoff  time,  which  spatially.  p.  102)  computer  used  the  same  s u r g i n g at B r u a r j o k u l l , I c e l a n d .  contain  I  think  it  from  Mclnnes,  velocity  steady  profiles.  Mclnnes  (unpublished,  programs  to  simulate  The broken c u r v e s i n F i g u r e  (unpublished),  state  with  no  p . 58)  I discuss  show  preferentially  use  of  damp h i g h wavenumbers.  (1.4.34))  glacier  the  the  velocity  these  if  term  methods i s hard t o j u s t i f y .  the  depends  flow on  the  equation local  to  the  (1.3.5)  to  I conclude t h a t the use  i n t h i s n u m e r i c a l model are s u p e r i o r First,  the  method used by Budd and Jenssen ( 1 9 7 5 ) , and a l s o ,  a d d i t i o n of a p u r e l y n u m e r i c a l d i s s i p a t i o n  either  2.3  s l i d i n g , and no smoothing of  In Appendix 1 S e c t i o n A 1 . 4 . 4 ,  smoothing  was  profile  i n s t a b i l i t y which arose as he attempted t o b u i l d up the a  the  n o n l i n e a r i n s t a b i l i t y . Budd and J e n s s e n attempted  i t began to o s c i l l a t e  to  above  an i n s t a b i l i t y which  error.  whenever  (redrawn  nonlinear  glacier.  Jenssen to  At  t o cope w i t h the i n s t a b i l i t y by smoothing the  58; p .  the  the a l i a s e d s i g n a l at h i g h wavenumbers,  must  attenuation  then,  for  The two methods I  on p h y s i c a l the  grounds.  continuum  (e.g..  ice surface slope o(x),  l a r g e a m p l i t u d e bumps i n the s o l u t i o n  profile  of  should  tend  then to  52  Vatnajokull model Plots at 50 years 1500 200 years 250 years 1000 o UJ X  500 10  20  DISTANCE  AO  30  (km)  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 c u r v e s redrawn from Mclnnes ( u n p u b l i s h e d , 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 f o r 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 m i d p o i n t s of the mesh i n t e r v a l s f o r 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 c u r v e s are the 50 year p r o f i l e s from my computer model using parameters i n T a b l e 2 . 1 . The Budd-Mclnnes model a l s o appears t o c r e a t e mass (see Appendix 16).  diffuse  out  medium (See  rapidly,  to  the  p h y s i c a l p r o p e r t i e s of  Appendix 6, which r e l a t e s p e r t u r b a t i o n s  to a d i f f u s i o n should  due  also  instabilities.  equation  of  the  (1.3.5)  f o l l o w i n g Nye ( 1 9 6 0 ) ) . T h i s same p r o c e s s  efficiently  smooth  out  high  In Appendix 13, I have shown t h a t , when  wavenumber the  ice  53 flux  is  . i . e . at process  calculated J±1/2, is  as  at  the  midpoints  shown i n F i g u r e 2 . 2 ,  incorporated  this,  intervals,  physical  diffusion  i n t o the n u m e r i c a 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 For example,  of the mesh  wavenumbers.  the broken c u r v e s  i n F i g u r e 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) model (Budd,  1975; M c l n n e s , u n p u b l i s h e d ) .  The  original  on F i g u r e 4.3 of M c l n n e s , from which these p r o f i l e s  are  computer caption redrawn,  was "Profiles from the V a t n a j o k u l l model at fifty year i n t e r v a l s , showing the i n c r e a s e i n the magnitude of the oscillations with time, due t o the two p o i n t finite d i f f e r e n c e a p p r o x i m a t i o n . In t h i s c a s e , no smoothing was used." These a u t h o r s d i d not c a l c u l a t e their  mesh  intervals.  the i c e f l u x at the m i d p o i n t s of  As a r e s u l t ,  the d i f f u s i v e  G l e n ' s flow law was unable t o p r e v e n t Nyquist  wavelength  instability  could  (in have  this  serious  mechanism of  aliasing  c a s e , 2 km). The t r i g g e r  been  a  large  truncation  at  the  for  the  error  (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  difference  (Budd and J e n s s e n ,  n bar 2  at  A -n  .225  TABLE 2 . 1 .  a  -1  the boundary  s  1.0  g  Ax  p  ms" 2  kg m " 3  9.8  910.  m  At a  1000.  Parameters f o r V a t n a j o k u l l ( F i g u r e  c u r v e s are the 50 year p r o f i l e s  1975). The s o l i d  1.0  2.3)  u s i n g my computer model w i t h  parameters i n T a b l e 2 . 1 . As f a r as  I can t e l l ,  the  Mclnnes a l s o used  54 these v a l u e s .  The  wavelengths  due  nonlinear to  instability  is  removed  at  all  the c h o i c e of n u m e r i c a l scheme. As w e l l  as  the h i g h wavenumber o s c i l l a t i o n , the Budd-Mclnnes model  appears  to suffer  further  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  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 . stress  may  be  For example,  the  c a l c u l a t e d u s i n g an i n t e r m e d i a t e or l a r g e s c a l e  s l o p e (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 . 9 2 ) . remove  the  nonlinear  instability  without  wavenumber s i g n a l at a l l , b y , at the step,  gravitational  In t h i s  case,  I  d i s t o r t i n g the low  completion  of  each  time  t a k i n g the F o u r i e r t r a n s f o r m 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  filter  a t 2/3m^. I then  in Figure 2.4;  t h i s f i l t e r has a c u t o f f  p e r f o r m the i n v e r s e F o u r i e r t r a n s f o r m  to  obtain  which are b a n d l i m i t e d a t 2/3 m^, and u n a f f e c t e d nonlinear  instability  cannot  grow,  and  s o l u t i o n below 2 / 3 m v . By a s u i t a b l y s m a l l  the  the  lowpass  profiles  by a l i a s i n g .  cannot choice  affect of  the  The the mesh  increment A x , m.. can be made l a r g e enough so t h a t a l l p h y s i c a l l y interesting the c u t o f f  wavenumbers  ( 1 9 5 9 ) , who o r i g i n a l l y  instability,  suppressed  procedure was  quite  1965).  the g l a c i e r p r o f i l e are w e l l below  wavenumber.  Phillips  predated  in  it  costly  in to  identified this  manner.  implement,  the  nonlinear  However,  because  his  the work  the F a s t F o u r i e r Transform a l g o r i t h m (Cooley and Tukey, In my work, u s i n g the f i l t e r  F o u r i e r Transform  method  did  not  of F i g u r e 2.4 w i t h the F a s t dramatically  computer e x e c u t i o n time f o r the t e s t s I p e r f o r m e d .  increase  the  55  IT(m)l  m  AX  TT/2  2TT/3  FIGURE 2.4. F i l t e r To Suppress The Nonlinear The Nyquist wavenumber i s at mAx=ir.  Instability,  2.2.5 ACCURACY Having accurately Section  achieved it  a  solves  stable  the  scheme,  partial  I  must  differential  now  ask how  equation.  In  A1.5, I examine the accuracy of the numerical scheme by  two somewhat complementary methods. The  first  method f o l l o w s from  the  von  Neumann  stability  analysis  i n S e c t i o n 2.2.4. At a l l wavenumbers m, I compare both  amplitude  and phase of the t r a n s f e r f u n c t i o n  numerical  scheme with those of the p a r t i a l  Errors  i n amplitude  T(m)  of  the the  differential  eqation.  represent i n c o r r e c t a t t e n u a t i o n , and e r r o r s  i n phase r e p r e s e n t i n c o r r e c t propagation speeds and For the l i n e a r  e q u a t i o n , the c o n c l u s i o n i s t o use 6=1/2  for  dispersion.  (2.2.14)  the most a c c u r a t e amplitudes  Ax as s m a l l as f e a s i b l y  possible  (see F i g u r e A1.7), and to take to  get  accurate  phase  (see  56 F i g u r e 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 f o r the n o n l i n e a r m o d e l . The second method error'.  is  an  estimation  of  of the f i n i t e d i f f e r e n c e e q u a t i o n s . estimate,  solution  h(x,t)  {hM|j=l,J}  I t i s a s p a t i a l domain e r r o r  i . e . i t e s t i m a t e s the t o t a l e r r o r at each x p o s i t i o n ,  r a t h e r than e s t i m a t i n g the e r r o r i n each  sine  A f t e r assuming t h a t i c e t h i c k n e s s h ( x , t ) infinitely differentiable, expressed  as  truncated  and Q ( x , t )  expressed result  'truncation  T h i s i s the d i f f e r e n c e between an e x a c t 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 , and an e x a c t  h(x,t)  the  in  Taylor  of  the  expansions  are  s o l u t i o n can be  about the s o l u t i o n  equation.  first  component.  and i c e f l u x Q ( x , t )  the f i n i t e d i f f e r e n c e  of the d i f f e r e n t i a l  terms  wave  error  is  neglected d e r i v a t i v e s .  The  i s t h a t the t r u n c a t i o n e r r o r  The  has the form  J €  n j  =  At  (1-26) * 2h 2  at  2  1  + At2 — 6  a3h n at  3  j  + Ax2  1  — 24  b 3Q n  dx3 j (2.2.15)  The l e a d i n g term v a n i s h e s w i t h the c h o i c e 6 = 1 / 2  T h i s i s the same r e s u l t o b t a i n e d from (2.2.14).  the  wavenumber  analysis  The t r u n c a t i o n e r r o r i s a l s o m i n i m i z e d by keeping the  mesh increments as s m a l l as p o s s i b l e . estimated  (2.2.16)  The  coefficients  can  be  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 n u m e r i c a l  model t o get a q u a n t i t a t i v e e s t i m a t e  of the e r r o r .  57 2.3 TESTING THE CONTINUITY MODEL  2.3.1 Two  INTRODUCTION  tests  verifies  are  described  in  t h a t the model s a t i s f i e s  width,  balance,  this  continuity  programs a c c u r a t e l y  second  test  other  The  first  nonconstant and t h a t  verifies  the  that  the  s o l v e a r e a l i s t i c n o n l i n e a r problem i n which  the flow v e l o c i t y depends on both h ( x , t ) tests  of  and i t s g r a d i e n t  numerical  models are  ( e . £ . Waddington, 1979), the t e s t s p r e s e n t e d reasonably  with  bed, v e l o c i t y , and i c e t h i c k n e s s ,  t e r m i n u s moves c o r r e c t l y . The  Although  subchapter.  and s t r i n g e n t  available  here are a  simple,  trial  of any n u m e r i c a l  f l o w l i n e model based on the c o n t i n u i t y e q u a t i o n  (1.3.5). I think  that  comprehensive,  in x.  i t i s a reasonable proposal  minimum  standard  for  t h a t these t e s t s be used  verification  n u m e r i c a l models of g l a c i e r  and  comparison  flow p r i o r t o a t t e m p t s t o  as of  a all  use  them  to  move  t o model c o m p l i c a t e d i c e masses.  2.3.2 The  CONTINUITY TEST WITH TERMINUS MOTION  ability  of  the  glacier  correctly,  i . . e . at a r a t e c o n s i s t e n t  used,  critical  is  Nye(l963[a], and  model with  the  flow  law  t o a c c u r a t e s i m u l a t i o n of g l a c i e r  I 9 6 3 [ b ] ) ; Nye i n d i s c u s s i o n  of a paper  being  flow ( e . £ . by  Mahaffy  Andrews ( 1 9 7 6 ) ) . I f the model t e r m i n u s advances too s l o w l y ,  i t w i l l a c t as a dam, r e s u l t i n g i n a g l a c i e r too  terminus  t h i c k , s l o w , and s h o r t ,  everywhere.  If  solution  even though c o n t i n u i t y i s  which  is  satisfied  the t e r m i n u s of the model advances more r a p i d l y  58 •  h  s  0  w  w  0.01 -0.01  1 .0 1.0  1 .0 1 .0  0  -1 .0 -0.1  0.1 0.1  s  0  c  At  AX  -0.02 -0.02  0.5 0.5  0.01 0.01  TABLE 2 . 2 . Parameters f o r 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 i v e s the advancing model, and the second line gives a retreating model. See F i g u r e 2.5 and Figure 2.6. than the c o r r e c t  r a t e , the g l a c i e r  f a s t f l o w i n g , and No  simple  terminus. glacier  Nye  (1.4.8) (1967)  terminus.  horizontal  bed,  are  through  published  available (1.4.10)  and  perfectly  that  plastic  in Section 2.2.3 the  the  choose  its  own  steady  ice,  state,  i t i s not  terminus  suitable  model  continuity  equation  (A1.1.1)  this solution accurately. Other s o l u t i o n s below  is  a  could easily  So* s,  of  solution  to  the  the model reproduces  i s t o some e x t e n t However,  arbitrary.  the  one  used  good one because i t i s e v i d e n t at a g l a n c e whether  the model r e s u l t s are c o r r e c t , number  analytical  be f o u n d .  allowed  I use a c h a n n e l , mass  and show t h a t  My c h o i c e  being  model s a t i s f i e s c o n t i n u i t y  position,  b a l a n c e , and a flow law g i v i n g an  a  and A 1 . 3 .  numerical  terminus  the  glacier  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 to  thin,  to s i m p l i f y  near  S i n c e t h i s s o l u t i o n assumed  tested i s described check  too  a s o l u t i o n f o r the shape of a  f o r i n c l u s i o n i n my n u m e r i c a l m o d e l . The  To  be  ex tended. approximations  s t r e s s equations  solution w i l l  nonconstant  and i t t e s t s  and n o n t r i v i a l  the  model  input f u n c t i o n s .  with  a  Let h 0 ,  W 0 , W , and c be c o n s t a n t s , l e t the c h a n n e l w i d t h W(x) be  59  Distance  FIGURE 2 . 5 . C o n t i n u i t y Test w i t h Moving T e r m i n u s . Flow i s t o 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 t o 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 e q u a t i o n s ( 2 . 3 . 1 ) through ( 2 . 3 . 5 ) w i t h the c o n s t a n t s i n Table 2 . 2 .  W(x) = W +" W x o and the mass balance A(x)  (2.3.1)  be  A(x) = sx + c/W(x) I use a v e l o c i t y  (averaged over depth) g i v e n  has no p h y s i c a l s i g n i f i c a n c e ;  (2.3.2) by  (2.3.3).  i t i s merely a n u m e r i c a l  test.  This  60  V(x,t) = cx/[W(x)(h(x,t)-h  )] x*0  0  = c / [ s ( t )W.(x) ]  x=0 (2.3.3)  where the s u r f a c e s l o p e s i s g i v e n by s(t)  = s  + s t (2.3.4)  0  (During  the  numerical s o l u t i o n procedure,  c a r e s h o u l d be taken  t o 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 ) this  can  be done by a p p r o x i m a t i n g  analytical  function  h(x,t)>h0.) (A1.1.1)  in  the  Substitution  ( 2 . 3 . 3 ) by a s u i t a b l y  region of  approaches  of  h(x,t)  (2.3.1)  near  through  = h  + s(t)  smooth h0  x  is  x (2.3.5)  measured a l o n g the g l a c i e r  bed. Any bed p r o f i l e  may be used i n the n u m e r i c a l m o d e l . The t h i c k n e s s h 0 at constant times, This s.  f o r a l l t i m e , and h ( x , t )  varies  with  slope  the c o r r e c t g l a c i e r shows  u s i n g the c o n s t a n t s i n the f i r s t  line  rate  h(L(t),t)=0,  t e r m i n u s p o s i t i o n i s found t o be at L(t)=-h/s(t) o the n u m e r i c a l r e s u l t s  is  s(t).  changes l i n e a r l y w i t h time at the c o n s t a n t  By s o l v i n g ( 2 . 3 . 5 ) f o r the v a l u e of L such t h a t  F i g u r e 2.5  x=0  l i n e a r l y w i t h x at a l l  r e s u l t i n g i n a wedge-shaped " g l a c i e r " slope s(t)  into  is  0  Distance  and  (2.3.5)  v e r i f i e s t h a t the t h i c k n e s s s o l u t i o n h ( x , t ) h(x,t)  h0;  (2.3.6)  f o r the a d v a n c i n g model 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  right  and  is  shown  at  equal  n o n d i m e n s i o n a l i z e d u n i t s f o r the f i r s t  time  i n t e r v a l s of  40 u n i t s , and at 5  10  units  61  FIGURE 2 . 6 . C o n t i n u i t y Test With Moving T e r m i n u s . Results i n F i g u r e 2.5 w i t h bed e l e v a t i o n d i s t a n c e measured a l o n g the b e d .  removed  thereafter.  w i t h the bed  F i g u r e 2.6  elevation subtracted, bed  rather  h(x,t) part  than  shows  and d i s t a n c e  103.  but  results  x measured  in  the  correctly  under  duplicates  the  curves  reverse order. This v e r i f i e s  quite  model  one  i s o b t a i n e d w i t h the  n u m e r i c a l model s a t i s f i e s c o n t i n u i t y everywhere terminus  the  (2.3.5) to w i t h i n  The same degree of a c c u r a c y which  along  i t s h o u l d be the s o l u t i o n  i t reproduces  r e t r e a t i n g model (Table 2.2) Figure 2.6,  same  horizontally, i.e.  i n ( 2 . 3 . 5 ) . In f a c t , in  the  and  and  that  moves  of the the  g e n e r a l c o n d i t i o n s ; the w i d t h  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  velocity  v a r i e s w i t h h , x , and t .  2.3.3 The  CONTINUITY TEST WITH BURGERS' EQUATION  test  in  Section  2.3.2  verifies  c o r r e c t l y with general geometrical In  this  section  I  show  on both  equation ice  The  theory  diffusion  and  by c o r r e c t l y  solving  i n the  a  fully  form of v e l o c i t y depending  slope.  The  problem  solved  also  waves. of  propagation  has been i n v e s t i g a t e d  the l i t e r a t u r e of gas Rankine  terminus.  the i t e r a t i v e procedure  with a r e a l i s t i c  thickness  includes kinematic  i n p u t and a moving  that  n u m e r i c a l model works a c c u r a t e l y nonlinear  t h a t the model works  of  shock waves i n a gas  by many a u t h o r s .  dynamics  have  included  C o n t r i b u t o r s to Stokes  (1848),  (1870) and T a y l o r ( 1 9 1 0 ) . Methods d i s c u s s e d i n the  on n o n l i n e a r waves by Whitham (1974) are c l o s e l y f o l l o w e d One  standard  approach  space and time i n equation  a  shock  analogous  to  r i g h t hand s i d e ) t o g e t h e r law  relating  flux relations  gas  to  find  here.  is  to  solve  a  continuity  (but w i t h no source term on the  w i t h an e q u a t i o n analogous t o  f l u x Q t o gas  d e n s i t y H . One of the  a  thickness,  is  diffusion  term g i v e s the and  flow  simplest  Q = Q (H) - v oH  The f i r s t  text  the gas d e n s i t y as a f u n c t i o n  front  (1.3.5)  with  the  coefficient  tendency  second  term  (2.3.7)  of  flux  to  is  diffusive  increase  with  damping  with  i/>0. When the v a r i a b l e change  63  dQ (H) c(H) = — dH is  introduced to  balance result  A(x,t)  the set  (2.3.8)  continuity to  zero  equation  (1.3.5)  and w i d t h W(x) set  is a nonlinear diffusion  with  to u n i t y ,  (2.3.9)  bx 7  equation  the  equation  = v a2c(x,t) known as B u r g e r s '  mass  (Burgers,  1948). When  the  flux  law  ( 2 . 3 . 7 ) has the form Q(H)  = oH  2  (2.3.10)  + pH + r - vdH  dx 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 reduces  (2.3.9)  analytical  solution  transformation Burgers'  The  a  1951; Hopf,  linear diffusion equation,  is  well-known.  Applying  1950),  t o which the the  inverse  to the s o l u t i o n g i v e s the a n a l y t i c a l s o l u t i o n  equation,  (2.3.10). Since  to  (Cole,  and thus t o ( 1 . 3 . 5 ) w i t h  the  flux  to  relation  d e t a i l s of the s o l u t i o n are g i v e n i n Appendix 9.  i t i s rare to  find  nonlinear  partial  remarkable.  It provides  nontrivial  differential an e x c e l l e n t  analytical  equations,  solutions  this  opportunity  for  to  result an  is  exact  t e s t of the model w i t h a n o n l i n e a r flow l a w . Burgers'  equation  glaciological  literature  amplitude  (Johnson,  1968;  has  appeared  previously  i n papers on k i n e m a t i c waves of Lick,  When the i n i t i a l c o n d i t i o n  the  finite  1980).  is (2.3.11)  c ( x , 0 ) = A 6(x) and the boundary c o n d i t i o n s  1970; H u t t e r ,  in  are  64 c( c o , t ) the s o l u t i o n i s  = c ( - oo,t)  (see  (2.3.12)  Appendix .9)  V t  c(x,t) =  = 0  2  e +  J ?  X /(4l/t)  (eR  (e  J  - 1)  R  - 1 )  OO  (2.3.13)  _ 2 z  e  dz  xA/TTTF  The  parameter  single  i s e q u a l t o h/2v.  R  asymmetric  (2.3.13) i s a  d e c a y i n g hump which p r o p a g a t e s i n the  x d i r e c t i o n . A shock leading  This -solution  (vertical  front)  edge, but i s p r e v e n t e d  tends  to  positive  form  on  the  from d o i n g so by d i f f u s i o n . When  o=1/2 and 0 = 7 = 0 , we o b t a i n from (2.3.8) c(x,t) so t h a t  the s o l u t i o n  It  (2.3.14)  (2.3.13) i s a l s o the s o l u t i o n f o r  i s worth n o t i n g t h a t e q u a t i o n  a nonlinear t e s t , (Lighthill glaciers, that  = H(x,t)  and see  is also a test Whitham,  1955).  Nye (1960). I t  the d i f f u s i v e  ( 2 . 3 . 9 ) , as w e l l as  of • k i n e m a t i c For  wave  application  is readily  p r o p e r t y of  H(x,t).  apparent  behaviour t o waves on  from  i s c a r r i e d as a k i n e m a t i c wave at  (2.3.15)  velocity  dx = c ( x , t ) dt is  thus  the  is  i n c l u d e d i n the complete To  test  the  analytical solution  (2.3.16)  k i n e m a t i c wave v e l o c i t y .  (2.3.15) i s a t o t a l d e r i v a t i v e . )  (The l e f t  t h i s kinematic  solution  numerical  (2.3.9)  c(x,t)  dc(x,t) = i/d2c(x,t) dt dx"2"  which  being  wave  s i d e of behaviour  (2.3.13).  solution  of  (2.3.13) and ( 2 . 3 . 1 4 ) , I  (1.4.1) a g a i n s t set  the  the  material  65  1 .0  FIGURE 2 . 7 . N o n l i n e a r Test W i t h B u r g e r 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 solution (2.3.13) have been superimposed on the n u m e r i c a l model s o l u t i o n f o r a s i n g l e hump. Because of the c l o s e agreement (one p a r t in 103), the c u r v e s are indistinguishable in this Figure. 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 i v e n i n Table 2 . 3 .  velocity V(x,t)  t o be,  using  V(x,t)  (2.3.10),  =  Q(x,t) H(x,t)  H*0  = 0 I  s o l v e the f i n i t e d i f f e r e n c e  choosing  L  approximates  sufficiently  H=0  (2.3.17)  e q u a t i o n s on the i n t e r v a l  large  that  H(-L,t)=0  ( 2 . 3 . 1 2 ) . The c o n d i t i o n at +co  [-L,L],  adequately  i s not needed s i n c e  66 A  V  L  to  1.0  0.1  7.5  2.0  TABLE 2 . 3 . (1.4.1)  Parameters  is f i r s t  Since f i n i t e impulse initial  initial  equation  test.  order. difference condition  Figure 2.7,  superimposed parameters (generally  model  for Burgers'  0.05  0. 125  0.0  0.0  1/2  At  Ax  y  e  equations (2.3.11),  cannot I  represent  start  be  analytical  to>0.  solution  (2.3.13)  is  as dashed l i n e s on the n u m e r i c a l model r e s u l t s . are  given  to better  in  scheme  mass  with  The  Table 2 . 3 . The agreement i s so c l o s e  than t h r e e f i g u r e s )  distinguished.  conserves  iterative  the  the  i n s t e a d w i t h an  c o n d i t i o n ( 2 . 3 . 1 3 ) c ( x , t 0 ) at some s m a l l  In  cannot  o  This  test  fully  t h a t the  dashed  lines  shows t h a t the n u m e r i c a l  nonlinear  equations.  f o r 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 i n c e  The the  n u m e r i c a l s o l u t i o n f i n d s the c o r r e c t time response f o r the hump, the b e h a v i o u r of k i n e m a t i c waves i n the also  correct.  numerical  solution  is  67  2.4 THE ICE TRAJECTORY COMPUTER MODEL  2.4.1 I  INTRODUCTION  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 through detail  a  of s p e c i f i e d p a r t i c l e s of  time-varying  glacier.  ice,  constants  for  Glen's  flow  I have d e s c r i b e d t h e model i n  trajectory  the  bedrock  of the i c e p a r t i c l e s t o be  channel  geometry.  At  each  { h - | j = 1 , J } and the b a s a l J  by  the  tracked.  The  model uses t h e same mesh increments Ax and At as t h e  c o n t i n u i t y model ( S e c t i o n 2 . 2 ) , and the same  profile  topography,  flow 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 ,  used  they  i n Appendix 2 .  The i n p u t s t o the model a r e  the  as  the" c o n t i n u i t y  model  time sliding  assumptions  about  s t e p , the i c e t h i c k n e s s profile  { u . |j=1,J} J a r e a l s o used as i n p u t t o t h i s 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 displacement  vector  as  a  given  by  P(t), its  f u n c t i o n of t i m e . F o r a p a r t i c l e a t  p o s i t i o n P 0 a t time t 0 , t P(t)  P (2.4.1 )  o  where v ( x , t ) (2.4.1),  is  the  ice  velocity  field  (u,w,v).  To  solve  I s t a r t a t each time s t e p by f i n d i n g the v e l o c i t y  field  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) on  throughout  the  x-z  Appendix 2 ,  mesh  the v e r t i c a l p l a n e shown  in  Section A 2 . 1 .  component u ( x , z )  i n the g l a c i e r  Figure 2.8, First,  the  described  downslope  i s found at each meshpoint u s i n g the  form ( 1 . 4 . 3 4 ) of G l e n ' s flow l a w , w i t h given  and  centreline,  h(x),  o(x),  velocity integrated and  ug(x)  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  r a t e du/dx i s e s t i m a t e d ou(i,j) dx  by a f i n i t e =  u(i+1,j) DX i-1,j  difference - u(i-1,j) + DX  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 rate  in  dw/^y  is  given  strain  i n terms of u ( x , z ) and the c h a n n e l w i d t h  W(x) by ( 1 . 3 . 3 ) . The  incompressibility  condition  with  mass  conservation  69 gives dv = - du - dw 51 d"x 5y  Using 1972)  the a p p r o x i m a t i o n  (usually  (2.4.3)  v e r y good, e.g_.  t h a t b a s a l m e l t i n g i s n e g l i g i b l e as  concerned,  Rothlisberger,  f a r as mass b a l a n c e  i.e. (2.4.4)  v(x,0,t) = 0  (2.4.3)  is  Carnahan give  integrated  and  numerically  others,  1969,  by  p . 73)  Simpson's  (e.c[.  v(x,z,t).  flowline  average  is exactly  zero  on  down the c e n t r e of the c h a n n e l , and i s v e r y s m a l l  i n a narrow flow volume ( F i g u r e The  rule  from the bed t o l e v e l z to  The l a t e r a l v e l o c i t y component w ( x , y , z ) the  is  value  of  w  1.1)  c e n t r e d on  this  flowline.  i s z e r o 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 c o m p l e t e l y the  meshpoints  in Figure 2.8,  the d i s p l a c e m e n t  l e a v i n g each meshpoint P 0 and g o i n g t o i n t e r v a l At i s found by e s t i m a t i n g P(t+At) i.e.,  - P0(t)  a  the i n t e g r a l  = 1 v(P0,t) 2L  using  an  ice  P  scheme  in  a  time  ( 2 . 4 . 1 ) by  (2.4.5)  the b e g i n n i n g and end of is  estimated  the four s u r r o u n d i n g meshpoints  interpolation  Section A 2 . 3 . 1 .  f i e l d of the  At  the average of the v e l o c i t i e s at  v a l u e s of v a t  at  point  the time i n t e r v a l A t . The v e l o c i t y v ( P , t + A t ) the  determined  described  in  from  at time t+At Appendix 2,  70 2 . 4 . 3 THE ICE PARTICLE TRAJECTORIES At  each  time  s t e p , and f o r each of the N p a r t i c l e s  t r a c k e d , the c o o r d i n a t e s are  (6x,6z)  r e l a t i v e to a  r e c o r d e d . The d i s p l a c e m e n t s  are i n t e r p o l a t e d from the surrounding meshpoints, are then  of  program  checks  whether  the  at  the  of the  four  particles  new p o s i t i o n s are  w i t h i n the g l a c i e r mass, and saves the p o s i t i o n s at which i c e p a r t i c l e s using  At<0,  the  interpolated  reach the i c e  program  is  p a r t i c l e s backward i n time and upslope to  i n the time At  ice  and the new c o o r d i n a t e s  (i,j)  saved.  The  By  meshpoint  of the p a r t i c l e s  displacements  being  find  where  and  when  they  times  and  surface.  e a s i l y adapted (e.g.. from  entered  still  a  the  to track  borehole),  ice  mass  as  precipitation.  2.4.4  ACCURACY OF THE TRAJECTORY MODEL  The v e l o c i t y f i e l d that  shearing  i n t h i s model i s  parallel  to  the  obtained  glacier  bed  by  assuming  i s the dominant  component of d e f o r m a t i o n . The downslope v e l o c i t y u ( x , z ) found of  using  the s t r e s s e q u a t i o n s  and the m e c h a n i c a l  i c e . T h i s i s the o n l y p l a c e where  used.  The  other  velocity  the  components  force are  is  properties  equations  determined  k i n e m a t i c a l l y , u s i n g c o n t i n u i t y and geometric assumptions the l a t e r a l v a r i a t i o n of the flow When  the  assumption  of  are  purely about  field.  predominantly  h o l d s , t h i s approach works v e r y w e l l . The longitudinal strain rate,  then  shear d e f o r m a t i o n  fractional  error  and i n the v e r t i c a l v e l o c i t y , i s  in  quite  71 small.  The  unbalanced  leading  is  approximately  l o n g i t u d i n a l f o r c e s on  b a s a l shear f o r c e . (A7.5.9),  term  The d e t a i l s  a  e(x)  = 0  vertical  r a t i o of  column,  the  to  are g i v e n i n Appendix 7,  which i s r e p e a t e d here as  2h  the  the  equation  (2.4.6).  dtf' dtf' xx + h yy 6x 6x />gho  d v /ou o~x/ dz  6T  +  (n-1) a  xz  max (2.4.6)  When not  the  accurate,  relatively  assumption the  of p r e d o m i n a n t l y shear d e f o r m a t i o n  downslope  inaccurate.  The  velocity final  component  term  ( 2 . 4 . 6 ) may be of order u n i t y , or l a r g e r ,  u(x,z)  i n the e r r o r when d u / o z  or l a r g e r ,  when s t r e s s d e v i a t o r s  c o n t r i b u t e t o the e f f e c t i v e any  additional  the  ice,  so  approximations This  there  is  no  in this  term. may  as  Fortunately, that safe,  the  ice  of  order  other than the shear  o t h e r than * x Z  possibility  arise  at  predominant motions may be v e r t i c a l extension  small.  flows  an  ice  sinking,  downslope  in  always  XZ  of  softens  compensating  d i v i d e , where the and  longitudinal  both  directions.  the i c e v e l o c i t y i s v e r y s m a l l at an i c e d i v i d e ,  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 the t r a j e c t o r i e s  <J  s t r e s s . F u r t h e r m o r e , the presence of  s t r e s s components  situation  is  estimate  is  The second term i n the same e q u a t i o n may a l s o be unity,  is  is  not  large.  near an i c e d i v i d e s h o u l d be  To  so be  interpreted  o n l y 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 s e s are a l s o important  in  icefalls. The"  use  of  Glen's  flow  law,  and  the. inclusion  of  72 l o n g i t u d i n a l s t r a i n r a t e s ( t o the a p p r o x i m a t i o n of an  (2.4.6))  improvement over the t r a j e c t o r y models of Nagata  ( 1 9 7 7 ) , who  assumed t h a t du/dz was z e r o , and Dansgaard and Johnsen who set and  (I969[a])  du/dz to a c o n s t a n t at each x i n the lowest 400  zero  above t h a t l e v e l ,  for a given z,  in  a  metres,  and assumed t h a t du/dx was  model  for  the  Camp  is  Century,  constant Greenland  borehole. The  fact  that  my ' model  is  time  dependent  i s a l s o an  improvement over these models. One advantage of t h i s model i s t h a t the law  is s t i l l  not been  conservation  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  the  instance,  mass  case  with  some  other  trajectory  models.  the Weertman (1968) a n a l y t i c a l model f o r i c e v e l o c i t y  and temperature at Camp C e n t u r y , G r e e n l a n d , assumed vertical  strain  rate,  but  used  a  a  constant  h o r i z o n t a l v e l o c i t y given  i n d e p e n d e n t l y by i n t e g r a t i n g G l e n ' s flow law ( 1 . 4 . 2 2 ) t o result  For  like  (1.4.34).  Dansgaard  and Johnsen  t h a t 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  a  (1969[b]) showed of  a  d i s c r e p a n c y of 2°C between p r e d i c t e d and measured temperature  at  the  bed  at  Camp  Century,  Dansgaard and Johnsen but assumed the  form  the  get  Greenland.  likely  The flow model used by  (1969[a]) s a t i s f i e d of  the  cause  c o n t i n u i t y everywhere,  horizontal  velocity  component,  i n s t e a d of u s i n g G l e n ' s flow l a w . While  the  model  I  described  i n c l u d e temperature v a r i a t i o n s ,  that  i n S e c t i o n 2.2 does not yet is  a  simple  development  t h a t 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 c o n s t a n t v e r t i c a l  s t r a i n r a t e d v / d z independent of  depth,  at  each  position  x,  73 regardless  of  the  values  a l s o d e s c r i b e d another strain  rate  was  of the o t h e r terms i n ( 2 . 4 . 3 ) . They  model  weighted  (p. 43)  in  which  the  vertical  by the downslope v e l o c i t y u ( x , z )  each p o i n t . While t h i s makes the  vertical  strain  rate  curves  more  c l o s e l y resemble the expected shape i n a r e a l i c e mass,  still  does not s a t i s f y  incompressible  continuity locally  continuum  flowing  with  through such  at  it  (2.4.3).  a velocity  An field  would have t o l o c a l l y c r e a t e or d e s t r o y mass. Constant v e r t i c a l s t r a i n r a t e (or a s t r a i n r a t e s p e c i f i e d a priori)  may  analytical  be  a  reasonable  solutions  to  some  p h y s i c s i n v o l v e d (e.cj. R o b i n , on temperature p r o f i l e s ) , kinds layer  of f l o w ,  given  to  deriving  flow problems t o i l l u s t r a t e the  1955; f o r the e f f e c t  of  advection  and i t a r i s e s n a t u r a l l y i n some simple  1950, the  p . 233;  Nye,  sophistication  accepting otherwise q u i t e general appears  when  i_.e. h o r i z o n t a l s h e a r i n g r e s t r i c t e d to the  (e.c=. H i l l ,  However,  approximation  1951;  or  basal  Nye,  1957).  of c u r r e n t computer models  inputs,  this  assumption  now  be a n e e d l e s s l i m i t a t i o n . W h i l e the e r r o r s i n v o l v e d  may t u r n o u t ,  i n some s i t u a t i o n s ,  t o be s m a l l , i t i s  preferable  t o 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  the c o n t i n u i t y e q u a t i o n described  in  of which I am aware which a t t e m p t s t o in  a  Section 2.4.2  manner to  similar  find  the  to  vertical  component, and a l s o makes no a p r i o r i assumptions f i e l d o t h e r than preliminary  those  described  three-dimensional  in  model  Section  way  for  about the  flow  is  the  2.4.2,  the  f o r use i n the temperature dependent  I  velocity  of Jenssen (1977) f o r  G r e e n l a n d i c e c a p . T h i s model a l s o s o l v e s distribution,  the  use  the  temperature flow l a w . The  74 Jenssen model i s s t i l l  under development t o improve the boundary  t r e a t m e n t and and t o reduce i n a c c u r a c y r e s u l t i n g from the c o a r s e grid  imposed by computer l i m i t a t i o n s (Greenland was m o d e l l e d by  a 12x12 map g r i d , i s the  logical  numerical  10 p o i n t s d e e p ) ; next  solution  step of  that  the  however, the Jenssen approach must  complete  be set  taken  towards  of e q u a t i o n s  the  for  ice  sheets.  2.5 TESTING THE TRAJECTORY MODEL  2.5.1  INTRODUCTION  In S e c t i o n 2 . 5 ,  I describe  two  tests  of  the  trajectory  model. In  Section  and v e l o c i t y analytical  2.5.2,  field  I compare the steady s t a t e  calculated  by  my  computer  s o l u t i o n f o r a steady s t a t e i c e sheet  In S e c t i o n 2 . 5 . 3 ,  trajectories model  to  an  (Nagata,1977).  I use the b a l a n c e c o n d i t i o n (A5.8) at  i c e s u r f a c e t o show t h a t c o n t i n u i t y i s s a t i s f i e d  the  by the v e l o c i t y  field.  2.5.2 Nagata  NAGATA ICE SHEET TEST  (1977)  derived  an  surface p r o f i l e , v e l o c i t y f i e l d , two-dimensional state  case,  describe  the  analytical and  solution  streamlines  of  for a  the  steady  i c e sheet r e s t i n g on a f l a t b e d . For the steady  streamlines  and  Nagata  model  S e c t i o n A 1 5 . 3 . Nagata  (1978)  ice in used  trajectories detail this  ice  in  coincide.  I  Appendix 15,  sheet  model  to  75  E 0  CD O  c  Nagata (1977) Model Steady state ice sheet Basal sliding m=3 Velocity f i e l S and streamlines  o o CD  -2  •6723  1.0  0.5  0.5  1.0  x/L  FIGURE 2.9. N a g a t a S t e a d y I c e S h e e t . The analytical solution for the ice thickness, mass balance, particle paths, and velocity field u s i n g the p a r a m e t e r s i n e q u a t i o n ( 2 . 5 . 2 ) . The velocities have been multiplied by 250 y e a r s . The numbers a r e t h e t i m e i n y e a r s f o r i c e t o flow the l e n g t h of each p a r t i c l e p a t h .  explain Field ice  the  concentration  in Antarctica. velocity  depends  on  Weertman-type  u(x) the  The  of  relation  at the M e t e o r i t e  N a g a t a model assumes t h a t  i s constant  basal  meteorites  shear  throughout stress  (see Appendix  a vertical  forward  column  and  through  a  S e c t i o n A8.3.1) of  the  (1.4.25) 8,  the  Ice  76  T  X  (km)  FIGURE 2 . 1 0 . Growth Of Nagata Ice S h e e t . The i c e s u r f a c e elevation is shown at i n t e r v a l s of 500 years, starting from i c e - f r e e conditions. The flow parameters, are g i v e n i n T a b l e 2.4 and the mass balance i n F i g u r e 2 . 9 . The steady s t a t e 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) t o one p a r t i n 1 0 3 .  form m u(x)  = A T  (2.5.1)  A l t h o u g h i t i s not c l e a r i n the o r i g i n a l p a p e r , a l s o assumes t h a t the equal  vertical  velocity  component  t o a c o n s t a n t b a l o n g the upper s u r f a c e  The s c a l e of the i c e sheet model i s s e t  the Nagata model v(x,z)  is  of the i c e s h e e t .  by H , the  thickness  at  the i c e d i v i d e . While  the  form of the mass b a l a n c e 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 model  is  well  suited  sheets,  the  Nagata  t o t e s t i n g the a c c u r a c y of my n u m e r i c a l  77  Distance (km) FIGURE 2 . 1 1 . V e l o c i t y F i e l d For Nagata M o d e l . The v e l o c i t y v e c t o r s calculated by the Waddington trajectory model for the steady s t a t e p r o f i l e in F i g u r e 2.10 and u s i n g flow parameters i n Table 2 . 4 . The velocities have been m u l t i p l i e d by 250 y e a r s (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 ) .  trajectory  program. I compare my n u m e r i c a l r e s u l t s  t o the Nagata  model shown i n F i g u r e 2.9 f o r the c o n s t a n t s i n Table F i g u r e 2.10 shows the growth of the  Nagata  steady s t a t e u s i n g the p r o f i l e model d e s c r i b e d The  steady  state  profile  velocity  in  ice  sheet  Section  to 2.2).  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 t o w i t h i n one p a r t The s o l i d v e c t o r s  in  2.4.  in  103.  F i g u r e 2.11  show  the  steady  state  f i e l d c a l c u l a t e d by my n u m e r i c a l m o d e l . 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  factor  of  250 y e a r s .  The  velocity  78 A bar"  m  2  b m a"  a'  1  1.0  100.  2  1  H m  L km  kg m " 3  g m s"2  3000  454.6  910.  9.8  p  At  AX  a  km  10.0  7.215  TABLE 2 . 4 .  Parameters f o r Nagata i c e  sheet.  f i e l d a g r e e s w i t h the a n a l y t i c a l s o l u t i o n ( F i g u r e 2.9)  to w i t h i n  a few p a r t s  agreement  is  only  i n 10 3  to  two  (except near the t e r m i n u s where the parts  i n 1 0 2 , because the mass b a l a n c e i n the  n u m e r i c a 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 b a l a n c e , goes t o -oo ) . F i g u r e 2.12 shows the t r a j e c t o r i e s sheet at time t = 0 . ; streamlines  in  intervals.  The  f o r i c e e n t e r i n g the  these are the same f i v e p o i n t s shown f o r the  Figure 2.9. total  The  residence  arrowheads times  indicate  The good agreement between the n u m e r i c a l solution  indicates  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 f i e l d t o get  that  250 year  a l o n g each of the  s t r e a m l i n e s are compared t o the a n a l y t i c a l v a l u e s  analytical  ice  the  i n Table  results  velocity  integration  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 .  of  the  five 2.5.  and field  the is  velocity  79  0  200  100  300 X  400  500  (km)  FIGURE 2 . 1 2 . T r a j e c t o r i e s In Nagata M o d e l . The s o l i d c u r v e s a r e the p a r t i c l e p a t h 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 by i n t e g r a t i n g the velocity field in Figure 2.11. 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 r e s i d e n c e times i n the ice sheet are compared to values for the analytical solution.  2.5.3  SURFACE MASS CONSERVATION TEST  Conservation vector n interface  implies is  of  mass  that  equal  at  the  the g l a c i e r  normal  s u r f a c e w i t h normal  velocity  of  the  t o the sum of the normal i c e v e l o c i t y  the s u r f a c e a c c u m u l a t i o n r a t e ,  plus  i..e.  — •n = v •n + a • n  dt  This equation  ice-air  i s d e r i v e d i n Appendix 5, S e c t i o n A 5 . 2 .  (2.5.2) When  ice  80 A n a l y t i c a l model years  N u m e r i c a l model years  Streamline Number  6723. 4606. 3322. 2346. 1466.  6746. 4625. 3335. 2350. 1477.  1 2 3 4 -> 5  TABLE. 2 . 5 .  Residence  thickness  change  measured  times i n Nagata i c e s h e e t .  rate  dh/dt  and  mass  normal t o the bed, and the  velocity  are p a r a l l e l t o and normal t o the bed, oh = v ( x , h )  A(x)  derivation  is  v e c t o r a are  components  (u,v)  ( 2 . 5 . 2 ) becomes  - u ( x , h ) djh + A(x) 6x  St where  balance  now  a  scalar  (Section 2.4.2  (2.5.3)  giving  the  mass  balance.  and Appendix 2, S e c t i o n A 2 . 2 )  My  of  the  v e l o c i t y f i e l d uses o n l y i n c o m p r e s s i b i l i t y ( 2 . 4 . 3 ) and the b a s a l boundary  condition  (2.4.4);  (2.5.3)  independent check of the degree conserves velocity  to  mass. Models which use ( e . £ . Weertman, 1968;  not have t h i s c o n s i s t e n c y F i g u r e 2.13  (a)  can-  which  the  be  used  as  velocity  an field  ( 2 . 5 . 3 ) t o d e r i v e the v e r t i c a l  Budd and o t h e r s ,  1971, p . 42)  do  check.  shows the mass balance A(x) and the  surface  r i s e d h / d t f o r a t y p i c a l time s t e p , w i t h parameters d e s c r i b e d by Table 2.6,  during  the  growth of the Nagata i c e sheet m o d e l . I  ' a p p l y the t e s t a t a time (2500 y e a r s ) growing  vigorously  because  sheet  is  c o n d i t i o n s a t t h a t time impose  the  most s t r i n g e n t c o n d i t i o n s on a c c u r a c y . r e s i d u a l of right  (2.5.3),  when  the  F i g u r e 2.13  i . . e . the d i f f e r e n c e  between  ice  (b) shows the  s i d e s of the e q u a t i o n on s u b s t i t u t i n g the v a l u e s  left  the and  of h , A ,  81  FIGURE 2 . 1 3 . S u r f a c e Mass C o n s e r v a t i o n T e s t . Time=2500 y e a r s d u r i n g the growth of the Nagata i c e sheet model t o a steady s t a t e , (a) shows the mass b a l a n c e (solid curve) and the r a t e of s u r f a c e r i s e d h / d t (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 v e r y s m a l l (note s c a l e change of 1 0 ~ 2 , i n d i c a t i n g t h a t 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 i n c r e a s e s near the terminus because t h e r e are few meshpoints in each vertical column; vertical integration is inaccurate.  u , and v from the n u m e r i c a l  error  is  t h r e e o r d e r s of magnitude  l e s s than the average magnitude of  the  mass  near  are  balance,  except  i n s u f f i c i e n t mesh p o i n t s accurate  vertical  solution.  i n any  integration  the  The  terminus  vertical of  residual  where  column  to  there  guarantee  ( 2 . 4 * 3 ) . T h i s r e g i o n has no  82 Ax km  Time a  7.215  2500.  TABLE 2 . 6 . Nagata  DZ m  At a  150.  10.  i c e sheet s u r f a c e  effect  on the t r a j e c t o r i e s  tests  (not  shown)  on  boundary t e s t .  investigated in  this  mass  Similar  the S t e e l e G l a c i e r Model 1 ( F i g u r e  r o u t i n e l y g i v e r e s i d u a l s of the o r d e r of one p a r t average  work.  balance,  even  though  the  3.3)  i n 10 2 of  mass  the  balance,  is  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 w i d t h i s v a r i a b l e , and the flow law  includes  Results satisfies  of the  a  height-dependent  these  longitudinal  tests  indicate  that  the  continuity  equation  (2.3.3)  to  strain  rate.  velocity  field  a  very  good  a p p r o x i m a t i o n . T h i s t e s t v e r i f i e s the a c c u r a c y 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  stable  heavy i s o t o p e s O 1 8 of oxygen and D (deuterium)  of hydrogen i n g l a c i e r of c l i m a t i c change procedure  ( e . c i . Dansgaard and o t h e r s ,  requires  assumptions  flow.  In t h i s c h a p t e r ,  that  the  past  distribution  i c e have been w i d e l y used  about  I investigate  climate  is  known,  the  as  1969;  indicators 1971).  p a t t e r n of  a r e l a t e d problem; can  the  This  glacier assuming  stable  isotope  be used t o 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  glaciers?  The example I c o n s i d e r  Glacier,  Yukon T e r r i t o r y . T h i s g l a c i e r was observed t o surge i n  1966-1967 ( S t a n l e y , in  Chapter 2  throughout  to  the  trajectories.  1969). find  surge Knowing  i n t h i s chapter  the cycle, the  and  to  calculate  trajectories  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 the  isotopic  Steele  described  i c e s u r f a c e and the v e l o c i t y  field  the  ice  and  the  isotopic  place  the  material  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 ,  showing  the  I used the computer models  c o m p o s i t i o n (from c l i m a t e ) at the time and was  is  and  has a l l o w e d me  surface  profiles  d i s t r i b u t i o n at a s e r i e s of times d u r i n g  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 optimum  b o r e h o l e and s u r f a c e sampling l o c a t i o n s f o r an i s o t o p i c  study of p a s t flow In S e c t i o n 3 . 2 , history.  of  In S e c t i o n s  patterns. I describe  the S t e e l e G l a c i e r and i t s  3.3 and 3.4 I d e s c r i b e  the computer  surge models  I used t o s i m u l a t e the S t e e l e G l a c i e r . Model 1 i n S e c t i o n 3.3  is  84 a  detailed  model  based  on a l l the a v a i l a b l e d a t a .  l i m i t e d s l i d i n g d a t a and an a p p r o x i m a t i o n i n the it  is  not  at  present  it  may  become  s i m p l i f i e d version better sliding  observations.  isotopic calculations describe two  so.  matched  Model 2  is  in Sections  some  possible  isotopic  indicate  period  that  appropriate  a  surge  trajectory  Model 2  i n S e c t i o n 3.4  to  resolution  the  3.6 and 3 . 7 .  for  calculations of  and  in  approximately  f o r the S t e e l e G l a c i e r i f the p r e s e n t  a  surge  periodicity  of  100 y e a r s  the  pattern.  In  interface  which  for  DEL  hidden by n o i s e  a  Isotopic were  at  i n the a c c u m u l a t i o n r e g i o n when a surge  Even f o r the i s o t o p i c - p r e c i p i t a t i o n model most  one  is  balance  favourable  t o the f o r m a t i o n 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 exceed  3.6  of  97 y e a r s .  d i s c o n t i n u i t i e s occur i n the i c e a l o n g s u r f a c e s  began.  Steele  Section  mass  flow  I  I p r e s e n t 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  with  ice-air  and  present  s u r f a c e 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 )  the  the  In S e c t i o n 3 . 5 ,  p r e c i p i t a t i o n at  average l o n g - t e r m c l i m a t e and of the g l a c i e r  model  is a  of  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  Section 3.7,  model  used f o r the t r a j e c t o r y  relations profile  balance,  computer  the use of s t a b l e i s o t o p e s i n g l a c i o l o g y ,  G l a c i e r . The i c e s u r f a c e  and  mass  the best model f o r p a r t i c l e  c a l c u l a t i o n s . W i t h more complete d a t a and refinements,  Because of  not  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 i n the S t e e l e G l a c i e r e n v i r o n m e n t .  85 3.2 STEELE GLACIER 3.2.1  GENERAL DESCRIPTION  Steele Glacier glacier  on  F i g u r e 3.1) Prior  the  (61°10'N,  northeast  of the S t .  140°15'W) slopes  is  a  surging  of the I c e f i e l d Ranges  E l i a s M o u n t a i n s , Yukon T e r r i t o r y ,  t o 1963 i t was c a l l e d the Wolf G l a c i e r .  Bostock  (1948, p . 99) s t a t e t h a t t h i s i s more  often-used varies  notation  valley  Wolf  Creek  Canada.  (Sharp (1951) and proper  than  G l a c i e r . ) . The g l a c i e r  from 34 to 44 km, and the w i d t h of the  (see  main  the  length  channel  is  one t o two km. The main c h a n n e l flows down from 3000 m e l e v a t i o n on the n o r t h s i d e of M t . S t e e l e ( F i g u r e in  Steele  between  3.2)  t o 1200 m e l e v a t i o n  C r e e k , where the i c e i s s t a g n a n t and  moraine-covered  surges.  The c o n t i n e n t a l s l o p e of the I c e f i e l d Ranges i s Annual  precipitation  G u l f of A l a s k a (see (Wood,  1972).  drops  from  Figure 3.1),  The f i r n  300 cm a " 1  from  system  the n o r t h face ( l a b e l l e d  tributaries. F i g u r e 3.2) Glacier  of  also contribute a substantial balance.  In i t s  presently  accumulation  is  (0)  i n F i g u r e 3.2)  of  reaches,  accumulation  Lake  of  the  basins  These t r i b u t a r y i c e streams (see  mass  Kluane  l i n e on the S t e e l e G l a c i e r i s  M t . S t e e l e (5070 m). In i t s upper extensive  at Y a k u t a t on the  t o about 35 cm at  v e r y h i g h (2400 m or m o r e ) . A major source avalanches  semi-arid.  Steele and  (1)  is  an  converging  through  (5),  f r a c t i o n of the S t e e l e  lower r e a c h e s , the S t e e l e G l a c i e r  makes a 9 0 ° bend t o the e a s t , and e n t e r s a s t r a i g h t v a l l e y formed by the Wolf Creek monocline ( S h a r p ,  and 1943).  narrow  86  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 a r e a s are major g l a c i e r s and i c e f i e l d s . The triangles 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 M o u n t a i n s ) . The S t e e l e G l a c i e r i s on the n o r t h - e a s t ( c o n t i n e n t a l ) s l o p e (top c e n t r e ) . 3.2.2 A  glacier  GLACIER SURGES surge  (e.cj.  Meier  and P o s t ,  1969)  i s a short  p e r i o d of v e r y 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 an  ice  reservoir  area  from  t o an i c e r e c e i v i n g a r e a downstream. A  surge i s f o l l o w e d by a l o n g e r p e r i o d of s t a g n a t i o n and  ablation  87 in  the  receiving  reservoir;  area,  these areas do  accumulation  and  a  surge  ice displacement periodic  renewed  not  i c e b u i l d u p i n the  necessarily  correspond  like  the  Steele,  may be 500 m a " 1 of one  phenomenon.  to For  the  maximum  t o 10 km a " 1 ,  10 km. glaciers  Surging  For a  velocity  appears  to  be  quiescent  phase  from 20 t o 150 y e a r s . Examples  the T r a p r i d g e G l a c i e r ( i n the S t e e l e Creek watershed)  i s o n l y t h r e e km l o n g , and the Muldrow G l a c i e r on  which  Mt. McKinley,  A l a s k a which i s over 50 km l o n g . There has been s p e c u l a t i o n the  A n t a r c t i c i c e sheet may surge ( H o l l i n ,  Surging g l a c i e r s Alaska  are found i n many  (Tarr and M a r t i n ,  Yukon  Territory.  (Hattersley-Smith, 1969),  the  1914, p .  (Post,  168),  of  1969),  (Hewitt,  the  in  that 1969).  world,  e«2.  B r i t i s h Columbia and the  1969),  1964; L 0 k e n ,  Karakoram  parts  1969; W i l s o n ,  the  Arctic  Iceland  1969),  Islands  (Thorarinsson,  the P a m i r s , T i e n Shan,  Caucasus, and Kamchatka ( D o l g o u s h i n and O s i p o v a , occur  a  l i k e , the S t e e l e , the surge  G l a c i e r s of a l l s i z e s have been observed t o s u r g e . are  the  w i t h downstream  d u r a t i o n i s t y p i c a l l y one t o two y e a r s , w i t h a lasting  ice  to  a b l a t i o n zones d e f i n e d by mass b a l a n c e .  large v a l l e y glacier during  with  both temperate and c o l d or s u b p o l a r  does not appear t o be t r i g g e r e d by c l i m a t e  1975).  Surges  glaciers.  Surging  variations.  The h i g h v e l o c i t y d u r i n g s u r g i n g i s g e n e r a l l y a t t r i b u t e d rapid basal s l i d i n g . Various surging  have  thermal  regulation  I969[b]),  hypotheses  on  the  mechanism  to of  been put f o r w a r d . The more p l a u s i b l e ones i n c l u d e  stress  (Robin,  1955;  instabilities  and b a s a l water f i l m i n s t a b i l i t i e s  Clarke,  (e.g.  Post,  1976;  Lliboutry,  1960; R o b i n ,  (Weertman, 1962,  1969;  1969) Robin  88 and  Weertman,  1973;  Budd,  g l a c i e r s u r g i n g has yet i n v e s t i g a t e the e f f e c t distribution  to  1975). A concensus on the cause of be  reached.  In  this  chapter,  of p e r i o d i c s u r g i n g on the s t a b l e  isotope  i n the S t e e l e G l a c i e r . The amount of b a s a l s l i d i n g  i s i m p o r t a n t t o t h i s q u e s t i o n ; the mechanism which causes i t n o t . For t h i s r e a s o n , consistent 1966-1967,  I  I specify a basal  w i t h the o b s e r v a t i o n s and I use  this  as  sliding velocity  is  us(x,t)  of the S t e e l e G l a c i e r surge of  a  boundary  condition  for  the  computer m o d e l .  3 . 2 . 3 OBSERVATIONS OF THE STEELE GLACIER The  flow p a t t e r n b e f o r e and d u r i n g a surge i s b e t t e r known  f o r the S t e e l e G l a c i e r than f o r Scientific Society  expeditions  and  watershed  led  in  by  most  sponsored  1936,  surging  glaciers.  by the American G e o g r a p h i c a l  W. A . Wood  1935,  other  explored  1939,  and  the  Steele  Creek  1941. Work i n c l u d e d the  e x p e r i m e n t a l use of o b l i q u e a e r i a l photography f o r mapping areas of h i g h r e l i e f . These a i r photos and the panoramas c o n t r o l p o i n t s g i v e i n f o r m a t i o n on the i c e s u r f a c e s t a t e of flow d u r i n g the q u i e s c e n t Sharp  (1943)  described  from  ground  e l e v a t i o n and  phase.  the  geology  of  observed the g l a c i e r s of the S t e e l e Creek b a s i n 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,  Steele  Creek,  (Sharp,  1947),  1951).  The Surveys and Mapping Branch of the Department of E n e r g y , Mines  and  Resources,  Government  of Canada, o b t a i n e d 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 summer  of  1951.  Flight  lines  A13232/33  covered  in  the  the S t e e l e  89  61°20  61°05 140° 30  U0°  FIGURE 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 . 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 c h a n n e l w i d t h f u n c t i o n ( F i g u r e 3.3 (b)) f o r the main ice s t r e a m . 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 included i n the mass b a l a n c e f u n c t i o n ( F i g u r e 3.3 (a)). G l a c i e r . T h i s photography 115F at of  the  Survey, was  the s c a l e 1:250,000; i t was a l s o used to p r e p a r e Steele  Glacier  at  1967). O b l i q u e a e r i a l  obtained  periodically,  (American G e o g r a p h i c a l Survey).  was used t o p r e p a r e the government map  In  1960,  the  s c a l e 1:25,000  photography between  predicted  map  (Topographical  of the S t e e l e  Glacier  1960 and 1965 by W. A . Wood  S o c i e t y ) and by A . Post ( U . S .  Post  a  an  Geological  imminent surge f o r  the  90 Steele Glacier of  (Wood, 1972), and i n 1965, noted the f i r s t  i n c r e a s e d flow on the upper S t e e l e D u r i n g the 1966-1967  surge,  o b t a i n e d by Wood, pre-1941 the  Surveys  and  (1969)  1972).  h i g h o b l i q u e a e r i a l photos were  survey s t a t i o n s were r e - o c c u p i e d , and  Mapping  Branch  obtained  coverage on August 13 and September Stanley  (Wood,  signs  measured  v e r t i c a l a i r photo  15, 1966.  Wood  (1972)  and  d i s p l a c e m e n t s of i d e n t i f i a b l e s u r f a c e  m a r k i n g s . F e a t u r e s o r i g i n a l l y between the Hodgson c o n f l u e n c e and the 90° bend ( F i g u r e 3.2) velocity  of  5 km a " 1  all  advanced  at  of  same  1967.  There  was  no  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 m a r g i n s . Bayrock (1967) terminus  the  d u r i n g 1966, w i t h a t o t a l d i s p l a c e m e n t of  r o u g h l y 8 km when the surge ended l a t e i n significant  roughly  observed  the  details  advance, and the r e a c t i v a t i o n of stagnant  active  ice  (3650 m a " 1 ) .  30 m h i g h  moved  forward  at  of  the  i c e . A bulge  about  10 m d " 1  T h i s bulge sometimes has been c a l l e d the t e r m i n u s  i n the l i t e r a t u r e  on  the  1966-1967  surge.  Alford  from  the  W h i t e h o r s e o f f i c e of the Water Survey of Canada o b t a i n e d monthly air  photographs  1966-1967. September 1967,  the  (Thomson,  of  the  advancing  He observed an advance of 6000 f e e t 10,  1966 and January  active  bulge  had  15,  (1830 m)  1967 ( R o o t s ,  slowed  to  between  1967). By August  2 m d"1  (730  ma-1)  1972).  S t a n l e y (1969) i d e n t i f i e d elevation  bulge d u r i n g the w i n t e r of  changes  three  d u r i n g the s u r g e .  based  on  surface  An upper z o n e , e n t i r e l y i n  the a c c u m u l a t i o n  area,  in  the  1966-1967  In a m i d d l e zone (the i c e r e s e r v o i r zone)  from  surge.  apparently  zones  a p o i n t above the f i r n l i n e  was  not  (based on S t a n l e y ' s  involved  description,  I  91 locate  this  p o i n t at about x=7 km) t o a p o i n t about 3 km above  the 9 0 ° bend (at x=27 km i n my m o d e l ) , t h e r e was a net of  the  ice  surface,  with  a maximum v a l u e of  Hodgson-Steele c o n f l u e n c e . In r e c e i v i n g a r e a ) the s u r f a c e D u r i n g the w i n t e r of year-long  the  lower  130 m above the  remaining  zone  (ice  r o s e by up to 100 m.  1966-1967,  the Hodgson G l a c i e r began a  surge d u r i n g which i t pushed the s t i l l - s u r g i n g S t e e l e  i c e stream t o one t h i r d of i t s normal w i d t h at the i c e The  lowering  effect  of  the  Hodgson  surface.  surge on deep i c e i s unknown. The  Hodgson i c e formed a l a r g e l o b e e x t e n d i n g t h r e e km down the main S t e e l e c h a n n e l . T h i s t r i b u t a r y surge may have been t r i g g e r e d a  by  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 l o w e r i n g of  the S t e e l e 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 Steele  Glacier  earlier. most  below  observations,  the bend had been stagnant  Based on b i o l o g i c a l  recolonization  the  s i n c e 1916 or  rates  inside  the  r e c e n t t r i m l i n e , he e s t i m a t e d t h a t the l a s t advance ended  i n the  period  1966-1967 Wood  1840-1890,  j..e.  115  to  75 y e a r s  before  the  surge. (1972)  over two km  p o i n t e d out e v i d e n c e f o r i c e d i s p l a c e m e n t s of  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 t h a t r e g i o n was h e a v i l y crevassed  in  1941, and the c r e v a s s e s l o o k e d s e v e r a l y e a r s o l d .  The e v i d e n c e suggests t h a t the S t e e l e 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. S h a r p ' s  (1951)  e s t i m a t e of the time of the p r e v i o u s t e r m i n u s advance would then imply a b u i l d u p time of 50 t o 90 y e a r s . Sharp (1951) p o i n t e d out moraines 1941  100 m t o 150 m above  the  i c e s u r f a c e below the bend. An a b l a t i o n e s t i m a t e of 2 m a " 1  implies  the maximum i c e e l e v a t i o n o c c u r r e d 50 t o 75 y e a r s p r i o r  t o 1941 . There are no d i r e c t o b s e r v a t i o n s regular  looped  moraines,  on the S t e e l e  or of d i s r u p t e d o g i v e p a t t e r n s  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 e p i s o d e s . periodicity elsewhere others, computer  is  based  on  observation  ( e . £ . Variegated G l a c i e r , 1977)). model  The is,  assumption at  best,  of  Alaska, of  an  The  documented  (e.£.  Bryson  exact p e r i o d i c i t y for  and  Goodman,  and the  significant  surges has  been  1980; G r i b b o n ,  1979). I f the S t e e l e 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 n the range 50 t o 150 y e a r s .  of  glaciers  (Bindschadler  approximation;  of  which  assumption  surging  c l i m a t e v a r i a b i l i t y on the time s c a l e of g l a c i e r widely  Glacier  is  93 3.3 NUMERICAL MODEL 1 3.3.1  FLOW LAW CONSTANTS AND SHAPE FACTOR  The computer model i n i t s p r e s e n t is  isothermal.  Temperatures  in  Steele Glacier (Jarvis  and  1976)  1966-1967  -7°C.  following It  shear  the  melting  the  Clarke,  point.  takes  The  upper  1974;  100 metres of  Clarke  and  ice the  Jarvis,  surge were i n the range - 1 ° C t o  i s l i k e l y t h a t the b a s a l  deformation  form assumes t h a t the  i c e , w i t h i n which most of  place,  assumption  is  at  of  or  near the  temperate  ice  the  pressure is  not  unreasonable. I use the c o n s t a n t s A = 5.3 in  Glen's  (Paterson, 10-15  s"1  1981, p .  39)  kPa'3  n=3  flow law ( 1 . 4 . 2 2 ) f o r the S t e e l e G l a c i e r s i m u l a t i o n .  S i n c e most of the motion i n a surge i s due t o in  the  (3.3.1)  flow  law  constants  have  little  sliding, effect  on  changes 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 shape f a c t o r  (see  i s unknown. I  use  a  is raised  to  S e c t i o n A 7 . 4 ) of s = 0.8  independent the  nth  of  power  introduces velocity.  p o s i t i o n x . S i n c e the shape f a c t o r (see  equation ( 1 . 4 . 3 4 ) ) ,  substantial However,  deformation t o t a l motion.  of  the  as  uncertainty pointed  out  uncertainty  into above,  the the  in  s  deformation internal  S t e e l e G l a c i e r i s a s m a l l component of  the  94 3 . 3 . 2 BED TOPOGRAPHY There a r e no Glacier.  published  Jarvis  reached d e p t h s without  any  the  Steele  and C l a r k e (1974) and C l a r k e and J a r v i s  (1976)  exceeding indication  ice  depth  100 metres of  data  with  for  a  hot-point  bottom; t h e 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 r g e v a l l e y exceed  600 metres  (e.g..  Lowell  Glacier,  based on monopulse r a d i o - e c h o s o u n d i n g , personal  communication).  drill  A  rough  glaciers  S t . E l i a s Mountains,  1977;  G . K. C . C l a r k e ,  estimate  of  the  e.cj. above t h e c o n f l u e n c e of t h e S t e e l e and the Hodgson obtained  by  assuming  a  b a s a l shear s t r e s s  depth  can be  of one bar and  measuring t h e i c e s u r f a c e s l o p e t o be c=0.03 from t h e (Topographical  1951 map  S u r v e y , 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 i v e s t h e depth as ( a p p r o x i m a t e l y )  h =  xz  = 472 metres  s/>go Alternatively, and  using  can  (3.3.2)  assuming flow by s i m p l e shear  with  no  sliding,  t h e 1951 v e l o c i t y of 25 m a " 1 measured a t t h i s  by Wood (1972)  (1.4.38)  with  the  flow  law c o n s t a n t s  point (3.3.1)  gives the depth estimate (n+2)V h = Since  the  motion  of  2A( Syoga)  the  Steele  l/(n+D = 445 metres  n  (3.3.3) Glacier  i s o b v i o u s l y not an  example o f s t e a d y n o n s l i p f l o w i n a c y l i n d r i c a l  channel,  these  are merely rough e s t i m a t e s . The p r e s e n c e of b a s a l s t r e s s of l e s s than  one  bar,  or  of nonzero b a s a l s l i d i n g would r e s u l t i n an  95 overestimate  of  perturbations  the  and  true  thickness.  lateral  asymmetry  c h a n n e l , I assume the c e n t r a l broken  curve  topography  shown  in  i n F i g u r e 3.3  Neglecting  flowline  Figure 3.2. (c),  caused  I  (x  by  the  stress  bends i n the  axis)  follows  the  For the S t e e l e G l a c i e r bed  use  an  exponential  function  h a v i n g the form -bx h(x) = ae where  a,  + c  (3.3.4)  b , and c are c o n s t a n t s determined by f i t t i n g  the  three  points: ( 1 ) 2900 metres e l e v a t i o n at the bergschrund stream  of  the  main  (x=0).  (2)  1200 metres at the 1978 t e r m i n u s p o s i t i o n (x=42 km).  (3)  1650 metres  Glaciers  (X=18  400 metres The  ice  confluence  t h i s p o i n t i n 1951  longitudinal  F i g u r e 3.3 The  the  of  the Hodgson and S t e e l e  km); t h i s v a l u e would g i v e an i c e depth of  at  1951  at  (Topographical Survey,  surface  profile  is  about 1969).,  shown  in  Glacier  is  (c). actual  topography  beneath  the  Steele  undoubtedly more c o m p l i c a t e d . T h i s a p p r o x i m a t i o n means t h a t computer  model  the observed  cannot  ice surface  be expected elevations.  to q u a n t i t a t i v e l y  the  reproduce  96 3.3.3  CHANNEL WIDTH  The w i d t h of the main i c e stream of S t e e l e G l a c i e r estimated ridges  by  on  measuring  the  topographic  (Topographical and  115G  Survey,  115F(El/2),  F i g u r e 3.2  is  model  mass b a l a n c e thickness  not a v a i l a b l e Instead, ways  of the l a t e r a l  map  the  at  scale  at  the  is  scale  not  1:250,000,  a  reliable  1:25,000  tributaries  s h o u l d be  (from  the  included  the  and c o u p l e d t o the main  the  computer  model  in  its  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 the  width  a p p r o x i m a t i o n . In f a c t ,  the  Hodgson  present  above  this  point.  everywhere i n t h i s r e g i o n , t h e i r mass b a l a n c e differ, last  surged  their (1)  and the two g l a c i e r s  shallower  depths,  t h r o u g h ( 5 ) ( F i g u r e 3.2)  reduce  are  I use  t o i n c l u d e the minor  i n t h i s way would tend  the  a  simple  not  equal  functions  the S t e e l e ) .  the  probably  is  do not surge t o g e t h e r  i n 1967-68, one year a f t e r  much  are  This  t h e i r surface gradients  is  form.  f u n c t i o n s . Since  Glacier  by  i n approximate  t o those of the S t e e l e a t t h e i r c o n f l u e n c e , widths  of  mass balance  t h i s option  comparable  their  in  stream  and  of  depth  or  ice  discharge  sum  which  i n d i c a t i o n of  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 ; in  through  moraine  as s e p a r a t e i c e streams w i t h t h e i r own bed and  functions,  and  separation  be  1969). The p u b l i s h e d government mapsheet  drawn),  channel w i d t h . I d e a l l y , numerical  the  can  probably  (the Hodgson Because  of  tributaries to  grossly  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 addition r e g i o n s of  to  the  mass b a l a n c e  confluence.  is  i n c l u d e d i n Model 1  f u n c t i o n ( F i g u r e 3.3  (a))  by  an  at  the  97 3.3.4  MASS BALANCE  No comprehensive  mass b a l a n c e measurements have  been  made  19,  1974,  on the S t e e l e G l a c i e r . Between an average date of J u l y and  an  average  date  of  measured a b l a t i o n at e i g h t  J u l y 2,  1975, C o l l i n s  survey t a r g e t s near the  bend (x=30 km). The measurements ranged from 1.57 ice,  with  a mean of 2.15 m. A l l o t h e r evidence  the lower S t e e l e G l a c i e r , S t a n l e y to be 1.5 m a " 1  (unpublished) right  m t o 2.77 m of is  (1969) e s t i m a t e d  based on downwasting of  angle  i n d i r e c t . On the  ablation  i c e shown t o be stagnant  by Sharp ( 1 9 5 1 ) . Since  snowfall  extrapolating valleys  is  mass  can  balance  s t r o n g l y on l o c a l  information  at best a r i s k y p r o c e d u r e .  p a t t e r n s of mass b a l a n c e , the  depend  Icefield  Ranges  even  conditions,  from  adjacent  However, the  qualitative  and i t s o r d e r of magnitude  throughout  g i v e some c o n t r o l on r e a s o n a b l e e s t i m a t e s  f o r the S t e e l e G l a c i e r . Marcus and Ragle (1970) r e p o r t e d accumulation  measurements  on  a  traverse  a c r o s s the  Ranges from the lower Seward G l a c i e r t o the Kaskawulsh The  values  for  the  Kaskawulsh  are  east s i d e of the I c e f i e l d Ranges ( F i g u r e 3 . 1 ) . the  elevation,  decreasing  to  ice  the  on  the  The p r e c i p i t a t i o n  equivalent  0.35 metres  1615 metres e l e v a t i o n . These a u t h o r s  of  Glacier.  Kaskawulsh i n c r e a s e s w i t h e l e v a t i o n . The 1964-65 w i n t e r  a c c u m u l a t i o n was 1.7 metres of  summer  Icefield  i n s t r u c t i v e because both  Kaskawulsh and S t e e l e 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  winter  melting  on the i c e p l a t e a u  Kaskawulsh  Glacier.  of  also  at  ice  2640 metres  equivalent  reported  at  negligible  at D i v i d e S t a t i o n 30 km west  Divide  Station  at  2620 metres  98  1  ol  10001 0  1  1 10  1  r  i  L  1 20  30  40  50  x (km) FIGURE 3.3..Model 1 For S t e e l e G l a c i e r . (a) Mass balance. The dashed curve i s balance f o r the main i c e stream. The s o l i d curve i n c l u d e s mass contributions from the tributaries (0) through (5) i d e n t i f i e d i n F i g u r e 3.2. (b) G l a c i e r width. The broken curve i s the width of the main i c e stream o n l y . The s o l i d curve i n c l u d e s the width of the Hodgson G l a c i e r . (c) Bed topography. The i c e s u r f a c e p r o f i l e i n 1951 i s a l s o shown.  99 elevation  receives  about  2.2 metres  (ice  equivalent)  a c c u m u l a t i o n each y e a r ; the r e c o r d t h e r e i n d i c a t e s t h a t was  a  low  net  1965-66  s n o w f a l l y e a r by about 40% i n the Kaskawulsh r e g i o n  (Marcus and R a g l e , (1970, F i g u r e 7 ) ) . Keeler  (1969) r e p o r t e d t h a t , on Mount Logan, 60 km s o u t h of  the S t e e l e 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 above 2500 metres  elevation.  The  net  on p r e c i p i t a t i o n  .accumulation  is  about  0.8 metres of i c e a n n u a l l y . Stanley  (1969,  F i g u r e 1)  2400 metres on the S t e e l e photography.  Wood  located  Glacier  (1972)  put  the f i r n l i n e at about  based  the  on  firn  the  line  1951 in  aerial  the h i g h e r  e l e v a t i o n range of 2750 t o 2900 m e t r e s . T h i s would l e a v e no  a c c u m u l a t i o n a r e a . Sharp (1947) gave the e s t i m a t e of 8000 t o  9100 f e e t  (2440 t o 2775 m e t r e s ) .  The broken l i n e i n F i g u r e 3.3 balance  function  observations.  (unpublished)  The i c e f l u x  i s i n c l u d e d through a function.  In  (a) shows an  estimated  f o r the main i c e stream of the S t e e l e  consistent with C o l l i n s  and  with  the  6 b £ . The  perturbation  Appendix 19  solid  Glacier indirect  6b^ to  the  mass  balance  I d e s c r i b e two methods of e s t i m a t i n g  curve  in  F i g u r e 3.3  b a l a n c e w i t h the t r i b u t a r y terms 6 b ^ W h i l e F i g u r e 3.3 (a)  mass  from the i t h s m a l l t r i b u t a r y g l a c i e r  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 t o the  almost  (a)  estimate  shows the mass  included.  r e p r e s e n t s the best a v a i l a b l e e s t i m a t e  of the S t e e l e G l a c i e r mass b a l a n c e , i t s h o u l d be  kept  in  that  with  longterm  mass  balance  may  change  c l i m a t i c change. The decades data  for  the  Steele  significantly  1930-1950, upon which much  of  mind  the  G l a c i e r i s based, appear t o have been an  100 exceptionally  warm  period  (§_.£. Hansen  S c h n e i d e r and M e s i r o w , 1976, Chapter  3.3.5 Because  others,  1981;  3).  CYCLIC SURGE PATTERN FOR THE MODEL I am i n v e s t i g a t i n g consequences of s u r g i n g ,  than surge mechanisms, u5(x,t).  and  This  I specify a p r i o r i  aspect  of  surge  the  sliding  modelling  S e c t i o n 1 . 5 . 3 . I have chosen t o use a  is  sliding  rather  velocity  discussed  velocity  in  having  the form u (x,t) s because  it  can  = X(x,t) T(t)  represent  the  (3.3.5)  observed s l i d i n g of the S t e e l e  G l a c i e r r e a s o n a b l y w e l l . Other f u n c t i o n a l forms of u s ( x , t ) e q u a l l y w e l l f i t the o b s e r v a t i o n s (1972).  The  time  dependent  of  term  w e i g h t i n g f a c t o r between z e r o and the  form  quiescent gives  of  T(t)  for  T(t) unity.  (1969)  is  a  and  Wood  nondimensional  F i g u r e 3.4  (b)  shows  a surge c y c l e of l e n g t h t « . D u r i n g the  s t a g e , T has some s m a l l  constant  value  f  (£.£.  f=0  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 t o the peak v a l u e by time t 1 f t2,  Stanley  could  then  remains at the maximum u n t i l  f a l l s back t o the nonsurge l e v e l by time t 3 .  X ( x , t ) shown i n distribution  F i g u r e 3.4 of  the  (a)  gives  sliding  the  velocity  O b s e r v a t i o n s on the S t e e l e G l a c i e r ( S t a n l e y , surging Robin, sliding  glaciers, 1969;  and  Robin  starts  in  and a  some  The term  normalized at  a  spatial  time  (t-t0).  1969) and on  t h e o r e t i c a l work on s u r g i n g  Weertman,  1973)  suggest  that  s m a l l r e g i o n , and the b o u n d a r i e s of  other (e.g. rapid this  101  FIGURE 3 . 4 . S l i d i n g Model For S t e e l e G l a c i e r . (a) X ( x , t ) g i v e s the n o r m a l i z e d 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 w e i g h t i n g f u n c t i o n f o r 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 glacier. this  each  transition  distribution  surge  these  x,(t),  and  c , ( t ) , and c 2 ( t )  the  allow detailed for  x.,(t),  c.,(t), of  cosine.  of s l i d i n g can change d u r i n g the surge as  x.2(t),  Steele  four  move  respectively.  of the c - ( t ) .  velocities,  p o i n t s move a c c o r d i n g  x2(t)  at  to  so  that  I  use  The  the  four  velocities  The data from the  G l a c i e r are not s u f f i c i e n t l y  estimation  In  from a zone of r a p i d s l i d i n g t o a  zone of no s l i d i n g i s g i v e n by one h a l f c y c l e of a  c.2(t),  the  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 s u r g e .  model,  points  up)  complete to  constant  the v e l o c i t y  values  transition  1 02 x (t) 2  with s i m i l a r equations x.2(t0), trigger  x.,(t0),  = x (t 2  ) + (t-t  0  0  ) c  f o r the o t h e r t h r e e p o i n t s . The  x,(t0),  and x 2 ( t 0 ) d e f i n e the  by Raymond and o t h e r s  extent  of  the  (unpublished, Figure  on the V a r i e g a t e d G l a c i e r , A l a s k a , i n d i c a t e a sliding  velocity,  0.3 m d " 1  (110 m a " 1 )  glacier.  The  mid-1980's.  from in  0.05 m d " 1  1979  in  the  (18 m a " 1 ) upper  The  and i s expected simulations  to  it  increase  i n 1973,  reaches  surge  of  to the  in  the  which I have c a r r i e d out f o r  this  c h a p t e r do not i n c l u d e a p r e - s u r g e i n c r e a s e although  regular  9)  V a r i e g a t e d G l a c i e r appears t o have a surge p e r i o d  of about 20 y e a r s ,  probably  could  sometime  in  basal  computer  program.  However,  observed p r e - s u r g e v e l o c i t i e s  sliding,  be m o d e l l e d 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 the  constants  zone.  Observations  of  (3.3.6)  2  for  the  available  in  S t e e l e G l a c i e r , the  (Wood, 1972) are l o w , and are  not  s u f f i c i e n t l y d e t a i l e d t o 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  leading  to x 2 ( t ) )  edge of the zone of r a p i d s l i d i n g ( i . e .  i s a r e g i o n of compressive  flow.  The  instability  r e g i o n s of compressive  flow i s w i d e l y r e c o g n i z e d ( e . g .  1969,  nonsurging  p. 207).  For  r e s u l t i n g from compressive balanced  by  surface  glaciers,  surface  to  have  any  largely  appreciable  on s u r f a c e e l e v a t i o n ; t o 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 the  rise  A surge l a s t i n g o n l y one or two  y e a r s o c c u r s too q u i c k l y f o r m e l t i n g control  the  of  Paterson,  flow i n the a b l a t i o n a r e a i s  melting.  x,(t)  velocities  c,  and  sliding,  c 2 of the v e l o c i t y d i s t u r b a n c e must be  103 sufficiently velocity  greater  than the  material  velocity  that  the  change can move ahead of any bulge forming i n the i c e .  For the elementary case of i c e  with  sliding  t h i c k n e s s h 0 advancing i n t o t h i n n e r stagnant in  U0  velocity  U0  and  i c e of t h i c k n e s s h ,  a c h a n n e l of c o n s t a n t w i d t h , a s i m p l e c o n t i n u i t y argument i n  Appendix 18 shows t h a t when c , = c 2 , the v e l o c i t y d i s t u r b a n c e must move at l e a s t  at the speed  c  i  £ U  o  h -h o  (7. 7.  1  to p r e v e n t the development of a shock velocity  transition.  Equation  front  in  the  advancing  ( 3 . 3 . 7 ) can be used t o  the i c e t h i c k n e s s h 0 and h , ; f o r the S t e e l e G l a c i e r ,  7 )  estimate  I have used  i t o n l y to choose r e a s o n a b l e v a l u e s of c , and c 2 .  3.4 STEELE GLACIER MODEL 2  3.4.1  PROBLEMS WITH STEELE MODEL 1  I ran the S t e e l e G l a c i e r flow  law  parameters  parameters are g i v e n  (3.3.1) in  model 1 with  Table 3 . 1 .  ( F i g u r e 3.3)  no The  using  the  s l i d i n g . The n u m e r i c a l surface  profiles  at  50 y e a r 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 a r e shown i n Figure 3.5.  The  i n d i v i d u a l t r i b u t a r i e s c o a l e s c e d by 250 y e a r s .  The maximum v e l o c i t y 33 m a " 1 surface  at is  41 m a " 1 .  x=12.5 km.  (see The  (averaged over depth) The  S e c t i o n 1.4.4) final  steady  in  steady  state  is  c o r r e s p o n d i n g v e l o c i t y a t the i c e (n + 2 ) / ( n + l ) state  times  this,  j..e.  l e n g t h i s 35.5 km, and the  104 n  3  A - n -1 kPa s 5.3  10"15  5  0.8  g  Ax  fi  ms" 2  kg  9.8  900.  irr  TABLE 3 . 1 . Parameters f o r S t e e l e steady  3  At  m  a  500.  1 .0  state.  T  AOOOl  Steele Glacier Growth to steady state No sliding Profiles every 50 years 3000  2000h  1000 0  10  20  30 x (km)  FIGURE 3 . 5 . S t e e l e Model 1 Growth To Steady S t a t e . There was no s l i d i n g . The model parameters a r e g i v e n i n T a b l e 3 . 1 . The i c e s u r f a c e p r o f i l e s are shown at 50 year intervals 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 l o b e s 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 d e p t h i s 444 m at x=25 km. The  s u r f a c e p r o f i l e s from Model 1 appear t o 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  width functions i n Figure 3.3. steady  state  F i g u r e 3.6  streamlines  However, I c a l c u l a t e d some of  for  Model 1;  (note t h a t the a b s c i s s a i s  elevation).  The  vertical  1  mesh  these  are  shown  ice thickness rather  increment  for  the  and the in than  velocity  r  1  Steele Glacier Model 1 Steady streamlines No sliding Arrows every 100 years  FIGURE 3 . 6 . Steady S t a t e S t r e a m l i n e s For Model 1. This cross-section shows ice depth. Arrows on the streamlines i n d i c a t e each 100 y e a r s 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 t o the smooth s t r e a m l i n e s are due t o 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 t e r m s .  calculations  (see  F i g u r e 2 . 8 ) was DZ = 1 5 m  The p e r t u r b a t i o n s of  tributary  ice  i n the s t r e a m l i n e s  are caused by the  f l u x t h r o u g h the mass balance  looking for i r r e g u l a r i t i e s  in streamline  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 ;  shape  t o have  additions  f u n c t i o n . I am  and  in  features  isotope of  this  106  3 £  2  co  1  o c a a  0  CO.  -1  u) a  •2  3000h £  2000h  •g £ 1000 3000  Steele Glacier Model 2 Growth to steady state No sliding Profiles every 50 years 2000 CXI X  1000  0  10  20  (km)  30  AO  50  FIGURE 3 . 7 . Model 2 For S t e e l e G l a c i e r . (a) Mass b a l a n c e . (b) G l a c i e r w i d t h . (c) Bed topography and i c e s u r f a c e 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 t o s t e a d y s t a t e w i t h no s l i d i n g .  form  introduced  through  the  mass  balance  function  is  107 undesirable. rather if  Adding the t r i b u t a r y f l u x at  than at the c h a n n e l margins  only  the  trajectories  ice  surface  are  also  shape  the  glacier  surface  i s an adequate a p p r o x i m a t i o n is  desired,  desired. this  If  particle  approximation  is  unacceptable. 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 added.  The  width  representation distribution  of of  function the  valley  the  obtained  velocities.  of  3.3  (b))  Steele  is  Creek.  only  resolution;  Wood  spatially  and  (1972)  surface  configurations.  where the c h a n n e l w i d t h g r a d i e n t thicken  rapidly  d e t a i l e d d a t a on how  It  For  the  glacier  actually  width  function  averaged velocity in in  some  regions the  ice  surface slope.  More  changes  speed  to  f o r the 1966-67 s u r g e . which  has  the  same  d e t a i l as the s l i d i n g d a t a .  When a network computer model f o r t r i b u t a r i e s 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  amount of d e t a i l  is  Stanley  result  instance,  and o b t a i n a r e v e r s e d  i s n e c e s s a r y t o use a  when  spatial  and  i s l a r g e and n e g a t i v e ,  t h i s s i t u a t i o n are not a v a i l a b l e  degree of s p a t i a l  and  The  F o r c i n g the i c e t o s l i d e at n e a r l y c o n s t a n t  unrealistic  prevent  detailed  temporally  through a c h a n n e l of h i g h l y v a r i a b l e w i d t h can  can  a  is  s u r g i n g v e l o c i t y of the S t e e l e G l a c i e r  not known w i t h the same (1969)  (Figure  sliding  i n Model 1 can be  justified.  is  developed,  are a v a i l a b l e ,  the  108 3.4.2  SIMPLIFICATIONS  F i g u r e 3.7 Steele  (a)  Glacier.  features, previous  yet  and (b)  This  show  model  avoids  the  a  simplified  resembles  model  for  Model 1 ' i n  difficulties  I  its  described  the gross  in  the  section.  Figure 3.7(c)  shows  the growth of Model 2 t o steady  state  w i t h no s l i d i n g , u s i n g 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  the  i c e t h i c k n e s s i s 452 m at x=21  maximum  velocity  (averaged over depth)  v a l u e s are c l o s e t o the v a l u e s  i s 30.4 m a " 1  length  is  36.5 km,  km, and the maximum at  x=11  f o r Model 1 (see  km.  Section  These  3.4.1).  3.5 STABLE ISOTOPES IN GLACIOLOGY  3.5.1 The  DEFINITION OF THE DEL SCALE  s t a n d a r d method of d e s c r i b i n g the i s o t o p i c c o m p o s i t i o n  of oxygen and hydrogen i n water i s the DEL s c a l e ( 6 ) . R  The  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 s o t o p e s  and D/H can be measured w i t h a mass  b18/016  a  practical  c o n c e n t r a t i o n s c a l e s h o u l d be based on R. The 6 v a l u e  ( 3 . 5 . 1 ) of  an  ice  sample  sample and R  is  the  relative  of a r e f e r e n c e  Ocean Water (SMOW) ( C r a i g ,  spectrometer;  ratio  difference  sample  1961).  known  between R^ of as  Standard  the Mean  109  6  A  drawback  available. U.S.  of  p a s t by the  Isotope  true  at  Standards set  This  I c e f i e l d Ranges  and  the / 0  8  1  fact  spectrometer reported  Dansgaard's  that  isotopic  6  Energy  up  Intercalibration a standard  difference,  (West  s t a n d a r d s from  Agency,  laboratory  O  o  in  is  several  one  glacier.  are the  i n the  Vienna.  In  M e e t i n g on S t a b l e Hydrology  and  SMOW  is  for  my  of i s o t o p i c d a t a f o r  the  not  Vienna  significant  1972;  West,  unpublished;  t h i s change. Dansgaard  ± 0 . 1 2 ° /  measurements ± 0 . 0 3 ° /  1 9 6 1 ) and  Krouse,  predate  in  of  o  o  (per  mil)  for  ( 1 9 6 9 )  routine  6 . Ahern ( u n p u b l i s h e d  f o r samples from the S t e e l e Copenhagen  also  l e v e l . T h i s i s adequate f o r work on g l a c i e r s , by  no samples  sample c a l l e d V i e n n a SMOW. The  The r e p o r t s  and  u n p u b l i s h e d [b])  1 5 8 )  1  ( 3 . 5 . 1 )  V i e n n a , the C o n s u l t a n t s '  r e p o r t e d a r e p r o d u c i b i l i t y of  p.  3  SMOW -I is  0  1 0  Atomic  glaciological applications.  mass  SMOW  between t r u e SMOW ( C r a i g ,  - 0 . 0 5 ° / o o -  Ahern,  R  SMOW  International  Geochemistry  - R  Bureau of Standards have been d i s t r i b u t e d  1976,  difference  S  Samples of o t h e r  National  September  R  =  Glacier.  now a c h i e v e s since 6  [b],  may  this vary  p a r t s ° / 0 0 t o t e n s of p a r t s ° / o f o r samples from any 0  1 10 3.5.2  FACTORS AFFECTING DEL  The nonzero v a l u e of DEL (6) glacier  is  the  result  of  for  a  ice  sample  from  left  the  well-mixed  ocean  6 i s c l o s e t o z e r o ) . At 0 ° C , the vapour p r e s s u r e s of  t h r e e major i s o t o p i c forms of water have the approximate (e.g.  Dansgaard, 2  the  dependent  : 0.989 : 0.904  2  differences  resulting  the  ratios  1964)  H O 1 6 : H O 1 8 : HDO = 1.000 and  a  l o n g s e r i e s of p r o c e s s e s i n the  h y d r o l o g i c a l c y c l e s i n c e the water (where  an  (3.5.2)  i n c r e a s e w i t h d e c r e a s i n g t e m p e r a t u r e . The  differences  in  volatility  lead  to  temperature  i s o t o p i c f r a c t i o n a t i o n i n e v a p o r a t i o n and c o n d e n s a t i o n  processes.  Under  fast  evaporation  or c o n d e n s a t i o n 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 concentrations  of  a  particular  isotopic  (the  ratio  species  of  the  i n the two  phases) i s c o m p l i c a t e d . The p r o c e s s which tends t o c o n t r o l the 6 values from  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 c o n d e n s a t i o n of cloud  vapour;  fortunately,  this  c a n , i n most c a s e s ,  a d e q u a t e l y m o d e l l e d by a R a y l e i g h c o n d e n s a t i o n slow  condensation  (quasi-equilibrium  phases) w i t h immediate removal 1964). factor  For  of  process,  .i.e. a  condensate  (Dansgaard,  slow c o n d e n s a t i o n 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  species  of  water  at  r a t i o s are w e l l known above 0°C have  be  of the vapour and l i q u i d  i s j u s t a r a t i o of the vapour p r e s s u r e s of the  isotopic  and  the  droplets  different  the ambient t e m p e r a t u r e . laboratory  measurements,  been e x t r a p o l a t e d t o - 2 0 ° C (Dansgaard,  1964, T a b l e 1)  u s i n g a formula of Zhavoronkov  from  These  and  others,  (1955).  When  the  111 Rayleigh  condensation  model  precipitation  i s p r i m a r i l y an  temperature.  In  general  i s a p p l i c a b l e , the 6 v a l u e of indication  of  site,  to  be  i n w i n t e r than i n summer. A c o n t i n e n t a l e f f e c t  sometimes  observed  precipitation  the  (Dansgaard,  1964);  the  6  more  is  also  values  of  at c o n s t a n t c o n d e n s a t i o n tejnperature may decrease  with distance from  condensation  t e r m s , 6 v a l u e s tend t o 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, a t any one negative  the  the  from the ocean due t o d e p l e t i o n of heavy  storm  systems  through  precipitation,  isotopes  and  due  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 f r e s h w a t e r  to  sources.  F a c t o r s o t h e r than temperature can i n f l u e n c e the 6 v a l u e of precipitation. evaporation  Dansgaard  from  (1964)  falling  droplets,  drops and a i r through which changes,  and  cause  they  variations  c o m p o s i t i o n of the source variations  discussed  in  effects  fall,  n o n - e q u i l i b r i u m phase frequency  While  and  these  year.  i n 6 from storm t o s t o r m , t h e i r e f f e c t  6(018/016)  and  6(D/H)  R a y l e i g h c o n d i t i o n s (Dansgaard, of  6(018/016)  isotopic  processes  average summer or w i n t e r 6 v a l u e tends t o be c o n s t a n t to  are  1964);  of  i s o t o p i c exchange between  the  storms.  the  linearly  can  on the  from  year  r e l a t e d under  simultaneous  measurement  and 6(D/H) can be used t o r e v e a l the presence  of  non-equilibrium condensation. • S e v e r a l p r o c e s s e s i n snow and f i r n tend t o distribution,  individual  storms,  difference.  In r e g i o n s 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  sometimes  percolating  the  meltwater  v e r t i c a l column of snowpack t o the  average  differences  the  isotopic  an"d  obliterating  homogenize  summer  can  bring  6-value  to  between winter  the whole (Dansgaard  1 12 and  others,  always  1973).  simple;  percolating  However,  Ahern  meltwater  the e f f e c t s of meltwater are not  (unpublished [ a ] , could  enhance  p.  109)  rather  found  than  that  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 . r e g i o n s w i t h no summer m e l t i n g , some smoothing of profile  occurs  v e r t i c a l gradients  of temperature  frequent  others,  barometric  pressure  in  i n the f i r n  l a r g e s e a s o n a l temperature v a r i a t i o n s ) , with  isotopic  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 high  the  In  or  in  regions  ( e . £ . due t o  stormy  changes  of  regions  (Dansgaard,  and  1973).  Below  the  550 kg n r 3 ,  depth  the " vapour  isotopic profile is process  at  which spaces  smoothed  firn in  only  reaches  a  density  of  the f i r n are i s o l a t e d . The by  solid  diffusion.  i s too slow t o have an a p p r e c i a b l e e f f e c t  This  on i c e i n the  Steele G l a c i e r .  3.5.3  PREVIOUS ISOTOPIC STUDIES  Assuming (1) similar  tracks  year t o y e a r ,  that  the  precipitating  (5)  the  follow  (2) t h a t R a y l e i g h c o n d e n s a t i o n  occurs,  (3)  that  can be s i m p l y r e l a t e d ,  (4)  6-temperature r e l a t i o n s h i p i s c o n s t a n t w i t h t i m e , and  t h a t the flow p a t t e r n of the i c e  then 6 v a l u e s of  masses  w i t h s i m i l a r frequency the year around and from  s u r f a c e and c o n d e n s a t i o n temperatures that  air  mass  can  be  calculated,  i n i c e c o r e s can be r e l a t e d t o c l i m a t e at the time  precipitation.  This  was  first  pointed  out  by Dansgaard  ( 1 9 5 4 ) . A thorough d i s c u s s i o n of i c e c o r e s t u d i e s can  be  found  1 13 i n Chapter 15 of P a t e r s o n  (1981).  The G r e e n l a n d i c e sheet p r o v i d e s the most s u i t a b l e and  meteorological  conditions  i n t e r p r e t a t i o n of a deep i c e c o r e The  first  major  drilling  S.I.P.R.E.  (U.S.  Army  Establishment,  now  for  called  simple  (Dansgaard  program  Snow,  a  Ice,  CRREL,  flow  climatic  and o t h e r s ,  1973).  was  undertaken i n 1956 by  and  Permafrost  Cold  Regions  Research  Research and  E n g i n e e r i n g L a b o r a t o r y ) ; a 411 m c o r e was r e c o v e r e d in  ice  at  Site  2,  northwest G r e e n l a n d . T h i s was f o l l o w e d by deep c o r e s at Camp  Century, Greenland Antarctica,  in  in  1968  1966  (1387 m),  (2164 m).  and  1969,  1971).  variations but  and J o h n s e n ,  is  Byrd  Station,  The Camp Century c o r e has  used t o 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 (Dansgaard  at  I969[a],  I969[b];  past  100,000  more  difficult  v a r i a t i o n s of 6 were not  t o date a b s o l u t e l y , preserved  (Robin,  climate  1970; Johnsen and o t h e r s ,  during  p r o c e s s . The l e n g t h of the f l o w l i n e , h o l e are i n d i s p u t e  years  Dansgaard and o t h e r s ,  The Byrd S t a t i o n c o r e a l s o shows l o n g term  ( E p s t e i n and o t h e r s ,  been  1972),  because the annual the  ice  formation  and the time s c a l e f o r  this  1977).  C o r i n g 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 G r e e n l a n d (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 , 1975,  1977),  (Dansgaard Merlivat,  Dome C ( L o r i u s and o t h e r s , and  others,  1977),  1977),  and  1979), Terre  Little  1974) ,  Devon  Ice  Clarke,  1978; F i s h e r ,  and  including  1973; Koerner and P a t e r s o n ,  Cap ( P a t e r s o n and o t h e r s , 1979),  America V  A d e l i e ( L o r i u s and  and at s i t e s i n the Canadian A r c t i c ,  Meighen Ice Cap (Koerner and o t h e r s ,  1974,  Agassiz  1977; P a t e r s o n and  Ice  Cap,  Ellesmere  11 4 Island  (D. F i s h e r ,  personal  communication).  Ice  cores  i s o t o p i c a n a l y s i s have been o b t a i n e d from the p l a t e a u at  for  5400 m  on Mount L o g a n , Yukon T e r r i t o r y by G . H o l d s w o r t h . During  p e r i o d s of e x t e n s i v e g l a c i a t i o n , deep sea  are e n r i c h e d i n O 1 8 ; the i s o t o p i c c o m p o s i t i o n of  sea  a l t e r e d because of the l a r g e volume of 0 1 8 - d e p l e t e d Measurements (e.g.  of  water  of deep i c e  1976)  have  complemented  sediments  the  climatic  cores.  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 c o r e s also  been  supported by o t h e r s t u d i e s .  (1977) found t h a t  the  temperature  isotopic  was  compatible  records  d i s t r i b u t i o n of temperature (1978)  used  the  history with  derived  the  in boreholes.  present  Paterson  linear  and o t h e r s ,  relationship  (i960)  and  between  flow model f o r  Picciotto  and  6(D/H)  and  others,  (1968),  to  measure  the  and  L o r i u s and o t h e r s ,  6(D/H)  in Antarctic precipitation.  ten  6  metre and  the  of the  in  the  St.  Elias  firn  L o r i u s and near-surface Antarctic  (1969) measured r e g i o n a l v a r i a t i o n s of  West and Krouse (1972) measured i s o t o p i c r a t i o s at in  on  6  recent accumulation rate at  sites.  sites  Clarke  demonstrated the e x i s t e n c e  (1970) used the a n n u a l v a r i a t i o n of  samples  vertical  E a s t A n t a r c t i c a , and L o r i u s and Delmas (1975) found a  temperatures. others,  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  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 of  from  boreholes.  Picciotto  coast  has  Robin (1976) and Johnsen  boundary c o n d i t i o n f o r a time-dependent heat Devon I s l a n d  is  i c e on l a n d .  the i s o t o p i c c o m p o s i t i o n of deep sea  Hayes and o t h e r s ,  studies  sediments  M o u n t a i n s , Yukon T e r r i t o r y ,  several including  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 , balance  estimates  patterns. small  surge-type  glacier  increase  demonstrating  mass  and r e l a t i n g i s o t o p i c c o m p o s i t i o n t o weather  A l o n g i t u d i n a l s u r f a c e s a m p l i n g of  systematic  obtaining  that  in  of 6  Rusty  Glacier  (a  the S t e e l e Creek b a s i n ) showed a  with  height  the g e n e r a l  in  the  ablation  i s o t o p i c p a t t e r n of the  zone, glacier  was not d e s t r o y e d by s u r g i n g . Ahern  (unpublished [b],  p.  162)  obtained  p r o f i l e t o a depth of 36 m i n a b o r e h o l e at x=15 on  the  Steele  an  isotopic  km ( F i g u r e  3.2)  G l a c i e r . T h i s p r o f i l e appeared t o 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°/0o  in  6(018/016). Isotopic  studies  on  glaciers  and i c e s h e e t s are  r e l a t e d t o s t u d i e s of i c e f l o w . The determined,  and  annual  layers  age  of  cannot  the  ice  always  i s o t o p i c a l l y . Time s c a l e s f o r i c e c o r e s  have  Dansgaard  Philberth  and  Johnsen  (I969[a]),  (1971), and by Hammer and o t h e r s steady  (1978),  (1977)  boreholes;  this  rate  or  measured allowed  temperature the them  deriving  and  l a r g e d i s t a n c e s from i c e the  time  time  t w o - d i m e n s i o n a l flow s i m u l a t i o n s , (1971,  used t o d a t e the  p.  148)  ice.  for  by  Federer  assumptions  gradient.  eliminate  boreholes at  influences  a  to  when  others  derived  of  forms of  Paterson  and  v e r t i c a l d e f o r m a t i o n i n the Devon  assumptions  velocity  using  be  detected  s t a t e , p r o x i m i t y t o an i c e d i v i d e , and s p e c i f i c  vertical strain others  by  must  be  been  closely  such  scale  the from  divides, scale as  strain  the f l o w . For the  horizontal  calculation, given  rate  by  Budd  and and  the B y r d s t a t i o n f l o w l i n e , must be  1 16 A conceptual  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  s t a t e flow models;  the e x i s t e n c e  of  which the i s o t o p i c p r o f i l e s attempt the  existence  surface  essential  For  may this  isotopic  trajectories ice  vary  in  reason,  the  absence  the e f f e c t  profile  in  flowline  Budd and Mclnnes (1978,  Vostok  Vostok  to  interested  flow models  interpretation; using  the  are  Jenssen  on  the  ice  Antarctic periodic  i c e s u r f a c e p r o f i l e s were computed by  1979). T h i s i s the o n l y p r e v i o u s work of  effect  of  in t h i s chapter,  isotopic  data.  c o r e . He c a l c u l a t e d  and i s o t o p i c d i s t r i b u t i o n  i n a t i m e - v a r y i n g i c e mass.  in  climatic  M i r n y , assuming  which I am .aware i n which t r a j e c t o r i e s are c a l c u l a t e d  ice  of i c e sheet e l e v a t i o n changes on  the  from  s u r g e s ; the time-dependent  of  of  time-dependent  and s i m u l a t e d i s o t o p i c p r o f i l e s  sheet  variations,  to determine, could preclude  t o 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  (1978) i n v e s t i g a t e d the  climatic  of steady s t a t e c o n d i t i o n s . In a d d i t i o n , the  elevation  variations.  the  steady  Jenssen  surging  on  I investigate  information  to  reveal  (1978) the  the  was  climatic  feasibility  the surge h i s t o r y .  J e n s s e n (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  equation  (3.5.3)  c a l c u l a t i o n s were not  3.5.4  below.  The  details  of  the  flowline  described.  DEL RELATIONS FOR THE MODEL  In the computer model, I assume t h a t : (1)  all given  precipitation elevation  has  the  or  location,  h o m o g e n i z a t i o n i n the f i r n ) ,  annual  and  average v a l u e f o r (!•£_• r a p i d  the  isotopic  1 17 (2)  isotopic diffusion I  use  two  6-precipitation  extremes of t o p o g r a p h i c which  I  call  i s n e g l i g i b l e i n the s o l i d  control  Model 61,  models, of  is  ice.  representing  precipitation.  based  opposite  The  on the assumption  first, that  the  i s o t o p i c c o m p o s i t i o n of s n o w f a l l i s c o m p l e t e l y c o n t r o l l e d by the s u r f a c e e l e v a t i o n of the g l a c i e r h ( x , t ) ,  and can be  by  range  a  linear  considered,  relationship  over  the  of  elevations  i.e. 6(x,t) = 6  Values  approximated  of  the  constants  6  + c h(x,t)  0  and  0  c  (3.5.3)  consistent  measurements i n the S t e e l e Creek b a s i n  (West and  w i t h the few Krouse,  and on the S t e e l e G l a c i e r ( A h e r n , u n p u b l i s h e d [b]) 6  = -15.0  0  are  °/oo  c = - 0 . 0 0 4 m- 1 % o The g r a d i e n t Sharp  and  U.S.A.  It  (3.5.4)  c i s c l o s e t o the v a l u e of - 0 . 0 0 5 m " 1 °/oo found by others, is  - 0 . 0 0 2 nr 1 ° /  (1960)  within  0 0  used c = - 0 . 0 0 6 irr flowline.  1972)  a  reported 1  Isotopic  for factor by  °/oo i n a values  c at Blue G l a c i e r , Washington, of  two  Dansgaard modelling on  d e s c r i b e d by Model 61, because  ice the  of  the  value  (1961). Jenssen study  of  an  of  (1978)  Antarctic  s h e e t s and i c e caps can be ice  sheet  surface  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,  called  Model 62,  assumes  that  the  i s o t o p i c c o m p o s i t i o n of the s n o w f a l l i s c o m p l e t e l y 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 G u l f of A l a s k a ; 6 i s assumed t o be a f u n c t i o n o n l y of p o s i t i o n x . A p p r o x i m a t i n g the  1951  (Topographical  by  Survey,  1969)  Steele a  Glacier  straight  line  surface in  the  1 18  FIGURE 3 . 8 . Reference S u r f a c e For 6-x F u n c t i o n . The a c c u m u l a t i o n a r e a of the S t e e l e G l a c i e r i s assumed to extend (see F i g u r e 3.7 (a)) t o x=13 km ( v e r t i c a l broken line). The c u r v e d l i n e is the observed ice surface elevation i n 1951. A p p l y i n g the 6 - e l e v a t i o n relation ( 3 . 5 . 3 ) t o the l i n e a r a p p r o x i m a t i o n (dashed l i n e ) to this s u r f a c e , and u s i n g the c o n s t a n t s ( 3 . 5 . 4 ) g i v e s 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 c o n s t a n t s i n ( 3 . 5 . 6 ) .  a c c u m u l a t i o n a r e a , as shown i n F i g u r e 3 . 8 , r e l a t i o n s h i p to t h i s l i n e , x-6  and a p p l y i n g the  6-h  w i t h the c o n s t a n t s ( 3 . 5 . 4 ) , g i v e s  the  relationship 6(x,t) = 60 + k x  (3.5.5)  w i t h the c o n s t a n t s 6o  =  —  26.6  °/oo  k = 0.00022 m- 1 ° / o o Model 62  is  more  glaciers.  The S t e e l e G l a c i e r ,  probably f a l l s Because  appropriate  than being  (3.5.6) Model 61 f o r s m a l l  a  large  valley  glacier,  between the two e x t r e m e s . Model 61  allows  the  s n o w f a l l t o v a r y w i t h both p o s i t i o n  isotopic and  time,  composition while  or  discontinuities  of  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 variations  valley  larger  i n 6 ( 0 1 8 / 0 1 6 ) . Models 61 and 62  119 can be c o n s i d e r e d t o g i v e respectively  to  the  the  isotopic  maximum  and  minimum  distribution  structure  w i t h i n the S t e e l e  Glacier.  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 t e e l e G l a c i e r i s unknown. Section 3.6.2,  I  Glacier  surge  using  150 y e a r s  In  show t h r e e computer s i m u l a t i o n s of the S t e e l e  suggested  periods by  field  spanning  the  range  observations.  I  of  50  to  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 t h r e e c a s e s . The p e r i o d (97 y e a r s ) which gave the most r e a s o n a b l e i c e t h i c k n e s s  profiles  at  isotopic  all  times  calculations. for  ice  was In  selected  for  trajectory  and  Section 3 . 6 . 3 , I present t y p i c a l  particles  in  this  model.  These  trajectories  trajectories  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 To  o b t a i n a p e r i o d i c a l l y r e p e a t i n g surge c y c l e , I used the  constants sliding  in  Table 3.2  pattern  Figure 3.9. periods;  PERIODICALLY REPEATING STATE  shown  at  intervals  of  0.25 y e a r s  The surge d u r a t i o n was two y e a r s . I used t h r e e  47,  97,  periods described steady  is  t o generate the s l i d i n g v e l o c i t y . T h i s  state  in  and  147 y e a r s  span  the  range of  in surge  possible  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 Figure 3.6,  I ran the S t e e l e G l a c i e r Model 2  t h r o u g h 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  time,  the  120 (a) (a) (a) (a) (a)  to t, t2 t3 t« f  U0  (m  x. 2 < t )  0.0 0.75 1.5 2.0 1 ) 47. 0 2) 97. 0 3) 147. 0 0.0 5000.  a-1)  0  X . Xi X2  i(t0)  (t0) (t0)  c .:(t0) c. i ( t ) C i (t0) c2 (t0)  (km) (km) (km) (km) (m (m (m (m  0  8.0 18.0 19.0 26.0  aa' a' a"  ) ) 1 ) 1 ) •  -1000.0 7500.0 15000.0 15000.0  1 1  TABLE 3 . 2 . V e l o c i t y p a t t e r n f o r S t e e l e s u r g e . The t h r e e v a l u e s of surge p e r i o d t « were used 3 . 1 0 , 3 . 1 1 , and 3.12 r e s p e c t i v e l y . model no l o n g e r "remembered" the i n i t i a l  for  Figures  steady s t a t e c o n d i t i o n ;  i t r e p e a t e d the same s u r f a c e p r o f i l e sequence t o one p a r t with  each  new  c y c l e f o r 97 and 147 year p e r i o d s ,  p a r t s i n 10 3 " w i t h a  47 year  period.  The  time  computer model must be very s m a l l when the g l a c i e r order  to  maintain  accuracy  summarizes  the n u m e r i c a l time  numerical  and  physical  Table 3 . 1 . F i g u r e s surface  (solid  l i n e ) at t 3  than  3.10)  step  constants  3.10 t h r o u g h 3.12  l i n e ) at t 0 ,  and t o a few  steps  sequence  used.  had  values  Table The  shown  show the p r e - s u r g e  3.3  other in  glacier  and the p o s t - s u r g e s u r f a c e  (broken  The  model  with  a  has a p r e - s u r g e p r o f i l e  47 year  surge  period  ( s o l i d l i n e ) which i s  less  100 m t h i c k beyond x=20 km. The p o s t - s u r g e p r o f i l e i s  less  than 100 m t h i c k beyond x=25 km. I t seems u n l i k e l y t h a t a of  the  is surging in  Section A1.5).  the  in  f o r the t h r e e surge p e r i o d s of 47, 97, and 147 y e a r s  respectively. (figure  (see  i n 10'  ice  surge  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 conclusion  shown i n t h i s diagram are u n r e a l i s t i c .  indicates  This  t h a t the S t e e l e G l a c i e r a c c u m u l a t i o n area  cannot p r o v i d e s u f f i c i e n t mass  in  just  47 y e a r s  to  generate  121  T  8000  Steele Glacier Model 2 Surge period 97 year* Sliding velocity at intervals of 0.25 years  6000 t=1.0  o  "54000 >  cn c CO  2000  0  10  20  30  x  40  (km)  FIGURE 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 M o d e 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 y e a r s d u r i n g the surge of two y e a r s d u r a t i o n . The end of each c u r v e i n d i c a t e s the p o s i t i o n of the a d v a n c i n g t e r m i n u s ( f o r 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 must  be  substantially  p e r i o d of times.  longer  147 y e a r s ( F i g u r e 3.12)  than 47 y e a r