SLOW SURGE OF TRAPRIDGE GLACIER, Y U K O N TERRITORY, 1951-2005 by T O M - P I E R R E F R A P P E - S E N E C L A U Z E B.Sc. physique, Universite Laval, 2002 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Geophysics) T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A February 2006 © Tom-Pierre Frappe-Seneclauze, 2006 Abstract Trapridge Glacier, a surging glacier located in the St. Elias Mountains, Yukon Territory, Canada, went through a complete surge cycle between 1951 and 2005. Air photos (1951-1981) and ground-based optical surveys (1969-2005) are used to describe the modifications in flow and geometry that occurred over this period. The acceleration of flow during the surge is detected by repeated measurement of poles drilled into the ice. Between 1974 and 1980, the median velocity in the lower basin went from 1 5 . 6 m a - 1 to 3 8 . 6 m a - 1 . Downstream from this zone of active flow, cold-based ice accumulated during the previous surge impeded the flow, and a steep front formed at the boundary between the two ice masses. Over the following ten years, this bulge propagated downglacier, advancing faster than the ice and integrating stagnant ice by continuous deformation. After it peaked in 1984, the flow in the lower basin remained above 2 5 m a - 1 for 12 years, but was on a slowing trend. The slowdown followed strangely regular 4-year pulses: 1-2 years of timid acceleration ( 2 - 3 m a - 1 ) , followed by 2-3 years of rapid deceleration ( 4 - 8 m a - 1 ) . The 1997-1999 acceleration was particularly vigorous, as the median velocity went from 2 0 . 3 m a - 1 to 2 8 . 5 m a - 1 . After this last pulse, the glacier gradually slowed down to pre-surge velocities. In 2005, the lower basin was flowing at less than 8 . 5 m a - 1 . Based on this flow history, I divide the 1969-2005 period into three parts: the activation phase (1969-1974), the surge (1980-1999), and the waning phase (1999-2005). To quantify the geometrical changes that occurred during each phase, digital elevation models are constructed from air photos and optical survey measurements. Optical and radar surveys are joined to photogrammetric measurements of the proglacial field to obtain a bed topography map. D E M s for 1951, 1970, 1972, 1977, and 1981 are generated by stereographic analysis of air photos. These background models are then updated year after year by ground-based survey data, using a Bayesian kriging algorithm. By assuming the ice surface to be similar in shape from one year to the next, only shifted vertically, this algorithm allows the propagation and diffusion of information on the ice surface topography. The interpolation obeys the data, and follows the prior model where no data are available. Uncertainty is attributed to each estimation, based on the estimated covariance structure of the field and on the uncertainty of the prior model. The assumption of stationarity and the parametrization of the covariance are tested by the analysis of orthonormal residuals. Despite the significant changes in surface area that occurred between 1969 and 2005, the total volume of ice remained relatively constant at 0.14 ± 0.03 k m 3 . The spread of the glacier during the surge was accompanied by a corresponding decrease in thickness. From the surface slope and ice thickness fields, the stress driving the flow is calculated. After filtering the surface slope to eliminate small-scale variations, the shallow ice approximation is used to estimate creep velocities from the stress and ice thickness fields. Deformation of the ice is shown to accounts for a negligible portion of the flow observed in the lower basin. Sliding or bed deformation must therefore be responsible for most of the motion in that area. In the upper basin, thicker and flowing more slowly, creep accounts for maybe as much as half of the recorded motion. iii Contents Abstract ii Contents iii List of Tables vi List of Figures vii Acknowledgments ix 1 Introduction 1 Literature review 4 2.1 Surge-type glaciers 4 2.1.1 Environmental controls on surge-type glacier distribution 6 2.1.2 Surge dynamics: overview of existing models 7 2.2 State of research on Trapridge Glacier 10 2.2.1 1940s surge 10 2.2.2 Thermal regime and flow evolution 1969-1989 10 2.2.3 Subglacial hydrology and ice-bed mechanics 11 Photogrammetric analysis: 1951-1981 ice surface topography, and bed topogra-phy . . . 14 3.1 Construction of D E M s from aerial photography 14 3.1.1 Testing of the D E M s and evaluation of uncertainty 16 3.1.2 Resulting D E M s 17 3.2 Bed topography 19 3.3 Results: Ice thickness maps, outlines, cross-sections 22 Ice flow calculation and estimation of glacier outlines 25 4.1 Description of the data set 25 4.1.1 Measurement accuracy 26 4.1.2 Expunging of blunders 28 4.2 Results: Changes in flow velocity 30 4.3 Estimation of glacier outlines 33 Kriging 36 5.1 Introduction: geostatistics and the interpolation problem 36 5.2 Ordinary kriging 38 5.2.1 Structural assumption 38 5.2.2 Unbiasedness . . 39 Contents iv 5.2.3 Estimation variance '. . 39 5.2.4 Normal equations 40 5.2.5 Optimal estimation error variance 40 5.3 Model selection: the experimental variogram 40 5.4 Model validation: analysis of orthonormal residuals 41 5.5 Application: interpolation of bed topography data 43 5.6 Including a priori information: Bayesian kriging 47 5.6.1 Structural assumptions 47 5.6.2 Bayesian block kriging algorithm 48 5.6.3 Model selection 50 6 Time-evolving D E M 52 6.1 Description of the algorithm 52 6.2 Incorporating information about the glacier limits 55 6.3 Selection of covariance parameters and orthonormal residuals analysis 59 6.4 Results: Changes in total ice volume, glacier profile, and ice distribution 64 6.5 Results: Driving stress 69 6.6 Results: Creep velocity and sliding 71 7 Conclusion 78 7.1 Flow , . . . 78 7.2 Geometry 79 7.3 Discussion and outlook 80 Bibliography 84 A Photogrammetric models of the Rusty Creek basin: technical issues 88 A . l Procedures for the orientation of air photos and collection of elevation data 88 A . 1.1 Orientation of air photos • • • • 88 A . l . 2 D E M surface collection 91 A.2 Other issues regarding the Trapridge D E M s : . . . 92 A . 2.1 Importing and gridding D E M s using M A T L A B 92 A . 3 Other D E M results 100 B Survey data import and filtering summary 102 B. l Accessing survey data 102 B.2 Filtering survey data 103 B. 2.1 By visual inspection 103 B.2.2 By flow analysis 104 B.3 Filter and segment flow-pole time series: makeallcontinuous .m 104 B.4 Calculate flow: calculatef low.m 107 B. 5 Results: Flow maps, 1970-2005 107 C Bayesian block kriging 116 C l Structural assumptions . . 116 C . 2 Preliminary definitions 116 C.3 Unbiasedness 118 C.4 Estimation error 119 C.5 Normal equations 120 Contents v C.6 Optimal estimation error variance 120 C.7 Algorithm summary and numerical implementation 121 D Stat is t ica l analysis of D E M t ime series 125 E Results : C o l o r plates 129 List of Tables 2.1 Duration of active phase of surge-cycle for various glaciers 6 4.1 Summary of flow history, 1969-2005 31 5.1 Orthonormal residuals statistics and decision rules 43 5.2 Orthonormal residuals analysis of two statistical models of the bed topography . . . 45 5.3 Bayesian block kriging system 49 6.1 Prior model used for each year of the time series 54 6.2 Orthonormal residuals analysis of conditional covariance models 63 7.1 Summary of geometrical changes, 1951-2005 83 A . l Ground control points used in the absolute orientation of air photos 90 A . 2 U T M coordinates of ground control points around Trapridge Glacier, in NAD27 . . . 93 B. l Decision parameters for the function i s p o s s i b l e f low 106 C . 2 Bayesian block kriging system, matrix form 124 D. l Orthonormal residuals statistics for all years, time-independent Bayesian model (best model) 125 D.2 Orthonormal residuals statistics for all years, time-independent second-order station-ary model 126 D.3 Orthonormal residuals statistics and parameters for all years, time-dependent Bayesian model 127 D.4 Orthonormal residuals statistics and parameters for all years, time-dependent second-order stationary model 128 List of Figures 1.1 Trapridge Glacier study site 2 2.1 Trapridge Glacier, 1941 11 3.1 Air photos of Rusty Creek basin, 1951-1981 15 3.2 Mean and standard deviation of the photogrammetric D E M s 16 3.3 Comparison of ground measurements and interpolated photogrammetric D E M eleva-tion, 1972 and 1981 18 3.4 1981 D E M s 19 3.5 Elevation difference between subsequent photogrammetric D E M s 21 3.6 Bed topography map 21 3.7 Ice-thickness maps, 1951-1981 23 3.8 Total ice volume and area, 1951-1981 23 3.9 Glacier outline and cross-section, 1951-1981 24 4.1 Flow pole grid in 2004, and partition of the glacier in zones 26 4.2 Distribution of daily displacement 27 4.3 Gradual elimination of blunders and its effect on annually-averaged flow vectors . . . 29 4.4 Zonally-averaged ice flow 32 4.5 Glacier area: 1951-2005 34 4.6 Estimated glacier outlines: 1974, 1981, 1991, 1996, 1999, and 2004 35 5.1 Data points for the kriging of the bed topography map 44 5.2 Statistical analysis of bed topography data . 45 5.3 Bed topography: detrended elevation and kriging standard deviation 46 6.1 One time step in the generation of the D E M time series 53 6.2 Ice thickness field and kriging standard deviation 56 6.3 Resolving the glacier boundary 58 6.4 Experimental conditional variograms 61 6.5 Bulge propagation, 1981-2005 65 6.6 Total ice volume and glacier area, 1951-2005 66 6.7 Profile evolution and flow, 1951-2005 67 6.8 Area under the profile • • 68 6.9 Time series of mean ice thickness, median slope, and median stress 69 6.10 Ice thickness, surface slope, and driving stress 70 6.11 Median creep rate time series 73 6.12 Observed surface velocity, creep velocity, and percentage of flow due to sliding . . . . 74 6.13 Percentage of longitudinal velocity due to sliding in zone U B , A , B , and C 75 6.14 Zonally-averaged sliding velocities, 1969-2005 76 List of Figures viii 7.1 Trapridge Glacier, 1941 and 2004 82 A . l Typical fiducial detail and layout 89 A.2 Location of control points used to reference the D E M s 94 A.3 Aerial photograph of Rusty Creek basin, 1951 95 A.4 Aerial photograph of Rusty Creek basin, 1970 96 A.5 Aerial photograph of Rusty Creek basin, 1972 97 A.6 Aerial photograph of Rusty Creek basin, 1977 98 A.7 Aerial photograph of Rusty Creek basin, 1981 99 A.8 Elevation difference from one Rusty Creek basin D E M to the next 100 A . 9 Elevation difference between each Rusty Creek basin D E M and the mean D E M . . . 101 B. l Flow maps: 1970-1982 108 B.2 Flow maps: 1983-1988 109 B.3 Flow maps: 1989-1994 . 110 B.4 Flow maps: 1995-2000 I l l B.5 Flow maps: 2001-2005 112 B.6 Glacier outlines: 1969-1986 113 B.7 Glacier outlines: 1986-1998 114 B.8 Glacier outlines: 1999-2005 115 E . l Profiles: 1951-1980 129 E.2 Profiles: 1991-1987 130 E.3 Profiles: 1988-1994 131 E.4 Profiles: 1995-2001 132 E.5 Profiles: 2002-2005 133 E.6 Ice thickness and elevation change: 1951-1974 134 E.7 Ice thickness and elevation change: 1980-1985 135 E.8 Ice thickness and elevation change: 1986-1991 136 E.9 Ice thickness and elevation change: 1992-1997 137 E.10 Ice thickness and elevation change: 1998-2003 138 E . l l Ice thickness and elevation change: 2004-2005 139 E.12 Kriging standard deviation: 1951-1985 140 E.13 Kriging standard deviation: 1986-1997 141 E.14 Kriging standard deviation: 1998-2005 142 E.15 Driving stress: 1951-1985 143 E.16 Driving stress: 1986-1997 144 E.17 Driving stress: 1998-2005 145 E . 18 Creep velocity: 1951-1985 146 E.19 Creep velocity: 1986-1997 147 E.20 Creep velocity: 1998-2005 148 ix Acknowledgments Rarely does one have the chance to print, distribute, and burn on microfiche one's gratitude for the generosity of others. In an effort towards concision, I have tried to limit my thanks to the people who have contributed directly to this thesis. However, as the wise keep reminding us, everything is in everything, and vice versa. Few people I know did not have an impact on this thesis, either directly of via butterfly wings. In 2000, after a few emails and reading an untranslated French C V , Garry invited me to Trapridge as a field assistant. In a typical example of Garry's breadth of vision, the fact that I once raised praying mantises and herded goats in Greece seemed to matter more than the fact that I had never touched a soldering iron, fixed a generator, or hooked up propane tanks in parallel. Over the following four years, Garry allowed me to take on increasing responsibilities in the logistics of our work at Trapridge, which made each field season a new challenge. He granted me the same trust to shape my degree. Even if sometimes frightening and solitary, this unbounded intellectual freedom was a blessing. Garry never asked me to restrict my curiosity to glaciology, and always encouraged me in the pursuit of various interests, whether within the department, at the school of music, or outside campus. For your guidance in all practical aspects of the academic world, for creating such a welcoming space, both in the field and on campus, and for the self-confidence I acquired by meeting challenges in such favorable conditions, I thank you, Garry. I am also grateful to the other members of my committee, Roger Beckie and Will iam Hsieh, for their encouragement and consultation. Their open door policy and dedication to the improvement of this thesis was much appreciated. Collaborating with Kelly Russell on problems common to volcanology and glaciology has been very stimulating, and I would like to thank him for these few months spent between hot rocks and cold hard places. Field work at Trapridge is a group effort, and lively evening discussions in the Weatherheaven are an essential part of the the success a field season. For their professionalism and congeniality I would like to thank: Fern Web, T i m Creyts, Sean Fleming, Neil Hall, Jessica Logher, David King, Rob Eso, Christian Schoofe, Natalie Jackson, Alex Babakov, and Andrew Scheaffer. The 2002 and 2003 seasons spent with J A L , D S K , and R A E hold a special place in my memories. The weeks of preparation for the field are spent going back and forth between our lab and the machine shop; Doug Poison, Ray Rodway, and Joern Unger have always supported us through their skillfulness, their creativity, and their over-tested tolerance for last minute requests. David Jones has been of great help in the domain of electronics, from chips to software. Every season in the Yukon starts and ends at the Kluane Lake Research Institute; thank you Andy, Carol, Cian, and Lance for your logistical support, whether for work or for play, and for making K L R S such a hospitable space. In the last year, I found myself spending more and more time in the Geophysics building, which became (sometimes quite literally) a second home. The recent arrival of Jenny Brauch in the office next to mine brightened my daily life, and I am very grateful for her generous support in the last stages of this thesis. I have enjoyed sharing ideas with the other members of the U B C Glaciology Group, and am indebted to T i m Creyts, Christian Schoofe, and Gwenn Elizabeth Flowers for their scientific and editorial advice. Acknowledgments x M y graduate studies were funded through the Natural Sciences and Engineering Research Coun-cil of Canada (NSERC) , with additional support from U B C Earth and Ocean Sciences. The excellent work of the administrative staff of the department made all bureaucratic and logistical issues simple, a quiet luxury than cannot be undervalued. Special thanks to Alex Allen and David Shorthouse for their good spirit and kind support. M y time in Vancouver has been brightened by the beautiful people I have lived with — Maya, Paul, Jim, Paula, Molly, Rebecca, Alev, Gina, Naava; you know how central you have been to my happiness. The guidance of Arnaud Desjardins, Ajahn Sona and the Birken Monastery Community, Chogyam Trungpa Rinpoche and my friends at the Shambhala Center has framed most of my evolution in recent years. M y work with Zoe Eakle, Radhasri (Rhonda Fogel) and Grant Couture has also been the catalyst of many changes, and I am most grateful for their teachings. — How should we call this? Gallicynism, the inescapable transboardment of French syntax into English. T P F S 1 Chapter 1 Introduction Fast-flowing ice plays a predominant role in the mass balance of ice sheets. It has been suggested that as much of 90% of the discharge from the Antarctic Ice Sheet occurs through a small number of ice streams and outlet glaciers (Bamber et al., 2000). Creep alone cannot account for these high discharges, and most of the displacement is thought to occur at the bed, either by sliding or by sediment deformation. Efforts to model the response of modern and past ice sheets to climate change have been hindered by the unrealistic treatment of the basal conditions controlling these localized flow pathways. Fast flow is also exhibited by a unusual type of mountain glacier known as surging glaciers. Glacier surges are a flow instability characterized by a short phase of fast flow followed by a long period of quiescence. Geological evidence has shown that surging is a cyclic phenomenon, and that the periods can be remarkably uniform. Budd (1975) suggested that the main difference between ice streams and surge-type glaciers lies in the size of their catchment area. Ice streams, tapping into large reservoirs of ice, can sustain continuously a large flux, while surging glaciers, unable to maintain the ice fluxes necessary for fast flow, oscillate between fast and slow flow. Important geometrical and dynamical differences exist between ice streams and surging glaciers, but the mechanisms of sliding and sediment deformation are common to both. The study of surge-type glaciers provides an opportunity to observe these phenomena more closely. Alpine glaciers are more accessible, warmer and considerably thinner than ice streams. Their thinness reduces drilling time and facilitates instrumentation of the bed. Detailed observation of the subglacial environment has been the main focus of the study of Trapridge Glacier, a surge-type glacier situated in the St. Elias Mountains, Yukon Territory, Canada (Figure 1.1). Over the last 31 years, numerous aspects of its dynamics have been studied: ther-mal regime, ice-bed coupling, sediment deformation, subglacial electrical phenomena, subglacial, englacial, and supraglacial hydrology, and meteorological forcing (see Chapter 2 for a review). Each research project has extended our knowledge of the subglacial environment of Trapridge Glacier, and has shed light on phenomena that are relevant to other soft-bedded glaciers. The case-by-case interpretation of field experiments and serendipitous events facilitated a conceptual understanding of the hydro-mechanical processes taking place under the glacier, but this knowledge has yet to be unified into a synoptic model allowing glacier-wide predictions on seasonal time scales. The monitoring of a surge is the perfect experimental support for the construction of such a model. The last synoptic description of the flow of Trapridge Glacier dates from the late 1980s, when a new surge seemed imminent. In spite of promising debuts, the advance did not evolve into a catastrophic surge similar to that observed in the 1940s. As the wait for a surge dragged on, Chapter 1. Introduction 2 F I G U R E 1 . 1 : Trapridge Glacier study site. The topographic map is from 1981. Shading represents ice or snow. The rectangle locate the area where most subglacial instruments are installed. The inset shows the location of the glacier in the western Yukon Territory, Canada. Chapter 1. Introduction 3 interest in modeling flow evolution on a glacier scale waned. The primary goal of this thesis is to reopen this branch of research, by drawing attention to the significant changes that occurred over the last 31 years. I will present evidence that a shift in the flow regime occurred between 1972 and 1981, leading to a ~20-year period of glacier advance and increased surface velocity. This phase of fast flow, admittedly not as dramatic as the previous surge, considerably modified the geometry of the glacier. The advance terminated gradually over a period of ~4 years, starting in 1999. The glacier has now returned to pre-surge velocities, and is entering a new phase of quiescence. To reconstruct the flow and geometrical evolution of the glacier over this complete surge cycle (1951-2005), I use optical survey measurements made in the field, and aerial photography. Air photos provide early information on the state of the glacier. Digital elevation models (DEMs) are obtained by stereographic analysis for 1951, shortly after the end of the previous surge, and for the first years of the field study (Chapter 3). The repeated surveying of flow poles inserted into the ice allows the calculation of mean annual surface velocities from 1969 to present. Using these data, important variations in flow rates are detected, and the advance of the terminus is measured (Chapter 4). Survey measurements also provide information on the ice surface topography. To generate D E M s , these data must be interpolated, and often extrapolated. Chapter 5 presents the theoretical development of ordinary and Bayesian kriging, two best linear unbiased estimators used for the analysis of spatial data. In a first application of ordinary kriging, optical and radar surveys are joined to photogrammetric measurements of the proglacial field to obtain a bed topography map. Bayesian kriging is used to generate a time-evolving D E M of the ice surface (Chapter 6). From this time series, changes in volume and aspect can be quantified, and the creep velocity can be calculated. By subtracting the contribution of creep to the observed velocity vectors, sliding velocities are estimated. 4 Chapter 2 Literature review 2.1 Surge-type glaciers Surge-type glaciers exhibit periodic changes in annually averaged flow rate, with alternating phases of slow flow (quiescence) and fast flow (surge). In their foundational 1969 paper, Meier and Post identified 204 surging glaciers in western North America, and outlined the following characteristics of surge behavior. (1) Surges occur repeatedly, and the duration of the cycle is generally constant for a single glacier. (3) Surges are short-lived (< 1 to 6 years, most commonly 2-3 years) and quiescent phases are much longer (~15 to > 100 years, commonly 20-30 years). (4) The surge-phase average ice velocity can be 10 times (or more) that of the quiescent phase. During the surge, ice drains rapidly from an upglacier reservoir area to a downglacier receiving area, causing thickness changes of 10-100 m and horizontal displacements of < 1km to ~ 1 0 k m ( 1 0 - 1 or more of the glacier length). Accumulated ice displacement is 10 to 100 times greater than the total quiescent phase displacement. (5) During the quiescent phase, ice accumulates in the reservoir area and melts in the receiving area, thus gradually reversing the changes that occurred during the surge. The glacier eventually returns to near its presurge state. This description of the surging phenomenon as observed in Alaska and the Yukon is considered as the classic surge-type behavior. However, information on polar surging glaciers is sparse, and some evidence shows that they could have surges of a different kind. Dowdeswell et al. (1991) examined 8 surge-type glaciers of the arctic archipelago of Svalbard and compared these to 36 glaciers from other areas of the world. The study was limited to the glaciers for which the active phase duration was known. They found several important differences. The glaciers in Svalbard have a much longer active-phase, from 3 to 10 years, compared to 1-2 years for most other glaciers. Ice velocities during the active-phase are also considerably lower. The termination of the surge is gradual, with velocities decreasing over several years, while most Alaskan glaciers have very abrupt surge termination, with dramatic drops in velocities over a few days (Raymond, 1987). The quiescent phase of Svalbard glaciers is also much longer, 50-500 years by comparison to a typical 20-30 years for most other surge-type glaciers. Table 2.1 reproduces the data gathered by Dowdeswell et al. (1991). The Svalbard glaciers are all polythermal (their temperature varies with depth), while most but not all the others are temperate (isothermal, at melting point). Chapter 2. Literature review 5 T A B L E 2.1: Duration of active phase of surge-cycle surge. Glacier Length Area Surge duration Thermal regi km k m 2 years Alaska, Yukon Territory, and British Columbia Variegated, A K 20 49 2 temperate Walsh, A K 89 830 4 ? Muldrow, A K 46 393 2 ? Peters, A K 27 120 2 temperate Tikke, B C 19 75 3 temperate Tyeen, A K 7 11 2 temperate West Fork, A K 41 311 1 temperate Rendu, A K 17 50 2 temperate Hazard, Y T 8 • - 2 polythermal Childs, A K 19 20 1 temperate Unnamed, A K 6 - 2 temperate Black Rapids, A K 45 341 1 temperate Carroll, A K 42 200 1 temperate Iceland Siduj0kull 40 350 1 temperate Dyngjuj0kull - - 2 temperate Tungnarj0kull 25 120 1 temperate Bruarj0kull 50 1500 2 temperate Eyjabakaj0kull 10 50 1 temperate Teigadalsj0kull 1 1 1 temperate Hagafellsj0kull eystri 17 110 1 temperate Hagafellsj0kull vestari 17 100 1 temperate Pamirs and Caucasus Medvezhiy 13 25 1 ? Byrs 2 10 '2 ? Shini-bini 10 16 2 ? Sugran 22 47 3 ? Muzgazy 11 16 2 ? Abramova - - 3 ? Ravak 3 2 1 ? Tanyrnas 10 61 2 ? Burokurmas 7 • 8 2 ? Kolka 6 3 1 ? Chapter 2. Literature review 6 T A B L E 2.1: continued Tien-shan Mushketov 21 70 2 ? Bezymyanny 6 11 2 ? Karakoram Kutiah 12 - 1 ? Baltbare 8 - 2 ? Andes Grande del Junca 10 9 1 ? Grande del Nevado 6 5 1 ? Svalbard Bakaninbreen 17 57 iot polythermal Bodleybreen 15 87 '5+ polythermal Fyrisbreen 3.5 3 8 polythermal Hessbreen 5 6 5 polythermal Hinopenbreen 68 1250 4+ polythermal Hyllingebreen 3 3 10 polythermal Osbornebreen 20 152 3+ polythermal Usherbreen 12 29 8 polythermal T A B L E 2 .1 : Duration of active phase of surge-cycle for various glaciers. Length, area, and surge dura-tion data reproduced from Dowdeswell et al. (1991), and (t) from Murray et al. (1998). Thermal regime information from Clarke (personal communication). 2.1.1 E n v i r o n m e n t a l controls on surge-type glacier d i s t r i b u t i o n About 1% of contemporary glaciers are thought to be of surge type (Jiskoot et al., 1998), the majority of which are concentrated in a few regions. In North America, most surge-type glaciers are situated in the northwestern mountain ranges of Alaska, Yukon Territory and northwestern British Columbia (St. Elias Mountains, Alaska Range, and Wrangell Mountains), and none exist in the southern ranges such as the Coast Mountains, the Cascade Mountains, or the Rockies. Worldwide, surge-type glaciers are found in the Pamirs, the Caucasus, the Tien-Shan, the Karakoram, the Andes, Svalbard, Iceland, and East Greenland (Paterson, 1994). These regions include many different climatic, tectonic, and geologic settings. Within each area, some drainage basins may contain many surging glaciers, while others contain none. This markedly non-random geographical distribution has inspired research on the environmental controls that, distinguish surge-type glacier from their non-surging counterparts (Post, 1969; Clarke et al., 1986; Wilbur, 1988; Clarke, 1991; Hamilton and Dowdeswell, 1996; Jiskoot et al., 2000). Different factors have been investigated. Among them, glacier slope, elevation, orientation, and the presence of tributaries were shown not to have a significant effect on the probability of surging (Clarke et al., 1986; Clarke, 1991; Hamilton and Dowdeswell, 1996). From a population of 146 Chapter 2. Literature review 7 glaciers in the Alaska and the Yukon, Wilbur (1988) found that "bottom-heavy" glaciers, those with concave profiles and a large portion of their area at low elevation, had a higher probability to surge. However, this geometry could well be a result of surging, rather than a cause of it. From a population of 1754 glaciers in the eastern St. Elias Mountains, Clarke (1991) showed that long glaciers are more likely to surge. Hamilton and Dowdeswell (1996) came to the same conclusion for a population of 615 glaciers in Svalbard. Hamilton and Dowdeswell (1996) also detected that glaciers with an internal reflecting horizon in radio-echo sounding profiles have a higher probability to be surge-type. Internal reflections can be caused by a transition from cold to warm ice within the glacier, and a reflecting horizon is often interpreted as the signature of a polythermal glacier. Obviously, statistical analysis alone cannot say whether polythermal glaciers are more prone to surging, or if surging modifies the thermal regime of glaciers. Lithology also affects the probability of surging. Post (1969) noted that surging glaciers are more often associated with fault-shattered valleys, and that none were found in the granitic Coast Mountains. He speculated that permeability of the bedrock might be a factor that influences surging. In Svalbard, the probability of surging is greater for glaciers underlain by sedimentary bedrock (Hamilton and Dowdeswell, 1996; Jiskoot et al., 2000). Clarke et al. (1984) and Fowler et al. (2001) suggested that surge-type behavior occurred over easily eroded material, where subglacial till is abundant and easily regenerated. Alone, none of these factors suffices to predict the distribution of surging glaciers within a region, or to predict the location of surge clusters in the world. Either the conditions for surging depend on an interplay of many of these factors, or the measurements that might reveal the environmental controls are lacking. If subglacial dynamics are responsible for the surges, it is not certain that macroscopic geometrical parameters (e.g. accumulation-ablation area ratio, slope, orientation, etc.) can reflect the dynamical conditions necessary for surging. The paucity of data on the subglacial environment of most glaciers makes it difficult to corroborate the suspected link between surge behavior and lithology. 2.1.2 Surge dynamics : overview of ex is t ing models The occasional periodicity of the phenomenon and the asynchronous surging of neighboring glaciers suggest that surges are triggered by internal instabilities rather than external forces such as earth-quakes or climate (Meier and Post, 1969). Various physical models able to produce flow oscillations have been proposed, with different degrees of support by field evidence. Three main factors have been considered as plausible controls on surge behavior: the thermal regime of the glacier, the na-ture of the subglacial hydrologic system, and the character of the subglacial bed (Sharp, 1988). A brief survey of the main models under these three groups is presented below. Thermal mechanisms Ice temperature has an effect on deformation rate, and thus on creep velocity. However, the most important effect of temperature is to enable or inhibit the presence of liquid water at the bed. A cold-based glacier (with basal temperature below melting point) is strongly coupled to the bed and Chapter 2. Literature review 8 can only flow by deformation; a warmed-based glacier (with basal temperature at the melting point) can be decoupled from its bed by hydromechanical effects and increase its motion by sliding or bed deformation. The critical distinction between these two states has inspired the hypothesis that flow instabilities can be explained in terms of a transition between a cold-based to a warm-based state (Clarke, 1976; Fowler et al., 2001). Clarke (1976) proposed that the timing and initiation of the surge of Trapridge Glacier was the result of the competition between (a) the generation of frictional heat at the bed (which by melting the basal ice enables rapid sliding) and (b) changes in glacier geometry (which have an important influence on the distribution of temperature within the glacier). During quiescence, the progressive thickening of the glacier insulates the base and traps geothermal heat, thus yielding an increase in basal ice temperature. At a critical thickness, the melting point is reached, sliding is enabled and the surge begins. The transport of ice from reservoir to receiving area and the concomitant thinning result in the downward advection of cold ice. The thermal gradient at the base is increased, and heat is conducted upward more easily. O n the other hand, the accelerated motion creates frictional heat. The competition between these two factors controls the duration of the surge, which terminates when the basal temperature returns to below melting point. Fowler et al. (2001) proposed a mathematical representation of a similar idea, where the surge mechanism is explained in term of a positive feedback between the creation of meltwater by enhanced motion and the failure of the subglacial sediments due to an increase in water storage. The argument goes as follow. As the build up of ice brings parts of the bed to the melting temperature, the accumulation of meltwater in the till increases pore-water pressure towards flotation pressure, thus reducing effective pressure and weakening the till. Failure of the till increases basal motion and heat production, thus creating more meltwater, which will raise pore-water pressure further and propagate the failure of the till. The surge is made possible by this cascade effect, and terminates when thinning reduces shear stress below the sediments shear strength. Fowler et al. (2001) suggest that the development and propagation of a wave front, such as the one observed during the surge of Bakaninbreen (Svalbard) and the one observed at Trapridge Glacier (Clarke and Blake, 1991) are associated with these thermal oscillations. Hydrological mechanisms Observations made during the 1982 to 1983 surge of Variegated Glacier provide evidence that surge motion occurred almost entirely by basal sliding made possible by high subglacial water pressure (Kamb et al., 1985). Based on these observations, Kamb (1987) proposed that the surge trigger was a transition of the subglacial drainage system from a normal tunnel configuration to a distributed "linked-cavity" configuration. Rothlisberger (1972) showed that the diameter of a subglacial tunnel depends on the competition between frictional melting of the walls and closure of the cavity by ice deformation. Creep closure rate is governed by the difference between the ice overburden pressure and the water pressure in the tunnel; if an increase in water flux accelerates the melting of the walls, the rate at which the tunnel is closed by creep must increase accordingly for the tunnel to maintain its diameter, and therefore the Chapter 2. Literature review 9 water pressure must drop. Therefore, in a channelized drainage system, the steady-state pressure in the conduit is inversely proportional to the water flux. This inverse relation between water flux and pressure makes large tunnels (low pressure) grow at the expense of smaller ones (high pressure), thus creating a dendritic drainage system. Kamb (1987) describes an alternative mode of drainage, where water is distributed across the bed through a network of linked cavities. These water pockets are formed behind protuberances in the bed (cavitation) and are connected by narrow orifices that constrict the flow. Due to the throttling effect of these orifices, an increase in water flux generates an increase in water pressure. This direct relationship between pressure and flux insures the stability of this distributed drainage system, because bigger pathways do not tend to attract water from smaller ones as in the channelized drainage system. According to Kamb (1987), surges are caused by a temporary transition from tunnel to linked-cavity drainage system. Under these conditions, an increase in water flux at the bed generates an increase in basal pressure, and therefore enhances sliding. In turn, sliding counteracts the enlargement of the connecting orifices (by melting) and stabilizes the linked-cavity system. If the water pressure in the system grows beyond a critical value, the growth and connection of the orifices will reverse the pressure-flow relation, rapidly reestablishing a tunnel-type drainage system, bringing the surge to an end. This is consistent with the release of basal water that commonly accompanies the termination of surges, and with the observation at Variegated Glacier of a peak in water pressure preceding the 1982-1983 surge termination. A similar transition from an efficient discrete basal drainage system to an inefficient distributed one has been observed during the 1991 surge of Skeioararjokull, Iceland (Bjornsson, 1998). Deformable bed mechanisms When a glacier rests on a bed of soft sediment, it is possible for a significant proportion of its motion to be accomplished by deformation of the sediments, rather than by creep or sliding over the bed. Clarke et al. (1984) suggested that this could be a possible explanation of the surges of Trapridge Glacier. Similarly to the model of Kamb (1987), the increase in basal water pressure resulting from disruption of the subglacial hydraulic system is the surge trigger. In the "soft" bed theory however,1 the rapid motion is explained in terms of failure of the bed rather than sliding over the bed. The disruption of the subglacial drainage system causes an increase in basal water pressure, which diminishes the effective stress on the sediment bed thus reducing its shear strength. This local failure of the sediments must be propagated for a surge to occur. Shearing of the subglacial sediments would normally be accompanied by dilatancy and would therefore increase the substrate permeability, allowing for extra storage of basal water, reducing the pore-water pressure and returning the system to an equilibrium. Clarke et al. (1984) suggested, however, that if the drainage in the subglacial till was primarily through pipes within the till or channels at the t i l l -'Thus named by contrast with the "hard" bed theory of Kamb (1987), where the drainage system is carved in the ice, as it would be for glaciers resting on impermeable bedrock. It is important to note that Kamb's (1987) model does not explicitly suppose an impermeable bed; in fact the model is corroborated using observations from the surge of the mainly soft-bedded Variegated Glacier. Its characterization as a "hard-bed" theory mainly comes for the fact that water routing in the substrate is not accounted for. Chapter 2. Literature review 10 ice interface, the deformation of the till could destroy these preferred pathways and reduce the permeability of the substrate. As the permeability of the substrate decreases and basal water pressure increases, the zone of weakened bed extends until a sufficient area of the bed reaches the failure point and rapid motion occurs. The surge terminates when thinning of the glacier reduces the stress to values below the shear strength of the sediments. Substantial layers of subglacial till have been found under most surging glaciers, and for this reason the "hard-bed" theories are losing ground (Harrison and Post, 2003). In both the "hard-bed" model of Kamb (1987) and the "soft-bed" models, the rapid release of stored water is necessary for the surge to occur. Both the switching between hydraulic regimes and the extended weakening of the till require a large input of meltwater at the base of the glacier. Englacially stored water is the most likely source for this release (Lingle and Fatland, 2003). Much still needs to be discovered, however, on the required volumes, mechanisms of storage, and causes of release. 2.2 State of research on Trapridge Glacier Study of Trapridge Glacier began in the late 1960s and has been pursued to this day. Many aspects of the glacier have been examined, with a particular emphasis on basal processes. This section presents an overview of the current state of knowledge of Trapridge Glacier. 2.2.1 1940s surge Scientific expeditions organized by W . A . W o o d from 1935 to 1941 provide what little information we have about the state of the glacier prior to 1951. The expedition geologist R.P. Sharp published his observations about the Wolf Creek Glaciers. "Glacier 13 [Trapridge Glacier], draining from the east slope of Mount Wood, advanced eastward at least 200 feet between 1939 and 1941. As observed from the air in August, 1949, this glacier was advanced far beyond its terminal position of 1941" (Sharp, 1951). World War II interrupted the field work and there is an observational hiatus between the 1941 expedition and a 1951 photogrammetric aerial survey by the Canadian Government. Comparing Figure 2.1 and Figure A.7 shows the dramatic advance that occurred over this period. 2.2.2 T h e r m a l regime and flow evolu t ion 1969-1989 The thermal structure has been studied by insertion of thermistor strings from 1972 to 1989. The 1972 program revealed the thermal map of the bed: a central zone at the melting point surrounded by a ring of cold-based ice. Ice temperature increases with depth, with 15m-depth temperature averaging - 7 ° C (Jarvis and Clarke, 1975). In the early 1970s, the glacier appeared to be in quies-cence, with the receiving area of the 1940s surge gradually being ablated (see air photo, Figure A.4). Surveying of flow poles revealed a zone of strongly emergent flow, resulting from the collision of active and inactive ice near the terminus (Collins, 1972). It was only in 1980, after a hiatus of 6 years, that the morphologic consequences of these two flow regimes became obvious: a wave-like bulge had formed at the boundary between these two zones (Figure A.7)'. A new drilling program was conducted in 1980 and 1981 for the installation of thermistor strings. It established that the "apron" of thin inactive ice downflow from the bulge was cold-based, while the thick active ice upflow from the bulge was warm-based (Clarke et al., 1984). The cold-based apron acted as a mechanical dam and caused the upflow thickening of the glacier. From 1980 to 1989, the bulge advanced at a fairly constant rate, averaging 3 0 m a - 1 , eventually overriding the apron. In 1987 and 1988, new thermistor strings were installed at the same geographical position where thermistors were inserted in 1980-1981 . Comparison with the previous temperature profiles showed that a dramatic change in thermal structure had occurred. Water penetration into surface crevasses had warmed the 15-m ice temperature by ~ 2 ° C . The transition zone between cold- and warm-based ice migrated downglacier, but. at a slower rate than the bulge. Clarke and Blake (1991) argue that the activation and incorporation of the apron was accomplished primarily by creep (see Figures 6 and 7 in Clarke and Blake, 1991). Drilling near the terminus was impeded 15 to 22m above the bed by an internal layer of sediments too thick to be penetrated by the drill. This internal structure was not present in 1980-1981, and is therefore assumed to have been created by the advance of the bulge. It is thought to result from the entrainment of basal sediments by blind thrusting at the warm-bed cold-bed interface (Clarke and Blake, 1991). No further measurements of temperature were taken after 1989, but it is probable that the warm-cold boundary continued to migrate downglacier. 2.2.3 Subglac ia l hydro logy and ice—bed mechanics The glacier bed consists of thick sequences of unlithified sediments, which overlie a flood basalt, underlain by highly fractured low-grade metamorphic carbonates (Sharp, 1943). In the central area of the ablation region, the till has been shown to be actively deforming to a depth of at least 0.3 m, and has an estimated viscosity of 3 x 10 9 - 1.8 x 10 1 0 Pas (Blake et al., 1992). Sedimentary studies Chapter 2. Literature review 12 of the forefield of Trapridge Glacier reveal transmissive layers of sand and gravel of a few meters in thickness with hydraulic conductivity between 10~ 5 and 1 0 _ 3 m s - 1 . These layers are capped by 0.1-1 m of till with conductivities in the range 1 0 ~ 1 2 - 1 0 - 7 m s - 1 (Stone and Clarke, 1993). Water pressure and conductivity are measured year-round beneath the warm-based central area of Trapridge Glacier. During the summer, diurnal oscillations are observed in some pressure records, implying that surface meltwater reaches the bed (e.g. Stone and Clarke, 1996; Flowers and Clarke, 2002). Upon completion of a drilled hole three different behaviors are possible: the hole is uncon-nected to the subglacial drainage system and therefore does not drain, the hole is connected to a low-pressure system and the water level drops, or the hole is connected to a high pressure system and artesian outflow occurs. The state of connection of a hole can also be detected by variations in water pressure. In the summer, connected holes respond to diurnal input of meltwater, and the subglacial water pressure fluctuates. The presence of a diurnal signal is evidence of "vertical connectivity", i.e. connectivity between basal and supraglacial water systems. In the winter, the input of meltwater from the surface is limited and the diurnal signal disappears. "Horizontal con-nectivity", i.e. connectivity of points on the subglacial horizon irrespective of their connection to the supraglacial system, can be determined by searching for events detected by sensors at various locations. In unconnected holes, pressure sensors measure the sediment pore-water pressure, which is generally higher than the pressure in connected holes, approaching the overburden pressure. If anything, the pressure in unconnected holes is anticorrelated with the input of superficial meltwater at the base: a pressure rise in the connected system can lift the glacier, thus reducing the effective pressure over the sediments and the pore water pressure (Murray and Clarke, 1995). The state of connectivity at a given location changes over time and the water pressure in the subglacial environment can change over very short length scales. To examine such variability we installed two arrays of 16 sensors arranged as a 4 x 4 matrix having 2 m spacing. The 16 sensors generally showed a high degree of coherence over the first summer. After a year, however, this coher-ence was lost and each sensor behaved independently, despite their close proximity. This small-scale variability of the hydraulic system makes the synoptic description of subglacial conditions difficult. Point measurements do not represent areally averaged properties of the subglacial environment. The connectivity of boreholes, the pressure sensor records, and the observation of release events lead to the development of a conceptual model for the basal hydrology of Trapridge Glacier (Stone and Clarke, 1993, 1996). Water escapes from the warm-region of the bed via an aquifer passing beneath the marginal permafrost, and emerges in the forefield as springs. This is the main drainage route for both connected and unconnected zones. The input of surface meltwater and generation of basal meltwater creates local "blisters" of water, which form the connected system. As the water pressure in these blisters is generally below that of the unconnected system they must be connected to some form of outlet, possibly preferential drainage routes through the basal sediments. Input of surface water can exceed the capacity of the connected system, causing a release event. These events are rare, observed to occur less than once a year. Before the release, water pressure in the connected blisters can exceed overburden. Pressure drops when an exit route through the cold marginal ice is found. The release event of 22 July 1990, for example, resulted in the expulsion of dyed basal water Chapter 2. Literature review 13 over a period of 2-3 days. During winter, as the input of water diminishes, the areal coverage of the connected system decreases, filled by ice or sediment creep (Stone and Clarke, 1993, 1996). Subglacial water pressure measurements portray a basal drainage system that can be dramati-cally inhomogeneous, both spatially and temporally. Since ice-bed coupling is strongly influenced by basal water pressure, the stresses at the bed will be also heterogenous. Several subglacial instru-ments have been especially designed to study the effect of hydromechanical coupling at Trapridge Glacier. "Slidometers" allow the continuous measurement of basal sliding over a period of a month or so, tilt cells estimate bed deformation, and "ploughmeters" allow the assessment of sediment mechanical properties (Kavanaugh and Clarke, 2001; Fischer and Clarke, 2001). The more recent "bed-stylus", composed of a tilt-cell attached to the tip of a thin and flexible rod of titanium (the stylus), itself mounted on a P V C body for support, gives an estimate of glacier sliding year-round. To study the dynamics of motion at the bed, these instruments are deployed jointly over a small distance (< 20 m), together with water pressure sensors. The first of these arrays, installed in 1992, revealed that basal processes were surprisingly inhomogeneous at such small length scales. Fischer and Clarke (2001) summarized these results and the conceptual models of ice—bed coupling dynamics that they inspired. The bed is seen as a time- and space-varying patchwork of sticky-spots and zones of stick-slip motion; the flow field of ice immediately above the bed is not uniform, bed deformation is not a steady and continuous process, and sliding is not temporally smooth. At places and times of low subglacial water pressure (increased ice-bed coupling) there is resistance to sliding, causing the downglacier shear deformation of the bed and the accumulation of elastic strain in the ice and sediments. In contrast, during rising water-pressure (which decouples the ice from the bed) at other places or times, the glacier slips forward and the bed relaxes elastically in shear upglacier (Fischer and Clarke, 2001). This model can explain the complex and sometimes counterintuitive response of the instruments in the 1992 bed array. However, hydromechanical modeling of the glacier bed as a whole has not been attempted yet. Flowers (2000) created a basin-scale hydrological model capable of quantifying water routing in the complex hydrological setting of Trapridge Glacier. Meteorological data and digital elevation models of the ice surface and the bed are used as input, and quantitative comparison of observed and simulated pressure records allow parameter selection. The model comprises four coupled com-ponents: surface ablation and runoff, englacial water transport and storage, drainage through a subglacial water sheet, and subsurface groundwater flow. Ablation of snow and ice feeds a surface runoff that transfers water to the glacier interior via crevasses or exits the system supraglacially. The englacial system delivers water to the bed, where drainage takes place in a porous sheet composed of water and dilatant sediments capable of porosity adjustment in response to changes in water flux. Beneath this conductive layer a low-permeability till horizon acts as an aquitard. Subsurface flow in a underlying aquifer completes the hydrological system. The Flowers (2000) model allows basin-scale hydrologic description of glaciated environments: volume of water storage, volume and routing of water runoff. While the spatial resolution is too coarse to predict local pressure measure-ments, the general behavior accompanying the development and collapse of drainage modes at the spring and fall transitions can be reproduced (Flowers and Clarke, 2002). 14 Chapter 3 Photogrammetric analysis: 1951-1981 ice surface topography, and bed topography 3.1 Construction of DEMs from aerial photography The extensive characterization of the ice surface topography by ground surveys was not a regular field activity until 1980. To infer the changes in ice distribution that occurred before that time, I rely on aerial photography. Several air photo missions were conducted by the Canadian Government in the St. Elias Mountains, all of them between 1948 and 1981. Five pairs of air photos (1951, 1970, 1972, 1978, 1981) were found to have sufficient overlap and texture to enable stereographic reconstructions of Trapridge Glacier and its surroundings. The orientation of the photographs and the collection of elevation data were performed by a subcontractor. From each stereo pair, two models were created: a fine-meshed model, at 20-m (x, y) resolution, for Trapridge Glacier and a coarser one, at 100-rn (x, y) resolution, extended to cover other glaciers of the Rusty Creek basin. The outlines of glaciers and snow patches were also digitized. Figure 3.1 shows the air photos for 1951, 1972, 1977, and 1981 overlayed by the digitized outlines. Details on the collection of photogrammetric data, including full page reproduction of the photos, are presented in Appendix A . From geological evidence and direct observations, we know that most glaciers in this region are surge-type. In that, regard, certain features of these photographs are worthy of note: 1. In 1951, Trapridge Glacier is in a post-surge stage: highly-crevassed and extended far beyond its 1941 limit (compare to Figure 2.1). 2. Rusty Glacier surged prior to 1951. The ice of its receiving area was not completely melted in 1951, but was mostly gone by 1970. Since that time, the glacier has continued to retreat. 3. Backe Glacier surged sometime after 1951 and before 1970, overriding the down-wasting ter-minus of Rusty Glacier. Between 1970 and 1981 there is little evidence for retreat of its terminus. 4. It is difficult to judge whether the apparent increase in the number of crevasses at the surface of the Hodgson Glacier in 1972 is due to an increase in activity or merely due to differences in illumination. Inspection of the D E M s show that there was a significant variation in elevation over the course of these years (see Figures 3.2 and A.9). Chapter 3. Photogrammetric analysis: 1951-1981 ice surface topography, and bed topography 15 FIGURE 3 . 1 : A i r photos of Rusty Creek basin, 1951-1981. Digit ized glacier outlines are superposed. O n this image, the outline does not always match the glacier contour because these photographs were not corrected for parallax, lens distortion, and film deformation. These corrections were made before the outlines were digitized, but the corrected photos could not, be outputted due to software limitations (see Appendix A for more details). Unfortunately, the post-surge receiving area of Rusty Glacier was not included in the 1951 outline. The debris-covered tongue was probably mistaken for the ground surface by the operator. The southwest corner of these frames, and of all subsequent maps of the Rusty Creek basin, is located at 533000 m East, 6874000 m North ( U T M coordinates). Chapter 3. Photogrammetric analysis: 1951-1981 ice surface topography, and bed topography 16 500m N F I G U R E 3 . 2 : (a) Mean and standard deviation of the 1951, 1970, 1972, 1977, and 1981 photogrammetric D E M s . The contours mark the mean elevation in increments of 20 m and the color the standard deviation. L imi t s of the 1951 and 1970 photographs, which do not cover the whole area, are marked by a dashed line, (b) Same, except that the 1951 D E M was excluded from the population. Note the dramatic change in variance. 3.1.1 Tes t ing of the D E M s and evaluat ion of uncer ta in ty Relative uncertainty I test these D E M s in two different ways. First, to insure that the models are consistent, I calculate the mean Z{x, y) = i £ " = 1 zi(^, y) and variance d2{x, y) = YA=\ (zi(x, y) ~ Z(x, y))2 of the elevation of each cell in the five models (n = 5). The standard deviation field (Figure 3.2) allows me to verify that stable elements of the landscape align from one year to the next. During the early stages of the photogrammetric analysis, this test identified models that were poorly registered and thus required revision. Figure 3.2 shows the results of these local statistics for all years (left side) and all years except 1951 (right side). Several conclusions concerning the quality of the D E M s can be reached: 1. The difference between Figures 3.2a and 3.2b indicates that much of the variance in Figure 3.2a is due to the inaccuracy of the 1951 model. As explained in Appendix A , the collection of points for the 1951 model was made difficult by the large overlap (small parallax) of the photo pair. I therefore consider 1951 to be a special case, distinct from the remaining four D E M s . 2. As seen in Figure 3.2b, mountains and non-glaciated zones down valley do not change sig-nificantly among these four models, with an average standard deviation of ~ 1 0 m . This is a reassuring result. The standard deviation of these unchanging features can be considered as representative of the relative uncertainty of the D E M s . This criterion cannot detect system-atic errors in the registration. Because absolute positioning is of no interest for this work, this is the uncertainty value I attributed to these D E M s . 3. As expected, glaciated areas have higher variances. These are elements of the landscape Chapter 3. Photogrammetric analysis: 1951-1981 ice surface topography, and bed topography 17 that have changed significantly over the 30 year photographic record. For example, note in Figure 3.2a the effect of the surges of Trapridge Glacier and Backe Glacier. 4. Not all regions of high variance can be attributed to glacier activity. The large elevation deviation in the upper basin of Backe Glacier (Figure 3.2b) is due to anomalous values in the 1970 D E M . The absence of ground control points in the vicinity makes this zone vulnerable to distortion. Errors in the estimated elevation can also be explained by the lack of contrast in this area, which makes it difficult to evaluate the position of the ground surface. This anomaly does not affect the zone around Trapridge Glacier, which is the area of interest. A b s o l u t e u n c e r t a i n t y The second test compares the 20-m (x, y) resolution D E M s for 1972 and 1981 to field measurements taken over these two years. This allows me to estimate the accuracy of the 20-m resolution D E M s , and to verify that the two data sets are based on the same datum. Figure 3.3 plots the difference between the survey measurements and the (interpolated) D E M values. In 1972, the D E M elevations are on average 5 m below the surveyed elevations. The D E M has a slight positive bias in 1981, but it underestimates the height of the bulge. In both years, the D E M tends to underestimate the upper basin elevation. Because of the high reflectance of the snow, the photos lack contrast in this area. This makes it difficult to determine precisely the location of the ground surface when viewing the photo in stereo. The difference between most other poles and the D E M s are within the 10 m uncertainty bounds. 3.1.2 Resulting DEMs Figure 3.4 shows the 1981 coarse and fine D E M s draped by the original registered air photos. As concluded in the previous section, the elevation uncertainty for the 1970, 1972, 1977, and 1981 models is roughly ± 1 0 m. I estimate the 1951 elevation uncertainty to be ± 1 5 m. These bounds suffice to explain the observed changes in topographical features that should have remained stable over such a short time period. By themselves, these D E M s are of limited use to characterize the surge. The digitized outlines describe the retreat and advance of the glacier. A n absolute year-to-year volume change could be evaluated by assigning an arbitrary reference bed, but the relative year-to-year volume change requires knowledge of the actual glacier volume, and hence of the bed topography. Ice thickness information is also necessary to evaluate basal stress and hydrological potential. Serendipitously, the last series of air photos of Trapridge Glacier (1981) were taken when it was close to its position of maximal retreat. Much of the ground that was once covered with ice was then exposed. Therefore, the 1981 D E M gives not only the ice surface topography for that year, but also a significant portion of the bed topography for other years. Radar soundings taken during the 1994-1997 field seasons complete this portrait. Chapter 3. Photogrammetric analysis: 19511981 ice surface topography, and bed topography 18 J • 1 1 ' I L ^^^••lllllllllllllllllllllllllll I I 4 12 Elevation difference (m) F I G U R E and 1981 3 . 3 : Comparison of ground measurements and interpolated photogrammetric D E M elevation, 1972 . Circles mark the location of flow poles. Their color gives their elevation relative to the D E M . Chapter 3. Photogrammetric analysis: 1951-1981 ice surface topography, and bed topography 19 F I G U R E 3 .4 : 1981 DEMs, coarse (100-m resolution) and fine (20-m resolution), draped with the original (right-side) air photo. The vertical exaggeration is x4, and the axis ticks for the coarse model are separated by 500 m. 3.2 Bed topography Before attempting to interpolate a bed topography map from these two data sets, it is reasonable to question whether one can suppose the bed topography to be constant in time. By splicing radar data from the mid-1990s to a 1981 D E M , and using the resulting bed map to calculate the 1951 ice thickness, am I piecing together parts of different puzzles? Can I assume the 1981 exposed ground surface to be the same as what in was in 1951, when it was covered by ice? Two questions need to be addressed: (1) does the topography of the bed change significantly when it is covered by ice? and (2) how quickly does the bed erode when it is uncovered? The main processes by which sediments are moved under a glacier are entrainment by subglacial flow, and sediment advection by ploughing or regelation. During the quiescent phase of Trapridge Glacier, it is unlikely that either of these could move a significant mass of sediment. Subglacial water flow at the bed is too weak, both in discharge and velocity, to carry a significant sediment Chapter 3. Photogrammetric analysis: 1951-1981 ice surface topography, and bed topography 20 load. W i t h the exception of the rare release of englacially stored water, most of the basal water escapes by ground-water flow through an aquifer (Stone and Clarke, 1996). Transport of sediment by ice advection is limited by ice flow velocities. Given the low velocities observed, I do not expect this to significantly alter the bed topography over the 30 years studied. The strongest evidence for a relatively stable bed is the absence of sediment moraines in the proglacial field, which should be observed if there was significant evacuation of bed sediments. The rate of erosion in the proglacial area between 1951 and 1981 is harder to estimate. If erosion was significant, by using the 1981 ground surface as a proxy for the 1951 bed, I will over-estimate the 1951 ice thickness. Figure 3.5 maps the elevation difference between subsequent D E M s . As one can see, it is a dynamical landscape, and both erosion and accumulation occur. The area below the 1951 terminus is ~ 5 0 m above what it was in 1981; inspection of the 1951 photograph reveals that this is due to ice remaining from a previous surge of Rusty Glacier (Figure A.3). The 1977 — 1972 difference shows that there is a bias in the absolute positioning of the 1977 model. I therefore used the 1972 model to compare with 1981. The hot spot on the bottom-left is probably due to the accumulation of sediments transported by Rusty Creek, which flows nearby. Surprisingly, the ground level seem to have risen in the proglacial area of Trapridge Glacier between 1972 and 1981. Could it also be the result of the accumulation of sediments transported by the proglacial streams? It is difficult to ascertain. The lack of precision of the D E M s is probably responsible for part of that variability. Some field evidence suggests a more stable ground surface. In the 1970s and 1980s, several years after the retreat of the ice, delicate subglacial features could still be observed in the proglacial area (Clarke, personal communication). The preservation of these features indicates that erosion or accumulation were limited. If there was erosion in the lower reaches of the valley, it was most likely along proglacial streams. The stream channels visible on the 1981 were not necessarily present under the 1951 ice. However, these changes are localized, and will not greatly affect the 1951 ice thickness estimate. Face with contradicting evidence, I come to no real conclusion. By lack of a better option, I use the 1981 model to infer the 1951 bed. If the D E M s have correctly captured a rise of the ground level between 1972 and 1981, this will result in an underestimation of the ice thickness. The parts of the 1981 D E M showing the ground surface and the radar data are used jointly to obtain a bed topography map. Because these two data sets are sampled irregularly, I use ordinary kriging to interpolate the data on a regular 20-m grid. To avoid interrupting the current discussion by technicalities, the details of the analysis are presented in Chapter 5, after after an introduction to the theory of kriging. The result of this interpolation is used to infer ice thickness maps from the D E M s . Chapter 3. Photogrammetric analysis: 1951 1981 ice surface topography, and bed topography 21 F I G U R E 3.5: Elevation difference between subsequent photogrammetric DEMs. F I G U R E 3.6: Interpolated bed topography. The contour lines give the bed elevation, and the color gives the detrended bed topography. Notice the two pinning points (elevation of 2250 m) that are responsible for the characteristic shape of the terminus (see Figure 7.1). The details of the geostatistical analysis are presented in Section 5.5. Chapter 3. Photogrammetric analysis: 1951-1981 ice surface topography, and bed topography 22 3.3 Results: Ice thickness maps, outlines, cross-sections The ice thickness field is obtained by subtracting the 20-m gridded bed map from the 20-m resolu-tion photogrammetric D E M s (Figure 3.7). The glacier outlines are defined by visual inspection of the air photos and digitized separately from the D E M s . Their concordance with the ice thickness map (Figure 3.7), especially in the thin apron region, inspires confidence in the capacity of the photogrammetric analysis to resolve elevation differences with precision. The period 1951—1981 is characterized by the melting of the ice accumulated in the receiving area during the 1940s surge. Between 1971 and 1981, the area covered by left-over ice (the "apron") is rapidly decreasing. Up-glacier, a thicker zone of active ice is advancing. Figure 3.9a shows the evolution of the glacier outline. In its 1951 post-surge state, the glacier has a bottom-heavy profile, with a ~90-m thick terminus. Between 1970 and 1981, the obstruction of the flow leads to a thickening of the ice up-stream from the apron. By 1981, a ~ 4 0 m ice front is formed at the contact zone between active and stagnant ice. From these D E M s and the digitized outlines, the area and volume of the glacier can be estimated (Figure 3.8). Because the limits of the upper basin cannot be clearly established from the D E M s , I limit for now the discussion to changes that occurred in the lower basin. A n estimate of the upper basin area is established in Section 4.3 from survey data, and the variation in upper basin volume is calculated for all years in Chapter 6. The apparent increase in volume between 1972 and 1977 is mainly due to the positive bias of the 1977 model, as illustrated in Figure 3.5. Because no ground survey measurements were taken that year, I cannot correct the elevation bias of this D E M . I will disregard this model for the rest of my analysis. Figure 3.8a indicates that the volume evolution in the lower basin is non-monotonic. Two opposed processes affect the lower basin mass balance: the melting of the receiving area, and the transfer of ice from the upper basin. The near disappearance of the apron between 1970 and 1981 is responsible for a rapid decrease in glacier area. This dramatic retreat, however, did not change significantly the total ice volume, because the apron is on average only ~10-m thick. Its ablation yielded only a slight decrease in volume, barely detectable between 1970 and 1972. After 1981, the melting of the apron and the advance of the active ice have continued. Whether these two inverse tendencies have affected significantly the total volume of ice is unknown. To investigate this further, D E M s for the subsequent years are necessary. Aerial photography missions in the St. Elias Mountains stopped in 1981, and no further information can be gleaned from this source. After 1980, Trapridge Glacier was visited annually, and survey measurements of flow poles and of the ice surface were taken. The coverage of the survey data increased sufficiently to give a fair idea of the geometrical changes that were happening. In Chapter 5, I present a Bayesian algorithm for kriging the ice thickness field that takes advantage of the continuity of this data set. Before this can be achieved, however, the data must be filtered to eliminate the inescapable measurement blunders and data entry errors. The following chapter explains how ice flow calculation can be used to this effect. Once the data are refined, they are used to quantify the changes in ice flow that accompanied the current surge. Chapter 3. Photogrammetric analysis: 1951-1981 ice surface topography, and bed topography 23 0.16 220 0.13 200 on ^ .5 1 8 0 | Z 160 0 0 8 S j u o 120 100 70 72 0.05 0.03 70 72 F I G U R E 3.8: Total ice volume (a ) and glacier area ( b ) , relative to the 1981 values (left axis) and in absolute values (right axis). Chapter 3. Photogrammetric analysis: 19511981 ice surface topography, and bed topography 24 F I G U R E 3.9: (a) Evolution of the glacier outline from 1951 to 1981. The ice accumulated in the receiving area during the 1940s surge is being ablated. In 1970, only about ~10m of stagnant ice remain. Given its thinness, the apron rapidly shrinks from 1970 to 1981. (b) Longitudinal D E M cross-sections. The width of the line represents the uncertainty of the surface elevation, The 1951 profile is shown as a reference (dashed line). The dotted line marks the limit of the lower and upper basin. The lower 1.5 km of the lower basin is melting away, while the upper 1 km is thickening. The upper basin has not changed significantly over this period. 25 Chapter 4 Ice flow calculation and estimation of glacier outlines 4.1 Description of the data set Two types of ice topography measurements were made over the 36 years of field work at Trapridge Glacier. First and foremost is the surveying of flow poles drilled into the glacier, which allows local estimation of the speed and position of the ice surface. Second is the surveying of free moving targets on the glacier surface (line measurements). During the survey of a "line", a person walks a prism along the glacier, stopping regularly for the surveyor to take a reading. These measurements are not attached to a specific marker on the ice surface, and therefore cannot be used to calculate ice velocity. Their main purpose is to characterize the ice surface, and to detect changes in geometry. The flow poles used at Trapridge are 3 m (10 feet) in length and generally drilled ~ 2 . 5 m into the ice. Poles are tracers, marking the trajectory of a particle of ice located at their base. Their lifespan depends on the ablation rate; the average melt out time for a pole is 3.6 years. A prism is fixed to the top of every pole and its location measured daily throughout the four-week summer field season. The average speed of the ice particle, in the horizontal plane as well as in the vertical, can be estimated from repeated measurements. The position of the ice surface, distinct from the position of the ice particle the pole is tracking, can also be measured relative to the top of the pole. This is done at the beginning and the end of every field season, when the prisms are installed and removed. Figure 4.1 shows the pole network as it stood at the end of the 2004 season. The extent of the flow pole array has varied significantly over the years. Hot summers and low-precipitation winters increase the net ablation and leave the ice surface covered with fallen poles that must be redrilled. Years of high accumulation cause some poles to disappear under the snow; these are temporarily lost and cannot be surveyed. When a pole has fallen, the position of the ice particle it was following is lost and the time series ends. If a pole is redrilled before it has fallen, or shortly after, its new position is known with respect to its old one and the continuity of the time series can be preserved. The presence of such "synchronous resets" (by my definition, any reset occurring within three days of the previous measurement) in the middle of a time series does not prevent the calculation of displacement; care must simply be taken to subtract the artificial movement caused by the pole reallocation. The flow that occurred between the time of the previous measurement and the resetting of the pole is neglected. Chapter 4. Ice flow calculation and estimation of glacier outlines 26 F I G U R E 4 . 1 : Extent of the pole network al the end of the 2004 season. This is the maximal number of poles we had on the glacier. To analyze changes in flow, the glacier was divided into five zones: four in the lower basin (A, B, C, D) and two in the upper basin (UB1, UB2). 4.1.1 Measu remen t accuracy O n several occasions during the 2005 field season, the flow pole array was surveyed twice in the same day. The difference between these repeated readings gives an estimate of the measurement error (Figure 4.2 a). Measurements taken on successive days can be added to increase the size of the population. If the average glacier velocity is ~ 2 0 m a _ 1 , the displacement expected from one day to the next is roughly 5 cm. If one neglects the flow that occurred between the two measurements, their difference gives another estimate of the measurement error. Because the daily surveying of the pole array is standard practice at Trapridge, including these measurements in the analysis increases dramatically the size of the population. Figure 4.2b shows the distribution of daily-displacements in the easting, northing and elevation, for 3316 pairs of measurements taken in the central area of the glacier (zone B in Figure 4.1). If the measurements were perfectly accurate, the mean displacement would indicate the average daily flow, and the displacement variance would reflect the spatial and temporal variability of the flow. The magnitude of the mean displacement ( 5 . 4 c m d _ 1 , or 1 9 . 6 m a - 1 ) is consistent with the velocity history of the glacier, and its orientation respects the general E N E direction of the flow. The variance of this population, however, is larger than that of the 2005 repeated measurements. Some of this difference can be explained by natural flow variability. These displacements were calculated using the 30 years of data available, including periods when the glacier was much more active than it was in 2005. Subglacial events causing the release of stress accumulated at the bed have been shown to cause rapid glacier motion (Fischer and Clarke, 2001). However, extreme displacement always occurs in the direction of the flow; if natural variability was the only source of variance, we would expect the dx distribution to be highly skewed towards positive values. Its symmetry suggests that Chapter 4. Ice flow calculation and estimation of glacier outlines 27 dx (m) dy(m) dz(m) FIGURE 4.2: Difference in easting (dx), northing (dy), and elevation (dz) between repeated measurements of flow poles, (a) To minimize displacement due to flow, only pairs of measurements taken on the same day are considered (N = 140). (b) To increase the size of the population, measurements taken one day apart are included (TV = 3316; note the logarithmic scale). More outliers are visible, and the standard deviation in the direction of flow (x) is larger than across the flow (y). The standard deviation in the vertical is larger because of numerous outliers. Most of them are due to negative displacements, and could be caused by the tipping of poles or the sliding of prism mounts. The data used for this analysis have been previously refined to eliminate outliers. This suggests that the sieving method presented in the next section is not 100% efficient. measurement errors are responsible for a significant portion of the variance in that direction. For this population, the variance in the vertical is much larger than in the other directions. The extreme values in vertical displacement are probably caused by the tipping of poles or the sliding of the prism mount, two common occurrences. The fact that most extreme values are negative corroborate this hypothesis. While the data used for this analysis were previously treated to expunge blunders (in the way explained below), some might remain. Blunders contribute extreme values to the distributions, irrespective of the instrument precision. Based on these two populations, and bearing in mind the contribution of blunders and natural variability, I estimate the survey measurements to be accurate to ~10cm. Chapter 4. Ice flow calculation and estimation of glacier outlines 28 4.1.2 Expunging of blunders Errors in readings or in data entry are inescapable; the survey site is a cold and lonely place where fingers get numb easily and where laptop screens are hard to read. In early years, the surveyor was accompanied by a notetaker, but this proved to be an expensive luxury. Figure 4.3a shows the superposition of flow vectors calculated by differentiating the first and last measurement of each pole for each year. The effect of blunders is obvious. It is necessary to exclude these outliers before the data can be used to derive flow maps, kriging being an exact interpolator, it does not smooth out spikes due to outliers; elimination of blunders is therefore also necessary for a correct estimation of the ice surface topography. The filtering was done in two steps. First, obvious outliers were eliminated by visual inspection. The criteria for the flagging of blunders depend on the type of measurement. Successive points on a profile line are named sequentially and should follow a walking path. The trajectory of flow poles should be monotonic and aligned with the general flow direction. Inconsistencies in these sequences are attributed to blunders. The 328 flow pole time series, and 211 lines were individually inspected and outliers weeded out based on these criteria. However, as Figure 4.3b shows, this visual inspection was not sufficient to eliminate all blunders. The second filtering step is appropriate for flow pole measurements only. The basic idea is simple: ice velocity vectors are calculated from the data, and are judged as being physically plausible or the result of a measurement error. However, this is not as straightforward as it may seem. Two measurements are needed to calculate a displacement; if the result is deemed to be unreasonable, one or both measurements could be responsible. The probability of two successive blunders being very low, the task is to determine which of the two points must be rejected. Eliminating both would be simple but wasteful; despite a certain degree of redundancy (poles are surveyed several times over a field season), the elimination of a measurement can come with a significant cost if it breaks the continuity between a year and the next. This would occur for example if the measurement prior to a resetting of the pole was deleted. The first years of the data set do not offer much redundancy, and the rejection of a measurement generally means the loss of all information available on that pole for that year, and therefore the rupture of the time series. I designed an algorithm to evaluate which of the two end-points of the trajectory is responsible for the erratic displacement. It is briefly outlined here, and details can be found in Appendix B. First, the time series of each pole is segmented into continuous sequences (i.e. uninterrupted by the melting of the pole). Each measurement within that sequence, is tested by comparing it to other measurements within that sequence. Velocity vectors are calculated using the tested measurement as one of the end-points. The plausibility of each of these velocity vectors is judged, and a weight is assigned based on the quality of the second end-point. Measurements that have been already tested are trusted, while those that are separated by more than one year, or less than ten days are penalized. (The averaging associated with a long time separation diminishes the effect of blunders and reduces the power of the test; for short time separation, the displacement is comparable to the accuracy of the measurement and the calculated velocity is uninformative.) Finally, the weighted Chapter 4. Ice Bow calculation and estimation of glacier outlines 29 I 1 500 m F I G U R E 4.3: To evaluate the method used to detect blunders, it is useful to consider the effect of each step on the calculation of annually-averaged flow vectors, (a) Flow vectors (for all years) obtained using raw data. (Arrow heads are not drawn to avoid cluttering further the figure.) Total number of measure-ments: 10609: number of annually-averaged flow vectors: 1416. (b) Flow vectors calculated from handpicked data. 140 measurements were discarded, resulting in the loss of 18 flow vectors, (c) Flow vectors calculated from numerically-filtered handpicked data. 251 more measurements were discarded by the automated rou-tine, and another 41 flow vectors were lost. In total, 391 measurements were rejected (3.7% of total), causing the loss of 59 flow vectors (4.2% of total). No more outliers are visible. A Af Chapter 4. Ice flow calculation and estimation of glacier outlines 30 results of the velocity tests are summed, and a decision is reached on the quality of the measurement based on whether it yielded mainly plausible flow vectors, or mainly implausible ones. The effect of these two weeding steps on the calculation of annually-averaged flow vectors is shown in Figure 4.3c. No more outliers are visible, and less than 5% of the measurements were lost. 4.2 Results: Changes in flow velocity Given the frequent melting of poles and their tendency to disappear and reappear depending on snow conditions, the flow information extracted from these data varies greatly in coverage from year to year. This variability makes it difficult to represent changes in the structure and the intensity of the flow in a consistent manner. Sacrificing spatial resolution for temporal continuity, I separated the glacier into six zones (see Figure 4.1) and considered the flow vectors in each zone as an ensemble. Figure 4.4 shows the median value and spatial variability in these zones for all years. For more spatial resolution, one should refer to the complete set of flow maps presented in Appendix B (Figures B.1-B.3). Figure 4.4 highlights these features of the recent flow history of the glacier: 1. From 1969 to ca. 1991, the glacier is divided into two regions with very distinct flow rate. Upglacier the ice is flowing, while downglacier the ice is stagnant. As discussed earlier, the in-teraction between these two flow modes is responsible for many of the significant morphological changes that occurred during the last 30 years. 2. Sometime between 1974 and 1981, the glacier started a new surge. Before 1974, the speed in the region of active ice (zones UB2, A , and B) was ~ 1 0 m a _ 1 . In 1981, the median velocity in the lower reaches of the upper basin (UB2) reached 1 9 . 8 m a - 1 . Immediately downglacier, in zone A , the ice was flowing at 3 8 . 6 m a - 1 ; a 4-fold increase from its 1974 value. Zone B took three years to reach the same velocity, but thereafter behaved synchronously with zone A . In 1974, zone C was fully comprised of stagnant ice; in 1981, some of the poles in that region were flowing at more than 25 m a - 1 . The coexistence of the two flow modes explains the high variability observed in zone C between 1981 and 1986. Sometime between 1991 and 1995, the bulge reached zone D, and the median ice flow in that region reached ~ 2 0 m a _ 1 . 3. The higher reaches of the upper basin (UB1) do not seem to have reacted to the surge, as the median velocity remained around 1 0 m a - 1 . The apparent slowdown between 1981 and 1987 is a consequence of the lack coverage in that region. During these years, only poles close to the ice margin where surveyed, and the rate closer to the flow axis is unknown. 4. The lower part of the upper basin (UB2) is more active. It nearly doubled velocity between 1974 and 1981 (from 11.1 to 19 .8ma _ 1 ) . From 1994 to 2001, median ice speed remained relatively constant at around 17 m a - 1 , then started to decrease. In 2005, the flow was reduced to 1 0 m a - 1 — below the pre-surge values. 5. From 1985 to 1997, the flow in the lower basin is on a slowing-down trend. The slowdown followed strangely regular 4-year pulse: 1-2 years of timid acceleration ( 2 - 3 m a - 1 ) , followed Chapter 4. Ice flow calculation and estimation of glacier outlines 31 Zone Activation Surge Waning Now J969-1972 1980-1999 1999-2005 2005 mean entry peak exit mean UB1 5.4 0.3 13.2 (1990) 11.5 11.4 10.7 UB2 11.4 19.8 23.0 (1981) 19.4 14.9 9.9 A 14.4 38.6 41.8 (1983) 28.5 17.2 8.4 B 7.6 26.8 35.1 (1983) 25.7 16.2 7.6 C 0.2 5.4 33.8 (1990) 23.9 16.5 7.6 D 0.3 0.0 31.7 (1993) 24.9 24.9 T A B L E 4 . 1 : Summary of flow history, 1969-2005. Based on the flow record, I divide the last 31 years in three periods: activation, surge, and waning. The median velocity in each zones is averaged over the activation and waning phase. The median velocities at the beginning (1980) and end (1999) of the surge are given, with the peak velocity attained in each zone (year specified in parenthesis). In 2005, there was no flow pole in zone D. The zones are defined in Figure 4.1. by 2-3 years of rapid deceleration ( 4 - 8 m a - 1 ) . The 1997-1999 acceleration was particularly vigorous, as the velocity in zone A went from 2 0 . 3 m a - 1 to 2 8 . 5 m a - 1 , an increase of 8 . 2 m a - 1 . After this last pulse, the glacier gradually slowed down to pre-surge velocities. In 2005, the lower basin was flowing at less than 8.5 m a - 1 . Based on this flow history, I divide the 1969-2005 period in three parts: the activation phase (1969-1974), the surge (1980-1999), and the waning phase (1999-2005). Table 4.1 summarizes the flow characteristics in each zone during these three periods. Chapter 4. Ice flow calculation and estimation of glacier outlines 32 0 l — l l I l I ' l I i T . . I -r Q -r- J. J l . l I l l . l I i 1 1 1 1 1970 1975 1980 1985 1990 1995 2000 2005 year F I G U R E 4.4: Zonally-averaged ice flow. The zones are defined in Figure 4.1. Zone U B 1 is the furthest upglacier, and D is the furthest downglacier. For each year, the box indicates the upper and lower quaxtiles of the speed distribution. The whiskers are set to 1.5 times the interquartile range, and any value outside this range is plotted as a circle. The solid lines follow the median values. Chapter 4. Ice Bow calculation and estimation of glacier outlines 33 4.3 Estimation of glacier outlines Information on the glacier limits is provided by air photos for the few years where they are available (1970, 1972, 1978, and 1981; see Section 3.3). In 2004, a walk around the glacier with a handheld G P S provided an outline with a ± 10 m uncertainty. For all other years, the glacier limits must be estimated from the survey measurements available. Starting ca. 1993, the end of the profile line surveyed along the glacier centerline ends at the glacier terminus, marking explicitly the position of the snout. Before that year, it is unclear if the center profile stops at the limit of the active ice front, at the tip of the apron, or anywhere between these two. Wi th the exception of a few measurements in 1996, no other survey points were taken to mark explicitly the limits of the ice. Measurements taken on the ice give a minimum value for the glacier extent, but simply defining the outline by connecting the most outward measurements would highly underestimate the glacier area. Instead, I draw a "most likely" outline for each year based on the known position of the snout and on other points surveyed near the glacier boundary, assuming a constant progression. The 1970 photogrammetric outline starts this time series, and the 2004 G P S outline terminates it, giving me two known end-points. The tracing of the outlines is somewhat subjective, but a human mind is often the most efficient way to include soft information in an inference. In any case, an algorithm to infer the glacier outline numerically would have had to be tuned until it gave "satisfactory" results, the definition of which would be, in the end, based on the same subjective analysis. To quantify the uncertainty, I bounded the nominal outline by minimal and maximal limits. These bounds are traced for each year, based on hard evidence (measurements made on the ice in the adjacent years, position of the snout, reference outlines marking the definite minimal (1981) and maximal (2004) limits) and on a subjective appreciation of what rate of advance (or retreat) is plausible. Some of the resulting outlines, with their uncertainty bounds, are presented in Figure 4.6; the complete set can be found in Appendix B. The most contentious areas are the southern portion of the terminus and the apron. The southern terminus region is highly crevassed and cannot be accessed. In addition, a significant portion of it cannot be seen from the survey site. The 2004 G P S outline is the only information available on this front. The melt rate of the inactive ice left over from the previous surge is also impossible to estimate. As it is not a zone of great dynamical interest, it was not surveyed often. Flow poles indicate that ice is present, but the extent of the apron is impossible to determine. Even aerial photography gives ambiguous results, as distinguishing the debris-covered apron and the ground can be difficult. Sometime between 1974 and 1981, a portion of the apron was disconnected from the main body of the glacier. I did not include it in the nominal outline, but tried to simulate, very subjectively, its disappearance when tracing the maximal outlines (see Figure B.6). The southwestern end of the glacier is bound by snow-covered mountains and visual inspection of air photo cannot establish where the glacier ends and when the mountain begins. In 1995 the upper basin was extensively surveyed to locate radar soundings; I used these points to estimate the glacier extent. A rapid break in slope in the photogrammetric D E M at the outer edge of these data marks the most likely boundary between the mountain and the glacier. The southern and northern limits of the glacier Chapter 4. Ice flow calculation and estimation of glacier outlines 34 were the same in 1981 and 2004, and are unlikely to have changed in-between. Figure 4.5 summarizes the change in glacier area that occurred since 1951, as captured by these outlines. Chapter 4. Ice Row calculation and estimation of glacier outlines 35 F I G U R E 4.6: Estimated glacier outlines: 1974, 1981, 1991, 1996. 1999, and 2004. The shading gives the uncertainty of the outline: dots are flow poles and +'s are profile line measurements. The 1981 outline was digitized from an air photo (Figure A.7). A significant portion of the apron got disconnected from the glacier between 1974 and 1981; it was not included in the nominal 1981 outline, but the maximal outline was traced to account for its presence. In 1996, the northern end of the terminus was surveyed, providing a uniquely explicit source of information. The southern end of the terminus was not surveyed because most of it is not visible from the survey site. The 2004 GPS outline gives the final limit of the glacier advance. The snout, position has barely changed in the last five years. The 10 m uncertainty bound around the 2004 GPS outline is not distinguishable from the line width at this scale. 36 Chapter 5 Kriging This chapter defines the geostatistical tools used in Chapters 3 and 6. The first section is a heuristical presentation of geostatistics and the interpolation problem. Ordinary kriging, a best linear unbiased estimator, is then introduced. Before any stochastic estimation is possible, a statistical model must be assumed. The experimental variogram (Section 5.3) provides a means to parameterize the covariance function, and the analysis of orthonormal residuals (Section 5.4) a way to test the model. In Section 5.5, these concepts are applied to obtain a gridded bed topography map. Finally, a Bayesian extension of ordinary kriging is presented. This tool is used in Chapter 6 to generate a time-evolving D E M of the glacier. 5.1 Introduction: geostatistics and the interpolation problem Given a field Z(x) sampled at a finite number of sites {x\,... , X J V } 1 and some other prior infor-mation, how can one estimate the value of the field at an unsampled location XQ, and quantify the estimation error? This is known as the interpolation problem. The geostatistical approach considers the data to be the partial realization of a random field. The distribution of this random field is constrained by simplifying assumptions, and an estimate can be obtained by averaging the ensemble of realizations that obey the data. Geostatistics, like time series analysis, must depart from classical statistics to address the fact that the data used are autocorrelated. The classical assumptions of random sampling and inde-pendence of the variates does not hold in spatial analysis; geological variables have deterministic sampling locations that are especially selected to insure that results are not independent from one another, but rather expressions of a common geological object. Whereas classical statistics would normally assign one random variable per attribute for all sites, geostatistics considers one random variable per attribute per site. Because the geographic location space is continuous, the number of variates is essentially infinite and the regionalized variable is represented by a random function, or random field. Let T> be the spatial domain in a n-dimensional euclidian space and let Cl be the permissible domain of variation of the attribute, i.e. the sample space. The random function representing the regionalized variable is the collection of random variables {Z(X,LJ) : x G T>,LJ G Cl}. Z(xi,u>), the value of the attribute at the location x^ is a random variable dependent on the parameter LJ. Given a specific LJJ G Cl, Z(x, LJJ) is a deterministic function of the spatial coordinates, a realization of the 1 Through the text x symbolizes the general position vector, in two or three dimensions. Because the components of the vectors are never used separately, I preferred the simpler notation x over the more explicit x. Chapter 5. Kriging 37 random function. Z(xi,Uj) is the value of the attribute at the location Xi for that specific realization, a number. For convenience, it is customary to refer to {Z(X,LJ) : x G V, u e fl} as simply Z(x) and to denote the observed realization Z(X,OJO) by the uncapitalized z(x). Geostatistical analysis can be subdivided into two steps: structural analysis and best linear un-biased estimation. The goal of structural analysis is to limit the population of plausible realizations. When asked to guess the value of a given attribute at a site, a practitioner will intuitively limit the scope of his estimation by steps: the value must be physically plausible, it must be within the observed range of variation, it should be similar to results obtained at neighboring sites, or if no nearby data are available, it can be assumed to follow regional trends. When asked to qualify the reliability of his guess, he will consider how effectively this selection process narrowed down the possible answer. Structural analysis follows the same sort of process. It attempts to answer ques-tions such as: can a central "most plausible" value be assumed, or estimated from the data? How much does the field vary from one location to the next? How quickly does the similarity between two readings decrease with distance? Does this depend on the orientation of the two sites? The result of the analysis is a statistical model of the random field, specifying the first two moments of the population probability function. But how can one infer the statistical moments of a whole population based on a single realization? Underlining the geostatistical analysis is an often unstated ergodic assumption (literally: work assumption), under which one presupposes that one realization is sufficient to determine properties of the ensemble. In effect, the spatial mean and covariance of the population are set to mimic that of the (partial) realization. Insofar as the population is expected to represent the natural setting, the validity of the ergodic assumption depends greatly on the quality of the sampling. In linear geostatistics, the estimate at an unsampled site is described as a linear combination of the neighboring data points. Once a statistical model is defined, the optimal weights can be calculated by minimizing the estimation error. The interpolation weights are dependent only on the estimation location and its distance to the data points; they do not depend on the actual value of the data points, only on their position. The kriging weights represent the expected similarity between the estimation location and the neighboring sites. To summarize, the kriging procedure estimates the regionalized variable by narrowing the en-semble of plausible solutions in two steps: first the ensemble is limited to functions sharing a certain structure; second, the ensemble is limited by requiring that all solutions obey the data. The final es-timate is obtained by averaging this reduced population, and its variance is the estimation variance. In point-estimation kriging, the data are fitted perfectly, even if a very complex field is required to do so. Unlike methods commonly used in geophysical inverse problems, the objective function does not include a penalty for complex solutions; the minimization of the misfit is not limited by the minimization of structure. Chapter 5. Kriging 38 5.2 Ordinary kriging Ordinary kriging is by far the most common form of optimal linear estimation for spatial data. Its mathematical development is straightforward, and the computation of the kriging weights is inexpensive. The kriging estimator is a linear combination of the field values at the sampling locations. The coefficients are chosen so that the estimator is unbiased, and so that the variance of the estimation error is minimal. The minimization problem is well-posed, but requires strong assumptions about the first two moments of the random function. 5.2.1 Structural assumption Ordinary kriging assumes the random function to be second order stationary, i.e. to have a spatially constant mean and a covariance function that is invariant under translation. E[Z{x)] = u. (5.1) Cov[Z(x), Z(x + h)] = C(h), (5.2) where h = x — x' is the distance vector between two locations. The model can be constrained further by assuming an isotropic field, in which case the two-point covariance depends only on the scalar distance h = \\x — x'\\. The requirement for second-order stationarity can be somewhat relaxed without affecting the kriging system by assuming the random function to honor what is known as the intrinsic hypothesis: E[Z{x)] = u. (5.3) Var[Z(x) -Z(x + h)} = 27(h), (5.4) that is, the mean of the random field is constant and the variance of the difference between two loca-tions depends only on the distance separating them. The function 7 ( - ) is known as the variogram 2, and is defined as 1(x,x') = ^E[(Z(x)-Z(x'))2]. (5.5) The variogram is essentially a measure of how dissimilarity between values of the attribute increases with the sampling distance. The intrinsic hypothesis is less restrictive than the assumption of second order stationarity; a stationary model honors the intrinsic hypothesis, but the converse is not necessarily true. Under second-order stationarity, the variogram and the covariance function are linearly related: 7(h) = C(0) - C(h). 2Originally, 7 as defined in (5.5) was referred to as the semivariograrn while the term variogram was reserved for 27. This convention has fallen into disuse and 7 is now commonly referred to as the variogram, despite the factor 1/2. Chapter 5. Kriging 39 5.2.2 Unbiasedness Consider the estimator Z(x$) to be a linear combination of the value of the field at the measurement location 3 N (5.6) i=l The interpolation problem reduces to the selection of weights that will warrant the desired properties of unbiasedness and minimal estimation variance. The unbiasedness of the estimator requires E Z(x0) - Z(x0) 0 N ^\iZ(xi) .t=i - E[Z(x 0)] JV which will only hold for any value of the mean if Y17=i A* = 1. 5.2.3 Estimation variance For the estimation location XQ, the variance of the estimation error is <T2(X0) = Var \z(x0) - Z(x0) . By (5.6) a2(xo) = Var N - Z(x0) (5.7) (5.8) From (5.8) and the assumption of stationarity (5.1) it can be shown that the estimate variance depends on the kriging weights and on the covariance function such that N N N a\x0) = C(0) - 2 Xi C ^ o - Xi) + XiXi °(Xi ~ x^ • i=l j = l (5.9) 2=1 The intermediate algebra is straightforward, and follows closely the derivation of the Bayesian estimate variance presented in Section C.4. It can also be found in any standard geostatistics textbook (e.g. Kitanidis, 1997; Olea, 1999). ''Note that the estimator Z(xo) = Z(XQ,UI) is a random function. As will be shown, the weights A* depend only on the location of the measurements, {.TI, . . . , xjv}, not on their values. Therefore, Z(xo) is a general expression, independent of the particular realization from which the data are sampled. The desired estimate z(x, UJO) is obtained when the data values are substituted in (5.6) i.e. z(xo) = X ^ i l i >^ z(xt). z(xo) is a deterministic function. Chapter 5. Kriging 40 5.2.4 Normal equations The minimization of the estimation error under the unbiasedness constraint is solved by Lagrange's method of multipliers. The objective function is L(\1,X2,...,Xi,P) = a2(x0) + 2/3^J2\i - 1 j . (5.10) By (5.9), this function is quadratic in term of the kriging weights Ai's, with a positive-definite coefficient (the covariance function). Consequently, the solution to the minimization problem is unique, and the optimum is a global minimum. Differentiation of this lagrangian with respect to the weights A '^s and the Lagrange multiplier (3 yields the following system of equations JV Y^KC(xi-xj)+/3 = C(x0-Xj), f o r j = l . . . J V (5.11) N 5~jAi=l. (5-12) j=i 5.2.5 Optimal estimation error variance Multiplying each equation of the system (5.12) by Xj and summing the N equations gives N N N' 1=1 j=l j=l Inserting the left-hand side of (5.13) in (5.9), and substituting the optimal weights obtained from (5.12) gives the optimal estimate variance N a*2{x0) = Cov[0] - ^ Xi C(x0 - xi) - (3. (5.14) i=i C(0) is the variance of the random function Z(x); it is the variance of the estimation error if no data are available to reduce the spread of the population. Each data point within correlation range of the estimation location will affect the estimate value, and lower the estimate variance by AiC(xrj — xi)-5.3 Model selection: the experimental variogram In the preceding derivation of the kriging estimate the covariance (or variogram) of the random field is supposed to be known exactly. In practice, it is estimated from the data. It is important to keep in mind that the random function Z(x, w) does not represent a natural "random process"; it is a mental construction that attempts to qualify our knowledge (and ignorance) of a deterministic field. A set of fictional realizations is constructed in the hope that their average behavior will resemble the natural field, and that their dispersion will cover its plausible range of variation. In this sense, the covariance function is a prescription more than an estimation. It should not, however, be arbitrary; Chapter 5. Kriging 41 the more closely the covariance function resembles the "natural variability" of the field, the more effectively the random population will qualify our knowledge of the natural setting. If one assumes that the field is second-order stationary, its covariance can be estimated from the data. A l l possible pairs of measurement are compared, and grouped into bins according to the distance separating them. The variance within each bin gives an estimate of the covariance for that distance value: C ^ = AT E *te)*(y.) - E E (5-15) fc i = l i=l i = l where hk is a distance value representative of the fc-th bin, and Nk the number of pairs in the bin. If the field is non-stationary but obeys the intrinsic hypothesis, the experimental variogram can be calculated instead: ^ k i = i The experimental variogram of a second-order stationary field will level to a maximal value, the sill, representing the field variance. The distance at which this plateau is reached is called the range, and is a measure of the autocorrelation length of the field. The j/-axis intersect of the variogram is a measure of the microscale variability of the field. It can represent measurement errors, or a small-scale signal unresolved by the current sampling intervals. In geostatistical literature, this microscale variability is known as a nugget effect, an artifact of the early use of geostatistics as a tool for the estimation of ore deposits. The function used to fit the experimental variogram is selected from a family of functions that have the property to yield positive-definite covariance matrices. Two such functions are the spherical model (a2(¥-&) iorO<h<r , s •y(h) = { ~ ~ , (5.17) (cr for h > r (where a2 is the variance and r the range), and the gaussian model 7 W = ^ ( l - e x p ( ^ ) ) . (5.18) The power model is an example of a non-stationary model: -y(h) = 0hs, (5.19) where the two parameters to be fitted are the coefficient 0 > 0 and the exponent 0 < s < 2. More such functions can be found in any geostatistics textbook (e.g. Kitanidis, 1997; Olea, 1999). 5.4 Model validation: analysis of orthonormal residuals Section 5.3 describes how an analytic expression for the covariance can be estimated from the data. In the derivation of the estimation error variance, the covariance function is treated as a deterministic function and is supposed to be known exactly. Consequently, the variance of the estimate does not Chapter 5. Kriging 42 represent the uncertainty related to the selection of an appropriate covariance model. To be precise, two distinct sources of uncertainty are overlooked: the choice of the analytical function used to describe the covariance, and 9, the set of parameters used to fit this model to the experimental variogram. The probability assigned to a realization depends on the fitting function, the parameters, and the data. The effect of the parameter selection on the prediction and its variance can be accounted for by considering the compound (or Bayesian) distribution (Kitanidis, 2005) p(Z(x,u)\F,y) = Jp(Z(x,aj)\6,F,y)p(9\F,y)dO, (5.20) where y denotes the data and p{Z{x,oj)\9,Jr,y) is the probability of a realization given a fitting function !F and a set of parameters 9. The compound distribution p(Z(x,u>)\F,y), on the other hand, includes the uncertainty attributed to the model parameters 9. For this reason, it is the compound distribution that should be used for prediction Z(x) = E[Z(x, u)\7, y] = J E[Z(x, u)\9, T, y] p(9\f, y)d9. (5.21) Assuming the parameters to be known exactly is equivalent to declaring their posterior pdf p(9\J-,y) to be a Dirac delta function, and assuming E[Z(x,cj)\F,y] ~ E[Z(x,u)\9,F,y}. (5.22) Kitanidis (2005) argues that this is often an adequate approximation, as E[Z(X,UJ)\9:T,y] is gen-erally not very sensitive to 6 over the range of variability of 9. Without considering a fully Bayesian approach to kriging based on the compound distribu-tion, the quality of the statistical model can be assessed by the analysis of orthonormal residuals. Residuals are differences between observations and model predictions. They can be calculated by predicting the value of the field at a sampled location while ignoring the corresponding data point. Orthonormal residuals are calculated in a similar fashion, except that the cross-validation follows a strict procedure (Kitanidis, 1991): 1. Organize the data points in an ordered sequence: { z(xi), 2(2:2)1 • • • , z(xn) }. 2. Using the model to be tested, obtain the kriging estimate at the location of the second data point X2 using the first measurement z(xi) as data. 3. Calculate the associated estimation error, the residual 82 = z(x2) — Z(x2), and normalize it by the estimation standard deviation to obtain the first orthonormal residual, 62 = 82/0-2-4. Estimate the field at the location of the third data point, using the first and second measure-ments as data. 5. Calculate the residual £ 3 = z(xz) — Z{x%) and the orthonormal residual 63 = 83/0-3. 6. Repeat steps 4-5 for every other data point z(xk), each time using the k — 1 preceding measurements to calculate the orthonormal residual e^ . Chapter 5. Kriging 43 Statistic Decision rule: Reject the model if ... Mean Qi = 1 " k=2 | Q i | > 2/y/n-l Variance Q2 = 1 n \ E 4 \Q2 - 1| > 2.8/>/n - 1 fc=2 T A B L E 5.1: Orthonormal residuals statistics and decision rules, for a 5% significance level, (n — 1)Q2 follows a x2 distribution with parameter (n— 1). The decision rule presented here is an approximation valid for large populations (n > 40): for smaller samples, the upper and lower percentiles must; be obtained from X 2 tables. Just as the estimator Z(XQ) was a random function until the actual values of the data points were substituted in the linear combination 2~2u=i ^» Z(xi), the residuals can be considered as random variables, and their ensemble statistics can be computed. If the statistical model is correct, the innovations e/t should have a mean of 0, a variance of 1, and be uncorrelated; that is to say, they should be orthonormal. Kitanidis (1991) proposes a series of statistics to test that distribution of the experimental orthonormal residuals follow these characteristics (Table 5.1). The orthonormal residuals should have a gaussian distribution, since minimum-variance linear unbiased estimation implicitly assumes that the estimation errors are approximately normal. Tests for the goodness of fit such as the Lilliefors test can be used to refute this hypothesis. By these tests, the computation of the orthonormal residuals allows the user to verify how well the statistical model represents the data. The rejection of the model by any of these tests does not imply that the estimate is inexact or invalid, but that the statistical model does not represent faithfully the complexity of the data. Given the specificity of the assumptions behind the model, this is a likely occurrence. It can also warn the practitioner against a poor parametrization of the covariance function, which is easier to remedy. Using these statistics, the performance of competing models can be objectively compared, and the best model can be chosen. 5.5 Application: interpolation of bed topography data As mentioned in Section 3.2, two data sets are joined to obtain a complete bed map: the 1981 20-m resolution D E M , covering the area situated outside the 1981 glacier outline, and the 1994-1997 radar data, covering the area situated inside the 1997 glacier limits (Figure 5.1). The procedure for collecting and processing the radar data is described in Flowers and Clarke (1999). I used two independent methods to evaluate the covariance parameters: least-squares fitting of an omnidirectional variogram, and estimation by approximate maximal likelihood ( A M L ) (Pardo-Iguzquiza and Dowd, 1997). The experimental variogram was calculated using my own routine, Chapter 5. Kriging 44 F I G U R E 5 . 1 : Data points for the kriging of the bed topography map. The x's mark radar sounding locations, and the dots mark elevation data collected from the 1981 air photos (1981 outline in red). For computational efficiency, not all the D E M points were used; the 1951 glacier area, being the zone of interest, was sampled more densely. inspired by the function VARIOGRAM of Moacir Cornetti & Armando Remacre posted on the Math-works file exchange. The approximate maximal likelihood estimates were obtained using a (slightly modified) F O R T R A N routine published by Pardo-Igiizquiza (1998). I use the COKRI function of Marcotte (1991) in my kriging algorithm and for the calculation of orthonormal residuals. Figure 5.2a shows the omnidirectional experimental variogram and the fitting function that provided the first set of parameters. The bed elevation values were detrended before the analysis. Using a spherical model (5.17), the least-squares fit give a correlation range of 890m and a variance of 957 m 2 . F i gure 5.2b and Figure 5.2c present two directional variograms in the two principal directions detected by the A M L algorithm. The data have been detrended based on the linear drift coefficients calculated by A M L . The solid lines give the estimated variogram functions. The fact that the A M L estimates do not fit the experimental variograms is interesting. Fitting experimental variograms might be a more intuitive approach to the determination of covariance parameters, but it is not an optimal one. The fact that estimates maximizing the likelihood of the data given the proposed model do not fit the experimental variogram reveals the limitations of this tool. The superiority of the A M L estimates is not only theoretical, but can be verified by orthonormal residuals analysis. Table 5.2 shows the mean, variance and the Lilliefors test for normality of the orthonormal residuals for the two models. The A M L parameters fare better on all fronts. The resulting bed map is presented in Figure 5.3. Chapter 5. Kriging 45 1200 1000 800 €. 600 400 800 1 7 5 . 7 200 AMLE range x AMLE ' range y 500 1000 distance (m) 500 1000 distance (m) 1500 F I G U R E 5 . 2 : Statistical analysis of bed topography data, (a) Omnidirectional variogram of detrended bed topography. The solid line gives the covariance function adopted in the "best fit" scenario, (b) Directional variogram, 175.7°. (c) Directional variogram, 85.7°. The solid line give the anisotropic covariance function obtained by approximate maximal likelihood. The A M L estimates do not fit the experimental variogram. Nevertheless, this statistical model gives better results when tested by orthonormal residuals analysis. Mean Variance Normality Model N Qy Qlmax Test Q2 Q2min Qimax Test L Lmax Test selection least- 389 0.08 0.10 0 0.59 0.86 1.14 1 0.14 0.05 1 squares A M L 388 -0.03 0.10 0 0.98 0.86 1.14 0 0.11 0.05 1 T A B L E 5 . 2 : Results of orthonormal residuals analysis for two statistical models of the bed topography. A positive test value implies the rejection of the statistical model. The A M L residuals have a mean closer to 0 and a variance closer to 1 than the least-squares fitted omnidirectional model. The residuals are also closer to a normal distribution, even if they fail the Lilliefors test for normality. F I G U R E 5 .3 : Kr iged bed topography map. In perspective, looking North-West, w i th vertical axis exagger-ated four times and elevation contour 20m apart (left). Top view of detrended bed elevation, wi th elevation contours 50m apart (middle). Top view of kriging standard deviation (right). Chapter 5. Kriging 47 5.6 Including a priori information: Bayesian kriging This section is based on the Bayesian kriging algorithm developed by Omre (1987). This extension of ordinary kriging incorporates qualified guesses in the inference by assuming the mean of the random field to be known within a specified range of uncertainty. The estimate obeys the data wherever they are available, and resembles the prior guess when they are not. The degree of certainty attributed to the guess is reflected in the final estimate variance. This algorithm is described as Bayesian because it includes a prior model of the field, and accounts for its uncertainty. However, it does not account for the uncertainty related to the selection of covariance parameters. Just as does ordinary kriging, it assumes the parameters 0 to be known exactly, and approximates the compound distribution p(Z(x,oj)\J-,y) ~ p(Z(x,oj)\9,Jr,y). Other Bayesian algorithms have been designed to quantify the effect of model selection on the prediction variance (Kitanidis, 2005). The algorithm of Omre's (1987) has the advantage of being similar to ordinary kriging, which allowed me to implement it by building on published ordinary kriging routines. I did not see the justification for a fully Bayesian approach, and was not inspired to attempt to program one. 5.6.1 S t r u c t u r a l assumptions Consider a random function {M(x, 4>); x G V, (f> G $} (thereafter referred to simply as M(x)), where T> is the spatial domain of variation and <E> is the sample space of the random variable <j>. Suppose its first two moments known a priori: E[M(x)]=m{x) Cov[M{x),M{x')} = CM{x,x') x, x' G V. xev (5.23) (5.24) Assuming that the random field Z(x) obeys x eT> (5.25) (5.26) Cov [Z{x), Z{x')\M{x"); x" G P] = Cz]M{x - x') x, x' € P, its first two moments are E[Z(x)] = E[E[Z{x)\M(x"); x" G = CLQ + p,M{x) V]] (5.27) Cav[Z(x), Z{x')] = E[Cov[Z(x), Z(x')\M(x"); x" G V]] + Cov[E[Z(a;)|M(a:"); x" £ V] ,E[Z(x)\M{x"); x" G V]] = CZ\M(X - x') + CM{x, X'). (5.28) Equation (5.27) shows that the field U-M{X), the mean of the random function M(x), can be in-terpreted as a qualified guess on the expected value of Z(x). The constant ao is unspecified, and Chapter 5. Kriging 48 the prior is only meant to represent the shape of the field, with no control on its level. Equa-tions (5.25) and (5.26) mean that Z(x) is second-order stationary around its expected function {an + U-M{X)\ % € T>}. The fluctuation is the variation of the field Z(x) around the prior estimate U-M{X)> and is noted {Y(x) = Z(x) — ^ M ( ^ ) ; X £ V}. 5.6.2 Bayesian block kriging algorithm The introductory section on ordinary kriging was limited to point kriging, that is the estimation of the field at a point location. Block kriging on the other hand, estimates the average value of the field over a given volume support. It was historically the first form of kriging to be developed, designed to estimate the volume of ore deposits. To construct a glacier D E M that represents faithfully the ice volume, it is preferable to estimate the average value of ice thickness over a grid cell rather than at a unique point. Block-estimation is a straightforward generalization of point-estimation, and the derivation of the Bayesian block kriging system follows the step outlined in Section 5.2: definition of a linear estimator, establishment of conditions for unbiasedness, minimization of estimation error with unbiasedness constraint, derivation of corresponding optimal variance, solution of kriging system. The details of the derivation are presented in Appendix C and the results are summarized in Table 5.3. Chapter 5. Kriging 49 Estimator JV Z(S0) = J2Xi (Z(xi)-fiM(xi))+nM(So). (5.29) 2 = 1 Normal equations A' [Cz\M(xi -Xj) + CM(xi,Xj)] +0 = [CZ\M(So,Xj) + CM{So,Xj)] , for j = 1...N i=l (5.30) N 5 > = 1 (5.31) i=i Optimal estimation error variance TV a*2(S0) = [Cz\M(P) + CM(S0,So)] [Cz{M(xi,So) + CM(xi:So)]-0. (5.32) i=l where VM{SO) = 7 T / fiM(x)dx (5.33) CZ\M(xI,SJ) = j - f Cav[Z{xi),Z(x)\M\dx (5.34) 3 J Sj CZLM(SI,SJ) = ^ -- f I Cov[Z(x),Z(x')\M]dxdx' (5.35) CM(xi,Sj) = ±- [ Cov[M{Xi),M(x)}dx (5.36) b3 JSj CM(Si,Si)=eV / / C o v ^ r r / ) ^ ' ) ] ^ ^ ' (5.37) T A B L E 5.3: The Bayesian block kriging system. See Appendix C for details. As one can see, the system of equations for the Bayesian kriging weights has the same general form as that for ordinary kriging. Furthermore, rewriting the Bayesian block kriging estimator (5.29) in term of the fluctuations, Y(x) = Z(x) - HM(X), shows that it is identical to the ordinary Chapter 5. Kriging 50 block kriging estimator for Y(So) N Z(So) - VM(SO) = YXi (z(xi) -i=l N Y(s0) = j2^y(xi)-i=i What is, then, the difference between including the a priori guess by Bayesian kriging, and simply estimating the difference between the data and the guess by ordinary kriging? Bayesian kriging can account for guesses with spatially variable uncertainty. When constructing a prior model, the analyst might know very well some areas of the domain, while others are more uncertain. This varying degree of trust is accounted for in the inference by the location-specific prior covariance (CM(XI, xj)). If the precision of the a priori guess is constant throughout the domain, Bayesian kriging reduces to the ordinary kriging of the random fluctuation function {Y(x);x 6 D}. In this case, removing the guess from the data can be seen as a complex form of detrending. The resulting field is second-order stationary and can be estimated by ordinary kriging. If the precision of the a priori guess varies in space, however, the fluctuation field is not second-order stationary, and Bayesian kriging must be used. The covariance function is then divided into a location-specific term and a translation-symmetric term Cz(xi,Xj) = CZ\M{XI ~ xj) + CM(xi,Xj). The final estimation variance (5.32) represents the uncertainty of the guess and that due to the distance of the neighboring data points. 5.6.3 Model selection To guarantee a positive estimation variance, the covariance function Cz{x,x') = CM{X,X') + CZ\M(X~~X') must be conditionally positive-definite. This will trivially be achieved if both CM{X, X') and CZ\M{X ~ X') are conditionally positive-definite. In the case of CZ\M{x — x'): this is insured by choosing the function used to fit the experimental variogram 7z|M(h) (defined below) from the standard set of elementary models used in geostatistics. The spherical and gaussian models pre-sented in Section 5.3 are two examples. The set of acceptable functions for the prior covariance function CM(X, X') can be written as CM(x, x') = O-M{X)O-M{X')PM{X - x') (5.38) where o-\I(x) = Var[M(x)] is a non-zero variance field defined over the domain T>, and PM(X — x') is a location-invariant positive-definite correlation function. Omre (1987) proposes an unbiased estimator for the conditional variogram Jz\M'-7z|Jtf(h) = 2^f Yl {[Z{xi)-Z(xj)}2-[fiM(xi)-^M{xj)}2-2'yM{xi,xj)}, (5.39) were is an ensemble of Nh pairs such ~ h. Under the assumptions (5.27) and (5.28), the conditional variogram is equivalent to the conditional covariance CZ\M(X,X') = CZ\M(®) ~ Chapter 5. Kriging 51 1Z\M{X,X>) (Omre and Halvorsen, 1989). As discussed above, the function used to fit the ex-perimental variogram must be positive-definite. However, nothing insures the positivity of the estimator (5.39). If the prior model has more spatial variance than the data, it is likely that Ah^M^i) - M M ( Z J ) ] 2 will be greater than Y^(i,j)<=AhlZ(xi) ~ z(xj)}2 for certain values of h, yielding a negative JZIMO3)- This possible failure of the variogram estimator is a consequence of the weakness of the ergodic assumption. When using data to calculate an experimental variogram, the spatial characteristics of the field are used to estimate the variance of its underlying random variable. While it is clear that the variance of the target field, being the sum of the prior and con-ditional variances, should be greater than the prior variance, nothing guarantees that the a priori guess will have less spatial variance than the data. Despite this possible source of inconvenience, the experimental variogram (5.39) remains a useful tool for the structural analysis of the target field. This concludes the chapter on kriging theory. In the next chapter, Bayesian kriging is used in an innovative approach to time—space interpolation. A time-evolving D E M of the glacier is generated, starting from the photogrammetric D E M s , and updated year after year by the ground survey data. 52 Chapter 6 Time-evolving D E M 6.1 Description of the algorithm I use the Bayesian kriging technique presented in the previous chapter to iteratively generate a time series of D E M s for Trapridge Glacier. To my knowledge, this is the first time this technique has been adapted for time-space analysis. The basic idea is simple: at each time step, the target field from the previous time step is used as prior model, and the available data are used to update the field. The resulting estimate is then passed on as prior model for the next time step, with its estimation variance used as model uncertainty. Figure 6.1 illustrates this procedure. When used in this time series context, the first statistical assumption of Bayesian kriging (5.25) states that the expected field at a given time step t is the value of the field at the previous time step t — l , 1 plus a constant shift: E[Zt\Zt-i(x"); x" £ P ] = a 0 + Zt^(x) x € V. (6.1) Given that the glacier thickness is not expected to change drastically from year to year, this is a reasonable assumption. The "DC component" an accounts for a mean change in ice thickness. A more specific model could have added a linear term to this drift, to represent the different rates of thinning between the lower and higher regions of the glacier. What degree of structure should be accounted for by the drift, and what part should be left to the random component of the field is always an issue for geostatisticians. The principle of parsimony suggests that the drift should be kept as simple as possible, and that extra parameters should be added only when necessary. In this case, using a linear drift would not only be extraneous but also unwise; because the data are unevenly sampled between the lower and upper basin, assuming a linear trend based mainly on lower basin data could lead to unrealistic behavior in the undersampled upper basin. I therefore limited the drift to a constant value. It does not account for the variation in the annual change along the length of the glacier, but it does insure that the distribution of the fluctuations will be centered. Note that the parameter an, unlike the mean drift in ordinary kriging, does not appear in the expression for the estimator (equation 5.29). It affects the kriging system through the unbiasedness constraint Y^i=\ x i — L which is only necessary because an ^ 0. Assuming a drift value of 0 makes 1 Strictly speaking, Bayesian kriging does not require a causal relationship between the prior and the field. The sequence should rather follow information gradients; going from a state where the field is well known, to a state where it is undetermined. As such, this algorithm can be used to iterate in the two time directions. On the few occasions when the elevation field for the next year is better constrained than for the previous one, I reversed the chronological sequence. For example, the 1970 photogrammetric model is used as prior for the inference of the 1969 elevation field (see Table 6.1). Chapter 6. Time-evolving DEM 53 PRIOR < — RESULT KRIGED FIELD Z(x) F I G U R E 6 . 1 : One time step in the generation of the D E M time series. The 1995 D E M (a) is updated by survey data from 1996 ( b ) . The variance of the resulting D E M (c) depends on the confidence attributed to the prior estimate and on the data coverage. In the next iteration, the 1996 D E M wi l l be used as prior for the kriging of the 1997 data, its standard deviation giving the prior uncertainty. Note that the standard deviation in the upper basin is low, despite the fact that very few measurements were taken there in 1996. This is because the prior was very well defined in that area. The 1996 model "remembers" the 1995 surveying of the upper basin, ( d ) shows the change in elevation between 1995 and 1996. Contours give the ice thickness field, in 20 m increments. Chapter 6. Time-evolving DEM 54 Year Prior model used 1969 1970 1971 1972 1974 1980 1981 1982 1983 1984 photogrammetric photogrammetric kriged photogrammetric kriged photogrammetric photogrammetric kriged kriged kriged 1970 1970 1970 1972 1972 1981 1981 1981 1982 1983 y kriged [y-i) 2005 kriged 2004 T A B L E 6 . 1 : Prior model used for each year of the time series. The 1970 photogrammetric D E M , being the closest available, was used as prior for the kriging of the 1969 and 1970 survey data. For 1971, the kriged 1970 D E M was used rather than the photogrammetric model because the kriged field includes the information from the survey data, and is therefore more precise. For the 1980 prior model, 1 hesitated between the kriged 1974 D E M , or the photogrammetric 1981 D E M . I finally chose the photogrammetric 1981 D E M , preferring proximity in time over precision. the estimator unbiased irrespective of the kriging weight values, and no extra constraint has to be imposed during their evaluation. In addition to propagating information in time, this Bayesian algorithm allows the inclusion of secondary sources of information. Here, the primary data come from ground surveys. They are the most precise data available, and yield the most consistent time series. O n the other hand, the photogrammetric data, even if only available for a few years and less precise, have the advantage of covering the whole domain. They are included in the inference as prior estimates of the glacier surface. Survey data are used to correct the elevation field, and to update this initial model in subsequent years. The model used as prior for each time step is specified in Table 6.1. I had several options as to what the kriged field Z could be. The desired final result for a D E M is an elevation field, but the primary quantity of interest for glaciology is the ice thickness, kriging the elevation field uses the survey measurements directly, without intermediate processing. The ice thickness field requires subtraction (and interpolation) of the bed topography map. Omre (1987) mentions that the Bayesian kriging development is exact if the prior and kriged random fields are jointly gaussian. Since only one realization is available, the joint normality of the random fields is assumed, rather than verified. However, the estimation of the covariance parameters is based on the assumption that the spatial variability of this unique realization represents the variability of the random function population. To strengthen this ergodic assumption, it is preferable to use a data Chapter 6. Time-evolving DEM 55 set that is spatially gaussian. By subtraction of the bed topography, the downward trend of the elevation data is naturally filtered, and the ice thickness data have a more symmetrical distribution than the elevation data. A trend persists, as the glacier tapers downstream, but is much weaker. I therefore chose to krige the ice thickness field. Flowers and Clarke (1999) kriged the natural logarithm of the ice thickness field. While this has the big advantage of insuring positivity, the non-linearity of the transformation can lead to uneven weighting when ice thickness values are small. This was not an issue for Flowers and Clarke (1999) because no areas of thin ice remained in the years considered (1994-1997). In earlier years, however, surveying of the proglacial apron yielded small ice thickness values. When using log-transformed data, these locations have a large numerical impact on the kriging due to the non-linearity of the transformation. Using the ice thickness field directly yields better results. Negative estimation values are rare, and only occur outside the glacier outline. I did not detrend the ice thickness data before kriging, because it can generate biases in the estimation (Olea, 1999). It is also not as important for Bayesian kriging as it is for ordinary kriging. The prior and the target field have to be jointly gaussian; their individual normality is a sufficient, condition, but not a necessary one. It suffices that the fluctuation of the target field around the prior be normally distributed. The ice thickness trend is filtered out by subtraction when calculating the fluctuations, and therefore its presence is not an inconvenience. 6.2 Incorporating information about the glacier limits Figure 6.2 presents a few D E M s obtained by Bayesian kriging (the complete set can be found in Appendix E) . The estimated glacier outlines (Section 4.3) are superposed for reference. As one can see, the kriged ice thickness field does not reproduce the glacier boundaries; ice thickness beyond the outline can be considerably larger than zero. Because measurements are only taken on the ice surface, nothing controls the ice thickness field beyond the limits of the glacier. Data taken near the glacier boundary influences the thickness field over a distance equal to the correlation range (500m; see next section). This can cause the ice thickness field to advance beyond its supposed limits. A remedy commonly used in geostatistics is to incorporate fake data to represent "soft" information not contained in the original data set. In this case, the ice thickness field could be forced to zero at the outline. However the presence of these zero-thickness data points would make the ice thickness decrease gradually from the nearest data point to the glacier boundary, resulting in a tapered terminus wherever data points are not close to the boundary (i.e. along most of the outline). The glacier terminus over the last twenty years has been cliff-like rather than tapered, and this approach would fail to reproduce this important feature of the glacier geometry. Bayesian kriging was developed precisely to incorporate this sort of soft information in the inference process. Information concerning the glacier limits could be included in the inference by setting the prior field to zero outside the outline. Unfortunately, this promising procedure induced erratic behavior at the glacier boundary. Data points situated between two subsequent outlines have a prior estimated value of zero, as they are situated beyond the previous outline, and a current ice thickness of ~50m. The fluctuation values at these location are therefore abnormally large Chapter 6. Time-evolving DEM 56 0 20 40 60 80 100 120 0 6 12 18 24 Ice thickness (m) Standard deviation (m) F I G U R E 6 . 2 : Ice thickness field and kriging standard deviation. No information on the glacier boundary is included in the interpolation, and therefore nothing stops the ice thickness field from growing beyond the supposed limits of the glacier. The cliff that formed at the glacier terminus during the 1990s is a statistical singularity that wi l l not be reproduced unless it is captured by the data. The scarcity of data, at the glacier boundary makes its resolution problematic. The issue is solved by setting a posteriori the ice thickness field beyond the outline to zero (see Figure 6.3). Note how the survey of the upper basin in 1991 and 1995 decreased the standard deviation in that area in 1992 and 1996. This is a good il lustration of the difference between Bayesian and ordinary kriging. In the iterative Bayesian algorithm, information is passed from one year to the next, and the estimate variance accounts for this prior knowledge. Chapter 6. Time-evolving DEM 57 compared to the rest of the field, and create peaks in the kriged ice surface. The transition from ice to no-ice at the glacier boundary is a statistical anomaly, a singularity. The softest way to incorporate this information is simply to ignore it during the estimation, letting the ice thickness field evolve smoothly according to the data as if no singularity existed, and to impose the boundary after the fact by forcing all ice thickness values outside the outline to zero. This avoids the numerical irregularities associated with the boundary, and replicates the sharp cliffs that are observed at the glacier terminus (Figure 6.3). This approach was also adopted by Flowers (2000). Chapter 6. Time-evolving DEM 58 F I G U R E 6 . 3 : In the absence of data beyond the glacier outline, the kriged ice thickness decreases gradually, creating a tapered terminus. To render the cliffs that formed at the terminus in the 1990s, the field is forced to zero outside the glacier outline. The results are illustrated here for the 2004 model. Chapter 6. Time-evolving DEM 59 6.3 Selection of covariance parameters and orthonormal residuals analysis The assumption of second order stationarity around the prior model (equations 5.25 and 5.26) defines the general structure of the statistical model. The specifics of the model, the covariance parameters, remain to be determined. Bayesian kriging has the disadvantage of requiring analytical expressions for two covariance functions: the covariance of the field conditional on the prior CZ\M(X,X'), a n d the covariance of the prior CM{X,X')- When generating a time series, the prior and the target field are the same physical field, and their spatial characteristics are expected to be similar. The variance field, describing how much is known about the target field, might vary from one year to the next, but the autocorrelation lengths, depending primarily on the physical properties of the field, should not change. If the prior and target fields share the same correlation function, the conditional covariance must also have the same spatial behavior: Cz]M(x - x') = Cz(x,x') - CM(x,x') by (5.28) Ct^1(x-x') = Ct(x,x')-Ct-i(x,x') (M = Zt-i) = [at{x)at(x')pt{x-x')] - [at-i{x)at-i{x')pt-i{x - x')} = o-t(x) at(x') - CTt-i(x) o-t-i{x') p(x - x'). (pt = pt-i) (6.2) Assuming the covariance spatial structure to be constant in time greatly simplifies the analysis, re-ducing the number of parameters to estimate from 93 (3 per year) to 33: the 31 conditional variances C t | t _i (0) , the correlation range, and the microscale variability 2 . A n analytical expression for the correlation function must also be chosen (spherical, gaussian, exponential, etc.). These parameters can be estimated by fitting a function to the experimental conditional variograms (equation 5.39) calculated from the data 3 . Before doing so, it is useful to consider the significance of the conditional variance cr 2^_ 1 = C t | t _ 1 (0 ) . Separating the conditional variance and the correlation function in (6.2) and dividing both side by the correlation functions yields „2 at\t-i at(x)at(x') - at-i{x)at-i(x') . (6.3) The variance field a2(x) gives the spread of the random function population {Z(x, cu); x € D, UJ € Cl} before it is restricted to obey the data; it is the maximal value of the estimation variance a^2(xQ) — o~2t{xo) - ^ZiLi XiCt{xi,xo) — Q (5.32). Equation (6.3) shows that the conditional variance o-2^t_1 represents the loss of information (or gain in variance) for every iteration where no data are available. It is the increase in variance that will occur at all locations that are more than a correlation range away from the nearest data point. In other words, it represents the rate at which information on the glacier geometry is dissipated. As the correlation range, it is an intrinsic characteristic of the 2 O r , in traditional geostatistical terms, the range, and the nugget. The variance is also known as the sill 3Note that the covariance of the target field Ct (or of the prior Ct-i) cannot be estimated using an experimental variogram because the field is not second-order stationary; its variogram depends on the location of the points, not only on their distance -ft(x.x') ^ j(h). Chapter 6. Time-evolving DEM 60 field. Instead of adjusting this value every year based on the experimental conditional variogram, I selected a central value that made physical sense and fitted most variograms. Figure 6.4 shows a selection of experimental conditional variograms with the selected spherical model. A gaussian model would have fitted the experimental variogram more closely, but repeatedly led to poorly conditioned matrices and negative variance values. The range and sill parameters were chosen by trial and error. The quality of the fit was used as a first estimate, but the final choice was based on the analysis of the orthonormal residuals. To compare model performance over the 31 years of the time series, the statistics Q\ and \Q2 — 1| for all time steps are averaged, and the test results are summed. More than 30 different combinations of parameters were compared (see Table 6.2 for a subset). The best results were obtained using a range of 500m and a sill of 100m 2 . The complete set of orthonormal residuals results for these parameters are presented in Table D . l . As seen in Figure 6.4, this model captures well the central tendency of the time series variograms. However, a conditional variance of 100 m 2 implies that there will be a standard deviation increase of 10 m per year for each location where no data are available within a correlation range. Given that the standard deviation is used as uncertainty estimate, this is a rather high figure. As an alternative to the fitting of the experimental variograms by a time-independent model, I tried using the program AMLE2D (Pardo-Iguzquiza, 1998) to infer annual covariance parameters from the fluctuation data. Because covariance is invariant in translation, the data covariance Cov[Z(x), Z(x')] is identical to the fluctuation data covariance Cov[Z(x) — U-M{X), Z(X') — HM{^')\-AMLE2D estimates covariance parameters that maximize the likelihood of the data, if they are as-sumed to be second-order stationary. As discussed in Section 5.6, if the prior variance is assumed to be constant in space, the data are second-order stationary and Bayesian kriging reduces to ordinary kriging. The parameters estimated by AMLE2D are therefore those that would be the most probable, if the variance of the prior was constant. Despite the fact that this is not the case for this Bayesian estimate, it might nevertheless be able to capture the spatial attributes of the field. This approach also has the advantage of being automated, which is a prerequisite when 31 estimations have to be made for each run. The range and sill values estimated by AMLE2D from the fluctuation data were extremely variable from year to year. The range varied between 0.01 and 1758 m, with jumps of up to one order of magnitude from one year to the next. This is a conspicuous behavior, difficult to justify physically. To add insult to the injury, the analysis of the orthonormal residuals has showed that this time-varying model did not perform as well as the stationary one (see Table D.3 in Appendix D for the parameters and the orthonormal residuals statistics calculated for each year). After this failed attempt at adapting to Bayesian kriging a tool designed for ordinary kriging, I tried using it in its legitimate context. I estimated the fluctuation field by ordinary kriging, using AMLE2D to infer statistical parameters for each year. To compare with my Bayesian time-independent model, Ialso tried ordinary kriging with a fixed covariance. The results are presented in Table 6.2. Two conclusions can be reached from these various attempts. First, comparing orthonormal residuals calculated from a Bayesian and a ordinary model with equal variances shows that Bayesian kriging can predict the data more accurately, maintaining the estimates within the limits allowed by the model variance. Second, time varying models do not necessarily improve the estimate. While it Chapter 6. Time-evolving DEM 61 100 300 300 200 400 distance (m) F I G U R E 6.4: Experimental conditional variograms. As discussed in Section 5.6. the conditional variogram estimator is not always positive. An extreme case of this can be seen in 1995, where the whole variogram is negative. The upper basin was thoroughly surveyed in 1995 to locate radar soundings. The relatively flat topography of that region caused the data of that year to have less spatial variance than that of the preceding year, and as a consequence, the conditional variogram estimates are negative (see equation 5.39). The proposed time-independent variogram function (superposed) is a spherical model with range 500 in and sill 100 m 2 . Chapter 6. Time-evolving DEM 62 is true that I have only attempted it for ordinary kriging, or by bending the assumptions of Bayesian kriging, I doubt that the estimates would have been improved if I had fitted experimental variograms for each year. As can be seen in Figure 6.4, many time steps have experimental variograms that are partially (or entirely) negative, and therefore uninformative. Before 1987, when data are in limited numbers, the experimental variograms have only a few points, each inferred from small populations. Trying to deduce covariance parameters for each year from is therefore difficult. The time-independent model is physically justifiable, more considerate of the limitations of the data set, and give good cross-validation result for all years. Before moving on to the description of the results, I should say a few words on the estimation of uncertainty. The common usage is to consider the standard deviation of the kriging estimates to be a measure of their uncertainty; this is, however, a debatable interpretation of this statistical quantity. The main interest of geostatistics lies in its capacity to quantify "how well we know what we know." The standard deviation of the field represents the envelope around the expected model that includes most realizations of the random function. The width this envelope needs to have in order to pass the cross-validation tests depends greatly on the assumptions made on the expected value of the field. A bad prior model will require a large variance to account for the difference between the estimates and the data. When certain areas of the glacier change rapidly, the assumption of a stable shape is violated, and the variance of the field as a whole must be adjusted to account for this local failure. As a consequence, the uncertainty attributed to interpolations based on this assumption will increase. How well we know the ice level at point C between data points A and B depends on how quickly the ice thickness is expected to change, and how much. It is this "gradient of knowledge" that the variogram attempts to quantify. If the covariance function models correctly the natural variability of the field, the standard deviation is a good measure of uncertainty. However, two factors prevent this. First, the estimation of the "natural variability" of the field is greatly sensitive to the sampling. The conditional variogram for 1995 is a good example of this. The dense sampling of a region with less spatial variance than the preceding data set yielded a negative variogram, a statistical anomaly. In essence, this first limitation is due to the fallibility of the ergodic assumption. Second, to pass the cross-validation tests, the estimated variance has to be inflated to account for the failure of the assumption of stationarity. The shape of the glacier is not constant from one year to the next, and the variance has to account for this discrepancy. This will not only affect sites that are far from data, but also site well surrounded by data points, although to a lesser extent. In any case, interpreting the field variance as its "natural variability," neglects this influence. As consequence, the standard deviation does not only reflect the rate of dissipation of the information granted by data points, but also the quality of the assumption on the prior, irrespective of whether the expected model has an impact or not on the estimate. As a consequence, the standard deviation tends to underestimate our knowledge of the field. Chapter 6. Time-evolving DEM 63 Model Mean Variance Normality Range Variance (Qi) Rejections < | i - Q 2 | ) Rejections (L) Rejections variable* 0.03 1 0.41 23 0.15 26 100 50 -0.02 0 0.84 30 0.15 25 200 50 0.00 2 0.70 26 0.15 28 250 50 -0.03 2 0.64 23 0.15 27 300 25 -0.01 2 0.70 26 0.16 26 300 50 0.03 2 0.62 24 0.15 27 300 100 -0.02 0 0.71 27 0.16 27 350 50 0.00 1 0.65 26 0.15 28 400 25 -0.02 3 0.93 27 0.16 27 400 50 0.01 1 0.62 26 0.16 29 400 100 -0.01 0 0.62 25 0.15 27 500 10 0.02 11 4.36 28 0.15 26 500 25 -0.01 2 1.21 27 0.15 26 500 50 -0.01 0 0.81 26 0.15 26 500 100 -0.01 0 0.58 24 0.15 26 500 180 -0.01 0 0.67 27 0.15 26 550 100 -0.01 1 0.58 26 0.16 26 600 25 0.00 2 1.55 25 0.16 26 600 50 -0.01 2 0.81 30 0.15 27 600 100 0.02 1 0.59 26 0.16 27 700 25 -0.01 4 1.93 27 0.16 27 700 50 -0.03 3 0.95 28 0.15 25 700 100 -0.00 0 0.62 24 0.16 27 tvariable* 0.15 6 21.47 29 0.13 24 t500 100 0.11 7 3.15 22 0.14 27 T A B L E 6 . 2 : Orthonormal residuals statistics for various sets of conditional covariance parameters. For each set, the statistics qualifying the orthonormal residuals mean Q±, variance |1 — Q2I, and normality L have been averaged over the 31 time steps. The number of times each test has rejected the model is indicated. The chosen model is in bold. * The correlation range and conditional variance are estimated annually from the fluctuation data using the program AMLE2D. The time series of parameters can be found in Appendix D. t Results obtained by applying ordinary kriging on the fluctuation data. This is equivalent to considering the prior variance <r2_1(a;, x') to be constant in space. Chapter 6. Time-evolving DEM 64 6.4 Results: Changes in total ice volume, glacier profile, and ice distribution Block kriging estimates the value of the field averaged over a grid cell; the product of the ice thick-ness by the cell area gives the estimated ice volume in the cell. The total glacier volume is calculated by adding the contribution of all cells whose southwest corner falls within the glacier outline (Fig-ures B.6-B.8). The kriging standard deviation is taken as the uncertainty of the ice thickness; the integration of this uncertainty over the glacier area amounts to 20-25% of the total ice volume, de-pending on the year. Because the outlines are used as integration bounds, their uncertainty affects the ice volume calculations. This other source of uncertainty can be quantified by calculating the ice volume contained in the maximal and minimal outlines, as defined in Section 4.3. In Figure 6.6a, the error bars represent the volume uncertainty caused by the ice thickness standard deviation, and the triangles represent the cumulated effect of the ice thickness and outline uncertainty. The relative volume uncertainty associated with the definition of the outlines is small, less than 4%. As one can see in Figure 6.6a the total volume of ice did not change significantly over the last surge cycle. This suggests that the advance of the glacier was due to a redistribution of mass, rather than to a positive mass balance. The lack of precision of the estimate, however, makes it impossible to reject completely the possibility of a mass-balance contribution. The redistribution of mass can also be detected by following the glacier profile evolution. Figure 6.7 gives the centerline profile for 8 years between 1951 and 2005. The ice velocity inferred from flow poles drilled along the centerline is represented by arrows. The bulge becomes visible between 1974 and 1980 (see the complete set of profiles in Appendix E) . From 1981 to 2001 the bulge advanced at an average speed of ~ 2 5 m a _ 1 (Figure 6.5). After 1999, the advance slowed, and came nearly to a stop in 2001. The existence of a zone of strongly emergent flow at the boundary of active and stagnant ice, precursor of the formation of the bulge, was detected as early as 1970 (Collins, 1972). Figure 6.7 shows that emergent flow is present upstream from the bulge for most of the subsequent years. The advance of the terminus was accompanied by a significant thinning of the glacier. If the changes in geometry are mostly longitudinal, the area under the center profile is a good proxy for the total ice volume. The similarity between Figure 6.8e, which gives the total area under the profile (relative to 1981), and Figure 6.8f, which gives the total ice volume (also relative to 1981), shows that this is the case. Based on this assumption, the distribution of ice glacier-wide can be summarized by the distribution of ice along the centerline. Dividing the profiles in sections (as illustrated in Figure 6.7), one can capture the redistribution of ice more precisely than by integrating the ice volume. Figure 6.8a-d presents the contribution of each section to the total profile area, expressed as a percentage of the 1981 profile section. From 1980 to ca. 2000 mass is propagated downglacier, from section B to D. The behavior of the upper basin is harder to qualify because of the lack of data in this region. Between 1980 and 1995, the sparse data available did not signal any significant change in level; after the thorough surveying of the upper basin in 1995, however, it became clear that the total ice mass in the upper basin had decreased by more than 4% (Figures 6.6c-6.6e). This suggests that the transfer of mass from section B to D might have also drained part of the upper Chapter 6. Time-evolving DEM 65 2350 h 2100 -1 1 1 1 L _ 0 100 200 300 400 Distance (m) F I G U R E 6.5: The bulge is drawn at three year intervals from 1981 to 2005 (1981, 1984, 1987, 1990, 1993, 1996, 1999, 2002, 2005). The 2002 and 2005 profiles are nearly identical, and cannot be distinguished at this scale. Notice the progressive thinning of the terminus front. The 1987 profile (first dashed line) lacks data between the bulge and the tip of the apron. In the absence of data, kriging interpolates a smooth surface between these data points, which explains the unrealistic thickness of the apron for that year. basin reservoir. If this was the case, the apparent increase in total mass between 1985 and 1995 could be explained by the failure to detect the decrease in mass in the upper basin before 1995. The fact that the "sudden" 1995 mass loss balances almost perfectly the last few years of increase in the lower basin supports this hypothesis. However, given the large uncertainty bounds in the ice volume calculations, this remains speculative. Chapter 6. Time-evolving DEM 66 70 60 IS 50 CO vP —l fL-30 20 0.09 0.08 0.07 „ ~ E 0.05 * 0.04 0.03 1970 1980 1990 2000 F I G U R E 6.6: Total ice volume (a) and glacier area (b), 1951-2005, relative to the 1981 values (left axis) and in absolute values (right axis). The total ice volume is broken down in upper basin volume (c) and lower basin volume (d). The rapid decrease in upper basin elevation detected in 1995 returns the ice volume to near 1981 values (e). This suggest that the apparent increase in total ice volume between 1984 and 1995 could be due to a transfer of mass between the upper and lower basin. Because of the delay in detecting the resulting thinning of the upper basin, this transfer creates the illusion of a net gain in mass. Chapter 6. Time-evolving DEM 67 F I G U R E 6 .7 : Profile evolution and flow, 1951-2005. The arrows represent flow vectors calculated from poles in a 200 m corridor along the centerline. To avoid cluttering, I kept only one vector every 25 m, the closest to the centerline. Chapter 6. Time-evolving DEM 68 60 5) 50 40 40 Ui 30 -H-oh 5— 150 "D cn 1^ 100 ™ "~ 50 >~ 150 h 100 i = ~ 50 7 - r Hfi HHgfifiHHiiii l iHHII 1970 Hi 1980 1990 year 2000 F I G U R E 6 .8 : Area under centerline profile as a proxy for total ice volume. Because the centerline is densely surveyed, the ice thickness is well constrained in that area. To illustrate the changes in ice distribution that .accompanied the surge, the glacier profile is divided in four sections (SI, S2, S3, and S4, as defined in Figure 6.7), and the area under each section is integrated, shows that Chapter 6. Time-evolving DEM 69 6.5 Results: Driving stress Glacier flow is caused by the difference in lithostatic pressure between adjacent columns of ice. The stress that results from this imbalance is proportional to the ice thickness h and to the surface slope dZ Tdx = PiSh-^, (6.4) where pj is the density of ice, g the gravitational acceleration, h is the ice thickness, Z the ice surface elevation, and x represents either of the horizontal coordinates. The driving stress depends solely on the glacier geometry, and can be calculated from the bed topography and ice surface D E M s . Figure 6.10 gives the ice thickness, surface slope, and resulting driving stress for 1951, 1970, 1981, 1995, and 2005. (The complete time series can be found in Appendix D.) As one can see, the stress field is very heterogenous. High frequency variations are due to small-scale features in surface topography; variations in the ice thickness are more gradual, and affect the stress on longer length scales. Figure 6.9 presents time series of the mean thickness, median slope, and median stress. I use the median as a measure of central tendency for the slope and stress because it is more robust to extreme values. This precaution was not necessary for the ice thickness field. To eliminate the influence of the apron on these statistics, I have only considered cells with an ice thickness > 20 m. In the years preceding the onset of the surge (1969-1974), the mean ice thickness in the active area increased by ~15m, and the median driving stress by ~20kPa . The median surface slope decreased by ~ 1 ° . With the exception of the area near the terminus, the stress pattern did not change significantly between 1969 and 2005. High frequency features appear to have diffused over the years, but it is difficult to tell whether this is an artifact of the interpolation, or the result of a natural process. A natural explanation for the disappearance of these short-scale variations could be the gradual "healing" of crevasses opened during the previous surge. However, the fact that this diffusion operates at the same rate in the lower and upper basins, despite their very distinct flow rates, contradicts this hypothesis. M y opinion is that the diffusion is numerical rather than natural; the resolution of fine details in the photogrammetric D E M s is dissipated as the model is propagated in time, updated by more sparsely sampled survey measurements. 51 70 75 80 85 90 95 00 05 51 70 75 80 85 90 95 00 05 51 70 75 80 85 90 95 00 05 F I G U R E 6 . 9 : Time series of mean ice thickness, median slope, and median stress. To eliminate the effect of the apron on the statistics, only cells with an ice thickness > 20 m were considered. Chapter 6. Time-evolving DEM 7 0 F I G U R E 6 . 1 0 : Ice thickness (top), easting component of surface slope (middle), and easting component of driving stress (bottom). The background color in the slope map gives the absolute value of the median slope. Note that the color scheme for the stress map is not symmetrical. Chapter 6. Time-evolving DEM 71 The contact zone between active and stagnant ice is marked by an arc of high stress. Upstream from this perturbation, there generally exists a zone of relaxation, with low surface slope and low stress. This front can be detected as early as 1972, and becomes clearly defined in 1980. The cold-based apron forms a mechanical dam that impedes the advance of the active ice Clarke et al. (1986). Upstream from this obstruction, ice accumulates and the glacier thickens, eventually forming a bulge. As established by Clarke and Blake (1991), the bulge propagates faster than the ice, incorporating downstream stagnant ice by continuous deformation. The increased thickness and steepened slope of the bulge significantly enhance the creep, facilitating the advance of this kinematic wave. It is however improbable that the bulge would have been maintained if creep was the only source of flow. Sliding must have played an important role in bringing active ice forward to maintain the ice thickness and sustain the advance of the bulge. The increase in driving stress caused by the thickening of the glacier during the activation period (1969-1980) certainly contributed to the speed-up that occurred between 1974 and 1980. This is not to say, however, that enhanced creep is responsible for the surge. The fact that the glacier was still in a period of rapid flow between 1988 and 1996 (see Figure 4.4), despite the rapid thinning that occurred during that interval indicates that sliding must account for a significant part of the motion. Once activated, sliding can move ice downstream rapidly; if the progression of the ice is impeded, this will cause a thickening of the ice and an increase in driving stress. The observed increase in driving stress is therefore not necessarily a condition for surging, but may simply be a consequence of obstructed flow. The presence of the apron is therefore responsible for an interesting positive feedback between sliding and creep, which could have an significant impact on the dynamics of the surge. To investigate further the role of creep in the development of the surge, I need a flow model relating the driving stress to surface strain rates and velocities. 6.6 Results: Creep velocity and sliding Because glacier velocities are low, acceleration can be neglected and the second law of Newton reduces to an equilibrium of forces. The driving stress, generated by the gravitational force, is resisted by the basal drag, the lateral drag, and the longitudinal stress. When modeling glacier flow, a typical simplification is to assume that: (1) the basal drag is the only resistance to flow, and (2) the strain rate tensor has only one non-zero component, exz. Under this lamellar flow model, the vertical gradient in ice velocity u(z) is where A(T) is a temperature-dependent viscosity parameter, and n is the exponent of Glenn's flow law, the stress-strain relation for ice (Van der Veen, 1999, p. 104). A value of n = 3 is typical. Velocity at any depth can be obtained by integrating (6.5). The temperature at Trapridge Glacier varies with depth (Jarvis and Clarke, 1975). Given that most of the deformation occurs in the basal ice, which for most of Trapridge Glacier is at the melting point, one could take the value of A for temperate ice as a first approximation. In reality, the cold interior ice stiffens the glacier, and (6.5) Chapter 6. Time-evolving DEM 72 reduces the deformation rate. To obtain a flow parameter that accounts for this stiffening, I simplify the thermal structure of the glacier to a self-similar copy of a unique profile f(z/h), stretched or contracted to fit any ice thickness. Based on this assumption, the velocity at the glacier surface (z = h) is u{h) = 2AThT2x + us, (6.6) where us is the sliding velocity,4 and AT is the vertically-integrated flow parameter rh AT = j A(f(z/h))(l-^jndz. (6.7) For f(z/h), I used the profile from hole 11C (1988) from Clarke and Blake (1991, Figure 9). The temperature dependency of the flow parameter A follows the Arrhenius relation A = A0exp(-Q/RT), (6.8) where AQ is independent of temperature, Q is the activation energy for creep and R is the uni-versal gas constant (Paterson, 1994, p.86). Above — 1 0 ° C , Paterson (1994) suggest a value of <3= 139kJ/mol. The value of A(T(z)) at each point of the profile was extrapolated from the tabu-lated value A(0°C) = 6.8 x 1 0 " 1 5 s _ 1 k P a - ^ , (Paterson, 1994, p. 96) using the Arrhenius relation (6.8). The resulting temperature-integrated flow parameter is AT — 5.75 x 1 0 - 1 5 s _ 1 k P a - 3 . As expected, this is quite close to the value for a temperate glacier. Because it neglects longitudinal and transverse stresses, this lamellar flow model is also known as the shallow ice approximation. It is most suited for glaciers or ice caps that are much larger than they are thick. It is not a particularly good model for valley glaciers, where lateral and longitudinal stresses can be important. Nye (1965), in one of the first instances of computational glaciology, calculated numerical solutions to the flow of a glacier in a channel of rectangular, semi-circular, or parabolic cross-section. For different half-width to thickness ratios (W), he obtained a "shape factor," f, accounting for the transverse stress caused by the valley walls. Multiplying the shallow ice estimate for velocity by this shape factor reproduces the results from his numerical simulations. The cross section of Trapridge Glacier is somewhere between a rectangle and a parabola, with a half-width to thickness ratio of ~500m/60m = 8.3. The closest shape factors tabulated by Nye (1965) are for a rectangle with W = 3 (f = 0.884), and for a parabola with W = 4 (f=0.806). The effect of transverse stresses decreases rapidly as the half-width to thickness ratio is increased, and a ratio of 8.3 should have a shape factor very close to 1. Neglecting the stress caused by the valley walls is therefore a reasonable assumption in the case of Trapridge Glacier. Important longitudinal stresses occur wherever there are rapid variations in surface slope. At glacier scale, these short scale stresses cancel each other and only the average slope has an impact on the flow. To eliminate these high frequency variations, I used a two-dimensional median filter, with a square window of side 4 In this section, unless stated otherwise, I use the term "sliding" to describe any motion happening at the base. Mathematically, it represents the boundary condition of the velocity field, u(z = 0) = us. Physically, it is explained by sliding of the ice over the bed, deformation of the bed, or the cumulative effect of both. Chapter 6. Time-evolving DEM 73 51 70 75 80 85 90 95 00 05 51 70 75 80 85 90 95 00 05 F I G U R E 6 . 1 1 : Time series of the median creep rate in the upper basin (a) and in the lower basin (b). To eliminate the effect of the apron on the median, only cells with an ice thickness > 20 m were considered. = 2n(h) ~ 440m (22 cells). Of course, other sources of stress, such as protuberances in the bed, breaks in bed slope, throttling of the flow at the intersection of the upper and lower basin, 5 should reduce the creep rate. The shallow ice approximation is nevertheless expected to give a reasonable estimate. The time series of median creep rate in the lower and upper basin are presented in Figure 6.11. The difference between the calculated creep velocity and the observed flow gives a minimum value for the motion occurring at the bed. Creep velocity is parallel to the surface gradient; because this is not necessarily the case for the observed velocity, the sliding vectors may not be aligned with the flow. Figure 6.12 compares the predicted creep speed to the observed flow speed. The information obtained from flow poles is averaged over cells of 100 mx 100 m. In each cell, the percentage of motion in the flow direction due to sliding is indicated. Where the predicted creep velocity is larger than the observed flow, the sliding vectors point in the direction opposed to the flow. This is an obvious indication that, in these locations at least, the shallow ice approximation is not valid. The longitudinal or transverse stresses may be significant, or the ice may be much colder than estimated. Figure 6.14 summarizes the evolution of the sliding velocity in each flow zones (as defined in Figure 4.1). Figure 6.13 gives the percentage of motion due to sliding. As one can see, creep accounts for a negligible portion of the motion in the lower basin (< 10%). Over the years, the median creep velocity in the lower basin varies between 1 and 7 m a - 1 , following the changes in ice thickness (Figure 6.11). The apparent decrease of sliding velocities in the early 1980s (Figure 6.14), first in zone B, then in zone C , follows the propagation of the bulge. Because of the steep slope, the creep velocity predicted by the lamellar flow model at the bulge is quite large. However, this is certainly an overestimate. Such a rapid variation in level would generate significant longitudinal stresses, which are not accounted for by the shallow ice approximation. Including the effect of these stresses would reduce the creep velocity, and consequently increase the percentage of motion attributed to sliding. In the upper basin, which is considerably thicker (70-110 m), ice deformation plays a significant 5note that the half-width to thickness ratio in that region is closer to 4 Chapter 6. Time-evolving DEM 74 F I G U R E 6.12: Observed surface velocity (top), creep velocity (middle), and percentage of flow due to sliding (bottom). The ice speed deduced from flow poles is averaged over cells of 100mx 100m. Subtracting the predicted creep velocity from the observed flow yields sliding velocity vectors. The creep velocity is parallel to the surface topography gradient; since this is not necessarily the case for the observed velocity, the deduced sliding vectors may not be aligned with the flow. The bottom row gives the percentage of motion in the flow direction due to sliding. Where the predicted creep velocity is larger than the observed flow, this is a negative quantity (marked by crossed squared). Obviously, these are regions where the shallow ice approximation is not valid. Where the predicted creep vectors point in a direction opposed to the direction of flow, the percentage of motion due to sliding is greater than 100%. This is not an impossible occurrence. Chapter 6. Time-evolving DEM 75 F I G U R E 6 . 1 3 : Percentage of longitudinal velocity due to sliding in the upper basin (the union of zones UB1 and UB2), and in the lower basin (zones A, B. and C). The zones are defined in Figure 4.1. Chapter 6. Time-evolving DEM 76 F I G U R E 6 . 1 4 : Zonally-averaged sliding velocities. 1969-2005. The zones are defined in Figure 4.1. Chapter 6. Time-evolving DEM 77 role. The median creep velocity went from ~ 8 m a - 1 in 1969, peaked at 1 5 . 2 m a - 1 in 1989 and dropped rapidly to roughly 7 m a - 1 after the survey of the upper basin in 1995. As I discussed earlier, it is likely that the ice thickness (and therefore the creep velocity) was overestimated in the years 1984 to 1995, due to a lack of data. It is also probable that the ice in the upper basin is colder than that of the lower basin. The flow parameter was adjusted for a thermal profile typical of the lower basin; using it in the colder upper basin might overestimate the creep velocity. Irrespective of the creep estimate, the flow in the upper basin is much slower, indicating lower sliding velocities. Low sliding rates can be caused by the lack of water at the base. Only a small quantity of surface meltwater is produced each summer in the upper basin, and most of it is probably channeled to the lower basin before it can reach the bed. If only the meltwater produced by basal friction and geothermal heat is present to lubricate the bed, sliding will be hindered. If the basal ice in the upper basin is below freezing, sliding over the bed will be prohibited, leaving sediment deformation as the only possible source of basal motion. Although speculative, these two scenarios could explain the important difference in basal velocities between the upper and lower basin. Chapter 7 Conclusion 78 This thesis summarizes the geometrical and flow evolution of Trapridge Glacier, a surge-type glacier situated in the St. Elias Mountains, Yukon Territory, Canada. Over the period of study, 1951-2005, the glacier has been through a complete surge cycle. The previous surge occurred sometime between 1941 and 1949. It caused an advance of the terminus of ~1 km (an increase in length of about 40%) and left the glacier surface highly crevassed (Figure A.3). When field work began at Trapridge Glacier in 1969, the glacier was still recovering from that surge. The ice accumulated in the receiving area was stagnant and slowly melting away. Upstream, thicker "active ice" was flowing at ~10 m a - 1 . The installation of thermistor strings during the 1971 and 1972 field seasons established that the ice in the receiving area was below freezing temperature throughout its thickness, while the region of active ice was mostly warmed-based. The cold-based ice, acting as a mechanical dam, was impeding the flow of the warm-based ice. This contrast in mobility is responsible for the important morphological changes that occurred during the surge. 7.1 Flow The acceleration of flow during the surge is detected by repeated measurement of poles drilled into the ice. The regular surveying of these markers between 1969 and 2005 yielded 328 pole trajectories, from which annually averaged velocities can be calculated. Measurement blunders and typing errors, inescapable in this difficult environment, were eliminated by visual inspection and by a numerical algorithm designed to detect unrealistic pole displacements. From this expunged data set, 1339 annually averaged flow vectors are obtained for the periods 1969-1974 and 1981-2005. A regular grid of flow poles covers the central area of the lower basin; the upper basin and the marginal zones of the lower basin are sampled more sparsely. This is more than sufficient to detect annual variations in flow with a reasonable degree of spatial resolution. To characterize the changes in flow on a synoptic scale, I divided the glacier into six zones (moving downglacier: UB1, UB2, A , B , C , and D; see Figure 4.1) and considered the velocity vectors in each as an ensemble. At the beginning of the time series, the boundary between active and stagnant ice was situated in zone B . Between 1974 and 1980, the flow rate in the region of active ice increased dramatically. The me-dian velocity in zone A went from 1 5 . 6 m a - 1 to 3 8 . 6 m a - 1 , while zone B accelerated from 7 . 3 m a - 1 to 2 6 . 8 m a _ 1 . The flow peaked in 1984, at 4 1 . 8 m a _ 1 for zone A , 35.1 m a - 1 for zone B, and 33.2 for zone C , which by then was mostly comprised of active ice. As the flow accelerated between 1974 and 1980, a steep front formed at the boundary between active and stagnant ice. Over the following ten years, this bulge propagated downward, advancing faster than the ice, at an average speed of Chapter 7. Conclusion 79 30 m a . The cold-based ice was incorporated to the moving front by continuous deformation, and the limit of the warm-based ice migrated downglacier, although more slowly than the bulge (Clarke and Blake, 1991). At the same time, the stagnant ice downstream continued to melt, and the ter-minus retreated. Between 1981 and 1991, The bulge advanced ~300m. Sometime after 1991, the apron disappeared, its ice having been ablated or integrated into the bulge, and the active ice formed the terminus. After it peaked in 1984, the flow in the lower basin remained above 2 5 m a - 1 for 12 years, but was on a slowing trend. The slowdown followed a strangely regular 4-year pulse: 1-2 years of timid acceleration (2-3 m a - 1 ) , followed by 2-3 years of rapid deceleration (4-8 m a - 1 ) . The 1997-1999 acceleration was particularly vigorous, as the velocity in zone A went from 2 0 . 3 m a - 1 to 2 8 . 5 m a - 1 , an increase of 8 . 2 m a - 1 . After this last pulse, the glacier gradually slowed down to pre-surge velocities. The terminus advanced ~150m between 1991 and 1999, and only ~ 5 0 m from 1999 to present. In 2005, the lower basin was flowing at less than 8 . 5 m a - 1 . 7.2 Geometry Based on this flow history, I divide the 1969-2005 period in three parts: the activation phase (1969-1974), marking the tail end of the period of quiescence that followed the 1940s surge, the surge (1980-1999), and the waning phase (1999-2005). To quantify the geometrical changes that occurred during that period, I used air photos and optical survey measurements to construct 31 digital elevation models. First, D E M s for 1951, 1970, 1972, 1977, and 1981 were generated by stereographic analysis or air photos. These background models were then updated year-by-year by ground-based surveys data, using my implementation of a Bayesian kriging algorithm. This technique allows for the propagation and diffusion of information in time. It assumes the ice surface to be similar in shape from one year to the next, only shifted vertically. The interpolation obeys the data wherever they are available, and follows this assumption everywhere else. Uncertainty is attributed to each estimate, based on the distance to the nearest measurements and on the uncertainty of the prior model. When no data are available within a correlation range of a given location, the standard deviation of the ice thickness field at that location increases by 10 m. Based on the data, the correlation range was estimated to be 500 m. These statistical assumptions on the expected value and covariance structure of the ice thickness field were tested by analysis of orthonormal residuals, an advanced form of cross-validation. Tests on the mean of the residuals, which should be close to zero if the model can predict the data, were passed for each year at the 95% level. The variance of the residuals, on the other hand, was often too large or too small, and the residuals were not always normally distributed. The model was nevertheless judged reasonable in its simplicity; over 30 other statistical models were tested, but did not yield better results. The standard deviation of the estimates was used as a measure of uncertainty; as noted in Section 6.3, this is not strictly valid. The width of the statistical population required to predict the data is not only affected by the natural variability of the field, but also by the assumption of stationarity underlying the estimation. To pass the cross-validation tests, the standard deviation must compensate for the stringency of the statistical assumption. Using the standard deviation as uncertainty, although a Chapter 7. Conclusion 80 common practice, underestimates the accuracy of the results. Table 7.1 summarizes the geometrical characteristics of the glacier over the three phases of the surge. During the activation phase and the first half of the surge, the apron was melting and the terminus was retreating. When the active ice reached the terminus ca. 1991, the glacier entered a phase of advance. During the 8 last years of the surge phase, the snout moved ~125m downstream (15 .6ma - 1 ) . The advance during the 6 years of the waning phase was only of ~ 5 0 m , an average of 8 . 3 m a - 1 . The total area of the lower basin increased from 0.99 to 1.16km 2 between 1991 and 2005 (an increase of 17%). Despite these changes in area, the total volume of ice remained relatively constant from 1969 to 2005, at 0 . 1 4 ± 0 . 0 3 k m 3 . The apron was too thin for its disappearance to affect significantly the total ice volume. The spread of the glacier during the second half of the surge was accompanied by a corresponding decrease in mean thickness. Although the thinning of the upper basin was not detected until 1995, it may have started much earlier. The thickening of the lower basin in the 1980s must have required a transfer of mass from the upper basin. Assuming a constant ice volume, this should have caused a decrease of the ice level in that reservoir. A delay in the detection of this thinning could explain the apparent increase in total ice volume from 1984 to 1995. The fact that in 1995 the total ice volume suddenly returns to the 1981 level corroborate this interpretation, and suggest that the 1984-1995 positive mass balance is an artifact. From the surface slope and ice thickness fields derived from the D E M s of the surface and the bed, the stress driving the flow can be calculated. After filtering the surface slope to eliminate short scales variations, the stress and ice thickness fields are used to estimate creep velocities, based on a lamellar flow model. This model, based on the shallow ice approximation, supposes that the basal stress is the only resistance to flow, neglecting transverse and longitudinal stresses. Given the large half-width to thickness ratio of Trapridge Glacier (W~8.3), it is reasonable to neglect the transversal stresses caused by the valley walls. Stress due to irregularities of the bed, however, could be significant. Despite its simplicity, this model shows that creep accounts for a negligible portion of the flow in the lower basin, and that sliding or bed deformation are responsible for most of the motion. In the upper basin however, creep can account for more than half of the recorded motion, and sliding velocities are low. The creep velocities are possibly slightly overestimated, as they are calculated using a flow parameter adjusted for a temperature profile typical of the lower basin. It is likely that the upper basin is colder, thus composed of stiffer ice. Longitudinal stresses generated by several ice cliffs are also non-negligible. The fact that the calculated creep velocity surpasses the observed velocity in some steep areas of the upper basin indicates that, at least in these locations, longitudinal stresses restrict the flow significantly. The lower sliding velocity in the upper basin could be explained by the lack of meltwater at the base, or by a basal temperature below freezing point. 7.3 Discussion and outlook These modifications in flow and geometry, even if significant, are modest compared to the typical Alaskan surge behavior; they also occurred on a much longer time interval. The 1951 air photo Chapter 7. Conclusion 81 indicates the consequences of the previous surge of Trapridge Glacier, which clearly exceed the changes that occurred in the last 31 years. It is not surprising that the investigators, waiting for a surge of that amplitude, did not recognize the modification in flow that occurred between 1974 and 1980 as the onset of a new surge. Whether the instability responsible for this acceleration constitutes a surge, or a "failed attempt at surging" is largely a semantic question. What is interesting, however, is whether this instability could have degenerated into a 1940s-type catastrophic advance, and if so, what prevented it. Comparing the glacier in 1941 and 2004 shows a surprising resemblance. As reported by Sharp (1951), the terminus advanced by 200 feet between 1939 and 1941. This implies an advance of ~ 3 0 m a _ 1 , twice the rate measured between 1991 and 1999. So far, it has been thought that the 1940s surge had begun sometime after the 1941 expedition, and before the 1951 air photo. This overlooked comment of Sharp, however, suggests that the glacier was already in a state of enhanced flow prior to 1939. Using the 2004 outline as a proxy for the 1941 terminus position (which based on Figure 7.1 seems reasonable), the glacier advanced at least 1km between 1941 and 1951, implying advance rates greater than 1 0 0 m a - 1 . Was the flow instability responsible for the advance of the terminus between 1939 and 1941 similar to that responsible for the current surge? In which case, what would explain their different destinies? From the 1951 photogrammetric model, the post-surge volume was estimated to be 0.18 ± 0 . 0 5 k m 3 . The uncertainty of the 1951 D E M is quite large, but if taken at face value, this volume estimate represents 1.4 times the volume available for the current surge. Could this deficit in mass explain why the recent surge did not continue on the path of its predecessor? What sort of barrier was overcome between 1941 and 1951 to allow such rapid descent? Could the thickness of the glacier be a significant control? The flow, as we know, depends mainly on the rate of sliding and bed deformation, which are controlled by the subglacial environment. Given the absence of information on the basal conditions during the 1940s surge, it is difficult to answer these questions. Modeling the glacier is one way to investigate this further. To be able to predict the response of the glacier, a model would have to account for the heterogeneity of the basal conditions at Trapridge Glacier. It would also have to predict the transition between hydraulic states, and the response of the sediment bed to these different hydrological settings. Compared to these significant challenges, the production of digital elevation models describing the glacier is a modest contribution, but hopefully one that will inspire others to pursue the work. Even if it would not do well as a coffee-table book, the current surge remains an interesting flow anomaly. It has the unique characteristic of having been observed as carefully from underneath as from the surface. The year-round recording of subglacial properties at Trapridge Glacier started in 1989, and will end in 2007. The onset of the surge was missed, but the transition from the surge to the waning phase occurred when the full instrument array was deployed. The rapid increase in velocity between 1997 and 1999 should be closely observed, as it marks the end of the surge phase. Could a significant disruption of the subglacial setting have occurred over that period of acceleration? Can any other synoptic measurements correlate with the 4-year cycles in surface velocity? These features of the flow history are leads for the investigation of the links between the highly heterogeneous subglacial properties and the integrated response of the glacier. In addition to Chapter 7. Conclusion 82 F I G U R E 7.1: Trapridge Glacier, 1941 (left, photo by R.P. Sharp) and 2004 (right, photo by A . J . Schaeffer). synoptic history, the D E M s and flow maps presented in this thesis offer a second interesting avenue for the characterization of the subglacial environment. Using the flow maps and the surface and bed D E M s , one can obtain by inversion the basal stress and sliding velocity fields. This spatial information on the basal environment can then be compared to the records of instruments inserted at the bed. The time series of subglacial instruments are complex and very different from one site the the next. To be interpreted, it is very helpful to know what to look for. This thesis provided a first key to their interpretation: a time series of global events, specifying distinct dynamical phases and important periods of transition. By defining regions of high stress and other "hot-spots," the inversion of the flow field would provide a spatial structure to facilitate the interpretation of subglacial records. Knowing when and where to look for shifts in the subglacial environment, the basal conditions controlling the flow and geometry of the glacier could be determined. Chapter 7. Conclusion 83 Activation Surge Waning Now 1951 1969-1972 1980-1999 1999-2005 2005 Median velocity U B l , m a - 1 - 8.2 11.7 11.4 10.7 Median velocity U B 2 , m a - 1 - 11.5 18.3 14.9 9.9 Median velocity L B , m a - 1 - 10.6 29.1 16.5 8.1 Total area, k m 2 3.08 ± 0 . 0 7 2.62 ± 0 . 1 6 1.98 ± 0 . 1 1 2.10 ± 0 . 0 7 2.11 ± 0 . 0 5 Total volume, k m 3 0.18 ± 0 . 0 5 0.14 ± 0 . 0 3 0.14 ± 0 . 0 3 0.14 ± 0 . 0 3 0.13 ± 0 . 0 3 Mean ice thickness, m 62 67 74 65 64 Median slope, 0 11 12 10 10 10 Median driving stress, kPa 98 107 113 100 99 • Median creep ve loc i ty ,ma - 1 1.9 5.7 6.4 3.4 3.4 L O W E R B A S I N Total area, k m 2 2.12 ± 0 . 0 7 1.67 ± 0 . 1 6 1.03 ± 0 . 1 1 1.16 ± 0 . 0 7 1.16 ± 0 . 0 5 Total volume, 1cm3 0.10 ± 0 . 0 3 0.06 ± 0 . 0 2 0.06 ± 0 . 0 1 0.06 ± 0 . 0 1 0.06 ± 0 . 0 1 Mean ice thickness, m 53 49 64 55 54 Median driving stress, kPa 89 82 106 89 88 Median creep velocity, m a - 1 1.2 1.9 4.8 2.2 1.9 U P P E R . B A S I N Total area — 0.94 k m 2 (fixed) Total volume, k m 3 0.08 ± 0 . 0 1 0.08 ± 0 . 0 1 0.08 ± 0 . 0 2 0.07 ± 0 . 0 2 0.07 ± 0 . 0 2 Mean ice thickness, m 82 87 83 76 77 Median driving stress, kPa 132 140 123 116 113 Median creep velocity, m a - 1 9.2 9.7 10.4 7.1 6.8 T A B L E 7.1: Summary of geometrical changes, 1951-2005. 84 Bibliography Bamber, J. L., Vaughan, D. G. and Joughin, I. (2000), 'Widespread complex flow in the interior of the antarctic ice sheet', Science 287(5456), 1248-1250. Bjornsson, H. (1998), 'Hydrological characteristics of the drainage system beneath a surging glacier', Nature 395(6704), 771-774. Blake, E., Clarke, G. K. C. and Gerin, M . C. (1992), 'Tools for examining subglacial bed deforma-tion', Journal of Glaciology 38(130), 388-396. Budd, W. (1975), 'A first simple model for periodically self-surging glaciers', Journal of Glaciology 14, 3-21. Clarke, G. K. C. (1976), 'Thermal regulation of glacier surging', Journal of Glaciology (16), 231-250. Clarke, G. K. C. (1991), 'Length, width and slope influences on glacier surging', Journal of Glaciol-ogy 37(126), 236-246. Clarke, G. K. C. and Blake, E. W. (1991), 'Geometric and thermal evolution of a surge-type glacier in its quiescent state - Trapridge Glacier, Yukon Territory, Canada, 1969-89', Journal of Glaciology 37(125), 158-169. Clarke, G. K. C , Collins, S. and Thomson, D. (1984), 'Flow, thermal structure, and subglacial conditions of a surge-type glacier', Canadian Journal of Earth Sciences 21(2), 232-240. Clarke, G. K. C , Schmok, J. P., Ommanney, C. S. L. and Collins, S. G. (1986), 'Characteristics of surge-type glaciers', Journal of Geophysical Research - Solid Earth and Planets 91(B7), 7165-7180. Collins, S. G. (1972), 'Survey of the Rusty Glacier area, Yukon Territory, Canada, 1967-1970', Journal of Glaciology 11(62), 235-253. Dowdeswell, J. A. , Hamilton, G. S. and Hagen, J. O. (1991), 'The duration of the active phase on surge-type glaciers - contrasts between Svalbard and other regions', Journal of Glaciology 37(127), 388-400. Fischer, U . H. and Clarke, G. K. C. (2001), 'Review of subglacial hydro-mechanical coupling: Trapridge Glacier, Yukon Territory, Canada', Quaternary International 86, 29-43. Bibliography 85 Flowers, G . E . (2000), A multicomponent coupled model of glacier hydrology, P h D thesis, The University of British Columbia. Flowers, G . E . and Clarke, G . K . C . (1999), 'Surface and bed topography of Trapridge Glacier, Yukon Territory, Canada: digital elevation models and derived hydraulic geometry', Journal of Glaciology 45(149), 165-174. Flowers, G . E . and Clarke, G . K . C . (2002), 'A multicomponent coupled model of glacier hydrology-2. application to Trapridge Glacier, Yukon, Canada', Journal of Geophysical Research - Solid Earth 107(B11). Fowler, A . C , Murray, T . and Ng, F . S. L . (2001), 'Thermally controlled glacier surging', Journal of Glaciology 47(159), 527-538. Hamilton, G . S. and Dowdeswell, J . A . (1996), 'Controls on glacier surging in Svalbard', Journal of Glaciology 42(140), 157-168. Harrison, W . and Post, A . (2003), How much do we really know about glacier surging?, in 'Annals of Glaciology, Vol 36', Vol. 36 of Annals of Glaciology, International Glaciological Society, Cambridge, pp. 1-6. Jarvis, J . and Clarke, G . (1975), 'The thermal regime of Trapridge Glacier and its relevance to glacier surging', Journal of Glaciology 14(71), 235-250. Jiskoot, H . , Boyle, P. and Murray, T . (1998), 'The incidence of glacier surging in Svalbard: E v i -dence from multivariate statistics', Computers & Geosciences 24(4), 387-399. Jiskoot, H . , Murray, T . and Boyle, P. (2000), 'Controls on the distribution of surge-type glaciers in Svalbard', Journal of Glaciology 46(154), 412-422. Kamb, B . (1987), 'Glacier surge mechanism based on linked cavity configuration of the basal water conduit system', Journal of Geophysical Research - Solid Earth and Planets 92(B9), 9083-9100. Kamb, B. , Raymond, C . F . , Harrison, W . D. , Engelhardt, H . , Echelmeyer, K . A . , Humphrey, N . , Brugman, M . M . and Pfeffer, T . (1985), 'Glacier surge mechanism - 1982-1983 surge of Variegated Glacier, Alaska', Science 227(4686), 469-479. Kavanaugh, J . L . and Clarke, G . K . C. (2001), 'Abrupt glacier motion and reorganization of basal shear stress following the establishment of a connected drainage system', Journal of Glaciology 47(158), 472-480. Kitanidis, P. K . (1991), 'Orthonormal residuals in geostatistics - model criticism and parameter-estimation', Mathematical Geology 23(5), 741-758. Kitanidis, P. K . (1997), Introduction to Geostatistics, Cambridge University Press, Cambridge. Kitanidis, P. K . (2005), Applied Stochastic Inverse Methods, Course Notes. Bibliography 86 Lingle, C . S. and Fatland, D. R. (2003), 'Does englacial water storage drive temperate glacier surges?', Annals of Glaciology 36, 14-20. Marcotte, D . (1991), 'Cokriging with Matlab', Computers and Geosciences 17(9), 1265-1280. Meier, M . F . and Post, A . (1969), 'What are glacier surges', Canadian Journal of Earth Sciences 6(4P2), 807 - 823. Murray, T . and Clarke, G . K . C. (1995), 'Black-box modeling of the subglacial water-system', Journal of Geophysical Research - Solid Earth 100(B6), 10231-10245. Murray, T . , Dowdeswell, J . A . , Drewry, D. J . and Frearson, I. (1998), 'Geometric evolution and ice dynamics during a surge of Bakaninbreen, Svalbard', Journal of Glaciology 44(147), 263-272. Myers, D. E . (1984), Cokriging: new developments, in G . Verly et al., eds, 'Geostatistics for natural resources characterization', Vol. v. 122, No 1, N A T O - A S I , Dordrecht, pp. 295-305. Nye, J . (1965), 'The flow of a glacier in a channel of rectangular, elliptic or parabolic cross-section', Journal of Glaciology 5(41), 661—690. Olea, R. A . (1999), Geostatistics for Engineers and Scientists, Kluwer, Dordrecht. Omre, H . (1987), 'Bayesian kriging-merging observations and qualified guesses in kriging', Mathe-matical Geology 19(1), 25-39. Omre, H . and Halvorsen, K . B. (1989), 'The Bayesian bridge between simple and universal kriging', Mathematical Geology 21(7), 767-786. Pardo-Iguzquiza, E . (1998), 'Maximum likelihood inference of spatial covariance parameters of soil properties', Soil Science 163(3), 212-219. Pardo-Iguzquiza, E . and Dowd, P. A . (1997), ' A M L E 3 D : A computer program for the inference of spatial covariance parameters by approximate maximum likelihood estimation', Computers & Geosciences 23(7), 793-805. Paterson, W . (1994), The Physics of Glaciers, 3rd edn, Pergamon, Oxford. Post, A . (1969), 'Distribution of surging glaciers in North America', Journal of Glaciology 8(53), 229-240. Raymond, C . F . (1987), 'How do glaciers surge - a review', Journal of Geophysical Research -. Solid Earth and Planets 92(B9), 9121-9134. Rothlisberger, H . (1972), 'Water pressure in intra- and subglacial channels', Journal of Glaciology 11, 177-203. Sharp, M . (1988), 'Surging glaciers — behavior and mechanisms', Progress in Physical Geography 12(3), 349-370. Bibliography 87 Sharp, R. (1943), 'Geology of the Wolfe Creek area, St. Elias Range, Yukon Territory, Canada', Geological Society of America Bulletin 54, 182-189. Sharp, R. (1951), 'The glacial history of Wolfe Creek, St. Elias Range, Canada', Journal of Geology 59, 97-117. Stone, D. B. and Clarke, G . K . C . (1993), 'Estimation of subglacial hydraulic-properties from induced changes in basal water-pressure - a theoretical framework for borehole-response tests', Journal of Glaciology 39(132), 327-340. Stone, D . B. and Clarke, G . K . C . (1996), 'In situ measurements of basal water quality and pressure as an indicator of the character of subglacial drainage systems', Hydrological Processes 10(4), 615-628. Van der Veen, C . (1999), Fundamentals of Glacier Dynamics, A. A. Balkema, Rotterdam. Wilbur, S. (1988), Surging versus non-surging glaciers: a comparison using morphology and bal-ance, P h D thesis, University of Alaska. 88 Appendix A P h o t o g r a m m e t r i c m o d e l s o f t h e R u s t y C r e e k b a s i n : t e c h n i c a l i s s u e s A . l Procedures for the orientation of air photos and collection of elevation data, by E . Saczuk, ISYS Consultants A.1.1 Orientation of air photos Interior orientation A n interior orientation must be performed on each photo using parameters and values found in the camera calibration report corresponding to each photo. This is to minimize possible distortions introduced by optical imperfections of the camera lens and to reduce the atmospheric and earth curvature effects on the photos. The first step is to create a virtual camera using the following parameters found in the camera calibration report: • Calibrated focal length • Principal point of best symmetry • Principal point of auto collimation • F i lm width and length • Fiducial coordinates • Lens distortion values • Orientation of camera • View geometry The interior orientation is then carried out by measuring the coordinates of the fiducial marks on each photo in the appropriate sequence. Figure A . l below shows details and typical layout of fiducials used in two of the commonly used aerial mapping cameras, namely the Wi ld R C 20/30 and the Zeiss R M K . The orientation is complete when the maximum residual values between the fiducial coordinates as computed based on lens and film distortion and those actually measured on the air photo by the operator, fall below a predefined limit. The predefined limit is typically set according to the scan resolution (14 microns for this project) so that the maximum residual is smaller than one pixel (i.e. 14 microns) and this is also the case here. Appendix A. Photogrammetric models of the Rusty Creek basin: technical issues 89 Camera Configurat ions for Strip Wizard tmm Imp r^*jg F I G U R E A . l : Typical fiducial detail aud layout with the numbers corresponding to the sequence in which the fiducials are to be measured. Re la t ive or ientat ion Following a successful interior orientation, a relative orientation is performed to minimize y-parallax in order to recreate the relative geometry of each photo at the time its exposure. This allows an air photo pair to be viewed in stereo. A pattern of relative control points is displayed (typically 3 to 5 on each photo) and the operator measures each point on both the left and right photos, carefully picking the same location in order to maintain the desired level of accuracy. Through the process of measuring the points, the software applies a mathematical transformation (typically polynomial) to each photo in order to recreate their relative alignment at the time of exposure. Once the maximum residual value falls below a predefined limit (again based on scan resolution, as with interior orientations), the relative orientation is complete and the stereo model can be used to accurately measure the position and elevation of objects on the ground. A b s o l u t e or ientat ion The goal of an absolute orientation is to register each stereo model to a ground coordinate system and datum. This is carried out by measuring locations on the model which correspond to ground control points at which the actual coordinates and elevations are known either through G P S or traditional surveying techniques or derived from large-scale maps. A stereo model is considered to be absolutely orientated when the maximum residual in both the horizontal and vertical position falls below a predefined limit. As with relative orientations, the software applies a mathematical transformation (linear, quadratic or polynomial) in order to fit the stereo model to a ground coordinate system. Appendix A. Photogrammetric models of the Rusty Creek basin: technical issues 90 Point ID Easting Northing Elevation m m m L N 536085.9 6788770 1 2348 8 P T 535904.3 6788657 2 2367 8 T N 534745.1 6788305 7 2545 8 A N 534776.7 6788250 2 2536 3 J L 534890.4 6788194 1 2510 2 U R 535604.2 6787281 3 2464 2 L R 536084.5 6787297 7 2357 2 M S 536160.1 6787794 0 2261 2 M K 536095.9 6788231 3 2268 2 R L 537868.3 6788033 3 2126 2 R M 537754.7 6787701 4 2163 8 R U 537712.6 6787546 3 2187 3 T A B L E A . l : Ground control points used in the absolute orientation of the air photos of Trapridge Glacier, using the WGS 84 ellipsoid and N A D 83 datum. Typically, the horizontal limit is determined by multiplying the scan resolution by the scale of the photo. For example, the length of one side of a pixel on a 1:50,000 scale photo scanned at 14 microns is equal to 0.7 m. The vertical limit is often set at 2 to 3 times the horizontal limit to account for the lower precision associated with horizontal measurements. The control points used to orient the stereo models are listed in Table A . l below, and mapped in Figure A.2. The residual values for the interior, relative and absolute orientations for each photo and model are found in the following files which are located on the accompanying D V D : 1951 Report.txt 1970 Report.txt 1972 Report.txt 1977 Report.txt 1978 Report.txt 1981 Report.txt . These points are visible on the 1981 air photos because L-shaped aerial targets were placed on the ground prior to the acquisition of the photos. It was necessary to compare the locations of the control points on photos taken in other years with the locations of the points on the 1981 photos in order to ensure that the appropriate control point was being measured. In most cases topographic features helped to positively identify the control points however shadows, snow cover and rapidly changing drainage patterns hindered the identification process and directly affected the precision with which some of the control points were measured. Appendix A. Photogrammetric models of the Rusty Creek basin: technical issues 91 A.1.2 DEM surface collection Elevation points are collected from each stereo model using a regularly-spaced grid of points, each of which is placed on the ground surface and measured by the operator. The vertical accuracy of each point is affected by a number of factors associated with the subjective decision of the operator who ultimately decides if the measuring mark is in fact on the ground surface. The ability of the operator to accurately place the measuring mark on the surface depends on lighting and contrast variations in the photo as well the type of terrain. The operator may skip points in featureless areas devoid of shadows with flat lighting and poor contrast or where the viewing geometry prohibits the positive identification of the ground surface due to excessive overlap of air photos or the obstruction of the ground by dense trees for example. In addition to errors in the vertical position introduced through the processes of orientation, operator experience also affects the uncertainty of terrain heights. In this case, the maximum vertical uncertainty of measured points as a result of terrain conditions and operator experience was deemed to be in the order of ± 2 0 m. Two separate elevation models were collected from each stereo model; a fine grid and coarse grid. The fine grid was spaced at 20 m intervals between points and encompassed only the Trapridge Glacier (533820 m E - 537420 m E , 6787320 m N - 6788580 mN). The coarse grid was spaced at 100 m intervals between points and encompassed an area bounded by 533000 m E - 539000 m E , 6784000 m N - 6789900 mN. Once all of the elevation points are measured or skipped due to poor viewing geometry or inability to locate the ground surface, the point file can be imported into a software package capable of creating a T I N , grid, D E M or other elevation surface. L i m i t a t i o n s The surface model derived from the 1951 imagery is deemed to be less accurate than the other surface models mainly due to the poor stereo geometry between the stereo pair. Instead of the standard 60% sidelap between adjacent air photos, the 1951 photos overlapped by over 80% which reduced the angle between light rays from a common point on both photos resulting in a larger area of focus thus degrading the precision with which the measuring mark could be placed on the ground, especially in shadowless, featureless areas like glacier surfaces. Although an additional two sets of air photos were available (1967 and 1978), they were not used in the analysis for the following reasons. The 1967 air photos did not provide sufficient stereo overlap of the Trapridge glacier and therefore could not be used to setup a stereo model and generate the elevation models. The 1978 air photos could not be absolutely oriented with sufficient precision because more than half of the ground control points fell outside of the air photo coverage or in areas covered by snow. This greatly reduced the horizontal and vertical accuracy of the stereo model and a decision was made to eliminate this dataset due to these problems which could not be easily rectified. The software package used to perform the orientations (i.e. Integraph Softcopy Workstation) has no option to output the oriented and transformed air photos. The software performs a mathematical transformation on the photos which does not result in a transformed output image. Thus the Appendix A. Photogrammetric models of the Rusty Creek basin: technical issues 92 photos used in the drape operation were not corrected for lens distortions, atmospheric and ground curvature effects or film distortion and were not oriented. The raw air photos were georeferenced in Idrisi GIS software by first assigning U T M coordinates based on the N A D 83 datum to each photo then registering each corner to an actual U T M coor-dinate derived from the oriented and corrected air photos. A n alternate option is to assign corner coordinates to the raw photos derived from topographic maps, although given the type of terrain which lacks easily recognizable features that can be used as control points, this approach would most likely lead to less reliable georeferenced image. A.2.1 Impo r t i ng and g r idd ing D E M s us ing MATLAB The elevation points and digitized outlines where handed to me as a dxf-file. These files were im-ported into Matlab using a reader published on the Mathworks File Exchange website: dxf 2coord.m, by Lukas Wischounig. The D E M data points are only approximately gridded, I therefore used an interpolation routine (griddata.m) to evaluate the elevation at the exact grid-node locations. A l l Trapridge survey data were measured relative to a set of control points, which position were established by astronomical measurements by Sam Collins in the early days of the study. These measurements were made relative to the datum of the time, NAD27. Hence the survey data are also in NAD27. For the sake of compatibility, the D E M s and all the subsequent work presented in this thesis are defined in NAD83, the current datum. Hence the control points had to be transformed from NAD27 to NAD83. Table A . l gives the NAD83 coordinates, while Table A.2 gives the original values in the NAD27 datum. I used these control points to obtain by regression a linear transformation to convert Trapridge data to NAD83; the official transformation is of course more complex, but the linear approximation suffice given the small range of the data: A . 2 Other issues regarding the Trapridge DEMs, BY Tom-Pierre Frappe-Seneclauze XNADSS = 1.00002429558483 x X N A D 2 7 - 142.456651242439 YNAD83 = 0.99976140050240 x Y N A D 2 7 + 1789.70245712966 ZNADS3 = 1.00099751166935 x Z N A D 2 7 - 0.983565105288996 (A.1) (A.2) (A.3) Appendix A. Photogrammetric models of the Rusty Creek basin: technical issues 93 Point ID Easting Northing Elevation m m in A N 534905 801 6788080 175 2534 732 F X 536980 700 6786371 725 2409 937 H K 538361 800 6787497 165 2358 289 J L 535019 567 6788024 055 2508 676 L N 536216 223 6788600 271 2347 380 L R 536213 736 6787127 451 2355 551 M K 536225 290 6788060 928 2266 927 M S 536289 246 6787624 171 2260 401 O C 537836 326 6787294 571 2195 410 P O 534954 543 6788008 901 2521 955 P T 536033 543 6788487 061 2367 047 R L 537997 564. 6787863 249 2124 926 R M 537884 008 6787531 284 2162 526 R U 537841 879 6787376 180 2186 015 T N 534874 626 6788135 605 2543 952 U R 535734 219 6787110 559 2462 622 G K 536264 344 6788107 009 2263 156 BS 536383 576 6787673 685 2226 637 CS 535527 684 6788225 446 2377 558 T A B L E A . 2 : U T M coordinates of ground control points around Trapridge Glacier, in NAD27. The lower basin is surveyed from L N , and the upper from T N . UR is used as.reference point to close the rounds. MS and M K have also been used as survey sites. Appendix A. Photogrammetric models of the Rusty Creek basin: technical issues 91 F I G U R E A . 2 : Location of control points used to reference the D E M s . Appendix A. Photogrammetric models of the Rusty Creek basin: technical issues 95 F I G U R E A . 3 : Aerial photograph of Rusty Creek basin, 1951. The origin of the frames (southwest corner) is located at 533000m East, 6874000 m North in U T M coordinates. Roll number: A13136, photo number: 44-45, scale: 70000. July 9th 1951. Appendix A. Photogrammetric models of the Rusty Creek basin: technical issues 96 F I G U R E A . 4 : Georeferenced aerial photograph of Rusty Creek basin, 1970. The origin of the frames (southwest corner) is located at 533000m East, 6874000 m North in U T M coordinates. Roll number: A21523, photo number: 123-124, scale: 60000. August 5th 1970. F I G U R E A . 5 : Georeferenced aerial photograph of Rusty Creek basin, 1972. The origin of the frames (southwest corner) is located at 533000m East, 6874000m North in U T M coordinates. Roll number: A22998, photo number: 66 67, scale: 80000. August 11th 1972. F I G U R E A . 6 : Georeferenced aerial photograph of Rusty Creek basin, 1977. The origin of the frames (southwest corner) is located at 533000m East, 6874000m North in U T M coordinates. Roll number: A24762, photo number: 130-131, scale: 70000. August 4th 1977. F I G U R E A . 7 : Georeferenced aerial photograph of Rusty Creek basin, 1981. The origin of the frames (southwest corner) is located at 533000m East, 6874000m North in U T M coordinates. Roll number: A25841, photo number: 35-36, scale: 80000. August 24th 1981. Appendix A. Photogrammetric models of the Rusty Creek basin: technical issues 100 XT3 Other DEM results To investigate further the topographic changes that cause these zones of high variance I mapped for each year the elevation change relative to the previous available model (Figure A.8). I also compared each D E M to the population average D E M (Figure A.9). F I G U R E A . 8 : Elevation difference from one Rusty Creek basin D E M to the next. F I G U R E A . 9 : Elevation difference between each Rusty Creek basin D E M and the mean D E M . 102 Appendix B Survey data import and filtering summary The goal of this appendix is twofold: first, to describe the algorithm used to filter the survey data; second, to introduce the various Matlab functions that I created to handle this plethora. Hence a good part of this chapter is with naive optimism dedicated to the lost soul who might one day wish to reuse this data set. To the busy reader, I suggest going directly to Section B . 2 , where the filtering algorithm is presented with a minimal amount of information on its implementation. B.l Accessing survey data 1) Reading survey log-files: S = read_log(f ile) • outputs a structure array with one entry per line of the log-file. Each element of the "survey structure" contains information relating to one measurement, in the fields: — tsdn: time of the shot, in serial date number format — X : Easting position (m) — Y : Northing position (m) — Z : Elevation of the lower tip of the flow pole (an ice-particle) — E l : Elevation of the prism — E2 : Elevation of the prism - pole height (measured at the beginning of the season), ~elevation of the glacier surface. — reset ('Y'/'N'/'S'): if 'Y' the pole was redrilled, and continuity of the time series was lost. If 'S' the pole was redrilled but the continuity the time series preserved (see testreset.m) — height: position of the prism relative to the ground as measured at the beginning of the season. — offset: position of the prism relative to the top of the pole. — ID : Target ID. — line: line number if source file. This identifies the shot uniquely. — ref 1: Reference point 1. — ref 2 : Reference point 2 . — surveyor: indeed. Appendix B. Survey data import and filtering summary 103 — notetaker: ah, the good ole days. — note: all that's left at the end of the line. Notes. — source: source file. 2) Finding specific types of measurements: [ndx . l i s t ] = svy_f i l t e r ( I D , type) • From a given list of target-IDs svy_f i l t e r selects the targets of a certain type and return the corresponding indexes (ndx). See the paragraph on s p l i t . m for details on the types available. The IDs can also be selected from a list file passed on as argument. S(ndx) returns the corresponding data-structures. 3) Sorting survey data: s p l i t ( f i l e .options) • Using svy_f i l t e r .m, splits the lines in a log-files in various functional categories based on their ID string. Shots are assigned to one of these directories (in order of preference): — grid : R X X C Y Y pattern — lines : either C L X X X , N L X X X , S L X X X pattern or from the 'radar' or "oldlines" lists defined in s v y . f i l t e r . — holes : H Y Y pattern — loggers : # Y Y Y Y Y pattern — upperbasin : flow poles in the u.b. From a list in svy_f i l t e r . — oldgrid : flow poles that do not obey R X X C Y Y pattern. From list. — junk : based on an explicit list in svy_f i l t e r . — singles : shots with unique IDs. — others : whatever is left... • The "grid", "upperbasin" and "oldgrid" lists are combined to make an "allpoles" log-file. These data-points are the ones used for flow calculation. • For each category, a lst-file lists the target names, the number of times they were shots and their life-span. The lst-file can be sorted in various ways using the opt ion parameter or the function resort .m. • A multi-column, printer-friendly, sum-file gives the ID of all targets in the directory. B.2 Filtering survey data B.2.1 By visual inspection The data are first filtered by visual inspection using p u r g e a l l .m • purgeall(fi le,type,dirout) Appendix B. Survey data import and filtering summary 104 — type is either "poles" or "lines": poles: are sorted by ID names, ordered by time. l i n e s : generally of the form A B C X X X , they are sorted by the letters A B C , then ordered by the numbers X X X , separately for each year. This sorting organizes the targets so that each of them can be easily compared to their neighbor (for lines) or previous shots (for poles). Each line or time series is plotted separately and the user can navigate the sequence and eliminate outliers. This is done using the function purge.m. • purge(Ssubset ,d irout , pnow_, pole_ , plotoptns) : Plots the data section passed as argument (as a survey-structure) using p lot_svy .m. The user is invited to delete targets (press'd') or mark them as interesting or suspicious (press'm'). p lot_svy returns the corresponding shots in the global structures Mark and Delete. If in "pole_" mode, purge calls check4reset .m to insure that the reset flag of a deleted reset-shot is passed to the next non-deleted measurement of that target. That way the information that the pole was reset is not lost when the faulty data point is taken away. When pnow_ is set high, the Mark and Delete structure arrays are written to a mrk- and del-file into the d i rout directory, which is created the first time the function is called. The deleted lines are taken away from a copy of the original log-file in d irout . The original f i l e is not modified, however, if d irout is explicitely set to the directory of f i l e , the log-file is edited directly, and the marked and deleted lines are appended to the existing mrk- and del-files. This allows to start the hand-filtering job one day, and go back to it later! B .2.2 By flow analysis Incorrect measurement of the position of a pole, used to track the trajectory of an ice-particle, will yield an incorrect flow vector. "Unreasonable" flow-vectors can therefore be used to detect outliers. This a posteriori quality assessment was implemented as a part of the routine used to calculate the flow-field. It is described in the following section. B.3 Filter and segment flow-pole time series: makeallcontinuous .m (1) makeallcontinuous .m calls makecontinuous .m for each target in a log-file, makecontinuous.m separates the target's time series into continuous segments (i.e. sequences not interrupted by the melting of the pole) and filter each segment individually, looking for outliers by calculating ice-velocities. • [ A l l , D e l e t e ] = makea l l cont inuous ( f i l e , l eve l ) — A l l : A structure array of length equal to the number of targets. It has two fields: ID is the target name. Data is a structure array outputted by makecontinuous. Each element Appendix B. Survey data import and filtering summary 105 of Data contains the information on one continuous sequence of the given target: * Data( i ) .Shots Is a survey-structure of length equal to the number of shots in the continuous sequence. * Data( i ) .per iod contains the serial date number of the first and last shots of the sequence * Data( i ) .year contains the year of the first and last shots of the sequence. — Delete: Survey-structure array of the shots that were rejected by the filter. • make continuous: separates the time series for each pole into continuous segments. To calcu-late the flow that occurred between two measurements of a pole, I need to insure that the pole was not reset meanwhile. In the log-file, a reset status is specified with each measurement; either 'N' or 'Y' . A continuous segment is from one reset to the next, makecontinuous will declare a reset "synchronous" if the shot previous to the reset is less than maxgap (= 3) days away. The reset status is then set to 'S'. S-resets are not considered a rupture of continuity, since the location of the ice-particle the pole was attached too is known relative to the new position of the pole. When calculatespeed is called, it verifies if a S-reset occurred between the two end-points of the trajectory and subtract the artificial movement due to the reset. (2) Test shots in each continuous sequence: test_shot.m and i sposs ib le f low.m — [yn_,data,mssg,summary] = t e s t _ s h o t ( n d x , d a t a . l e v e l , f a c t o r ) : This function is the core of this filtering algorithm. It compares each shot to be tested (Data(ndx)) to selected points in the same continuous time series, and calculates the flow that occurred for each combinations. The function i s p o s s i b l e f low (see below) is called for each of these pairs to decide if the amount of flow that occurred between the two end-points of the trajectory is reasonable. From this series of binary tests a conclusion is reached on the quality of the tested shot. Annually averaged flow is used for testing to smooth short term flow instability and measurement errors. Shots less that dtmin = 10 days away are definitively not used as end-points, unless no other shots are available. Shots more that one year away are not used either, since the averaging over a long time period reduces the power of the test. The preferred comparison points are a year away. The flow from the tested shot to every other point in this time-window is tested using i s p o s s i b l e f low. The many test are then averaged, to give the final test result. Generally speaking, a positive result is more credible then a negative one, as it is unlikely that the two shots are consistent blunders leading to a plausible displacement. A negative result is not necessarily sufficient to affirm that the tested shot is a blunder, because either end-point could be responsible for the unreasonable flow-vector. This is why I use many Appendix B. Survey data import and filtering summary 106 Vector component (spherical coordinates) Lower basin Upper basin min max min max V (m/s) 0 45 0 40 6 (deg) -40 30 -60 60 4> (deg) -35 10 -35 20 T A B L E B . l : Decision parameters for ispossiblef low with the level set to "generous". end-points to test a shot: if the tested point is indeed a blunder, each flow calculation will lead to unreasonable result, even if some of the other end-points are themselves blunders. When averaging the i s p o s s i b l e f low results, some weighting can be applied to account for what is known about the data, f a c t o r is a preferential weighting that can be applied to favor the result of tests carried with the end-point either before (f ac tor> 1) or after ( factor< 1) the tested shot. In the case of makecontinuous .m, f a c t o r was set to 1.5 because since shots are tested in sequence, measurements taken before the tested points have already been tested and their chance of being a blunder is low. Weighting is also applied in the following cases. If the interval between the two end-point is less than a year, the impact of negative results is diminished by a factor 0.8 (= weight_short, set in script). If the interval between the two end-point is two years or more, the impact of positive results is diminished by a factor 0.8 (= weight_long,set in script). If the test is positive and the shot is considered not-a-blunder, the result of each ispossibleflow-test carried with a yet-to-be-tested shot is compiled in the conf idence field of the element of the survey-structure corresponding to that shot. That way I have an estimation of the quality of the shots after the currently tested shot. A shot with a "bad reputation" (when most of the flow calculation carried using it yielded unreasonable results) is not used to test subsequent shots. — [ tes t , message] = i sposs iblef low(speed, l o c , l e v e l , dt) Explicitly enough, this function decides if a calculated flow-vector is realistic. The de-cision is based on predefined bounds for the velocity, horizontal direction and vertical direction of the flow, l e v e l ("very generous," "generous", "tight") is a qualitative appreci-ation of the significance level of the test for various setting of these bounds. The bounds are also dependent on whether the pole is in the upper basin or not, and on the delay between the two end-points; shots that are very close in time are allowed more latitude, as the displacement expected is comparable to the accuracy of the instrument. • S = checkqmreset(S) Before 1989, the resetting of poles was not recorded explicitly in the log-files. The reset status of shots taken before that year are therefore unknown, checkqmreset uses test_shot to detect the resets. In two series of tests, each measurement is compared to shots taken before Appendix B. Survey data import and filtering summary 107 and after it. The hypothesis is that a reset will appear to be a blunder when compared to the previous shots, but not when compared to the following ones. So the shot is declared a reset if test_shot returns a negative value when using antecedent measurements and a positive value when using measurements taken after the shot. B.4 Calculate flow: calculateflow.m Now that the data are filtered and sorted by continuous sequences, flow vectors are calculated by differentiation. • Flow = c a l c u l a t e f l o w ( A l l , i n t e r v a l , i n i t i a l j d ) : i n t e r v a l specify over what period of time the flow should be calculated. So far, I only considered annual averages, i n i t i a l j d is the Julian day that serves as reference point for the annual average flow: the function will try to calculate the flow that occurred between that day and the same day the next year, or the closest available point. The output is an array structure, with each element corresponding to a flow vector, described by the fields: — Flow.ID : target ID — Flow.speed : speed vector — Flow. loc : location of vector (midway between the two end-points) — Flow.tsdn : time of vector (time of last end-point) — Flow.Ipt : survey-structure of the initial end-point — Flow.Fpt : survey-structure of the final end-point — F l o w . t s t r : date-string of end-points B.5 Results: Flow maps, 1970-2005 The following pages present the full set of maps of (1) annually-averaged flow maps, and (2) data and outlines. Appendix B. Survey data import and filtering summary 108 Appendix B. Survey data import and filtering summary 109 Appendix B. Survey data import and filtering summary 110 • Surface speed (m/yr) F I G U R E B . 3 : Flow maps: 1989-1994. The length of the flow vectors is five times the animal flow. Appendix B. Survey data import and filtering summary 111 Appendix B. Survey data import and filtering summary 112 F I G U R E B.6: Glacier outlines: 1969-1986. The shading gives the uncertainty. Dots are flow poles and +'s are profile line measurements. Appendix B. Survey data import and filtering summary 114 F I G U R E B . 7 : Glacier outlines: 1990-1998. The shading gives the uncertainty. Dots are flow poles and +'s are profile line measurements. F I G U R E B . 8 : Glacier outlines: 1999—2005. The shading gives the uncertainty. Dots are flow poles and +'s are profile line measurements. Appendix C Bayesian block kriging 116 C l Structural assumptions The field {m(x);x G D} is a realization of a random function {M(x,<f>);x G D,<f> G $} (thereafter referred to as simply M(x)), with its first two moments known a priori: E[M(x)] = u.M(x) xeD ( C l ) Cov[M(x), M{x')} = CM(x,x') x,x' e D. (C.2) The first moments of the random field Z(x) are supposed to obey the following relations E [Z(x)\M(x")\ x" eD]=a0 + M(x') x' G D (C.3) Cov [Z{x), Z(x')\M{x"); x" e D] = Cz{M(x - x') x, x' e D, (C.4) which entails that {Z(X,UJ);X € D,OJ G fl} is second-order stationary around its expected function {ao + m(x); x G D}. C.2 Preliminary definitions The block kriging estimator for the control volume So is N Z(So) = YXi (Z(xi)-m{xi)) + u.M{S0). (C.5) i=i where PM(SQ) = ^ / g u-M(x)dx. It is known from linear Bayesian theory that ElYJ = ElElY^Yi}} (C.6) Coy[YuY2} = E[Coy{YuY2\Y3}} + Coy[E[Y1,Y3},E[Y2,Y3}], (C.7) hence the random field point-to-point covariance, expressed in term of the model covariance (known a priori) and the covariance of the field conditional on the model (estimated from the data) is Cov[Z(x),Z(x')] = E [Cov[Z{x),Zp\M{x"); x" G D\] + Cov [E[Z(x)\M(x"); x" G D],E[Z(x')\M(x"); x" G D]] = E[Cz\M(x - x')) + Cov[a0 + M(x), a0 + M{x')} = CzlM(x - x1) + CM(x, x1). (C.8) Appendix C. Bayesian block kriging 117 The point-to-block and block-to-block covariances are Cz{xUSJ) = ±- f Cov (Z(Xi),Z(x))dx (C.9) 3 J Sj Cz{Si,Sj) = ^— f [ Cov(Z{x),Z(x'))dxdx'. (CIO) and similarly, by mere separation of the integral, the conditional and model point-to-block and block-to-block covariances are obtained Cz\M(xi,Si) = ±- ! Cov (Z(Xi),.Z(x)\M)dx ( C . l l ) B3 JSj Cz\M (Si, Sj) = / / Cov (Z(x),Z(x')\M)dx dx' (C.12) t 3 " Si J Sj CM(xi,Sj) = ±- [ Cov (M(Xi),M(x))dx (C.13) B3 JSj CM(Si, Sj) = c ^ T / / Cov (Z(Xi), Z(x')) dxdx'. (C.14) i 3 JSiJSj From these definition, ( C l ) and (C.3) the following expression is obtained for the expectation of the block average E[ZS] - E[E[ZS\M}} = | J E[M(x)] + a0]dx = ^ J^nM(x)dx + a0 = u.M(S) + a0. (C.15) Since the integrand Cov (Z(x), Z(x')\M) is by assumption stationary, the integrals (C.11)-(C.12) are functions of the euclidian distance and not of position, and therefore CZ\M(xU SJ) and CZ\M{SI, Sj) are stationary. Note that by this result and by (C.15), it is shown that if the point random function Z(x) is second-order stationary around ao + m(x), then so is its block average Zs(x). In the early days of kriging, the analytical solution of the integrals of the point-to-block and block-to-block covariances was critical to the applicability of the method. Nowadays, with the increase in numerical power and the concomitant disregard for analytical solutions, it is customary to approximate these integrals by sums by dividing the volume support into k = pD subsections, Appendix C. Bayesian block kriging 118 where p is the planar division and D the dimensionality of the support: 1 k CZ\M(xi,Sj) -Ycz\M(xi,xJ0) (C.16) 0=1 CzMSuSj) =i r ^ E E ^ I M ^ . ^ ) ( c - 1 7 ) a=l/?=l k CM{xi,S0) =i ^ ^CMixuXjp) (C.18) 0=1 C M { S u S j ) ^ ± Y , Y , C » ' ( x ^ x i t > ) (°-19) fc2 a=l0=1 where the X i a ' s and s arc the center points of the subdivisions. C.3 Unbiasedness Unbiasedness of the estimator requires E[Z(S0) - Z(x0)} = 0 t E[Z(S0)] = E[Z(x0)] (C.20) N E [E[Z(S0)\M]} = YXi E l z ( x i ) ~ t*M(xi)} + HM(SQ) I=I N o-o + PM(SQ) = Y,AI a° + M M ( S O ) 1=1 AT a o = A i a o * = i and hence as it was the case for ordinary kriging, 2~liLi xi = 1 i s a necessary and sufficient condition for unbiasedness. Appendix C. Bayesian block kriging 119 C.4 Estimation error The estimation error variance is o-2(S0) = Var[z(S0)-Z(So) = Var N = Var Z(S0)-Y/\lY(xi) + u.M(S0) 1=1 TV J ^ A i { Z ( 5 0 ) - F ( ^ ) - M M ( 5 0 ) } JV N = J ^ A l A J C o v ( Z ( 5 0 ) - y(xi) - HM(SO), Z(SO) - Y(XJ) - nM{S0)) (C.21) i=i j=i where the relations X ^ i = 1 a n d Var Eil i (^^ )] - Ei l i Ef=i ^ Cov (y(xi), V-(^)) were used. By the unbiasedness property of the estimator the covariance reduces to the second moment Cov (Z(S0) - Y{Xi) - M M ( S O ) , Z(SO) - Y(Xj) - u.M(S0)) = E [(Z(S0) - Y(Xi) - M M ( S O ) ) • (Z(S0) - Y(Xj) - M M ( S O ) ) ] . (C.22) Expanding the product and adding and subtracting two times O Q , E[(Z(So) - Y(Xi) - M M ( S o ) ) • (Z(S0) - Y(Xj) - M M ( S O ) ) ] = + E{{Z{S0)-u.M{SQ))2]-a2 -E[(Z(S0)-m(SoW(Xi)} + a2 -E{(Z(S0)-u.M(S0))Y(Xj)] + a2 + E[Y(xi)Y(xj)]+a20, (C.23) but since E[(Z(SQ) — U,M(SQ))] = E[Y(xi)] = E[Y(xj)] = ao and by the definition of covariance E [(Z(S0) - Y(Xi) - HM(SQ)) • (Z(S0) - Y(Xj) - M M ( S 0 ) ) ] = + V a r [ Z ( 5 0 ) - M M ( 5 0 ) ] -Cov(Z(S0),Y{xi)) -Cov(Z(S0),Y(Xj)) -CaviYixj^Yixj)). (C.24) As covariance is independent of shift of the variables the expectation is E [ ( Z ( S 0 ) - Y{Xi) - HM(S0)) • (Z(S0) - Y(XJ) - M M ( S O ) ) ] = Var [Z(S0)} - Cov (Z(S0), Z(XI)) - Cov (Z(SQ), Z(XJ)) - Cov {Z(XI), Z(XJ)). (C.25) Appendix C. Bayesian block kriging 120 Inserting back (C.25) and (C.22) into (C.21) the block estimation error variance is N N N cr2(5o) - Var[Z(50)] - 2 ^ A i C o v ( Z ( 5 o ) , ^ i ) ) + E ^ A i A J - C o v ( Z ( a : i ) , Z ( x j ) ) , (C.26) i=l i=l j=l which by (C.8) and definitions (C.11)-(C.14) can be rewritten N <T2(S0) = [Cz[M(0) + CM(S0,S0)} - 2 ^ A i [ C z | M ( S 0 , ^ ) + C M ( S o , Z i ) ] *=i N N + YYtxiX3 [CZ\M{xi-Xj) + CM{xi:Xj)] . (C.27) i=l 3 = 1 C.5 Normal equations The objective function for the minimization of the estimation error with unbiasedness constraint is the lagrangian L(Ai , ...,\N,3)= a2(S0) - 28 ^ K - lj • (C.28) Derivation with respect to the Aj's and the lagrange multiplier B yields the following system of equations Y Xi [Cz\M(xi - x3) + CM{xi, Xj)] + 3= [Cz\M(So, x5) + CM{S0, Xj)] , for j = 1... N (C.29) i=l N $3^=1. (C30) i=i Comparison with the ordinary kriging set of normal equations reveals two differences. First, the right-hand-side of (C.29) is dependent on an estimation-block to data-point covariance, instead of the usual point-to-point covariance. This is characteristic of block kriging, and occurs similarly in ordinary block kriging. Second, the covariance terms reflect the presence of two sources of variance: that of the random field conditional on the prior and that of the prior itself. Note that the prior field M(x) is not assumed to be second-order stationary or intrinsic: its mean is not constant, and its covariance depends on the specific locations of the points, not only on their relative distance. Hence, the random field Z(x) itself is not second-order stationary. C .6 Optimal estimation error variance Multiplying the j th equation (C.29) by Xj and summing the N equations one obtains N N N ^ ^ A i A j [CZ\M(XZ - XJ) + CM{xi,Xj)] +8 = Yxi [CZ\M(SO,XJ) + CM(S0,Xj)] . (C.31) i=l j=l j = l Appendix C. Bayesian block kriging 121 Substituting Y J i = 1 J2j=i [CZ\M(X% - Xj) + CM{xi,Xj)] in (C.27) yields N N a*2(S0)= [Cz\m(0) + CM(SQ,SO)] - 2 ^ ^ A J A J [Cz\M(xi - x3) + CM{xi,Xj)] i=l j=\ N N + ^ ^ A i A j [Cz\M(xi - Xj) + CM(xi,Xj)] -3 i=l j=l hence the optimal estimation error variance is N N <T*2(S0) = [CZ]M(0) + CM(S0,SO)] A I A J [cwixi - x5) + CM{xuXj)] - 3. (C.32) i=l j = l C.7 Algorithm summary and numerical implementation Table C.2 summarizes the kriging procedure in matrix form. Myers (1984) shows how a matrix ordinary kriging system can easily be extended to allow multivariate attributes (cokriging), poly-nomial mean-functions (universal kriging) and to solve for more than one location at the time. The universal cokriging code by Marcotte (1991), which I used as a base for the programming of the block Bayesian kriging, implements these result. At the risk of burdening the notation, I present in Table C.2 the multiple-locations matrices instead of the traditional single-location ones, hoping that it will facilitate the deciphering of the code (should it ever be attempted). Note that the system presented is essentially identical to that of ordinary kriging; only a few elements in the definition of the matrices are different. Among them, the non-stationarity of the CM matrices is most problematic. Whereas the rest of the code depends only on euclidian distances, the calculation of the data points fluctuation Y(xi) and that of the prior covariance matrices neces-sitates the evaluation of the prior mean and variance functions ((J-M> O~M) a f definite locations. This rupture of symmetry complicates the code. In addition, since no analytical form for these functions are derived, the evaluation must be carried by linear interpolation of the gridded fields. Appendix C. Bayesian block kriging 122 Estimator JV Z(S0) = ]PAi (Z{xi) - nM{xi)) + U.M{SQ). i=l Normal equations YXi [CZ\M{XI - Xj) + CM(xi,Xj)] + 8 = [CZ\M(S0,Xj) + CM(SQ,XJ)] , for j = 1 ...N i=l JV i=l Optimal estimation error variance JV c r * 2 ( S 0 ) = [ C Z | M ( 0 ) + C M ( 5 O , 5 O ) ] [ C ^ f o , 5 0 ) + CM{xu S0) ] - 8. i=l MATRIX EQUATIONS Estimator Zs = TTz + n Normal equations KT = K0 Optimal estimation error variance cr2 = C 5 Z | M + C S M - d i a g ( r T K 0 ) Matrix expression Notes and associ-ated v a r i a b l e s ( Zsixo1) \ block estimation. ZS Zs = (raxl) \ Zs(xom) ) Z (mxl ) / z{x{) - HM(XI) \ z = z(xn) - HM{xn) \ o / fluctuation of data around the prior. Appendix C. Bayesian block kriging 123 (mx l ) r ((n+l)xm) B ( m x l ) A (nxm) K ((n+l)x(n+l)) Cz\M (n X n) C M (nxn) ( ( n + l ) x m ) CZ\M0 (nxm) CM0 (nxm) a ( m x l ) M M ^ O 1 ) \ flM(S0m) J r = A BT B A = i f = / P1 \ \0m J A i 1 . . . A i m \ A n 1 • • • A n m y C Z | M + CM 1 1 0 /Cz]M(xi-xi) . . . C W O E I - z„) \ CZ\M — \Cz\M{xi ( CM(xi,xi) CM = C Z \ M ( x n - Xn)J CM(xi,Xn)\ \CM(xi,xn) K0 = CM(Xu, xn) J CZ\M0 + C M 0 1 c Z\M0 CM0 — Cz^ixuSo1) ... czlM(xi,s0my CZ\M(xn,So1) ... Cz\M(xn, Som) j I CMixuSo1) ... CM(xuS0m) \ \ CM{xn,SQl) ... CM{xn,S0m) j a2 = \ < ^ 2 ( § o m ) prior block average at es-timation locations. lagrange parameter for the m estimation locations. kriging weights. conditional covariance matrix. Dependent only on the. distance. prior covariance matrix. Dependent on the location. conditional data-point to estimation-block covariance matrix. prior data-point to estimation-block covari-ance matrix. estimation error variance. Appendix C. Bayesian block kriging 124 ( m x m ) CSM ( m x m ) CSM = fCz\M(Sol,Sol) ... CZ\M{Bo1,Som)^ \Cz\M{Son,So1) ... Cz\M(Son,Som)J (CM{SQl,Sj) ... CM(S0\S0m)\ \CM(Son,S01) ... CM(S0n,S0m)J conditional estimation-block to estimation-block covariance matrix. prior estimation-block to estimation-block covariance matrix. where HM{SQ) = --j- / u.M{x)dx So JSo Cz\M(xi, Si) = / Cov {Z(x)h Z(x)\M) dx CZ\M(Si,Si) = -±- [ f Cov (Z(x),Z(x')\M)dxdx' JSiJSj CM(xi, S^ = I Cov (M(xi), M(x)) dx CM(Si,Si) = ^ - f f Cov(Z(x),Z(x'))dxdx' T A B L E C . 2 : The Bayesian block kriging system, n is the number of data points and m the number of estimation locations. 125 Appendix D S t a t i s t i c a l a n a l y s i s o f D E M t i m e s e r i e s Mean Variance Normality Year N Qi Test Q2 Q2min Test L L m a x Test 1969 25 -0.08 0.41 0 1.13 0.52 1.64 0 0.16 0.18 0 1970 23 -0.02 0.43 0 1.09 0.50 1.67 0 0.22 0.18 1 1971 22 0.08 0.44 0 0.08 0.49 1.69 1 0.17 0.19 0 1972 23 -0.06 0.43 0 0.29 0.50 1.67 1 0.11 0.18 0 1974 20 -0.06 0.46 0 0.07 0.47 1.73 1 0.08 0.20 0 1980 57 -0.00 0.27 0 0.25 0.63 1.37 1 0.10 0.12 0 1981 66 0.05 0.25 0 0.28 0.65 1.35 1 0.13 0.11 1 1982 51 0.01 0.28 0 0.96 0.60 1.40 0 0.28 0.13 1 1983 100 -0.01 0.20 0 1.70 0.72 1.28 1 0.13 0.09 1 1984 66 0.00 0.25 0 0.84 0.65 1.35 0 0.21 0.11 1 1985 113 -0.08 0.19 0 1.42 0.74 1.26 1 0.15 0.08 1 1986 137 0.10 0.17 0 o:87 0.76 1.24 0 0.18 0.08 1 1987 144 -0.03 0.17 0 0.55 0.77 1.23 1 0.13 0.07 1 1988 156 0.07 0.16 0 0.54 0.78 1.22 1 0.15 0.07 1 1989 178 -0.04 0.15 0 0.32 0.79 1.21 1 0.13 0.07 1 1990 277 0.05 0.12 0 0.25 0.83 1.17 1 0.12 0.05 1 1991 174 -0.05 0.15 0 0.35 0.79 1.21 1 0.14 0.07 1 1992 64 -0.06 0.25 0 0.14 0.65 1.35 1 0.15 0.11 1 1993 197 -0.01 0.14 0 0.21 0.80 1.20 1 0.09 0.06 1 1994 280 -0.01 0.12 0 0.22 0.83 1.17 1 0.09 0.05 1 1995 507 -0.03 0.09 0 0.28 0.88 1.12 1 0.13 0.04 1 1996 393 -0.01 0.10 0 0.20 0.86 1.14 1 0.09 0.04 1 1997 399 0.03 0.10 0 0.22 0.86 1.14 1 0.11 0.04 1 1998 272 0.08 0.12 0 0.52 0.83 1.17 1 0.20 0.05 1 1999 262 0.00 0.12 0 0.13 0.83 1.17 1 0.10 0.05 1 2000 204 -0.04 0.14 0 0.20 0.80 1.20 1 0.11 0.06 1 2001 139 0.01 0.17 0 0.08 0.76 1.24 1 0.11 0.08 1 2002 187 -0.04 0.15 0 0.14 0.79 1.21 1 0.15 0.06 1 2003 209 -0.06 0.14 0 1.13 0.81 1.19 0 0.31 0.06 1 2004 221. -0.00 0.13 0 1.15 0.81 1.19 0 0.29 0.06 1 2005 202 -0.02 0.14 0 0.61 0.80 1.20 1 0.20 0.06 1 T A B L E D . l : Orthonormal residuals statistics for all years, Bayesian model. The correlation range is 500m, and the variance is 100m2. This is the model that was selected. Appendix D. Statistical analysis of DEM time series 126 Mean Variance Normality Year N Qi Test Q 2 Test L Lrnax Test 1969 25 0.89 0.41 1 13.38 0.52 1.64 1 0.18 0.18 1 1970 23 0.85 0.43 1 13.52 0.50 1.67 1 0.21 0.18 1 1971 22 -0.07 0.44 0 14.25 0.49 1.69 1 0.20 0.19 1 1972 23 0.11 0.43 0 13.13 0.50 1.67 1 0.16 0.18 0 1974 20 0.19 0.46 0 13.66 0.47 1.73 1 0.20 0.20 1 1980 57 0.33 0.27 1 4.43 0.63 1.37 1 0.11 0.12 0 1981 66 0.18 0.25 0 5.60 0.65 1.35 1 0.11 0.11 0 1982 51 0.45 0.28 1 6.21 0.60 1.40 1 0.06 0.13 0 1983 100 0.46 0.20 1 6.33 0.72 1.28 1 0.16 0.09 1 1984 66 0.02 0.25 0 2.83 0.65 1.35 1 0.17 0.11 1 1985 113 -0.06 0.19 0 3.58 0.74 1.26 1 0.16 0.08 1 1986 137 -0.20 0.17 1 2.39 0.76 1.24 1 0.12 0.08 1 1987 144 -0.07 0.17 0 1.88 0.77 1.23 1 0.15 0.07 1 1988 156 0.11 0.16 0 2.34 0.78 1.22 1 0.15 0.07 1 1989 178 0.00 0.15 0 1.14 0.79 1.21 0 0.12 0.07 1 1990 277 -0.01 0.12 0 1.03 0.83 1.17 0 0.14 0.05 1 1991 174 -0.12 0.15 0 1.34 0.79 1.21 1 0.16 0.07 1 1992 64 0.02 0.25 0 0.74 0.65 1.35 0 0.13 0.11 1 1993 197 0.02 0.14 0 0.98 0.80 1.20 0 0.11 0.06 1 1994 280 -0.02 0.12 0 0.91 0.83 1.17 0 0.11 0.05 1 1995 507 -0.01 0.09 0 1.46 0.88 1.12 1 0.11 0.04 1 1996 393 0.01 0.10 0 0.82 0.86 1.14 1 0.11 0.04 1 1997 399 0.06 0.10 0 1.58 0.86 1.14 1 0.13 0.04 1 1998 272 -0.07 0.12 0 1.88 0.83 1.17 1 0.14 0.05 1 1999 262 0.06 0.12 0 0.55 0.83 1.17 1 0.06 0.05 1 2000 204 -0.08 0.14 0 1.06 0.80 1.20 0 0.12 0.06 1 2001 139 0.15 0.17 0 0.85 0.76 1.24 0 0.13 0.08 1 2002 187 0.07 0.15 0 1.12 0.79 1.21 0 0.14 0.06 1 2003 209 0.20 0.14 1 3.07 0.81 1,19 1 0.23 0.06 1 2004 221 0.03 0.13 0 3.38 0.81 1.19 1 0.23 0.06 1 2005 202 -0.09 0.14 0 1.04 0.80 1.20 0 0.16 0.06 1 T A B L E D . 2 : Orthonormal residuals statistics for all years, time-independent second-order stationary model interpolated by ordinary kriging. The correlation range is 500 m, and the variance is 100 m 2 . PM Cz\M Mean Variance Normality Year N Range Nugget Range Variance Qi Test Q2 Q2MIN Q2 Test L Lmax Test 1969 25 260.00 0.00 213.82 127.45 0.08 0.41 0 0.98 0.52 1.64 0 0.15 0.18 0 1970 23 260.00 0.00 229.76 139.79 -0.10 0.43 0 0.98 0.50 1.67 0 0.22 0.18 1 1971 22 260.00 0.00 356.13 8.02 0.14 0.44 0 0.21 0.49 1.69 1 0.15 0.19 0 1972 23 260.00 0.00 515.01 53.51 0.28 0.43 0 0.86 0.50 1.67 0 0.14 0.18 0 1974 20 260.00 0.00 0.01 5.88 -0.03 0.46 0 0.42 0.47 1.73 1 0.10 0.20 0 •a C6 1980 57 260.00 0.00 275.07 50.17 0.34 0.27 1 0.78 0.63 1.37 0 0.10 0.12 0 a a. 1981 66 260.00 0.00 0.01 29.20 . -0.08 0.25 0 0.67 0.65 1.35 0 0.19 0.11 1 1982 51 260.00 0.00 60.71 28.35 -0.02 0.28 0 0.61 0.60 1.40 0 0.22 0.13 1 1983 100 260.00 0.00 122.51 47.64 0.10 0.20 0 0.79 0.72 1.28 0 0.13 0.09 1 . Statistic 1984 66 260.00 0.00 142.32 29.99 0.03 0.25 0 0.69 0.65 1.35 0 0.20 0.11 1 . Statistic 1985 113 260.00 0.00 136.51 59.92 0.08 0.19 0 0.59 0.74 1.26 1 0.20 0.08 1 . Statistic 1986 137 260.00 0.00 154.41 35.25 0.04 0.17 0 0.66 0.76 1.24 1 0.18 0.08 1 . Statistic 1987 144 260.00 0.00 240.20 41.60 -0.02 0.17 0 0.55 0.77 1.23 1 0.11 0.07 1 1988 156 260.00 0.00 986.00 142.58 -0.03 0.16 0 0.56 0.78 1.22 1 0.15 0.07 1 i" 1989 178 260.00 0.00 1272.00 79.86 -0.00 0.15 0 0.64 0.79 1.21 1 0.12 0.07 1 .alysis 0 1990 277 260.00 0.00 1406.00 88.89 -0.10 0.12 0 0.65 0.83 1.17 1 0.09 0.05 1 .alysis 0 1991 174 260.00 0.00 215.07 19.17 -0.05 0.15 0 0.47 0.79 1.21 1 0.11 0.07 1 .alysis 0 1992 64 260.00 0.00 221.01 8.32 0.10 0.25 0 0.28 0.65 1.35 1 0.15 0.11 1 to 1993 197 260.00 0.00 1758.00 122.46 0.00 0.14 0 0.50 0.80 1.20 1 0.11 0.06 1 1994 280 260.00 0.00 124.88 6.95 -0.03 0.12 0 0.56 0.83 1.17 1 0.08 0.05 1 1995 507 260.00 0.00 317.63 32.10 0.03 0.09 0 0.33 0.88 1.12 1 0.15 0.04 1 S' 1996 393 260.00 0.00 138.13 5.33 0.05 0.10 0 0.62 0.86 1.14 1 0.09 0.04 1 S 1997 399 260.00 0.00 260.82 20.22 0.04 0.10 0 0.43 0.86 1.14 1 0.08 0.04 1 0) 1998 272 260.00 0.00 16.81 5.99 0.04 0.12 0 0.68 0.83 1.17 1 0.18 0.05 1 SA' 1999 262 260.00 0.00 106.51 3.22 0.06 0.12 0 0.49 0.83 1.17 1 0.09 0.05 1 2000 204 260.00 0.00 103.18 5.93 -0.05 0.14 0 0.60 0.80 1.20 1 0.10 0.06 1 2001 139 260.00 0.00 76.20 1.72 -0.04 0.17 0 0.56 0.76 1.24 1 0.11 0.08 1 2002 187 260.00 0.00 75.84 4.03 0.02 0.15 0 0.37 0.79 1.21 1 0.11 0.06 1 2003 209 260.00 0.00 36.10 25.64 0.01 0.14 0 0.73 0.81 1.19 1 0.32 0.06 1 2004 221 260.00 0.00 39.38 23.61 -0.02 0.13 0 0.42 0.81 1.19 1 0.25 0.06 1 2005 202 260.00 0.00 117.57 14.60 0.01 0.14 0 0.55 0.80 1.20 1 0.20 0.06 1 T A B L E D . 3 : Orthonormal residuals statistics and covariance parameters; time-dependent Bayesian model. to Year JV PM Range Cz\M Mean Variance Normality Nugget Range Variance Qi ' Test Q2 Q 2 m , „ Test L Lmax Test 1969 25 213.82 0.00 213.82 127.45 0.12. 0.41 0 15.39 0.52 1.64 1 0.16 0.18 0 1970 23 229.76 0.00 229.76 139.79 -0.09 0.43 0 14.12 0.50 1.67 1 0.17 0.18 0 1971 22 1295.00 0.00 1295.00 107.23 -0.45 0.44 1 16.81 0.49 1.69 1 0.16 0.19 0 1972 23 515.01 0.00 515.01 53.51 -2.03 0.43 1 28.84 0.50 1.67 1 0.18 0.18 0 1974 20 131.01 0.00 131.01 24.35 1.55 0.46 1 84.77 O.47 1.73 1 0.16 0.20 0 1980 57 275.07 0.00 275.07 50.17 0.11 0.27 0 9.23 0.63 1.37 1 0.12 0.12 1 1981 66 0.01 0.00 0.01 29.20 4.81 0.25 1 385.20 0.65 1.35 1 0.10 0.11 0 1982 51 2072.00 0.00 2072.00 1325.64 -0.11 0.28 0 1.12 0.60 1.40 0 0.07 0.13 0 1983 100 131.26 0.00 131.26 73.74 0.01 0.20 0 7.99 0.72 1.28 1 0.13 0.09 1 1984 66 1564.00 0.00 1564.00 577.06 -0.04 0.25 0 0.90 0.65 1.35 0 0.14 0.11 1 1985 113 113.76 0.00 113.76 39.30 -0.09 0.19 0 9.36 0.74 1.26 1 0.12 0.08 1 1986 137 549.01 0.00 549.01 181.82 0.05 0.17 0 1,45 0.76 1.24 1 0.13 0.08 1 1987 144 343.13 0.00 343.13 44.27 -0.06 0.17 0 3.49 0.77 1.23 1 0.14 0.07 1 1988 156 1104.00 0.00 1104.00 195.20 -0.04 0.16 0 2.30 0.78 1.22 1 0.14 0.07 1 . 1989 178 1262.00 0.00 1262.00 99.64 0.09 0.15 0 2.34 0.79 1.21 1 0.12 0.07 1 1990 277 1928.00 0.00 1928.00 115.43 0.07 0.12 0 2.89 0.83 1.17 1 0.11 0.05 1 1991 174 242.82 0.00 242.82 19.28 -0.02 0.15 0 5.59 0.79 1.21 1 0.18 0.07 1 1992 64 337.26 0.00 337.26 23.29 -0.08 0.25 0 2.20 0.65 1.35 1 0.16 0.11 1 1993 197 1105.00 0.00 1105.00 131.22 0.02 0.14 0 1.47 0.80 1.20 1 0.11 0.06 1 1994 280 263.63 0.00 263.63 18.98 0.05 0.12 0 3.57 0.83 1.17 1 0.15 0.05 1 1995 507 686.01 0.00 686.01 148.51 -0.05 0.09 0 1.36 0.88 1.12 1 0.10 0.04 1 1996 393 804.01 0.00 804.01 50.45 -0.09 0.10 0 2.62 0.86 1.14 1 0.10 0.04 1 1997 399 802.51 0.00 802.51 158.97 0.05 0.10 0 1.67 0.86 1.14 1 0.13 0.04 1 1998 272 • 18.49 0.00 18.49 8.91 0.10 0.12 0 46.17 0.83 1.17 1 0.14 0.05 1 • 1999 262 1704.00 0.00 1704.00 148.02 -0.01 0.12 0 1.21 0.83 1.17 1 0.06 0.05 1 2000 204 252.95 0.00 252.95 20.94 0.09' 0.14 0 3.61 0.80 1.20 1 0.14 0.06 1 2001 139 1696.00 0.00 1696.00 58.81 0.20 0.17 1 3.69 0.76 1.24 1 0.11 0.08 1 2002 187 273.76 0.00 273.76 25.63 0.11 0.15 0 . 3.42 0.79 1.21 1 0.16 0.06 1 2003 209 26.26 0.00 26.26 19.61 0.18 0.14 1 11.34 0.81 1.19 1 0.15 0.06 1 2004 221 735.01 0.00 735.01 298.82 0.05 0.13 0 1.54 0.81 1.19 1 0.22 0.06 1 2005 202 52.26 0.00 52.26 5.16 0.12 0.14 0 20.78 0.80 1.20 1 0.13 0.06 1 T A B L E D . 4 : Orthonormal residuals statistics and covariance parameters; second-order stationary model interpolated by ordinary kriging. to 0 0 Appendix E Results: Color plates 129 i i 1 1 1 1 1 = — 1 1.5 2 2.5 3 3.5 4 Eas t i ng ( k m ) F I G U R E E . l : Profiles: 1951-1980. The length of the flow vectors is five times the annual flow. Appendix E. Results: Color plates 130 Easting (km) F I G U R E E . 2 : Profiles: 1981-1987. The length of the flow vectors is five times the annual flow. Appendix E. Results: Color plates 131 11988 1989* 11990' 1991I 11992" 19931 100 0 -H994* 1.5 2.5 3 Easting (km) 3.5 F I G U R E E . 3 : Profiles: 1988-1994. The length of the flow vectors is five times the annual flow. Appendix E. Results: Color plates 132 - 1 r 2.5 3 Easting (km) 1.5 3.5 F I G U R E E . 4 : Profiles: 1995-2001. The length of the flow vectors is five times the annual flow. Appendix E. Results: Color plates 133 Easting (km) F I G U R E E . 5 : Profiles: 2002-2005. The length of the flow vectors is five times the annual flow. Appendix E. Results: Color plates 134 F I G U R E E . 6 : Ice thickness and elevation difference from previous model: 1951-1974. Appendix E. Results: Color plates 135 0 20 40 60 80 100 120 -12 -8 -4 0 4 8 12 Ice thickness (m) Elevation change (m) F I G U R E E.7: Ice thickness and elevation difference from previous model: 1980-1985. Appendix E. Results: Color plates 136 ^ H i : : t : ii [ | i :- ( J^BBBBBBBBBBBBBBBBBBBBBH 0 20 40 60 80 100 120 -12 -8 -4 0 4 8 12 Ice thickness (m) Elevation change (m) F I G U R E E . 8 : Ice thickness and elevation difference from previous model: 1986-1991. Appendix E. Results: Color plates 137 Ice thickness (m) Elevation change (m) F I G U R E E . 9 : Ice thickness and elevation difference from previous model: 1992-1997. Appendix E. Results: Color plates 138 0 20 40 60 80 100 120 -12 -8 -4 0 4 8 12 Ice thickness (m) Elevation change (m) F I G U R E E.10: Ice thickness and elevation difference from previous model: 1998-2003. Appendix E. Results: Color plates 139 40 60 80 100 120 -12 -8 -4 0 4 8 12 Ice thickness (m) Elevation change (m) F I G U R E E . l l : Ice thickness and elevation difference from previous model: 2004-2005. Appendix E. Results: Color plates 140 0 6 12 18 24 5 00 m ftf Standard deviation (m) F I G U R E E . 1 2 : Kriging standard deviation: 1951-1985. Appendix E. Results: Color plates 141 Standard deviation (m) F I G U R E E . 1 3 : Kriging standard deviation: 1986-1997. Appendix E. Results: Color plates 142 0 6 12 18 24 5 00 m ftf Standard deviation (m) F I G U R E E . 1 4 : Kriging standard deviation: 1998-2005. Appendix E. Results: Color plates 143 Driving stress (kPa) F I G U R E E . 1 5 : Driving stress: 1951-1985. Appendix E. Results: Color plates 144 ^ , , , , -40 0 40 80 120 160 200 500 m J\J Driving stress (kPa) F I G U R E E . 1 6 : Driving stress: 1986-1997. Appendix E. Results: Color plates 145 FIGURE E .17: Driving stress: 1998-2005. Appendix E. Results: Color plates 146 I , L I 1 0 6 12 18 24 30 36 500 m Creep (m/yr) F I G U R E E .18 : Maximal creep velocity: 1951-1985. Appendix E. Results: Color plates 147 Appendix E. Results: Color plates 148 Creep (m/yr) F I G U R E E . 2 0 : Maximal creep velocity: 1998-2005.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Slow surge of Trapridge Glacier, Yukon Territory, 1951-2005
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Slow surge of Trapridge Glacier, Yukon Territory, 1951-2005 Frappe-Seneclauze, Tom-Pierre 2006
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Slow surge of Trapridge Glacier, Yukon Territory, 1951-2005 |
Creator |
Frappe-Seneclauze, Tom-Pierre |
Publisher | University of British Columbia |
Date Issued | 2006 |
Description | Trapridge Glacier, a surging glacier located in the St. Elias Mountains, Yukon Territory, Canada, went through a complete surge cycle between 1951 and 2005. Air photos (1951-1981) and groundbased optical surveys (1969-2005) are used to describe the modifications in flow and geometry that occurred over this period. The acceleration of flow during the surge is detected by repeated measurement of poles drilled into the ice. Between 1974 and 1980, the median velocity in the lower basin went from 15.6ma⁻¹ to 38.6ma⁻¹. Downstream from this zone of active flow, coldbased ice accumulated during the previous surge impeded the flow, and a steep front formed at the boundary between the two ice masses. Over the following ten years, this bulge propagated downglacier, advancing faster than the ice and integrating stagnant ice by continuous deformation. After it peaked in 1984, the flow in the lower basin remained above 25ma⁻¹ for 12 years, but was on a slowing trend. The slowdown followed strangely regular 4-year pulses: 1-2 years of timid acceleration (2-3ma⁻¹) , followed by 2-3 years of rapid deceleration (4-8ma⁻¹) . The 1997-1999 acceleration was particularly vigorous, as the median velocity went from 20.3ma⁻¹ to 28.5ma⁻¹. After this last pulse, the glacier gradually slowed down to pre-surge velocities. In 2005, the lower basin was flowing at less than 8.5ma⁻¹. Based on this flow history, I divide the 1969-2005 period into three parts: the activation phase (1969-1974), the surge (1980-1999), and the waning phase (1999-2005). To quantify the geometrical changes that occurred during each phase, digital elevation models are constructed from air photos and optical survey measurements. Optical and radar surveys are joined to photogrammetric measurements of the proglacial field to obtain a bed topography map. DEMs for 1951, 1970, 1972, 1977, and 1981 are generated by stereographic analysis of air photos. These background models are then updated year after year by ground-based survey data, using a Bayesian kriging algorithm. By assuming the ice surface to be similar in shape from one year to the next, only shifted vertically, this algorithm allows the propagation and diffusion of information on the ice surface topography. The interpolation obeys the data, and follows the prior model where no data are available. Uncertainty is attributed to each estimation, based on the estimated covariance structure of the field and on the uncertainty of the prior model. The assumption of stationarity and the parametrization of the covariance are tested by the analysis of orthonormal residuals. Despite the significant changes in surface area that occurred between 1969 and 2005, the total volume of ice remained relatively constant at 0.14 ± 0.03 km³. The spread of the glacier during the surge was accompanied by a corresponding decrease in thickness. From the surface slope and ice thickness fields, the stress driving the flow is calculated. After filtering the surface slope to eliminate small-scale variations, the shallow ice approximation is used to estimate creep velocities from the stress and ice thickness fields. Deformation of the ice is shown to accounts for a negligible portion of the flow observed in the lower basin. Sliding or bed deformation must therefore be responsible for most of the motion in that area. In the upper basin, thicker and flowing more slowly, creep accounts for maybe as much as half of the recorded motion. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-05 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052585 |
URI | http://hdl.handle.net/2429/17532 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2006-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2006-0042.pdf [ 33.54MB ]
- Metadata
- JSON: 831-1.0052585.json
- JSON-LD: 831-1.0052585-ld.json
- RDF/XML (Pretty): 831-1.0052585-rdf.xml
- RDF/JSON: 831-1.0052585-rdf.json
- Turtle: 831-1.0052585-turtle.txt
- N-Triples: 831-1.0052585-rdf-ntriples.txt
- Original Record: 831-1.0052585-source.json
- Full Text
- 831-1.0052585-fulltext.txt
- Citation
- 831-1.0052585.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052585/manifest