MULTIPLE GEOCHRONOLOGICAL APPROACHES TO CONSTRAIN LATE QUATERNARY GLACIAL HISTORIES IN WESTERN CANADA by Adam Christopher Hawkins B.Sc., Western Washington University, 2014 THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN NATURAL RESOURCES ENVIRONMENTAL STUDIES UNIVERSITY OF NORTHERN BRITISH COLUMBIA December 2023 © Adam C. Hawkins, 2023 Abstract My dissertation investigates late Pleistocene and Holocene glacier change in western Canada and is presented in five chapters. In chapter one, I discuss the importance of the cryosphere, the techniques used for reconstructing past glacier behavior, previous research that has contributed to our current understanding of past glacier change, and the objectives and structure of my dissertation. In chapter two, I complete a multi-proxy investigation of Holocene glacier change at Gilbert Glacier, in the Southern Coast mountains. I use radiocarbon dating of lateral moraine stratigraphy and cosmogenic nuclide surface exposure dating to constrain the timing and duration of late Holocene advances at Gilbert Glacier to 2.0–1.8, ~1.5–1.3, ~0.9–0.8, and 0.4–0.1 ka. Organic matter associated with glacier advance within the north-lateral moraine at Gilbert Glacier records advances at 4.8–4.6, 4.5–4.3, 4.0–3.9, 3.8–3.6, 3.4, 3.2–2.9, 2.7, and 0.5–0.3 kilo calendar years BP (ka; 2-sigma age range). I advocate for the tandem use of multiple glacial geochronologic tools to better constrain the onset and termination of glacier advances. Chapter three applies cosmogenic surface exposure dating to a previously understudied region in the Mackenzie and Selwyn mountains of the eastern Yukon and Northwest Territories to develop a late Holocene glacier chronology. Surface exposure ages from 27 moraine boulders across nine glaciers show that glaciers reached their greatest Holocene extent between 1600-1850 CE. I use this glacier chronology to tune a glacier model forced by models of past climate to estimate regional changes in ice volume over the past millennium. Additionally, I use the same glacier model to estimate that glaciers in the region will decline in ice volume by 85% 97% by 2100 CE. Chapter four returns to the southern Coast mountains where I use cosmogenic dating on moraine boulders, bedrock surfaces, and shallow bedrock cores to investigate the history of ii deglaciation and Holocene glacier behavior of a small cirque glacier. Bedrock surfaces within and outside of the late Holocene maximum extent of the glacier record complex burial and exposure histories. I employ a Monte Carlo approach to evaluate the most likely glacier history that can be explained by our data. I discuss how our results compare with the work of previous chapters and past research, and the limitations of this technique to constrain the many variables that impact surface exposure ages. Chapter five provides a summary of my study’s major findings and further discusses their limitations and broader implications. I conclude with recommendations for future research that will expand on the work presented in this dissertation. iii TABLE OF CONTENTS Abstract-------------------------------------------------------------------------------------------------------- ii TABLE OF CONTENTS---------------------------------------------------------------------------------- iv LIST OF TABLES------------------------------------------------------------------------------------------vii LIST OF FIGURES--------------------------------------------------------------------------------------- viii Acknowledgements------------------------------------------------------------------------------------------ xi 1. Introduction------------------------------------------------------------------------------------------------ 1 1.1 The cryosphere and society------------------------------------------------------------------------- 1 1.2 Reconstructing past glacier behavior-------------------------------------------------------------- 3 1.3 Glacier response to Holocene climate in western Canada-------------------------------------10 1.4 Dissertation objectives----------------------------------------------------------------------------- 12 1.5 Dissertation structure------------------------------------------------------------------------------- 12 2. Tandem dating methods constrain late Holocene glacier advances, southern Coast Mountains, British Columbia---------------------------------------------------------------------------- 13 2.1 Abstract---------------------------------------------------------------------------------------------- 13 2.2 Introduction------------------------------------------------------------------------------------------14 2.3 Study Area------------------------------------------------------------------------------------------- 16 2.4 Methods---------------------------------------------------------------------------------------------- 17 2.4.1 Site selection----------------------------------------------------------------------------------17 2.4.2 Lateral moraine stratigraphy----------------------------------------------------------------17 2.4.3 Moraine exposure dating-------------------------------------------------------------------- 20 2.4.4 Regional radiocarbon record--------------------------------------------------------------- 21 2.5. Results----------------------------------------------------------------------------------------------- 21 2.5.1 Moraine morphology------------------------------------------------------------------------ 21 2.5.2 Stratigraphy, Gully A------------------------------------------------------------------------ 22 2.5.2.1 Unit descriptions---------------------------------------------------------------------- 22 2.5.2.2 Radiocarbon samples----------------------------------------------------------------- 24 2.5.4 Moraine exposure dating-------------------------------------------------------------------- 26 2.6. Discussion-------------------------------------------------------------------------------------------30 2.6.1 Neoglacial expansion------------------------------------------------------------------------ 30 2.6.2 Late Neoglacial advances------------------------------------------------------------------- 32 2.6.3 Implications----------------------------------------------------------------------------------- 35 2.7 Summary and concluding remarks--------------------------------------------------------------- 36 3. Late Holocene glacier and climate fluctuations in the Mackenzie and Selwyn Mountain Ranges, Northwest Canada-------------------------------------------------------------------------------38 3.1 Abstract---------------------------------------------------------------------------------------------- 38 3.2 Introduction------------------------------------------------------------------------------------------39 3.3 Study Area------------------------------------------------------------------------------------------- 41 iv 3.4 Methods---------------------------------------------------------------------------------------------- 42 3.4.1 Field site selection--------------------------------------------------------------------------- 43 3.4.2 Mapping of former and present glacier extents------------------------------------------ 43 3.4.3 10Be field sampling---------------------------------------------------------------------------45 3.4.4 10Be laboratory procedures and AMS measurements----------------------------------- 46 3.4.5 ELA reconstructions------------------------------------------------------------------------- 47 3.4.6 Glacier modeling----------------------------------------------------------------------------- 50 3.4.6.1 Open Global Glacier Model--------------------------------------------------------- 50 3.4.6.2 Equilibrium run------------------------------------------------------------------------50 3.4.6.3 Transient run--------------------------------------------------------------------------- 51 3.4.6.4 Future glacier simulations------------------------------------------------------------53 3.5 Results------------------------------------------------------------------------------------------------53 3.5.1 Glacier chronology---------------------------------------------------------------------------53 3.5.2 Climate reconstructions since the late Holocene-----------------------------------------55 3.5.3 Past millennium glacier change------------------------------------------------------------ 58 3.5.4 21st Century glacier projections----------------------------------------------------------- 59 3.6 Discussion------------------------------------------------------------------------------------------- 60 3.6.1 Holocene glacier fluctuations-------------------------------------------------------------- 60 3.6.2 ELA and climate reconstruction----------------------------------------------------------- 63 3.6.3 GCM evaluation------------------------------------------------------------------------------65 3.6.4 Future response of glaciers to climate change------------------------------------------- 67 3.7 Conclusions------------------------------------------------------------------------------------------67 4. Terrestrial Cosmogenic Nuclide depth profiles used to infer changes in Holocene glacier cover, Vintage Peak, Southern Coast Mountains, British Columbia-----------------------------69 4.1 Abstract---------------------------------------------------------------------------------------------- 69 4.2 Introduction------------------------------------------------------------------------------------------70 4.3 Study Area------------------------------------------------------------------------------------------- 72 4.4 Methods---------------------------------------------------------------------------------------------- 74 4.4.1 Field methods---------------------------------------------------------------------------------74 4.4.2 Laboratory processing----------------------------------------------------------------------- 76 4.4.3 Bedrock exposure history modeling------------------------------------------------------- 77 4.5 Results and Discussion-----------------------------------------------------------------------------79 4.5.1 Late Holocene moraine ages----------------------------------------------------------------79 4.5.2 Paired 14C and 10Be on erratic and bedrock----------------------------------------------- 85 4.5.3 14C bedrock depth profiles----------------------------------------------------------------- 88 4.5.4 Monte Carlo simulations-------------------------------------------------------------------- 89 4.5.5 Results synthesis and regional context---------------------------------------------------- 92 4.6 Conclusions------------------------------------------------------------------------------------------93 v 5. Conclusion------------------------------------------------------------------------------------------------ 95 5.1 Summary of major findings----------------------------------------------------------------------- 95 5.1.1 Holocene glacier change at Gilbert Glacier---------------------------------------------- 95 5.1.2 Mackenzie and Selwyn mountains---------------------------------------------------------96 5.1.3 Vintage Peak---------------------------------------------------------------------------------- 96 5.2 Study limitations------------------------------------------------------------------------------------ 97 5.3 Broader Implications------------------------------------------------------------------------------- 99 5.4 Recommendations---------------------------------------------------------------------------------101 Bibliography----------------------------------------------------------------------------------------------- 103 Appendix----------------------------------------------------------------------------------------------------114 Appendix A. Author contributions------------------------------------------------------------------ 114 Appendix B. Chapter 2 Supplement-----------------------------------------------------------------115 B.1 TCN field sample collection---------------------------------------------------------------- 115 B.2 Cosmogenic laboratory analysis------------------------------------------------------------115 B.3 Regional radiocarbon compilation--------------------------------------------------------- 116 B.4 Electron microprobe analysis--------------------------------------------------------------- 116 Appendix C. Chapter 3 Supplement-----------------------------------------------------------------123 Appendix D. Chapter 4 Supplement---------------------------------------------------------------- 133 vi LIST OF TABLES 2.1 Radiocarbon ages from this study 26 2.2 Gilbert Glacier 10Be sampling information 29 3.1 10 Be sampling information for all boulders in this study 47 4.1 Vintage Peak 10Be sample information 84 B.1 Cosmogenic nuclide input data 115 B.2 Cosmogenic nuclide sample blank information 117 B.3 Radiocarbon ages from this study 117 B.4 Regional radiocarbon compilation 118 C.1 Airphoto and satellite information 123 C.2 Cosmogenic sample blank information 123 C.3 10 Be ages with snow correction 124 D.1 Vintage Peak sample blank information 144 D.2 Radiocarbon AMS results 145 vii LIST OF FIGURES 1.1 Applications for Quaternary dating methods 5 2.1 Gilbert Glacier study area map 17 2.2 North lateral moraines of Gilbert Glacier 20 2.3 Gilbert Glacier stratigraphic sections 27 2.4 Late Holocene glacier advance records 32 3.1 Study area map of 10Be sampling locations in YT/NWT 40 3.2 Glaciers from which 10Be samples were collected 44 3.3 Box and whisker plots of 10Be ages 55 3.4 Changes in ELA and estimated temperature change 56 3.5 Modeled ice volume change 59 3.6 Fractional glacier volume change until 2100 CE 60 4.1 Bedrock and boulder sampling locations on Vintage Peak 73 4.2 Bedrock core from site 17-VPC4 75 4.3 Lockie’s Glacier Holocene extent change 80 4.4 Comparison of Vintage moraine record with regional glacier records 81 4.5 Snow depth variations on Vintage Peak 83 4.6 Paired 14C-10Be diagram 86 4.7 Elevation-normalized 10Be and 14C concentrations in bedrock cores 88 4.8 Bedrock 14C profiles with depth 89 viii 4.9 Summary of Monte Carlo results 91 B.1 Oxide plots of regional tephra 121 B.2 Oxide plots of tephra sampled at Gilbert Glacier and Mt Mazama 122 C.1 CRU TS 4.06 annual surface temperature anomalies 125 C.2 Comparison between different glacier model outputs 125 C.3 Example of temperature and precipitation bias calibration results 126 C.4 Example of glacier volume and length evolution 126 C.5 Modeled glacier profiles during non-transient 1000-year run in OGGM 127 C.6 Glacier hypsometries from non-transient equilibrium experiments in OGGM 128 C.7 Looking SE at North Moraine Hill Glacier 129 C.8 Temp, precipitation, and modeled glacier evolution with CCSM4 GCM climate 129 C.9 Temp, precipitation, and modeled glacier evolution with MRI GCM climate 130 C.10 Temp, precipitation, and modeled glacier evolution with MPI GCM climate 130 C.11 Temp, precipitation, and modeled glacier evolution with MIROC GCM climate 131 C.12 Manually digitized glacier margins from aerial photos and satellite imagery 132 D.1 17-VPC1 core site 133 D.2 17-VPC2 core site 133 D.3 17-VPC3 core site 134 D4 17-VPC4 core site 134 D.5 17-VP-07 erratic boulder 135 ix D.6 Example late Holocene moraine boulder (17-VP-04) 135 D.7 Monte Carlo runs within measured nuclide concentrations 136 D.8 Monte Carlo results for Core 1, no inheritance 137 D.9 Monte Carlo results for Core 2, no inheritance 138 D.10 Monte Carlo results for Core 3, no inheritance 139 D.11 Monte Carlo results for Core 4, no inheritance 140 D.12 Monte Carlo results for Core 1, 3 kyr inheritance 141 D.13 Monte Carlo results for Core 2, 3 kyr inheritance 142 D.14 Monte Carlo results for Core 3, 3 kyr inheritance 143 D.15 Monte Carlo results for Core 4, 3 kyr inheritance 144 x Acknowledgements I am grateful for the guidance, support, and care so many people gave me, which was critical to the success of my PhD. I first want to thank my supervisor, Brian Menounos. Brian gave generously his ideas and time to greatly improve the quality of my research and my development as a scientist. His mentorship pushed me to consider the importance and broader implications of my research. My writing has been heavily influenced by Brian’s teaching, though I continue to sneak in the odd nominalization to keep him on his toes. I became interested in glacial geology first for the opportunity to spend time in the mountains, immersed in the raw and powerful beauty of the alpine. Brian’s research support has allowed me to work in absolutely breathtaking locations in western Canada and I will forever be grateful for the opportunity to visit these places. Lastly, I want to thank Brian for his patience and motivation as I finished my PhD. This has been a long, faltering, and at times uncertain journey. Brian was stalwart in his support as I slowly made progress towards finishing my dissertation. I have also benefited from the experience and guidance from my supervisory committee: Brent Goehring, Gerald (Jerry) Osborn, and Peter Jackson. I want to especially thank Brent for his exceptional mentorship and hospitality during my laboratory work at Tulane University in New Orleans. He brought me into his lab and his home while I worked there and I am forever grateful for his kindness. I am honored to have had the chance to work with the likes of Jerry Osborn out in the field. His experience, wit, and good humor made fieldwork all the more enjoyable and memorable. This dissertation was financially supported by grants from UNBC, the Geological Society of America, the Hakai Institute (Tula Foundation), the National Sciences and Engineering Research Council (NSERC) of Canada, and the Canada Research Chairs program. I am deeply grateful for the opportunity to research and live on the traditional lands of many First Nations in xi western Canada, including those of the Lheidli T’enneh, Denendeh, Klahoose, and Tla’amin, among many others I have been lucky to visit. Interestingly enough, people don’t seem to say “no” when you offer them the chance of a free helicopter ride into a beautiful section of mountains where someone else (me) cooks for them. Little did they know I would make them hike up precariously steep moraines with concrete saws or have them scrambling around on icy bedrock slopes while gathering bedrock cores. My many thanks to those who helped with field work on this project: Rebecca Lerch, Chris Darvill, Jane Markin, Hunter Gleason, and (the) John Clague. I also want to thank Keir Nichols, Rachel Sortor, and the rest of the Tulane University Cosmogenic Nuclide Laboratory crew for their mentorship and friendship during my time there. Though I do not miss the cockroaches, I do miss you fine folks (and Friday free donuts). I also give my heartfelt appreciation for our research group at UNBC, whose knowledge and interest in the cryosphere was a steady source of inspiration. Many thanks to Ben Pelto, Alex Bevington, Caleb Mathias, Charlie Bourque, Maria Kroeger, and many others for their friendship. An unorthodox but important acknowledgement is to the volunteer Search and Rescue community, which has been a vital outlet for me during my graduate studies. This group of outstanding individuals has taught me how to be a better leader, mentor, and friend. I am so proud of the hard work these people do in training and response to help those in need. A few sentences can’t begin to address the foundational support my parents and siblings have provided that has gotten me where I am today. Suffice it to say, I love you all and yes, I am finally done! Last but certainly not least, I thank my partner, Jane, for her patience and support throughout my studies. Our little family with “the kids” Barrik, Nahanni, and James is a refuge from the anxiety of graduate school. Thank you for all of your love. xii 1. Introduction 1.1 The cryosphere and society The cryosphere includes seasonal snow, mountain glaciers, ice sheets, permafrost, and ice formed on lakes, rivers, and oceans. Snow and ice influence Earth’s climate by reflecting much of the high latitude incoming solar short-wave radiation back into space and cooling the surrounding landscape, thereby helping to regulate the surface temperature of the Earth. At longer timescales, changes in the extent of ice sheets influence global moisture and heat transport. Anthropogenic warming over the past several decades increased the average temperature in mountain areas in western North America, European Alps, and High Mountain Asia by about 0.3 °C per decade (Masson-Delmotte et al., 2018). This warming contributed to negative mass balance of glaciers in western North America and an estimated 7.4 Gt annual ice loss between 2000 and 2019 (Hugonnet et al., 2021; Menounos et al., 2019). Globally, glaciers are expected to lose 26 - 41% of their ice mass by 2100, compared to 2015 levels (Rounce et al., 2023). Nearly 2 billion people live downstream of mountain snowpacks and glaciers, and rely on runoff for consumption (Immerzeel et al., 2020). Glacier runoff is used for hydropower and irrigation for agriculture, especially during the summer/dry seasons (Biemans et al., 2019; Carey et al., 2017). Diminishing glacierized area changes the amount of meltwater contributed to streams (Moore et al., 2020). Though some recent global ice thickness estimates show that mountain glaciers may contribute less to sea level rise than previously thought, many mountain communities may have less time to adapt to future change (Millan et al., 2022). As glaciers retreat, they debustress and expose oversteepened valley walls, leading to an increase in mass 1 movement events and other glacier-related hazards (Geertsema et al., 2022; Hock et al., 2019; Moore et al., 2009). Glaciers are enigmatic features in mountain landscapes with deep cultural significance to communities that live around them (Jurt et al., 2015). Rapidly melting glaciers in New Zealand pose challenges to the glacier tourism industry, making access increasingly difficult (Purdie, 2013). Yet, in some locations, glacier retreat exposes new landscape elements such as lakes that provide additional tourism opportunities. Visiting retreating glaciers is an important component of “last chance tourism”, as the public travels to visit vanishing ecosystems (Lemelin et al., 2010). Glaciers represent important, abiotic features in glaciated mountain ecosystems. These ice masses serve as “water towers” to their downstream catchments, storing snow during the winter and releasing glacial meltwater during the summer. This meltwater helps to stabilize river temperature and discharge through the summer and early autumn months, which is vital to the health of aquatic ecosystems (Moore et al., 2020; O’Neel et al., 2015). Retreating glaciers in western Canada will have complex impacts on Pacific salmon populations, in some cases reducing viable habitat and in other areas creating new habitat (Pitman et al., 2020). Communities that rely on the influence of glaciers on their freshwater resources need to have an understanding of how continued glacier retreat will influence them. Our understanding of glacier response to ongoing climate change can be informed, in part, by how glaciers have responded to past climate perturbations. This dissertation aims to improve our understanding of past glacier and climate change in western Canada with several dating techniques that constrain the timing of past glacier fluctuations. This work gives greater 2 context to modern glacier retreat than previous research and it is used to improve estimates of future glacier fluctuations. 1.2 Reconstructing past glacier behavior Reconstructions of past glacier length, area or volume constrain the timing and magnitude of glacier fluctuations, which in turn can provide proxy evidence for past changes in climate (Oerlemans, 2005). Glacier mass balance is largely controlled by the local climate. For land-terminating glaciers, the primary drivers of mass balance are winter accumulation of snow and ice through snowfall and avalanching onto the glacier and mass loss during the summer months mostly by melting of snow and ice (Benn and Evans, 2014). As the climate cools and/or becomes wet, glacier mass balance becomes positive, eventually leading to glacier advance. Glaciers do not respond immediately to changes in climate, but have a lagged response time that can vary depending on glacier size and hypsometry, thermal regime, and previous mass balance history (Beedle et al., 2009; Cuffey and Paterson, 2010; Pelto and Hedlund, 2001). In contrast, increasing summer temperatures and/or less precipitation causes greater mass loss and subsequent glacier retreat. In order for a glacier reconstruction to be used as a paleoclimatic proxy, we must first develop a sufficiently resolved chronology of glacier fluctuations. Evidence of glacier fluctuations before photographs or eyewitness accounts includes landforms such as moraines or trimlines, proglacial lake sediments, and polished and eroded surfaces. Useful glacier chronologies provide tight constraints on the timing of an event, such as when a moraine was emplaced, or a glacier crossed a drainage divide, as well the duration a glacier was in an advanced or restricted position. Historical methods of “relative” dating, including schmidt-hammer tests, weathering rind thickness, pedogenesis, and moraine crest 3 morphology have largely been eschewed for absolute or numerical dating techniques, which assign specific ages to landforms or events (Haeberli et al., 2003; Schmidt, 1951; Wayne, 1984). The work presented in this dissertation uses numerical dating techniques to understand past glacier activity so that our results can be consistently compared across time and space. Below, I provide a brief synopsis of potential dating techniques pertinent to Holocene glacier studies (Fig. 1.1). 4 Figure 1.1: Cartoon diagram of applications for various Quaternary dating methods and the timelines overwhich the dating methods can be applied. Figure from (Davies, 2022). 5 Lichenometry is a biological dating method that uses the diameter of lichen thalli to estimate the age of a surface (Calkin and Ellis, 1980; Denton and Karlen, 1973; Osborn et al., 2015). The technique relies on researchers in the field measuring the largest observed lichen on rock surfaces (such as moraine boulders), correctly identifying the lichen species (most commonly Rhizocarpon Geographicum), and applying an independently-calibrated lichen growth curve to estimate lichen age. Lichenometry has been used throughout western Canada to date Holocene moraines, including around Nahanni National Park Reserve (Dyke, 1990), Yukon Territory (Denton and Karlén, 1977), in the southern Coast Mountains at Bridge Glacier (Allen and Smith, 2007), Manatee Valley (Koehler and Smith, 2011), Mount Waddington area (Larocque and Smith, 2004, 2003), and the Canadian Rockies (Luckman, 2007, 1977; Osborn and Taylor, 1975). Lichenometry is an inexpensive dating method, but it suffers from significant methodological shortcomings. As discussed thoroughly by Osborn et al. (2015), lichen misidentification, poorly constrained growth curves with limited ability to extrapolate to other areas, inconsistent treatment of errors, and likely erroneous assumptions of the preservation of lichens on rock surfaces combine to seriously question the utility of the method in Quaternary studies. Radiocarbon dating of organic material in formerly glaciated areas can be used to constrain past glacier advances (Luckman, 2007; Maurer et al., 2012; Ryder and Thomson, 1986; Yi et al., 2008). Materials useful for this dating technique include overridden stumps and trees, detrital wood fragments or preserved leaf litter. This organic material may be found within composite lateral moraines, within glacier forefields, or in proglacial lake sediments. The geomorphic position of materials collected for radiocarbon dating is crucial to its interpretation (Ryder and Thomson, 1986). For example, a stump in growth position, overlain by till, can be 6 interpreted with confidence to closely record the advance of a glacier. Organic material overlying till in a composite moraine, on the other hand, can provide a minimum limiting age of the underlying till. Dendrochronology can be used to date the timing of glacial advances and help constrain the duration of ice free periods (Luckman, 1988; Luckman et al., 1993; Schweingruber, 2012). Samples for dendrochronology can be collected from living trees growing on moraine surfaces to determine a minimum-limiting age for the emplacement of a moraine (Luckman, 2007). In some cases, detrital wood in glacier forefields can be cross-dated with living trees to develop a master chronology that may provide precise calendar dates for when the tree was killed (Lowe and Walker, 2014; Luckman and Masiokas, 2017; Reyes et al., 2006a). The age of the killed trees can provide a minimum duration of ice-free cover prior to ice advance. While dendrochronology, in combination with radiocarbon dating, can provide closely constrained glacier chronologies, suitable material for radiocarbon dating is not present in many formerly glaciated alpine areas in western Canada. In the Yukon and Northwest Territories, many forefields exist above treeline and so trees and shrubs were uncommon prior to late Holocene glacier advances. Additionally, the context of samples collected for radiocarbon dating is crucial, and often ignored in larger review studies (Menounos et al., 2009; Mood and Smith, 2015a). For example, the radiocarbon age of an overridden tree in growth position in a glacier forefield may date when a glacier was advancing at that extent on the forefield, but it does not indicate the timing of the glacier reaching its maximum extent. There is a need for a dating method that constrains the timing of specific events without the need for datable organic material and can be widely applied throughout high mountain regions. 7 When applied to glacier forefields and moraines, the dating methods described above provide discontinuous “snapshots” of glacier activity. Proglacial lacustrine sediments have the potential to provide a near-continuous record of sediment influx through the Holocene. Variations in glaciogenic sediment delivery to the lake can be used as a proxy for glacier size (Bakke et al., 2005; Karlen and Matthews, 1992; Maurer et al., 2012; Menounos et al., 2008; Røthe et al., 2018). Pro-glacial lake sediment records are most useful in cases where lake sedimentation is dominated by glacial influence. Non-glaciogenic sediment delivery can introduce noise to the sedimentary record and confound the glaciogenic signature. Lake sediments can also provide evidence of when a glacier expands beyond a flow divide; careful analysis of those sediments can demarcate times the receiving lake was fed by the glacier (e.g. Maurer et al., 2012). While pro-glacial lake sediments may provide relative glacier size information, the actual past size and extent of the glacier being studied can be difficult to determine. Finally, lake sediment records must be chronologically constrained. An age-depth relationship along the sediment core can be developed using identified tephra layers, terrestrial macrofossils dated with radiocarbon dating, and 210Pb and 137Cs dating. In areas without identifiable macrofossils or tephras, age control on sediment cores can be difficult to establish. Terrestrial cosmogenic surface exposure dating is a numerical dating technique that constrains the duration a rock surface has been exposed at or near the surface of the Earth and solves many of the limitations of the other dating methods described above (Balco, 2010; Dunai, 2010; Gosse and Phillips, 2001). High energy cosmic rays from deep space continuously bombard Earth’s surface at a rate that is essentially constant over time. High energy particles collide with atoms in the Earth’s atmosphere, creating a cascade of secondary cosmic-ray particles, primarily neutrons and muons (Dunai, 2010). The geomagnetic field and thickness of 8 the atmosphere above the Earth’s surface both influence the amount of secondary particles that reach the surface (Gosse and Phillips, 2001). The greatest number of secondary cosmic-rays reach the surface at high latitudes, where particles travel close to parallel with magnetic field flux lines, and at high elevations, where the atmosphere is thin. The bombardment of minerals within rocks at the Earth’s surface produces rare isotopes in the mineral lattice. These rare isotopes are created through several mechanisms, the most dominant being spallation by collision of a high-energy neutron with the target nucleus in an element, creating an isotope of a new element. Additional nuclide production occurs as a result of the interaction of low-energy muons. The most common rare isotopes used in cosmogenic dating are stable 3He and 21Ne and radioactive 10 Be, 14C, 26Al, and 36Cl. Nuclide production occurs in just the upper few meters of Earth’s surface, as the high-energy secondary cosmic rays are rapidly attenuated with depth. As depth below the surface increases, the relative production of nuclides by muons increases, as muons are less attenuated with depth compared to neutrons (Braucher et al., 2013; Dunai, 2010; Gosse and Phillips, 2001). A surface exposure age can be calculated with a modelled cosmogenic nuclide production rate, scaled to the latitude and elevation of the sample in question (Balco et al., 2008; Dunai, 2010). In cases of simple exposure, where the rock surface does not contain inherited cosmogenic nuclides from previous exposure, the measured isotope concentration allows for the exposure age to be directly calculated. In scenarios where the rock surface does have inherited nuclides, then the measured nuclide concentration is the sum of periods of accumulation, less any isotope loss through radioactive decay or erosion. Improvements in accelerator mass spectrometry (AMS) and isotope extraction techniques now allow samples as young as ~100 9 years to be dated. By using stable or long-lived isotopes, cosmogenic dating can be applied to surfaces throughout the Quaternary Period. Cosmogenic dating has been widely used in glacial environments to date the timing of moraine emplacement, bedrock erosion rates, and ice cover duration (Rand and Goehring, 2019; Schaefer et al., 2009; Young et al., 2021). While cosmogenic dating can answer varied questions about glacier history or earth-surface processes, the geomorphic context of each sample is required. A common dating scenario in glacial geology, and heavily featured in this dissertation, is constraining when boulders were emplaced on a moraine or forefield surface. An assumption in this approach is that the boulders on a moraine were sourced by subglacial plucking of the underlying bedrock surface, followed by transport to the lateral or frontal margin of the glacier. An additional assumption is that the bedrock surface has been either covered for a sufficient duration, allowing any inherited radionuclides to decay, or eroded deeply enough to remove any nuclide inventory from the bedrock. Once the boulder is deposited on the surface of the moraine, the boulder begins to accumulate cosmogenic isotopes. As described above, once sampled, the measured isotope concentration can be used to calculate an exposure age using a scaled production rate. This scenario represents a simple exposure history, with no prior inheritance. The measurement of multiple cosmogenic isotopes can elucidate complex periods of exposure and is discussed in greater detail in Chapter 4. 1.3 Glacier response to Holocene climate in western Canada Glaciers in western Canada and Alaska have a rich history of study since the mid-1900’s (Barclay et al., 2009; Calkin and Ellis, 1980; Mathews, 1951; Menounos et al., 2009; Mood and Smith, 2015a; Porter, 1989; Reyes et al., 2006b). The Cordilleran Ice Sheet underwent a complex 10 pattern of decay at the end of Pleistocene, punctuated by a brief period of advance or stillstand during the Younger Dryas cold interval (Darvill et al., 2022; Menounos et al., 2017). Following the Younger Dryas, most alpine glaciers retreated to minimal extents during the warm early Holocene (Menounos et al., 2009). Spatially limited evidence suggests glaciers may have begun small advances around 7 ka (Menounos et al., 2009, 2004). Largely in concert with decreasing summer insolation in the northern Hemisphere through the Holocene, glaciers periodically advanced multiple times since the mid-Holocene. Glaciers in the southern Coast and southern Interior mountains reached positions close to their Holocene maxima by at least 2 ka, with most glaciers reaching their greatest Holocene extents during the climatic advances of the Little Ice Age between 1300 to 1850 CE (Hawkins et al., 2021; Maurer et al., 2012; Solomina et al., 2015). In Alaska, glaciers advanced from retracted early Holocene positions around 4.5 ka. Glaciers later advanced at 3.3-2.9, 2.2-2.0 and 1.4-1.2 ka, culminating in two advances during the Little Ice Age at 1540-1710 CE and 1810-1880 CE (Barclay et al., 2009). While the majority of Holocene moraines in western North America date to the Little Ice Age (LIA), several exceptions exist. Mid-to-late Holocene terminal moraines lie beyond LIA moraines in northern and interior Alaska (Badding et al., 2013) and Tiedemann Glacier in the southern Coast Mountains (Menounos et al., 2013; Osborn and Luckman, 1988). Though there has been significant research on Holocene glacier fluctuations in western Canada, there are several knowledge gaps and limitations to our current understanding of Holocene glacier change. In the eastern Yukon and Northwest Territories, there has hitherto been comparatively little study of the region’s 1,235 glaciers (Pfeffer et al., 2014) likely due to the remoteness of the field sites (Menounos et al., 2009) and lack of organic material for dating. Some lichenometric ages of Holocene moraines and rock glaciers (Dyke, 1990) exist, and one 11 study inferred upvalley glacier fluctuations from proglacial sediments (Tomkins et al., 2008). Additionally, there is much less preserved evidence of Holocene glacier activity prior to the Little Ice Age period. There is a need for: i) robust dating methods to develop glacier chronologies in previously understudied regions of western Canada; ii) more techniques to elucidate climate and glacier change prior to the late Holocene maximum. 1.4 Dissertation objectives My dissertation seeks to improve our understanding of late Quaternary glacier and climate change in western Canada. To accomplish this goal, I completed field work, laboratory analysis and numerical modeling in the Yukon and Northwest Territories, as well as the southern Coast Mountains of British Columbia. My work focuses on answering the following questions: 1) How have glaciers in western Canada responded to climate change throughout the Holocene? 2) Can cosmogenic surface exposure dating be integrated with other existing dating methods to develop robust glacier chronologies? 3) Can cosmogenic surface exposure dating be used in novel ways to better understand glacier behavior since the late Pleistocene? 1.5 Dissertation structure This dissertation is paper-based and chapters 2, 3, and 4 represent contributions in peer-reviewed manuscripts. The chapters that have been published are presented with identical text as is in the published manuscript but may differ in formatting to conform the stylistic requirements required for the thesis. Chapter 5 provides a summary and conclusion of my dissertation, including a discussion of limitations within my studies and recommendations for further research. 12 2. Tandem dating methods constrain late Holocene glacier advances, southern Coast Mountains, British Columbia This chapter has been peer reviewed and is published in Quaternary Science Reviews. Please see Appendix A: Authorship Statements for details of the contributions of each author. Hawkins, A.C., Menounos, B., Goehring, B.M., Osborn, G.D., Clague, J.J., Jensen, B., 2021. Tandem dating methods constrain late Holocene glacier advances, southern Coast Mountains, British Columbia. Quat. Sci. Rev. 274, 107282. 2.1 Abstract Combined use of radiocarbon-dated subfossil wood within lateral moraines and surface exposure ages on moraine boulders provides an approach to better constrain times of glacier advance and onset of retreat. We test this method at Gilbert Glacier in the southern Coast Mountains of British Columbia where units of sediments associated with glacier expansion date to 4.8–4.6, 4.5–4.3, 4.0–3.9, 3.8–3.6, 3.4, 3.2–2.9, 2.7, and 0.5–0.3 kilo calendar years BP (ka; 2-sigma age range). Surface exposure (10Be) ages reveal times of moraine stabilization at 1.83–1.78, 1.38–1.28, 0.85–0.76, and 0.13–0.06 ka (interquartile range). Analysis of both datasets, as well as previously published regional advance records, narrows the age range of four late Holocene advances to 2.0–1.8, ~1.5–1.3, ~0.9–0.8, and 0.4–0.1 ka. We advocate for widespread use of our tandem approach at other sites throughout Earth’s high mountains to narrow the uncertainties associated with glacier expansion and better understand how glaciers respond to climate change. 2.2 Introduction Terrestrial evidence of Holocene glacier activity in western North America is temporally biased because most glaciers in this region underwent long-term, Holocene expansion, 13 destroying evidence of earlier activity (Menounos et al., 2009; Mood and Smith, 2015a; Solomina et al., 2015). Alpine glaciers in western Canada retreated to minimum extents during the early Holocene and then grew episodically until about 170 years ago (1850 CE) during the “Little Ice Age” (LIA). This general expansion was characterized by discrete periods of glacier advance separated by intervals of stabilization or glacier retreat (Clague et al., 2009; Koch et al., 2014; Menounos et al., 2009; Mood and Smith, 2015a). Many of these advances, however, remain poorly resolved. The timing of Holocene glacier advances in western Canada has traditionally been informed by radiocarbon dating subfossil plant material in lateral moraines and in glacially overridden forefields beyond present glacier termini (Menounos et al., 2009; Osborn et al., 2012; Osborn and Luckman, 1988). While useful, these ages provide only maximum limiting ages for glacier advances. When relying on a radiocarbon record alone, the duration of a prehistoric advance can only be roughly constrained by a sequence of radiocarbon ages of overridden subfossil wood associated with composite moraines. Dendrochronology can be used to precisely date trees damaged or killed at the distal limit of an advance (Koch, 2009), but sites containing such trees are rare (Luckman et al., 2020; Luckman and Masiokas, 2017). An exception is the work of Luckman et al. (2020), who combined composite moraine radiocarbon ages and dendrochronology on moraine crests to constrain the late Holocene glacier activity of several outlet glaciers of the Columbia Icefield in the Rocky Mountains of Alberta. Lichenometry has been shown to be problematic for a number of reasons, including the misuse of lichen growth curves, confounding lichen mortality, and a lack of standardized approaches when using the method (Osborn et al., 2015). 14 Another approach that has been widely used to numerically constrain times of Holocene glacier activity is cosmogenic surface exposure dating of moraine boulders (Kaplan et al., 2011; Schaefer et al., 2009). Exposure ages on moraine boulders provide minimum limiting ages for moraine stabilization and thus the termination of an advance. Cosmogenic surface exposure dating can be applied in areas lacking subfossil wood for radiocarbon dating and where no trees are available for dendrochronology. In this study, we used both surface exposure dating of moraine-crest boulders and radiocarbon ages from wood in a composite lateral moraine to constrain the Holocene activity of Gilbert Glacier in the southern Coast Mountains of British Columbia (Fig. 2.1). Combining these two dating methods offered us the opportunity to bracket the timing of a given glacier advance, an important prerequisite to more closely link climate drivers and glacier response. Previous studies used radiocarbon and surface exposure dating to develop glacier chronologies, but these studies employed sparse radiocarbon dating from lateral moraines or lacustrine sediments to develop regional glacier reconstructions (Ivy-Ochs et al., 2006; Kronig et al., 2018; Le Roy et al., 2017; Munroe and Laabs, 2017). To our knowledge, this is the first study to combine extensive surface exposure ages with detailed radiocarbon ages from a lateral moraine to narrow the age range for advances at a single location. 15 Figure 2.1: A: Map of southwest British Columbia showing the locations of Gilbert Glacier and other glaciers and icefields cited in this paper. B: Lidar-derived bare earth hillshade of Gilbert Glacier, with the modern (2019 CE) extent of Gilbert Glacier shown in teal. Lower left inset shows the extent of panel A within British Columbia, Canada. 2.3 Study Area Gilbert Glacier is a 9-km long valley glacier located in mainly granitic rocks of the southern Coast Mountains, 200 km northwest of Vancouver, British Columbia (Fig. 2.1). Large, composite lateral moraines flank both sides of the lower part of the glacier (Fig. 2.2), with the north lateral moraine being the more continuous and better exposed. This moraine was previously studied by Ryder and Thomson (1986). We examined this north lateral moraine over its full extent from 2100 to 900 meters above sea level (m a.s.l.). A series of nested minor moraine ridges is associated with the lower part of the north lateral moraine (Fig. 2.2). Ryder and 16 Thomson (1986) interpreted the three most proximal ridges to be latest Holocene (LIA) in age. They also described four units of till in an exposure through the north lateral moraine and interpreted them to record the late Pleistocene Fraser Glaciation (Marine Isotope Stage 2), Neoglacial expansion of Gilbert Glacier, and a Little Ice Age (LIA) advance. 2.4 Methods 2.4.1 Site selection We used satellite imagery and aerial photographs to locate exposures in the moraines and accessed the glacier forefield by helicopter. The moraine ridges sampled for surface exposure dating were located using air photos and ground-truthed in the field. An airborne lidar survey of the study area was completed while we were in the field and informed our geomorphic interpretations after fieldwork was completed. 2.4.2 Lateral moraine stratigraphy The proximal side of the north lateral moraine at Gilbert Glacier is incised by numerous gullies, exposing the stratigraphy of the composite moraine (Fig. 2.2). We accessed a gully, herein referred to as Gully “A” (Fig. 2.2A), documented the stratigraphy over a continuous vertical distance of 100 m, and collected wood samples from in situ overridden stumps, detrital logs, roots, and peat for radiocarbon dating. Samples of both detrital and in situ wood had preserved bark. Numerous wood mats contained large detrital logs (Fig. 2.3), but we could not safely access them. We determined stratigraphic unit and wood sample elevations with a laser range finder, using a lidar-derived elevation at the moraine crest as a datum. We also documented stratigraphy and collected radiocarbon samples at two smaller gullies down-valley of Gully A 17 (Gullies B and C, Fig. 2.2A). The stratigraphic sections were photographed, described in the field, and divided into lithostratigraphic zones using the procedure and nomenclature of Eyles et al. (1983). Wood mats containing a mass of detrital logs and roots, as well as significant organic horizons were counted as separate units. We submitted radiocarbon samples to the W.M. Keck Carbon Cycle Accelerator Mass Spectrometry Laboratory at the University of California Irvine. The samples were physically and chemically prepared with an acid-base-acid (1N HCl and 1N NaOH) treatment prior to combustion, and isotopic concentrations were measured in an accelerator to determine radiocarbon ages. We calibrated our radiocarbon ages and those reported by others using the Calib 8.1 software and the IntCal20 calibration curve (Reimer et al., 2020; Stuiver et al., 2021). An organic-rich sediment sample collected near the bottom of the Gully A moraine sequence contained glass shards representing a re-worked tephra. This tephra was identified by electron microprobe analysis at the University of Alberta to confirm provenance. The methods, microprobe data, and results are described in the Supplementary Materials accompanying this paper (Fig. B.1 and B.2). 18 Figure 2.2: North lateral moraine of Gilbert Glacier, showing dated boulder locations and median 10Be ages ± interquartile range for each moraine set. Panel A base imagery is from Google Earth; panel B is a base map created from hillshaded bare-earth lidar image acquired in 2019 and without any interpretive overlay; panel C shows the same section of moraine in panel B with additional interpretations. 19 2.4.3 Moraine exposure dating We collected 33 rock samples from granodiorite boulders resting on eight nested moraine crests (Fig. 2.2C; Table B.1) using the sampling methodology described by Darvill (2013). Large (generally >1 m3) boulders on the crests of the moraines were selected to reduce the chance of post-depositional movement and the impact of shielding by snow. We determined sample locations using a handheld GPS with approximate horizontal accuracy of 3-5 m and vertical accuracy of 5 m. We measured the dimensions of each boulder with a tape measure and took photos of each boulder from different aspects. At each site, we measured topographic and boulder self-shielding using a Brunton compass and inclinometer. We cut the surface of the sampled boulders with a concrete rock saw, then used a hammer and chisel to remove approximately 3x3x2.5 cm “brownies” of rock. We collected approximately 1.5 kg of rock from each boulder with consideration for the low abundance of quartz in the boulders and an assumed low isotope ratio due to the short exposure duration of the samples. All sub-samples from each boulder were placed into labeled cloth bags and transported in five-gallon buckets from the field. We processed the samples at Tulane University’s Cosmogenic Nuclide Laboratory, where the samples were physically and chemically purified to quartz and 10Be was extracted from each sample following standard procedures (Ditchburn and Whitehead, 1994). Three samples yielded insufficient quartz for further analysis. Pure BeO was loaded into cathodes and sent to Purdue University’s PRIME laboratory for AMS measurements, with isotope ratios normalized to standard 07KNSTD (Nishiizumi et al., 2007). We calculated exposure ages using version three of the University of Washington-hosted exposure age calculator, formerly known as CRONUS-Earth(Balco et al., 2008), with the default production rate dataset (Borchers et al., 20 2016). We report exposure ages using Lifton-Sato-Dunai (LSDn) scaling, assuming no significant snow cover or erosion (Lifton et al., 2014). Additional information used to calculate surface exposure ages can be found in the Supplementary Materials. 2.4.4 Regional radiocarbon record To provide regional context for our study at Gilbert Glacier, we compiled previously published radiocarbon ages on trees and stumps in growth position from glacier forefields in British Columbia, Canada, and Washington State (Table B.3). We restrict the regional record to sites within 350 km of Gilbert Glacier in order to extract a regional, coherent climate signal. The dataset comprises published radiocarbon ages from the Coast Mountains (Mood and Smith, 2015a), Mount Baker (Osborn et al., 2012), and the Cariboo Mountains (Maurer et al., 2012) that closely limit times of glacier advance (i.e. in situ stumps and logs with preserved bark resting on paleosols). 2.5. Results 2.5.1 Moraine morphology During its late Holocene climatic advance, Gilbert Glacier reached 4.5 km down-valley of its current (2019 CE) terminus. Airborne lidar-derived bare-earth hillshades reveal the complex structure of the north lateral moraine (Fig. 2.2B). The outermost dated moraine is a subdued, forested ridge with large (>3 m3) boulders along the sampled portion of the moraine crest. The intermediate moraines are forested, boulder-rich ridges, 5-10 m wide, that extend 100-200 m in length before terminating against the north valley wall slope or the eroded moraine wall slope. At the east end of the sample area, the north lateral moraine consists of three sparsely 21 treed ridges; two of these extend to the forested valley bottom (Fig. 2.2B). The terminal moraine at Gilbert Glacier is poorly preserved and was not sampled in this study. Moraine “sets” are defined by their surface exposure ages and geomorphic position (Fig. 2.2C). 2.5.2 Stratigraphy, Gully A Gully A is located about 600 m up valley from the site studied by Ryder and Thomson (1986). Recent erosion exposed six matrix-supported bouldery diamictons, separated by stratified sand and gravel units and wood mats, some of which contain stumps in growth position (Fig. 2.3). We subdivided the section into 22 lithostratigraphic units, numbered from the bottom of the gully exposure, about 105 m above the 2019 CE surface of Gilbert Glacier (Fig. 2.3; Fig. B.3). The location and material description of samples collected for radiocarbon dating are listed in Table 2.1. 2.5.2.1 Unit descriptions Units 1, 2, 4, 7, 16, and 22 are massive to weakly stratified, matrix-supported diamictons that differ from one another largely in terms of average or maximum grain size. Clasts range up to boulder size and are subangular to rounded; many of the clasts are striated and faceted. The matrix of these units consists primarily of sand, granules, and small pebbles. With the exception of unit 2, we interpret these units to be till on the basis of their lateral separation from the adjacent valley wall and the lack of valleyward-dipping stratification. Unit 2 is interpreted to be colluvium, due to its weak, valley-ward stratification. The lowest diamicton (unit 1) contains partially grüssified sub-angular to sub-rounded cobbles and boulders and is locally covered by a thin (~10 cm) layer of oxidized sand and silt. This unit overlies the steep bedrock slope of the 22 valley wall. The uppermost diamicton (unit 22) is the topmost unit in the section and forms the crest of the moraine, parallel to the valley axis. Units 5, 9, 11, 13, and 19 are horizontally stratified sandy pebble gravel units with some cobbles and boulders. Units 15, 17, and 21 are horizontally stratified sand units; the highest of these (unit 21) is backstopped by unit 22 till. We interpret the gravel and sand units to be ice-proximal glaciofluvial sediments deposited during periods of glacier thickening (Ryder and Thomson, 1986). Units 3, 6, 8, 10, 12, 14, 18, and 20 are wood mats consisting largely of detrital logs, branches, and roots, and, in some cases, in situ stumps and vertical tree stems in growth position. The wood mats range in thickness from a few tens of centimetres to nearly 5 m. Prominent wood mat layers are composed almost entirely of woody debris, but in some cases they grade laterally into organic-rich sandy silt. We interpret stumps in growth position within sand and gravel units to be the remains of trees killed when buried by glaciofluvial sediments in ice-proximal environments (Menounos et al., 2013; Ryder and Thomson, 1986). Wood mats beneath tills provide maximum limiting ages for times when the glacier thickened to the elevation of the mats. Most of the wood mats contain rooted material and are interpreted to be previously vegetated surfaces that were overrun as Gilbert Glacier advanced or were buried by glaciofluvial sediments. With exception of in situ stumps, we treat the radiocarbon ages from detrital branches, and roots, and logs as maximum limiting ages, as they may have been killed prior to ice expansion. 23 2.5.2.2 Radiocarbon samples Radiocarbon ages on 14 wood samples collected from Gully A range from about 5.60 to 0.43 ka (Table 2.1). The thick lowermost diamicton (unit 1) is likely till deposited during the Fraser Glaciation. Sandy peat at the top of unit 2 (Fig. 2.3) contains reworked Mazama ash, which has been dated to ~7.6 ka (Egan et al., 2015; Jensen et al., 2019). Because the tephra is reworked, its presence at the top of unit 2 provides no constraints on the deposition of unit 2. The unit 3 wood mat (Fig. 2.3) contains a detrital branch and ~3-cm thick root that returned radiocarbon ages of 4.52–4.42 ka and 5.60–5.49 ka, respectively. Unit 4 till is overlain by unit 5 sandy gravel (unit 5) containing two subfossil wood samples dated to 4.52–4.42 ka and 3.96–3.84 ka. Given the radiocarbon age on detrital wood from underlying unit 3, unit 4 till was emplaced 4.5 ka ago or later. Unit 7 till overlies a wood mat (unit 6) containing a detrital branch dated to 3.82–3.64 ka and is overlain by another wood mat (unit 8) containing a detrital log dated to 3.69–3.57 ka. A higher wood mat (unit 14) contains a stump in growth position that returned an age of 3.00–2.87 ka. It is overlain by a till (unit 16), which, in turn, is covered by sandy gravel (unit 17) containing stumps that were dated to about 2.7 ka. Still higher, a 3-5-m thick wood mat (unit 20) underlying the uppermost till (unit 22) yielded a radiocarbon age of 0.48–0.32 ka. In summary, maximum ages for periods of expansion of Gilbert Glacier based on the stratigraphy in Gully A are 5.3–5.0, 4.8–4.6, 4.5–4.3, 4.0–3.9, 3.8–3.6, 3.4–3.3, 3.2–2.9, 2.7, and 0.5–0.3 ka (2-sigma ranges). Wood in growth position with ages of 3.37–3.25, 3.00–2.87, and 2.76–2.72 ka provide limiting ages that are thought to closely date times when Gilbert Glacier 24 was advancing, whereas ages on detrital wood are less constraining, although still maximum-limiting ages. Table 2.1: Radiocarbon ages from this study. 2.5.3 Stratigraphy, gullies B and C Detrital wood in Gully B overlies a diamicton that, based on its physical properties (e.g. texture and inclusion of oxidized cobbles and boulders) and its stratigraphic position, correlates with the lowermost till exposed in Gully A (Fig. 2.3). Wood resting on top of this diamicton within an organic-rich layer yielded an age of 4.85–4.66 ka, similar to the radiocarbon ages above and below unit 4 in Gully A (Fig. 2.3). A detrital wood sample collected from Gully C 25 returned an age of 4.08–3.93 ka. The sand unit containing the detrital wood may correlate with unit 5 in Gully A, but the limited exposure at gullies B and C precludes definitive correlation. Figure 2.3: Stratigraphic sections in three gullies on the north lateral moraine of Gilbert Glacier. Lithostratigraphic units referred to in the text are labeled with bold numbers; wood mat units are labeled with brown numbers. Two sigma ranges of radiocarbon ages (ka) are shown in brackets. The photo at the right is a view towards the southwest of Gully A; note wood mats separated by till, gravel, and sand units. Unit widths in the stratigraphic section approximate lateral exposure, rather than average grain size. 2.5.4 Moraine exposure dating The 10Be surface exposure ages record four intervals of moraine abandonment, 1.83–1.78 (n=5), 1.38–1.28 (n=6), 0.85–0.76 (n=5), and 0.13–0.06 ka (n=14) (interquartile ranges). The 26 three youngest moraine sets are each composed of multiple crests, which we have grouped according to their surface exposure ages. The ages of the moraines are consistent with their geomorphic positions; that is, the oldest moraine crest is the most distal, and the youngest crest(s) is the most proximal to the edge of Gilbert Glacier (Fig. 2.2B). The three most proximal moraine crests, treated as a singular moraine set due to overlapping exposure ages, date to the classic LIA (Menounos et al., 2009). 27 Table 2.2: Gilbert Glacier 10Be sample information. 2.6. Discussion Gilbert Glacier provides an exceptional record of mid-to-latest Holocene glacier fluctuations, elucidated through the tandem use of radiocarbon and surface exposure dating. We demonstrate the tandem use of the stratigraphic and moraine surface exposure record to constrain two late Holocene advances of Gilbert Glacier. We then incorporate previously published radiocarbon data from nearby glaciers to place the entire surface exposure record in context. 2.6.1 Neoglacial expansion The radiocarbon record from the north lateral moraine of Gilbert Glacier details glacier expansion during the period 5–2 ka. The earliest advance for which we have evidence occurred after 5.6 ka (earliest possible age of unit 3, Gully A) and before 4.2 ka (unit 5). This advance correlates with advances of Castle Creek, Helm, and Franklin glaciers in British Columbia (Koch et al., 2007b; Maurer et al., 2012; Mood and Smith, 2015b). Several radiocarbon ages on detrital wood samples at Gilbert Glacier may provide maximum limiting ages for glacier thickening around 4.5–4.4 ka, similar to ages of a reported advance of Tiedemann Glacier (Menounos et al., 2008; Ryder and Thomson, 1986). Unit 7 till is underlain by a wood mat containing a branch that returned a radiocarbon age of 3.82–3.64 ka and is overlain by a wood mat containing a detrital log that yielded an age of 3.69–3.58 ka. Ice expansion around 3.7 ka coincides with an advance at Confederation Glacier (Coulthard et al., 2013). Unit 17 till is underlain by an in situ tree stem that dates to 2.99–2.89 ka (unit 14) and is overlain by three ages on the outer rings of in situ tree stems in growth position ranging from 2.76 to 2.72 ka (unit 17), coincident with advances of Lillooet, Bridge, Decker, and Tiedemann glaciers (Allen and Smith, 2007; Menounos et al., 2009; Osborn et al., 2007; Reyes and Clague, 2004). Beluga Glacier advanced shortly after 2.7 30 ka, near the end of this interval (Koehler and Smith, 2011; Menounos et al., 2009). Additionally, Castle Creek Glacier in the Cariboo Mountains was within several hundred metres of its LIA limit as early as 2.7 ka and fluctuated around this advanced position several times until its climatic LIA advance (Maurer et al., 2012). In their investigation at Gilbert Glacier, Ryder and Thomson (1986) reported a litter layer with an age of 2.11–1.88 ka. We use this age to bracket an advance of Gilbert Glacier described below (Fig. 2.4). The youngest radiocarbon age from Gully A (0.48–0.33 ka) is from a detrital log within a large wood mat (unit 20) underlying the uppermost till (unit 22). This log provides a maximum limiting age for the onset of the classic LIA advance of Gilbert Glacier, in agreement with other records throughout western Canada (e.g. Menounos et al., 2009; Mood and Smith, 2015a; Solomina et al., 2015). 31 Figure 2.4: Late Holocene glacier advance records for Gilbert Glacier (red line) and regional advance records from other sites in western Canada and northwest Washington State (blue line). Blue vertical bars indicate the times of regional advances (see text for details); the width of the bars are 2-sigma uncertainties from single dated samples or the published age range for the advance. The two horizontal black bars are the 2-sigma range of radiocarbon-dated material at Gilbert Glacier over this time interval. Red vertical bars are times of moraine stabilization and retreat of Gilbert Glacier, with bar widths set by the interquartile range of the 10Be dated moraine. Taken together, the pairs of advance and retreat dates constrain individual advances of Gilbert Glacier. The black temperature line is the mean of all temperature reconstructions for 30-60°N from Kaufman et al. (2020). Temperature anomaly (100 year bins) is relative to the 1800-1900 CE mean. 2.6.2 Late Neoglacial advances We use radiocarbon ages from this study (Table 2.1) and those from our compilation (Table B.2) to bracket the onset and termination of late Holocene advances at Gilbert Glacier. Of the four dated moraines at Gilbert Glacier, two are bracketed by correlative radiocarbon ages. We use the age of the litter layer from Ryder and Thomson (1986) mentioned above and the outermost surface exposure dated moraine at Gilbert Glacier to constrain an advance of 32 Gilbert Glacier commencing at ~2.0 ka and terminating by ~1.8 ka. Given the uncertainties in the radiocarbon and surface exposure ages, this advance may have lasted as long as 450 years, or as short as several decades. This advance accords with well dated advances of Bridge and Deming glaciers (Allen and Smith, 2007; Osborn et al., 2012). We tentatively correlate a group of moraine crests formed at Gilbert Glacier and dated by 10 Be surface exposure to ~1.3 ka (Fig. 2.2B) with the regional “First Millennium Advance”, documented at Lillooet, Bridge, Franklin, and Deming glaciers (Allen and Smith, 2007; Mood and Smith, 2015b; Osborn et al., 2012; Reyes et al., 2006b). If the onset of glacier expansion at Gilbert Glacier is synchronous with regional glacier activity, the duration of this advance is ca. 300 years, having initiated by ~1.6 ka and ended by 1.3 ka. A third advance is recorded by a group of moraine crests on the north lateral moraine with a median surface exposure age of 0.8 ka (1150 CE). A correlative advance, centered around 0.9 ka (1050 CE), has been documented at Klinaklini, Warren, Deming, and Castle Creek glaciers (Koch et al., 2007a; Maurer et al., 2012; Osborn et al., 2012; Ryder and Thomson, 1986). If Gilbert Glacier advanced at 0.9 ka, the advance lasted for about a century. The age of the final advance of Gilbert Glacier is bracketed by a radiocarbon age of 0.48–0.33 ka (1470–1620 CE) on wood from a large wood mat (unit 20) underlying the uppermost till (unit 22) in Gully A and surface exposure ages of 0.13–0.06 ka (1820–1890 CE) on three closely spaced moraines. The three closely spaced ridges represent the latest Holocene terminal and recessional positions of Gilbert Glacier. The regional glacier record shows two distinct intervals of latest Holocene glacier advance, one around 0.5 ka and a second at about 0.3 ka (Fig. 2.4). At Gilbert Glacier, however, we observe only one set of moraine crests within this period, taking into account uncertainties in the surface exposure ages (Fig. 2.2). Assuming 33 Gilbert Glacier also experienced two advances in this period, we attribute the lack of a second, older moraine crest to be the result of either poor moraine preservation or, more likely, to the second advance having been more extensive than the first, which would have destroyed or obscured evidence of the first advance. Tandem use of radiocarbon and surface exposure dating refines the timing of advances of Gilbert Glacier. Without the radiocarbon record, we would have no information on the mid-Holocene Neoglacial activity, and without the surface exposure ages, we would have no evidence of glacier activity between ~2.0 ka and ~0.4 ka. Additionally, we are able to closely constrain the ages of two local advances of Gilbert Glacier, while providing reasonable estimates of advance timing for two more moraine sets based on regional radiocarbon records. Terrestrial glacier records are imperfectly preserved, making it necessary to use multiple types of evidence to reconstruct past behavior. Kronig et al. (2018) employed surface exposure dating on moraine crests and radiocarbon dating on organic material within composite moraines on two glaciers in the Swiss Alps. However, in their study only one new radiocarbon age was reported from a wood sample in till. They used previously published radiocarbon ages from overridden paleosols within composite moraines as maximum limiting ages from times of glacier advance in order to place their 10Be exposure ages in context. Le Roy et al. (2017) similarly obtained surface exposure ages on Holocene moraines in the French Alps, while using previously published radiocarbon evidence to place their data into a regional context. Our study at Gilbert Glacier demonstrates the advantage of employing tandem dating methods at a single site. Ambiguities introduced by combining data from glaciers that may have different response times to climate forcing are explicitly avoided. 34 Developing a detailed glacier chronology provides greater insight into the influence of past climate on glacier activity. Our glacier chronology is in good agreement with Northern Hemisphere temperature reconstructions (Kaufman et al., 2020), notably with moraine abandonment and glacier retreat at the onset of century-scale intervals of warming (Fig. 2.4). Although both local temperature and precipitation controlled the behavior of Gilbert Glacier, our chronology suggests that Gilbert Glacier generally responded in concert with regional temperature variations. 2.6.3 Implications Our study refines understanding of the duration of late Holocene glacier advances in western Canada. Gilbert Glacier was at or near its maximum Holocene extent by 2.0 ka and reached this limit several times prior to the classic Little Ice Age. Gilbert Glacier is only the second site in western Canada where dated late Holocene moraines lie outside moraines constructed during the climatic advances of the LIA. The other locality is Tiedemann Glacier, although modelling suggests a large rock avalanche event onto the surface of the glacier likely contributed to its exceptional pre-LIA advance (Menounos et al., 2013). There are numerous other sites throughout western Canada where previous work indicates the presence of pre-LIA moraines (Allen and Smith, 2007; Larocque and Smith, 2003; Osborn et al., 2007; Reyes and Clague, 2004), but these moraine surfaces remain undated. Expanding spatial coverage of aerial lidar in high mountains (i.e. Johnson et al., 2015) is revealing previously unrecognized glacial landforms obscured by vegetation and thus providing many opportunities for research of the type reported here. Where moraine complexes are found, surface exposure dating and, if possible, dating of subfossil plant material will help to improve our understanding of past glacier fluctuations in this and other regions. 35 Comparison of radiocarbon-based records of glacier advances is typically difficult, as the ages are maxima or minima and relate only to the thickening and passage of a glacier at a specific site. This contextual information is commonly lost in regional comparisons. For example, a glacially overridden tree on the valley floor close to the modern terminus of a glacier may have died several hundred years before the death of an overridden tree at the distal limit of the advance or near the moraine crest, which might be 100 m or more above the valley floor. In such cases, the radiocarbon ages record the same advance, but might be interpreted to record two separate advances (e.g. Osborn et al., 2012). In contrast, moraine surface exposure ages more closely date the termination of an advance, providing additional information that complements a radiocarbon chronology. In situations where moraine boulders have little to no inheritance and minimal post-depositional movement, surface exposure ages allow for confident dating of a moraine. This method can be applied broadly and does not suffer from the limitations faced by lichenometry discussed in Osborn et al. (2015). Robust chronologies are needed to investigate regional and hemispheric climate and glacier linkages. Priority should be given to sites with composite moraines containing overridden organic material and lithologies suitable for cosmogenic exposure dating, for example Patagonia (Kaplan et al., 2016; Menounos et al., 2020), New Zealand (e.g. Schaefer et al., 2009), and sites throughout western North America, including those with dated composite moraines but poorly dated moraine crests (e.g. Allen and Smith, 2007; Barclay et al., 2009; Reyes and Clague, 2004). 2.7 Summary and concluding remarks Tandem radiocarbon dating of subfossil plant material in lateral moraines and exposure dating of boulders on moraine crests more tightly constrain Holocene glacier fluctuations than is 36 possible with either technique alone. We targeted Gilbert Glacier as a case study of our tandem dating technique for several reasons. First, the north lateral moraine has several gullies that expose a composite stratigraphy of tills and glaciofluvial sediments separated by wood mats. Second, the crest of the moraine is ornamented with numerous smaller morainal ridges that are suitable for surface exposure dating. Third, the local granitic bedrock is quartz-bearing and thus appropriate for 10Be exposure dating. And fourth, a foundational earlier study (Ryder and Thomson, 1986) showed the potential of the site for testing our method. Eighteen radiocarbon ages from the north lateral moraine of Gilbert Glacier show that the glacier thickened and advanced shortly before 4.8–.6, 4.5–4.3, 4.0–3.9, 3.8–3.6, 3.4, 3.2–.9, 2.7, and 0.5–0.3 ka. Thirty 10 Be exposure ages on nested morainal ridges place the onset of glacier retreat at 1.83–1.78, 1.38–1.28, 0.85–0.76, and 0.13–0.06 ka. We used these data to constrain several advances: an advance that began by 2.0 ka and terminated around 1.8 ka; another advance around 1.5 ka that terminated by 1.3 ka; a short-lived advance between 0.9 ka and 0.8 ka; and one of two Little Ice Age advances that began about 0.5 ka (1450 CE) and terminated about 150 years ago. Our study adds to the evidence that many glaciers in western Canada were at or near their LIA maximum extents by 2.0 ka or earlier, after which they reached similar positions several times prior to the advances of the past several centuries. Surface exposure dating is a useful tool for dating Holocene moraines, especially in areas where complex moraine morphologies are preserved. Further application of this tandem dating method will better constrain Holocene glacier chronologies in western Canada and other high mountain areas, leading to further insights into the influence of past and future climate on glacier fluctuations. 37 3. Late Holocene glacier and climate fluctuations in the Mackenzie and Selwyn Mountain Ranges, Northwest Canada This chapter has been peer reviewed and is published in The Cryosphere. Please see Appendix A: Authorship Statements for details of the contributions of each author. Hawkins, A.C., Menounos, B., Goehring, B.M., Osborn, G., Pelto, B.M., Darvill, C.M., Schaefer, J.M., 2023. Late Holocene glacier and climate fluctuations in the Mackenzie and Selwyn mountain ranges, northwestern Canada. The Cryosphere 17, 4381–4397. 3.1 Abstract Over the last century, northwestern Canada experienced some of the highest rates of tropospheric warming globally, which caused glaciers in the region to rapidly retreat. Our study seeks to extend the record of glacier fluctuations and assess climate drivers prior to the instrumental record in the Mackenzie and Selwyn Mountains of northwestern Canada. We collected 27 10Be surface exposure ages across nine cirque and valley glacier moraines to constrain the timing of their emplacement. Cirque and valley glaciers in this region reached their greatest Holocene extents in the latter half of the Little Ice Age (1600-1850 CE). Four erratics, 10-250 m distal from late Holocene moraines, yielded 10Be exposure ages of 10.9-11.6 ka, demonstrating that by ca. 11 ka, alpine glaciers were no more extensive than during the last several hundred years. Estimated temperature change obtained through reconstruction of equilibrium line altitudes show that since ca. 1850 CE, mean annual temperatures rose 0.2-2.3 °C. We use our glacier chronology and the Open Global Glacier Model (OGGM) to estimate that since 1000 CE, glaciers in this region reached a maximum total volume of 34-38 km3 between 1765-1855 CE and have lost nearly half their ice volume by 2019 CE. OGGM was unable to produce modeled glacier lengths that match the timing or magnitude of the maximum glacier extent indicated by the 10Be chronology. However, when applied to the entire Mackenzie and 38 Selwyn Mountain region, past-millennium OGGM simulations using the Max Planck Institute Earth System Model (MPI-ESM) and the Community Climate System Model 4 (CCSM4) yield late Holocene glacier volume change temporally consistent with our moraine and remote sensing record, while the Meteorological Research Institute Earth System Model 2 (MRI-ESM2) and the Model for Interdisciplinary Research on Climate (MIROC) fail to produce modeled glacier change consistent with our glacier chronology. Finally, OGGM forced by future climate projections under varying greenhouse gas emissions scenarios predict 85% to over 97% glacier volume loss by the end of the 21st century. The loss of glaciers from this region will have profound impacts to local ecosystems and communities that rely on meltwaters from glacierized catchments. 3.2 Introduction Between 1990-2020 CE, northwestern Canada warmed by 1.1 °C above the 1961-1990 CE average (Muñoz-Sabater, 2019, 2021), which contributed to the loss of an estimated 0.429 ± 0.232 km3 of ice in the Mackenzie and Selwyn Mountains of eastern Yukon and Northwest Territories between 2000 and 2020 CE (Fig. 3.1; Hugonnet et al., 2021). Glaciers in this region are clearly responding to recent climate warming, but proxy evidence of past climate change is scarce (Tomkins et al., 2008; Dyke, 1990). Reconstructions of when and how glaciers responded to past climate change provide one method for estimating paleoclimatic conditions, while also placing the rate of modern glacier change into a geologic context. 39 Figure 3.1: Study area map of 10Be sampling locations. Panel (a) is the temperature trend from ERA5land between 1950 and 2021 CE. Few glacier change studies exist for the Mackenzie and Selwyn Mountains as compared to other mountainous regions in SW Yukon, British Columbia, and Alaska. Previous Quaternary research in this region focused on Pleistocene glacial deposits and Holocene rock glaciers (i.e. Duk-Rodkin et al., 1996; Fritz et al., 2012; Menounos et al., 2017; Dyke, 1990). The remote location and related logistical challenges of conducting fieldwork in this area are likely reasons this region is underrepresented in Holocene climate reconstructions (e.g. Marcott et al., 2013). The timing and magnitude of the most extensive Holocene glacier expansion in the eastern Yukon and Northwest Territories, which places modern glacier retreat in context, remains uncertain. Research in northern and interior Alaska indicates that glaciers reached their 40 maximum Holocene extents around 3.0-2.0 ka (Badding et al., 2013) while nearly all glaciers in southern Alaska and western Canada reached their greatest Holocene positions around 1600-1850 CE, at the culmination of the Little Ice Age (LIA, ~1300-1850 CE) (Menounos et al., 2009; Barclay et al., 2009; Hawkins et al., 2021). The primary objectives of our study are to develop a Holocene glacier chronology in the Mackenzie and Selwyn mountains of eastern Yukon and Northwest Territories and use our glacier chronology to estimate changes in climate responsible for these glacier fluctuations. We then deepen our understanding of glacier activity in this area by estimating glacier volume change using multiple models of past climate to force a glacier flowline model. Finally, we briefly evaluate future glacier change in this region in response to various greenhouse gas emissions scenarios. 3.3 Study Area The Mackenzie and Selwyn ranges extend over 600 km from north of the Liard River in northwestern British Columbia to the Stewart River and northern extent of the Mackenzie Range in northern Yukon (Fig. 3.1). This region is covered by 650 km2 of ice from nearly 1200 glaciers situated among peaks that rise as high as 2952 m above sea level (Pfeffer et al., 2014). Bedrock consists of faulted and folded Paleozoic sedimentary rocks with Early Cretaceous granitic intrusions (Pfeffer et al., 2014; Cecile and Abbott, 1989). A portion of our study area is situated in the Nahanni (Nááts'ihch'oh) National Park Reserve, which was expanded in 2009 to >30,000 km2 (Demuth et al., 2014). Glacier runoff within the Nahanni National Park Reserve flows into the Liard River watershed which later joins the Mackenzie River, eventually draining north to the Beaufort Sea. Two of our nine field sites are located nearly 200 kilometers to the northwest of Nahanni National Park Reserve and are situated on or adjacent to the Keele Peak massif, which 41 is similarly composed of Early Cretaceous granitic rock. Meltwater from our study sites on and near the Keele Peak massif flows into the Stewart River, which flows west to the Yukon River and eventually to the Bearing Sea. The watersheds in our study area are culturally and ecologically important for the numerous First Nations communities who have lived on this land for millennia, including the Dënéndeh, Kaska Dena, and Na-Cho Nyak Dun First Nations, among others. 3.4 Methods Our glacier chronology originates from digitized glacier margins of aerial photos and satellite imagery and constraining the age of late Holocene moraines using cosmogenic 10Be surface exposure dating. Cosmogenic surface exposure dating relies on the accumulation of rare isotopes, in this case 10Be, in the bedrock surface during periods of exposure at or near the surface of the Earth (Gosse and Phillips, 2001). We use this chronology to estimate paleoclimate conditions in the late Holocene using several methods. First, we estimate past and present equilibrium line altitudes (ELA) using the maximum elevation of lateral moraines (MELM, LIA maximum only), toe-to-headwall altitude ratio (THAR), and accumulation area ratio (AAR) and infer changes in temperature and precipitation from estimated ELA changes (Braithwaite and Raper, 2009; Meier and Post, 1962; Ohmura and Boettcher, 2018). We then estimate the temperature decrease needed to grow glaciers to their late Holocene positions using a flowline glacier model. Additionally, we perturb monthly temperature and precipitation from several General Circulation Model (GCM) simulations of climate since 1000 CE to produce modeled glacier extents that most closely match the terrestrial and remotely sensed record (Taylor et al., 2012) before evaluating past modelled glacier volume change for all glaciers in the Mackenzie 42 and Selwyn mountains. Finally, we model future glacier change in this region under various Representative Concentration Pathways (RCPs; Moss et al., 2010). 3.4.1 Field site selection We selected sampling locations within the Mackenzie and Selwyn Mountain ranges using satellite imagery, aerial photos, and digital elevation data to identify purported late Holocene moraines. We consulted bedrock geologic maps of the area to locate sites that likely contained quartz-bearing lithologies suitable for 10Be surface exposure dating (hereafter 10Be dating), which was then confirmed in hand-samples in the field (Cecile and Abbott, 1989; Gordey, 1992). Helicopters and floatplanes during late summer in 2014, 2016, and 2017 provided access to the field sites. 3.4.2 Mapping of former and present glacier extents We manually digitized past glacier outlines for six of the nine glaciers sampled for 10Be dating. Those glaciers represent sites with multiple dated moraine boulders and morphologies better suited for glacier flowline modeling. It is the author’s understanding that only two of the glaciers included in this study, North Moraine Hill and Butterfly glaciers, have formal names. The remaining glaciers are referred to with informal names below. The resulting glaciers used in paleoclimate reconstructions are Anderson, Mordor, North Moraine Hill, Butterfly, Keele Peak, and Arrowhead glaciers (Fig. 3.2). We used imagery from airphotos between 1949 and the mid-1970’s CE and satellite imagery from 1985 CE, onward (Table C.1). Air photos represent digitally scanned negatives housed at the Canadian National Airphoto Library (NAPL). We georeferenced each airphoto by manually selecting 40-60 ground control points (GCPs) on the air photographs and high-resolution satellite imagery (e.g. large boulders, peaks, and ridges). We 43 subsequently performed a thin-plate spline transformation in GIS software (QGIS), visually inspecting the georeferenced image for any obvious distortions. Portions of glacier outlines further from GCPs have positional errors smaller than 20 m. Figure 3.2: Glaciers from which 10Be samples were collected. Sample locations are shown with green circles. Moraine crests are depicted as black dashed lines. Exposure ages ± analytical errors for individual boulders are in text boxes, with erratic boulders ages shown in italics. Grey insets show sampling sites at larger scale. Imagery is from PlanetLabs, acquired between July and August, 2021 and 2022. 44 We used Landsat 5, 7, and 8 satellite imagery to delineate glacier margins at roughly 5-10 year intervals from the mid-1980’s onward (Fig. C.12). To aid in the manual digitization, we made false color composites for each Landsat scene to highlight the glacier surface relative to non-glaciated terrain. The surfaces of most glacier termini are debris free, which facilitated glacier mapping. We mapped late Holocene glacier margins using high resolution satellite imagery from Mapbox and PlanetLabs to delineate glacier trimlines and moraine crests. In areas with cloud cover or snow-covered terrain, we used hillshades from ArcticDEM to help identify moraine ridges (Porter et al., 2018). 3.4.3 10Be field sampling We targeted samples from large (generally taller than 1 m), granitic boulders on or near moraine crests (Fig. 3.2). It is commonly assumed that large boulders on moraine crests are windswept such that snow cover is minimal, and their large size limits the chance of being previously covered by moraine material or moving following deposition (Heyman et al., 2016). Recent work by Tomkins and others (2021) provides evidence that sampling from the crests of moraines may not reduce the chance of geomorphic exposure age scatter, however at the time of sampling in this study, we followed the common practice of targeting boulders on moraine crests. Several erratic boulders directly overlying bedrock and distal to the moraine crests were sampled as well. We measured topographic shielding of the incoming cosmic ray flux and boulder self-shielding using a Brunton compass and inclinometer, and then determined the location and elevation of each sample with a handheld GPS receiver with barometric altimeter. Samples were collected from the top surfaces of boulders using a concrete saw and hammer and chisel to collect approximately 1 kg of rock. 45 3.4.4 10Be laboratory procedures and AMS measurements The Lamont-Doherty Earth Observatory Cosmogenic Nuclide Laboratory processed samples collected in 2014, and we analyzed the remaining samples in the Tulane University Cosmogenic Nuclide Laboratory. All samples were crushed, milled, and sieved to 250-750 μm. Physical and chemical isolation of quartz was completed following the procedures of Nichols and Goehring (2019). We isolated Be using standard chemical isolation procedures, including anion and cation exchange columns (Ditchburn and Whitehead, 1994; Schaefer et al., 2009). We included a process blank with every batch of ~eight samples (Table C.2). We sent sample aliquots of extracted Be to either the Purdue Rare Isotope Measurement (PRIME) Laboratory or the Lawrence-Livermore National Laboratory (LLNL_CAMS) for AMS measurements, which were normalized to the standard KNSTD dilution series (Nishiizumi et al., 2007). We calculated the exposure ages for all samples using version 3 of the online exposure age calculator formerly known as CRONUS-Earth, hosted by the University of Washington (https://hess.ess.washington.edu/). We used the default 10Be reference production rates from the “primary” calibration dataset (Borchers et al., 2016) and report individual sample ages using the Lifton-Sato-Dunai (LSDn) scaling scheme and 1-sigma analytical errors (Table 3.1). No corrections for burial by snow or surface erosion are applied to the moraines as snow depth and its variation and rates of surface erosion are poorly constrained. We do, however, provide estimates of how exposure ages may be influenced by snow cover (Table C.3). Moraine ages are reported as the median exposure age ± interquartile range to avoid the issue of using statistics that assume an underlying distribution of the ages of the moraine boulders, a key requirement of parametric approaches to characterize central tendency and dispersion (Menounos et al., 2017; Darvill et al., 2022). 46 Table 3.1: 10Be sample information for all boulders sampled in this study. 3.4.5 ELA reconstructions Variations in the equilibrium line altitude of a glacier relate to long term changes in climate. Such variations have been used to estimate changes in either temperature or precipitation (Dahl and Nesje, 1992; Moore et al., 2022; Oien et al., 2022). Commonly used 47 methods to reconstruct past ELAs include the maximum elevation of lateral moraines, toe-to-headwall altitude ratio, and accumulation area ratio, among others. Each method offers advantages and limitations in reconstructing past ELAs (Benn et al., 2005; Nesje, 1992; Porter, 2001; Osmaston, 2005). We use the MELM, THAR, and AAR methods of ELA reconstruction to estimate glacier ELAs between the Little Ice Age (ca. 1300-1850 CE) and modern time (2000-2021 CE). To record the MELM for each glacier, we used high resolution satellite imagery and elevation data from ASTER GDEM version 3 (NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team, 2019) to identify the highest elevation of preserved lateral moraines. The THAR method assumes a glacier’s ELA is positioned at a fixed ratio between the maximum and minimum elevation of the glacier, shown in Eq. (1): = + ( ) (1) Work by Meirding (1982) and Murray & Locke (1989) found that ratios of 0.35 to 0.4 yielded satisfactory estimates of alpine glacier ELAs. Here, we use the mean ELA from a THAR of 0.35 and 0.4. The accumulation area ratio assumes a fixed ratio of the accumulation area to the total area of a glacier in equilibrium (Braithwaite and Raper, 2009; Meier and Post, 1962). Here, we assume the AAR for glaciers in this region to be 0.6, which is generally considered to be the ratio of steady state cirque and valley glaciers in NW North America (Porter, 1975). We generated LIA and modern glacier hypsometries by clipping the ASTER DEM to the digitized glacier extents. In this case, the modern glacier extents are from the latest satellite imagery used for each glacier (imagery from 2017-2021 CE). We acknowledge that the modern 48 DEM does not account for the paleo surface of the glacier during the LIA and may negatively bias the paleo-ELA (Porter, 2001). For each ELA reconstruction method, we inferred the change in average temperature (dT) from the Little Ice Age to present as a function of changing ELA by assuming an environmental lapse rate of -6.5 °C km-1. The ELA of a glacier is also influenced by changes in precipitation. Ohmura et al. (2018; 1992) empirically derive an equation (Eq. 2) to estimate the annual precipitation, P, in millimeters water equivalent (mm w.e.) at the ELA of a glacier, given a mean summer (JJA) temperature T: = + 2 + (2) , where, a = 966, b = 230, and c = 5.87. We estimated changes in precipitation at the ELA of each study glacier by assuming a modern (1986-2015 CE mean) JJA temperature (T) at the modern ELA from the fifth generation European Centre for Medium-Range Weather Forecasts (ECMWF) global climate atmospheric reanalysis (ERA5). We use our dT estimate from our ELA reconstructions to yield Eq. 3: = + ( − ) + ( − ) 2 (3) We selected ERA5 2 m surface temperatures (Hersbach et al., 2020) from the grid cell nearest to the study glacier and used the same -6.5 °C km-1 lapse rate to approximate T at the modern ELA. 49 3.4.6 Glacier modeling 3.4.6.1 Open Global Glacier Model Our final method of ELA reconstruction uses the Open Global Glacier Model (OGGM; Maussion et al., 2019) which is a modular, open-source model framework with the capacity to model glacier evolution for all glaciers on Earth. The glacier model within OGGM is a depth-integrated flowline model that solves the continuity equation for ice using the shallow ice approximation (Cuffey and Paterson, 2010). Multiple flowlines for each glacier are calculated using a DEM clipped around the glacier polygon using the routing algorithm of Kienholz et al. (2014). The default mass-balance model used in OGGM begins with gridded monthly climate data, here the Climatic Research Unit gridded Time Series (CRU TS) version 4.04 (Harris et al., 2020). The climate data feeds a temperature index model described in Marzeion et al. (2012), incorporating a temperature sensitivity parameter that is calibrated using nearby glaciers with observations of specific mass balance (Zemp et al., 2021). Ice thickness is estimated by assuming a given glacier bed shape (parabolic, rectangular, or mixed) and applying a mass-conservation approach that employs the shallow-ice approximation. OGGM assumes that the “modern” glacier outline, sourced from the Randolph Glacier Inventory (RGI), is from the same date as the DEM. Users are also able to supply their own glacier outlines. More information on OGGM can be found on OGGM.org, or in publications on the model (Maussion et al., 2019; Eis et al., 2021). 3.4.6.2 Equilibrium run In our first experiment using OGGM, we started with the RGI polygons for the six of our study glaciers targeted for surface exposure dating (Anderson, Mordor, Butterfly, North Moraine Hill, Keele Peak, and Arrowhead glaciers). We then ran a 1000-year simulation under a constant climate, iteratively adjusting a temperature bias relative to the average CRU TS climate centered 50 around 2000 CE (close to the RGI polygon date of most glaciers in the region) until the modeled glacier reached equilibrium at or very near the glacier length indicated by the moraine record. From these equilibrium run experiments, we produce three different estimates of ELA and temperature change. First, the temperature lowering required to expand a glacier to its LIA length was interpreted as the approximate temperature change from the LIA to 2000 CE. Second, we then extracted the hypsometry of the modeled glacier at t=0 (modern extent) and t=1000 (LIA extent) and estimated the modeled ELA using the same AAR method as described in section 3.5, again assuming an AAR of 0.6. We can again apply the -6.5 °C km-1 lapse rate to estimate the apparent temperature change from modelled glacier extents between the two time periods. Third, for the modern glacier extent, we extracted the elevation at which the modeled surface mass balance of each glacier is equal to zero without any temperature bias. This represents the modern climatic ELA and is not based on glacier morphology. 3.4.6.3 Transient run In our next experiment with OGGM, we simulate changes in glacier volume in the Mackenzie and Selwyn mountains using our glacier chronology to tune the climate model input. We used OGGM to simulate the response of our five glaciers driven by monthly temperature and precipitation variability from four Coupled Model Intercomparison Project Phase 5 (CMIP5) GCM runs (CCSM4, MIROC-ESM, MPI-ESM-P, and MRI-ESM2; Taylor et al., 2012). All GCMs incorporate volcanic, total solar irradiance, summer insolation in both hemispheres, aerosol and greenhouse gas emission, and land use change forcings over the period 850-2005 CE (Landrum et al., 2013; Sueyoshi et al., 2013; Yukimoto et al., 2019). We omitted the glacier on Keele Peak, as its RGI outline includes several cirque glaciers separated from the main glacier, which causes OGGM to produce a problematic flowline that 51 crosses several flow divides. We set the mass balance gradient for each glacier to 5.2 mm w.e. m-1 based on the mass balance gradient for Bologna Glacier in Nahanni National Park Reserve for the 2014-2015 CE balance year (Ednie and Demuth, 2019). For each GCM, we ran 300-500 simulations incrementally perturbing the temperature bias (Tbias) and unitless precipitation factor (Pbias) to determine which combination of temperature and precipitation bias produces a modeled glacier length time series that best fits our glacier chronology. Tbias values ranged from -5 to +2 °C and Pbias between 1.0 and 4.0. Initial testing prior to running the larger simulations showed that Tbias and Pbias values beyond the above range produced glacier extents that far exceeded the late Holocene maximum extent of the glacier or made them disappear entirely. When the glacier flowline exceeded 80 grid points beyond the modern glacier extent, the simulation was discarded. For each simulation, we calculated the summed root mean squared error (RMSE) of modeled glacier length versus the moraine and remotely sensed glacier length at multiple timesteps. The combination of Tbias and Pbias that produced the lowest RMSE was selected as the “optimized” set of parameters for each glacier and GCM. The exact values of Tbias and Pbias are not meant to convey specific information about past climate. These values allow for regional tuning of the OGGM model to better fit the reconstructed and observed glacier response. Finally, we averaged the set of Tbias and Pbias from each glacier that produced the lowest RMSE for each GCM and applied those corrections before running simulations of the past millennium for all (1,235) glaciers in the eastern YT/NWT, forced by each “calibrated” GCM. The past millennium climate is of interest as it covers the onset and termination of Little Ice Age cooling. We start all past millennium runs at 1000 CE. We then compared the modeled glacier volume change over the past millennium to our chronology as well as what is already known 52 about late Holocene glacier change in this region to evaluate if the modeling results were reasonable. 3.4.6.4 Future glacier simulations To predict the fate of glaciers in this region, we use OGGM to project 21st-century glacier change for all 1235 glaciers in the eastern Yukon and Northwest Territories, forced by four different CCSM4 projection runs under different representative concentration pathways (RCPs). We use the default model parameters of OGGM v1.5.3 and rely on OGGM’s pre-processed glacier directories, which already contain glacier geometry and climate data. The historical climate data is CRU TS version 4.04 (Harris et al., 2020). We then download the CMIP5 (CCSM4) climate model output from four different RCP’s and run OGGM’s bias correction against the CRU calibration data, which in turn calculates anomalies from the CRU reference climatology (1961-1990 CE). Finally, we run OGGM for all 1235 glaciers forced by the calibrated climate scenarios from 2020 to 2100 CE and analyze the projected change in glacier area and volume. 3.5 Results 3.5.1 Glacier chronology Glaciers in the Mackenzie and Selwyn mountains deposited moraines fronting cirque and valley glaciers 0.7 to 2 km beyond their ca. 2020 CE extents. These moraines are typically devoid of vegetation other than widespread lichen cover. The moraines we sampled are commonly boulder-rich, with pebble-cobble matrices. Many alpine cirques preserve two nested moraines within tens of meters of each other. We observed nested moraine crests at Keele Peak, Arrowhead, North Moraine Hill, and Mordor 53 glaciers. There is also a partially-nested crest preserved at Anderson Glacier. We did not sample both crests at most locations since our focus was to date the outermost moraines. Erratic boulders 10-40 m beyond cirque moraines at Anderson and Mordor glaciers date to 10.9-11.1 ka (Table 3.1). An erratic sampled ~250 m beyond the late Holocene moraine fronting Keele Peak glacier dates to 11.6 ± 0.3 ka. Erratic boulders directly overlaid bedrock and had abundant lichen cover. We did not observe any obvious signs of boulder surface erosion, such as grüssification, solution pitting, or enhanced relief of resistant minerals. In the Nahanni National Park region, the median 10Be age on moraine boulders is 610 ± 850 a (ca. 1405 CE, n = 19). Adjacent to Keele Peak, the median moraine exposure age is 370 ± 110 a (ca. 1650 CE, n = 8). Together, the sampled moraines in this study date to 460 ± 415 a (ca. 1560 CE). We sampled both the inner and outer crest of the moraine couplet at Arrowhead and Mordor glaciers. At Anderson Glacier, the outer moraine dates to 390 ± 50 a (1620 CE, n = 3) and the inner moraine to 150 ± 20 a (1860 CE, n = 2). At Mordor Glacier, the outer moraine dates to 1260 ± 295 a (760 CE, n = 3) and the inner moraine dates to 390 ± 20 a (1630 CE, n = 2). There is notable scatter in the exposure ages on many of the sampled moraines (Table 3.1, Fig. 3.3). At Nahanni 01, Butterfly, Anderson, Mordor, and North Moraine Hill glaciers, there is at least one sample from each moraine that returned ages older than 1 ka. This scatter gives individual moraine ages large errors, however when we analyze all moraine boulder ages together, there is a distinct peak in exposure ages between ~800 to 100 a exposure (ca. 1200 to 1900 CE), with the greatest peak around 480 to 280 a (1540-1740 CE, Fig. 3.3). 54 Figure 3.3: Box and whisker plots of 10Be surface exposure ages for each glacier, showing the interquartile range and median age of each moraine surface and the normalized probability density function (red line) for all 10Be samples and kernel density plot (grey lines) for each individual 10Be sample. 3.5.2 Climate reconstructions since the late Holocene ELA reconstruction using the different methods described above yield a range of estimated changes in ELA between the LIA and modern time (Fig. 3.4). We use ELAs from the AAR method using mapped former and modern glacier extents as the “standard” ELA against which we compare our other ELA estimates. Any ELA reconstruction method could serve as the “standard”; the AAR method was selected due to its common usage in glacier reconstructions (Benn et al., 2005; Dahl and Nesje, 1992; Oien et al., 2022). When comparing ELA change within a single method, “dELA” is the change in reconstructed ELA between the LIA and modern time using the method in question. As discussed more below, we assume that 55 precipitation remains constant between the LIA and modern time for ELA reconstructions using the MELM, THAR, and AAR methods. Figure 3.4: Changes in ELA and estimated temperature change between the Little Ice Age maximum to modern (ca. 2015) for six glaciers in this study. Each bar represents a different ELA reconstruction method as described in text. OGGM TLower is the temperature lowering from ca. 2000 CE climatology required to allow the modeled glacier to reach their late Holocene maximum extent. OGGM Smb is the change in ELA where the modeled surface mass balance on the glacier equals zero between the late Holocene maximum and ca. 2000 CE. OGGM AAR is the difference in AAR-derived ELA from the modeled glacier extent at the late Holocene maximum and ca. 2000 CE. Green squares with capped error bars are the mean and 1-sigma standard deviation for all ELA reconstruction methods for each glacier. The modern ELA derived from the AAR method is +12 m to +171 m (average 107 m) higher than the LIA ELA using the maximum elevation of lateral moraines method, corresponding to a +0.1 to +1.1 °C (average 0.9 °C) increase in temperature (Fig. 3.4). Using the THAR method, the dELAs range from +47 m to +240 m (average 138 m), corresponding to a dT of +0.3 to +1.6 °C (average 0.9 °C) since the LIA. 56 ELAs reconstructed from LIA and modern glacier extent mapping, assuming an AAR of 0.6, indicate a rise in ELA since the LIA of +60 to +182 m, corresponding to a +0.4 to +1.2 °C (average 0.8 °C) increase in annual average temperature (Fig. 3.4). Using OGGM, we include three estimates of ELA change. Non-transient simulations on glaciers in the Nahanni National Park region using OGGM require +2.3 °C of warming, relative to the 30-yr average climate centered around 2000 CE, to retreat from their LIA extents to modern positions. Keele Peak and Arrowhead glaciers require nearly +3.2 °C average warming since the LIA relative to their modern temperature (Fig. 3.4). This warming is equivalent to a dELA since the LIA of +354 m in Nahanni National Park and +492 m in the Keele Peak area. Applying the AAR method, but with OGGM-derived glacier hypsometries at the LIA and modern time, indicates much less warming since the LIA, with rises in ELAs between +7 m and +76 m, corresponding to a rise in temperature of <0.1 to 0.5 °C. We interpret this minimal change in ELA to be the result of glacier surface thickening in the OGGM model when the glacier expands to LIA extents, which reduces the apparent ELA change as the lower portion of the modeled glacier surface thickens (Fig. C.5 & C.6). The third variation of ELA reconstruction using OGGM estimates the modern ELA not from modeled glacier hypsometry, but rather the elevation at which the modeled surface mass balance on the glacier is equal to zero. In a warming climate, this estimate of glacier ELA is expected to be higher than the AAR-derived ELA, as a glacier undergoing rapid retreat has a morphometry that lags behind the climate signal. Changes in ELA using the modern mass balance-derived ELA and the AAR-derived LIA ELA range from +277 m to +535 m. Estimated temperature change indicates a rise in temperature since the LIA of +1.6-3.5 °C. 57 Using the equation of Ohmura et al. (2018) and temperature change estimates from our AAR-derived ELAs, we estimate that compared to modern values, there was -117 to -339 mm w.e. yr-1, or 5-15% (average 10%), less precipitation at the ELA of our study glaciers during the LIA. 3.5.3 Past millennium glacier change Estimates of glacier evolution in the YT and NWT over the past millennium vary among the four GCMs (Fig. 3.5). The MPI simulation shows steady glacier volume until 1600 CE, while MRI, MIROC, and CCSM4 indicate a reduction in glacier volume until ca. 1250 CE, afterwards CCSM4 and MRI (and to a lesser degree MPI) show an increase in glacier volume until ca. 1400 CE before a period of stable ice volume until ca. 1600 CE. MRI, MPI and CCSM4 all indicate glacier expansion ca. 1600 CE, with MPI reaching a maximum ice volume of 38.1 km3 at 1765 CE and CCSM4 producing a maximum ice volume of 34.7 km3 at 1855 CE (Fig. 3.5). MRI appears to largely miss 20th century glacier retreat and continues to show glacier expansion until 1980 CE, followed by volume loss. Glacier volume simulated by MIROC decreases through the past millennium, in contrast to the other GCM simulations. Projections of future glacier loss (below) using CCSM4 climate simulations begin with an initial regional ice volume of 18.1 km 3 in 2019 CE. Compared to the maximum modeled ice volume in the CCSM4 past millennium simulations, this represents a 48% loss in ice volume since ca. 1850 CE. 58 Figure 3.5: Modeled ice volume change for all glaciers in the eastern YT and NWT produced by OGGM using four different GCMs. Dashed lines from 1000 CE to 1250 CE are used to indicate spin up duration of the model. Dots and vertical lines respectively denote average and standard deviation (1-sigma) of normalized mean glacier length binned by decade. 3.5.4 21st Century glacier projections Under all CCSM4 21st century emissions scenarios, glacier volume in the eastern YT and NWT significantly declines throughout this century (Fig. 3.6). Glacier volume is projected to decrease by 85% under RCP2.6 and 97% under RCP8.5, compared to 2019 CE values. The greatest rate of ice loss is projected to be between present day and ca. 2040 CE, then the rate of volume decline slowly decreases through to the end of the century. 59 Figure 3.6: Fractional glacier volume change until 2100 CE under various representative concentration pathways (RCPs) for all glaciers in the eastern YT and NWT. 3.6 Discussion 3.6.1 Holocene glacier fluctuations Early Holocene erratic boulders just beyond moraines dating to the last millennium, as well as a lack of moraines down valley of the latest Holocene moraines, implies that since ca. 11 ka glaciers in this region were no more extensive than during the latest Holocene. These results accord with records from southern Alaska and western Canada (Menounos et al., 2009; Mood and Smith, 2015; Barclay et al., 2009) that show most alpine glaciers within these regions reached their greatest Holocene positions during the last several hundred years. We interpret the erratic boulders of latest Pleistocene age to record local deglaciation associated with the termination of the Younger Dryas cold interval (Menounos et al., 2017; Seguinot et al., 2016; Braumann et al., 2022). Similar erratic boulders that lie beyond late Holocene cirque moraines 60 were dated by Menounos et al. (2017) and were also interpreted to record local deglaciation. The erratic boulders sampled in the present study were not part of a moraine, so their ages are interpreted to reflect deglaciation at those sites; the absence of an associated moraine precludes us from drawing conclusions about the size of the up valley glaciers. The most parsimonious explanation for coeval ages of erratic boulders and end moraines is the complex decay of the Cordilleran Ice Sheet; some cirques were still covered by the ice sheet while others were ice free prior to the Younger Dryas and so were able to form an end moraine (Menounos et al., 2017). Our moraine chronology generally accords with the limited previous work in this region. Moraine ages from this study suggest glaciers reached their LIA maximum closer to 1560 CE, with a possible readvance or standstill in the mid-1800’s. Tomkins et al. (2008) used varve and tree ring records near Tungsten, YT to infer periods of glacier growth around the late 1300s to 1450 CE, 1600 to 1670 CE, 1730 to 1778 CE, and an apparent Little Ice Age maximum 1778-1892 CE. Dyke (1990) completed an extensive lichenometric survey of rock glaciers and late Holocene moraines directly west and south of Tungsten, dating most late Holocene moraines to within the past 400 years. Our moraine chronology is in general agreement with the lichenometric ages of Dyke (1990) and suggests an earlier Little Ice Age maximum than interpreted by Tomkins et al. (2008). The significant scatter in our 10Be moraine dataset complicates our interpretations of decadal-to-century scale glacier fluctuations, however. Several scenarios could yield moraine exposure ages that are either older or younger than the true depositional age of the moraine. Inherited nuclides from episodes of previous exposure would result in exposure ages older than the true depositional age. One source of inherited nuclides could be from rockfall followed by supraglacial transport before deposition on the moraine. It is also possible that there was insufficient resetting of the 10Be inventory in the local 61 bedrock during the Last Glacial Maximum (LGM) as these sites sit at the periphery of the LGM extent of the Cordilleran Ice Sheet. A third possibility is that the inclusion of old outliers reflects the incorporation of previously exposed boulders within the glacier forefield. A review of Holocene glacier fluctuations in western Canada revealed a progressive expansion of ice that culminated with climatic advances during the Little Ice Age (Menounos et al., 2009). Given what is known about Holocene glacier activity, the most likely explanation for our pre Little Ice Age boulder ages is that these boulders contain inherited nuclides from previous moraine building events and were subsequently reincorporated into the late Holocene moraines during the advances of the Little Ice Age. A final possibility to explain the scatter in our moraine ages is that many boulder ages are too young. Mass shielding by previous burial within a moraine followed by exhumation of a sampled boulder, or from snow cover, would reduce the nuclide production rate and result in erroneously young exposure ages. Exhumation and post-depositional movement would be more likely if our moraines were originally ice cored (Crump et al., 2017). Snow cover results in younger apparent ages on moraine boulders, however unrealistic quantities of snow cover are required to meaningfully impact the exposure age of our moraines. One meter of 0.25 g cm-3 snow on the surface of our boulders for four months of the year would decrease the calculated age by 15-27% (Table C.3). This decrease in age does not significantly impact our interpretations, as the moraines would still predominantly date to the Little Ice Age. The timing of glacier fluctuations in the eastern Yukon and Northwest Territories agrees with records of late Holocene glacier advance in Europe (Braumann et al., 2020, 2021; Ivy-Ochs et al., 2009). Though Europe has different climate forcings than western North America, the similar timing of late Holocene glacier response suggests that lower temperatures associated with 62 decreasing summer insolation in the Northern Hemisphere played an important role in the timing of glacier advance in the late Holocene in both regions. 3.6.2 ELA and climate reconstruction In this study, we reconstructed and estimated past and present glacier ELAs through several methods, inline with recommendations by Benn et al. (2005) that multiple ELA reconstruction methods be used to provide a more robust estimation of past ELAs and uncertainty with each reconstruction method. An important limitation to the AAR and THAR methods is that they do not account for modern glaciers being out of equilibrium with modern climate. If the modern ELA is not accurately known and the glacier is retreating or advancing in response to climate perturbations, then comparisons in ELA change between modern and other time periods will under- or over-estimate ELA departures (Porter, 2001). Additionally, the assumption that a glacier’s ELA only fluctuates due to changes in temperature is an oversimplification (Ohmura et al., 1992). Increased (decreased) precipitation will lead to a higher (lower) mass balance and may obscure the impact of temperature change on glacier response (i.e. Shea et al., 2004). Anderson et al. (2011) presents lacustrine δ18O records from the central Yukon that suggest a wet, early Little Ice Age, then dry conditions until modern day, in response to the changing position and strength of the Aleutian Low. If glaciers in the Mackenzie and Selwyn Mountains received greater snowfall during the LIA, then less cooling would be needed to grow glaciers to their LIA extents. Tomkins et al. (2008) developed a July mean temperature reconstruction from tree rings and varved lake sediments close to Tungsten, near the northern end of Nahanni National Park Reserve. Their amalgamated temperature reconstruction demonstrates the differing signals of varved lacustrine sediment and tree ring records but does 63 suggest cooler temperatures in the early 1800’s, a warm interval at the end of the 1800’s to early 1900’s, followed by cooling until at least the 1940’s before warmer than average July temperatures until modern time. Our non-transient experiment using OGGM provides another estimate for temperature change since the LIA, though it still ignores the effect of precipitation variability. By determining the temperature lowering from the present climate needed to grow a modeled glacier to LIA extents, we remove the likely erroneous estimation of the modern glacier ELA based on current glacier hypsometry and more directly compare modern temperatures with the inferred temperature during the LIA maximum, when the glacier was in equilibrium with climate. Both the non-transient (“OGGM TLower” in Fig. 3.4) and surface mass balance (“OGGM Smb” in Fig. 3.4) incorporate modern climatology and as a result indicate generally greater temperature change since the LIA compared to glacier geometry-based reconstruction methods. A bedrock borehole temperature reconstruction (62.47° N, 129.22° W) between Nahanni National Park and Keele Peak indicates around +3 °C of surface warming since 1500 CE (Huang et al., 2000), consistent with our temperature change estimates comparing past ELAs to modern climatology. A similar study design as presented in this manuscript would be improved by selecting a site with a multi-year in situ mass balance record to compare the modelled modern ELA estimate with the ELA derived from in situ measurements. OGGM is built to perform best at regional to global scales and may produce problematic results at the scale of individual glaciers (Maussion et al., 2019). Differences between the year of DEM acquisition and RGI glacier extent, erroneous glacier margins, and lack of nearby mass balance calibration information can all have significant impacts on the evolution of individual modeled glaciers. To help give confidence that the modeling results from OGGM were 64 producing reasonable glacier evolution, we ran a simple flowline glacier model modified from Jarosch et al. (2013), which was able to grow glaciers to similar extents as OGGM (Fig. C.2). The similar glacier evolution between the two models indicates that modeled glacier response is the result of climate inputs, rather than unique properties of each model. As mentioned above, regular mass balance data from in situ mass balance measurements or remote sensing on glaciers in remote areas will help improve the performance and validation of global glacier models like OGGM (Eis et al., 2021). A similar study design as is presented in this paper may be successfully implemented in areas with robust glacier chronologies from the late Holocene to present from many more glaciers than are included in our study. Well-constrained glacier chronologies would serve to extend the calibration or validation dataset for large scale glacier modeling efforts (i.e. Rounce et al., 2023). 3.6.3 GCM evaluation Of the four different CMIP5 GCM simulations tested, glacier model runs forced by CCSM4 and MPI yield glacier fluctuations that best match our general understanding of latest Holocene glacier expansion and glacier retreat over the past millennium (Menounos et al., 2009; Luckman, 2000; Fig. 3.5). We consider the results from MRI to be unreasonable due to the continued ice expansion through most of the 20th century, and similarly discount the results from MIROC due to the modeled steady glacier volume decline over the entire past millennium. Our 10Be chronology suggests glacier advance and moraine formation earlier than what the modeling results show. At Arrowhead Glacier, the outer and inner moraine 10Be ages (1620 and 1860 CE, respectively) are comparable with the modeled glacier evolution under the CCSM4 climate, however. MRI suggests a period of glacier retreat shortly before 1600 CE, which is consistent with our moraine chronology, however MRI, CCSM4, and MPI all suggest further ice 65 expansion which would have overridden previously deposited moraines. If the exposure age of a moraine is interpreted to more closely record the onset of glacier retreat, rather than advance, then our moraine chronology further indicates that glaciers reached their LIA maximum extents prior to when OGGM suggests. The four GCMs used in our study simulate varied temperature and precipitation time series over the past millennium, which results in differing modeled glacier responses (Fig. C.8-11). Modeled glaciers forced by CCSM4 and MPI reach late Holocene maxima between 1765 and 1860 CE, coincident with other late Holocene glacier records (Menounos et al., 2009; Barclay et al., 2009; Mood and Smith, 2015). Our moraine and remote sensing record allowed for four GCM’s to be calibrated for a small selection of glaciers in the region prior to being run for all 1235 glaciers. Without a well-dated moraine chronology, we would be unable to assess how to model performs beyond the remote sensing record. Further research is needed to evaluate why the existing GCM simulations fail to grow glaciers at the same time as our moraine chronology suggests in northwestern Canada. The moraine record offers an important method of validating glacier models beyond the remote sensing record, however moraine chronologies must be tightly constrained in order to confidently evaluate model results. Additional cosmogenic surface exposure dating in this region, especially in areas where there is an unambiguous lack of post-depositional movement may help to produce moraine chronologies with less scatter. Measuring multiple nuclides on moraine boulders (such as using paired 14C-10Be) would allow potential inheritance to be investigated (i.e. Goehring et al., 2022). Finally, as mentioned above, consistent mass balance records from glaciers in this region would help to better constrain the influence of local climate 66 on glacier response in the Mackenzie and Selwyn Mountains (Pelto et al., 2019; Ednie and Demuth, 2019). 3.6.4 Future response of glaciers to climate change The Mackenzie and Selwyn mountains are almost certain to experience profound glacier mass loss throughout the 21st century. The estimated magnitude of ice volume decline agrees with modeling results by Clarke et al. (2015) who estimate a 70-95% reduction in glacier volume in the Canadian Rocky Mountains by 2100 CE. Additionally, recent work by Rounce et al. (2023) estimates 93-100% deglaciation in the Mackenzie and Selwyn Mountains by 2100 CE, depending on the magnitude of global temperature change. Under SSP3.7 and SSP5.85, this region is predicted to be fully deglaciated by 2080 CE (Rounce et al., 2023). By 2019 CE, approximately half of the ice volume was lost in the Mackenzie and Selwyn Mountains in the CCSM4 run compared to the glacier maximum in 1860 CE (Fig. 3.5). The loss of glaciers in this region will cause greater fluctuations in streamflow and temperature that may have negative impacts on thermally stressed species, including fish that are important food sources for local communities (Babaluk et al., 2015; Clason et al., 2023; Moore et al., 2009). 3.7 Conclusions Based on geomorphic mapping, surface exposure ages, and numerical modeling, the following conclusions can be drawn from our study. (1) The probability distribution of 10Be ages suggests that most glaciers in eastern YT and NWT reached their greatest Holocene extents during the latter half of the Little Ice Age [1600-1850 CE]; (2) The uncertainty ascribed to some moraines is high, given the presence of some boulders that yielded 10Be ages that predate the Little Ice Age, and future work utilizing multi-nuclide approaches would allow this scatter to be 67 further investigated; (3) We find no evidence of glaciers extending beyond LIA limits since at least 10.9-11.6 ka, in accord with most other Holocene glacier records in the Northern Hemisphere; (4) Our ELA reconstructions suggest warming of 0.2-2.3 °C since the LIA, with morphology-based ELA reconstructions likely underestimating the modern ELA of glaciers undergoing retreat; and (5) Projections of future glacier change estimate a further 85% - 97% loss of glacier volume in the Mackenzie and Selwyn mountains by 2100 CE, in agreement with recent global modeling efforts. Glacier chronologies from late Holocene glacier fluctuations can provide important sources of validation of GCM simulations beyond the instrumental record, especially given the variety between individual GCM simulations of past climate. Nearby in situ mass balance records and well-constrained late Holocene glacier chronologies are needed to help validate past millennium GCM simulations and highlight important feedbacks between the arctic and the global climate system. Modern tropospheric warming will continue to dramatically reduce glacier volume in this region, with significant impacts to the local ecosystem that relies on glacier-fed rivers and streams through the summer months. 68 4. Terrestrial Cosmogenic Nuclide depth profiles used to infer changes in Holocene glacier cover, Vintage Peak, Southern Coast Mountains, British Columbia 4.1 Abstract In North America, most glaciers reached their maximum Holocene downvalley positions during the Little Ice Age (1300-1850 CE), and in most cases, this expansion destroyed earlier evidence of glacier activity. Substantial retreat in the 20th and early 21st centuries exposed bedrock that fronts many glaciers that may record early-to-mid Holocene exposure and later burial by ice which can be elucidated using multiple-nuclide cosmogenic surface exposure dating. Furthermore, cores of bedrock allow the measurement of cosmogenic nuclide depth profiles to observe the slope of the concentration curves and thus better constrain likely exposure and burial histories. To investigate this technique, we collected four bedrock surface samples for 10 Be and 14C surface exposure dating, as well as bedrock cores up to 55 cm deep on Vintage Peak, in the southern Coast Mountains of British Columbia, Canada. We apply a Monte Carlo approach to generate combinations of exposure, burial and subglacial erosion that can be explained by our data. We find that Vintage Peak became deglaciated between 14.5 and 11.6 ka, though alpine sites within late Holocene maximum extents were not exposed until 10-11 ka. Glaciers on Vintage Peak advanced to near late Holocene maximum positions around 5-6 ka. Nine 10Be ages on late Holocene moraines reveal that glaciers reached their greatest Holocene extents ca. 1300 CE. Our results agree with other regional glacier records and demonstrate the utility of surface exposure dating applied to deglaciated bedrock as a technique to help construct a record of Holocene glacier activity where organic material associated with glacier expansions may be absent or poorly-preserved. Further work to increase exposure/burial history modeling complexity may help to better constrain complex exposure histories in glaciated alpine areas. 69 4.2 Introduction The most unequivocal evidence to gauge the magnitude of a glacier advance is to measure dimensional changes of a glacier such as its length, as is recorded in end moraines. In many cases, these end moraines are well preserved, can be directly dated, and thus provide important proxies for glacier expansion. Unfortunately, the terrestrial record of Holocene alpine glacier fluctuations in North America is largely incomplete due to most glaciers reaching their greatest Holocene extents during the latest Holocene, overriding or destroying evidence of previous glacier positions (Barclay et al., 2009; Menounos et al., 2009; Solomina et al., 2015). Alternative methods are needed to determine the duration of ice cover at alpine sites in order to place current rates of glacier retreat into a long-term perspective. Cosmogenic surface exposure dating offers a useful method to document past glacier activity (Balco, 2010). Cosmogenic surface exposure dating uses the accumulation, and in some instances decay, of rare isotopes in rock at the surface of the Earth to calculate the duration the surface has been exposed and/or buried (Balco, 2010; Gosse and Phillips, 2001). For late Pleistocene and Holocene studies, the combination of 14C (t1/2 5,730 yr) and 10Be (t1/2 1.38 Myr), provide several advantages over the measurement of a single nuclide alone. Due to the relatively short half-life of 14C, times of previous exposure, such as during Pleistocene interglacials are quickly “forgotten” through radioactive decay. Conversely, 10Be decays slowly enough that it is effectively stable over Holocene timescales. In situations of negligible erosion, 10Be acts to record total exposure duration at a site, while 14C can elucidate periods of burial. At locations where a well-constrained glacier chronology has been established, the cosmogenic nuclide inventory within previously glaciated bedrock surfaces has been used to constrain subglacial erosion rates (Balter-Kennedy et al., 2021; Goehring et al., 2013). Hippe (2017) provides a clear 70 overview of the concepts behind using paired 14C and 10Be to investigate Holocene surficial processes. In formerly glaciated alpine sites in western Canada, bedrock surfaces that lie beyond limits of late Holocene glacier expansion would have been exposed to cosmogenic nuclide flux since retreat of the Cordilleran Ice Sheet following the termination of the Last Glacial Maximum (Menounos et al., 2017, 2009; Seguinot et al., 2016). The CIS, large decaying masses of the CIS, or even cirque glaciers may have experienced a re-advance during the Younger Dryas cold interval, afterwhich the ice sheet retreated rapidly (Menounos et al., 2017). Alpine glaciers in western Canada were likely smaller than their 2000-2010 CE extents for much of the early Holocene (Koch et al., 2014, 2007b; Menounos et al., 2004) and may have begun periodically advancing as early as 8 ka, with more evidence for glacier advances of generally increasing magnitude by ca. 5 ka. At most sites in western Canada, glaciers reached their greatest Holocene extent during the Little Ice Age, between 700-150 years before 2000 CE (a b2k) (Menounos et al., 2009; Mood and Smith, 2015a). What remains unclear is the size of alpine glaciers prior to Neoglacial expansion in the mid-Holocene and the amount of time glaciers were near their late Holocene maximum extents or retreated to within their modern (ca. 2017 CE) extents. We seek to investigate the record of late Pleistocene through Holocene glacier change at an alpine site in the southern Coast Mountains of British Columbia, Canada, through several different applications of cosmogenic nuclide surface exposure dating. We aim to constrain the duration and spatial variability of exposure, cumulative burial by ice, sediment, and/or snow, and subglacial erosion rates at several locations in an alpine cirque. This study highlights the utility of paired nuclide dating and bedrock cores to more closely constrain complex exposure histories, 71 the limitations of bedrock exposure dating, and provides new numerical ages for end moraines in the southern Coast Mountains. 4.3 Study Area Vintage Peak (1876 m above sea level - asl) is located at the head of Powell Lake in the southern Coast Mountains of British Columbia, Canada (50.25° N, 124.30° W, Fig. 4.1). The mountain is composed of a mid-Cretaceous granodiorite pluton that forms broad sloping surfaces and steep rock walls that extend down to the valley bottom (Bellefontaine et al., 1994). A small (<0.05 km2) north-facing glacierette referred to in this paper as “Lockie’s Glacier” is located at 1780 m asl in a cirque west of Vintage Peak, below a sub-peak called Lockie’s Table. Another north-facing glacierette, referred to as Vintage Glacier, resides on the north side of Vintage Peak. The cirque headwalls are variably steep, though large portions of the headwall give way to rolling, convex ridgelines (Fig. 4.1). 72 Figure 4.1: Bedrock and boulder sampling locations on Vintage Peak. 10Be surface exposure ages are listed in text boxes with analytical errors. Green star on inset marks the study location. Base hillshade is from a 2017 aerial LiDAR survey. Contour interval is 100 meters. The study site is 50 kilometers east of the Strait of Georgia, with the coastal climate producing high precipitation and moderate temperatures. Between 1998 and 2022, the mean annual air temperature and mean annual precipitation at Toba Camp (52 m asl), 38 km NE of Vintage Peak was 10.1 °C and 1401 mm, respectively (Env. and Climate Change Canada, 2022). Assuming a lapse rate of 6.5 °C km-1, the average annual air temperature at 1600 m on Vintage Peak is close to 0 °C. There are two Province of British Columbia snow survey sites with near annual observations since 1938 CE located just 2 km away from Lockie’s Glacier at 910 and 1040 m a.s.l. (Fig. 4.1), that together have a mean annual snow water equivalent (SWE) at the 73 end of March or early April of 847 ± 399 mm. Snow survey data used in this study can be found on the BC Government’s Snow Survey website. 4.4 Methods 4.4.1 Field methods To constrain the age of the inferred late Holocene moraines in the study area, we sampled nine large boulders on or near the crest of moraines or those resting directly on exposed bedrock. Six boulders were from the Lockie’s Glacier moraine and three were from the Vintage Peak moraine (Fig. 4.1). In addition, we sampled one erratic (17-VP-07), approximately 60 m beyond the Lockie’s Glacier end moraine. At each site, we recorded boulder dimensions, GPS position and elevation, local horizon and boulder self-shielding using a Brunton compass and inclinometer. We used a gas-powered concrete saw and hammer and chisel to collect ~1 kg of rock from each boulder. Seasonal snow cover on moraine boulders partially shields the boulder surface from incoming cosmic rays, lowering the production rate of cosmogenic isotopes. The distribution of snow on and adjacent to moraines has not been well-documented, though common surface exposure sample collection procedures selectively target boulders on moraine crests, assuming these features are windswept and therefore are less impacted by snow shielding (Gosse and Phillips, 2001; Heyman et al., 2016). To test if moraine crests in this area have less snow than the surrounding landscape, the Hakai Airborne Coastal Observatory completed an aerial LiDAR and orthoimagery survey of Vintage Peak on September 15, 2017 and again on May 16, 2023. Snow is retained in May on Vintage Peak, while the September survey was largely snow-free. The elevation difference between the co-registered, 1-m spatial resolution digital elevation models 74 (DEMs) gives an approximation of the snow depth in May 2023 across the survey area, notably around the sample moraines and erratic. We collected four bedrock core samples within or just beyond the late Holocene extent of Lockie’s Glacier (Fig. 4.2). Core sample locations were chosen to capture a variety of potential bedrock exposure histories. At each location, we used a portable drill to collect a pair of 0.35-0.55 m deep, 41 mm outer diameter bedrock cores. Cores were not taken deeper than 0.55 m due to the slow penetration rate of the drill (supplied drill bits for the rental drill were for sedimentary rock opposed to igneous, intrusive rock). A pair of cores within 0.10-0.15 m of each other were taken at each site in order to have enough material at each depth for dating. Similar to boulder samples, we recorded sample coordinates, elevation, surface orientation, and horizon shielding at each bedrock site. Samples were transported in plastic tubes, with core orientation labeled and maintained throughout transport and laboratory sampling (Fig. 4.2). Below, we refer to the pairs of cores at each site as a single core. Figure 4.2: Bedrock core from site 17-VPC4. Photo scale converted to digital representation. One core (17-VPC4) was collected from a bedrock knob, approximately 150 m beyond the presumed late Holocene extent of Lockie’s Glacier. Cores 17-VPC3 and 17-VPC2 are within late Holocene glacier limits, with 17-VPC3 ~50 m inside the late Holocene moraine and 17-VPC2 ~40 m beyond the modern (2017 CE) ice limit (Fig. 4.1). 17-VPC1 is located just off a 75 broad ridgeline on the NE side of the Lockie’s Table cirque. We attempted to sample bedrock with minimal glacier erosion by sampling from areas with oxidized surfaces and preserved striations. We avoided surfaces on the lee of bedrock steps that may be the result of subglacial plucking (Graham et al., 2023; Rand and Goehring, 2019). 4.4.2 Laboratory processing All samples were processed at the Tulane University Cosmogenic Nuclide Laboratory. We followed the procedure of Nichols and Goehring (2019) to physically and chemically isolate quartz from each sample. Extraction of 10Be from the sampled boulders and the top 5 cm of each bedrock core followed standard chemical isolation procedures (Ditchburn and Whitehead, 1994). Each batch of approximately eight samples (Table D.1) included a sample blank. Each bedrock core was divided into five, 0.05 m-long sections. While both 10Be and 14C were measured at the top 0.05 m of each core, we only measured 14C along the remaining length of each core. We measured both 10Be and 14C on the one erratic boulder, 17-VP-07. Extraction of 14C followed the procedures of Goehring et al. (2019). Concentrations of 10Be were measured by accelerator mass spectrometry (AMS) at Purdue University's PRIME Laboratory, with all 10 Be measurements normalized to the 07KNSTD standard (Nishiizumi et al., 2007). We measured 14C concentrations via AMS at the Woods Hole Oceanographic Institution’s National Ocean Sciences Accelerator Mass Spectrometry (NOSAMS) facility. Sample blank corrections for the 14C samples followed the procedure outlined in Balco et al. (2022), wherein blank corrections are treated as a lognormal distribution. Exposure ages were calculated using version 3 of the online exposure age calculator formerly known as CRONUS-Earth (Balco et al., 2008; Borchers et al., 2016). Topographic and self-shielding corrections were applied to all moraine and bedrock sample surfaces. Corrections 76 for snow cover were applied to all samples separately and are discussed further below. Moraine ages are presented as a median ± interquartile range. This choice of summary statistics does not rely on the assumption of a normal distribution in our dataset. Median ages and individual sample ages are presented in years ago (a) before the year of sampling (2017 CE). 4.4.3 Bedrock exposure history modeling The measured nuclide concentration in each of the bedrock cores represents an integrated signal of exposure, burial by ice and snow, and erosion. While this integrated signal precludes us from being able to constrain the precise timing of Holocene glacier advance and retreat, we are able to apply reasonable assumptions to limit the most likely burial and exposure histories at each sample location. Monte Carlo simulations for each sample location varied the duration of exposure, burial, and erosion rate at each site. We used the exposure age of an erratic boulder outside of late Holocene glacier limits to constrain the maximum possible exposure duration for all sites, assuming the exposure age of the erratic represents the timing of late Pleistocene deglaciation on Vintage Peak. The minimum exposure duration at each sample location was constrained by the 14C exposure age of bedrock surface samples. Over Holocene timescales, 10Be is effectively stable compared to the much shorter half-life of 14C. We plotted the measured 14C-10Be ratios on a two-isotope diagram with burial isochrons to determine the apparent burial history of each sample. All sites are assumed to have been seasonally covered by snow throughout their exposure history. We used the average snow depth and density from nearby snow monitoring sites to apply a snow cover correction for 6 months of the year. The average (± standard deviation) snow 77 density is 0.38 ± 6.6 g cm-3 and average depth is 220 ± 102 cm. During periods of burial by ice, we assume the ice was sufficiently thick to cease nuclide production. We assume that whenever the glacier was overlying bedrock at each site, the glacier was thick enough to effectively attenuate all incoming cosmogenic particles. Subaerial erosion rates in alpine environments are typically low (Elkadi et al., 2022; Lehmann et al., 2020), rarely more than 0.001 cm yr-1. As such, we assume that during periods of exposure, there was negligible subaerial erosion of bedrock surfaces. Subglacial erosion rates have been shown to be highly variable across regions and within single glacier forefields (i.e. Koppes et al., 2015; Magrani et al., 2022; Rand and Goehring, 2019). We do not include results with subglacial erosion rates higher than 0.2 cm yr-1, as very few solutions can be found with greater erosion rates. For each site, our Monte Carlo simulation was run 10,000 times per discrete erosion rate ( ) of 0.0, 0.05, 0.1, or 0.2 cm yr-1, with each run picking a uniform, randomly distributed value of exposure duration (te) and burial duration (tb) in years. This results in 40,000 runs per site. The exposure duration is a value between the maximum possible exposure duration of 15 ka and the minimum exposure duration, which is the 14C age of the bedrock surface for each site. The burial duration is a value between zero and the difference between the maximum exposure duration and the minimum exposure duration at each site. The theoretical nuclide concentration at each depth in the model space was calculated using Eq. 1: = λ (1 − (− λ ) )( −λ ρϵ )( Λ ), (1) where Ni is the nuclide concentration (atoms g-1) of the isotope in question (14C or 10Be), Pi is the surface production rate (atoms g-1 yr-1), i is the nuclide-specific decay constant (s-1), is 78 the rock density (g cm-3), and Λ is the attenuation length (160 g cm-2). We assumed the bedrock has an average density of 2.65 g cm-3. We ran the model setup described above assuming no inherited 14C or 10Be. Since inherited 10Be is more likely because of its long half-life, we ran our model again for each site assuming a duration of possible inheritance as indicated by the relative 14C-10Be concentrations at Core 4. This produced another 40,000 total simulations per site. We tested the influence of a 0.001 cm yr-1 subaerial erosion rate on our simulations and found that for Cores 1-3 there was no significant difference in the modeled burial or exposure ages, and at Core 4, inclusion of subaerial erosion lead to no simulations producing exposure and burial histories that matched our observations. Thus, we present our results with no subaerial erosion. The iterations we retained were those that yielded modelled 14C profiles within the 2 error bounds of the measured 14C profiles and a modeled surface 10Be concentration within the 2 errors of the measured 10Be concentration. We then used a 2 minimization to evaluate how closely the modeled concentrations match the measured nuclide concentrations with depth. The solution with the lowest 2 value was selected as the most likely burial and exposure history for each site. 4.5 Results and Discussion 4.5.1 Late Holocene moraine ages At its late Holocene maximum, Vintage Glacier covered a maximum area of nearly 1.9 km2, while Lockie’s Glacier covered approximately 1.5 km2. By 2020 CE, both glaciers were less than 0.1 km2 and are now likely too thin to flow. 79 Figure 4.3. Holocene extent change of Lockie’s Glacier. The boundary of Lockie’s Glacier when Core 1 was exposed is uncertain, but must have at least covered the site of Core 3. Four time periods are shown on this figure: 1) The Little Ice Age maximum extent of Lockie’s Glacier; 2) The intermediate extent of Lockie’s Glacier that must have covered the site of Core 3, but left Core 1 exposed; 3) the 2006 CE extent; and 4) the 2017 CE extent of Lockie’s glacier taken from satellite imagery. Moraines sampled in this study are bouldery and draped over bedrock, with crests that were no more than 5 m above the bedrock surface. Collectively, the nine moraine boulders sampled on Vintage Peak 10Be date to 705 ± 219 a (1090-1530 CE, Table 4.1, Fig. 4.4). A single boulder on the Lockie’s Glacier moraine dates to 2849 ± 88 a, while the remaining boulders on the two moraines date between 934 and 538 a. 80 Figure 4.4: Comparison of Vintage moraine record with regional glacier records. Top panel: Probability density function of 10Be moraine ages on Vintage Peak (blue line) with kernel density plot showing individual sample ages. Black markers are the median ± interquartile range 10Be ages from moraines fronting Gilbert Glacier, 65 km north of Vintage Peak (Hawkins et al. 2021). Bottom Panel: Normalized probability from organic in situ 14C samples from glacier forefield and composite morines in western Canada, interpreted to closely constrain periods of glacier advance (modified from Hawkins et al. (2021)). Black line is the mean reconstructed temperature anomaly with respect to the 1800-1900 CE average for the 30 - 60° N latitude band from Kaufman et al. (2020). Blue shaded region is the Little Ice Age period from 1300 - 1850 CE. Differencing a late spring LiDAR survey from one obtained during autumn reveals the depth and distribution of snow across Vintage Peak (Fig. 4.5). In our survey data, moraine crests were either snow-free or were covered by up to several meters less snow than over adjacent areas. This supports the assumption that sampling boulders for surface exposure dating on moraine crests will reduce the influence of snow cover on moraine exposure ages. The LiDAR surveys provide a snapshot of snow distribution and seasonal melt patterns, which are likely consistent each year. Snow depth measured at the nearby snow survey sites was ~25% of normal on March 30th, 2023. An automated snow survey site (3A25P, 1340 m asl) 80 km SE of Vintage 81 Peak similarly recorded the snowpack to be 25% of normal in Mid-May. The influence of snow cover on the young moraine boulders does not significantly impact our interpretations, however, we include exposure ages assuming the moraine boulders were covered with the average snow depth (2.20 m) and density (380 kg m-3) from two nearby snow pillow sites for 6 months of the year in Table 4.1. Glaciers on Vintage Peak reached their Holocene maxima at the beginning of the classic LIA. The timing of glaciers reaching their late Holocene maximum on Vintage Peak accords with other nearby glacier records (Fig. 4.4). Nearby glaciers advanced multiple times from 1.4 ka to 0.3 ka (Hawkins et al., 2021; Koehler and Smith, 2011; Mood and Smith, 2015b; Reyes and Clague, 2004; Ryder and Thomson, 1986). The regional radiocarbon record provides evidence of glacier advance, but not the timing of glaciers reaching their maximum positions and depositing moraines. While there is good correspondence among records of glacier advance preceding the age of moraines on Vintage Peak, the majority of glaciers in the region continued to advance after 1300 CE. Though many glaciers in western Canada reached their greatest Holocene extents between 1600 and 1850 CE (Menounos et al., 2009), there is increasing evidence that glaciers were near their late Holocene maximum positions several times as early as 2 ka (Hawkins et al., 2021; Maurer et al., 2012). Stochastic variability in glacier response to climate forcing, or differing response times to climate perturbations, may explain why glaciers on Vintage Peak reached their greatest Holocene extents slightly earlier than many other glaciers in western Canada (e.g. Roe, 2011). 82 Figure 4.5: Snow depth on Vintage Peak. Top panel is an orthoimage collected on May 16, 2023. Bottom panel shows elevation change between aerial LiDAR-derived DEMs collected May 16, 2013 (snow covered) and September 15, 2017 (snow free). 83 Table 4.1 Vintage peak 10Be sample information. 4.5.2 Paired 14C and 10Be on erratic and bedrock The erratic boulder (17-VP-07), some 65 m beyond the late Holocene moraine that fronts Lockie’s Glacier, returned a 10Be age of 11.65 ± 0.2 ka and a 14C age of 9.74 ± 0.9 ka. When plotted on a paired isotope diagram, the erratic falls within the field of continuous exposure (Fig. 4.6). The discrepancy between the apparent 10Be age and 14C age could be the result of inheritance. The erratic may contain the equivalent of ~2 kyr of inherited 10Be, which is more likely than the boulder having any significant quantity of inherited 14C. Another way to produce the observed difference in 10Be and 14C ages on the erratic is by mass shielding under seasonal snow cover. Mass shielding, in this case by snow, increases the calculated age of a surface. 14C is more sensitive to thin mass shielding, which means that the apparent 14C age increases faster than the apparent 10Be age with increasing shielding correction. With a shielding factor that would be produced by 190 cm of 0.38 g cm-3 snow for 6 months of the year, both nuclides yield an equivalent exposure age of ~14.5 ka. Given the average snow depth of 220 ± 100 cm and average density of 0.38 ± 0.67 g cm-3 (n=123, mean ± standard deviation) at nearby snow monitoring sites which are nearly 600 m a.s.l. lower than 17-VP-07, this snow cover value is reasonable. Work by Menounos et al. (2017) and Seguinot et al. (2016) provide cosmogenic ages and ice sheet modeling that suggest the area around Vintage Peak was likely deglaciated around 14 ka, giving further confidence that the snow-corrected age of the erratic is an accurate deglaciation age for Vintage Peak. However, as we can not rule out the possibility of inheritance, we constrain the age of deglaciation on Vintage Peak to 14.5 ka - 11.6 ka. 85 Figure 4.6: Paired 14C-10Be nuclide diagram with y-axis normalized to the production ratio of 14C-10Be and x-axis (atoms/g 10Be) normalized to the site specific production ratio of 10Be. Dashed lines represent burial isochrons, the uppermost indicating 2 kyr burial, with each descending line representing an additional 2 kyr of burial. The uppermost black solid line is the line of continuous exposure with zero erosion, and lower black line is the line of steady-state equilibrium; samples plotting between these two lines were continuously exposed and eroding at a given rate. The erratic boulder and lack of apparent Younger Dryas moraine(s) indicates that glaciers on Vintage Peak did not extend beyond their late Holocene positions since late Pleistocene deglaciation that may have commenced as early as 14.5 ka. Some cirque glaciers in western Canada that were above the decaying CIS advanced during the Younger Dryas deposited moraines on the landscape (Menounos et al., 2017). However, we did not observe any moraines beyond the late Holocene moraines at Vintage Peak. Our findings on Vintage Peak further 86 support that, on average, glacier advances during the Little Ice Age were the most extensive since the termination of the LGM. Both 10Be and 14C bedrock surface concentrations at each core site follow expected trends based on geomorphic position relative to the modern glacierettes (Fig. 4.7). The greatest nuclide concentrations are found at Core 4 (17-VPC4), outside of the late Holocene moraines, with a surface 10Be concentration of 1.669 ± 0.0167 x 105 atm g-1 and 14C concentration of 1.886 ± 0.445 ×105 atm g-1. Nuclide concentrations at Core 3 (17-VPC3, just inside of the late Holocene moraine) and Core 2 (17-VPC2, just distal from the 2017 CE glacier terminus) are equivalent within errors (Table 4.1). This suggests that, within our measurement uncertainty, the bedrock surface at Cores 2 and 3 experienced nearly equivalent exposure histories over the Holocene. Cores 2 and 3 have 63-66% less 10Be and 55-57% less 14C than Core 4. We interpret equivalent nuclide concentrations at Core 2 and 3 to indicate that Lockie’s Glacier responded to climate perturbations in a roughly binary manner; existing near the glacier’s late Holocene maximum or retracted close to modern extents. 87 Figure 4.7: 10Be and 14C concentrations in bedrock cores with elevation. Nuclide concentration is normalized relative to Core 4, with the lowest production rate of the sample locations. The 10Be and 14C surface concentrations of Core 1 (17-VPC1), which is assumed to have experienced variable cover by minimally erosive ice, but less total ice cover than Core 2 or 3, are 1.179 ± 0.0351 ×105 atm g-1 and 1.401 ± 0.1795 ×105 atm g-1, respectively. Core 1 has 29% less 10 Be and 26% less 14C than Core 4. 4.5.3 14C bedrock depth profiles Five measurements of 14C concentration along the depth of each core show decreasing or consistent, within uncertainties, nuclide concentrations. At Core 4, the 14C concentration is 1.886 ± 0.445 ×105 atm g-1 in the upper 0.05 m of the core, and decreases to 1.464 ± 0.225 ×105 atm g-1 at the bottom 0.40-0.45 m of core (Fig. 4.8, Table D.2). There is no clear pattern to measurement 88 uncertainty with depth in each of the four cores, with measurement uncertainty ranging from 10% to 47% (average 25%) of the mean concentration. Cores 2 and 3 have quasi-vertical 14C concentration profiles of equivalent magnitude, within errors (Fig. 4.8). Core 1 has a slight decrease in 14C concentration in the upper section of core, but within uncertainties is quasi-vertical as well. Figure 4.8: Bedrock 14C profiles with depth. Solid lines are the mean 14C concentration, error bounds are 1-sigma measurement uncertainties at each sample depth. 4.5.4 Monte Carlo simulations The number of Monte Carlo simulations that reproduced our measured nuclide concentration in the bedrock cores are strongly influenced by the prescribed subglacial erosion rate (Fig. D.8-D.15). For all four cores, assuming no subglacial erosion, roughly 50% of 89 simulations produced 14C profiles that fit within error bounds at all measurement depths. For Core 4, assuming no inheritance, no simulations produced acceptable nuclide profiles for any subglacial erosion rate. With 3 kyr of inheritance in the model, 2-6% of simulations could be explained with subglacial erosion rates between 0.05 and 0.2 cm yr -1 for Core 4 (Fig. D.7). For Cores 3 and 2, our modeling could reproduce 14C depth profiles for roughly 10% of simulations with a subglacial erosion rate of 0.05 cm yr-1, though only ~1% of simulations could produce both 14C and 10Be concentrations within measurement errors. Simulations of exposure and burial history at Core 1, when assuming zero erosion, only produced around one acceptable solution for a subglacial erosion rate of 0.05 cm yr-1 and no solutions for greater subglacial erosion rates (Fig. 4.9). With 3 kyr of 10Be inheritance, simulations at Core 1 produced 1% or fewer simulations at subglacial erosion rates of 0.05 to 0.2 cm yr-1. 90 Figure 4.9: Summary of Monte Carlo Results with no inheritance. Each panel shows results assuming different subglacial erosion rates ( ). Where there are no symbols for a given site, there were no results within measured nuclide concentration errors. In summary, our results show that at Core 1 and 4, there is little to no subglacial erosion apparent in the nuclide concentrations. At Cores 2 and 3, a small increase in subglacial erosion would have to be explained by greater total exposure duration and less burial duration by ice. If we assume that Core 2 must have experienced at least as much burial by snow and ice as at Core 1, then we can constrain the burial duration at Core 2 to ≥5 kyr, which then limits the total exposure duration at Core 2 to around 5 kyr. Similarly, any period of burial by ice at Core 3 must have covered Core 2, yielding nearly equivalent total burial and exposure durations. Preserved striations on the bedrock surfaces of Cores 2 and 3 demonstrate that some amount of subglacial erosion did occur, however the subglacial erosion rate was most likely less than 0.05 cm yr-1. 91 Work by Koppes et al. (2015; 2009) determined erosion rates in glaciated catchments to be upwards of 1 cm yr-1 at tidewater glaciers in Alaska and Patagonia, though there is significant subglacial erosion variability dependent on variations in climate, lithology, and glacier bed conditions (Cowton et al., 2012; Graham et al., 2023; Koppes, 2022). Subglacial erosion rates in crystalline bedrock have been constrained from alpine glaciers in the Alps by Goehring et al. (2011) and Wirsig et al. (2017). 4.5.5 Results synthesis and regional context Based on our results and comparison to published glacier records from the region, we conclude that Vintage Peak was fully covered by erosive Cordilleran Ice Sheet during the Last Glacial Maximum. As early as ~14.5 ka, the Cordilleran Ice Sheet had wasted enough to expose most of the bedrock surface around Vintage Peak, though alpine glaciation kept bedrock within the late Holocene glacier maxima covered. Following deglaciation, the site of Core 4 may have been covered by till for ~3 kyr before the till was eroded away. By ca. 11 ka, the site around Core 3 became ice free, with Lockie’s Glacier retreated to modern extents by ca. 10 ka. The Vintage Peak massif remained mostly ice free until on the onset of the Neoglaciation Period around 5-6 ka, at which point Lockie’s Glacier began to advance past the site of Core 2 and Core 3. Lockie’s Glacier fluctuated around its maximum late Holocene extent for nearly 6 kyr until the culmination of the Little Ice Age, which locally peaked around 0.7 ka (ca. 1300 CE), though ice likely stayed in an advanced position until the mid-1800’s CE before beginning to retreat to its minimal modern extent. Decreasing solar insolation in the northern Hemisphere strongly influenced the growth of glaciers in western Canada from the mid-Holocene until the culmination of the Little Ice Age (Menounos et al., 2009; Solomina et al., 2015). Our data on Vintage Peak accords with our 92 current understanding of Holocene glacier change in this region and highlights that glaciers advanced near their late Holocene maxima soon after the onset of Neoglaciation around 6 ka. Holocene temperature reconstruction by Kaufman et al. (2020) shows a peak in Holocene temperatures at around 6.8 ka, followed by a consistent decrease in temperature until around 1850 CE. Our results are at odds with the interpretations of other surface exposure dating studies, namely Menounos et al. (2017) and Darvill et al. (2022). Without correcting for snow cover, the 10 Be ages for the bedrock surface at Core 4 and the erratic accord with the median age of cirque moraines across western Canada formed in response to Younger Dryas cooling (Menounos et al., 2017). However, our paired 14C-10Be dating of Core 4 and the erratic reveal a more complex history of exposure that is best explained by significant influence by snow and in the case of Core 4, inheritance and/or transient mass shielding following deglaciation. Measurements of 14C on boulders presented in Menounos et al. (2017) would show whether the boulders have experienced a complex exposure history. 4.6 Conclusions Our study used paired 14C-10Be dating on alpine bedrock surfaces and bedrock cores to elucidate glacier activity prior to the late Holocene maximum extent of an alpine glacier in the southern Coast Mountains of British Columbia. We estimate late Pleistocene deglaciation at Vintage Peak to have occurred around 14.5 ka. Bedrock outside of Holocene ice extents suggests limited burial that cannot be explained by snow cover alone and may be the result of previous cover by till or inherited 10Be. A broad ridgeline that would have been covered by minimally erosive ice during maximum ice extents indicates around 9.2 ka of exposure and ~4.9 ka of burial since late Pleistocene glaciation. Bedrock exposure histories close to the late Holocene 93 maximum extent and the much less extensive 2017 CE ice margin are nearly equivalent, and indicate that the bedrock at this site may not have become exposed following late Pleistocene deglaciation until 10-11 ka. For most of the Holocene, ice was more extensive than the 2017 CE glacier extent for around 5.6 kyr, likely attributed to glacier advance during Neoglaciation around 6 ka in response to cooling temperatures. Bedrock nuclide concentrations are sensitive to changes in subglacial erosion rates. We can constrain the subglacial erosion rate at three of our sites to likely less than 0.05 cm yr-1. Nine 10Be dates on boulders from two moraines on Vintage Peak yield a median exposure age of 700 ± 220 a (1090-1530 CE), dating the late Holocene maximum extent of glaciers on Vintage Peak to the earlier portion of the Little Ice Age. Future work to test additional simulated burial/exposure histories to investigate the impact of transient till cover and variable pre-LGM inheritance could further constrain glacier change throughout the Holocene. 94 5. Conclusion My dissertation focused on incorporating cosmogenic surface exposure dating to improve our understanding of Holocene glacier behavior in western Canada. First, I combined radiocarbon dating of organic material within a lateral moraine with 10Be surface exposure dating on moraine boulders to develop a refined chronology of Holocene glacier change at Gilbert Glacier in the southern Coast Mountains. Second, I used 10Be surface exposure ages on moraines in the Mackenzie and Selwyn Mountains to develop glacier chronologies on six glaciers, then used several methods to estimate changes in climate since the latest Holocene. Third, I investigated glacier change throughout the Holocene using paired 10Be-14C concentrations from bedrock cores fronting a glacierette in the southern Coast Mountains. This concluding chapter summarizes the major findings of this dissertation and its limitations, and makes recommendations for further research. 5.1 Summary of major findings 5.1.1 Holocene glacier change at Gilbert Glacier My work at Gilbert Glacier demonstrates the utility of applying a multi-proxy approach to glacier reconstructions, specifically using radiocarbon evidence attributed to glacier advance and 10Be exposure dating of moraine boulders to date the onset of glacier retreat. Ryder and Thomson (1986) visited Gilbert Glacier and completed an initial study of the exposed lateral moraine stratigraphy at the site. I add significant detail and chronologic control to the stratigraphic record at Gilbert Glacier, and date the emplacement age of several nested moraines farther down-valley. This study is one of the first to attempt to narrow the age range ascribed to individual advances using tandem dating techniques. I developed a robust glacier chronology from the mid-Holocene to latest Holocene that provided strong evidence that Gilbert Glacier was 95 at or near its late Holocene maximum position by at least 2 ka. I constrain four individual late Holocene advances to 2.0-1.8, ~1.5-1.3, ~0.9-0.8, and 0.4-0.1 ka. This result is corroborated by several other glacier records in western Canada, including our results from Vintage Peak. 5.1.2 Mackenzie and Selwyn mountains In chapter three, I developed late Pleistocene and Holocene glacier chronologies for nine different glaciers in the Mackenzie and Selwyn mountains. This region was understudied and the timing and magnitude of glacier response to Holocene climate change was uncertain (Dyke, 1990; Tomkins et al., 2008). We confirmed that glaciers in the region reached their greatest Holocene extent during the Little Ice Age period (1300-1850 CE), not ca. 2-3 ka, as has been shown in central and northern Alaska (Badding et al., 2013). Our dataset contains notable outliers, which we suggest may record previous periods of moraine development. We use several equilibrium line altitude reconstruction methods to estimate that the climate warmed in northwestern Canada by as much as 2.3 °C since the Little Ice Age. Glacier modeling forced by four different models of past climate estimate that since their Little Ice Age maximum extent, glaciers in this region have already lost nearly half their ice volume by 2019 CE. Projections of future glacier change predict that nearly all glaciers in this region will disappear by the end of the century, with the highest rates of ice loss occurring during the next two decades. Our study provides numerical constraints to landforms and ice limits in a region where few existed. 5.1.3 Vintage Peak In chapter four, we investigated moraine and bedrock exposure histories on Vintage Peak in the southern Coast Mountains of British Columbia. Paired nuclide surface exposure dating on an erratic boulder indicated that Vintage Peak became deglaciated between 14.5 and 11.6 ka, in 96 general agreement with deglaciation ages from other studies (Menounos et al., 2017; Seguinot et al., 2016). Nine moraine boulders yielded 10Be exposure ages that date the maximum Holocene extent of two glaciers to ca. 1300 CE. We also examined if cosmogenic nuclide depth profiles in bedrock could more tightly constrain bedrock exposure histories than surface samples alone. We found that additional isotope measurements in a bedrock profile limit the possible exposure histories that can explain measured nuclide concentrations. Bedrock exposure dating on surfaces within and outside Holocene glacier limits show that glaciers were larger than their 2017 CE extent for much of the last 5 kyr until the onset of post-Little Ice Age retreat. 5.2 Study limitations Several limitations to the work presented in this dissertation exist. The 10Be ages on moraine boulders included notable scatter within most moraines, especially in our dataset presented in Chapter Three. In most cases, there were one or two boulders per moraine that have 2-3 ka of exposure, while the remaining boulders dated to within the last millennium. Collectively, the exposure ages from our entire dataset gives a median exposure age that accords with the LIA maximum extent of many other glaciers in western Canada. However, it is less clear when individual glaciers reached their maximum Holocene extent. In this dissertation, I report the age of individual moraines using the median ± interquartile range of the moraine boulder exposure ages. While this is appropriate given we cannot assume a normally distributed population of exposure ages on the moraine, the conservative error estimate gives larger error bounds than a mean ± standard deviation, complicating the comparison of different moraine ages across regions. Furthermore, tightly constrained chronologies are beneficial when using moraine ages to tune or validate glacier models forced by models of past climate, as I did in Chapter Three. Sampling additional boulders on each moraine to have a larger sample size from which 97 outliers could be identified and using paired nuclide dating on the moraine boulders (i.e. 14C) would allow us to test whether boulders record a complex exposure history. A related limitation of cosmogenic dating is that the dating technique can quickly become cost prohibitive, limiting the number of samples collected and dated. As a general rule, the cost per sample for 10Be surface exposure dating is around $1,000 USD. When applied to moraine dating, an average of 3-5 samples per moraine is commonplace due to this practical limitation. However, this limited sample size per moraine complicates our ability to investigate scatter within a dataset. In Chapter Four, we measured only 14C along the length of the bedrock cores because of both limited quartz abundance and due to the significant additional cost of measuring 10 Be in another 16 samples. Measurements of 10Be along the length of the core would have further constrained the possible burial and exposure histories that could explain our data. Unless the cost of laboratory materials and AMS measurement are significantly reduced, the cost of cosmogenic dating is unlikely to change. I include several recommendations in section 5.4 below that may be implemented, especially if cost is not a concern. We do not attempt to model or constrain every possible influence on the observed nuclide concentration. Namely, we do not directly evaluate the impact of transient till cover on the observed nuclide concentration at Core 4 (outside of Holocene glacier limits). We assume a consistent snow depth overlying each bedrock surface for 6 months of the year during the entire exposure period. If snow depth was significantly different than our assumed depth and density, the inferred exposure and burial durations could be impacted. We provide a brief investigation of variable inheritance from pre-LGM exposure. Ultimately, we were unable to rule out the possibility of inheritance at all sites on Vintage Peak. Similarly, the exposure age on the moraine boulders is several hundred years older than we might expect for moraines deposited during the 98 Little Ice Age in this region. One explanation could be that the moraine boulders contain some inherited 10Be. Dating the boulders with 14C as well would help to investigate the possibility of inheritance in the moraine boulders. 5.3 Broader Implications The work presented in this dissertation provides a meaningful contribution to our understanding of late Pleistocene and Holocene glacier and climate fluctuations in western Canada by exploring the application of multiple geochronologic techniques on glaciated landforms. We developed firm chronologic control of the termination of late Holocene glacier advances in the southern Coast Mountains and in the Mackenzie and Selwyn mountains. The utility of widespread cosmogenic surface exposure dating is clearly demonstrated as an important tool for developing glacier chronologies that can be compared across regional and hemispheric scales. Given the limitations of lichenometry, consideration should be given to apply cosmogenic surface exposure dating in areas that previously have only been dated with lichenometry. Furthermore, using multiple geochronological approaches to reconstruct glacier growth and decay can better constrain past glacier change, which in turn can provide useful constraints to glacier-climate models. Though no single method presented in this dissertation is novel, few previous studies combine tandem approaches to better constrain past glacier activity. Our work in the Yukon and Northwest Territories used prehistoric evidence of past glacier change (i.e. the moraine record) to evaluate GCM simulations of past climate, developing a significant regional glacier chronology in the process. This extends our ability to tune and validate models of past climate beyond the instrumental record. We also showed that there are significant differences between different GCM simulations of past climate. We project future glacier change in the Mackensie and Selwyn mountains and highlight that this region will 99 experience significant continued glacier mass loss over the next several decades, with 85 to >97% ice volume loss by the end of the 21st century. Dating moraines and developing glacier chronologies provides an interesting record of glacier change, but these datasets become useful when applied to reconstruct climate and validate paleoclimate models. The glacial geology community has many geochronologic tools and modeling resources at their disposal; I advocate for researchers to expand their use of these tools to develop more comprehensive records of past glacier and climate change, and projections of future glacier fluctuations. In Chapter Four, I used paired 14C-10Be dating on bedrock surfaces and an erratic, and 10 Be dating on moraine boulders to investigate and constrain complex burial histories and erosion rates in an alpine environment. Measuring cosmogenic nuclide concentrations with depth provides additional constraints on possible exposure histories than surface nuclide concentrations alone. I highlight the utility of paired nuclide dating on moraine and erratic boulders to quantify complex exposure histories. As in Chapter 2, we found strong evidence for glaciers being at or near their late Holocene maximum positions for at least the last 2 kyr and possibly as early as 6 kya. Our findings provide greater context to modern glacier retreat in response to tropospheric warming over the past several decades. Modern tropospheric warming, especially over the past several decades, has increased temperatures to values similar to that of the mid-Holocene (Kaufman et al., 2020). Studies such as those presented in this dissertation provide vital context to the response of the cryosphere to warm temperatures. Linking records of Holocene glacier activity with other evidence of ecosystem change, such as fish abundance and reconstructed hydroclimate, may help us predict future impacts from climate change (Hudson et al., 2021; Lasher et al., 2021; Tunnicliffe et al., 2001). While the rate of modern climate warming is unprecedented in the Holocene, the 100 comparable magnitude of modern warming to conditions during the Holocene Thermal Maximum makes records that include this time period pertinent to addressing current challenges. 5.4 Recommendations Below, I make recommendations for future research that will help to expand on the findings presented in this dissertation. 1. Expand LiDAR coverage in mountainous areas. LiDAR allows geomorphic landforms such as moraines to be identified even in areas of dense vegetation cover. Greater geospatial mapping of glacial geomorphic features can help identify other sites where multiple moraine crests are preserved in mountain valleys (i.e. Johnson et al., 2015). These sites would be promising targets for cosmogenic dating to better understand glacier variability during the mid-to-late Holocene. Making aerial LiDAR data freely available to researchers would allow more users to investigate areas of interest and identify future study locations. 2. Use cosmogenic nuclide exposure dating on previously studied moraine surfaces. There are several moraines in western Canada that have been previously dated using lichenometry, which has numerous inherent methodological limitations (Osborn et al., 2015). These sites suggest complex moraine assemblages that, once robustly dated, may provide useful insight into Holocene glacier fluctuations. Sites such as Tiedemann, Lillooet, and Bridge glaciers preserve moraines that may predate the classic Little Ice Age (Allen and Smith, 2007; Larocque and Smith, 2003; Osborn et al., 2007; Reyes and Clague, 2004). Dating these moraines and comparing moraine ages across multiple sites would give further insight to the nature of late Holocene glacier change in western Canada. 3. Develop mass balance gradients on glaciers in northwestern Canada. Glacier modeling tools like OGGM are profoundly useful in assessing past and future glacier change (Eis et al., 101 2021; Rounce et al., 2023). However their accuracy is limited by a lack of nearby in situ mass balance records that measure the mass balance gradient of individual glaciers (Maussion et al., 2019). Widespread application of remotely-sensed mass glacier mass balance campaigns (Pelto et al., 2019) may offer a logistically feasible method for assessing multi-year mass balance change on remote glaciers. 4. Use terrestrial glacier records to test and validate models of past glacier change. As discussed in Chapter 3, well-constrained glacier chronologies can be used as a proxy of past climate change beyond the instrumental record (e.g. Doughty et al., 2017). In areas where past glacier activity has been confidently reconstructed, a version of our methodology used in Chapter 3 can be used to test whether models of past climate change are supported by the response of glaciers in that region. 5. Standardize multi-nuclide measurements on moraine boulders and bedrock surfaces. Many cosmogenic age studies assume boulders contain no inheritance or effects by erosion or snow cover. Measuring paired nuclides would allow researchers to confirm whether or not sampled boulders have a complex exposure history (Goehring et al., 2022). For Pleistocene and late Holocene studies, regularly measuring 14C and 10Be would provide greater insight to the exposure ages attributed to moraines. 102 Bibliography Allen, S.M., Smith, D.J., 2007. Late Holocene glacial activity of Bridge Glacier, British Columbia Coast Mountains. Can. J. Earth Sci. 44, 1753–1773. Badding, M.E., Briner, J.P., Kaufman, D.S., 2013. 10Be ages of late Pleistocene deglaciation and Neoglaciation in the north-central Brooks Range, Arctic Alaska. J. Quat. Sci. 28, 95–102. Bakke, J., Lie, Ø., Nesje, A., Dahl, S.O., Paasche, Ø., 2005. Utilizing physical sediment variability in glacier-fed lakes for continuous glacier reconstructions during the Holocene, northern Folgefonna, western Norway. Holocene 15, 161–176. Balco, G., 2010. Review of Cosmogenic Nuclides: Principles, Concepts, and Applications in the Earth Sciences. Quaternary Geochronology. https://doi.org/10.1016/j.quageo.2010.08.002 Balco, G., Brown, N., Nichols, K., Venturelli, R.A., Adams, J., Braddock, S., Campbell, S., Goehring, B., Johnson, J.S., Rood, D.H., Wilcken, K., Hall, B., Woodward, J., 2022. Reversible ice sheet thinning in the Amundsen Sea Embayment during the Late Holocene. The Cryosphere. https://doi.org/10.5194/tc-2022-172 Balco, G., Stone, J.O., Lifton, N.A., Dunai, T.J., 2008. A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10Be and 26Al measurements. Quat. Geochronol. 3, 174–195. Balter-Kennedy, A., Young, N.E., Briner, J.P., Graham, B.L., Schaefer, J.M., 2021. Centennial‐ and orbital‐scale erosion beneath the Greenland ice sheet near Jakobshavn Isbræ. J. Geophys. Res. Earth Surf. 126. https://doi.org/10.1029/2021jf006429 Barclay, D.J., Wiles, G.C., Calkin, P.E., 2009. Holocene glacier fluctuations in Alaska. Quat. Sci. Rev. 28, 2034–2048. Beedle, M.J., Menounos, B., Luckman, B.H., Wheate, R., 2009. Annual push moraines as climate proxy. Geophysical Research Letters 36. https://doi.org/10.1029/2009GL039533 Bellefontaine, K., Alldrick, D., Desjardins, P.J., 1994. Mid Coast (all or parts of 92F, G, J, K, L, M, N; 93D; 102P; 103A) (No. Open File 1994-17). British Columbia Ministry of Energy, Mines and Petroleum Resources, British Columbia Geological Survey. Benn, D., Evans, D.J.A., 2014. Glaciers and Glaciation, 2nd edition. Routledge. Biemans, H., Siderius, C., Lutz, A.F., Nepal, S., Ahmad, B., Hassan, T., von Bloh, W., Wijngaard, R.R., Wester, P., Shrestha, A.B., Immerzeel, W.W., 2019. Importance of snow and glacier meltwater for agriculture on the Indo-Gangetic Plain. Nature Sustainability 2, 594–601. Borchers, B., Marrero, S., Balco, G., Caffee, M., Goehring, B., Lifton, N., Nishiizumi, K., Phillips, F., Schaefer, J., Stone, J., 2016. Geological calibration of spallation production rates in the CRONUS-Earth project. Quat. Geochronol. 31, 188–198. 103 Braucher, R., Bourlès, D., Merchel, S., Vidal Romani, J., Fernadez-Mosquera, D., Marti, K., Léanni, L., Chauvet, F., Arnold, M., Aumaître, G., Keddadouche, K., 2013. Determination of muon attenuation lengths in depth profiles from in situ produced cosmogenic nuclides. Nucl. Instrum. Methods Phys. Res. B 294, 484–490. Calkin, P.E., Ellis, J.M., 1980. A Lichenometric Dating Curve and its Application to Holocene Glacier Studies in the Central Brooks Range, Alaska. Arct. Alp. Res. 12, 245–264. Carey, M., Molden, O.C., Rasmussen, M.B., Jackson, M., Nolin, A.W., Mark, B.G., 2017. Impacts of Glacier Recession and Declining Meltwater on Mountain Societies. Ann. Assoc. Am. Geogr. 107, 350–359. Clague, J.J., Menounos, B., Osborn, G., Luckman, B.H., Koch, J., 2009. Nomenclature and resolution in Holocene glacial chronologies. Quat. Sci. Rev. 28, 2231–2238. Coulthard, B., Smith, D.J., Lacourse, T., 2013. Dendroglaciological investigations of mid- to late-Holocene glacial activity in the Mt. Waddington area, British Columbia Coast Mountains, Canada. Holocene 23, 93–103. Cowton, T., Nienow, P., Bartholomew, I., Sole, A., Mair, D., 2012. Rapid erosion beneath the Greenland ice sheet. Geology 40, 343–346. Cuffey, K.M., Paterson, W.S.B., 2010. The Physics of Glaciers. Academic Press. Darvill, C.M., 2013. Cosmogenic nuclide analysis. Geomorphological techniques 1–25. Darvill, C.M., Menounos, B., Goehring, B.M., Lesnek, A.J., 2022. Cordilleran ice sheet stability during the last deglaciation. Geophys. Res. Lett. 49. https://doi.org/10.1029/2021gl097191 Davies, B.J., 2022. Cryospheric Geomorphology: Dating Glacial Landforms II: Radiometric Techniques, in: Treatise on Geomorphology. Elsevier, pp. 249–280. Denton, G.H., Karlén, W., 1977. Holocene glacial and tree-line variations in the White River Valley and Skolai Pass, Alaska and Yukon Territory. Quat. Res. 7, 63–111. Denton, G.H., Karlen, W., 1973. Lichenometry: Its Application to Holocene Moraine Studies in Southern Alaska and Swedish Lapland. Arctic and Alpine Research. https://doi.org/10.2307/1550128 Ditchburn, R.G., Whitehead, N.E., 1994. The separation of 10Be from silicates, in: 3d Workshop of the South Pacific Environmental Radioactivity Association. pp. 4–7. Doughty, A.M., Mackintosh, A.N., Anderson, B.M., Dadic, R., Putnam, A.E., Barrell, D.J.A., Denton, G.H., Chinn, T.J.H., Schaefer, J.M., 2017. An exercise in glacier length modeling: Interannual climatic variability alone cannot explain Holocene glacier fluctuations in New Zealand. Earth Planet. Sci. Lett. 470, 48–53. Dunai, T.J., 2010. Cosmogenic Nuclides: Principles, Concepts and Applications in the Earth Surface Sciences. Cambridge University Press. Dyke, A.S., 1990. A lichenometric study of Holocene rock glaciers and neoglacial moraines, 104 Frances Lake map area, southeastern Yukon Territory and Northwest Territories (No. Bulletin 394). Geological Survey of Canada. Egan, J., Staff, R., Blackford, J., 2015. A high-precision age estimate of the Holocene Plinian eruption of Mount Mazama, Oregon, USA. Holocene 25, 1054–1067. Eis, J., van der Laan, L., Maussion, F., Marzeion, B., 2021. Reconstruction of past glacier changes with an ice-flow glacier model: proof of concept and validation. Front. Earth Sci. 9, 77. Elkadi, J., Lehmann, B., King, G.E., Steinemann, O., Ivy-Ochs, S., Christl, M., Herman, F., 2022. Quantification of post-glacier bedrock surface erosion in the European Alps using 10 Be and optically stimulated luminescence exposure dating. Earth Surface Dynamics. https://doi.org/10.5194/esurf-10-909-2022 Env. and Climate Change Canada, 2022. . Environment and Climate Change Canada Historical Climate Data. URL https://climate.weather.gc.ca (accessed 6.1.22). Eyles, N., Eyles, C.H., Miall, A.D., 1983. Lithofacies types and vertical profile models; an alternative approach to the description and environmental interpretation of glacial diamict and diamictite sequences. Sedimentology 30, 393–410. Geertsema, M., Menounos, B., Bullard, G., Carrivick, J.L., Clague, J.J., Dai, C., Donati, D., Ekstrom, G., Jackson, J.M., Lynett, P., Pichierri, M., Pon, A., Shugar, D.H., Stead, D., Del Bel Belluz, J., Friele, P., Giesbrecht, I., Heathfield, D., Millard, T., Nasonova, S., Schaeffer, A.J., Ward, B.C., Blaney, D., Blaney, E., Brillon, C., Bunn, C., Floyd, W., Higman, B., Hughes, K.E., McInnes, W., Mukherjee, K., Sharp, M.A., 2022. The 28 November 2020 landslide, tsunami, and outburst flood – A hazard cascade associated with rapid deglaciation at Elliot creek, British Columbia, Canada. Geophys. Res. Lett. 49. https://doi.org/10.1029/2021gl096716 Goehring, B.M., Menounos, B., Osborn, G., Hawkins, A., Ward, B., 2022. Reconciling the apparent absence of a Last Glacial Maximum alpine glacial advance, Yukon Territory, Canada, through cosmogenic beryllium-10 and carbon-14 measurements. Geochronology 4, 311–322. Goehring, B.M., Muzikar, P., Lifton, N.A., 2013. An in situ 14C–10Be Bayesian isochron approach for interpreting complex glacial histories. Quat. Geochronol. 15, 61–66. Goehring, B.M., Schaefer, J.M., Schluechter, C., Lifton, N.A., Finkel, R.C., Timothy Jull, A.J., Akçar, N., Alley, R.B., 2011. The Rhone Glacier was smaller than today for most of the Holocene. Geology 39, 679–682. Goehring, B.M., Wilson, J., Nichols, K., 2019. A fully automated system for the extraction of in situ cosmogenic carbon-14 in the Tulane University cosmogenic nuclide laboratory. Nucl. Instrum. Methods Phys. Res. B 455, 284–292. Gosse, J.C., Phillips, F.M., 2001. Terrestrial in situ cosmogenic nuclides: theory and application. Quat. Sci. Rev. 20, 1475–1560. 105 Graham, B.L., Briner, J.P., Young, N.E., Balter-Kennedy, A., Koppes, M., Schaefer, J.M., Poinar, K., Thomas, E.K., 2023. In situ 10Be modeling and terrain analysis constrain subglacial quarrying and abrasion at Jakobshavn Isbræ, Greenland. EGUsphere. https://doi.org/10.5194/egusphere-2023-836 Haeberli, W., Brandova, D., Burga, C., Egli, M., Frauenfelder, R., Kääb, A., Maisch, M., 2003. Methods for absolute and relative age dating of rock-glacier surfaces in alpine permafrost, in: Phillips, Springman, Arenson (Eds.), Proceedings of the Eighth International Conference on Permafrost. Swets & Zeitlinger, Lisse, pp. 343–348. Hawkins, A.C., Menounos, B., Goehring, B.M., Osborn, G.D., Clague, J.J., Jensen, B., 2021. Tandem dating methods constrain late Holocene glacier advances, southern Coast Mountains, British Columbia. Quat. Sci. Rev. 274, 107282. Heyman, J., Applegate, P.J., Blomdin, R., Gribenski, N., Harbor, J.M., Stroeven, A.P., 2016. Boulder height – exposure age relationships from a global glacial 10Be compilation. Quat. Geochronol. 34, 1–11. Hippe, K., 2017. Constraining processes of landscape change with combined in situ cosmogenic 14C-10Be analysis. Quat. Sci. Rev. 173, 1–19. Hock, R., Rasul, G., Adler, C., Caceres, B., Gruber, S., Hirabayashi, Y., Jackson, M., Kääb, A., Kang, S., Kutuzov, S., Milner, A., Molau, U., Morin, S., Orlove, B., Steltzer, H., Allen, S., Arenson, L., Baneerjee, S., Barr, I., Bórquez, R., Brown, L., Cao, B., Carey, M., Cogley, G., Fischlin, A., de Sherbinin, A., Eckert, N., Geertsema, M., Hagenstad, M., Honsberg, M., Hood, E., Huss, M., Jimenez Zamora, E., Kotlarski, S., Lefeuvre, P., Ignacio López Moreno, J., Lundquist, J., McDowell, G., Mills, S., Mou, C., Nepal, S., Noetzli, J., Palazzi, E., Pepin, N., Rixen, C., Shahgedanova, M., McKenzie Skiles, S., Vincent, C., Viviroli, D., Gesa, A.W., Yangjee Sherpa, P., Weyer, N., Wouters, B., Yasunari, T., You, Q., Zhang, Y., 2019. High mountain areas 131–202. Hudson, A.M., Emery-Wetherell, M.M., Lubinski, P.M., Butler, V.L., Grimstead, D.N., Jenkins, D.L., 2021. Reconstructing paleohydrology in the northwest Great Basin since the last deglaciation using Paisley Caves fish remains (Oregon, U.S.A.). Quat. Sci. Rev. 262, 106936. Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss, M., Dussaillant, I., Brun, F., Kääb, A., 2021. Accelerated global glacier mass loss in the early twenty-first century. Nature 592, 726–731. Immerzeel, W.W., Lutz, A.F., Andrade, M., Bahl, A., Biemans, H., Bolch, T., Hyde, S., Brumby, S., Davies, B.J., Elmore, A.C., Emmer, A., Feng, M., Fernández, A., Haritashya, U., Kargel, J.S., Koppes, M., Kraaijenbrink, P.D.A., Kulkarni, A.V., Mayewski, P.A., Nepal, S., Pacheco, P., Painter, T.H., Pellicciotti, F., Rajaram, H., Rupper, S., Sinisalo, A., Shrestha, A.B., Viviroli, D., Wada, Y., Xiao, C., Yao, T., Baillie, J.E.M., 2020. Importance and vulnerability of the world’s water towers. Nature 577, 364–369. Ivy-Ochs, S., Kerschner, H., Kubik, P.W., Schlüchter, C., 2006. Glacier response in the European Alps to Heinrich Event 1 cooling: The Gschnitz stadial. J. Quat. Sci. 21, 115–130. 106 Jarosch, A.H., Schoof, C.G., Anslow, F.S., 2013. Restoring mass conservation to shallow ice flow models over complex terrain. The Cryosphere. https://doi.org/10.5194/tc-7-229-2013 Jensen, B.J.L., Beaudoin, A.B., Clynne, M.A., Harvey, J., Vallance, J.W., 2019. A re-examination of the three most prominent Holocene tephra deposits in western Canada: Bridge River, Mount St. Helens Yn and Mazama. Quat. Intern. 500, 83–95. Johnson, M.D., Fredin, O., Ojala, A.E.K., Peterson, G., 2015. Unraveling Scandinavian geomorphology: the LiDAR revolution. GFF 137, 245–251. Jurt, C., Brugger, J., Dunbar, K. w., Milch, K., Orlove, B., 2015. Cultural values of glaciers, in: Huggel, C., Carey, M., Clague, J.J., Kaab, A. (Eds.), The High-Mountain Cryosphere. Cambridge University Press, Cambridge, pp. 90–106. Kaplan, M.R., Schaefer, J.M., Strelin, J.A., Denton, G.H., Anderson, R.F., Vandergoes, M.J., Finkel, R.C., Schwartz, R., Travis, S.G., Garcia, J.L., Martini, M.A., Nielsen, S.H.H., 2016. Patagonian and southern South Atlantic view of Holocene climate. Quat. Sci. Rev. 141, 112–125. Kaplan, M.R., Strelin, J.A., Schaefer, J.M., Denton, G.H., Finkel, R.C., Schwartz, R., Putnam, A.E., Vandergoes, M.J., Goehring, B.M., Travis, S.G., 2011. In-situ cosmogenic 10Be production rate at Lago Argentino, Patagonia: Implications for late-glacial climate chronology. Earth Planet. Sci. Lett. 309, 21–32. Karlen, W., Matthews, J.A., 1992. Reconstructing Holocene Glacier Variations from Glacial Lake Sediments: Studies from Nordvestlandet and Jostedalsbreen-Jotunheimen, Southern Norway. Geografiska Annaler. Series A, Physical Geography. https://doi.org/10.2307/521430 Kaufman, D., McKay, N., Routson, C., Erb, M., Dätwyler, C., Sommer, P.S., Heiri, O., Davis, B., 2020. Holocene global mean surface temperature, a multi-method reconstruction approach. Sci Data 7, 201. Koch, J., 2009. Improving age estimates for late Holocene glacial landforms using dendrochronology – Some examples from Garibaldi Provincial Park, British Columbia. Quat. Geochron. https://doi.org/10.1016/j.quageo.2008.11.002 Koch, J., Clague, J.J., Osborn, G., 2014. Alpine glaciers and permanent ice and snow patches in western Canada approach their smallest sizes since the mid-Holocene, consistent with global trends. Holocene 24, 1639–1648. Koch, J., Clague, J.J., Osborn, G.D., 2007a. Glacier fluctuations during the past millennium in Garibaldi Provincial Park, southern Coast Mountains, British Columbia. Can. J. Earth Sci. 44, 1215–1233. Koch, J., Osborn, G.D., Clague, J.J., 2007b. Pre-`Little Ice Age’ glacier fluctuations in Garibaldi Provincial Park, Coast Mountains, British Columbia, Canada. Holocene 17, 1069–1078. Koehler, L., Smith, D.J., 2011. Late Holocene glacial activity in Manatee Valley, southern Coast Mountains, British Columbia, Canada. Can. J. Earth Sci. 48, 603–618. 107 Koppes, M., Hallet, B., Rignot, E., Mouginot, J., Wellner, J.S., Boldt, K., 2015. Observed latitudinal variations in erosion as a function of glacier dynamics. Nature 526, 100–103. Koppes, M.N., 2022. 4.09 - Rates and Processes of Glacial Erosion, in: Shroder, J. (jack) F. (Ed.), Treatise on Geomorphology (Second Edition). Academic Press, Oxford, pp. 169–181. Koppes, M.N., Montgomery, D.R., 2009. The relative efficacy of fluvial and glacial erosion over modern to orogenic timescales. Nat. Geosci. 2, 644–647. Kronig, O., Ivy-Ochs, S., Hajdas, I., Christl, M., Wirsig, C., Schlüchter, C., 2018. Holocene evolution of the Triftje- and the Oberseegletscher (Swiss Alps) constrained with 10Be exposure and radiocarbon dating. Swiss J. Geosci. 111, 117–131. Larocque, S.J., Smith, D.J., 2004. Calibrated Rhizocarpon spp. Growth Curve for the Mount Waddington Area, British Columbia Coast Mountains, Canada. Arct. Antarct. Alp. Res. 36, 407–418. Larocque, S.J., Smith, D.J., 2003. Little Ice Age glacial activity in the Mt. Waddington area, British Columbia Coast Mountains, Canada. Can. J. Earth Sci. 40, 1413–1436. Lasher, G.E., Abbott, M.B., Anderson, L., Yasarer, L., Rosenmeier, M., Finney, B.P., 2021. Holocene hydroclimatic reorganizations in northwest Canada inferred from lacustrine carbonate oxygen isotopes. Geophys. Res. Lett. 48. https://doi.org/10.1029/2021gl092948 Lehmann, B., Herman, F., Valla, P.G., King, G.E., Biswas, R.H., Ivy-Ochs, S., Steinemann, O., Christl, M., 2020. Postglacial erosion of bedrock surfaces and deglaciation timing: New insights from the Mont Blanc massif (western Alps). Geology 48, 139–144. Lemelin, H., Dawson, J., Stewart, E.J., Maher, P., Lueck, M., 2010. Last-chance tourism: the boom, doom, and gloom of visiting vanishing destinations. Curr. Issues Tourism 13, 477–493. Le Roy, M., Deline, P., Carcaillet, J., Schimmelpfennig, I., Ermini, M., 2017. 10Be exposure dating of the timing of Neoglacial glacier advances in the Ecrins-Pelvoux massif, southern French Alps. Quat. Sci. Rev. 178, 118–138. Lifton, N., Sato, T., Dunai, T.J., 2014. Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes. Earth Planet. Sci. Lett. 386, 149–160. Lowe, J.J., Walker, M., 2014. Reconstructing Quaternary Environments, Third Edition. ed. Routledge, 2 Park Square, Milton Park, Abingdon, Oxon. Luckman, B.H., 2007. Reconstruction of Little Ice Age events in the Canadian Rocky Mountains. Géogr. phys. quat. 40, 17–28. Luckman, B.H., 1988. Dating the Moraines and Recession of Athabasca and Dome Glaciers, Alberta, Canada. Arctic and Alpine Research. https://doi.org/10.2307/1551697 Luckman, B.H., 1977. Lichenometric dating of holocene moraines at Mount Edith Cavell, Jasper, 108 Alberta. Canadian Journal of Earth Sciences. https://doi.org/10.1139/e77-154 Luckman, B.H., Holdsworth, G., Osborn, G.D., 1993. Neoglacial Glacier Fluctuations in the Canadian Rockies. Quat. Res. 39, 144–153. Luckman, B.H., Masiokas, M.H., 2017. Neoglacial history of Robson Glacier, British Columbia. Canadian Journal of. Luckman, B.H., Sperling, B.J.R., Osborn, G.D., 2020. The Holocene history of the Columbia Icefield, Canada. Quat. Sci. Rev. 242, 106436. Magrani, F., Valla, P.G., Egholm, D., 2022. Modelling alpine glacier geometry and subglacial erosion patterns in response to contrasting climatic forcing. Earth Surf. Processes Landforms 47, 1054–1072. Masson-Delmotte, V., Zhai, P., Pörtner, H.-O., Roberts, D., Skea, J., Shukla, P.R., Pirani, A., Moufouma-Okia, W., Péan, C., Pidcock, R., Others, 2018. Global warming of 1.5 C. An IPCC Special Report on the impacts of global warming of 1. Mathews, W.H., 1951. Historic and prehistoric fluctuations of alpine glaciers in the mount Garibaldi map-area, southwestern British Columbia. J. Geol. 59, 357–380. Maurer, M.K., Menounos, B., Luckman, B.H., Osborn, G., Clague, J.J., Beedle, M.J., Smith, R., Atkinson, N., 2012. Late Holocene glacier expansion in the Cariboo and northern Rocky Mountains, British Columbia, Canada. Quat. Sci. Rev. 51, 71–80. Maussion, F., Butenko, A., Champollion, N., Dusch, M., Eis, J., Fourteau, K., Gregor, P., Jarosch, A.H., Landmann, J., Oesterle, F., Recinos, B., Rothenpieler, T., Vlug, A., Wild, C.T., Marzeion, B., 2019. The Open Global Glacier Model (OGGM) v1.1. Geosci. Model Dev. 12, 909–931. Menounos, B., Clague, J.J., Clarke, G.K.C., Marcott, S.A., Osborn, G., Clark, P.U., Tennant, C., Novak, A.M., 2013. Did rock avalanche deposits modulate the late Holocene advance of Tiedemann Glacier, southern Coast Mountains, British Columbia, Canada? Earth Planet. Sci. Lett. 384, 154–164. Menounos, B., Clague, J.J., Osborn, G., Luckman, B.H., Lakeman, T.R., Minkus, R., 2008. Western Canadian glaciers advance in concert with climate change circa 4.2 ka. Geophys. Res. Lett. 35. https://doi.org/10.1029/2008GL033172 Menounos, B., Goehring, B.M., Osborn, G., Margold, M., Ward, B., Bond, J., Clarke, G.K.C., Clague, J.J., Lakeman, T., Koch, J., Caffee, M.W., Gosse, J., Stroeven, A.P., Seguinot, J., Heyman, J., 2017. Cordilleran Ice Sheet mass loss preceded climate reversals near the Pleistocene Termination. Science 358, 781–784. Menounos, B., Hugonnet, R., Shean, D., Gardner, A., Howat, I., Berthier, E., Pelto, B., Tennant, C., Shea, J., Noh, M., Brun, F., Dehecq, A., 2019. Heterogeneous Changes in Western North American Glaciers Linked to Decadal Variability in Zonal Wind Strength. Geophys. Res. Lett. 46, 200–209. 109 Menounos, B., Koch, J., Osborn, G., Clague, J.J., Mazzucchi, D., 2004. Early Holocene glacier advance, southern Coast Mountains, British Columbia, Canada. Quat. Sci. Rev. 23, 1543–1550. Menounos, B., Maurer, L., Clague, J.J., Osborn, G., 2020. Late Holocene fluctuations of Stoppani Glacier, southernmost Patagonia. Quat. Res. 95, 56–64. Menounos, B., Osborn, G., Clague, J.J., Luckman, B.H., 2009. Latest Pleistocene and Holocene glacier fluctuations in western Canada. Quat. Sci. Rev. 28, 2049–2074. Millan, R., Mouginot, J., Rabatel, A., Morlighem, M., 2022. Ice velocity and thickness of the world’s glaciers. Nat. Geosci. 15, 124–129. Mood, B.J., Smith, D.J., 2015a. Holocene glacier activity in the British Columbia Coast Mountains, Canada. Quat. Sci. Rev. 128, 14–36. Mood, B.J., Smith, D.J., 2015b. Latest Pleistocene and Holocene behaviour of Franklin Glacier, Mt. Waddington area, British Columbia Coast Mountains, Canada. Holocene 25, 784–794. Moore, R.D., Fleming, S.W., Menounos, B., Wheate, R., Fountain, A., Stahl, K., Holm, K., Jakob, M., 2009. Glacier change in western North America: influences on hydrology, geomorphic hazards and water quality. Hydrol. Process. 23, 42–61. Moore, R.D., Pelto, B., Menounos, B., Hutchinson, D., 2020. Detecting the Effects of Sustained Glacier Wastage on Streamflow in Variably Glacierized Catchments. Front Earth Sci. Chin. 8. https://doi.org/10.3389/feart.2020.00136 Munroe, J.S., Laabs, B.J.C., 2017. Combining radiocarbon and cosmogenic ages to constrain the timing of the last glacial-interglacial transition in the Uinta Mountains, Utah, USA. Geology 45, 171–174. Nichols, K.A., Goehring, B.M., 2019. Isolation of quartz for cosmogenic in situ 14C analysis. Geochronology 1, 43–52. Nishiizumi, K., Imamura, M., Caffee, M.W., Southon, J.R., Finkel, R.C., McAninch, J., 2007. Absolute calibration of 10Be AMS standards. Nucl. Instrum. Methods Phys. Res. B 258, 403–413. Oerlemans, J., 2005. Extracting a climate signal from 169 glacier records. Science 308, 675–677. O’Neel, S., Hood, E., Bidlack, A.L., Fleming, S.W., Arimitsu, M.L., Arendt, A., Burgess, E., Sergeant, C.J., Beaudreau, A.H., Timm, K., Hayward, G.D., Reynolds, J.H., Pyare, S., 2015. Icefield-to-Ocean Linkages across the Northern Pacific Coastal Temperate Rainforest Ecosystem. Bioscience 65, 499–512. Osborn, G., Luckman, B.H., 1988. Holocene glacier fluctuations in the Canadian Cordillera (Alberta and British Columbia). Quat. Sci. Rev. 7, 115–128. Osborn, G., McCarthy, D., LaBrie, A., Burke, R., 2015. Lichenometric dating: Science or pseudo-science? Quat. Res. 83, 1–12. 110 Osborn, G., Menounos, B., Koch, J., Clague, J.J., Vallis, V., 2007. Multi-proxy record of Holocene glacial history of the Spearhead and Fitzsimmons ranges, southern Coast Mountains, British Columbia. Quat. Sci. Rev. 26, 479–493. Osborn, G., Menounos, B., Ryane, C., Riedel, J., Clague, J.J., Koch, J., Clark, D., Scott, K., Davis, P.T., 2012. Latest Pleistocene and Holocene glacier fluctuations on Mount Baker, Washington. Quat. Sci. Rev. 49, 33–51. Osborn, G., Taylor, J., 1975. Lichenometry on Calcareous Substrates in the Canadian Rockies. Quat. Res. 5, 111–120. Pelto, B.M., Menounos, B., Marshall, S.J., 2019. Multi-year evaluation of airborne geodetic surveys to estimate seasonal mass balance, Columbia and Rocky Mountains, Canada. Cryosph. Discuss. 1–30. Pelto, M.S., Hedlund, C., 2001. Terminus behavior and response time of North Cascade glaciers, Washington, U.S.A. J. Glaciol. 47, 497–506. Pfeffer, W.T., Tad Pfeffer, W., Arendt, A.A., Bliss, A., Bolch, T., Graham Cogley, J., Gardner, A.S., Hagen, J.-O., Hock, R., Kaser, G., Kienholz, C., Miles, E.S., Moholdt, G., Mölg, N., Paul, F., Radić, V., Rastner, P., Raup, B.H., Rich, J., Sharp, M.J., The Randolph Consortium, 2014. The Randolph Glacier Inventory: a globally complete inventory of glaciers. Journal of Glaciology. https://doi.org/10.3189/2014jog13j176 Pitman, K.J., Moore, J.W., Sloat, M.R., Beaudreau, A.H., Bidlack, A.L., Brenner, R.E., Hood, E.W., Pess, G.R., Mantua, N.J., Milner, A.M., Radić, V., Reeves, G.H., Schindler, D.E., Whited, D.C., 2020. Glacier Retreat and Pacific Salmon. Bioscience 70, 220–236. Porter, S.C., 1989. Late Holocene fluctuations of the fiord glacier system in Icy Bay, Alaska, USA. Arct. Alp. Res. https://doi.org/10.1080/00040851.1989.12002750 Purdie, H., 2013. Glacier Retreat and Tourism: Insights from New Zealand. mred 33, 463–472. Rand, C., Goehring, B.M., 2019. The distribution and magnitude of subglacial erosion on millennial timescales at Engabreen, Norway. Ann. Glaciol. 60, 73–81. Reimer, P.J., Austin, W.E.N., Bard, E., Bayliss, A., Blackwell, P.G., Ramsey, C.B., Butzin, M., Cheng, H., Lawrence Edwards, R., Friedrich, M., Grootes, P.M., Guilderson, T.P., Hajdas, I., Heaton, T.J., Hogg, A.G., Hughen, K.A., Kromer, B., Manning, S.W., Muscheler, R., Palmer, J.G., Pearson, C., van der Plicht, J., Reimer, R.W., Richards, D.A., Marian Scott, E., Southon, J.R., Turney, C.S.M., Wacker, L., Adolphi, F., Büntgen, U., Capano, M., Fahrni, S.M., Fogtmann-Schulz, A., Friedrich, R., Köhler, P., Kudsk, S., Miyake, F., Olsen, J., Reinig, F., Sakamoto, M., Sookdeo, A., Talamo, S., 2020. The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0–55 cal kBP). Radiocarbon 62, 725–757. Reyes, A.V., Clague, J.J., 2004. Stratigraphic evidence for multiple Holocene advances of Lillooet Glacier, southern Coast Mountains, British Columbia. Can. J. Earth Sci. 41, 903–918. Reyes, A.V., Luckman, B.H., Smith, D.J., Clague, J.J., Van Dorp, R.D., 2006a. Tree-Ring Dates 111 for the Maximum Little Ice Age Advance of Kaskawulsh Glacier, St. Elias Mountains, Canada. Arctic 59, 14–20. Reyes, A.V., Wiles, G.C., Smith, D.J., Barclay, D.J., Allen, S., Jackson, S., Larocque, S., Laxton, S., Lewis, D., Calkin, P.E., Clague, J.J., 2006b. Expansion of alpine glaciers in Pacific North America in the first millennium A.D. Geology 34, 57–60. Roe, G.H., 2011. What do glaciers tell us about climate variability and climate change? J. Glaciol. 57, 567–578. Røthe, T.O., Bakke, J., Støren, E.W.N., Bradley, R.S., 2018. Reconstructing Holocene Glacier and Climate Fluctuations From Lake Sediments in Vårfluesjøen, Northern Spitsbergen. Front Earth Sci. Chin. 6. https://doi.org/10.3389/feart.2018.00091 Rounce, D.R., Hock, R., Maussion, F., Hugonnet, R., Kochtitzky, W., Huss, M., Berthier, E., Brinkerhoff, D., Compagno, L., Copland, L., Farinotti, D., Menounos, B., McNabb, R.W., 2023. Global glacier change in the 21st century: Every increase in temperature matters. Science 379, 78–83. Ryder, J.M., Thomson, B., 1986. Neoglaciation in the southern Coast Mountains of British Columbia: chronology prior to the late Neoglacial maximum. Canadian Journal of Earth Sciences. https://doi.org/10.1139/e86-031 Schaefer, J.M., Denton, G.H., Kaplan, M., Putnam, A., Finkel, R.C., Barrell, D.J.A., Andersen, B.G., Schwartz, R., Mackintosh, A., Chinn, T., Schlüchter, C., 2009. High-frequency Holocene glacier fluctuations in New Zealand differ from the northern signature. Science 324, 622–625. Schmidt, E., 1951. A non-destructive concrete tester. Concrete 59, 34–35. Schweingruber, F.H., 2012. Tree Rings: Basics and Applications of Dendrochronology. Springer Science & Business Media. Seguinot, J., Rogozhina, I., Stroeven, A.P., Margold, M., Kleman, J., 2016. Numerical simulations of the Cordilleran Ice Sheet through the last glacial cycle. Cryosphere 10, 639–664. Solomina, O.N., Bradley, R.S., Hodgson, D.A., Ivy-Ochs, S., Jomelli, V., Mackintosh, A.N., Nesje, A., Owen, L.A., Wanner, H., Wiles, G.C., Young, N.E., 2015. Holocene glacier fluctuations. Quat. Sci. Rev. 111, 9–34. Stuiver, M., Reimer, P.J., Reimer, R.W., 2021. CALIB 8.1 [WWW program] [WWW Document]. URL http://calib.org (accessed 3.1.21). Tomkins, J.D., Lamoureux, S.F., Sauchyn, D.J., 2008. Reconstruction of climate and glacial history based on a comparison of varve and tree-ring records from Mirror Lake, Northwest Territories, Canada. Quat. Sci. Rev. 27, 1426–1441. Tunnicliffe, V., O’Connell, J.M., McQuoid, M.R., 2001. A Holocene record of marine fish remains from the Northeastern Pacific. Mar. Geol. 174, 197–210. 112 Wayne, W.J., 1984. Relative dating techniques to distinguish late Pleistocene-Holocene continental sediments. Chem. Geol. 44, 337–348. Wirsig, C., Ivy-Ochs, S., Reitner, J.M., Christl, M., Vockenhuber, C., Bichler, M., Reindl, M., 2017. Subglacial abrasion rates at Goldbergkees, Hohe Tauern, Austria, determined from cosmogenic10Be and36Cl concentrations. Earth Surface Processes and Landforms. https://doi.org/10.1002/esp.4093 Yi, C., Chen, H., Yang, J., Liu, B., Fu, P., Liu, K., Li, S., 2008. Review of Holocene glacial chronologies based on radiocarbon dating in Tibet and its surrounding mountains. J. Quat. Sci. 23, 533–543. Young, N.E., Lesnek, A.J., Cuzzone, J.K., Briner, J.P., 2021. In situ cosmogenic 10Be–14C–26Al measurements from recently deglaciated bedrock as a new tool to decipher changes in Greenland Ice Sheet size. Clim. Past. 113 Appendix Appendix A. Author contributions Authorship contributions are listed following the CRediT Authorship Guidelines. Chapter 2: A. Hawkins contributed to all authorship components except resources and funding acquisition. B. Menounos contributed to all components. J. Clague and G. Osborn both contributed investigation, and review and editing. B. Goehring contributed to significant laboratory investigation, supervision, resources and review and editing of the manuscript. B. Jensen contributed to investigation and writing of portions of the original and final manuscript. Chapter 3: A. Hawkins contributed to all 14 authorship components except resources and supervision. B. Menounos was involved in all authorship components. B. Goehring contributed to formal analysis, investigation, resources, supervision, validation, and review/editing. G. Osborn was involved in conceptualization, investigation, supervision, and review/editing. B. Pelto contributed to data curation, methodology, and software. C. Darvill was involved in investigation, visualization, and review/editing. J. Schaefer was involved in conceptualization, funding acquisition, investigation, and review/editing. Chapter 4: A. Hawkins contributed to all authorship components except funding acquisition, administration, resources, supervision, and validation. B. Goehring contributed to all components except funding acquisition, investigation, and original draft. B. Menounos contributed to all components except data curation, formal analysis, software, and validation. 114 Appendix B. Chapter 2 Supplement B.1 TCN field sample collection Samples for TCN dating were collected using a gas-powered concrete saw and hammer and chisel from prominent boulders on or very near the crest of moraines, where there is assumed to be less chance of post-depositional movement or exhumation and less burial by snow. The late Holocene moraines were poorly to moderately matrix supported, with abundant large (>1 m 3) clasts at the moraine crest surface. Sampling methodology largely followed the recommendations of Darvill (2013). Sample locations were marked with a handheld GPS with approximate horizontal accuracy of 3-5 m and vertical accuracy of 5 m. We measured the dimensions of each boulder with a tape measurer and took photos from various aspects of each boulder. At each site, using a Brunton inclinometer, horizon shielding was measured at approximately 45-degree increments (i.e N, NW, W…), with measurements aimed at major topographic inflections (valley bottoms, mountain tops, etc.). Approximately 1.5 kg of rock was collected from each boulder, due to low modal abundance of quartz in the boulders and assumed low isotope ratio due to the short exposure duration of the samples. All sub-samples from each boulder were collected into a labeled cloth bag and transported in 5 gallon buckets for laboratory analysis. B.2 Cosmogenic laboratory analysis Field samples had been cut into 3 x 3 x 2.5 cm (L x W x thickness) prisms. We measured each sub-sample thickness with calipers before they were crushed and ground at the start of the physical separation process. The average thickness of all pieces is listed in Table D.1. Physical and chemical preparation of each rock sample to quartz followed procedures outlined in (Nichols 115 and Goehring, 2019). Each quartz sample was spiked with Be carrier and underwent standard Be extraction procedures (Ditchburn and Whitehead, 1994). Final BeO samples were loaded into cathodes and sent to Purdue University’s PRIME lab for AMS measurement. Exposure ages are reported using LSDn scaling, assuming no significant snow cover or erosion (Table D.1). Snow cover is poorly constrained at this location and erosion is likely minimal, given the relatively short exposure age of the samples and resistant rock type. Analytical errors are reported for individual samples, while the reported age of each separate landform (i.e. distinct moraine ridge) is the median ± interquartile range (IQR). B.3 Regional radiocarbon compilation The compiled radiocarbon record from sites within 350 km of Gilbert Glacier in-cludes over 100 radiocarbon dates from overridden, in situ wood (Table D.4). Much of the record has been previously compiled by Mood and Smith (2015a), though our compilation is more restrictive, focusing on samples that closely limit the timing of glacier advance. Data from this study and the regional compilation were combined to generate the probability distribution in Fig. 2.4. Both individual samples, such as a leaf litter (sample “S-1572”) layer in the lateral stratigraphy at Gilbert Glacier from Ryder and Thomson (1986) and regional trends, such as the “First Millennium Advance” originally proposed by Reyes et al. (2006b), are used to best approximate the timing of late Holocene advances at Gilbert Glacier. B.4 Electron microprobe analysis A bulk sample of peat from the top of unit 2 in Gully A (see main text) had rare glass shards and was sent to the University of Alberta for isolation of glass shards and electron probe microanalysis. The sample was sieved to two size fractions: >150 µm and 150-25 µm. The finer 116 fraction was used for the identification of glass shards. The sample was then density separated and a 5 µm beam used during electron probe microanalysis, as detailed in Jensen et al. (2019). Geochemical comparison between tephra found in the sampled peat and known tephras in the area allowed for the identification of the sampled tephra as Mt Mazama (Fig. B.2). 117 Table B.1: Cosmogenic nuclide input data Table B.2: Cosmogenic nuclide sample blank information Table B.3: Radiocarbon ages from this study 117 Table B.4: Regional radiocarbon compilation Table B.4 (continued) 118 119 Table B.4 (continued) 120 Figure B.1: Oxide plots showing geochemical differences between Mazama, Bridge River (Mt Meager) and Mt St Helens tephras found in western Canada. 121 Figure B.2: Oxide plots of tephra sampled from peat at Gilbert Glacier and Mt Mazama representative tephra. The tephra at Gilbert Glacier matches well with Mt Mazama and are interpreted to come from the same event. 122 Appendix C. Chapter 3 Supplement Table C.1: Airphoto and satellite information Table C.2: Cosmogenic sample blank information. 123 Table C.3: 10Be ages with snow correction. 124 Figure C.1: CRU TS 4.06 annual surface temperature anomalies from our field sites around Nahanni National Park Reserve and Keele Peak. Anomalies are with respect to the average annual temperature between 1981 and 2010 CE. Note that since the early 1900’s, there has been over 2 °C of warming in this region. Figure C.2: Comparison between glacier model outputs using a glacier model modified from Jarosch et al.(2013) and OGGM outputs. Both simulations suggest a delayed onset of LIA glacier advance inconsistent with the moraine record. 125 Figure C.3: Example of temperature and precipitation bias calibration results from Anderson Glacier forced by CCSM4 climate. Colored contours represent the root mean squared error (RMSE) between modeled glacier length and observed length in the moraine or remote sensing record. Note that changes in precipitation can theoretically be compensated for by changes in temperature (and vice versa). Figure C.4: Example of glacier volume and length evolution in a non-transient 1000-year equilibrium run for Butterfly Glacier with a temperature bias of -2.28 °C. 126 Figure C.5: Modeled glacier profiles during non-transient 1000-year runs in OGGM. Note the glacier thickening, which serves to reduce the apparent difference in ELA when using the AAR method of ELA reconstruction. 127 Figure C.6: Modeled glacier hypsometries from non-transient equilibrium experiments in OGGM. Dashed lines mark the intersection of the estimated ELA assuming an AAR of 0.6, and the elevation of the ELA. 128 Figure C.7: Looking Southeast at North Moraine Hill Glacier, with sampled boulders highlighted in green. Note the person (red) on the left side of the image for an approximate scale. Figures C.8-11: 30-year rolling average precipitation and temperature for each GCM and resulting OGGM regional glacier volume over the past 1 kyr. 129 130 131 Figure C.12: Manually digitized glacier margins from aerial photos and satellite imagery. Margins are colored by the decade from which they are digitized. Note that not all margins cover the entire glacier extent, as the focus was to map the terminal position of the glacier through time. 132 Appendix D. Chapter 4 Supplement Figure D.1: 17-VPC1 core site, looking South. Summit of Lockie’s table in the background. Figure D.2: Sampling location at 17-VPC2 with remnant ice of informally named Lockie’s Glacier in the background. 133 Figure D.3: 17-VPC3 core site. Core collected on bedrock surface next to cobble under author’s knee in photo, between trekking poles indicating paleoflow direction (right to left in photo) of Lockie’s Glacier. Figure D.4: Sampling location of 17-VPC4. Note the cirque headwall of Lockie’s Table at the top of the image and a portion of the left lateral moraine in the top right corner of the image. 134 Figure D.5: 17-VP-07 erratic boulder. Figure D.6: Example boulder (17-VP-04) on the late Holocene moraine deposited by Lockie’s Glacier. Note the boulder sits atop a blocky, bouldery moraine draped over bedrock. 135 Figure D.7: Monte Carlo runs within measured nuclide concentrations for each site. Note the y-axis is logarithmic. Columns denoted ‘surface’ and cross-hatched are the number of solutions that are within error bounds for surface 14 C and 10Be concentrations. Columns denoted ‘all’ and dotted are the number of solutions where the entire modeled profile is within error bounds at the surface and along the 14C depth profile. 136 Figure D.8: Monte Carlo results for Core 1, showing runs that produced modelled nuclide concentrations within errors for surface 14C-10Be and 14C along the length of bedrock core. Points are colored by their χ2 value, with cooler colors denoting smaller χ2 values. The exposure and burial duration for the run with the lowest χ2 value is displayed on the figure. 137 Figure D.9: Monte Carlo results for Core 2. Please see figure caption for Fig. D.8 for an explanation of the figure. 138 Figure D.10: Monte Carlo results for Core 3. Please see figure caption for Figure D.8 for an explanation of the figure. 139 Figure D.11: Monte Carlo results for Core 4. Please see figure caption for Figure D.8 for an explanation of the figure. 140 Figure D.12: Monte Carlo results for Core 1, with 3 kyr of pre-LGM inheritance. Please see figure caption for Figure D.8 for an explanation of the figure. 141 Figure D.13: Monte Carlo results for Core 2, with 3 kyr of pre-LGM inheritance. Please see figure caption for Figure D.8 for an explanation of the figure. 142 Figure D.14: Monte Carlo results for Core 3, with 3 kyr of pre-LGM inheritance. Please see figure caption for Figure D.8 for an explanation of the figure. 143 Figure D.15: Monte Carlo results for Core 4, with 3 kyr of pre-LGM inheritance. Please see figure caption for Figure D.8 for an explanation of the figure. Table D.1: Vintage Peak sample blank information 144 Table D.2: Radiocarbon AMS results 145