A THERMAL INSTABILITY MECHANISM FOR GLACIER SURGES BY JOSEPH WALTER HOFFMANN B.A., Washington University (Missouri), 1970 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE i n the Department of GEOPHYSICS We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA March, 1972 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f the requirements f o r an advanced degree a t the U n i v e r s i t y o f 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 i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r 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 c o p y i n g o f 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 g r a n t e d by the Head o f my Department or by h i s r e p r e s e n t a t i v e s . I t i s understood t h a t c o p y i n g or p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a llowed without my w r i t t e n p e r m i s s i o n . Department o f GEOPHYSICS The U n i v e r s i t y o f B r i t i s h Columbia Vancouver 8, Canada Date March 24, 19 72 (i) ABSTRACT A theory of.thermal conduction i n s t a b i l i t i e s that could account for many of the properties of surging glaciers is presented. A one-dimensional diffusion equation with advection is solved numerically using a variation of the Crank-Nicholson method. Geothermal and f r i c t i o n a l heating are shown todominate the melting of ice from the base of the glacier, and the water that i s formed is assumed to flow in an approximately uniform sheet between glacier and bedrock. The sliding theory developed by Weertman is used to determine the effects of parameters such as geothermal heat flow, basal shear stress, and bedrock roughness on the behavior of surging glaciers. A cyclic evolution of the temperature profile near the bed of a surging glacier is discovered and is thought to be a principal factor in starting the surge and in determining the surge period. Some good correlation with observation is obtained when the theory i s applied to models which resemble known surging glaciers. An attempt to predict the future surge behavior of the Rusty Glacier, Yukon (previously known as the "Fox" Glacier) from the temperature data of Classen and Clarke produced uncertain results. ( i i ) TABLE OF CONTENTS ABSTRACT (*) CONTENTS ( i : L) LIST OF TABLES ( i v) LIST OF FIGURES ( v) ACKNOWLEDGMENTS ( v i) CHAPTER 1 Introduction 1 1.1 The Rusty G l a c i e r Project 1 1.2 G l a c i e r Surges 4 1.3 Summary of Previous Theories 9 CHAPTER 2 Theory 1 4 2.1 Gl a c i e r S l i d i n g 1 6 2.1.1 S l i d i n g Theory 16 2.1.2 Water Flow at the Base of a 2 2 G l a c i e r 2.2 Heat Conduction 2 8 2.2.1 Temperatures i n Glaciers 2 8 2.2.2 Mathematical Theory f o r the 30 Conduction of Heat i n Solids 2.2.3 Accumulation, Ablation, and Advection 32 2.2.4 Heat Sources at the G l a c i e r Bed 36 2.2.5 Internal Heat Generation 4 2 .2.3 The Surge Cycle 4 5 ( i i i ) CHAPTER 3 Results 3.1 The Modeling Procedure 3.2 Numerical Methods 3.3 Results f o r Surging Glaciers 3.3.1 Model.Glaciers 3.3.2 Rusty. G l a c i e r CHAPTER 4 Summary and Discussion BIBLIOGRAPHY APPENDIX: L i s t i n g s of Computer Programs (iv) LIST OF TABLES TABLE I Comparison of the s t r a i n heating and f r i c t i o n a l heating f o r surging g l a c i e r s . TABLE II Results f o r Model G l a c i e r No. 103. TABLE III Results f o r Model G l a c i e r No. 104. TABLE IV Results f o r Model G l a c i e r No. 105. (v) LIST OF FIGURES FIGURE 1 G l a c i e r surface and bedrock p r o f i l e s f o r Robin's n stress i n s t a b i l i t y theory. FIGURE 2 T o t a l melted i c e thickness versus rate of water 25 formation f o r several values of the combined g l a c i e r surface area. FIGURE 3 Flow of a water f i l m of uniform thickness beneath 26 a g l a c i e r . FIGURE 4 . . The boundary conditions f o r a g l a c i e r i n the s o l u t i o n 33 of the one-dimensional heat conduction equation. FIGURE 5 The replacement of the moving boundary temperature due 37 to advection during a surge by a var i a b l e temperature at a f i x e d point within the g l a c i e r . FIGURE 6 I l l u s t r a t i o n of the pre-surge temperature p r o f i l e . 46 FIGURE .7 . .. I l l u s t r a t i o n of the post-surge temperature p r o f i l e . 48 FIGURE 8 . Temperature p r o f i l e f o r the Rusty G l a c i e r , Yukon 62 (vi) ACKNOWLEDGMENTS I thank my supervisor, Dr. Garry Clarke, for a l l of h i s e x c e l l e n t guidance and e n t h u s i a s t i c i n t e r e s t during the course of t h i s work. I also thank Dr. R. D. R u s s e l l and the Department of Geophysics for the f i n a n c i a l support I received while 1 was completing the t h e s i s . I am also g r a t e f u l to Dr..J. Weertman of Northwestern Uni v e r s i t y , I l l i n o i s f o r the opportunity to read a preliminary d r a f t of his paper on s u b g l a c i a l hydrology, and to Dr. Stan Paterson and Dr. Ron Goodman, both from Ottawa, for useful discussions during t h e i r v i s i t s to Vancouver. Frequent use was made of the University's Computing Centre, and I g r a t e f u l l y acknowledge the services rendered by i t s s t a f f . . . The. computing costs and other expenses were financed by grant.#67-4327 from the. National Research Council of Canada to Dr. Clarke. I also thank Mrs. Corinne McAdam for a rapid and accurate performance i n the typing of the manuscript. F i n a l l y , and e s p e c i a l l y , I thank a l l of my friends among the f a c u l t y , s t a f f , and students of the Department, for each one.has contributed much to whatever I have achieved at this u n i v e r s i t y . CHAPTER 1 INTRODUCTION This thesis i s concerned wita developing a theory for g l a c i e r surges based on thermal fl u c t u a t i o n s at the i c e -rock i n t e r f a c e and incl u d i n g such parameters as i c e thickness, surface slope, geothermal heat flow, bedrock roughness, and mean annual accumulations and surface temperatures. The work was developed as one part of a study i n t o the subject of surging g l a c i e r s that began four years ago on the Rusty G l a c i e r , Yukon T e r r i t o r y , Canada. 1.1 The Rusty G l a c i e r Project The Rusty G l a c i e r (previously known as the "Fox" Glacier) i s a small v a l l e y g l a c i e r which l i e s west of the Steele G l a c i e r and north of the Hodgson G l a c i e r on the i n t e r i o r side o f the I c e f i e l d Ranges. Examination of the mcrainal patterns i n the va l l e y revealed a h i s t o r y of past advances of the kind generally associated with surging g l a c i e r s . Further i n d i c a t i o n s that another surge might occur i n the near future resulted i n the beginning o f an extensive and long-term.research program on the g l a c i e r . Dr. Garry Clarke and h i s students from the University o f B r i t i s h Columbia were among the p r i n c i p a l i n v e s t i g a t o r s f o r the project. During the summer of 1968 Clarke and David Crossley 2 undertook the mapping of the thickness of the glacier [Crossley, 1969]. The standard seismic reflection technique was abandoned when they discovered that the glacier thickness did not exceed 150 meters and was therefore too shallow for reliable seismic record interpretation. Successful determinations of ice depths were obtained however by conducting a detailed gravity survey over the glacier surface. Crossley also undertook a limited near-surface temperature program which identified the Rusty as a sub-polar glacier. The next summer Clarke and David Classen conducted a program of thermal d r i l l i n g and deep ice-temperature measurements on the glacier [Classen, 1970]. Their d r i l l i n g results verified the ice thickness values obtained by Crossley and Clarke the previous year, and the temperature measurements confirmed that the Rusty was below the pressure-melting point throughout. Furthermore they discovered two fascinating facts about the Rusty's thermal properties. First the observed temperature profiles were not linear near the glacier base but instead the gradient values decreased with depth. They believed that this unusual curvature was possibly caused by the glacier's surge history. Second they calculated a geothermal heat flow from the basal temperature gradients and obtained an anomalously high value. Upon applying this information to the construction of basal temperatures models - 3 -for the Rusty, they concluded that c e r t a i n regions of the base of the g l a c i e r might be at the pressure-melting point. In the spring of 1971, Dr. Clarke discussed these re s u l t s with me and i t was decided at that time that we should pursue a general theory of g l a c i e r surging based upon thermal i n s t a b i l i t i e s that would be induced by the surge process i t s e l f . - 4 -1.2. G l a c i e r Surges A surging g l a c i e r i s defined as one which p e r i o d i c a l l y experiences a sudden and rapid displacement of a large amount of i t s i c e at v e l o c i t i e s that are one or two orders of magnitude greater than the normal flow rates for the g l a c i e r between surges. Those g l a c i e r s which have such high v e l o c i t i e s constantly, such as some large g l a c i e r s i n West Greenland, are not c l a s s i f i e d as surging g l a c i e r s because they do not appear to have the same pe r i o d i c i n s t a b i l i t y . A few examples of g l a c i e r surges w i l l i l l u s t r a t e what amazing na t u r a l phenomena these g l a c i e r s are. The Steele G l a c i e r i n southwestern Yukon T e r r i t o r y began i t s rapid advance sometime during 1965. By August, 1966 a w a l l of i c e t h i r t y meters high had moved twelve kilometers down the lower part of the g l a c i e r , and the press had begun to r e f e r to the Steele as a "galloping g l a c i e r " . However, when a e r i a l photographs were taken i n early 1968, the rapid movement had almost stopped [Stanley, 1969] . T y p i c a l of surging g l a c i e r s , the Steele's surface had been severely crevassed and the thickening of the g l a c i e r i n the lower region was compensated by a slumping of the surface i n the upper zone. The surface lowering f o r the Steele amounted to more than one hundred meters and was responsible for releasing a surge of the Hodgson G l a c i e r , i t s l a r g e s t t r i b u t a r y , i n the following year. In Spitsbergen, the Negri G l a c i e r surged during the winter of 1935-36. The g l a c i e r advanced twelve kilometers along a f i f t e e n kilometer front at an average speed of 35 meters per day [ L i e s t ^ l , 1969]. Perhaps the f a s t e s t v e l o c i t y observed f o r a large surging g l a c i e r was achieved during the 1953 surge of the Kutiah G l a c i e r i n the Karakoram Mountain Range i n the Himalayas. Desio [1954] has reported a sustained speed during this surge of 113 meters per day. Smaller g l a c i e r s can also surge. Tyeen G l a c i e r , at G l a c i e r Bay, Alaska has a t o t a l area of only 6.4 square kilometers, but i t has surged twice i n modern times, i n 1946-1948 and 1964-1966 [ F i e l d , 1969]. Two aspects of surging g l a c i e r s deserve s p e c i a l a t t e n t i o n . One i s t h e i r rather s p e c i a l geographic d i s t r i b u t i o n and the other i s t h e i r p e r i o d i c i t y . Although the s c a r c i t y of reports of surges i n some places probably r e f l e c t s a lack of observations i n those areas, there i s enough evidence to conclude d e f i n i t e l y that the d i s t r i b u t i o n of surging g l a c i e r s i s not random. In North America, a l l but a few of the surging g l a c i e r s are located i n a r e l a t i v e l y small area near the Alaska-Yukon border co n s i s t i n g of the Alaska Range Mountains, eastern Wrangell - 6 -Mountains, eastern Chugach Mountains, I c e f i e l d Ranges, the St. E l i a s Mountains of south-central Alaska near G l a c i e r and Yakutat Bay, and the St. E l i a s Mountains i n southwestern Yukon T e r r i t o r y and northwestern B r i t i s h Columbia [Post, 1969; Meier and Post, 1969; Horvath and F i e l d , 1969]. The p r i n c i p a l exception to th i s trend i s the Otto G l a c i e r on Ellesmere Island.in the Canadian A r c t i c which surged sometime between 1950 and 1959 [Hattersley-Smith, 1964,1969]. No surging g l a c i e r s have yet been i d e n t i f i e d i n the Brooks Range, Kenai Mountains, west and ce n t r a l Wrangell and Chugach Mountains, Rocky Mountains, or S i e r r a , Nevada [Post, 1969]. The main region containing surging g l a c i e r s outside of North America i s the Karakoram Range of the Himalayas. Perhaps 90 per cent of a l l observed surge events occur i n e i t h e r t h i s area or the Alaska-Yukon region [Hewitt, 1969]. Surging g l a c i e r s are also found i n Iceland and Spitsbergen, and i n selected portions.of the Pamir Mountains i n U.S.S.R. and of the Chilean Andes [Paterson, 1969]. There have been no i d e n t i f i c a t i o n s of surging g l a c i e r s i n the other major mountain ranges of the world, inc l u d i n g the Alps. Some doubt about the p e r i o d i c i t y of g l a c i e r surges s t i l l e x i s t s , but.many authors have come to believe that a l l surging g l a c i e r s surge repeatedly at approximately equal i n t e r v a l s . - 7 -At least three observed surges would be required to give any d i r e c t support to the idea that surges are p e r i o d i c ; however, only.a few g l a c i e r s have a c t u a l l y been observed to surge more than once. Tarr and Martin [1914] recorded the dramatic changes i n the Variegated G l a c i e r caused by i t s surge i n 1905-1906. The Variegated was also observed to have surged j u s t p r i o r to 1948 and again i n 1964-1965.[Post, 1969]. Although no observations were made i n the mid-1920's which would v e r i f y whether the g l a c i e r a c t u a l l y surged during that time, this information does seem to suggest a surge period f or the Variegated of approximately twenty years. Thorarinsson [1969] has uncovered h i s t o r i c a l evidence to document four surges of the Bruarjb'kull G l a c i e r i n Iceland which preceded, the .one , that he observed i n , 1963-1964. These e a r l i e r surges occurred i n 1625, around 1730, i n 1810, and i n 1890. Thus the B r d a r j b k u l l surges r e g u l a r l y at i n t e r v a l s of 70-100 years. Indirect,evidence for the p e r i o d i c i t y of g l a c i e r surges-has been cited.by Meier and Post [1969] i n a report of a study of 204 . l i k e l y surging g l a c i e r s which have been i d e n t i f i e d i n North....Ameri,ca [Post, 1969]. A e r i a l photographs of these g l a c i e r s have shown.a s e r i e s , o f o ld medial moraine loops associated with each of the large g l a c i e r s i n the group.. Such loops are formed from the flow of .a-.tributary onto the main g l a c i e r . The loops are c a r r i e d down the main, g l a c i e r when i t surges, and the t r i b u t a r y then begins to form.a new loop behind the old ones. I t i s t h i s s e r i e s of loops, - 8 -often nearly i d e n t i c a l and generally evenly spaced, which suggests that the surges are pe r i o d i c . An example of this s i t u a t i o n i s the case of the Muldrow G l a c i e r , which surged most recently i n 1956-1957, and i t s large t r i b u t a r y , the T r a l e i k a G l a c i e r [Post,I960]. In some instances, as many as ten of these loops can be seen on a sing l e g l a c i e r , thereby giving evidence f o r surges which date centuries i n t o the past. I t should be noted that the absence of these morainal patterns among the smaller g l a c i e r s i n the study can be explained by the ablation of the g l a c i e r between surges [Meier and Post, 1969]. - 9 -1.3 Summary of Previous Theories Theories regarding the causes of g l a c i e r surges have an i n t e r e s t i n g h i s t o r y . Tarr and Martin [1914] proposed that the surges they observed were caused by avalanches of snow and rock r e s u l t i n g from a seri e s of earthquakes that occurred i n the Yakutat Bay region i n 1899. In fa c t most of the surging g l a c i e r s i n the Alaska-Yukon region are located along the same active f a u l t zone [Post, I960]. Nevertheless there i s l i t t l e c o r r e l a t i o n i n space or time between surge events and large earthquakes or avalanches [Post, 1965]. An o r i g i n a l thermal i n s t a b i l i t y theory for g l a c i e r surges was proposed by Robin [1955]. He suggested that, following years of more or le s s steady accumulation, the increased shear stress would cause increased movement along the g l a c i e r bed. This, i n turn, would r a i s e .the temperature of the basal i c e and .allow the g l a c i e r to flow even more rap i d l y u n t i l "conduction and the downward 'advection' of cold i c e again s t a r t to lower the temperature of the ice-rock i n t e r f a c e with consequent (deceleration and) thickening of the g l a c i e r . " At the time, very l i t t l e of the theory of g l a c i e r s l i d i n g had been developed, nor had a s u f f i c i e n t number of observations of surging g l a c i e r s been made to allow Robin to test h i s theory q u a n t i t a t i v e l y . In the l a t e 1960's, he might have f e l t encouraged to reevaluate h i s ideas and explore new d e t a i l s , but, by that time, he had come to believe i n another theory of surging g l a c i e r s based upon stress i n s t a b i l i t i e s and to which thermal properties were mostly i n c i d e n t a l . - 10 -This l a t e r theory [Robin, 1969; Robin and Barnes, 1969] proposed that a surge w i l l occur when (1) the g l a c i e r changes from compressive to extending flow i n the ablation area, or (2) from extending to compressive flow i n the accumu-l a t i o n zone. These two d i f f e r e n t modes of flow are explained t h e o r e t i c a l l y by the choice of the sign of the r a d i c a l i n the r e c t i l i n e a r s o l u t i o n to Glen's flow law for i c e . In order to maintain a constant rate of discharge, Robin observes, a region of the g l a c i e r undergoing extensive flow w i l l be about t h i r t y meters thicker than a section experiencing compressive flow. The cause of the sudden switch can then be explained i n terms of p a r t i c u l a r v a r i a t i o n s i n the bedrock p r o f i l e . Figure 1 i l l u s -trates two such p r o f i l e s f o r the two si t u a t i o n s proposed by Robin. In support of t h e i r theory, Barnes and Robin give two arguments based upon observations. They note that, f o r many surging g l a c i e r s , approximately the same volume of i c e i s involved i n successive, surges of the same g l a c i e r , e.g. Muldrow Gl a c i e r [Post 1960] and Walsh G l a c i e r [Post, 1966]= They also c o r r e c t l y state that the surging of t r i b u t a r y g l a c i e r s caused by the surface lowering which r e s u l t s from the surge advance of the main g l a c i e r f i t s generally.well i n t o stress i n s t a b i l i t y theory of surges. However, they note further that these facts are b a s i c a l l y incon-c l u s i v e i n o f f e r i n g a preference to a surge mechanism based on . , - 1 1 -Ca) Suggested bedrock p r o f i l e i n the accumulation zone of a g l a c i e r which could cause a surge. (b) Suggested, bedrock p r o f i l e i n the ablation zone of a g l a c i e r which could cause a surge. Figure 1 . G l a c i e r surface and bedrock p r o f i l e s f o r Robin's stress i n s t a b i l i t y theory. ( £ i s the s t r a i n rate f o r the l o n g i t u d i n a l stress.) - 12 -stress or thermal i n s t a b i l i t i e s , since a l l of the theories i m p l i c i t l y assume that part of the g l a c i e r w i l l be t r y i n g to advance before the surge begins and w i l l have the freedom to continue to advance as long as the surge l a s t s . In addition, some surging g l a c i e r s seem to be e x h i b i t i n g surge behavior i n s p i t e of a very unfavourable mass balance. The Fisher G l a c i e r i n the Yukon, which has a negative mass balance, i s an example of a g l a c i e r of this type [Post,1971]. There i s a t h i r d type of i n s t a b i l i t y which might be the cause of, g l a c i e r surges. In a theory proposed by Weertman [1962, 1969], the thickness of the water layer between i c e and bedrock i s the c r i t i c a l parameter i n the surge mechanism. When th i s thickness exceeds the s o - c a l l e d c o n t r o l l i n g obstacle s i z e (see Section 2.1), s u f f i c i e n t stress can be released from the bed to reduce s t a t i c f r i c t i o n and allow the g l a c i e r to accelerate. Once the surge has begun the f r i c t i o n a l heating produced by the s l i d i n g movement would melt enough i c e to continue l u b r i c a t i n g the bed. Weertman believed several s p e c i a l c r i t e r i a would have to be met before this mechanism could produce a dramatic surge event: (1) the g l a c i e r should be at l e a s t ten kilometers long, (2) the g l a c i e r should be at the pressure-melting point at the bed, (3) the water beneath the g l a c i e r should flow mainly as a sheet with l i t t l e flow through streams or channels, (4) the g l a c i e r bed should be smoother for large obstacles than smaller ones, and (5) an unusually large shear - 13 -stress should act at the bed, such as one might expect from a big kinematic wave on the glacier. Even though this theory may be too specific to be applied to some surging glaciers, i t does offer a direct solution to the problem of those large glaciers which have very great rates of advance during a surge. Recently Weertman and Robin have been working together on a new modified theory of water-film i n s t a b i l i t i e s with consideration for variations in the basal.shear stress [Weertman, 1971b]. - 14 -CHAPTER 2 THEORY Because each i n d i v i d u a l surge event i s ultimately dependent upon the most i n t r i c a t e considerations of hydrology and g l a c i e r mechanics as w e l l as upon the p a r t i c u l a r t e r r a i n and climate, no theory can explain every aspect of the behavior of surging g l a c i e r s . What has been achieved i n this thesis i s the development of a fundamental thermal i n s t a b i l i t y which can be a property of most surging g l a c i e r s and which o s c i l l a t e s i n d i r e c t coordination with the active and dormant phases. Although we believe that t h i s thermal mechanism i s the p r i n c i p a l n a t u r a l control for surging g l a c i e r s , we admit that f o r some surges, e x p e c i a l l y those which may .occur i n temperate g l a c i e r s , various stress and melt water considerations may be of importance. * The theory can be divided i n t o three parts. The f i r s t part i s concerned with how the g l a c i e r moves during the surge. While some of the increase i n v e l o c i t y i n the e a r l i e r period of the advance i s due to the enhancement of viscous flow which occurs as the base of the g l a c i e r warms toward the pressure-melting point [Robin, 1955], most of the acceleration of a surging g l a c i e r i s assumed to r e s u l t from the sudden i n i t i a t i o n of rapid s l i d i n g . A summary of theory concerning g l a c i e r s l i d i n g i s presented as well - 15 -as a discussion of water flow at the base of a glacier which includes the recent study by Weertman [19 71b]. The second part of the theory introduces the several aspects of heat conduction that have been considered in the proposed thermal instability mechanism for glacier surges. Among these topics are geothermal and f r i c t i o n a l heating, advection of colder ice toward the bed as the glacier thins during a surge, and accumulation at the surface during the dormant phase. Also support is given for the approximation that one can ignore the contribution of internal sources of heat generation, such as those due to strain heating. Finally, in the third section, the complete surge cycle is developed and the effects of different terrains and climates are considered. 2.1 G l a c i e r S l i d i n g A major breakthrough i n the science of glaciology occured when Weertman f i r s t published the fundamental p r i n c i p l e s of h i s theory of g l a c i e r s l i d i n g . Previously there was no theory that had presented a d e t a i l e d explanation.of the p h y s i c a l processes that control the s l i d i n g v e l o c i t y of a g l a c i e r , nor did a theory e x i s t that would enable a quantitative estimate to be made of this v e l o c i t y . Weertman's theory has been c r i t i c i z e d by some, perhaps most severely by L l i b o u t r y , and generally accepted and then r e f i n e d by others, .notably Nye and Kamb. S t i l l i t s s i g n i f i c a n c e i n glaciology i s acknowledged by everyone, and i t s applications extend w e l l beyond i t s considerable importance for surging g l a c i e r s . 2.1.1 S l i d i n g Theory Weertman's theory i s based upon two complementary s l i d i n g mechanisms [Weertman, 1957]. Given that a g l a c i e r bed i s composed of many protuberances of various s i z e s , there e x i s t s a greater hydrostatic pressure on the u p h i l l side of the obstacles than on the downhill side. When the g l a c i e r i s assumed to be everywhere at the pressure-melting point, such a pressure difference w i l l produce a temperature difference across the obstacles giving r i s e to a flow of heat up the g l a c i e r . This u p h i l l flow of heat i s balanced by the downhill flow of g r a v i t a t i o n a l p o t e n t i a l energy as.the g l a c i e r moves forward. That i s , the heat conducted from the downhill side of the obstacles melts i c e on the u p h i l l side and the water that i s formed flows around the obstacles before r e f r e e z i n g . This process of melting and re-freezing i s known as regelation. The second means by which a g l a c i e r may s l i d e along i t s bed r e s u l t s from the enhancement of creep rates due to stress concentrations i n the i c e surrounding the obstacles. Weertman assumed that the i c e would quickly close any cavity which formed on e i t h e r side of the obstacles since the hydrostatic pressure at that point would be several times greater than the shear stress at the bottom of the g l a c i e r . . Evidence f o r both of the s l i d i n g mechanisms proposed by Weertman was found by Kamb and La Chapelle [1964] i n d i r e c t tunnel observations of the Blue G l a c i e r i n Washington*, For a given obstacle s i z e , L, and obstacle spacing, L , the s l i d i n g v e l o c i t i e s due to each mechanism can then be derived [Weertman,, 1964]. In. the case of the regelation process, the v e l o c i t y of s l i d i n g - i s given by the volume of water melted per unit time divided by the area of the obstacle. Following Weeertman's notation, the temperature gradient across the obstacle i s 2 3 equal.to Cx(L') /L , where T i s the basal shear stress and C i s the constant which contains the pressure dependence of the .melting.point of ice (C = AT/AP). This corresponds 2 3 to a heat f l u x through the obstacle of CDx(L') /L where D i s the thermal conductivity of the rock. Thus the t o t a l - 18 -2 heat transferred per unit time i s CDx(L') /L and the rate of melting 2 (and refreezing) i s C D T ( L ' ) /(HpL) where p and H are the density and the la t e n t heat of fusion of i c e . Hence the regelation s l i d i n g v e l o c i t y i s given as S r = (CDx/HpL)(L'/L) 2 = (CDxs 2/HpL) (2 ,1) where the "smoothness", s = L'/L, i s a c h a r a c t e r i s t i c property of the bed. The s l i d i n g v e l o c i t y due to creep-rate enhancement i s found i n a s i m i l a r manner by assuming that the distance over which the stress concentrations act i s of the order of the s i z e of the obstacle i t s e l f . The compressive stress on the u p h i l l side of the obstacle i s given by a = ( T / B ) ( L ' / L ) 2 = T S 2 / 6 ( 2 , 2 ) where g = 1 i f there i s cavity formation and 3 = 2 i f there i s not. The creep rate i s given by Glen's flow law of i c e £ = Bo 1 1 ( 2 . 3 ) where B and n are experimentally determined cons t a n t s s a n c j the s l i d i n g v e l o c i t y i s the creep rate m u l t i p l i e d by the distance over which i t i s applied; S = eL = (Bs^J'/^L. ( 2 . 4 ) c Examining equations (2.1) and ( 2 . 4 ) , one discovers that i s proportional to 1/L which means that the regelation s l i d i n g v e l o c i t y i s l a r g e s t f o r smaller obstacles, and that varies d i r e c t l y as L which means that the s l i d i n g v e l o c i t y due to stress concentrations i s greatest for larger obstacles. Hence one finds that, - .19 -fo r small obstacles, the s l i d i n g motion i s diminished by the l i m i t e d effectiveness of the creep enhancement mechanism, and that, f o r large obstacles, s l i d i n g i s hindered by the poor e f f i c i e n c y of the regelation process. Thus the actual s l i d i n g v e l o c i t y i s p r i n c i p a l l y determined, by those obstacles whose s i z e gives equal values for and S^. This p a r t i c u l a r value f o r L i s c a l l e d the c o n t r o l l i n g obstacle s i z e and i s given by * ;„ . .. X • - . • ( C D S n / B H p ) 1 / 2 t a - n ) / 2 B 1 - n (2.5) The t o t a l s l i d i n g v e l o c i t y , S = S + S , can now r c be expressed using this equation f o r X; , S = 2 ( B C D / H p 3 n ) 1 / 2 x ( n + 1 ) / 2 S n + 1 (2.6) This i s the fundamental r e s u l t , although Weertman [1964] does introduce several other undetermined constants, s i m i l a r to i n order., to generalize his theory. .He has also shown {Weertman, 1967a] that h i s theory can be. applied.to non-temperate g l a c i e r s which have t h e i r base at the melting point i f the temperature gradient at the bed i s not too large. I f the temperature gradient i s very large,., s l i d i n g w i l l not occur. The next step i n Weertman's, theory i s to consider the e f f e c t of a water layer of thickness u> upon the s l i d i n g v e l o c i t y . The v a l i d i t y of the assumption that this water flows as a sheet which can be regarded as having a uniform thickness i s an unresolved - 20 -question, cf. Weertman [1964], Lliboutry [1968a], Kamb [197P], and Weertman [1971b]., However, because we are convinced by Weertman's arguments and because there could be no tractable alternative in terms of .our modeling routines, we accepted this assumption in our theory. This point wi l l be discussed in more detail later in this section. If such a water layer is assumed, i t can have an important effect on the.sliding velocity. For some sufficiently large value of io= to', the contribution of the regelation process wil l become negligibly small and. the sliding velocity wi l l be totally expressed by equation (2.4). At this point;, the velocity is controlled by those obstacles which are just slightly bigger than the height of the. water layer. Hence, for o i > w ' , .... S .= (B/3 n)T ns 2 nu> (2.7) Note the. difference in the power dependence of the parameters T and s between equations (2.6) and (2- 7) . Such a change in the sliding behaviour.of a glacier is believed to contribute sub-stantially to the acceleration of a surge. The choice of the value to' is s t i l l fairly arbitrary at the present stage of the theory. Weertman [1964] suggests 1/10A, while Kamb [1970] would prefer 10A. For the models: presented in this thesis, we have chosen the simple equality, co' = X. . Major improvements and refinements have been made in Weertman's theory as a result of the excellent work by Nye[1969, 1970] and Kamb [1970]. Both authors employ Fourier analytical - 21 -methods to the study of the bedrock topography. Nye further assumes a l i n e a r (Newtonian), viscous approximation and no c a v i -t a t i o n around the obstacles i n order to derive i n t e r e s t i n g analyt-i c a l r e s u l t s for the f r i c t i o n a l drag. Kamb uses a nonlinear theory to derive a b e a u t i f u l l y rigorous expression f o r the s l i d i n g v e l o c i t y which agrees with the power dependence of the shear stress given by Weertman (equation 2.6). His approximations are based principally.upon the assumption that the second i n v a r i a n t of the s t r a i n - r a t e tensor depends only upon the distance from the ice-rock i n t e r f a c e . He also considers i n d i r e c t l y the e f f e c t of a water layer and again ob-tains, r e s u l t s s i m i l a r to Weertman's theory. He c r i t i c i z e s Weertman, however, on the choice of a c o n t r o l l i n g obstacle s i z e , and demonstrates i n his theory that a p a r t i c u l a r obstacle spacing i s the actual c r i t i c a l parameter for g l a c i e r s l i d i n g . For normal.valley g l a c i e r s , . Kamb's re s u l t s i n d i c a t e there should be l i t t l e .cavitation at the base of the g l a c i e r unless melt water with a pressure equivalent to one-half the g l a c i e r thickness has access to the bed. Weertman[1971a] recently has responded to some of the c r i t i c i s m s of his theory o f f e r e d by Nye and Kamb. An alternate theory for g l a c i e r s l i d i n g has been proposed by L l i b o u t r y [1965, 1968a, 1969] and was reviewed by Weertman [1967a]. The differences between, the theories centjer on the issues of s u b g l a c i a l c a v i t a t i o n and the l o c a t i o n of the melt water formed beneath the g l a c i e r . L l i b o u t r y ' s theory a l l o y s c a v i t i e s of l i m i t e d extent to be formed and maintained by the equ i l i b r i u m between the hydrostatic overburden pressure and the water pressure of the s u b g l a c i a l . h y d r a u l i c system. I f co r r e c t , t h i s displacement of water i n t o c a v i t i e s would c e r t a i n l y have an e f f e c t upon the r e l a t i o n s h i p s between g l a c i e r s l i d i n g and theories f o r surging g l a c i e r s . 2.1.2 Water Flow at the Base of a G l a c i e r Weertman [1971b] has now introduced a theory which, at l e a s t q u a l i t a t i v e l y , describes the c h a r a c t e r i s t i c s of water flow, at the.base.of a g l a c i e r or i c e sheet. The bed i s assumed to be impermeable to water, and the discussion turns to the types of drainage channels which may.exist.beneath the g l a c i e r : RtSthlisberger. channels, which are i n c i s e d upwards in t o the ic;e, and Nye channels which are i n c i s e d downwards into the bedrock r In both cases, the destruction of these channels occurs and the water tends to flow approximately as a sheet. For Rb'thlisberger channels the destruction occurs because of melting from the base of the g l a c i e r due to the heat of s l i d i n g , while, for Nye channels, a l l but the largest channels disappear because the rate of g l a c i a l erosion of the bedrock usually exceeds the rate at which the channels.can penetrate into the rock. For surging g l a c i e r s , those drainage channels i n the i c e which remain a f t e r the surge has ended w i l l be destroyed by the refreezing of the water and the subsequent deformation by viscous flow of the g l a c i e r . Hence the melt water generated by the new surge, when i t begins, w i l l be characterized by sheet flow. . . T e c h n i c a l l y , of course, the sheet pf water cannot be of uniform thickness, and t h i s i s demonstrated by Weertman 11971b], Even to a f i r s t approximation, the thickness of the water generally increases.from the, upper accumulation zone to the areas down-stream since the amount of water flowing from one l e v e l i s added to the water produced at.the next lower l e y e l . Also, fo r t y p i c a l v a l l e y g l a c i e r s , the water thickness i s greatest along the centre and decreases toward the edges. However, for a .glacier surge these e f f e c t s involve, only the r e l a t i v e p a r t i c i p a -t i o n of various parts.of,the g l a c i e r in. the surge advance, and such features are beyond the scope of the present theory. For.example, Weertman [1969] has shown from such f i r s t considera-tions that the thickest water layer occurs at a distance upstream from the terminus equal to one-third the t o t a l length of the g l a c i e r , and hence t h i s zone may be the most l i k e l y to f i r s t e x h i b i t the e f f e c t s of a surge. His c a l c u l a t i o n assumes that the g l a c i e r has a p r o f i l e determined by perfect p l a s t i c i t y theory, that the bed i s v i r t u a l l y f l a t , and that the rate at which i c e i s melted from the base of the g l a c i e r i s constant everywhere. In a more d e t a i l e d consideration of t h i s question, Weertman [1971b] - 24 -shows how the v a r i a t i o n of the water l e v e l can be rel a t e d to the obstacles and i r r e g u l a r i t i e s of the bed. Assuming a sheet of uniform thickness, there are several ways to calculate t h i s thickness, u, from the amount of melt water produced per unit area per unit time from the bottom of the g l a c i e r , W. The method used i s based on the work of Lamb [1945] and s l i g h t l y d i f f e r e n t from that suggested by Weertman [1971b]. In Figure 2, the average thickness n which must be melted from one or both g l a c i e r surfaces to produce a given value of W i s shown to be a l i n e a r function of W with a slope that i s inversely proportional to the area of the g l a c i e r . This water then forms a layer between i c e and bedrock as i l l u s -t rated i n Figure 3. The v e l o c i t y d i s t r i b u t i o n i n the layer i s given by f(y) and the g l a c i e r has a slope equal to tan a and a s l i d i n g v e l o c i t y equal to U. The Navier-Stokes equations i n the gap reduce to _9^f _ dp_ 2 d x 9y 0 = ^ (2.9) dx where u i s the v i s c o s i t y of water at the melting point and p represents the s c a l a r pressure. The boundary conditions are f(0) = 0 (2.10) f(ui) = U. (2.11) Figure 2. T o t a l melted', ice thickness n. versus rate of ' water formation W f o r several values of the combined g l a c i e r surface area A. - 26 -Figure 3. Flow of a water f i l m of uniform thickness beneath a g l a c i e r . - 27 -The s o l u t i o n to the system of equations (2.8-2.11) i s given by f ( y ) „ I D . yiazzi dp_ w co 2y dx Lamb [1945] points out that this expression i s approximately v a l i d even i f co varies slowly with the x-coordinate. The pressure de r i v a t i v e i s set equal to a constant which i s given by -pgtana, where g i s the acceleration due to gravity. Also, i t i s required that W = /j!jf(y)dy. (2.13) This corresponds to the assumption that the bed i s everywhere impermeable to water. Hence W = ^ coU + ^ — co3 (2.14) 2 2y where y i s a constant equal to 6y/pgtana. Solving equation (2.14) f o r the water l e v e l y i e l d s co — [YW + ( Y V + Y 3 U 3 / 2 7 ) 1 / 2 ] 1 / 3 + [yW - (y 2W 2 + y 3 U 3 / 2 7 ) 1 / 2 ] 1 / 3 . (2.15) Hence, when W can be known, equations (2,6), (2.7), and (2.15) w i l l completely determine the s l i d i n g v e l o c i t y of the g l a c i e r s i n the surge theory. - 28 -2.2 Heat Conduction In this section, results are obtained regarding the conduction of heat in ordinary and surging glaciers; these shall be used in the development of the thermal instability mechanism for glacier surges. The rigorous analysis of thermal conduction in glaciers presents a particularly d i f f i c u l t problem, The normal thermal regime is often close to the phase transition point between water and ice. Some mass is continually being lost, while new ice is usually being added. The glacier can slide along i t s bed as well as creep under non-linear viscous forces. Inso doing, it. may open crevasses and channels through which air and water can.circulate. Because i t behaves as a ..viscous solid, each internal element of the glacier can be regarded as a minute source of heat. Simplifying assumptions are therefore required. The assumptions we have chosen give the most simplicity for quantitiative models, without destroying the .. essential properties which.relate, to the theory. 2.2.1 Temperatures in Glaciers .. Glaciers have been traditionally classified into three groups according to their temperatures [Ahlmann, 1935]. Polar glaciers; are below the pressure-melting point everywhere including at their base. Sub-polar glaciers are below the pressure-smelting point throughout most of their thickness but may reach the melting point in a thin layer near the bed. Temperate glaciers ...are at the pressure-melting point everywhere, except for seasonal - 29 -f l u c t u a t i o n s near the surface. Sub-polar g l a c i e r s are of p r i n c i p a l i n t e r e s t f or the mechanism of g l a c i e r surging we present i n t h i s t h e s i s . Polar g l a c i e r s are of no i n t e r e s t since they probably do not surge. The question of whether temperate g l a c i e r s can surge has yet to be decided. Our opinion i s that they do not, although we b e l i e v e that i n some instances, surging g l a c i e r s can e x i s t i n regions where ordinary g l a c i e r s would be temperate. This i s because the i n t e r v a l s of strong advection which these g l a c i e r s experience during a surge, has the e f f e c t of b r i n g i n g cold surface i c e to deep l e v e l s within the g l a c i e r . For g l a c i e r s which have a high accumulation rate and frequent surges, t h i s e f f e c t could maintain a sub-polar temperature d i s t r i b u t i o n . The f i r s t assumption that i s made regarding the temperature d i s t r i b u t i o n i s that the annual f l u c t u a t i o n s of 'temperatures near the surface can be ignored. Since the amplitude of a a " . v s , t e m p e r a t u r e wave of frequency co i n c i d e n t at the' 1/2 upper surface.of a g l a c i e r diminishes as exp (-y ( . 013co) ) where co i s measured i n 2ir»yaerB " and y i s measured i n meters, an annual d i f f e r e n c e i n maximum and minimum temperatures of twenty degrees would correspond to changes of only three-tenths of a degree at a depth of t h i r t y meters. Therefore the approx-imation of the surface temperature as a constant (mean annual) value has been employed i n the theory. - 30 -The second assumption is that the percolation of surface melt water through the glacier w i l l not significantly affect the glacier temperatures. Each gram of water at the freezing point.contains enough latent heat to raise the temperature of 160 grams' of ice by one degree. Most of this water w i l l drain completely out of.the glacier. Furthermore, the warming effect of the water which penetrates into crevasses w i l l be approximately balanced by the additional surface area of the glacier that is exposed to cold winter temperatures when crevasses exist. 2 . 2 . 2 Mathematical Theory, for the Conduction of Heat, in Solids The principal mechanism for the transfer of heat in solids is conduction. This is in contrast to fluids where convection and radiation.usually have the dominant role. The most general.equation for thermal, conduction.in an isotropic solid which may be in, motion and have internal sources of heat generation i s given by, [Carslaw and Jaeger, 1 9 5 9 ] . . : . . v x(x,y,z,t ) f + V y ( x , y , z , t ) ^ + V 2 ( x , y , Z , t ) | | 2 2 2 -K + + " 9 " > ~ - ~ A(x,y,z,t) = 0 . ( 2 .16 ) 9x2 9y2 9z 2 P ° or in vector notation, | | + V(x,y,z,t)-vT -K V 2 T - ^ A(x,y,z,t) = 0 (2 .17 ) where T(x,y,z,t) is the scalar temperature f i e l d , V(x,y,z,t) is a vector that represents the particle velocities, and A(x,y,z,t) is the heat generation term. The quantities K,p, and C refer respectively to the thermal d i f f u s i v i t y , density, and s p e c i f i c heat, which are p h y s i c a l properties of the s o l i d and have been assumed constant. They are r e l a t e d by the thermal conductivity,K - Kpc, I f the coordinate axes are chosen such that the x-axis points in, the d i r e c t i o n of g l a c i e r flow, the y-axis i s perpendicular to the bed, and the z-axis i s mutually perpendicular to the other two axes, then to a good approximation the temperature var i a t i o n s can be r e s t r i c t e d to the s i n g l e v e r t i c a l coordinate, y. Equation (2.16) then becomes 2 • .||.+ . v y x , y , z , t ) f^-K -Mf -^ : A(x,y,z,t) = 0 (2.18) By fu r t h e r assuming that .V and A.also,do not vary appreciably i n the l a t e r a l d i r e c t i o n s , the following s i m p l i f i e d equation i s obtained - f + V ( y , t ) f - A - L A ( y , t ) . 0 (2.19) 9y where the subscript of the v e l o c i t y term has now been dropped. In order to obtain a p a r t i c u l a r s o l u t i o n for a physi-c a l l y r e a l i z a b l e problem, equation (2.19) must be solved i n conjunction with several, boundary conditions. For a g l a c i e r , the usual.boundary conditions a r e . ( l ) the i n i t i a l temperature d i s t r i b u t i o n , T(y,0); (2) the temperature at the bed, T ( 0 , t ) , or the gradient,*}) k -g— at the base of the i c e , (0,t); and (3). the temperature at the surface, T ( Y , t ) , where the thickness of the g l a c i e r , Y = Y ( t ) , may be changing with time. These - 32 -boundary conditions, together with.the functions V and A, completely specify the problem of thermal conduction i n a g l a c i e r under the stated assumptions (Figure 4). In the fi g u r e , X(.t) i s . t h e length of th e . g l a c i e r , ,U(t) i s the v e l o c i t y of the g l a c i e r down the v a l l e y , and a i s again the slope of the surface of the g l a c i e r . Because the thickness has been assumed constant along the length of the g l a c i e r . i n agreement with the previous .assumption made, regarding the .velocity V, i t .follows that a w i l l also be the slope of the bedrock. 2.2.3 Accumulation, Ablation, and Advection The changes i n the thickness of the g l a c i e r , Y(t), correspond, to a r a i s i n g or. lowering.of.the g l a c i e r surface. Two d i f f e r e n t mechanisms contribute to the cause of these changes. F i r s t i c e may be added to the g l a c i e r surface by winter..snowfalls (accumulation) or melted away by summer heat ( a b l a t i o n ) . The e f f e c t of these two processes.on the temperature d i s t r i b u t i o n i n a g l a c i e r has been considered by Be n f i e l d [1948,1951]. As one would expect, Benfield's r e s u l t s show that these e f f e c t s require a considerable time to manifest any su b s t a n t i a l f l u c t u a t i o n s i n .temperatures at depths.well below the g l a c i e r surface. Hence, accumulation and ablation can be neglected f o r the short time corresponding to a sing l e surge-type advance, but not for the much longer i n t e r v a l s between surges. The second, and more important, means by which a glacier.may change i t s thickness i s i n response to an extension - 33 -gure 4. The boundary conditions f o r a g l a c i e r i n the s o l u t i o n of the one-dimensional heat conduction equation. (or compression) of the g l a c i e r i n the l o n g i t u d i n a l d i r e c t i o n , as would occur during a surge. This process i s known as advection. H e r e the moving boundary condition at the surface of the g l a c i e r i s of the form T [ Y ( t ) , t ] = T q (2.20) where T q i s a constant equal to the mean annual temperature and where Y(t) i s determined by the v e l o c i t y function V(y,t) i n equation (2.19). Because the function V(y,t) represents the motion of the g l a c i e r perpendicular to the bed, i t must vanish f o r y = 0. Furthermore, i n order for the i c e within the g l a c i e r to remain continuous, V(y,t) must vary l i n e a r l y with distance from the bed. Since V(Y,t) = •— , then v" yfe * f ( 2 - 2 1 ) or . i f sct)4 dt 3 V(y,t) = S(t). (2.22) The function S(t) can be determined from the p r i n c i p l e of conservation of mass for a g l a c i e r of the type shown i n Figure 4, i f i t i s assumed that the f r a c t i o n of mass l o s t by melting during a surge i s n e g l i g i b l e . The mass, of the g l a c i e r i s proportional to X ( t ) , so that i t must be true that . ^ (X(t)-Y(t)) = 0. (2.23) Thus, U(t)-Y(t) + X(t)- S(t) = 0 (2.24) or S(t) = - U(t). (2.25) Hence V(y,t) = - 7 ^ " (2^26) Also, of course, Y ( t ) = Y q + /QS(t')dt' ( 2 . 2 7 ) and X ( t ) = X q + / ^ j ( t ' ) d t ' ( 2 . 2 8 ) where Y = Y ( 0 ) and X = X ( 0 ) . K•-•.->.ce.. i f the v e l o c i t y U(t) i s o o known, as from Weertman's s l i d i n g theory, the other functions X ( t ) , Y ( t ) , S ( t ) , and V(y,t) can be determined. I t should be noted that f o r surging g l a c i e r s , which are thinning over t h e i r active length, S(t) and V(y,t) w i l l be negative. In p r a c t i s e , the boundary condition ( 2 . 2 0 ) i s d i f f i c u l t to apply, both a n a l y t i c a l l y and numerically, except i n very s p e c i a l problems. However, i n the case of surging g l a c i e r s , where the p r i n c i p a l amount of heat production occurs at the bed, an a l t e r n a t i v e boundary condition i s a v a i l a b l e . Consider a point which.is i n i t i a l l y at the center of a. surging g l a c i e r . A temperature d i s t r i b u t i o n which increases monotonically with depth,is assumed. During the surge, conduction of the heat generated at the ice-rock i n t e r f a c e i s not l i k e l y to have much e f f e c t beyond a narrow region near the base of the g l a c i e r , and therefore the center point w i l l have i t s .temperature a l t e r e d only as a r e s u l t of i t s r e l a t i v e movement towards the colder surface i c e , as shown i n Figure 5. I f Y T^(t)=T(™,t) i s defined as t h i s temperature, then T h ( t ) = TC *° -/QV(|^,t')dt',0] (2.29) which can be written, from equation (2.26), as T h ( t ) = T[ +f tQ dt'),0]. (2.30) Since U(t) = , T (t) = xt I°(l + log X(.t)),0]. (2.31) h 2 Defining, H(t) = I°(l + log X(t)) (2.32) one obtains T h ( t ) = T(H(t),0). (2.33) Equation (2.33) i s a boundary.condition .for the surging g l a c i e r which corresponds to equation (2.20) and can be substituted i n t o the theory given t h e . v a l i d i t y of the assumptions under which i t was derived. . . , I t i s obvious that this process of advection will, -tend to make the temperature of a l l the i n t e r i o r points of the g l a c i e r become colder as the g l a c i e r thins during a surge, and not merely f o r the case described. 2.2 0.4 Heat Sources at the Gl a c i e r Bed There are two sources of heat at the .glacier bed, the geothermal heat flow from the earth's i n t e r i o r and the f r i c t i o n a l heat produced by the movement of the g l a c i e r over the bedrock. - 37 -SCt) Figure 5. The replacement of the moving boundary temperature due to advection during a surge by a v a r i a b l e temperature at a f i x e d point w i t h i n the g l a c i e r . - 38 -Although the geothermal heat flow for a given area w i l l remain r e l a t i v e l y constant for a long period of time, c r i t i c a l changes i n the amount of heat a c t u a l l y flowing i n t o the g l a c i e r can occur i n response to,small f l u c t u a t i o n s of the temperatures i n the upper few meters of bedrock. Therefore i t i s not s u f f i c i e n t f o r the theory to assume a constant value for this heat flow, but instead a continuous c a l c u l a t i o n of the bed temperatures must be done, so ..the, precise value of the bedrock temperature, gradient can be found. The geothermal heat flow i s given by % " -Vb ( 2- 3 4 ) where and cj^ are the thermal .conductivity and temperature gradient of the bedrock and ^ .<0.for the case where temperatures increase with depth. The temperature gradient ^ i s found from the s o l u t i o n of the one-dimensional heat equation f o r the temperatures i n the bedrock, where the problem i s s i m p l i f i e d because, the functions U(y,t) and V(y,,t) of equation (2.19) vanish; T h a t - i s , for ...the, bedrock, 2 | f - K ^ = 0 (2.35) 3y and * b ( t ) * If ly-0 • ( 2 ' 3 6 ) Again, boundary conditions are required to obtain a s o l u t i o n to equation (2.35). The i n i t i a l temperature d i s t r i b u t i o n T(y,0) i s assumed to be known. At some depth w e l l below the surface, the temperatures and gradient of the bedrock w i l l not be affected by the changes i n g l a c i a l a c t i v i t y that, are being considered., One i s then j u s t i f i e d - i n assuming a second boundary condition of the form T(Y f,t) = T f (2.37) where Y^ and are both constants. A reasonable value f o r Y^ would be on the order of one or two hundred meters. The boundary condition at the ice-rock contact (y=0) may be of two types,, depending on.whether there e x i s t s a f i l m .of water between.the two surfaces. If no f i l m i s present, then the heat flow across the i n t e r f a c e must be constant. That i s . • K''d> = K d> (2.38) g yg b Tb where K and are the thermal conductivity and temperature gradient at .the base of the g l a c i e r , or 4 g = 0 b (2.39) where the .approximation = has now been made. Equation (.2.39) i s then, the appropriate boundary condition which must be s a t i s f i e d by-both the g l a c i e r and. bedrock temperatures. I f a water f i l m does e x i s t between i c e and rock, then the heat flows i n the two media need not be equal and equation (2.38) would not n e c e s s a r i l y apply. In f a c t , i t can be seen that f.§ simply corresponds to a change i n the amount of water present i n the f i l m , with a d d i t i o n a l i c e being melted when cj> > - 40 -and some water refreezing to the base, of.the g l a c i e r when (j) <(j>^..- The appropriate boundary condition f o r this case i s c l e a r l y that the temperature of both g l a c i e r and bedrock along th i s surface be fi x e d at the pressure-melting point. Thus T(0,t) = T (2.40) m where i s very close .to 0°C, i s the f i n a l boundary condition that i s required.for a complete s o l u t i o n of the g l a c i e r and bedrock temperatures. This, d u a l i t y of the boundary conditions at the ice-rock i n t e r f a c e , as represented by equations (2.39) and.(2.40), forms an important part of the proposed thermal i n s t a b i l i t y mechanism. Consider.now f r i c t i o n at the.base of the g l a c i e r . F r a c t i o n a l heating w i l l .only, be s i g n i f i c a n t when the g l a c i e r i s able to s l i d e along i t s bed. The.work done per unit area upon the g l a c i e r i n moving a.distance A * i s .given by. TA x, where T i s ..the .basal shear stress.. Some of th i s work w i l l go in t o s l i g h t l y r a i s i n g the temperature, of the melt water r e l a t i v e to . the i c e ; .hence .only ascertain percentage.determined by an e f f i c i e n c y . f a c t o r which s h a l l be defined as v , w i l l be available to melt i c e . From this discussion, the amount of heat produced ...per unit., time and per. unit area by f r i c t i o n can now be written as Q - f . M (2.41) . where J i s the mechanical equivalent of heat (J = 4.18 jo u l e s / c a l o r i e ) , or, f i n a l l y , Q(t) = ^ u(t). (2.42) - 4.1 -The e f f i c i e n c y f a c t o r v can have a d d i t i o n a l p h y s i c a l i n t e r p r e t a t i o n besides the heating of water. The surge may come to a mechanical,, rather than a thermal, stop by advancing onto an area of less slope and greater roughness or by s l i d i n g into i t s own terminal moraine from previous surges.. These conditions can be approximated, in.our.model.by the choice of the fa c t o r v„ Since i t has been assumed that the g l a c i e r i s s l i d i n g and i s therefore at the .pressure-melting point at the bed, the. heat which i s produced.will only melt a d d i t i o n a l i c e from the base .of ..the g l a c i e r . The ..additional, water w i l l further l u b r i c a t e the,bed, allowing the glacier., to. s l i d e faster.. This, i n turn, produces more heat, according to equation (2.42) which w i l l melt more i c e , and therefore i t appears that the high v e l o c i t i e s of ..surging g l a c i e r s can. be. obtained i n this manner. ; -.. The, f i n a l step .is, to, write, the expression, for the t o t a l amount.of melt water .produced per.,unit time and per unit area based, upon.,the, considerations that, have .been presented i n this s e c t i o n : • W(t) = (Q(t) - K(b(t) -• 4> (t>))/L . . . (2.43) In .the equation, L i s . t h e l a t e n t heat of fusion of i c e , and K i s a mean value of the thermal c o n d u c t i v i t i e s of i c e and bedrocks - 42 -2.2.5 Internal Sources of Heat Generation The previous discussion has obviously neglected the heat, generation term, A ( y , t ) , which appears i n equation (2.19). This term corresponds to i n t e r n a l sources of heat production which, i n a g l a c i e r , r e s u l t from the d i s s i p a t i o n of viscous forces i n the.ice. It w i l l be shown that this s t r a i n heating w i l l have l i t t l e e f f e c t upon the thermal regime of a surging g l a c i e r . The equation for s t r a i n heating i n a g l a c i e r i s given by Paterson. [1969, P-179 ], and reduces, f o r the one-dimensional case, to • Qg(y) = e ( y W y ) / J . (2,44) The s t r a i n rate, e(y) can be expressed i n terms of the s t r e s s , x(y), by ..Glen's flow, law (equation 2.3), and making this sub-s t i t u t i o n gives Q s(y) = B ( c ( y ) ) N + 1 / J . (2.45) The parameter B depends on temperature i n the following way B = B e." q / R T (2.46) o where R i s the gas constant (R = 8.31 joules/mole°K) and q i s the. s o - c a l l e d .activation energy f o r creep (q = 58,520 joules/mole)„ Taking.the temperature to be 0°C (=273°K), the values for B vary considerably depending .on the choice of n and B q . T y p i c a l values are.those of Weertman [1964] who uses n = 3 and B = 0.017 -3 -1 bar y r . , Substituting- the numerical values for the constants, equation (2.45) becomes Q s(y) = 0.004 ( a ( y ) ) 4 . (2.47) - 43 -For the one-dimensional case, i t is also true that o(y) = p g:Y-y)sina. (2.48) Arbitrarily choosing a = 3.0° and evaluating the constants, one has a(y) = .00469(Y-y) (2.49) and Q s(y) = ( l o 9 2 x l 0 " 1 2 ) ( Y - y ) 4 (2.50) where Y and y. are measured in meters and Qg is given in calories/ 3 meter . Q s(y).is now integrated over the deepest 10 meters of the glacier and. compared, with the,heat produced by f r i c t i o n . : ..Qs„ -yj° ( 1 . 9 2 x l 0 " 1 2 ) ( Y - y ) 4 d y (2.51) = 3 . 8 4 x l 0 ~ 1 3 ( ( Y - 1 0 ) 5 - Y 5 ) . (2.52) Recall that the f r i c t i o n a l heating was given by • Q = j xU (2.53) = (vpgsina/J)YU (2.54) Evaluating the.constant (letting v= 0.10),,one obtains Q = (4.7xl0~ 4)YU (2,55) Table I compares Q and Qg for various values of Y and U of surging glaciers. It can be seen that the f r i c t i o n a l heating exceeds the.integrated strain heating by two or three orders of magnitude in. each instance. Also, according to equation ( 2 . 5 0 ) , one.would expect the strain heating to diminish rapidly as one moved, toward shallower depths in the glacier. Y(m) U(m/yr) 100. 100 500 1000 200 100 500 1000 300 100 500 1000 400 100 500 1000 - 44 -TABLE I — 2 Q s(cal/cm -yr) 0,012 0,169 0.829 2.580 2 Q(eal/cTO -yz) 4=70 23.50 47o00 9.40 47.00 9^.00 14.10 70.50 141.10 18.80 94.00 188.00 COMPARISON OF THE STRAIN HEATING AND FRICTIONAL HEATING FOR SURGING GLACIERS - 45 -2.3 The Surge Cycle We now unify the previous discussions of g l a c i e r s l i d i n g , thermal conduction, advection,.and basal melting to form a model of the thermal i n s t a b i l i t y mechanism,, The g l a c i e r i n Figure 6 could be representative of a surging g l a c i e r i n the proposed theory. I t i s a sub-polar g l a c i e r with-.mean, annual surface tempetature, T <0°C, and with i t s bed temperature j u s t s l i g h t l y o below the pressure-melting point. I t receives geothermal heat according to the geothermal- gradient, G<0, which i s well-established at .depth. The equ i l i b r i u m p r o f i l e , determined by. thermal conduction theory i s i l l u s t r a t e d .by the broken l i n e . However, th i s g l a c i e r i s c l e a r l y too thick to be. i n equilibrium. Perhaps i t has recently undergone a.heavy period of accumulation, or perhaps i t has experienced an increase in..the mean, annual surf ace temperature caused by a gradual warming of the l o c a l climate. In any case, it;, i s assumed that, by some means, the g l a c i e r i n h e r i t e d this present temperature d i s t r i b u t i o n . To achieve the new l i n e a r e q u i l i b r i u m curve (shown by the dashed-and-dotted l i n e i n Figure . 6 ) , the g l a c i e r w i l l begin to warm under, the influence of the geothermal. heat... Before the basal temperature can a t t a i n i t s new equilibrium value, i t w i l l reach the pressure-melting point, beyond which . it.can. warm no further. When the melting point i s reached,' the g l a c i e r w i l l begin to s l i d e and a thi n layer of water w i l l form.between the base of the g l a c i e r and the bedrock. I f the bed.is s u f f i c i e n t l y smooth and the basal shear stress s u f f i c i e n t l y large, then the theory which was developed e a r l i e r I l l u s t r a t i o n of the pre-surge temperature p r o f i l e . predicts that this s l i d i n g motion could melt enough a d d i t i o n a l water to allow very large s l i d i n g v e l o c i t i e s to be reached. Assuming these.conditions do e x i s t , and that a surge advance does occur, at l e a s t two processes w i l l eventually cause the surge to decelerate and f i n a l l y stop. These are the creation of a s u b g l a c i a l drainage system.which accelerates the flow of water from beneath the g l a c i e r , and the advection of..colder i c e towards the bed. which .decreases the rate of .melting. This l a s t f a c t follows d i r e c t l y from equation (2;43). As colder i c e moves toward the bed, the g l a c i e r temperature gradient (t) , w i l l become larger negative. When the equality , * gft) = 4>bCt) " ^ (2.56) i s . s a t i s f i e d , W(t) = 0 and the production of melt.water ceases. However, the .glacier continues to s l i d e and then one would expect the g l a c i e r temperature gradient, to continue increasing in. magnitude due to advection. In this case, W(t)<0 corresponding to melt water refreezing to the base of the g l a c i e r . This refreezing can continue u n t i l nearly a l l of the remaining water disappears, and the g l a c i e r f i n a l l y refreezes to the bed, ending the surge. The advection which has occurred during the surge s i g n i f i c a n t l y . c o o l s the i n t e r i o r of the g l a c i e r . A d d i t i o n a l cooling of the g l a c i e r i s accomplished by the addition of new surface i c e at the mean annual temperature. The s i t u a t i o n i s i l l u s t r a t e d i n Figure 7. - 48 -Figure 7. I l l u s t r a t i o n of the post-surge temperature p r o f i l e . - 49 -If the surging glacier i s now allowed to accumulate at a moderate to heavy rate during i t s dormant phase, i t w i l l again reach the pre-surge condition of Figure 6, and another surge w i l l be inevitable when the basal ice temperature reaches the pressure-melting point„ The time required before the second surge may occur is determined from the solution to the one-dimensional heat equations in the glacier and bedrock, and w i l l partially depend upon the size of the last surge advance. Since long, deep glaciers on f l a t terrain tend to have large advances, their intervals between surges w i l l be greater than the intervals between surges for shallow, steep glaciers. Regarding climate, we would expect glaciers in regions of relatively warm temperatures .and high accummulation rates to surge more,frequently than glaciers in cold climates with low accumulation rates.. .It. i s required that the climate be sufficiently warm so that the base of the glacier w i l l occasionally reach the pressure-melting point, and sufficiently cold so that thermal conduction and advection processes w i l l have a significant effect. It is believed that this last statement, more than any other.explains the special geographic.distribution of surging glaciers„ In addition to offering explanations for the distribution and apparent periodicity of glacier surges, the proposed thermal - 50 -i n s t a b i l i t y mechanism has another compelling feature. That i s , since the o s c i l l a t i o n between the active and dormant phases i s controlled by thermal conduction at the ice-rock interface, i t i s not necessary to recharge any reservoir by accumulation between surges. Thus this theory allows several surges of glaciers with zero or negative mass balance, although i t i s true that eventually the glacier thins i t s e l f so much that the mechanism becomes inoperable. CHAPTER 3 RESULTS In order to test our proposed theory, we created a computer modeling procedure for surging g l a c i e r s . Although not many quantitative r e s u l t s were obtained, the modeling approach was responsible f o r discovering several early errors i n the development of the theory.. As more data become ava i l a b l e from surging g l a c i e r observations.and experiments, computer programs s i m i l a r to those, created f or this thesis should be very useful . 3.1... The Modeling Procedure We are aware of no other previous attempts to model the thermal regime of a surging g l a c i e r . Models, of g l a c i e r surges .based upon flow .considerations.were obtained by Campbell and Rasmussen [1969] by assuming a l i n e a r r e l a t i o n between volume transport and basal shear s t r e s s . The .following parameters constitute input f o r the model: the length of.the g l a c i e r , the depth, the surface slope, the mean annual temperature, the bedrock smoothness, the geothermal gradient, the e f f i c i e n c y : o f f r i c t i o n a l heating, and w the f r a c t i o n of the water produced per unit time which leaves the g l a c i e r . An i n i t i a l temperature d i s t r i b u t i o n f o r the g l a c i e r and bedrock i s also part of the input. For the model g l a c i e r s the temperature at the bed was chosen.to be at the pressure-melting point, and the temperatures decreased with the distance from the bed. Curvatures for the g l a c i e r temperature p r o f i l e were chosen to resemble Figure 6. For the Rusty G l a c i e r model, the i n i t i a l temperatures were taken from the r e s u l t s of the deep i c e temperature measurements of Classen and Clarke [1971]. In each case, the i n i t i a l bedrock temperatures were chosen to be roughly l i n e a r with dep th. The program f i r s t determines whether the base of the g l a c i e r i s at the pressure-melting point. If not, the model g l a c i e r i s not allowed to s l i d e and c o n t r o l i s transferred to the part of the program.which solves the heat conduction equations. When the bed.reaches the melting point, control w i l l return.to the s l i d i n g routines. These routines calculate the heat generation and water production terms f o r that time increment; then they solve the.conduction equations with..the d i f f e r e n t boundary conditions explained i n section 2.2.4. The s l i d i n g v e l o c i t y i s then calculated and the model g l a c i e r advances. The procedure i s then repeated u n t i l the s l i d i n g (or surging) stops and the g l a c i e r refreezes to the bed. At that point, .the dormant phase begins again. The consideration of numerical techniques i n the o r i g i n a l program forced us to restore the model to .the same i n i t i a l thickness a f t e r each advance. This corresponded to assuming i r r e g u l a r accumulation rates, and eliminated accumulation as a possible parameter i n the model. As t h i s thesis was being finished,, we were w r i t i n g a new program which did not contain this l i m i t a t i o n . . . . - 5 3 -3o2 Numerical Procedures The solution to the system of equations [ ( 2 . 1 9 ) , ( 2 . 2 2 ) , ( 2 . 3 3 ) , ( 2 , 3 5 ) , ( 2 . 3 7 ) , ( 2 . 3 9 ) , and. (2 .40 ) ] which together with the i n i t i a l temperatures, completely specifies the one-dimensional thermal.conduction problem for surging glaciers does not exist analytically and'.must be obtained numerically. The general approach used in the solution of linear partial differential equations such as: (2.19). and (2 .35) i s the method of f i n i t e differences, and the particular technique used.in our modeling program is a variation of the classical Crank-Nicholson methods [Forsythe and Wasow, 1965] . The f i r s t step to be taken is to non-dimensionalize the variables. Several schemes are possible;.the one we employ is as follows: e = | (3 .1) o £ = (3 .2) o T = —j (3 .3) Y o Si = ^ o V (3 .4 ) d With these substitutions equation (2„19) becomes 12 + n 39 - ^ = 0 (3 5) where the internal heat generation term has been dropped. Similarly equation (2 .35) becomes The boundary conditions also transform: eqn.(2.33)-> G^x) = 0(H(x),O) (3.7) eqn. (2.39)-* g b (3.8) and 35 eqn.(2.40)+ 9(0,x) - 0 (3.9) T(y,0) -> 9 ( 5,0). (3.10) Equation (3.5) and (3.6) can now be expressed i n terms of f i n i t e differences [McCracken and Dora,. 1964; Carnaha, Luther, and Wilkes, 1969J. Assuming that 0 ( 5 ,x) has a s u f f i c i e n t number of p a r t i a l derivatives the value of 0 ( 5 ,x) can be expressed as a Taylor expansion 0(C+A5,x+Ax) =0(5,x)+(A?|-+Ax - ~ ) 0 ( 5,x) + 7^1)7 ( A 5 i 5 + ^ h n ~ l Q a > T ) + Y ( 3 - U ) where the remainder, R^, i s given by R n = ~ 7 ( A ? | Y +Ax |^) n 9(5+pA5,x+pAx) (3.12) with 0 1 5 ) 9' 2 „2_ 0 . .. .-20. „+0 . i _ | - _ l z i J X>J 1 + 1 »J- + 0 ( A C 2 ) (3.16) a r ( A C ) where the symbol 0(x) means there e x i s t s a p o s i t i v e constant M such that |B ! - M « X as x-K). We define 0 . , . - 0 . , . V ° i , j ; 2 A ? (3.17) and 0. - 0. . , <5 ( 0 . .) = 1 , 3 ~ l (3.18) T ^ 2Ax as the c e n t r a l difference operators. The values of the p a r t i a l d erivatives can then be approximated as [Carnahan, Luther, and Wilkes, 1969] a e = 6 (a. .) (3.19) 30 = 6 ( 0 . .) (3.20) and 3x 1 X » J 2 3 9 ^ ^ ( Q . .) +hr2(Q. . , ) . (3.21) 3^2 2 5 i , j ' 2 5 i , j + i ' The l a s t three equations are the basis of the Crank-Nicholson method of s o l u t i o n of para b o l i c p a r t i a l d i f f e r e n t i a l equations. [Crank and Nicholson, 1947]. - 56 -Using the central difference notation, we rewrite equations (3.5) and (3.6) as 6 (0. ... fi5£(0. .) - 6^ 2(6. .) + hc2 (9. ...) = 0 (3.22) and 6 T ( Q i j j ) - fs ? 2(ei,3) - | « 5 2 ( 9 l f = 0 (3.23) These two equations are represented, at each point in time, as tri-diagonal systems of simultaneous linear equations at the M-l grid-points ( i = 1, 2, 3, . .. , M, and 0. „ and 0„, . both being known). These are then solved using the standard Gaussian elimination technique [Carnahan, Luther, and Wilkes, 1969]. - 57 ~ 3.3 Results for Surging Glaciers 3.3.1 Model Glaciers Tables I I , I I I , and IV show the r e s u l t s f or three of the models that were calculated during the course of the work for this t h e s i s . Model No. .103, which had parameters corresponding to the Muldrow G l a c i e r i n the Alaska Range yie l d e d sporadic surges with an average length of 5.15 years and an average dormant period of 55.15 years. Meier and Post [1969] give 2.0 and 50±10 years. Model No. .104, which corresponded to no p a r t i c u l a r surging g l a c i e r , also had sporadic surge behavior with averages of 2.22 and 41.8 years. Model No.. 105 was the most successful from the point of view of p e r i o d i c i t y . The model, which was created to .approximate.possible conditions of the Tikke G l a c i e r i n the Alsek River area, of B r i t i s h Columbia, had nineteen nearly i d e n t i c a l surges during f i v e centuries which l a s t e d 1.75 years each and were separated by. dormant i n t e r v a l s of 24 years. The observed, values for Tikke G l a c i e r are 3.0 and 20±10 years [Meier and Post, 1969]. Problems were encountered i n the models for the Walsh Gl a c i e r (No. 101) and the .Steele G l a c i e r (No. 113). Some surges did occur, but the bed temperatures were quick to return to the melting point. This d i f f i c u l t y i s believed to e x i s t because of the uneven accumulation rates that constitute a systematic err o r i n the.present modeling procedure. - 53 -TABLE II Input Parame t e r s : X. 46 km. o Y 290 in. o a 2.2° T -8.0 °C o G = -0.03 °C/m. s = 20.0 w = 0.0 V = 0.001 Numerical Parameters: Ax = 0.05 A? = 0.04 Output.Parameters (Averages): .Length of Surge =5.15 yrs. Advance = 720 m Ve l o c i t y = 141.6 m/yr .Advection Rate = -0.89 m/yr Dormant I n t e r v a l s 55.15 yrs. RESULTS FOR MODEL GLACIER NO. 103 - 59 ~ TABLE III Input Parameters: X = 16 km. o Y = 200 m. o a = 4.5° T = -3.75 °C o G = -0.02 °C/m. s = 20.0 w = 0.50 v = 0.001 Numerical Parameters: AT = 0.05 A5 = 0.04 Output Parameters (Averages): Length of Surge = 2.22 yrs. Advance = 602 m. Ve l o c i t y = 261.9 m/yr. Advection Rate = -3.04m/yr. Dormant In t e r v a l = 41.8 yr. RESULTS FOR MODEL GLACIER NO. 104 - 60 -TABLE IV Input Parameters: X o — 20 km. Y o = 250 m. a = 3.6° T o = -6.7°C G = -0.028°C/m. s 15.0 w = 0.0 V — 0.001 Numerical Parameters: Ax = 0.05 A? = 0.10 Output Parameters (Averages): Length of Surge = 1.75 yrs. Advance = 155 m. Velocity = 87.9 m/yr. Advection Rate = -1.08 m/yr. Dormant Interval = 24.0 yr. RESULTS FOR MODEL GLACIER NO. 105 S i m i l a r l y , the smaller g l a c i e r s that were modeled, the Variegated (No. 106) and the Tyeen (No. 107), took excessively long to recover from each surge, because i t was not possible to incorporate t h e i r high accumulation rates into the program. The g l a c i a l parameters were taken from Meier and Post [1969], with the exception of Steele G l a c i e r [Stanley, 1969]. The geothermal gradients f or northwestern North America were obtained from Clark [1965J. 3.3.2 Rusty G l a c i e r The Rusty G l a c i e r was also modeled (No. 109). The i n i t i a l temperatures measured by Classen and Clarke [19 71] are shown i n Figure 8. The model was found to remain frozen to the bed f o r at l e a s t 150 years. However, when the Rusty model was allowed to increase i t s thickness to 120 meters, i t surged within 30 years. Hence, i f the depth.(or temperature) measurements have appreciable e r r o r s , our p r e d i c t i o n that the Rusty w i l l not surge again soon could be i n c o r r e c t . - 62 -Depth (m) 10 20 30 4 -40 50 4 -60 4-70 80 90 100 Temperature (°C) -7 -6 -5 -4 -3 H 1 1 1 -f-Hole 2B -2 -1 0 H 1— Hole 2C Gravity depth: 88 m Figure 8. Temperature profile for the Rusty Glacier, Yukon. CHAPTER 4 SUMMARY AND DISCUSSION This thesis has presented the theoretical foundation for a mechanism of glacier surges which.is characterized by thermal conduction i n s t a b i l i t i e s at the interface between glacier and bedrock. The theory encompasses most of.what is presently known about surging glaciers and .suggests many other aspects of their.behavior. . A.limited program.of computer modeling was completed in-order to test the proposed mechanism. The results, which are given in Tables .II, III, and IV, were encouraging, i f not conclusive. Several .problems s t i l l remain in the modeling routines, and their rapid solution.can only serve to clarify the situation regarding the theory. With one exception, the results of the modeling study should be regarded as generally a r t i f i c i a l , since parameters and i n i t i a l temperature profiles for these models were chosen arbitrarily. In the case of the Rusty Glacier, Yukon, an effort was made to predict the surge behavior of this glacier from i n i t i a l temperatures given by the thermal d r i l l i n g experiments of Classen and Clarke. Unfortunately, the result of this attempt was rather uncertain, although there were indica-tions that the Rusty would not surge again in the very near future. As a f i n a l point, we offer the suggestion that the present geographic distribution of surging glaciers might be correlated to the long-term climatic changes over the earth during the Present and Pleistocene periods. I t was mentioned i n the discussion of Figure 6 that the pre-surge temperature p r o f i l e might be obtained by a gradual warming of the l o c a l climate. This l o c a l v a r i a t i o n might w e l l be part of the gradual warming on a global scale that accompanies the end of an i c e age. As the warming occurs, some polar g l a c i e r s w i l l become sub-polar and the thermal regime of the g l a c i e r and bedrock system might be s u f f i c i e n t l y disturbed to allow thermal conduction i n s t a b i l i t i e s to slowly develop. I f this hypothesis were true, we would expect that, i n the past, the g l a c i e r s of the Alps and the Western United States, f o r example, could have been surging g l a c i e r s . We would further expect that, i f the present warming trend continues i n t o the future, the present d i s t r i b u t i o n of surging g l a c i e r s w i l l no longer surge, and w i l l be replaced by another group of g l a c i e r s which would be located nearer the Pole. In this connection, we note that Wilson [1969] has proposed a theory for the beginning of new i c e ages that i s a c t u a l l y based upon the surge of the Antarc-t i c Ice Sheet i t s e l f . Of course, these suggestions are tremendously o v e r - s i m p l i f i e d and have not been supported by e i t h e r theory or observation. However, such p o s s i b i l i t i e s should help to make further research of surging g l a c i e r s even more e x c i t i n g i n the future. - 65 -BIBLIOGRAPHY - 66 -Ahlmann, H.W. (1935) Geographical Journal, 86, 97. B e n f i e l d , A.E. (1948) A problem of the temperature d i s t r i b u t i o n in a moving medium. Quarterly of Applied Mathematics, 6_, 439. Benfi e l d , A.E. (1951) The temperature i n an accumulating snow-. f i e l d . Monthly Notices of the Royal Astronomical Society, Geophysical Supplement, _6, no. 3, 139. Campbell, W.J. and Rasmussen, L.A. (1969) Three-dimensional surges and recoveries i n a numerical g l a c i e r model. Canadian Journal of Earth Sciences, 6_, 9 79. Carnahan, B., Luther, H.A., and Wilkes, J.O. (1969) Applied Numer-i c a l .Methods. John Wiley and Sons, New York. Carslaw, H.S. and Jaeger, J.C. (1959) Conduction of Heat i n So l i d s . Oxford.at the Clarendon Press. Clark, S.R. (1965) Handbook of Physical Constants. Geological . . . Society of America, New York. Classen, D.F..(1970) Thermal d r i l l i n g and deep ice-temperature measurements on the "Fox" G l a c i e r , Yukon. M.Sc. Thesis, University of B r i t i s h Columbia. Classen, D.F. and Clarke, G.K.C. (1971) Basal hot spot on a.surge type g l a c i e r . Nature, 229, 481. Crank, J . and Nicholson, P. (1947) A p r a c t i c a l method for numerical evaluation of solutions of p a r t i a l d i f f e r e n t i a l equations of the heat conduction type. Proceedings of the Cambridge P h i l o -s o p h i c a l Society, 4_3, 50. Crossley, D.J. (1969) Gravity and temperature measurements on the "Fox" G l a c i e r , Yukon. M.Sc. Thesis, University of B r i t i s h Columbia. Desio, A. (1954) An exceptional advance i n the Karakoram-Ladakh region. Journal of Glaciology, 2^ , no. 16, 383. F i e l d , W..0. (1969) Current observations on three surges i n G l a c i e r Bay, Alaska, 1965-1968. Canadian Journal of Earth Sciences, 6, 831. Forsythe, G.E. and Wasow, W.R. (1960) F i n i t e - D i f f e r e n c e Methods for P a r t i a l D i f f e r e n t i a l Equations, John Wiley and Sons, New York. - 67 -Hattersley-Smith, G. (1964) Rapid advance of g l a c i e r i n northern Ellesmere Island. Nature, 201, 176. Hattersley-Smith, G. (1969) Recent observations on the surging Otto g l a c i e r , Ellesmere Island. Canadian Journal of Earth Sciences, _6, 883. Hewitt, K. (1969) G l a c i e r surges i n the Karakoram Himalaya . (Central A s i a ) . Canadian Journal of Earth Sciences, 1009. Horvath, E.V. and F i e l d , W.O. (1969) References to g l a c i e r surges i n North America. Canadian Journal of Earth Sciences, 6_, 845. Kamb, B. (1970) S l i d i n g motion of g l a c i e r s : theory and observation. Reviews of Geophysics and Space Physics, _8, no. 4, 6 73. Kamb, B. and LaChapelle, E. (1964) Direct observation of the mech-. anism.of g l a c i e r s l i d i n g over bedrock. Journal of Glaciology, ,_5, no. 38, 159. Lamb,. H. ;(1945) Hydrodynamics (6th ed.). Dover P u b l i c a t i o n s , New York. L i e s t ^ l , 0. . (1969). G l a c i e r surges i n West Spitsbergen. Canadian Journal of Earth Sciences, 6, 895. Ll i b o u t r y , L. (1965) Traite' de. g l a c i o l o g i e . P a r i s , Masson et Cie (2 vols.,) . Ll i b o u t r y , L. (1968a) General theory of s u b g l a c i a l c a v i t a t i o n and s l i d i n g of temperate g l a c i e r s . Journal of Glaciology, J7, no. 49, 21. Lliboutry., L. (1968b) Steady-state temperatures at the bottom of ice sheets and computation of the bottom i c e flow law from the surface p r o f i l e . Journal of Glaciology, _7, n o ' 51, 363. Ll i b o u t r y , L. (1969) The dynamics of temperate g l a c i e r s from the det a i l e d viewpoint. Journal of Glaciology, J3, no. 53, 185. McCracken, D.D. and Dorn, W.S. (1964) Numerical Methods and For-tran Programming. John Wiley and Sons, New York. Meier, M. (1969) Seminar on the causes and mechanics of g l a c i e r surges, St. H i l a i r e , Canada, September 10-11, 1968: a sum-mary. Canadian Journal of Earth Sciences, 6^, 987. Meier, M. and Post, A. (1969) What are g l a c i e r surges? Canadian Journal of Earth Sciences, 6_> 807. - 68 -Nye, J.F. (1969) A calculation, of the s l i d i n g of i c e over a wavy surface using a Newtonian viscous approximation. Proceedings of the Royal Society of London, A, 311, 445. Nye, J.F. (19 70) G l a c i e r s l i d i n g without c a v i t a t i o n i n a l i n e a r viscous approximation. Proceedings of the Royal Society of London, A, 315, 381. Paterson, W.S.B. (1969) Physics of G l a c i e r s . Pergammon Press, London. Post,. A. (1960) The exceptional advan ces of the Muldrow, Block Rapids, and Sustina G l a c i e r s . Journal of Geophysical Research, 65, no. 11, 3703. Post, A. (1965) Alaskan g l a c i e r s : recent observations i n respect .to the earthquake-advance theory., Science, 148, 366. Post, A. (1966) The recent surge of Walsh G l a c i e r , Yukon and Alaska. Journal of Glaciology, 6_, no. 45, 375. Post, A. (1969) D i s t r i b u t i o n of surging g l a c i e r s i n western North America. Journal of Glaciology, j5, no. 53, 229. Robin, G. de Q. (1955) Ice movement and temperature d i s t r i b u t i o n i n . g l a c i e r s and i c e sheets. Journal of Glaciology, 2_, 523. Robin, G. de Q. (1969) I n i t i a t i o n of g l a c i e r surges. Canadian Journal of Earth Sciences, 6^, 919. Stanley, A.D. (1969) Observations of the surge of Steele G l a c i e r , Yukon T e r r i t o r y , Canada. Canadian Journal of Earth Sciences, 6, 819. Tarr, R.S. and Martin, L. (1914) Alaskan G l a c i e r Studies. National Geographic Society, Washington, D.C.. Thorarinsson, 3., (1969) G l a c i e r surges i n Iceland, with s p e c i a l reference to the surges of BruayHkull. Canadian Journal of Earth Sciences, §_y 875. Weertman, J . (1957) On the s l i d i n g of g l a c i e r s . Journal of G l a c i o l -ogy, _3, 33. Weertman, J . (1962) Catastrophic g l a c i e r advances. Publications of the International Association of S c i e n t i f i c Hydrology, 58, 31. Weertman, J . (1964) The theory of g l a c i e r s l i d i n g . Journal of Glaciology, 5_, no. 39, 287. - 69 -Weertman, J . (1967a) S l i d i n g of nontemperate g l a c i e r s . Journal of Geophysical Research, 72, no. 2, 521. Weertman, J . (1967b) An examination of the L l i b o u t r y theory of g l a c i e r s l i d i n g . Journal of Glaciology, 6_, no. 46, 489. Weertman, J..(1969) Water l u b r i c a t i o n mechanism of g l a c i e r surges. Canadian Journal of Earth Sciences, 6^ 929. Weertman, J.. ..(19 71a) In defense of a simple model of g l a c i e r s l i d i n g . Journal of Geophysical Research, 76', no» .26, 6485. Weertman,. J . (1971b) General theory of water flow at the base of a g l a c i e r or i c e sheet. (preprint) Reviews of Geophysics and Space Physics, 1Q, no. 1, 287. Wilson, A.T. (1969) The c l i m a t i c e f f e c t s of large-scale surges of i c e sheets. Canadian Journal of Earth Sciences, 6^, 911. Post, .A. (1971) personal communication to Dr. Garry Clarke. - 70 -APPENDIX: L i s t i n g s of Computer Programs - 71 -ROUTINE MASTER PBOSSAK *5iJ RGB* (E-DITIOH FOUR) I H P L I C i r R E i\ L * 8 (A - a,0 - Z ) HBAL*li FLOAT DIME S3 I OS Y 8 I E E ( 1 0 0 ) , r i»TCB{100) ,Y8BED(100) ,TSiBED (100) YICE (500) , T1CE (503) //BED (500) , IBED (500) XI (500) , THEPI (500) , ZBI (500) , THET3 (500) XY.T"E(500) , XTICE(50D) wi(IOOO) , x a i (1000) , i w i ( i o o o ) ,THETW (1 OOO) ,rBKP ( i G UPS (1000) , OH EG (1003) DIMENSION D E lij I Oil DIMENSION DIMSN5 ION DIMENSION COMMON DE-MS , DELTA, CO, C 1,C2 READ READ READ READ READ READ AD HEAD READ BEAD READ (5, 10 1 , EMD = &0) KCODE,XO, 12,, SL3PE, TSU&F,GE33 (5, 102) (5, 102) (5, 102) (5, 103) (5, 103) (5, 103) (5, 10*4} (5, 104) (5, 104) (5,104) R,ACCUM„ Y F I X / i P A R ^ T F I X DTIHE, DWAII E F F I C H I E E , R BED, l i l C E, :1 BSD IiST,IHTB,ISR,ISRLS IJJTW, RFR FJQ , H PRE 3 (YKICE (I) ,1=1 ,NICE) (TKICE (I) ,I=1,Hi:E) (Yi^BSD ( I ) ,1=1 j ~i\ B E D) ( T i B S D (I) ,1=1, S3 ED) n a AX r TIMEHK, rZ W EMQ WRITE (6 , 20 1) WRITE (6,2 02) 'ti Z 0 D E L-FFIC= E F F I C * 1 .02 WRITS (6 , 203) O , SLOPE, TSJ.RF, 31203, S , E F F i : EFFIC= EFFIC*1.D-2 W P A R EPSLN= 1.D-07 DSS5= 9.14D2 COSD=1.59D7 EOUIV-- '4 . 18D-5 G= 9.804D-5 SPUE AT = 5.D2 HI. A T = 8.0D4 ~ K= 3 . ' 4 DO E0= 10.3/(S*R) -1 = (1 . 03D-2) * { ( R * B / C K ) **2) C2= C1/C0 D.IFF= C0t'!D/(SP HEAT^DEMS) DI F F B - D l F F PI= 1. 14 1592653539793 - 72 -SLOPE 1= PI*SLOPE / (1 80. DO) DELI&= G*DSIN(SLOPE 1) GAMMA- G* DCOS{SLOPE 1 ) TFIX = 3 E 0 G * Y F I X + T F I X rSMPX= r FI:<-TS DRF MWAir= M I C E * K B E D KIP1 = MICE+1 HBP1 = M3EDM HWP1 = NWAIT+1 FLGI3= DBLf. { FLOAT (ISR) ) F L O T S B « DBLE(FLOAT (ISSB)) FLOPH= DBL R{FLOAT(MICE)) FLOTi1B= D B L E ( F L O A T (MB.ED) ) F LO I HV7= DBL E ( FLO AT ('K W A11) ) F L O T U P F D B L E ( F L O A T (MM?1) ) DXI = 1. DO/FLOTH DXBI = 1 . DO/FLGTiiB DXWI= 1.DO/FLOTMK D IFF V- = (FL GT H * DIF F * F LOT n 3 * DIF F B) /FLO'fHM D T AO - -4. DO*DTI*;E*DIFF/ ( Y O * Y 3 * F L : : T 5 ) DTAUB= H . D0*DTIME*DIFFI>/ (YFIX*YFIX*FLGrS3) D T A0 H-~ DR A IT* D I FF ! ' // { (I0*¥FIiC) **2> R A T01 = Di'AU/DXI RAT02-= R A P 0 1 / D1 I R& T O B 1 = DTAU3/DXBI R ATO B2 = R A10B1/DX B I RA?OW1= DTAUW/DXWI RATOl-/2= EATOW1/DXSI I P AS 5= 0 T I K E = O.DO TAU = O.DO TA OB- O.DO TAUW= O.DO - A L L TEMPO (8 I 'J E , H I C E , x i* I C E , 2? 3 ICE, Y O , XI , YICE, TICS) WRITE (6 ,320) WRITE (6,300) ( x ICE ( I ) , I=1 ,MIP1 , IK I ) ,Y ICE (H IP1) WRITE (6,3 10) TIMS 3RITE ( 6 , 3 0 0 ) ( T I C E (I) , I = 1 , H I P 1 , O T ) , r r . C E ( H I P 1 ) CALL GRAY (YICE,TICE,GRADG) CALL T E M P O ( H B E D , N S E D , Y N B E D , ' I S B E D , Y F l X, XBI , I B E D , T3SD) WRITE (6,3 2 1 ) - 73 -WRITE (6 ,300) (YBED (I) , 1 = 1 ,MB P1 , INT B) , YBED (HBP? ) WHITE (5,311) TIME WHITE (6,3 00) (TBED (I) , 1 = 1 ,H3 PI , IS I E) , I B E D ("d B P1) CALL GRAY(YBED,TBED,GR A DR) DO 20 1= 1,HIP1 20 XYICE ( I )= XI (I) *Y0/2. BO WRITE (5 , 330) WRITE (6,300) (XYICE (I) ,1=1 ,H1P1,IHT) , X Y I C £ ( K I P 1 ) 10 IPAS5= IP&SS+1 KSTEP= 0 MSTSP= 0 HI5TEP= 0 MBSTEP= 0 • x=xo 11=0. DO V=0.DO a=o.oo TSUS= r i H E TAU= O.DO 1 All B= O.DO .WLSVEL = O.DO WLEVEL= W L E V E L M . 02 WRITE (6,301) TIL1S, X, Y ,0 ,V ,WLEVEL fifLEVEL= WLEtfELM. D-2 2 ALL 3PACE(MICE, MICE,YICE,TIC E,5CY ICE , XTICE,TSORF,THETI) CALL SPACE(MBBD,MBED,YBED,TBED,Y BED,T BED,TFIX , THE! B) IF (IPA33 . L E . 1) GO TO 2 CALL GKAY (XYICE,XTICE,GRftDG) CALL GRAY (YB ED, TBED, GRADE) 2 £ S T E P = NSTEP+1 -~\ CALL SLIDE ( Y , i L E VEL, BSS,[J) V= - (Y*'J+W) /X 3 = EFFIC*U*BSS/EQITIV DGR AD= Gfi&DR-GRADG ' IISAWk Y=COND*DGRAD W= WPAB*H+(Q-23HD*D3HAD)*DTiaR*X*l.D-3/HLAT WLE VEL = W/X IF (¥, . G T . 0 . DO) GO TO 3 - 74 -END OF SURGE . . GLACIER QUIET PHASE WRITE (o ,800) TIME TSTP= TIME TSIN = FIME-TSUE WRITE (6,305) T I M E , X , Y , U , 7 WRITE (6,310) TIME WRITE (6,300) (XTICE (I) ,1=1 ,3IP1 ,11)11) ,XTICE(MIP1) WRITE (5,311) TIME WRITS (6,300) ('IBED (I) ,1=1 ,MBP1 ,INT 3) , PBED(ilBPI) WRITS-(S,805) IPASS, ISIS ESTABLISH N EH GRID SYSTEM FDR D } il il A N P GLftCIER DO 11 I = 1,P!WP1 i-JI(I)= DBL E( FLOAT (I) ) XVII (I) = (WI (I) -1 .DO) *DXKT y HI (I) = (YO+Y FIX) *XWI (I) -YO DHES (I) = O.DO 1 1 OPS (I)= 0 . DO WRITS (6,75 0) WRITE (6 ,820) WRITE (6,300) (YBI (I) ,1=1 , HWP1 ,INTW) , YHI (MS PI) IYW= 0 IYZ= 0 TAUW= O.DO CALL 3PLINE(SIP 1 , 0 , E P 5 L N , X ¥ I 3 E , X T I C E , A R G , S S , S S 1 , S S 2 ) DO 15 I="I,KWP1 Y ¥ = YVJI (I) IF (YY . G E . - 1 . D - 4 ) TEMP (I) = TBED (I- MICE) IF (YY . LT . - 1 . D - 4 . AND. YY . G E . -Y 0/2 . DO) GO TO 12 IF (Y Y . L T . -YO/2.D0 .AND. YY .GT. -Y) GO TO 13 IF (YX . L E . -Y) TEMP (I) = TSORF GO TO 15 12 CO NT I MLJ E A R G = -YY C ALL 3PLIMV (ARG, S P, S P 1 , S P 2) TEMP (I) = SP GO TO 15 13 IYZ= IYZ+1 "> IF (IYS . G E . 1) GO TO 1 '4 1 CALL 3PL INE(MIP1 ,0 ,EPSLN,Y I3E ,T ICS , A 8 3 , S S , S S I , SS2) 14 IYH= IYW+1 A R G = YO-Y-YY CALL SPLINV(&RS,SP,S?1,SP2) TEMP (I) = SP . IF (IYZ . G T . 1) GO TO 15 - 75 -HIYZ= I ' 15 THE TVs (I) = (TEMP (I)-TSURF)/TEHPX T8KP (MIYZ) = (TEMP (MIYZ- 1)>TE?1P (HIKZ+1) )/2.D0 TEMP (HIP1) = (TEMP (HI CE+ 2 ) + TEf'i P (MICE) ) /2. DO THE!iv (MIYZ) = (TEMP (I1IYZ) -TS3S?) /TS HPX THETW (NIP1) = (TEMP (MIP 1)-TSU RF) /TEMPX WRITE {S , 8 10) TIMS I PRkY TO GOD FOR THE SUCCESS OF I HIS COM PUT ES' PROGRAM. WRITS (6,300) ( T E E P (I) , I = 1, fit 3 P 1, IN I W) , T.EHP (51 HP 1) 21 HSTEP= HSIEP+1 TIMS- riME*DWftIT TAUW= rAO W + DTA UH THSTW (1) = O.DO rHETK(MWP1)= 1.D0 THSTW (MW&IT) = ( (TFIX- DXWI*GEOG) -T5 0RF) / I EMPX CALL CRANK(DTAUK, RAT031,RftTD»2, MKAIT, Pa0w,XWI, Td3TR,UPS,0HE3) DO 15 I=1,M£P1 16 TEI1P (I) = TEL'iPX*THET\i {.I) * T S U RF DO 17 I=1,frlIP1 17 TICS (I) = TEMP (MIP1+1-1) DO 16 1=1,MB?1 13 TBED (I)= TEM P ( I + KICE) • TIHB = TAD W* ( (Y0 +Y F l X ) * * 2 ) / D I F F W IF (IIMS .LT. TIHCHK) GO TO 21 TEH3ED= TEMP (HIP1 ) IF (T EMBED . LT. O.DO) GO TO 22 TSAR= TIME-TSTP ACCUM= (YQ-Y) /TSAR WRITE (6,850) TIME TSAR ACCOM TIME (TEMP (I) , 1= 1 , H« P1 , INT W) , r SE P (MH P1) WRITE (6 ,86 0) WRITE (6,87 0) WRITE (5 , 8 10) WRITE (6,3,00) G O T O 10 22 CONTINUE IF ( (MS T EP/K FREQ) *MFREQ WRITE (6,8 10) TIMS WRITE(6, 300) «RITE (6,500) IF (TIMS . LT. WRITE (6,84 0) GO TO 1 ,NE. 3 STEP) GO TO 21 (TEMP (I) TEMBED TIM END)' 30 TO 21 TIHS I=1,HvJP 1 , INT W) r EH p ( a w P 1) - 76 -3 CONTINUE. C CALCULATE BEDROCK TE MP ERATH RES DO n j=1,MBP1 3ME3 (J) = O.DO 4 OPS (J) = 0. DO f w 5 MBS T EP = MBSTEP+1 TAUB= TAO'B + DTAUB THSTB(1) = O.DO THETB (MBED) = 1 . D0-DX8I*GE()G/IFIX THET B (MBP 1) = 1.D0 CALL CR A W K (DTA U8, R AT OB 1 , RAT0B2 , MBSD, T A'J B ,XB I, T HET B,OPS, OM EG) IF { {MBS T EP/ISRB) * IS R 3 . N E. SBSTE?) 30 TO 5 DO 6 J=1,HBP1 6 TBED (3) = T F l X T HE IB {J) CALL GRAY(YBED,TBSD,GEADR) ..C CALCULATE GLACIER TEMPERATURES . CALL 5PLINE(HICE,0,EPSLN,Y8ICE,TNICE,ARG,SS,SS1, SS2) DO 7 J=1,MIP1 OMEG (3) = 0.DO 7 UPS(J)= YO*V*XI( J) / (H. DO*DIFF) C . 8 KISTEP= MISTEP+1 TAU= TAU+DTAU AR3= (1.DO-UPS(KIP1 ) *TAU)*Y0/2 . DO C ALL 3PLINV (ARG,SP,SP1,3P2) THETI (1)= 0.DO THETI(MIP1) = SP/TSORF CALL CRANK (DTA0,RATO1,RAT02,MICE,T AO, X I , THETI, UPS,OMEG) IF ( (MI3 I EP/ISR) * ISR . NE. MISTEP) 30 TO 8 DO 9 J=1,MIP1 9 X T I C E ( J ) = TSURF*THEn (J) CALL GRAY{XYICE,XTICE,GRA DG) C TIME= TIME+DTIME X= X + U*DriME •£= Y + V*DTIME IF ( (M3PEP/NFREQ) *NF8EQ .ME. SSTS?) 30 TO 2 WLE VEL= «LEVE.L*1 . D2 WRIT 2 (5,30 1) TIME, X, Y, 0,7,WL.E7EL ffLEVBL= WLEVSL*1.D-2 WRITE (6,390) T IM E, Q, 38 A DR , GR A DG, OS A WA Y , W BRITE (6,310) TIME WRITE (6,300) (XTICE (I) ,I=1,!JIP1,IST) , XTICE (HIP1) WRITE (6,311) TIME c - 77 -WRITE (6 , 300) WRITE (6 ,600) IF (TIKE . LT. WRITE (6,700) GO TO 1 (TBED (I) ,I=1,M3?1, INTB) ,TBED{MBpi) TIME T IM AX) GO TO 2 101 132 103 134 201 23 2 20 3 300 33 1 335 31 D 311 32 3 321 333 39 0 500 63 3 700 753 800 835 810 823 840 850 860 870 60 FORMAT STAT EM EiSTS FORMAT (I10,5D10.4) , FORMAT (5 D10. 4) FORMAT (414) FORMAT (16 D5. 2) FORMAT (1111 ,'PROGRAM *SURGE* '///) FORMAT ( 8 GLACIER CODE NO. = «,I5/) FORMAT (« I N I T I A L VALUES'/ 0 X= » , 6 X , *¥=* , 5 X , ' SI.OPE = 8 , 2X , 5 TSURF= 4 , 2 X, 1 » GTG = * , 4X , ' ROUGH= 8, 2X, ' E.FFIC = 1 1 , 2X, 3 WP AR= 8/2F8 . 0 , 6F8. 2///) FORMAT (8F8.2) TIM E= ', F 6. 2/ ' X= », 6X, 8 Z=» , 6X, "0= ', &X, 1 V=« , 6 X , 'HLE V= V 2 P FORMAT (// 18.0,2F8.2 FORMAT (» F8.2//) END OF SURGE V A L U E S : V 8 1, 6X,'V= 8,6X/2F8.0,2F8.2//) FORMAT (// FORMAT (// FORMAT (// FORMAT (// FORMAT (// FORMAT (// 15D11.3) FORMAT(// FORMAT (/« FORMAT (' FORMAT (// FORMAT (// FORMAT (/' FORMAT (// FORMAT (// FORMAT (// 1 r o . 2////) FORMAT (// 1 g. TIM E = r i M E = « , F 6 . 2 / ' X = « , 6 X , 8 Y = 8,6X, 'U = I C E TEMP 3/ 5 PI:1E= «,F6.2) BED TEMP'/' r i K E = ' , F 6 . 2 ) I C E TEMP ARE GIVES AT THE BED TEMP ARE GIVES AT I K E ICE TEMP ARE tl08 GIVSN AT TIM E= 8,F6.2/' Q = « , 9 X , ! 3 a = ' FOLLOWING L E V E L S : 5 ) FOLLOW IN G DEPTHS: s) THE FOLLOWIS3 LEVELS: , 8 X, e GG=',8X,'K*DG=» *) 6X,*fl= ICE-ROCK CONTACT TEMP=" >F5. GLACIER IS S L I D I N G . . . T I M E = » , GLACIER DOES NOT REFREEZE TO M EC-J GRID SYSTEM FOR D 3 3 MA MP 2) F6 .2/) BED...TIME=',F6.2//) GLACIER HAS BEE w CREATED 5) END OF SURGE T I E E = « , P 6 . 2//) SURGE WO. ',12,' HAS LAS T E D 4 , F 6 . 2 , 8 ICE & BED TEMP 5/' TIME= « , FS . 2) TEMP ARE GIVEN.AT THE F OLLQl-M hi G DEPTHS GLACIER DOES NOT 3 ARM TO MELTING POINT YEAUSV) ' ) AT BED. ..TIME=' FORMAT (/' 1SORGS') FORMAT (/' CONTINUE STOP GLACIER , F6. 2//) THE GLACIER REACHES MELTING POINT AT BED AND BEGINS TO SLID HAS HAITED « , F 6 . 2 , ' YEARS SINCE BSD OF LAST ACCUMULATION RATE FOR DOR MAST GL AC IER= 3 , F5 . 2 , 8 'A/IR"//) ' - 78 -C SUBROUTINE T3RP0 (GENERALIZED I N I T I A L VALOES) (*5 U RG E*) SUBROUTINE TEHPO(N,INIT,YINIT,TINIT,YMAX,X,Y,T)' ' IMPLICIT REAL* 8 (A-H,0-Z) R E AL*U FLOAT DIMENSION YINIT(IOO) , TINIT (100) , Y (500) ,X (500) ,1(500) NP1= N+1 FLO A"N = DBL E(FLO AT (N)) EPSLN= 1 . D - 0 7 CALL SP L I N E ( I N IT, 0, EPSLN/YIMIT, TINIT, ARG,SS, SS1 , SS2) C DO 11 1=1,NPI' DZ= DBLE ( F L O A T ( I ) ) Z- Y PI A X * (DZ- 1. DO) /FLOAT N . ARG = Z CALL SPLINV (ARG, SP,SP1,SP2) T (I) = SP X (I) = Z/YMAX 1 1 1 (I) = Z RETURN - 79 -SUBROUTINE GRAY: CALCULATE APPROX. TEKP. GRAD'. USING GRID VALUES SUBROUTINE GRAY(Y,T,GRADNT) IMPLICIT REALMS (A-H, O-Z) DIMENSION Y (1000) ,T (1000) GRADNT = (T (2) -T (1 ) ) / (Y (2) -Y ( 1) ) IF(GRADNT .GT. O.DO) GRADNT=-GRADNT RETURN END - 80 -SUBROUTINE S.^IDE (CALCULATES WEES! MAN-KAMB SLIDING VELOCITY) SUBROUTINE SLIDE (Y , WLEVEL ,BSS , U) • IMPLICIT REftL*8 (A-H, O-Z) COMMON DENS,DELTA,CO,C1,C2 BSS= DENS*DELTA*Y SLA K= CO/BSS U= C I * (BSS * * 2 ) IF (WLEVEL .GT. SLAM) U = C2*BLEVEL*(ESS**3) RETURN END - 81 -SUBROUTINE SPACE: RESETS GRID VALUES (*SUR,GE*) SUBROUTINE SPACE (MIN,MOOT, YIS,TIN,YOUT, TOUT, MAX,TH ETA) IMPLICIT REAL*8(A-H,O-Z) DIMENSION TIN ( 1000) , TOUT ( 1000) , KIH (1000) ,Y0UT(1000) DIMENSION THETA(1000) MOP1= M0UT+1 EPSLN= 1.D-07 CALL 5PLINE(MIN,0,EP5LN,YIN,TIN, A i G , SS,SS1,SS2) DO 22 1=1,HOP! ARG= YOUT (I) CALL SPLINV (ARG, SP,SP1,SP2) TOUT (I) = SP THETA (I )= TOUT (I)/TM AX 22 CONTINUE RETURN END SUBROUTINE CRANK-NICHOLSON (ABRIDGED VERSION) (*S0RGE*) SUBROUTINE CRARK (DTAU,RAT01 ,RAT02, M,TAU ,XI,THETA, II PS , OM EG) IMPLICIT REALMS(A-H,O-Z) DIMENSION A (1000) , B (1000) ,C (1003) , D (1000) DIMENSION THETA (1000) , XI ( 1000) , UPS (1000) , OS EG {1 000) MP1= M+1 COMPUTE ARRAYS A, B, C, AND D A (1) = O.DO B (1 ) = 0. DO C (1) = O.DO D (1) = 0. DO DO 3 3 .1=2 ,M A ( I ) = -RAT02 B (I) = 1 . DO-UPS (I) *RAT01 + 2. D0*RATO2 C (I) = RAT01*UPS (I) -RAT02 D ( I ) = RAT02*THETA (I- 1 ) + (1 . DO -2. DO*B AT 02) *THETA {I) *R AT02*THETA (1 + 1) D (I) = D (I) -DTAU*OMEG ( I - 1) D{2) = D (2) +RAT02*THETA (1) D ( M ) = D (M) + (RAT02-UPS (MP1) *RAT01) *THETA (S P 1) COMPUTE SOLUTION VECTOR CALL TRIDAG (2,N,A,B, C, D,THETA) RETURN END - 83 -C • SUBROUTINE SOLUTION FOR TRIDIAGONAL SYSTEM OF LINEAR EQUATIONS SUBROUTINE TRIDAG (II,IL,A,B,C,D,VECI) IMPLICIT REAL*8 (A-H, O-Z) DIMENSION A (500) , B (500) ,C (500) , D (500) , V ECT (500 ) DIMENSION BETA(500),GAMMA(500) C COMPOTE INTERMEDIATE ARRAYS: BETA AND GAMMA BETA (II) = B (II) GAMMA (II) = D (II) /BETA (II) I IP1= I I H DO 1 1= I I P 1 , I L BETA (I) = B (I) -A (I) *C ( I - 1J/BETA (1-1 ) 1 GAMMA (I) = (D (I) -A (I) *3AMMA (1-1) )/BETA (I-1) C COMPUTE SOLUTION "/ECTOR: VECT V E C T ( I L ) = GAMMA (IL) LAS T= I L - I I -DO 2 K= 1,LAST 1= IL-K ' 2 V E C T ( I ) = GAMMA (I) -C (I) *VECT ( I + 'l)/BETA (I) RETURN END - 84 -S U B R O U T I N E S P ' L I H E / S P L I M I N T E R P O L A T I O N PROGRAM (*SURGE«) S U B R O U T I N E S P L I N E (N, M, E P S L N , f t , ¥ , A R G , S S , S S 1 , S S 2 ) I M P L I C I T REA L * 8 (A-H,O-Z) DIMENSION X (50 0) , Y ( 5 0 0 ) , T (500) , H (50 0) , DEL Y (500) ,H2 (5 00) ,D E L S Q Z ( 5 0 0 1) ,S2 (500) ,S3 ( 5 0 0 ) , 8 (5 0 0 ) ,C (500) 50 N-1= N-1 H( 1 ) = X ( 2 ) - X ( 1 ) D E L Y ( 1 ) = (Y (2) -Y (1) ) /H (1) 51 C O N T I N U E DO 5 2 1= 2,NI H ( I ) = X ( I + 1 ) - X ( I ) H2 ( I ) = H ( I - 1) + H ( I ) B ( I ) = ( 5 . D - 0 1 ) *H (1-1 )/H2 ( I ) D £ L Y ( I ) = (Y (1+1) -Y ( I ) )/H ( I ) • D E L S 3 Y ( I ) = ( D E LY ( I ) - D E L Y (1-1) )/H2 ( I ) S2 ( I ) = (2. D 0*DELSQY ( I ) ) 52 C ( I ) = (3 .DO) * D E L 5 Q Y ( I ) S2 (1) = 3 2 ( 2 ) S2 (N) = S2 (N-1) DM EGk- 1 . 0 7 1 7 3 6 8 D 0 53 ETA= 0.D00 54 DO 60 1=2,N1 55 W= (C ( I ) - B (I) *S2 (1-1 ) ~ (5. D-01*B (I) ) *S2 (1+1) -S2 ( I ) ) *OMEGA 56 I F (DABS (W) . L E . ETA) GO TO 60 57 ETA= DABS(W) 60 S2 ( I ) = S2 ( I ) +W 6 1 I F (ETA .GE. E P S L N ) GO TO 5 3 6 2 DO 6 3 I=1,H1 6 3 53 ( I ) = (S2 (1 + 1 ) - S 2 ( I ) )/H ( I ) AS= DELY (HI) +H (H 1) * 5 2 (N 1) / 6 . D00 A L= DELY ( 1 ) - H (1) *S2 (2 )/6 . D00 I F (M . L E . 0) RETURN GO TO 64 ENTRY S P L I N V ( A R G , S P , S P 1 , S P 2 ) M= 1 61 I-- 1 65 I F (ARG - X (1) ) 69, 73, 66 66 I F ( A R G - X ( N ) ) 6 8 , 7 1 , 7 0 67 I F ( A R G - X ( I ) ) 72, 73, 6 8 68 1= 1+1 GO TO 67 6 9 SS= A L * ( A R G - X ( I ) ) + Y (1) SS1 = AL - 85 -SS2= O.DOO : i GO TO 74 70 SS= AB* (ARG-X (S) ) +Y (-N) " S31 = AR SS2= 0.DOO GO TO 74 7 1 1= N 7 2 I=. 1-1 7 3 KT1= ARG-X(I) HT2= ARG-X(I + 1) PRODfi= HT1*HT2 SS2= S2 (I) +HT1*S3 (I) DELS35= (S2(I) +S2 (1 + 1) +SS2) /6.D00 SS= Y (I) + HT1 *DELY (I) +PRODH*DELSQS SS1= DELY (I) * (HT1+HI2)*DELSQS+PR3DH*S3(I) /6.D00 7 4 CONTINUE SP= 33 SP1= SS1 RETURt-i END