A MODEL FOR HOLOCENE GLACIAL EROSION AT PEYTO GLACIER, ALBERTA by Robert James Vogt B.Sc., University of Northern British Columbia, 2010 THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE IN NATURAL RESOURCES AND ENVIRONMENTAL STUDIES UNIVERSITY OF NORTHERN BRITISH COLUMBIA DECEMBER 2014 © Robert Vogt, 2014 UMI Number: 1526489 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion. Di!ss0?t&iori P iiblist’Mlg UMI 1526489 Published by ProQuest LLC 2015. Copyright in the Dissertation held by the Author. Microform Edition © ProQuest LLC. All rights reserved. This work is protected against unauthorized copying under Title 17, United States Code. ProQuest LLC 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, Ml 48106-1346 Abstract Glaciers play an important role in the evolution of many mountain landscapes. The primary objective of this study was to model rates of contemporary glacier erosion through numerical modeling. The research uses a regional glaciation model and couples to it (i) a process-driven erosion model that includes abrasion and quarrying; and (ii) a sliding-based power relation erosion model. The coupled erosion models were then used to estimate erosion for Peyto Glacier. The models predict higher rates of erosion than those estimated from Peyto Lake sediments. Contrary to observed sediment yield, both models predict a decrease in erosion during glacial retreat. Changes in sediment storage are believed to account for the discrepancy seen. Over the past 2000 years, major differences exist between the models during times of abrupt climate change. These differences depend on the inter-decadal to century-scale variability inherent in the climate proxy data used to force the models. Acknowledgements This thesis could not have been completed without the support of many people. Fore­ most, I would like to thank my supervisor, Dr. Brian Menounos, for his patience and support through my transition into the earth sciences, my period of over-teaching, and the thesis journey itself. I also thank my committee members Dr. Peter Jackson and Dr. Matt Reid for their insightful contributions and feedback. Thank you to Dr. Gwenn Flowers for taking the time as an external examiner. Many thanks go out to Garry Clarke and Christian Schoof for their help in kickstarting my understanding of glacier physics. Thank you to Theo Mlynowski for helpful discus­ sions and results regarding Peyto Lake sediment yield observations. I would also like to thank NSERC for helping support this thesis through a CGS scholarship. My deepest thanks go to my family for their many years of love and support throughout my university journey. I couldn’t have made it to this point without the enduring support of my parents. Thanks to my wife Natalie for support through periods of long days at the office. Finally, thank you to my office-mates Joe Shea and Christina Tennant for engaging in useful discussion and keeping the mood light. Contents Abstract i Acknowledgements ii Table of Contents iii List of Tables vi List of Figures vii Table of Symbols x 1 Introduction 1 2 Development of a Distributed Glacier Erosion Model 3 2.1 A b stra c t............................................................................................................ 3 2.2 Introduction...................................................................................................... 3 2.3 Study a r e a ......................................................................................................... 5 2.3.1 Mass b a la n c e ....................................................................................... 8 2.4 Ice m o d e l ......................................................................................................... 11 2.4.1 2.5 Ice m echanics....................................................................................... 11 2.4.2 Existing ice m o d e l ............................................................................. 15 2.4.3 Glacial hydrology................................................................................ 16 2.4.4 Model optimization............................................................................. 18 Erosion m o d e ls ............................................................................................... 21 2.5.1 21 Process-driven erosion model iii .......................................................... 2.5.2 Sliding erosion model ................................................................. 28 2.6 R esults............................................................................................................... 28 2.6.1 Model O p tim izatio n .......................................................................... 28 2.6.2 Process model ................................................................................... 31 2.6.3 Sliding m o d el...................................................................................... 39 2.7 Discussion andConclusions 3 4 ........................................................................... 41 Evaluating the performance of a glacier erosion model applied to Peyto Glacier, Alberta, Canada 45 3.1 A b stra c t............................................................................................................ 45 3.2 Introduction...................................................................................................... 45 3.3 Methods and d a t a ............................................................................................ 47 3.3.1 Sediment s i n k s ................................................................................... 48 3.3.2 Sediment b u d g e t ................................................................................ 51 3.3.3 Uncertainty a n a ly s is .......................................................................... 54 3.4 R esults............................................................................................................... 55 3.4.1 Model comparison to lake records, 1 9 1 7 -2 0 1 0 .............................. 55 3.4.2 Modeled sediment yieldfor the past 2000 y e a r s ............................... 57 3.4.3 Sediment budget im p licatio n s.......................................................... 63 3.5 Discussion......................................................................................................... 64 3.6 C onclusions...................................................................................................... 70 Conclusion 72 References 75 A Study Area A-2 B Sensitivity Analysis A-5 v List of Tables 2.1 Ice model param eters.................................................................................... 20 2.2 Optimization of ice model p a ra m e te rs ......................................................... 31 2.3 Parameter values used for abrasion sensitivity analysis ............................. 33 2.4 Grain-size distribution parameters used for abrasion sensitivity analysis . 33 2.5 Parameter values used for quarrying sensitivity analysis............................. 35 3.1 Average sediment deposit thicknesses required by process m o d e l 63 A. 1 Peyto Glacier ice thickness m easurem ents..............................................A-3 vi List of Figures 2.1 Map of Peyto Lake drainage b a s i n ................................................................ 6 2.2 Map of model dom ain...................................................................................... 7 2.3 Observed net mass balance by elevation b a n d ............................................. 9 2.4 Temperature reconstructions considered for the s tu d y ................................. 10 2.5 Vector diagram of ice m ovem ent................................................................... 12 2.6 Vector diagram of the abrasion process.......................................................... 22 2.7 Vector diagram of the quarrying p ro c e s s ....................................................... 26 2.8 Peyto Glacier velocity m a p ............................................................................. 29 2.9 Sensitivity analysis of the mean abrasion r a t e ............................................ 34 2.10 Sensitivity analysis of the mean quarrying r a te ........................................... 36 2.11 Process and sliding model comparison with icea r e a ................................... 37 2.12 Erosion rate versus sliding velocity for each p i x e l ..................................... 38 2.13 Erosion rate versus ice thickness for each pixel .......................................... 39 2.14 Erosion rate map of Peyto G l a c ie r ............................................................... 40 3.1 Aerial view of the sediment transport p a t h .................................................... 48 3.2 Aerial view of isolated forefield deposits....................................................... 50 3.3 Sum of observed sediment yield from lake and delta with observed ice extent 51 3.4 Digitized regions of forefield d ep o sitio n ....................................................... 53 3.5 Mean modeled and observed sediment yield (1917-2010).......................... 56 3.6 Abrasion and quarrying components of the sediment yield (1917-2010) . 57 3.7 Process and sliding sediment yield 0-2000 .................................................... 59 3.8 Sediment yield compared to ice volume over 0-2000 ................................. 60 3.9 Relationship between quarrying rate and sliding speed 61 vii .............................. 3.10 Adjusted sliding model comparison to process model sediment yield. . . 3.11 Flow diagram of the influence of melt on sliding v e lo c ity ................ 62 67 A. 1 Location of Peyto Glacier mass balance stakes....................................A-2 A. 2 Observed winter mass balance by elevation band A.3 Observed summer mass balance by elevation b a n d ................................ A-4 .......................................A-4 B. 1 Sensitivity of the abrasion model to the parameter m .................................. A-5 B.2 Sensitivity of the abrasion model to the parameter kabr............................... A-6 B.3 Sensitivity of the abrasion model to the parameter G .................................. A-6 B.4 Sensitivity of the abrasion model to the parameter Zr ............................... A-7 B.5 Sensitivity of the abrasion model to the parameter bc ............................... A-7 B.6 Sensitivity of the abrasion model to the parameter f£ed............................... A-8 B.7 Sensitivity of the abrasion model to the mass balance adjustment factor . A-8 B.8 Sensitivity of the abrasion model to the parameter C .................................. A-9 B.9 Sensitivity of the abrasion model to the parameter p .................................. A-9 B.10 Sensitivity of the abrasion model to the parameter f l ed............................... A-10 B. 11 Sensitivity of the abrasion model to the parameter p r ............................... A -10 B.12 Sensitivity of the abrasion model to the sediment grain-size distribution . A-l 1 B. 13 Sensitivity of the quarrying model to the parameter Lq .............................. A -12 B.14 Sensitivity of the quarrying model to the parameter h ...............................A-13 B.15 Sensitivity of the quarrying model to the parameter m .............................. A -13 B.16 Sensitivity of the quarrying model to the parameter 50 kPa). Under its own weight, this deformation begins to occur approximately 50 m below the glacier surface [Cuffey and Paterson, 2010]. On the other hand, basal sliding occurs when gravity acting on the glacier exceeds friction at the bed causing the ice mass to slide across the bed. As the contact between ice and rock typically has a high coefficient of friction, sliding occurs predominately in the presence of meltwater which reduces the effective pressure at the bed [Eyles, 2006]. Temperate glaciers, as assumed here for Peyto, have a seasonal supply of meltwater allowing enhanced seasonal sliding to occur. The physics of internal deformation of ice is moderately well understood but predicting the timing and magnitude of basal sliding is far more uncertain [Cuffey and Paterson, 2010]. This uncertainty arises from the complex conditions at the glacier bed. There are several boundary conditions for the ice. At the surface, accumulation/melt must match the mass balance and the atmosphere exerts negligible stress. At the bed, ice does not penetrate the bedrock and there is friction between the ice and the bed. Conservation of mass, momentum, and energy must hold in all physical systems. Con­ servation of energy is not addressed here since ice is assumed to be temperate. Including thermodynamics would require a more sophisticated ice model that used in this study. As- 11 ( i n n ity Figure 2.5: Vector diagram for ice movement, suming that ice is incompressible (ie. density is constant), conservation of mass requires dXjVj = 0 (2.1) where dXj is the partial derivative with respect to Xj, v, [ma~x] is the velocity vector and j represents the Jc, y, and z directions. Conservation of momentum in the j direction yields Pid,Vj + p,dXk (vjVk) = dXkcrjk + fj (2.2) where rjk = crjk - p8jk [Pa] is the effective stress (total surface stress cr and pressure p), and fj [N m~3] is the force per unit volume from the body forces. Ice deformation is described by Glen’s Flow Law [Glen, 1955] which states that the 12 relation between strain rate, e [s '], and applied stress, r [Pa], follows a power law as: eij=AT*-xTij (2.3) where the exponent n » 3 is roughly constant, and A [Pa~n . r 1] depends mainly on temperature and ice impurities [Cuffey and Paterson, 2010]. For temperate ice (T « 0°C), A = 2.4 • 10"24 Pa~n s~l is assumed to be a constant. These values are well constrained through experiments and assumed constant in this study. The strain rate as defined above relates to fluid (ice) velocity derivatives as 1 dVj dxj dvj dx, (2.4) While the concept of sliding appears simple, mathematical descriptions are difficult due to complicated and largely unknown basal conditions. Deformation of subglacial sed­ iments can be incorporated into the description of sliding velocity [Schoof, pers. comm. 2011]. Due to variable conditions at the glacier bed and the difficulty of observation, slip is largely described using empirical relations [Clarke, 2005]. Sliding velocity should depend on the shear stress at the bed, and approach zero as the shear stress approaches zero. The simplest such relation is given by v, = Cim (2.5) where C [Pa~m a-1] is a constant dependent on bed roughness and the thermal and me­ chanical properties of the ice, r [Pa] is the basal shear stress, and m is a positive exponent. This equation was derived by Weertman [1964] using m = \{n + 1) where n v 3 is the Glen’s Flow Law exponent. Field experiments have shown that different glaciers require different exponents to successfully describe ice motion. No choice of exponent adequately 13 describes the variations observed; any positive exponent is generally accepted for model­ ing purposes [Cuffey and Paterson, 2010]. Different exponents show dramatically different flow dynamics, so the value will be chosen which best simulates observed conditions at the study site. The approach in equation 2.5 assumes no dependence on the presence of meltwater at the bed. Since the sliding velocity depends only on the shear stress, sliding is more or less constant throughout the year. Basal water pressure should impact sliding speed. As water pressure rises, the effective pressure at the bed Pe [Pa] (difference between ice pressure P, and water pressure Pw) decreases so sliding speed should increase. Thus the relation kl = C W P ^ (2.6) should be expected [Cuffey and Paterson, 2010]. Generally m = n ^ 3 and q = 1 are assumed based on simplified calculations. The relation, however, applies for any positive values of m and q since no optimal set of exponents has been identified through observa­ tions [Raymond, 1987]. Caution must be used with Equation 2.6 as it predicts an infinite increase in sliding velocity as the ice approaches flotation (Pe = 0). It is important to apply a critical effective pressure which yields the maximum sliding velocity. Note that by taking q = 0, Weertman’s sliding equation is obtained. The shallow ice approximation will be used to simplify governing equations. The shal­ low ice approximation assumes that the glacier length is much larger than the ice thick­ ness [Hutter, 1983]. Glacier thicknesses for Peyto are generally 100 - 200 m, while the glacier length is on the order of kilometers. Writing the boundary conditions and conserva­ tion equations as dimensionless allows identification of which terms will dominate when thickness is assumed to be much smaller than extent. Removing these terms allows sim­ 14 plification of the governing equations. Vertically integrating equation 2.1 and the surface boundary condition yields the continuity equation of ice (2.7) where H [m\ is the ice thickness, Q[ m2 a '] is the vertically-integrated ice discharge per unit width and b [m a-1] is the mass balance. Conservation of energy is not addressed since I assume that ice is temperate (0 °C). Using the velocity from equations 2.3-2.6 and vertically integrating yields the discharge ( 2 .8 ) where the stress is calculated by r = pigH |VS (2.9) t , = P i g H d XiS ( 2 . 10) where r, denotes the x, component of the driving stress r and S [m] is the surface elevation [Cuffey and Paterson, 2010]. The term in brackets from equation 2.8 is known as the diffusion coefficient. 2.4.2 Existing ice model The regional glaciation model developed by Christian Schoof and Garry Clarke of the UBC Glaciology Group is used for this study. Initially coded in Matlab, a Python version coded by Brian Menounos is used here. The vertically-integrated model uses the shallow 15 ice approximation. It is governed by equations 2.7-2.10 with Weertman sliding (q = 0 in equation 2.8). The model inputs include subglacial topography and mass balance gradient. Cell size and timestep are user-defined, depending on the application. The model tracks ice thickness, flux and velocity of each cell through time. The governing equations are solved by a finite difference approach. For the cell (/', j), diffusion coefficients of the discharge are calculated for (i + j), (i - -, j), (/, j + j), (i, j - \). Diffusion coefficients are converted to a sparse matrix in coordinate format. A conjugate gradient solver [Nazareth, 2009] is used to solve the sparse linear system. Parameter choices will be discussed in the model optimization process, section 2.4.4. 2.4.3 Glacial hydrology Since erosion depends strongly on the glacier sliding speed, it is desired to adapt the model from Weertman sliding to allow for the seasonal variability of sliding. Some representa­ tion of subglacial hydrology is required to determine the water pressure in equation 2.6. A complete hydrology model has been proposed by Flowers and Clarke [2002], which includes four modes of water transport (surficial runoff, englacial transport and storage, subglacial water networks, and subglacial aquifers) and exchanges between each. Fur­ ther developments have coupled this with a higher-order ice dynamics model [Pimentel and Flowers, 2010]. Unfortunately since water flows on a much faster timescale than ice, computing time suffers greatly by introducing a 4-part coupled hydrology model [Herman et al., 2011]. Herman et al. [2011] suggests incorporating only the subglacial component of the hy­ drology into a glaciation model. While exchanges between the systems are critical to represent the hydrologic system at a given time, subglacial hydrology is the critical com- 16 ponent for glacier mechanics. The subglacial hydrology is governed by ( 2 . 11) (2 . 12) where hw [m] is the water depth, Qj [m2 s~]] is the discharge rate per unit width, bh [m a~[] is the basal water source term, a = 3.5 is an empirical exponent, and hc [m] is a critical water depth (typically 0.1 m) [Flowers and Clarke, 2002]. The basal water source term consists of basal melt as well as surface meltwater that reaches the bed. The discharge, Q, is calculated by (2.13) (log (Kmax) - log (Kmin)) arctan kt (2.14) + j (log (Kmax) + log (/rmij ) where ka, kb, Kmax, and Kmin are all hydraulic conductivity tuning parameters. I adopt the values from Herman et al. [2011] due to a lack of field data for site optimization (values listed in section 2.4.4). Even including one component of the glacial hydrology remains computationally ex­ pensive due to the different timescales of flow. This difficulty provides the motivation to test a quasi-hydrological approach. The quasi-hydrological approach aims to improve upon Weertman sliding without suffering the computational slow-down of a hydrologic model. A simple assumption is that the amount of water available to the bed would be linearly dependent upon the surface melt rate given by the mass balance. Only water pres­ 17 sure, not discharge, is needed for the basal sliding velocity. Using the previous subglacial hydrology model for motivation, water pressure is estimated by (2.15) where bc [m a 1] is a critical mass balance parameter (which requires calibration) replacing the critical water depth hc. Not all surface melt will reach the bed, and this proportion is uncertain. Clason et al. [2012] suggests this proportion is relatively constant, and this proportionality is included in the bc parameter which will be optimized. Water pressure is taken to be zero if b is positive. To prevent numerical instability near flotation, water pressure is capped at Pw = 0.95P;. This approach simulates the annual variability in wa­ ter pressure through the variability of the mass balance. While water volume is expected to increase non-linearly down the glacier, summer melt in the ablation area reaches the threshold so little difference is seen. Since drainage system development is not included, this approach could overestimate water pressure in areas of high melt. This simplification of the hydrology is far from perfect, but aims only to be an improvement upon Weertman sliding without the substantial computational expense of a full hydrological model. 2.4.4 Model optimization Optimization of the glaciation model for a study area is required to determine values for some hydrology and sliding parameters. Optimization involves choosing free parameters that best simulate observations (ice velocity, extent, and thickness). Peyto Glacier will be used here, in preparation for the case study in chapter 3. Before considering the free parameters, the value of many parameters can be con­ strained by their physical meaning. Gravitational acceleration and ice density are taken to be constant in this study. Glen’s flow law is reasonably well understood and its parameters are assumed constant for a constant ice temperature near 0°C [Cuffey and Paterson, 2010]. Bed topography is obtained from UBC’s subglacial topography model which estimates bed elevation through surface inversion techniques [Clarke et al., 2013]. Radar surveys of ice depth have been used to verify this model in numerous locations. Output from the subglacial inversion work [Clarke et al., 2013] was a 200 m digital elevation model (DEM) which was interpolated to 50 m by kriging. A 200 m grid is too coarse for a glacier the size of Peyto Glacier, so the model was run with a spatial resolution of 50 m. The effect of using such a small grid size with the shallow ice approximation will be addressed in a future study. A smaller grid would push the limits of computing resources as well as the topography interpolation. A timestep of a year or more does not capture high melt events which lead to seasonal basal sliding and negate the main benefit of estimating subglacial water pressure. A one month timestep is chosen to provide a smooth variation in mass bal­ ance through the year. Lacking field measurements, hydrology parameters are assumed to be constant (Table 2.1). Finally, surface mass balance from 1966-2009 is used for climate input. Model calibration Measurable observations to use for optimization are ice extent, ice thickness, and ice surface speed. Ice extents are available for numerous years via digitization of air photos. The ablation region is used for extent comparisons, since this is the region of the most dramatic changes. Several ice thickness measurements were made for Peyto Glacier be­ tween 1983-1986. Ice surface velocity measurements were made in the 1940’s to 1960’s by measuring plates, but were limited to the glacier snout [Ommanney, 1972]. These ve­ locities ranged from 4 m a ' 1 to 12 m a ' 1 for a given year. Air photo pairs from 1982, 19 Table 2.1: Ice model (left) and hydrology (right) parameters Parameter Value n 3 .0 a A [Pa-3 a~l] Parameter Value h c [m ] 0 .1 b K min [ m a ~ x] 0.76 x 10"16 a Kmax [m a ~ l ] At [a] 0.083 c Ax [m] 50d g [m s~2] 9 .8 1 a Pi [kg m~3] 910a 1 x 10~3 b 1 x 10"6 b ka 15.0 b h 0.85 b p w [kg tfT3] 1000 b G [W m ~ 2] 0.05 a L [kJ kg-1] 334a “Value from [Cuffey and Paterson, 2010] *Value from [Flowers and Clarke, 2002] “Equivalent to one month J200 m used for subglacial hydrology model configuration 1986, and 1991 were analyzed to extract average surface velocities. Unfortunately, due to the length of time between photos common features could not be found to determine surface velocities. Recent satellite images did not yield a good pair of cloud-free images. Ice thickness measurements [Demuth et al., 2006] and ice extent [Theo Mlynowski, pers. comm] were readily available for 1986. Optimizing the model to these data max­ imized the amount of observations at a single date. The observed mass balance curve is adjusted to yield a zero net mass balance for 1986. The adjusted curve is used to evolve the modeled ice to a steady state. This zero net mass balance adjustment, denoted mbadd, has been estimated to be 0.495 m greater than the 1966 - 1995 average mass balance [Demuth et al., 2006]. The steady state solution for a variety of parameter sets is then compared to the measured ice thickness and extent. Using a steady state to optimize parameters means 20 proper response time of the glacier to changes in mass balance may not be captured. 2.5 Erosion models 2.5.1 Process-driven erosion model While glacial erosion can occur via numerous processes, the dominant two glacier specific processes are abrasion and quarrying [Boulton, 1982]. Fluvial erosion can also be an im­ portant contributor [Koppes and Montgomery, 2009], but the model’s focus is on glacier specific processes. The erosion rate from other processes, such as adhesion and dissolu­ tion, is thought to be orders of magnitude less than abrasion and quarrying [Rabinowicz]. These other erosive processes are thus excluded from the current study. Abrasion is the process of debris scraping against either bedrock or other debris [Anderson and Anderson, 2010], Erosion of bedrock is considered primary erosion, while erosion of debris is re­ garded as secondary erosion since it is the breakdown of previously eroded particles. Since sediment is assumed to have a certain grain-size distribution and is not tracked through the model, only primary erosion is explicitly considered in the present study. Abrasion typically produces sediments such as silts and clays. Quarrying is the process by which fractured rock is removed from the bed by freeze-on (freezing of bedrock to glacier sole) [Anderson and Anderson, 2010]. Quarrying produces coarse sediment with grain size as large as boulders [Huddart and Stott, 2013]. There has been much debate over the relative importance of quarrying versus abrasion. For a smooth bed that is not easily fractured, or a glacier with substantial subglacial till, abrasion can be a dominant source of erosion. For a heavily fractured bed surface, how­ ever, the quarrying rate can greatly exceed the abrasion rate [Duhnforth et al., 2010]. In general, abrasion is thought to account for at most 40% of erosion but typically much less 21 Figure 2.6: Vector diagram of the abrasion process. [Menzies, 2002]. Bed deformation can also be a substantial contributor to erosion for tillbased glaciers [Boulton and Hindmarsh, 1987, Clarke, 1987]. This deformation cannot be included explicitly without sediment tracking. Erosion from bed deformation is an exam­ ple of secondary erosion [Cuffey and Alley, 1996], and the movement is incorporated into the sliding velocity [Schoof, pers. comm. 2011]. Abrasion Abrasion was first proposed by Tyndall in 1864 [Linton, 1963]. The influencing fac­ tors have been greatly debated since Tyndall, though the existence of the process never experienced real opposition. The obvious signs after glacier retreat include striations in the direction of glacier movement by large debris and polishing by small debris [Menzies, 2002]. Two contrasting theories were finally developed in the 1970’s by Boulton [1974] 22 and Hallet [1979] to model the abrasion rate a [m a ']. Both take the general form: k yjh y F W V n a= hs A (2.16) where kabr is dependent on debris hardness and concentration, FN [kg m .T2] is the nor­ mal force at the bed, vp [m a~l] is the debris velocity, hs [Pa] is the bedrock hardness, and A [m2] is the unit area [Hildes, 2001]. The bedrock hardness and unit area will be included in the parameter kabr [Par1 m 2] henceforth. The difference between the theories lies in the estimation of FN. For erosion, it is the normal force between debris and bedrock that is of importance, rather than between ice and bedrock. Boulton’s theory assumes the ice pressure is entirely on the debris; it works well for very thin or cold ice [Boulton, 1982]. Hallet’s theory assumes dependence on vertical velocity and debris buoyancy and is suit­ able for temperate ice which can flow around debris by regelation and/or ice deformation. I will use Hallet’s theory since it is more widely accepted for warm-based temperate glaciers like Peyto [Hart and Boulton, 1991]. Approximating debris as spherical, the normal force is described by: 4 4nnR2 FN = ^ nR3 (Pr ~Pi) 8 COS 9 + fbed R2 + R2 V” (2‘1 where R [m] is the radius of the debris, p r and pt [kg m3] are the rock and ice densities, g [m s-2] is the acceleration of gravity, 0 is the ice surface gradient, f £ d is a bed modi­ fication factor determined to be 2.4 [Goran, 1970], tj [Pa s] is the effective viscosity of ice, v„ [m a -1] is the bed-normal ice velocity relative to the debris particle, and R, is the transition radius between regelation and creep [Hallet, 1979]. The first term represents buoyancy within the ice, and the second term represents viscous drag. Buoyancy is ne­ glected in modeling since it is small in comparison to viscous drag unless R is extremely 23 large [Hildes et al., 2004]. Ice movement normal to the bed is due to a combination of basal melting and ice sliding (with the effect of vertical strain negligible) giving v„ = bmelt + v* sin 9d (2.18) where bmel‘ is the basal melt rate, v* is the sliding velocity, and 9d is the local bed dip angle [Hildes et al., 2004], Note that v„ < 0 indicates debris moving away from the bedrock. Since no erosion results from movement away from the bed, v„ is taken as zero. Considering the force balance yields the debris velocity [Hildes et al., 2004] vp = vs (l - / i / " d sin 9d / f l d) - bmeh {nf£d / f l d) (2.19) where fi is the coefficient of friction and f£ed is the tangential bed modifier determined to be 1.7 [O’Neill, 1968]. Note that vp < 0 indicates deposition not movement up-glacier, and is taken to be vp = 0. The effective viscosity of the glacier can be found by l = \ [ A T n- 1]"’ ( 2 .20 ) with A, t and n previously defined [Cuffey and Paterson, 2010]. Finally, the transition radius R,, calculated by Watts [1974], is given by where Pc [AT Pa !] is the pressure melting coefficient, L [J kg *] is the latent heat of fusion, and Zr [m K VP”1] is the thermal resistivity of the debris. 24 Quarrying After initial debate over the existence of quarrying, it was widely accepted by the early 1900’s. However, debate continued until recently about the sustainability of quarrying [Lindstom, 1988]. Stress from the ice did not appear sufficient to fracture rock, but frac­ turing processes requiring subcritical stresses have now been found [Glasser and Warren, 1990, Hallet, 1996]. Quarrying is largely restricted to uneven bedrock regions. On the upvalley side of an asperity, high pressure results from the ice flowing downslope. This in­ creases melt from regelation which flows over or through cracksin the outcrop, refreezing on the down-valleyside where a cavity forms with much lower pressure. The refreezing can entrain rock fragments into the ice that have been sufficiently loosened by fractures [Liu et al., 2009]. The quarrying rate is governed by e = vc ( i - ^ ) where vc [m a~l] is the crack growth velocity and (2.22) is the cavity size S q [m] normalized by the distance between bumps Lq [m\ [Hallet, 1996]. Crack growth is expected to scale with stress intensity K, and is found to follow a power law relation vc = £Km (2.23) where m = 40 ± 10, and ( is a rock dependent constant. Hallet [1996] notes concern over this description due to the extreme sensitivity of crack growth to stress intensity. More recently, different approaches to calculating crack growth have been developed for modeling applications [Hildes, 2001, Iverson, 2012], but most follow the functional form of equation 2.22. Hildes’ method implements a combination of crack growth veloc- 25 Figure 2.7: Vector diagram of the quarrying process ities due to tension and shear propogation. While an improvement over Hallet’s descrip­ tion, this approach is still complicated and exhibits extreme sensitivity to stresses [Iverson, 2012]. Iverson [2012] introduces a simpler approach. While less explicitly derived, it does exhibit less sensitivity. For a bed with step-shaped bumps of height h and distance between bumps Lq, Iverson [2012] derives that the crack growth velocity is given by vskh 2ZT (2.24) where Jc is the probability of a bump fracturing. Further, the cavity size S q can be deter­ mined by i h \ i/2 S q = 4 A ~Vn ( i r ) {Pi “ Pwynl2 (2 '25) where A is from Glen’s Flow Law. For an assumed bed shape (h and Lq), the only factor 26 not determined in the ice model is the probability of breaking, k. Iverson estimates this as (2.26) where Vo N 2] is a characteristic rock volume per unit width, (rd [Pa] is the deviatoric stress in the bedrock, o-0 [Pa] is the Weibull scale parameter, and w is the Weibull modulus. The factor k ~ \ is introduced to account for the stresses required for slow crack growth over fast crack growth. The deviatoric stress is calculated by (2.27) where c is a scalar relating boundary to deviatoric stresses, and stress is capped at the stress where ice crushes, crn [Pa], Transport The erosion processes discussed so far estimate sediment production. Other than sub­ glacial deformation, the only methods of sediment transport are entrainment in the ice and movement by subglacial water. The entrainment of eroded material has been modeled pre­ viously, largely driven by regulation at the bed [Iverson, 1993, Alley et al., 1997, Iverson, 2000], However, it was deemed to be beyond the scope of this study to attempt inclusion of a sediment transport component to the model. With the inclusion of an entrainment model, sediment transport could be addressed in a future study. 27 2.5.2 Sliding erosion model Many modeling studies, particularly in landscape evolution, have used a single empirical mathematical equation to represent combined abrasion and quarrying [e.g. Braun et al., 1999, Jamieson et al., 2008], Erosion rate is driven by an ice dynamics quantity such as either sliding velocity or ice discharge. Using the more common sliding velocity, this equation takes the form E = Ki |v,| (2.28) where vs [m a -1] is the sliding velocity and K, is a constant that encompasses all other factors of erosion [Herman and Braun, 2008]. An exponent is sometimes included on the sliding velocity, but often assumed to be 1.0. There is no physical interpretation of Kt. The value is simply chosen to yield reasonable erosion rates in the basins it was tested on. The range of values found in the literature will be used here for comparison to the process model. However, I note that this is not a definite range and any number could legitimately be chosen. This approach to erosion modeling is particularly prevalent in landscape evo­ lution studies as it simplifies model mechanics to allow for shorter computation times. Comparisons of the computing times of the two models will be made. 2.6 Results 2.6.1 Model Optimization The first step was optimizing the parameter choices for the ice dynamics model. While velocities were not used in optimization (and measured velocities correspond to locations which are no longer ice covered in 1986), similarly low (< 2 0 m a~l) velocities are seen at near termini locations. Note that velocities are not this low in general. One region 28 0 0.5 1 1.5 2 km Figure 2.8: Map of modeled Peyto Glacier summer ice velocity. The region of high veloc­ ity (in red) is at the location of an icefall. near the equilibrium line altitude shows particularly high modeled summer velocities up to 200 m a-1 (Figure 2.8). Peyto has a steep icefall near this location which explains the anomalously high velocities. The free parameters of the model relate to the sliding law. These parameters (C, m and q) vary from glacier to glacier, but are commonly considered to remain constant through time for a particular glacier [Cuffey and Paterson, 2010]. Three configurations of the ice model are used to assess the impact of different approaches to estimating water pressure. Configuration 1 uses Weertman sliding, where hydrology has no influence on sliding (q = 0). Configuration 2 uses the subglacial hydrology as in Herman et al. [2011] to estimate 29 water pressure. Configuration 3 uses the quasi-hydrology approach (Equation 2.15) which has an additional free parameter, bc. Optimal values for these free parameters must be found for each model configuration. Optimization of free parameters for each model configuration proceeded similarly, fol­ lowing a uniform stratified sampling approach [Chaudhuri et al., 2007]. The input mass balance and topography remained constant. The subglacial hydrology configuration, how­ ever, required a grid size of 200 m to perform optimization runs with compute times in the same order of magnitude. In all cases, m and q were initially chosen, and then C (and bc for the quasi-hydrology) was varied to obtain an optimal solution compared to obser­ vations for the given exponents. Values of m ranged from 1.0 to 5.0, q ranged from 1.0 to 3.0, and bc ranged from 1.0 m to 10.0 m. Solutions diverge from observations towards the boundaries of these prescribed regions. C varies over many orders of magnitude depend­ ing on the values of m and q so its range is considerably large. Coarse sampling is initially taken across each parameter range. Subsequent rounds of finer sampling are completed in the regions of best performance. Values of individual parameters cannot be significantly changed from the optimized values without requiring changes to the other parameters as well. In all configurations that they appear, bc * 2.8 m, m ~ 3.0 and q ^ 1.0 were found to be the optimal combination while C values varied among the model configurations. Larger values of m resulted in thin ice, while smaller values resulted in thick ice, particularly in the accumulation area, for a given ice extent (Table 2.2). The subglacial hydrology configuration performed the best for ice thickness, but did not do well at estimating the glacier extent. The quasi-hydrology configuration outper­ formed the Weertman configuration in estimates of both extent and ice thickness, with low 30 computational expense. It also performed nearly as well as the hydrology configuration for ice thickness, and better for ice extent estimates. The subglacial hydrology configuration is expected to perform the best on both accounts if a 50 m grid was used, since it is the most physical representation. Unfortunately, using parameters determined on the 200 m grid does not transfer well to a model with a 50 m grid and computation times proved prohibitive on a 50 m grid. The computational cost (2-3 days for a 400 year simulation) is due to hydrology requiring a time step orders of magnitude shorter than ice. Due to the strong performance with limited computational demand of the quasi-hydrology model during the optimization phase, it will be used as the ice model of choice for use with the erosion models. 2.6.2 Process model Since the abrasion and quarrying components of the process model depend on different parameters, analysis of each will be undertaken separately. The sensitivity analysis can be completed by a variety of methods [Hamby, 1994, 1995], Here, it was completed by varying each parameter between the minimum and maximum likely physical values. All Table 2.2: Comparison of optimized ice model performance against observed glacier con­ ditions. The hydrology configuration used a 200 m pixel size, while the others employed a 50 m pixel size. Configuration Extent Difference {km2) Thickness Difference (m) Net Absolute Mean RMS Weertman 0.002 0.142 -11 33 Quasi 0.020 0.085 -7 29 Hydrology 0.04 0.2 -5 27 31 other parameters are fixed at their best estimated value from the literature. Models are run for a 2000 year period with varied climate, similar to the reconstruction that will be used at Peyto in the following chapter. Parameter sensitivities are partially dependent on the current ice dynamics, and sensitivities were not seen to transfer well between short and long sensitivity tests. The sensitivity of individual erosion parameters varies with changing ice dynamics. Sensitivities reported are thus averages of these various conditions. Abrasion component For the abrasion model, there are twelve parameters in equations 2.16-2.21: bc, C, fbed, fled' G, kabn mi P> Pn Zr, mbadd, and the subglacial grain-size distribution. Param­ eter values are chosen to be representative of a dolomitic basin. Values for Peyto are chosen as the middle of this dolomitic range, since specific details of the basin rock prop­ erties are lacking. Each parameter is varied individually over 10 levels across the range of dolomitic values, with the exception of grain-size (Table 2.3). The grain-size distribution was addressed separately to allow consideration of unimodal and trimodal distributions with varying means and standard deviations for each (Table 2.4). The grain-size distribu­ tions used are representative of previous estimates for various glacier basins [Bennett and Glasser, 2009, Hildes, 2001, Haldorsen, 1981]. Numerous levels were considered for each parameter to determine the linearity of erosion rate variations. Grain size distributions are reported as 0 values, with

C c 0) I 1000 X 500 500 1000 1500 2000 tea r (AD) Figure 3.7: Process model and sliding model sediment yield over 0-2000 compared to long-term average sediment yield of the lake and delta. Vertical line indicates 1917, the start of the 1917-2010 study period. The shaded regions represent ±1 standard deviation. the difference is not consistent, with ice volume leading sometimes (eg. end of the LIA) and ice volume lagging at other times (eg. MCA minimum). Many erosion studies discuss average erosion rate instead of sediment yield when com­ parisons to sediment deposits are not being made. Trends in erosion rate, while linked to sediment yield, do show differences during periods of extreme climate.The mean annual process erosion rate modeled here is 1.7 ± 0.3 mm. During both the Medieval climate anomaly and the late 1900’s, an increase in erosion rate is seen. The Medieval climate anomaly is a period of relative high temperature about 1000 AD. This corresponds to 59 3000 e — Process SY — Sliding SY — Ice Volume 1.4 2500 1.2 2000 ■o 5 1500 >c 01 £ I 1.0 1000 0.8 500 0.6 1500 500 2000 Figure 3.8: Sediment yield of the process and sliding models over 0-2000 compared to ice volume. The shaded regions represent ±1 standard deviation. when glacier area declines to near 10 km2. In both cases, the sediment yield continues to decline. An increase in sediment yield does not occur until the glacier retreat slows. Hold­ ing all parameters constant and varying only the sliding speed illustrates a sharp change from increasing to decreasing erosion rate which could be responsible for this negative feedback (Figure 3.9). Erosion rates do not increase as expected during the peak of the Little Ice Age. The increased sediment yield seen is due mainly to increased ice area rather than increased erosion rate. In this case, the increasing cavity size (which decreases quarrying rate) is 60 3.0 2.5 2.0 t K 1.5 g l.o o 0.5 o.o 10 20 30 Sliding speed (m a 1) 40 50 Figure 3.9: Relationship between quarrying rate and sliding speed with all parameters held constant. Curves represent the impact of changing the constant ice thickness used. dominant over the increasing sliding velocity. Given the overall correlation between ero­ sion rate and ice area, it appears that there is a transition between which factor (increasing ice velocity or increasing cavity size) is dominant both at small and large ice extents. Sliding model adjustment Since the sliding parameter value is not physically constrained, further adjustments are made to compare the models under similar mean sediment yields. To produce a mean sed­ iment yield over the past 2000 years equivalent to the process model, a sliding parameter Kj ~ 5.5 x 10~2 is required (Figure 3.10). Other than a sliding model lag, the models now 61 6000 Sliding Process 5000 J 4000 •< 2 3000 Figure A.2: Winter mass balance (in water equivalent) by elevation band from field obser­ vations taken between 1966-1995 with a polynomial best fit. -1 -6 2&00 2200 2400 2600 2800 3000 3200 Elevation (m) Figure A.3: Summer mass balance (in water equivalent) by elevation band from field observations taken between 1966-1995 with a polynomial best fit. A-4 B Sensitivity Analysis This appendix details time series plots for the sensitivity analysis performed on each pa­ rameter. The long-term average of each plot line is used in the sensitivity analysis section of Chapter 2. Ten levels are used across the range of parameter values detailed in Chapter 2. 05 0 .4 I c«cc £in 0.2 o ta 0.1 0.0 500 1000 1500 f e a r (AD) Figure B.l: Sensitivity of the abrasion model to the sliding law exponent m. The abrasion rate increases as m increases across 2.95-3.05 A-5 0,5 0 .4 0 .3 S cce 2 0 .2 : w* 2 01 0.0 500 1500 Figure B.2: Sensitivity of the abrasion model to the abrasive wear coefficient kabr. The abrasion rate increases as kabr increases across 5.0 x 10-1I-2.0 x 10-10 0 .5 0 .4 0 .3 a * UJ 0.1 0.0 500 1000 1500 \tear (AO) Figure B.3: Sensitivity of the abrasion model to the geothermal heat flux G. The abrasion rate increases as G increases across 0.0-0.1 W m~2 A-6 Figure B.4: Sensitivity of the abrasion model to the thermal resistivity of inclusion Zr. The abrasion rate increases as Zr increases across 0.15-0.63 m K W_1 0 .5 0 .4 |c 0 .3 a» c«c 1V) 02 O UJ 0.1 0.0 500 1000 \ f e a r (A D ) 1500 Figure B.5: Sensitivity of the abrasion model to the critical melt rate bc. The abrasion rate decreases as bc increases across 2.7-2.9 m a~] t e a r (AD) Figure B.6: Sensitivity of the abrasion model to the normal viscous bed drag factor f^ed. The abrasion rate increases as fE . increases across 22-2.6 0 .5 0 ,4 5| 03 2% t oC | 0.2 o (A UJ 0.1 0.0 500 10 0 0 t e a r (AD) 1500 Figure B.7: Sensitivity of the abrasion model to the net mass balance adjustment factor mbadd- The abrasion rate decreases as mbadd increases across 0.4-0.6 m a~l A-8 0.5 0 .4 0.0 u 1000 f e a r (AD) 1500 Figure B.8: Sensitivity of the abrasion model to the sliding law coefficient C. The abrasion rate increases as C increases across 2.25 x 10~9-2.65 x 10~9 Pa-3 a-1 1500 f e a r AD) Figure B.9: Sensitivity of the abrasion model to the coefficient of friction //. The abrasion rate decreases as ji increases across 0.60-0.85 A-9 0.5 0.2 0.1 0.0 500 1000 1500 t e a r (AD) Figure B.10: Sensitivity of the abrasion model to the tangential viscous bed drag factor fled- The abrasion rate increases as f l ed increases across 1.5-1.9 0 .5 0 .4 c 0 .3 a» cc aM 02 o UJ 0.1 0.0 500 1000 t e a r (AO) 1500 Figure B. 11: Sensitivity of the abrasion model to the bedrock density p r. The abrasion rate increases as p r increases across 2700-2900 kg m~3 A -10 0.5 500 1000 f e a r (A D ) 1500 Figure B.12: Sensitivity of the abrasion model to the subglacial sediment grain-size distri­ bution. The abrasion rate decreases as grain size $ values increase across -2.5-7.5. Note that as